AR過程（TIME SERIES,GRENANDER-ROSENBRAT, CHELSEA PUBLISHING COMPANY)

AR過程をpackageを使わず、Yule-Walker法を用いてRで作成してみました。

下記はsin(x)をサンプルデータとして、それをAR過程で予測してみたものです。

一人で書いているだけなのもさみしいので、投稿してみました！

ある意味での情報共有の場にしたいです！

間違っていたりしたら、情報共有頂けると助かります。

```{r}

#AR過程

library(dplyr)
library(pryr)
N=100
data=data.frame(num=1:N) %>% dplyr::mutate(value=num*sin(num)+cos(num)*num^2)
data=data %>% dplyr::mutate(normalize_value=value)
plot(data$num,data$normalize_value,type="p")
mean=mean(data$value);var=var(data$value)
vector=data$normalize_value

w<-ts(vector)
v=acf(w)
autocorrelation_matrix=array(0,dim=c((length(v$acf)-1),(length(v$acf)-1)))

for(j in 1:ncol(autocorrelation_matrix)){
  vec=v$acf[1:(length(v$acf)-1)]
  vec=vec[1:(length(vec)-(j-1))]
  autocorrelation_matrix[j,j:(ncol(autocorrelation_matrix))]=vec
}
 autocorrelation_matrix2=t(autocorrelation_matrix)
 diag(autocorrelation_matrix2)<-0
 autocorrelation_matrix=autocorrelation_matrix+autocorrelation_matrix2 
vec=v$acf[2:length(v$acf)]
 alpha=solve(autocorrelation_matrix,vec)
 data=data %>% dplyr::mutate(estimate_value=0)
error=0
for(j in (length(alpha)+1):nrow(data)){
  val=(data$normalize_value[j]-sum(alpha*data$normalize_value[(j-1):(j-length(alpha))]))^2
  error=error+val
  data$estimate_value[j]=sum(alpha*data$normalize_value[(j-1):(j-length(alpha))])
}
stn_error=error/(nrow(data)-length(alpha))
estimated_variance=v$acf[1]-sum(alpha*v$acf[2:length(v$acf)])
if(floor(stn_error*10)/10==floor(estimated_variance*10)/10){
  print("Result is right!")
}

plot(data$num,data$normalize_value,xlab = "time",ylab="value",main = "yule_walker",xlim=c(min(data$num),max(data$num)),ylim=c(min(data$normalize_value),max(data$normalize_value)),col=2)
par(new=T)
plot(data$num,data$estimate_value,xlab = "time",ylab="value",main = "yule_walker",xlim=c(min(data$num),max(data$num)),ylim=c(min(data$normalize_value),max(data$normalize_value)),col=3)

log_likelihood=-(nrow(data)/2)*log(2*pi)-(nrow(data)/2)*log(estimated_variance)-(1/(2*estimated_variance^2))*error
spectrum=function(w){
  y<-estimated_variance/(abs(1-sum(alpha[2:length(alpha)]*exp(-(1i)*w*c(1:(length(alpha)-1)))))^2)
return(y)
}
spectral=data.frame(num=seq(0,100,by=0.1)) %>% dplyr::mutate(value=0)
for(j in 1:nrow(spectral)){
  spectral$value[j]=spectrum(spectral$num[j])
}
plot(spectral$num,spectral$value,type="l")

spectral_density_func=function(x){
  
val=c()

A=c(-1,alpha)

for(j in 0:length(alpha)){

val=c(val,exp(1i*x)^(length(alpha)-j))  
  
}
  
a=Re(sum(A*val));b=Im(sum(A*val))

return(estimated_variance/(2*pi*Re((a+b*1i)*(a-b*1i)))  )

}

spectral_dens_data=data.frame(x=seq(-pi,pi,0.01),df=0)

for(j in 1:nrow(spectral_dens_data)){
  
spectral_dens_data$df[j]=spectral_density_func(spectral_dens_data$x[j])  
  
}

plot(spectral_dens_data$x,spectral_dens_data$df,type="o",xlab="lambda",ylab="spectral_density_func")

#simpson積分公式による積分値の計算(一変数関数の積分)

a=-pi;H=2*pi;n=100000

Y=(H/3)*(1/(n-1))*(spectral_density_func(a)+spectral_density_func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*spectral_density_func(a+(j*H)/(n-1))

}
#simpson積分値(spectral density functionが全体の積分で1になっている)
print(paste0("スペクトル密度関数全体の全区間の積分は",Y))

#entropy rate の計算

log_spectral_density_func=function(s){
  
z<-log(4*pi*pi*exp(1)*spectral_density_func(s))  
  
 return(z) 
}

#simpson積分公式による積分値の計算(一変数関数の積分)

a=-pi;H=2*pi;n=100000

Y=(H/3)*(1/(n-1))*(log_spectral_density_func(a)+log_spectral_density_func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*log_spectral_density_func(a+(j*H)/(n-1))

}

#simpson積分値(entropy rate)
print(paste0("エントロピーレートは",Y/(4*pi)))



```

```{r}

#AR過程

library(dplyr)

data=data.frame(no=1:48,values=window(UKgas,c(1975,1),c(1999,3)))

r=4;

values=data$values

X=array(0,dim=c(nrow(data)-r,r))

Y=c()

for(j in 1:(nrow(data)-r)){
  
X[j,]=c(values[j:(j+r-1)])

Y=c(Y,(values[r+j]))
  
}

A=solve(t(X)%*%X)%*%t(X)%*%Y

b=mean(Y)-sum(apply(X,2,mean)*A)

predict=X%*%A+b

result=data.frame(predict=predict,Y=Y)

N=nrow(X);

sigma_hat=sum(c(Y-X%*%A+b)^2)/N

AIC=log(sigma_hat)+2*r/N

```


