SARIMAXモデルによる時系列データの予測

時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装第3部 1章

この章で使うパッケージ

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(ggplot2)
library(ggfortify)

## Registered S3 methods overwritten by 'ggfortify':
##   method                 from    
##   autoplot.Arima         forecast
##   autoplot.acf           forecast
##   autoplot.ar            forecast
##   autoplot.bats          forecast
##   autoplot.decomposed.ts forecast
##   autoplot.ets           forecast
##   autoplot.forecast      forecast
##   autoplot.stl           forecast
##   autoplot.ts            forecast
##   fitted.ar              forecast
##   fortify.ts             forecast
##   residuals.ar           forecast

library(tseries)
library(urca)

分析の対象となるデータ

イギリスの交通事故死傷者数のデータを用いた

R組み込みのデータ

dfSeatbelts <- Seatbelts
head(dfSeatbelts)

##          DriversKilled drivers front rear   kms PetrolPrice VanKilled law
## Jan 1969           107    1687   867  269  9059   0.1029718        12   0
## Feb 1969            97    1508   825  265  7685   0.1023630         6   0
## Mar 1969           102    1507   806  319  9963   0.1020625        12   0
## Apr 1969            87    1385   814  407 10955   0.1008733         8   0
## May 1969           119    1632   991  454 11823   0.1010197        10   0
## Jun 1969           106    1511   945  427 12391   0.1005812        13   0

front:前席における死傷者数 PetrolPrice:ガソリンの値段 law:前席においてシートベルトを装着することを義務付けた法律の施行の有無を表したフラグ。0だと未施工、1だと施工済み。1983年の1月31日に施工された。

前方座席の死傷者数のみ抽出した

front <- Seatbelts[, "front"]
front

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1969  867  825  806  814  991  945 1004 1091  958  850 1109 1113
## 1970  925  903 1006  892  990  866 1095 1204 1029 1147 1171 1299
## 1971  944  874  840  893 1007  973 1097 1194  988 1077 1045 1115
## 1972 1005  857  879  887 1075 1121 1190 1058  939 1074 1089 1208
## 1973  903  916  787 1114 1014 1022 1114 1132 1111 1008  916  992
## 1974  731  665  724  744  910  883  900 1057 1076  919  920  953
## 1975  664  607  777  633  791  790  803  884  769  732  859  994
## 1976  704  684  671  643  771  644  828  748  767  825  810  986
## 1977  714  567  616  678  742  840  888  852  774  831  889 1046
## 1978  889  626  808  746  754  865  980  959  856  798  942 1010
## 1979  796  643  794  750  809  716  851  931  834  762  880 1077
## 1980  748  593  720  646  765  820  807  885  803  860  825  911
## 1981  704  691  688  714  814  736  876  829  818  942  782  823
## 1982  595  673  660  676  755  815  867  933  798  950  825  911
## 1983  619  426  475  556  559  483  587  615  618  662  519  585
## 1984  483  434  513  548  586  522  601  644  643  641  711  721

前方座席の死傷者数のデータを図示した

ggplot2を使用した

autoplot(
  Seatbelts[, "front"],
  main="イギリスの交通事故死傷者数(前席)",
  xlab="年",
  ylab="死傷者数"
)

・交通事故死亡者数には年単位の周期性がありそう
・1974年と1983年ごろに2回ほど死亡者数が減っている

差分系列の作成方法

単位根検定（KPSS検定）

データの差分を取る必要があるかどうかを調べる。
yが非定常過程でその差分系列が定常過程である時、yは単位根過程である。
xとyに回帰分析を行うと、まったく関係のないxとyの間に有意な相関を見出してしまう「見せかけの回帰」が行われることがよくある。
「見せかけの回帰」を避けるためには、あらかじめ２つの時系列データxとyが単位根過程に従っているかどうかを確認しなくてはならない。それが単位根検定。
ここで対数を取ってから検定を行った理由は、対数変換したデータを対象に分析するため

