## UKgas データのモデリング
# ガス消費量を四半期ごとに観測した時系列データ
data(UKgas)
u <- UKgas
plot(u) # 観測値のプロット
library("forecast")
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: timeDate
## This is forecast 5.6
# stepwise=F で網羅的に調べてくれる
# trace=T で実行内容出力
auto.arima(u,stepwise = F,trace=T) -> fit.arima
##
## ARIMA(0,1,0)(0,1,0)[4] : 1061.603
## ARIMA(0,1,0)(0,1,1)[4] : 1062.474
## ARIMA(0,1,0)(0,1,2)[4] : 1064.484
## ARIMA(0,1,0)(1,1,0)[4] : 1066.418
## ARIMA(0,1,0)(1,1,1)[4] : 1068.54
## ARIMA(0,1,0)(1,1,2)[4] : Inf
## ARIMA(0,1,0)(2,1,0)[4] : 1072.662
## ARIMA(0,1,0)(2,1,1)[4] : 1074.83
## ARIMA(0,1,0)(2,1,2)[4] : 1070.62
## ARIMA(0,1,1)(0,1,0)[4] : 991.3633
## ARIMA(0,1,1)(0,1,1)[4] : 993.4824
## ARIMA(0,1,1)(0,1,2)[4] : 995.3999
## ARIMA(0,1,1)(1,1,0)[4] : 996.8972
## ARIMA(0,1,1)(1,1,1)[4] : Inf
## ARIMA(0,1,1)(1,1,2)[4] : Inf
## ARIMA(0,1,1)(2,1,0)[4] : 1003.22
## ARIMA(0,1,1)(2,1,1)[4] : Inf
## ARIMA(0,1,1)(2,1,2)[4] : Inf
## ARIMA(0,1,2)(0,1,0)[4] : 991.8972
## ARIMA(0,1,2)(0,1,1)[4] : 994.0086
## ARIMA(0,1,2)(0,1,2)[4] : 995.9851
## ARIMA(0,1,2)(1,1,0)[4] : 997.606
## ARIMA(0,1,2)(1,1,1)[4] : 999.7985
## ARIMA(0,1,2)(1,1,2)[4] : 998.8268
## ARIMA(0,1,2)(2,1,0)[4] : 1004.351
## ARIMA(0,1,2)(2,1,1)[4] : 1003.061
## ARIMA(0,1,3)(0,1,0)[4] : 989.9013
## ARIMA(0,1,3)(0,1,1)[4] : 991.609
## ARIMA(0,1,3)(0,1,2)[4] : 993.8635
## ARIMA(0,1,3)(1,1,0)[4] : 995.3873
## ARIMA(0,1,3)(1,1,1)[4] : 997.6109
## ARIMA(0,1,3)(2,1,0)[4] : 1002.871
## ARIMA(1,1,0)(0,1,0)[4] : 1040.247
## ARIMA(1,1,0)(0,1,1)[4] : 1042.176
## ARIMA(1,1,0)(0,1,2)[4] : 1044.042
## ARIMA(1,1,0)(1,1,0)[4] : 1046.119
## ARIMA(1,1,0)(1,1,1)[4] : Inf
## ARIMA(1,1,0)(1,1,2)[4] : 1040.29
## ARIMA(1,1,0)(2,1,0)[4] : 1051.862
## ARIMA(1,1,0)(2,1,1)[4] : Inf
## ARIMA(1,1,0)(2,1,2)[4] : Inf
## ARIMA(1,1,1)(0,1,0)[4] : 993.9869
## ARIMA(1,1,1)(0,1,1)[4] : 996.0927
## ARIMA(1,1,1)(0,1,2)[4] : 998.0513
## ARIMA(1,1,1)(1,1,0)[4] : 999.5489
## ARIMA(1,1,1)(1,1,1)[4] : 1001.476
## ARIMA(1,1,1)(1,1,2)[4] : 1003.639
## ARIMA(1,1,1)(2,1,0)[4] : 1005.205
## ARIMA(1,1,1)(2,1,1)[4] : Inf
## ARIMA(1,1,2)(0,1,0)[4] : 994.2836
## ARIMA(1,1,2)(0,1,1)[4] : 996.4936
## ARIMA(1,1,2)(0,1,2)[4] : 998.6107
## ARIMA(1,1,2)(1,1,0)[4] : Inf
## ARIMA(1,1,2)(1,1,1)[4] : 1001.837
## ARIMA(1,1,2)(2,1,0)[4] : 1005.879
## ARIMA(1,1,3)(0,1,0)[4] : Inf
## ARIMA(1,1,3)(0,1,1)[4] : Inf
## ARIMA(1,1,3)(1,1,0)[4] : Inf
## ARIMA(2,1,0)(0,1,0)[4] : 1017.687
## ARIMA(2,1,0)(0,1,1)[4] : 1019.851
## ARIMA(2,1,0)(0,1,2)[4] : 1022.045
## ARIMA(2,1,0)(1,1,0)[4] : 1023.777
## ARIMA(2,1,0)(1,1,1)[4] : 1025.988
## ARIMA(2,1,0)(1,1,2)[4] : Inf
## ARIMA(2,1,0)(2,1,0)[4] : 1030.223
## ARIMA(2,1,0)(2,1,1)[4] : Inf
## ARIMA(2,1,1)(0,1,0)[4] : 992.0569
## ARIMA(2,1,1)(0,1,1)[4] : 993.3738
## ARIMA(2,1,1)(0,1,2)[4] : 995.6026
## ARIMA(2,1,1)(1,1,0)[4] : 997.0855
## ARIMA(2,1,1)(1,1,1)[4] : 999.3362
## ARIMA(2,1,1)(2,1,0)[4] : 1003.471
## ARIMA(2,1,2)(0,1,0)[4] : Inf
## ARIMA(2,1,2)(0,1,1)[4] : 994.6148
## ARIMA(2,1,2)(1,1,0)[4] : 992.4116
## ARIMA(2,1,3)(0,1,0)[4] : Inf
## ARIMA(3,1,0)(0,1,0)[4] : 1012.752
## ARIMA(3,1,0)(0,1,1)[4] : 990.4821
## ARIMA(3,1,0)(0,1,2)[4] : 991.402
## ARIMA(3,1,0)(1,1,0)[4] : 1007.507
## ARIMA(3,1,0)(1,1,1)[4] : 994.738
## ARIMA(3,1,0)(2,1,0)[4] : 1005.787
## ARIMA(3,1,1)(0,1,0)[4] : 995.5411
## ARIMA(3,1,1)(0,1,1)[4] : 988.7248
## ARIMA(3,1,1)(1,1,0)[4] : 998.8104
## ARIMA(3,1,2)(0,1,0)[4] : Inf
fit.arima # AR(3),MA(1),sma(1)
## Series: u
## ARIMA(3,1,1)(0,1,1)[4]
##
## Coefficients:
## ar1 ar2 ar3 ma1 sma1
## -0.6938 -0.9083 -0.6588 -0.2470 -0.7296
## s.e. 0.1070 0.0691 0.1151 0.1289 0.1080
##
## sigma^2 estimated as 1061: log likelihood=-506.42
## AIC=1024.84 AICc=1025.72 BIC=1040.65
# もちろん、自分で次数の設定もできる
# arima(u,order=c(3,1,1),seasonal = list(order=c(0,1,1))) -> fit.arima.own
# frequery を持つtimeseriesオブジェクトだとSARIMA(p,d,q)(pp,dd,qq)としてくれる。
forecast(fit.arima,level=c(95),h=8) -> predict.fc
plot(predict.fc)
(predict.fc)
## Point Forecast Lo 95 Hi 95
## 1987 Q1 1197.6804 1133.8501 1261.5106
## 1987 Q2 675.7573 611.8153 739.6993
## 1987 Q3 368.1837 303.0241 433.3432
## 1987 Q4 824.8829 758.7539 891.0118
## 1988 Q1 1235.5669 1141.8586 1329.2752
## 1988 Q2 724.7930 630.1863 819.3997
## 1988 Q3 399.2647 303.3632 495.1662
## 1988 Q4 861.0579 763.7638 958.3520
# SのないARIMAモデルの場合
library("magrittr")
u %>%
ts %>%
auto.arima(stepwise = F) -> fit.arima2
fit.arima2 # arima(u,order=c(2,1,0)) -> fit.arima.own
## Series: .
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## -0.0244 -0.9771 5.4433
## s.e. 0.0220 0.0160 2.5178
##
## sigma^2 estimated as 2667: log likelihood=-571.07
## AIC=1150.15 AICc=1150.54 BIC=1160.84
fit.arima2 %T>%
summary %>%
forecast(h=8) %>% plot
## Series: .
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## -0.0244 -0.9771 5.4433
## s.e. 0.0220 0.0160 2.5178
##
## sigma^2 estimated as 2667: log likelihood=-571.07
## AIC=1150.15 AICc=1150.54 BIC=1160.84
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.150833 51.15932 36.94186 -3.126375 13.07156 0.2103359
## ACF1
## Training set -0.6919714
