R Time Series Analysis 時系列解析(4) (HoltWinters法)

#http://www.esri.cao.go.jp/jp/sna/data/data_list/sokuhou/files/2013/qe131_2/gdemenuja.html
#より取得し、加工したデータを使用
SNA <- read.csv("GDP_def_extract.csv",header=TRUE)
View(SNA)
GDP <- ts(SNA$GDP, start=c(1980,1), frequency=4)
  
GDP.HW1 <- HoltWinters(GDP,seasonal="multiplicative")
GDP.HW1
#Smoothing parameters:最小逐次推定二乗誤差を与えるパラメータの値
#Coefficients:aは水準,bは傾向,sjは季節成分

#T期におけるk期先(第T+k期)の予測は、次のようになる
#seasonal="additive"なら yt(k) = a(T)+k*b(T)+sT+k-j (1≦k≦j),  a(T)+k*b(T)+sT+k-2j (j<k≦2j)
#seasonal="multiplicative"なら yt(k) = (a(T)+k*b(T))*sj+k-j (1≦k≦j), (a(T)+k*b(T))*sj+k-2j (j<k≦2j)
#のように予測する

GDPfit1 <- GDP.HW1$fitted
#xhatは予測値
#levelはa(T),trendはb(T),seasonはsT+k-jに対応。
#seasonal = "multiplicative"としているので
#xhat = (level+trend)*seasonの関係が成立
#GDP.WH1はmatrix形式であることに注意

title="ホルトウィンタース法による近似曲線(2000年4-6月=100)"
ts.plot(GDP,GDPfit1[,1],type="l",col=c(1:2),lty=c(1:2),main=title)

##季節要素を含めない場合

#decomposeにより、傾向成分(TC),季節成分(S),不規則成分(I)を分解
bunkai1 <- decompose(GDP,type="multiplicative")
GDPTCI <- GDP/bunkai1$seasonal
GDP.HW2 <- HoltWinters(GDPTCI,gamma=0) #gamma=0は季節成分を含めない設定
GDPfit2 <- GDP.HW2$fitted

ts.plot(GDP,GDPfit1[,1],GDPfit2[,1],type="l",col=c(1:3),lty=c(1:3),main=title)

##predict関数を用いて予測
oldpar <- par(no.readonly=TRUE)
par(mfrow=c(2,1),mar=c(2,4,4,2),par(new=F))

#積型
pre1 <- predict(GDP.HW1,n.ahead=12)
title="季節要素のあるHoltWintersモデルによる予測(積型)"
ts.plot(GDP,pre1,type="l",lty=c(1,2),col=c(1,2),main=title,xlab="")

#積型(季節要素を含めない)
pre2 <- predict(GDP.HW2,n.ahead=12)
title="季節要素のないHoltWintersモデルによる予測(積型)"
ts.plot(GDP,pre2,type="l",lty=c(1,2),col=c(1,2),main=title,xlab="")
par(mfrow=c(1,1)); par(oldpar)