# 季節とトレンド library(dlm) data <- log(AirPassengers) plot(data,type="o") # レベルも傾きも、季節変動も、システムノイズ0 観測誤差のみ # Step1 # モデル作成のための関数を作る build.4<-function(theta){ dlmModPoly(order=2,dV=exp(theta[1]),dW=c(0,0))+ dlmModSeas(fr=12,dW=c(0,rep(0,10)),dV=0) } # Step2 # MLEでパラメタ推定。 fit.4<-dlmMLE(data,parm=dlmMLE(data,parm=c(1),build.4,method="Nelder-Mead")$par,build.4,method="BFGS") fit.4 # 推定されたパラメタを使ってモデルを作り直す DLM.4<-build.4(fit.4$par) DLM.4 DLM.4$V DLM.4$F # Step3 # カルマンフィルター Filt.4<-dlmFilter(data,DLM.4) Filt.4$m # Step4 # スムージング Smooth.4<-dlmSmooth(Filt.4) # plot plot(data,col=1,type="o") lines(dropFirst(Filt.4$m)[, 1] + dropFirst(Filt.4$m)[, 3],col=2) lines(dropFirst(Smooth.4$s)[, 1] + dropFirst(Smooth.4$s)[, 3], col=4) legend("bottomright", pch=c(1,NA,NA), col=c(1,2,4), lwd=1, legend=c("data","Filter","Smooth")) par(mfrow=c(3, 1)) # 元データ plot(data, col=1, type="o", main="data") # レベル+トレンド成分 plot(dropFirst(Smooth.4$s)[, 1], col=4, main="level+trend") #季節成分 plot(dropFirst(Smooth.4$s)[, 3], col=4, main="seasonal") par(mfrow=c(1, 1)) #======================================================================= # システムノイズもパラメタ推定 # さらに、1959年以降を予測する # 1959年以降を切り落とす test.data <- window(data,end=c(1958,12)) test.data # Step1 # モデル作成のための関数を作る build.5<-function(theta){ dlmModPoly(order=2,dV=exp(theta[1]),dW=c(exp(theta[2]),exp(theta[3])) )+ dlmModSeas(fr=12,dW=c(theta[4],rep(0,10)),dV=0) } # Step2 # MLEでパラメタ推定。 fit.5<-dlmMLE(test.data,parm=dlmMLE(test.data,parm=c(0,1,1,1),build.5,method="Nelder-Mead")$par,build.5,method="BFGS") fit.5 # 推定されたパラメタを使ってモデルを作り直す DLM.5<-build.5(fit.5$par) DLM.5 # Step3 # カルマンフィルター Filt.5<-dlmFilter(test.data,DLM.5) # Step4 # スムージング Smooth.5<-dlmSmooth(Filt.5) # plot plot(test.data,col=1,type="o") lines(dropFirst(Filt.5$m)[,1]+dropFirst(Filt.5$m)[,3],col=2) lines(dropFirst(Smooth.5$s)[,1]+dropFirst(Smooth.5$s)[,3],col=4) legend("bottomright", pch=c(1,NA,NA), col=c(1,2,4), lwd=1, legend=c("data","Filter","Smooth")) par(mfrow=c(3,1)) # 元データ plot(test.data,col=1,type="o", main="data") # レベル+トレンド成分 plot(dropFirst(Smooth.5$s)[,1],col=4, main="level+trend") # 季節成分 plot(dropFirst(Smooth.5$s)[,3],col=4, main="seasonal") par(mfrow=c(1,1)) # 予測 Fore <- dlmForecast(Filt.5,nAhead=24,sampleNew=5) # 予測の答え合わせ plot(data, type="o") lines(dropFirst(Smooth.5$s)[,1]+dropFirst(Smooth.5$s)[,3],col=4) lines(Fore$f,col=2,lwd=2) legend("bottomright",pch=c(1,NA),col=c(1,2),lwd=c(1,2),legend=c("実測値","予測値")) # 乱数を使って予測 Line <- function(x){ lines(x, col=8, type="o") } plot(window(test.data,start=c(1957,1)), xlim=c(1957,1961),ylim=c(5.7,6.5), type="o") lapply(Fore$newObs, Line) lines(window(data,start=c(1959,1)), col=1,lwd=2) lines(Fore$f,col=2) legend("topleft",pch=c(NA,NA,1),col=c(1,2,8),lwd=c(2,1,1),legend=c("実測値","予測値","乱数を使った予測値"))