Z変換とその応用　日科技連　柳井浩先生

```r

#p.78 

#パラメータ
p=0.8

#100人から指数関数的に採用人数を減少
f=function(t){
return(exp(-t)*100)  
}

#t年度の在職確率
g=function(t){
return(p^t)  
}

#f(t)g(t)の逆変換により得られたt年度の在社人数
h=function(t){
return(100*(p^t)*sum((p*exp(1))^(-c(0:t))))  
}

#退職金の見積もり
u=function(t){
return(50*(p^(t-1)-p^t)*t)  
}




```

```r

#p.26 z変換の割り算
#p.84の例 書籍：z変換とその応用　柳井浩先生

h=100*c(1,0.8,0.6,0.6)
b=c(0.8,0.5,0.3)
K=length(h)-1
N=length(b)
H=array(0,dim=c(K+1,N+1))

for(j in 1:(K+1)){
hvec=h[j:1]
H[j,1:j]=hvec
}

a=H%*%c(1,b)

```

p.85 作成メモ(何となく自身が忘れないため(笑))：

$q(t)=\beta^{t}(1-\beta) \space \space　\space$        $q(t)$はトラックの非残存率
  
$p(t)=\alpha^{t}(1-\alpha) \space \space　\space$        $p(t)$はタイヤの破損率

$p(t)$のZ変換より
$$P(z)=\frac{1-\alpha}{1-\alpha z}$$

$r(t)$を$r(t+1)=\sum_{k=0}^{t}r(k)p(t-k), r(0)=1, t=0,1,2,...$とする。このZ変換より

$$R(z)=1+\frac{(1-\alpha )z}{1-z}$$

また、この逆Z変換より$r(t)=\delta_{0}(t)+(1-\alpha)u_{1}(t)$

$v(t)=q(t)r(t)$として、$$v(t)=\beta^{t}(1-\beta)\delta_{0}(t)+(1-\alpha)\beta^{t}(1-\beta)u_{1}(t)$$

$$(※V(z)=\sum_{t=0}^{\infty}\beta^{t}(1-\beta)\delta_{0}(t)z^{t}+(1-\alpha)(1-\beta)\frac{z}{1-\beta z}
=(1-\beta)+(1-\alpha)(1-\beta)\frac{ z}{1-\beta z} )$$

$s_{1},s_{2}, ..., s_{T}$を各時刻 $1, ... ,T$ のトラックの販売数のデータとしてAR過程で計算したパラメータ$\hat{a}=(a_{1},a_{2}, ...,a_{p})$をもとに販売計画台数$s(t)=\sum_{k=1}^{p}s_{t-k}a_{k}$とし、一台のトラックあたり12本のタイヤがあるとして、$t$年度におけるタイヤの総需要本数$w(t)$を$$W(z)=12S(z)V(z)$$の逆Z変換より計算してみる。$S(z)$は$s(t)$の、$V(z)$は$v(t)$のZ変換である。Z変換の定義式：$$F(z)=\sum_{t=0}^{\infty}f(t)z^{t}$$より計算を行ってみると、$$W(z)=12S(z)V(z)=12s(t)(\frac{1-\beta}{1-z}+(1-\alpha)\beta(1-\beta)\frac{z}{(1-\beta z)(1-z)})$$
$$ =12s(t)(\frac{1-\beta}{1-z}+(1-\alpha)\beta(\frac{1}{1-z}-\frac{1}{1-\beta z}))$$
よって、
$$w(t)=12s(t)(1-\beta+(1-\alpha)\beta(1-\beta^t)) $$となったので、コードを作ってみた。

$1-\beta^t=(1-\beta)(1+\beta+\beta^2+...+\beta^{t-1})$なので各時点におけるトラックの非残存率の確率が含まれており、$1-\alpha、1-\beta$に関しても幾何分布の構造を見るとわかるように、$\alpha^{t}、\beta^t$を事象を発生しない確率とした発現確率となっている。

