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

#はじめに
ARモデルの復習として、自身の実施の記録をQiitaに記録します。基本的には[こちらの内容](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)に沿って実施しているだけですが、一応、Nileデータについても追加でやってみました。
なお、上記のリンク先は、「[Rによるデータサイエンス](https://www.morikita.co.jp/books/book/3063)」で書籍化もされています。一般的な統計手法が網羅されており、辞書的にも使えるのでとても良い書籍と思います。

#ARモデル
時系列地点t-pからtまでの各データの関係式

```math
y_t = \sum_{i=1}^{p}\phi_iy_{t-i}+e_t
```
を自己回帰(AutoRegression AR)モデルと呼びます。ここで、φは自己回帰係数と呼びます。また、pを次数と呼びます。次数pのARモデルはAR(p)を表現します。また、eはホワイトノイズを示します。
現在のデータが過去のデータの線形和と誤差項で考えるもので、各過去データとの相関の強さをφで表現し、過去いつまでの影響を検討するかというものをpで示す、というものです。
ARモデルの推定とは、上記を仮定した上で、次数pと次数分の各φの値をデータより推定する問題を解くことになります。

#ARモデルの推定
RではARモデルのパラメータを推定する関数として、ar関数が存在します。
ar関数では推定アルゴリズムはmethodで指定できます。指定可能なアルゴリズムは４つ存在し、ユール-ウォーカー法、最小２乗法、最尤法、バーグ法があります。デフォルトはユール-ウォーカー法が選択されます。ユール-ウォーカー法の数理的な話は[こちら](https://math-eng.com/autoregressive-yule-walker/)が参考になります。

ここで、標準データのlh(Luteinizing Hormone in Blood Sample)とNile(Flow of the River Nile)の各データをARモデルに当てはめてみて各パラメータを推定してみます。

まず、plot関数でデータに定常性(平均や分散が時間変化せず時間的に変化しない一定の確率過程)がありそうかを確認します。（今回は省略）
次に、自己相関係数と偏自己相関係数(注目している時以外を無視した自己相関係数)をそれぞれacfとpacfで確認してみます。

```R
> acf(lh)
```

![plot (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/dd57a711-135a-fce6-ece3-727f3e0160fe.png)

```R
> pacf(lh)
```

![plot (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/9b68de60-8dc2-25bd-6c32-3daf523842ba.png)

```R
> acf(Nile)
```

![plot (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/ccc3a9b1-939c-929d-1c28-6f81a228a488.png)

```R
> pacf(Nile)
```

![plot (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/84e9d8c8-ce09-b33b-20aa-a5d9efc15348.png)

lhは自己相関としては１、２期前に相関が見られますが、偏自己相関としては１期前のみのようです。また、Nileは１〜４期前に相関が見られますが、偏自己相関としては１期前のみのようです。
ここでそれぞれのデータについてARモデルを仮定し、次元と係数を推定してみます。ここでは、推定アルゴリズムはデフォルトのユール・ウォーカー法を使ってみました。

```R
> lh.ar <- ar(lh)
> lh.ar

Call:
ar(x = lh)

Coefficients:
      1        2        3  
 0.6534  -0.0636  -0.2269  

Order selected 3  sigma^2 estimated as  0.1959

> Nile.ar <- ar(Nile)
> Nile.ar

Call:
ar(x = Nile)

Coefficients:
     1       2  
0.4081  0.1812  

Order selected 2  sigma^2 estimated as  21247
```

係数や次数は次のように取得することができます。

```R
> Nile.ar$order
[1] 2
> Nile.ar$ar
[1] 0.4081111 0.1811710
```

この結果、lhデータでは、次のAR(3)モデルが推定されたことになります。

```math
\hat{y_t} = 0.653y_{t-1}-0.063y_{t-2}-0.227y_{t-3}
```

また、Nileデータでは、次のAR(2)モデルが推定されたことになります。

```math
\hat{y_t} = 0.408y_{t-1}-0.181y_{t-2}
```

推定されたモデルはacfやpacfで見られた自己相関係数の状況に比べても、そこまで離れた結果ではなさそうに思います。

#ARモデルでの予測
predict関数を使って、10期先まで予測してみます。
まず、lhデータを予測します

```R
> lh.pr<-predict(lh.ar,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 (9).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/eac9ea82-daff-2228-b3cb-51b2179e4c9f.png)


```R
> Nile.pr<-predict(Nile.ar,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 (10).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/187459/e4dbb4de-cbfa-a903-4bab-03e5b88b4e4e.png)

#最後に
先に書いた通り、[Rと時系列(2)](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)に習って試し、追加でNileデータにもARモデルを適用してみました。モデルの診断など足りない箇所については後日、理解の範囲で記載できればと思います。

#参考文献
[Rと時系列(2)](https://www1.doshisha.ac.jp/~mjin/R/Chap_34/34.html)
[ARモデルとユール–ウォーカー方程式](https://math-eng.com/autoregressive-yule-walker/)