R機械学習/SARIMAXモデルで時系列予測

# Introduction
時系列モデリング手法の1つであるSARIMAXをRで実装する方法についてのメモ．

# パッケージインストールとインポート
```r
library(xts)  # xts型の提供
library(forecast)  # 時系列モデルの作成と予測
library(ggfortify)  # グラフ描画
library(tidyverse)  # データ操作
library(tseries)  # 正規性の検定
```

# データ操作
今回は，Rにデフォルトで用意されているSeatbeltsデータを使用．

- データ説明：英国における1969年～1984年の自動車運転手の死傷者数．期間中の1983年1月31日にシートベルトの着用が義務づけられた．
- 変数説明
  - DriversKilled：自動車運転手の死亡者数
  - drivers：死傷者数
  - front：前座席の同乗者の死傷者数
  - rear：後部座席の同乗者の死傷者数
  - kms：走行距離
  - PetrolPrice：ガソリン価格
  - VanKilled：バン（軽貨物車）運転者の死亡数
  - law：シートベルトの着用義務

```r
dataset = as.xts(Seatbelts)
dim(dataset)
head(dataset)
```
Out:

```
[1] 192   8
       DriversKilled drivers front rear   kms PetrolPrice VanKilled law
1 1969           107    1687   867  269  9059   0.1029718        12   0
2 1969            97    1508   825  265  7685   0.1023630         6   0
3 1969           102    1507   806  319  9963   0.1020625        12   0
4 1969            87    1385   814  407 10955   0.1008733         8   0
5 1969           119    1632   991  454 11823   0.1010197        10   0
6 1969           106    1511   945  427 12391   0.1005812        13   0
```

## データ分割
今回はdriversを自身の過去の値から予測する．説明変数（外生変数）としてkmsとPetroPriceを加える．
シートベルト着用義務の法律が施行される前の期間で訓練データとテストデータを作成する．

```r
dataset_law0 = dataset['::1982-12', c('drivers', 'kms','PetrolPrice')]
train = dataset_law0['::1979-12']
test = dataset_law0['1980-01::']

y_train = train[,'drivers']
X_train = train[,c('kms','PetrolPrice')]
y_test = test[,'drivers']
X_test = test[,c('kms','PetrolPrice')]
```

## データの特徴把握
まずは全変数の時系列グラフを確認

```r
autoplot(train)
```
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/f26cbbd4-2432-4b6a-6c9e-047db60e8200.png)

続いて時系列データ特有の特徴を確認

```r
ggtsdisplay(y_train, main = 'drivers')
```

- 上：データの時系列グラフ
- 左下：自己相関コレログラム；Lag=12で相関が強くなっているため，12ヶ月の周期性があると見て取る．
- 右下：偏自己相関コレログラム；Lag=1で相関が強くなっているため，1時刻前の値が強く影響していると見て取る．
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/c133332a-c1ae-4f8a-557b-f4747171cff4.png)

先ほどとは別の方法でデータの特徴を把握する．
STL分解をすることでデータを以下のように3つの要素に分解する． 
$$
元データ=季節性+トレンド+残差
$$

- seasonalグラフ：似た波形が繰り返し出現しているため，季節性ありと判断できる．
- trendグラフ：一部区間で上昇トレンドが見て取れる．
- remainderグラフ：値が極端に高いところなどは外れ値となるが，今回の例では外れ値は無いと判断する．

```r
plot(stl(x=y_train, s.window='periodic'), main='driver')
```
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/a06c4cdb-ff5d-f9a2-12f0-6fe38814b7ec.png)

# モデル作成
SARIMAXモデルは以下のような要素から成るモデル．
$$
SARIMAX = Seasonal(季節成分)+AR(自己回帰)モデル+Integrated(和分)+MA(移動平均)モデル+eXgenous(外因性)
$$