```r

#p.85

#タイヤの仕入れとトラック

#sのデータ
#サンプルデータ(とりあえずいいのがなかったのでこれ)
data("AirPassengers")
Y=as.numeric(t(AirPassengers))

p=6　#予測するのに使う過去のデータ数
y=c()
X=array(0,dim=c(length(Y)-p,p))
t=c()

for(j in 1:nrow(X)){
t=c(t,length(Y)-j+1)
val=Y[length(Y)-j+1]  
y=c(y,val)
vec=Y[(length(Y)-j):(length(Y)-j-p+1)]
X[j,]=vec  
}

#AR過程のパラメータ
as=solve(t(X)%*%X)%*%t(X)%*%y

#タイヤの破損数
wheelnum=c(3,5,8,9,15)

alpha=1-1/mean(wheelnum) #幾何分布の最尤推定
s0=10

r=function(t){
return(ifelse(t>0,1,0)+12*(1-alpha)*s0)  
}

#トラックの残存率計算のためのサンプル
age=c(5.8,17,15,16,14)
beta=1-1/mean(age)#幾何分布の最尤推定

#AR過程の予測値
predict=X%*%as

#w(t)の計算
wpre=c()
for(j in 1:length(t)){
w=(1-beta)+(1-alpha)*(1-beta^t[j])
wpre=c(wpre,w*12*predict[j])
}

cbind(t,wpre,predict)




```

```r


#Gauss法
#指数平滑化の予測子
#1次関数での指数平滑化 p.118
#結果は単回帰に同じ

data("AirPassengers")

Y=AirPassengers[1:9]
time=c(1:9)*30

f=function(x){
a=x[1];b=x[2];alpha=x[3]  
f0=function(t){
 return(a*t+b) 
}
f1=function(t){
return(a-b*(1-alpha)/alpha+b*t)  
}
f2=function(t){
return(a-2*b*(1-alpha)/alpha+b*t)
}
at0=a+b*time[1]
bt0=b
pthi=at0+bt0*c(time-time[1])
return(pthi)
}

h=0.001
lam=1
ite=10^3
sita=c(1,1,0.7)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}

prevec=f(sita)

plot(c(1:length(Y)),prevec,xlim=c(1,length(Y)),ylim=c(min(c(Y,prevec)),max(c(Y,prevec))),col=2,xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,prevec)),max(c(Y,prevec))),col=3,xlab="index",ylab="values")

cbind(Y,prevec)


```

```r


#Gauss法
#指数平滑化の予測子
#2次関数での指数平滑化 p.118
#結果は多項式回帰に同じ

data("AirPassengers")

Y=AirPassengers[1:12]
time=c(1:12)*30

f=function(x){
a=x[1];b=x[2];c=x[3];alpha=x[4] 
at0=a+b*time[1]+c*time[1]^2
bt0=b+2*c*time[1]
ct0=c

f1=(at0-bt0*(1-alpha)/alpha+(1-alpha)*(2-alpha)*ct0/(alpha^2))  
f2=(at0-2*bt0*(1-alpha)/alpha+2*ct0*(1-alpha)*(3-2*alpha)/(alpha^2))
f3=at0-3*(1-alpha)*bt0/alpha+3*ct0*(1-alpha)*(4-3*alpha)/(alpha^2)

pthi=at0+bt0*c(time-time[1])+ct0*c(time-time[1])^2
return(pthi)
}

h=0.001
lam=1
ite=3*10^3
sita=c(1,1,1,0.7)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
print(sum((Y-f(sita))^2))
}

prevec=f(sita)

plot(c(1:length(Y)),prevec,xlim=c(1,length(Y)),ylim=c(min(c(Y,prevec)),max(c(Y,prevec))),col=2,xlab="index",ylab="values")

par(new=T)

plot(c(1:length(Y)),Y,xlim=c(1,length(Y)),ylim=c(min(c(Y,prevec)),max(c(Y,prevec))),col=3,xlab="index",ylab="values")

cbind(Y,prevec)


```


```r

#p.138

a=0.5
b=0.6
r=0.4
delta=100
n=100
x0=c(100,200,100)
g=c(delta,0,0)
h=10^(-3)

A=matrix(c(0,b,0,a,0,r,0,0,0),nrow=3)


Ztrans=function(z){
G=g/(1-z)  
X=solve(diag(1,length(x0))-z*A)%*%x0+solve(diag(1,length(x0))-z*A)%*%(z*G)
return(X)
}


#極限定理
zpoint=0.99999
(1-zpoint)*Ztrans(zpoint)/zpoint

#微分定理：時点tでのx(t)
(Ztrans(0+10^(-4))-Ztrans(0))/(10^(-4)) #1次勾配でのx(1)
(Ztrans(0+10^(-4))-2*Ztrans(0)+Ztrans(0-10^(-4)))/(2*10^(-4*2)) #2次勾配でのx(2)


#時系列モデルの計算
xmat=array(0,dim=c(length(x0),n))
xmat[,1]=x0
for(j in 1:(n-1)){  
if(j==1){
x=A%*%x0+g
xmat[,2]=x
}else{
x=A%*%x+g 
xmat[,(j+1)]=x
}}



```