summary(ur.kpss(log(Seatbelts[, "front"])))

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 2.3004 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

有意水準5%の棄却点が0.463であり、今回のデータの統計量2.3004が棄却点を上回っていることがわかった
単位根がないという帰無仮説が棄却されたので、単位根があるとみなすことができる

差分をとる回数を調べる

本来はここで差分をとってからもう一度KPSS検定をという流れになるが、forecastのndiffs関数を使うと、差分を取る回数がすぐにわかる

ndiffs(log(Seatbelts[, "front"]))

## [1] 1

対数変換したデータを図示

個数や人数などのデータは対数変換してからモデル化すると上手くモデル化できる傾向がある。
例えば、法律の施行は交通事故死傷者数を0.5倍にした。など掛け算で表現するほうが自然でもある。
図示にforecastパッケージのggtsdisplay関数を使用した。

log_front <- log(front)
ggtsdisplay(log_front, main="対数系列")

ACF:相関係数
PACF:偏相関係数

対数変換をしたことで、年毎の変動の幅がおよそ均一となった。
自己相関は長期にわたって続いていることがわかる。また12か月周期で大きくなっている。
偏自己相関も、およそ12か月単位で大きくなっている。

対数差分系列を作成して図示

diff関数で差分系列を作成できる
デフォルトでlag=1となる

log_diff <- diff(log_front)
ggtsdisplay(log_diff, main="対数差分系列")

長期にわたって平均値が変化せず、単位根がなくなっていることがわかる。
短期の自己相関は減ったが、不規則に絶対値の大きな自己相関がある。12か月単位での自己相関も目立つ。
一年単位での自己相関が目立ったため、季節成分があることは疑いようがない。

季節差分系列の作成方法

forecastパッケージのggsubseriesplot関数を用いて、月毎に分けたグラフを描いてみる

ggsubseriesplot(front)

12月が最も交通事故死傷者数が多い月であることがわかった

次は、季節差分を取ってみる

対数差分系列をさらに季節差分をとって図示

lagの値をfrequency関数で求めた。frequency(front):12。lag=12でもよい。

seas_log_diff <- diff(log_diff, lag=frequency(log_diff))
ggtsdisplay(seas_log_diff, main="季節差分系列")

12か月単位の自己相関が大きく残ったままとなっている。
季節階差を取ったからといって、季節の影響をすべて取り除けるわけではないようだ。

訓練データとテストデータに分ける

予測の評価のために、訓練データとテストデータに分割する。
あらかじめ対数変換もしておく。
ARIMAモデルは内部で勝手に差分を取ってくれるため、差分系列は使わない。

データの準備

Seatbelts_log <- Seatbelts[,c("front", "PetrolPrice", "law")]
Seatbelts_log[,"front"] <- log(Seatbelts[,"front"])
Seatbelts_log[,"PetrolPrice"] <- log(Seatbelts[,"PetrolPrice"])
head(Seatbelts_log)

##             front PetrolPrice law
## Jan 1969 6.765039    -2.27330   0
## Feb 1969 6.715383    -2.27923   0
## Mar 1969 6.692084    -2.28217   0
## Apr 1969 6.701960    -2.29389   0
## May 1969 6.898715    -2.29244   0
## Jun 1969 6.851185    -2.29679   0

訓練データとテストデータに分ける

train <- window(Seatbelts_log, end=c(1983,12))
test <- window(Seatbelts_log, start=c(1984,1))

説明変数だけ切り出す

petro_law <- train[, c("PetrolPrice", "law")]

ARIMAモデルによる推定

forecastパッケージのArima関数を用いた
yは応答変数としてfrontを指定
orderはSARIMA(p,d,q)(P,D,Q)の次数を設定
seasonalは季節成分の次数(P,D,Q)を設定
xregは説明変数を指定

model_sarimax <- Arima(
  y = train[, "front"], 
  order = c(1, 1, 1),
  seasonal = list(order = c(1, 0, 0)), 
  xreg = petro_law
)
model_sarimax