```{r}

#状態空間モデル

library(dplyr)

data=data.frame(no=1:50,round(LifeCycleSavings))

mat=matrix(sample(1:100,nrow(data)*(ncol(data)-1),replace=T),ncol=ncol(data)-1,nrow=nrow(data))

w=4;t=nrow(data)


X=t(as.matrix(mat*data[1:(t-w),2:5]))


Z=t(as.matrix(data[1:(t-w),2:5]))

Z_plus=t(as.matrix(data[2:(t-w),2:5]))

Z_minus=t(as.matrix(data[1:(t-w-1),2:5]))

mat2=Z_minus%*%t(Z_minus)

A=Z_plus%*%t(Z_minus)%*%solve(mat2)

C=((X)%*%t(Z))%*%solve(Z%*%t(Z))

R=(X-C%*%Z)%*%t(X-C%*%Z)/(t-w)

Q=(Z_plus-A%*%Z_minus)%*%t(Z_plus-A%*%Z_minus)/(t-w-1)

log_likelihood=-(t-w-1)*(ncol(A)*log(2*pi)+log(det(Q)))/2-sum(diag(solve(Q)%*%(Z_plus-A%*%Z_minus)%*%t(Z_plus-A%*%Z_minus)))/2

cor(c(C%*%Z),c(X))

plot(1:length(c(X)),c(C%*%Z),type="o",col=2,xlim=c(1,t-w),ylim=c(min(c(C%*%Z,X)),max(c(C%*%Z,X))),xlab="index",ylab="values")

par(new=T)

plot(1:length(c(X)),c(X),type="o",col=3,xlim=c(1,t-w),ylim=c(min(c(C%*%Z,X)),max(c(C%*%Z,X))),xlab="index",ylab="values")



```

```{r}

#入門 機械学習による異常検知―Rによる実践ガイド　p.219

#7.4.3 部分空間同定法

#kalman-filter

library(dplyr)

data=data.frame(no=1:50,round(LifeCycleSavings))

mat=matrix(sample(1:100,nrow(data)*(ncol(data)-1),replace=T),ncol=ncol(data)-1,nrow=nrow(data))

w=4;t=nrow(data)

Z=t(as.matrix(data[1:(t-w),2:5]))

Z_plus=t(as.matrix(data[2:(t-w),2:5]))

Z_minus=t(as.matrix(data[1:(t-w-1),2:5]))

X=t(as.matrix(mat*data[1:(t-w),2:5]))


M=nrow(X);N=M*w+1;time=25

A=Z_plus%*%t(Z_minus)%*%solve(Z_minus%*%t(Z_minus))

C=((X)%*%t(Z))%*%solve(Z%*%t(Z))

R=(X-C%*%Z)%*%t(X-C%*%Z)/(t-w)

Q0=(Z_plus-A%*%Z_minus)%*%t(Z_plus-A%*%Z_minus)/(t-w-1)

z0=as.numeric(data[1,2:5]);mu0=solve(A)%*%c(z0)

Xp=c(X[,(time-1):(time-1+w-1)])

Xf=c(X[,(time):(time+w-1)])

for(j in 1:(N-1)){

Xp_sub=c(X[,(time-1-j):(time-1-j+w-1)])

Xf_sub=c(X[,(time+j):(time+j+w-1)])

Xp=cbind(Xp_sub,Xp)

Xf=cbind(Xf,Xf_sub)

}

Spf=Xp%*%t(Xf)

Spp=Xp%*%t(Xp)

Sff=Xf%*%t(Xf)

Sfp=Xf%*%t(Xp)

diag(solve(Spp)%*%Spf%*%solve(Sff)%*%Sfp)

W=solve(sqrt(Spp))%*%Spf%*%solve(sqrt(Sff))

u=svd(W)$u

d=svd(W)$d

v=svd(W)$v

#alphaと同じ

eigen(W%*%t(W))$vectors

#betaと同じ

eigen(t(W)%*%W)$vectors

m=6

alpha=u[,1:m]

#ラグランジュ未定乗数の条件
t(alpha)%*%alpha

beta=v[,1:m]

#ラグランジュ未定乗数の条件
t(beta)%*%beta

alpha_hat=(solve(sqrt(Spp)))%*%alpha

beta_hat=(solve(sqrt(Sff)))%*%beta

Zp_hat=t(alpha_hat)%*%Xp

Zf_hat=t(beta_hat)%*%Xf






for(j in 1:ncol(X)){

if(j==1){  

K=Q0%*%t(C)%*%solve(R+C%*%Q0%*%t(C))

mu=A%*%mu0+K%*%c(X[,j]-C%*%A%*%mu0)

V=(diag(1,m)-K%*%C)%*%Q0

Q=Q0+A%*%V%*%t(A)

}else{

K=Q%*%t(C)%*%solve(R+C%*%Q%*%t(C))

mu=A%*%mu+K%*%c(X[,j]-C%*%A%*%mu)

V=(diag(1,m)-K%*%C)%*%Q

Q=Q+A%*%V%*%t(A)  


}

#print("predict:")
#print((C%*%mu))

#print("variance:")
#print(V)

print("kalman matrix:")
print(K)

}

  

```