SARIMAXモデルを構築する．
`auto.arima()`でパラメータを最適化したモデルを自動選択できる．
選択されたパラメータはARIMA(p,d,q)(P,D,Q)[s]の形で表示される．

- p,(P)：ARの次数（季節成分に対するARの次数）
- d,(D)：和分の階数（季節成分に対する和分の階数）
- q,(Q)：MAの次数（季節成分に対するMAの次数）
- s：季節成分の周期

<details><summary>一部引数の詳細</summary><div>

```r
auto.arima(
  y,  # 学習する時系列データ
  max.p = 5,  # 探索する次数p(ARモデル)の最大値
  max.q = 5,  # 探索する次数q(MAモデル)の最大値
  max.P = 2,  # 探索する次数P(季節成分のARモデル)の最大値
  max.Q = 2,  # 探索する次数Q(季節成分のMAモデル)の最大値
  max.order = 5,  #(p+q+P+Q)の合計値．
  seasonal = TRUE,  # TRUE:SARIMAモデル; FALSE:ARIMAモデル
  ic = c("aicc", "aic", "bic"),  # モデル選択に使用される情報量基準
  stepwise = TRUE,  # True:ステップワイズで探索; False:全探索
  trace = FALSE,  # Trueの場合探索したモデルの情報がリスト表示される
  approximation = (length(x) > 150 | frequency(x) > 12),  # 指定した条件に該当する場合近似解を用いる．近似解を使用しない場合はFALSEにする．
  xreg = NULL,  # 外生変数を指定
  parallel = FALSE,  # TRUE:並列処理を行う
  num.cores = 2,  # 並列処理を行う際のコア数の指定
  ...
)
```
</div></details>

```r
model_sarimax = auto.arima(y_train, max.order=14, ic='aic', ,stepwise=F, parallel=T, num.cores=4, xreg=X_train)
model_sarimax
```
Out:

```
Series: y_train 
Regression with ARIMA(3,0,1)(0,1,1)[12] errors 

Coefficients:
          ar1     ar2     ar3     ma1     sma1      kms  PetrolPrice
      -0.5666  0.5975  0.4085  0.9443  -0.8597  -0.0027    -7755.019
s.e.   0.0988  0.0905  0.0885  0.0832   0.1723   0.0197     2259.103

sigma^2 estimated as 16913:  log likelihood=-759.8
AIC=1535.6   AICc=1536.9   BIC=1557.9
```

ちなみに上記のモデルを`Arima()`を使って構築すると以下のようになる．

```r
model_sarimax2 = arima(y_train, order=c(3,0,1), seasonal=c(0,1,1), xreg=X_train)
model_sarimax2
```
Out:

```
Call:
arima(x = y_train, order = c(3, 0, 1), seasonal = c(0, 1, 1), xreg = X_train)

Coefficients:
          ar1     ar2     ar3     ma1     sma1      kms  PetrolPrice
      -0.5666  0.5975  0.4085  0.9443  -0.8597  -0.0027    -7755.019
s.e.   0.0988  0.0905  0.0885  0.0832   0.1723   0.0197     2259.103

sigma^2 estimated as 15926:  log likelihood = -759.8,  aic = 1535.6
```


# モデル評価
## 残差チェック
`checkresiduals`で残差に自己相関があるかのチェックを行う．
帰無仮説は「自己相関がない」なので，今回の例では，p-value=0.09であるため，有意水準5%では自己相関があるとは言えないことになる．よってチェックはクリア．
同時に出力されるグラフに関する説明は以下の通り．

- 上：残差のの時系列グラフ；周期性が見られるとそのモデルは信頼性に欠ける
- 左下：残差の自己相関コレログラム；自己相関があるとそのモデルは信頼性に欠ける．
- 右下：残差のヒストグラム；正規分布から逸脱しているとそのモデルは信頼性に欠ける．

```r
checkresiduals(model_sarimax)
```
Out:

