LoginSignup
2
0

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

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

Last updated at Posted at 2018-01-13

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

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

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

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

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


#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)))





#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



#状態空間モデル

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による実践ガイド　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)

}
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0