はじめに
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モデルを適用してみました。モデルの診断など足りない箇所については後日、理解の範囲で記載できればと思います。