```
	Ljung-Box test

data:  Residuals from Regression with ARIMA(3,0,1)(0,1,1)[12] errors
Q* = 25.124, df = 17, p-value = 0.09198

Model df: 7.   Total lags used: 24
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/c6e0d482-f31c-4da4-97d6-5efd122cd852.png)

`jarque.bera.test()`で残差の正規性の検定を行う．  
帰無仮説は「正規分布に従う」なので，今回の例ではp-value=0.6125であるため，有意水準5%では正規分布に従わないとは言えないことになる．よってチェックはクリア．

```r
jarque.bera.test(resid(model_sarimax))
```
Out:

```
	Jarque Bera Test

data:  resid(model_sarimax)
X-squared = 0.98046, df = 2, p-value = 0.6125
```

## 予測
今回は未来のの外生変数の値が分かっている前提で予測する．  

```r
pred_sarimax = forecast(model_sarimax, xreg=X_test, h=length(y_test))
autoplot(pred_sarimax, y_test)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/df1a5068-e891-d9d9-cf4e-e66f940e7f87.png)

精度評価指標の算出

```r
accuracy(pred_sarimax, y_test)
```
Out:

```
                     ME     RMSE       MAE        MPE     MAPE      MASE         ACF1
Training set  -3.552711 120.3265  95.78231 -0.6330721 5.553774 0.6443118 -0.001547839
Test set     -17.228042 147.2290 110.48983 -1.3294291 6.744713 0.7432468           NA
```

訓練データの実測値と予測値を可視化

```r
test_meas_pred = data.frame(time=index(y_train), measure=as.vector(y_train), predict=model_sarimax$fitted)
test_meas_pred = gather(test_meas_pred, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'red'))
```
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/31607e9c-ee4e-c4fe-efd9-dacc458270e5.png)

テストデータの実測値と予測値を可視化

```r
test_meas_pred = data.frame(time=index(y_test), measure=as.vector(y_test), predict=pred_sarimax$mean)
test_meas_pred = gather(test_meas_pred, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'red'))
```
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/ce10dab9-55e3-e3b2-c0b6-e26341f61fcb.png)

# モデル比較
モデル単体では精度の良し悪しわ分からないので，今回は以下の3つのモデルで精度を比較する．

- SARIMAXモデル：先ほど作成したもの．
- SARIMAモデル：外生変数を含まないもの．
- ナイーブ予測：モデルを使わないシンプルな予測．今回は1周期前の値をそのまま予測値とする．

SARIMAモデル

```r
model_sarima = auto.arima(y_train, max.order=10, ic='aic', ,stepwise=F, parallel=T, num.cores=4)
model_sarima
```
Out:

```
Series: y_train 
ARIMA(1,1,4)(0,1,1)[12] 

Coefficients:
          ar1     ma1      ma2      ma3      ma4     sma1
      -0.7517  0.2348  -0.4050  -0.0433  -0.3153  -0.8690
s.e.   0.1194  0.1455   0.1136   0.0932   0.1028   0.1932

sigma^2 estimated as 17777:  log likelihood=-757.71
AIC=1529.43   AICc=1530.43   BIC=1548.88
```
```r
pred_sarima = forecast(model_sarima, h=length(y_test))
```

ナイーブ予測

```r
pred_naive = snaive(y_train, h=length(y_test))
```

テストデータで実測値と各モデルの予測値を比較

```r
test_meas_pred3 = data.frame(time=index(y_test), measure=as.vector(y_test), SARIMAX=pred_sarimax$mean, SARIMA=pred_sarima$mean, naive=pred_naive$mean)
test_meas_pred3 = gather(test_meas_pred3, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred3, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'blue', 'green', 'red'))
```
Out:
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/931058/7dd59d6c-0d03-74f0-7e71-e85ce04508c3.png)

RMSEを比較したところ，テストデータにおいてはSARIMAモデルが最良の結果となった．

