時系列解析の勉強3

forecast、seasonal パッケージ　X13-ARIMA-SEATS

(参考)
R {seasonal} パッケージで X13-ARIMA-SEATSを使う
Examples of X 13ARIMA SEATS in R

用いたデータ
牛乳等生産量2002.1_2014.9

• 牛乳乳製品の生産動向 月次データ(政府統計の総合窓口)

#時系列データに変換して保存しておく。
#m<-scan("milk.txt")
#milk<-ts(m, frequency =12, start = c(2002,1))
#save(milk,file="milk.dat")

milk.dat

load("milk.dat") #牛乳等生産量2002.1_2014.9
#library(TTR) #移動平均(ここでは使わない）
#png("milk01.png",  width=1000, height=800)
par(mar=c(5,5,5,5),las=1,family="Takao P明朝")
plot(milk,main="牛乳等生産量 (2002.1_2014.9)")
grid()
#lines(SMA(milk, 3),col="red",lty=2)
#dev.off()

library(forecast)
#時系列処理/季節調整プログラム X13-ARIMA-SEATSを使う
library(seasonal)
Sys.setenv(X13_PATH = "/home/user")
checkX13()  #使用可能状態になっているか確認
fit<- auto.arima(milk,trace=T,stepwise=F,approximation=F,start.p=0,
  start.q=0,start.P=0,start.Q=0)
fit

Series: milk
ARIMA(0,1,1)(2,0,0)[12] with drift

Coefficients:
ma1 sar1 sar2 drift
-0.4426 0.4598 0.5016 -9.5177
s.e. 0.0752 0.0777 0.0798 2407.6606

sigma^2 estimated as 29396153: log likelihood=-1527.02
AIC=3064.03 AICc=3064.44 BIC=3079.15

d.forecast <- forecast(fit, level = c(95), h = 36)
m <- seas(milk,arima.model = "(0 1 1)(2 0 0)",regression.variables = c("const", "td"))
d.season<-series(m, c('forecast.forecasts'))
#png("milk02.png",  width=1000, height=800)
par(mar=c(5,5,5,8),las=1,family="Takao P明朝")
plot(forecast(d.forecast,h=36),main="予測  (forecast の結果  &  seasonal の結果)")
grid()
lines(d.season[,1],col="red")
axis(4,d.season[nrow(d.season),1], "seasonal package",col.axis="red")
axis(4,d.forecast$mean[length(d.forecast$mean)],"forecast package",col.axis="blue")
#dev.off()