## Series: train[, "front"] 
## Regression with ARIMA(1,1,1)(1,0,0)[12] errors 
## 
## Coefficients:
##          ar1      ma1    sar1  PetrolPrice      law
##       0.2589  -0.9503  0.6877      -0.3464  -0.3719
## s.e.  0.0826   0.0303  0.0548       0.0955   0.0467
## 
## sigma^2 = 0.009052:  log likelihood = 165.33
## AIC=-318.66   AICc=-318.18   BIC=-299.54

計算結果はCoefficientsを確認する。s.e.は標準誤差。
石油価格(PetroPrice)が上がると交通事故死傷者数は減ることがわかった。
法律(law)が施行されると交通事故死傷者数は減ることがわかった。
AICなどを見てモデル選択をする。

次に、最適なモデルの次数を選ぶ作業に移る。

自動モデル選択auto.arima関数

モデルの次数を決定する作業をする
モデルの次数をひたすら変えたAICの計算を自動モデル選択関数であるforecastパッケージのauto.arima関数を使って行う
ic = “aic”はAICを使ってモデル選択をするという指定、他にaiccやbicがある
max.order = 7はSARIMA(p,d,q)(P,D,Q)におけるp+q+P+Qの最大値。これを増やすほど複雑なモデルを候補に入れてモデル選択することができる
stepwise = F,approximation = Fは共に計算量をケチらないという指定
stepwise=Tとすると候補となる次数の組み合わせが減る
approximation=Tの場合は毎回の計算において近似的な手法を使うことで計算速度を上げる。特にapproximatio=Tになっていると誤った結果が出やすくなる
定常性・反転可能性のチェックはauto.arima関数の中で行われている

sarimax_petro_law <- auto.arima(
  y = train[, "front"], 
  xreg = petro_law,
  ic = "aic",
  max.order = 7,
  stepwise = F,
  approximation = F,
  parallel = T,
  num.cores = 4
)
sarimax_petro_law

## Series: train[, "front"] 
## Regression with ARIMA(2,0,1)(0,1,1)[12] errors 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1  PetrolPrice      law
##       1.1225  -0.1322  -0.8690  -0.8183      -0.3748  -0.3431
## s.e.  0.0906   0.0876   0.0443   0.1129       0.1000   0.0473
## 
## sigma^2 = 0.007624:  log likelihood = 168.12
## AIC=-322.23   AICc=-321.53   BIC=-300.36

かなり次数の大きいモデルが選ばれた。

残差のチェック

残差の自己相関のチェック（検定）

帰無仮説は「残差に自己相関がない」

checkresiduals(sarimax_petro_law)

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(2,0,1)(0,1,1)[12] errors
## Q* = 20.99, df = 18, p-value = 0.2799
## 
## Model df: 6.   Total lags used: 24

p-value>0.05であるため有意な自己相関は見られなかった
これは明らかな問題がないことの確認であって、良いモデルであることの保証ではない

上のグラフが残差系列
左下が残差のコレログラム
右下が残差のヒストグラム
異常に突出した残差などの問題がないことを確認した

残差の正規性の検定

帰無仮説は「正規分布に従う」
tseriesパッケージのjarque.bera.testを用いた
残差の計算はresid関数を用いた

jarque.bera.test(resid(sarimax_petro_law))

## 
##  Jarque Bera Test
## 
## data:  resid(sarimax_petro_law)
## X-squared = 0.39938, df = 2, p-value = 0.819

こちらも正規分布と有意に異なっているとは言えないという結果となった
これでモデルの同定を終わる

ARIMAによる予測

forecastパッケージのforecast関数を用いて予測した。
h=12は12時点先まで予測。 level = c(95, 70)は95%,70%予測区間も併せて出力するという指定。

petro_law_test <- test[, c("PetrolPrice", "law")]
sarimax_f <- forecast(
  sarimax_petro_law, 
  xreg = petro_law_test,
  h = 12,
  level = c(95, 70)
)
sarimax_f