```r
rmse_df = rbind(
  t(as.data.frame(accuracy(pred_sarimax, y_test)[,'RMSE'])),
  t(as.data.frame(accuracy(pred_sarima, y_test)[,'RMSE'])),
  t(as.data.frame(accuracy(pred_naive, y_test)[,'RMSE']))
)
rownames(rmse_df) = c('SARIMAX', 'SARIMA', 'naive')
rmse_df
```
Out:

```
        Training set Test set
SARIMAX     120.3265 147.2290
SARIMA      123.3607 133.0183
naive       196.8749 182.3017
```

# Conclusion
今回はSARIMAXモデルの構築から精度の評価まで実装してみた．
時系列モデリング手法はVARや状態空間モデル，ディープラーニング等まだまだ気になる手法があるので，それらも随時投稿していこうと思う．

------------------------------------------
------------------------------------------
# code
```py
# パッケージインストールとインポート
library(xts)  # xts型の提供
library(forecast)  # 時系列モデルの作成と予測
library(ggfortify)  # グラフ描画
library(tidyverse)  # データ操作
library(tseries)  # 正規性の検定


# データ操作
dataset = as.xts(Seatbelts)
dim(dataset)
head(dataset)

## データ分割
dataset_law0 = dataset['::1982-12', c('drivers', 'kms','PetrolPrice')]
train = dataset_law0['::1979-12']
test = dataset_law0['1980-01::']

y_train = train[,'drivers']
X_train = train[,c('kms','PetrolPrice')]
y_test = test[,'drivers']
X_test = test[,c('kms','PetrolPrice')]

## データの特徴把握
autoplot(train)

ggtsdisplay(y_train, main = 'drivers')

plot(stl(x=y_train, s.window='periodic'), main='driver')


# モデル作成
model_sarimax = auto.arima(y_train, max.order=14, ic='aic', ,stepwise=F, parallel=T, num.cores=4, xreg=X_train)
model_sarimax

model_sarimax2 = arima(y_train, order=c(3,0,1), seasonal=c(0,1,1), xreg=X_train)
model_sarimax2

# モデル評価
## 残差チェック
checkresiduals(model_sarimax)

jarque.bera.test(resid(model_sarimax))

## 予測
pred_sarimax = forecast(model_sarimax, xreg=X_test, h=length(y_test))

autoplot(pred_sarimax, y_test)

accuracy(pred_sarimax, y_test)

test_meas_pred = data.frame(time=index(y_train), measure=as.vector(y_train), predict=model_sarimax$fitted)
test_meas_pred = gather(test_meas_pred, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'red'))

test_meas_pred = data.frame(time=index(y_test), measure=as.vector(y_test), predict=pred_sarimax$mean)
test_meas_pred = gather(test_meas_pred, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'red'))


# モデル比較
model_sarima = auto.arima(y_train, max.order=10, ic='aic', ,stepwise=F, parallel=T, num.cores=4)
model_sarima

pred_sarima = forecast(model_sarima, h=length(y_test))

pred_naive = snaive(y_train, h=length(y_test))

test_meas_pred3 = data.frame(time=index(y_test), measure=as.vector(y_test), SARIMAX=pred_sarimax$mean, SARIMA=pred_sarima$mean, naive=pred_naive$mean)
test_meas_pred3 = gather(test_meas_pred3, key='method', value='value', -time)

ggplot() +
  geom_line(data=test_meas_pred3, aes(x=time, y=value, colour=method)) +
  scale_color_manual(values=c('black', 'blue', 'green', 'red'))

rmse_df = rbind(
  t(as.data.frame(accuracy(pred_sarimax, y_test)[,'RMSE'])),
  t(as.data.frame(accuracy(pred_sarima, y_test)[,'RMSE'])),
  t(as.data.frame(accuracy(pred_naive, y_test)[,'RMSE']))
)
rownames(rmse_df) = c('SARIMAX', 'SARIMA', 'naive')
rmse_df

```