RStudio

packageを使用しないでやってみた(ARMA)

ARMAを勉強のため、packageを使わずに実装してみました。
よければ一度、動かしてみてください。
変なところがあったら、教えてください。

#ARMA(修正) 

library(dplyr)
library(pryr)

Data=read.csv("C:/Users/kozakai/Desktop/arima_sample/arima/data/unemploy_train.csv")

data=Data %>% mutate(num=1:nrow(Data))
colnames(data)=c("Month","value","num")
plot(data$num,data$value,type="p",xlab="num",ylab="value")



p=2;r=0
print(paste0("num(a)=",p-1," ;num(b)=",r))
P=1
val=c()
pthi=rep(0,p+r-1)
pthi_data=matrix(0,ncol=p+r-1,nrow=nrow(data)-p+1)
V=c()
V=matrix(0,ncol=p+r-1,nrow=nrow(data)-p+1)
data=data %>% dplyr::mutate(predict=0)

sigma_vec=c()

for(j in p:nrow(data)){

  val=c(val,P)

  if(p==2){

  sigma=0.01  

  }else{

  if(var(data$value[(j-p+1):(j-1)])>0){

  sigma=var(data$value[(j-p+1):(j-1)])

  }else{

  if(j>p){

  sigma=sigma_vec[length(sigma_vec)] 

  }else{

  sigma=1  

  }

  }

  }  

  sigma_vec=c(sigma_vec,sigma)

  vec=c(t(data$value[(j-p+1):(j-1)]),rnorm(r,mean=0,sd=sqrt(sigma)))
  S=(sum(vec^2)*(1/sigma)+P^(-1))^(-1)
  pthi=S*(pthi/P+(data$value[j]*vec)/(sigma))
  P=S;pthi_data[j-p+1,]=pthi
  V[j-p+1,]=vec
  data$predict[j]=sum(pthi*vec)

}
P=val

#learning(b0=1)

colnames(data)=c("Month","value","num","learning")

#iteration
plot(data$num,data$value,xlab="time",ylab="value",main="Iteration",xlim=c(min(data$num),max(data$num)),ylim=c(min(data$learning),max(data$learning)),col=2)
par(new=T)
plot(data$num,data$learning,xlab="time",ylab="value",main="Iteration",xlim=c(min(data$num),max(data$num)),ylim=c(min(data$learning),max(data$learning)),col=3)

#predict

data=data %>% dplyr::mutate(predict=0)

for(j in p:nrow(data)){

a=pthi[1:(p-1)];b=pthi[p:length(pthi)]  

data$predict[j]=sum(data$value[(j-p+1):(j-1)]*a)+sum(b*rnorm(r,mean=0,sd=sqrt(sigma)))+rnorm(1,mean=0,sd=sqrt(sigma))

}

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

cor(data$value,data$predict)

unemploy_train.csv ;

image.png