##          Point Forecast    Lo 70    Hi 70    Lo 95    Hi 95
## Jan 1984       6.140922 6.050389 6.231454 5.969719 6.312124
## Feb 1984       6.079469 5.986077 6.172861 5.902858 6.256080
## Mar 1984       6.167368 6.072965 6.261771 5.988846 6.345890
## Apr 1984       6.196306 6.101088 6.291523 6.016243 6.376368
## May 1984       6.293745 6.197758 6.389732 6.112227 6.475262
## Jun 1984       6.274492 6.177762 6.371222 6.091568 6.457416
## Jul 1984       6.394263 6.296812 6.491714 6.209976 6.578550
## Aug 1984       6.427685 6.329534 6.525836 6.242075 6.613294
## Sep 1984       6.354929 6.256100 6.453759 6.168036 6.541823
## Oct 1984       6.396352 6.296863 6.495840 6.208212 6.584491
## Nov 1984       6.340966 6.240838 6.441095 6.151616 6.530317
## Dec 1984       6.448938 6.348188 6.549689 6.258413 6.639464

autoplot(sarimax_f, predict.colour=1, main = "ARIMAによる予測")

この予測には問題がある。将来の石油価格が分かっている前提で予測を出している。
未来の価格は分からないので、何らかの値で代用する必要がある。

過去の石油価格の平均値を使う

rep(1,12)は1が12回連続で出てくる。

petro_law_mean <- data.frame(
  PetrolPrice=rep(mean(train[, "PetrolPrice"]),12),
  law=rep(1, 12)
)
sarimax_f_mean <- forecast(sarimax_petro_law, xreg=as.matrix(petro_law_mean))
sarimax_f_mean

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1984       6.194413 6.082470 6.306357 6.023211 6.365616
## Feb 1984       6.123359 6.007879 6.238839 5.946748 6.299970
## Mar 1984       6.214309 6.097580 6.331038 6.035787 6.392831
## Apr 1984       6.242017 6.124280 6.359754 6.061955 6.422080
## May 1984       6.337695 6.219007 6.456383 6.156177 6.519212
## Jun 1984       6.318318 6.198711 6.437926 6.135395 6.501242
## Jul 1984       6.438607 6.318108 6.559105 6.254320 6.622894
## Aug 1984       6.471575 6.350211 6.592939 6.285965 6.657184
## Sep 1984       6.396515 6.274312 6.518718 6.209621 6.583408
## Oct 1984       6.445650 6.322632 6.568668 6.257510 6.633789
## Nov 1984       6.388848 6.265039 6.512657 6.199498 6.578198
## Dec 1984       6.496951 6.372373 6.621529 6.306425 6.687477

autoplot(sarimax_f_mean, predict.colour=1, main = "ARIMAによる予測")

直近の石油価格を使う

また、直近の石油価格を使った予測値も出しておく。
tail関数で末尾の値を抽出できる。

petro_law_tail <- data.frame(
  PetrolPrice=rep(tail(train[, "PetrolPrice"], n=1),12),
  law=rep(1, 12)
)
sarimax_f_tail <- forecast(sarimax_petro_law, xreg=as.matrix(petro_law_tail))
sarimax_f_tail

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1984       6.141143 6.029200 6.253086 5.969940 6.312345
## Feb 1984       6.070088 5.954609 6.185568 5.893477 6.246699
## Mar 1984       6.161038 6.044309 6.277768 5.982516 6.339560
## Apr 1984       6.188746 6.071010 6.306483 6.008684 6.368809
## May 1984       6.284424 6.165736 6.403112 6.102906 6.465942
## Jun 1984       6.265048 6.145440 6.384655 6.082124 6.447971
## Jul 1984       6.385336 6.264837 6.505835 6.201049 6.569623
## Aug 1984       6.418304 6.296941 6.539668 6.232695 6.603914
## Sep 1984       6.343244 6.221041 6.465447 6.156350 6.530137
## Oct 1984       6.392379 6.269361 6.515397 6.204239 6.580519
## Nov 1984       6.335577 6.211768 6.459387 6.146227 6.524927
## Dec 1984       6.443680 6.319102 6.568258 6.253154 6.634206

