ARIMAモデルをRで試してみる

#はじめに
ARIMAモデルの復習として、以前の[ARモデルをRで試してみる](https://qiita.com/purple_jp/items/b5b36bc78a4ccfb42027)の記事同様に記録を残します。
前回同様、基本的には[Rと時系列(2)](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)に沿って実施しているだけですが、
Nileデータについても追加でやっています。
また、パラメータ推定については実装はせずに[forecastパッケージ](https://cran.r-project.org/web/packages/forecast/forecast.pdf)のauto.arimaを利用しました。

#ARMAモデルとARIMAモデルについて
ARはAutoRegressiveのことで自己回帰のことでした。
これに誤差の移動平均（Moving Average）$\sum_{j=1}^{q}{b_je_{t-j}}$を加えたものが、ARMAモデルです。
数式にすると次です。

```math
y_t = \sum_{i=1}^{p}{a_iy_{i-1}} + e_t + \sum_{j=1}^{q}{b_je_{t-j}}
```

さらにy_tのd階の差分$\Delta^dy_d$のARMAモデルを自己回帰和分移動平均
(ARIMA:AutoRegressive Integrated Moving Average)モデルと言います。
ARIMA($p$,$d$,$q$)で表し、式では次のようになります。

```math
y_t - y_{t-d} = \sum_{i=1}^{p}{a_iy_{i-1}} + e_t + \sum_{j=1}^{q}{b_je_{t-j}}
```

上式を$y_t$で整理すると、$y_t$の項に$y_{t-d}$に足される形なので、和分となっていると思います。
なお、差分で表現するのはトレンド除去の対応など時系列データの整形でよくあることなので、
差分でデータ整形しARMAモデル適用=ARIMAモデルという理解で多分良いのかなと思います。

#ARIMAモデルのパラメータ推定方法
ARモデルの時はユールウォーカー法で必要なパラメータを予測しましたが、
今回のARIMA($p$,$d$,$q$)は、各パラメータのすべての組み合わせの中から
AICなどの情報量基準を使い、最も小さい組み合わせを選ぶようです。
参考としている[Rと時系列(2)](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)では、別途実装しパラメータ推定を行っていますが、
[forecastパッケージ](https://cran.r-project.org/web/packages/forecast/forecast.pdf)のauto.arimaを使えば、次のように計算を簡単に実行できました。

```R
> library(forecast)
> auto.arima(lh, ic="aic", stepwise=T, trace=T) 

 ARIMA(2,0,2) with non-zero mean : 66.42642
 ARIMA(0,0,0) with non-zero mean : 82.09291
 ARIMA(1,0,0) with non-zero mean : 64.75832
 ARIMA(0,0,1) with non-zero mean : 68.10389
 ARIMA(0,0,0) with zero mean     : 224.6837
 ARIMA(2,0,0) with non-zero mean : 64.50375
 ARIMA(3,0,0) with non-zero mean : 64.18482
 ARIMA(4,0,0) with non-zero mean : 65.84091
 ARIMA(3,0,1) with non-zero mean : 64.47047
 ARIMA(2,0,1) with non-zero mean : 65.20321
 ARIMA(4,0,1) with non-zero mean : 66.41819
 ARIMA(3,0,0) with zero mean     : 80.68744

 Best model: ARIMA(3,0,0) with non-zero mean 

Series: lh 
ARIMA(3,0,0) with non-zero mean 

Coefficients:
         ar1      ar2      ar3    mean
      0.6448  -0.0634  -0.2198  2.3931
s.e.  0.1394   0.1668   0.1421  0.0963

sigma^2 estimated as 0.1949:  log likelihood=-27.09
AIC=64.18   AICc=65.61   BIC=73.54 
```

同じようにNileデータもやってみました。

```R
> auto.arima(Nile, ic="aic", stepwise=T, trace=T) 

 ARIMA(2,1,2) with drift         : Inf
 ARIMA(0,1,0) with drift         : 1298.645
 ARIMA(1,1,0) with drift         : 1283.346
 ARIMA(0,1,1) with drift         : 1270.309
 ARIMA(0,1,0)                    : 1296.697
 ARIMA(1,1,1) with drift         : 1267.637
 ARIMA(2,1,1) with drift         : 1269.132
 ARIMA(1,1,2) with drift         : 1269.134
 ARIMA(0,1,2) with drift         : 1268.544
 ARIMA(2,1,0) with drift         : 1279.282
 ARIMA(1,1,1)                    : 1267.255
 ARIMA(0,1,1)                    : 1269.091
 ARIMA(1,1,0)                    : 1281.48
 ARIMA(2,1,1)                    : 1268.896
 ARIMA(1,1,2)                    : 1268.923
 ARIMA(0,1,2)                    : 1267.957
 ARIMA(2,1,0)                    : 1277.484
 ARIMA(2,1,2)                    : Inf

 Best model: ARIMA(1,1,1)                    

Series: Nile 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.2544  -0.8741
s.e.  0.1194   0.0605

sigma^2 estimated as 20177:  log likelihood=-630.63
AIC=1267.25   AICc=1267.51   BIC=1275.04
```

ただし、[Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測](https://tjo.hatenablog.com/entry/2013/07/12/184704)によると、
このパラメータ推定ですが探索的に行なっていることから、
真のパラメータ以外のところをベストモデルと判定してしまうことがあるようです。
そのため、ある程度あたりがついている場合、探索範囲を狭めるのが良いようです。

#推定されたパラメータを元に予測してみる
せっかくauto.arimaでパラメータ推定したので予測をしたいと思います。
推定したパラメータでarimaモデルを作ってから、predict関数で予測します。前回と同じです。
ということで、まずはlhデータを。

```R
> lh.arima<-arima(lh,order=c(3,0,0))
> lh.pr<-predict(lh.arima,n.ahead=10)
> lh.SE1<-lh.pr$pred+2*lh.pr$se
> lh.SE2<-lh.pr$pred-2*lh.pr$se
> ts.plot(lh,lh.pr$pred,lh.SE1,lh.SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
```

![plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/cbc97f61-f242-38f6-4917-2830e079b310.png)

次に、Nileデータを。

```R
> Nile.arima<-arima(Nile,order=c(1,1,1))
> Nile.pr<-predict(Nile.arima,n.ahead=10)
> Nile.SE1<-Nile.pr$pred+2*Nile.pr$se
> Nile.SE2<-Nile.pr$pred-2*Nile.pr$se
> ts.plot(Nile,Nile.pr$pred,Nile.SE1,Nile.SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
```

![plot (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/38609bef-65ab-4898-307a-5903d21ccb31.png)

ということで、簡単ですがARIMAモデルのパラメータ推定と予測を行ってみました。
なお、ARIMAモデルの診断についてはtsdiagを使うそうですが、
これについては後日この記事に追記するか、別の記事として改めて取り上げたいと思います。

#参考文献
- [Rと時系列(2)](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)
- [ARモデルをRで試してみる](https://qiita.com/purple_jp/items/b5b36bc78a4ccfb42027)
- [forecastパッケージ](https://cran.r-project.org/web/packages/forecast/forecast.pdf)
- [Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測](https://tjo.hatenablog.com/entry/2013/07/12/184704)