はじめに

ARモデルの復習として、自身の実施の記録をQiitaに記録します。基本的にはこちらの内容に沿って実施しているだけですが、一応、Nileデータについても追加でやってみました。
なお、上記のリンク先は、「Rによるデータサイエンス」で書籍化もされています。一般的な統計手法が網羅されており、辞書的にも使えるのでとても良い書籍と思います。

ARモデル

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

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で指定できます。指定可能なアルゴリズムは４つ存在し、ユール-ウォーカー法、最小２乗法、最尤法、バーグ法があります。デフォルトはユール-ウォーカー法が選択されます。ユール-ウォーカー法の数理的な話はこちらが参考になります。

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

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

> acf(lh)

> pacf(lh)

> acf(Nile)

> pacf(Nile)

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

> 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

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

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

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

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

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

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

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

ARモデルでの予測

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

> 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)))

> 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)))

最後に

先に書いた通り、Rと時系列(2)に習って試し、追加でNileデータにもARモデルを適用してみました。モデルの診断など足りない箇所については後日、理解の範囲で記載できればと思います。

参考文献

Rと時系列(2)
ARモデルとユール–ウォーカー方程式