autoplot(sarimax_f_tail, predict.colour=1, main = "ARIMAによる予測")

ナイーブ予測

過去の平均値を予測値として使う

mean(train[, “front”])と同じ。

naive_f_mean <- meanf(train[, "front"], h = 12)
naive_f_mean

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Feb 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Mar 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Apr 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## May 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Jun 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Jul 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Aug 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Sep 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Oct 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Nov 1984       6.729989 6.468432 6.991546 6.328736 7.131243
## Dec 1984       6.729989 6.468432 6.991546 6.328736 7.131243

過去の最新の値を、予測値として使う

naive_f_latest <- rwf(train[, "front"], h = 12)
naive_f_latest

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1984       6.371612 6.183035 6.560189 6.083208 6.660015
## Feb 1984       6.371612 6.104924 6.638300 5.963748 6.779476
## Mar 1984       6.371612 6.044987 6.698237 5.872082 6.871141
## Apr 1984       6.371612 5.994458 6.748766 5.794805 6.948419
## May 1984       6.371612 5.949941 6.793283 5.726722 7.016502
## Jun 1984       6.371612 5.909695 6.833529 5.665170 7.078053
## Jul 1984       6.371612 5.872684 6.870539 5.608568 7.134656
## Aug 1984       6.371612 5.838236 6.904988 5.555884 7.187340
## Sep 1984       6.371612 5.805881 6.937343 5.506401 7.236822
## Oct 1984       6.371612 5.775279 6.967944 5.459600 7.283624
## Nov 1984       6.371612 5.746173 6.997051 5.415086 7.328138
## Dec 1984       6.371612 5.718362 7.024861 5.372553 7.370671

tail(train[, "front"], n = 1)

##           Dec
## 1983 6.371612

予測の評価

訓練データとテストデータの予測のRMSEを調べた（未来の石油価格がわかっている前提の予測）

accuracy(sarimax_f, x=test[, "front"])

##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set -0.004353509 0.08283297 0.06430968 -0.07784816 0.9593004 0.5848660
## Test set      0.071495376 0.09674572 0.07537019  1.10793155 1.1703728 0.6854562
##                    ACF1 Theil's U
## Training set -0.0040321        NA
## Test set      0.3376395   1.03686

テストデータの予測のRMSEが0.0967であることがわかった。
訓練データの当てはめ精度よりも予測精度が悪くなりすぎていたら問題。

sarimax①のRMSE 石油価格の平均値を使用した予測のテストデータ

accuracy(sarimax_f_mean, x=test[, "front"])["Test set", "RMSE"]

## [1] 0.06945114

sarimax②のRMSE 直近の石油価格を使用した予測のテストデータ

accuracy(sarimax_f_tail, x=test[, "front"])["Test set", "RMSE"]

## [1] 0.1018344

sarimax①の予測精度のほうが高いことがわかった。

ナイーブ予測①のRMSE 石油価格の過去の平均値を使用したナイーブ予測

accuracy(naive_f_mean, x=test[, "front"])["Test set", "RMSE"]

## [1] 0.3949872

sarimax①と比較して値が上回っているため、予測精度がsarimaxの方が良いことがわかった。

ナイーブ予測①のRMSE 直近の石油価格を使用したナイーブ予測

accuracy(naive_f_latest, x=test[, "front"])["Test set", "RMSE"]

## [1] 0.1498196

sarimax②と比較して値が上回っているため、予測精度がsarimaxの方が良いことがわかった。

結論

・短期の自己相関や季節性など様々な要因がある複雑なデータでもSARIMAXを使うことでモデル化することができた。
・単純に平均値で求めた予測よりもSARIMAXモデルで求めた予測のほうが高い精度が出ることがわかった。
・石油価格の直近の値を使用したSARIMAXモデルよりも過去の平均値を使用したほうが良い精度で分析できることがわかった。