【時系列】ARモデルをわかりやすく解説|Yule-Walker法や最尤法も
こんにちは、青の統計学です。
今回解説するのは、時系列モデルの基礎であるARモデルです。
まずは数式を見てみましょう。
ARモデル(autoregression model)
$$y_{n} = \sum_{j=1}^{m}a_jy_{n-j}+v_n$$
\(y_n\)：定常時系列
\(v_n\)：正規白色雑音(ホワイトノイズ)=系列相関がない
\(m\)：自己回帰の次数
\(a_j\)：自己回帰の係数
時系列モデルは、回帰係数以外にも次数も推定せねばならず、次数によって結果が大きく変わってしまいます。
前期の値が翌期の値に影響を与えているのがわかりますね。
株価指数や失業率などの経済指標によく使われ、社会科学の分野では頻出します。
ARモデルによる予測
さて、仮に1期先の予測をしてみます。
予測誤差とその性質を確認しましょう。
まずはARモデルです
$$y_{n+1} = a_1y_{n}+…a_my_{n-m+1}+v_{n+1}$$
最終項がノイズで、それ以外は現在までのデータで決まる部分ですね。
次に1期先の予測です。
ノイズを除いた形ですね。
$$y_{n+1|n} = a_1y_{n}+…a_my_{n-m+1}$$
残差と考え方は同じで、一期先の予測誤差は\(ε_{n+1}=y_{n+1}-y_{n+1|n}=v_{n+1}\)となります。
正規白色雑音を想定しているので、予測誤差の性質は以下のようになります。
$$E[ε_{n+1}]=E[v_{n+1}]=0$$
$$E[ε_{n+1}^2]=E[v_{n+1}^2]=σ^2$$
ちなみに移動平均項を加えたARMAモデルや季節性を考慮したSARIMAモデルの方が一般的なのに、なぜARモデルを学ぶ必要があるのか疑問に感じた方も多いかもしれません。
というのもARMAモデルなどは自由度が高すぎて不安定、つまりどのようなデータにもfitしてしまうので、予測タスクが困難になったりする場合があります。
一方、ARモデルの推定は比較的容易(OLSなど)で、多変量モデルの場合など実用性が高いのが特徴です。
どの状況でも一方が優れているということはないので、適切に使い分けたいですね。
VARモデル(ベクトル自己回帰モデル)
ちなみにARモデルの多変数版としてVARモデルがあります。
$$X_t=\sum_{i=1}^{p}A_iX_{t-i}$$
VARモデルは複数の時系列変数がお互いにどのように影響し合うかをモデル化します。
各変数は過去の自分自身の値と他の変数の値の関数として表されます。
よく状態空間モデルとどう違うのか、と言われますが全然違います。
VARは過去の値を用いて将来の値を予測しますが、状態空間モデルほどの動的な更新は行いません。
主に定常過程に適しており、ラグを用いて予測を行います。
状態空間モデルは、状態方程式を通じて時間と共に状態を更新するため、動的なシステムや非定常過程に適しています。
将来予測、つまり外挿において適しているのは状態空間モデルです。
【時系列】状態空間モデルをわかりやすく解説|カルマンフィルタの仕組み
Yule-walker法
さて、脱線してしまいました。
ARモデルの推定の方法としていくつかご紹介します
Yule-Waler法は、ARモデルのパラメータを推定する一般的な方法の一つで、以下のような自己相関係数\(C_k\)を作ります。
\(C_k=\sum_{j=1}^ma_jC_{j-k},(k=1,2,…\)
\(C_0=\sum_{j=1}^ma_jC_{j}+σ^2\)
k=0とそれ以外で分けています。
$$\begin{equation}\begin{bmatrix}C_0 & C_{-1} & C_{-2} & \cdots & C_{m-1} \\ C_1 & C_0 & C_{-1} & \cdots & C_{m-2} \\ C_2 & C_1 & C_0 & \cdots & C_{m-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ C_{m-1} & C_{m-2} & C_{m-3} & \cdots & C_0\end{bmatrix}\begin{bmatrix}a_1 \\ a_{2} \\ \vdots \\ a_{m} \end{bmatrix}=\begin{bmatrix}C_1 \\ C_2 \\ \vdots \\ C_m\end{bmatrix}\end{equation}$$
そして上の方程式を解いて\(a_j\)を求めます。
一番右はtoepritz行列と呼ばれ、一次方程式が簡単に解ける行列です。
ただ、これは何をやっているのでしょうか？
自己相関係数と予測誤差の分散の関係をみてみればわかります。
$$σ^2=E[v_n^2]=E[(y_n-\sum_{j=1}^ma_jy_{n-j})^2]$$
$$=C_0-2\sum_{j=1}^ma_jC_j+\sum_{j=1}^m \sum_{i=1}^ma_ja_iC_{i-j}$$
このように予測誤差の分散の期待値は自己相関係数\(C\)を使って、上のように表せます。
パラメータ\(a_j\)で微分します。
$\frac{∂σ^2}{∂a_{j}}=-2C_j+2\sum{i=1}^ma_iC_{i-j}=0$
予測誤差の分散の期待値を最小にするように係数\(a_j\)を決めていると言えます。
この方法は計算が比較的簡単で、効率的ですが、モデルの次数が高くなると精度が低下することがあります。
普通の時系列タスクで次数\(m\)は未知です。
導出は省きますが、\(\frac{M^4}{12}\)となります。
その代わりに、Levinson’s Algorithm を通すと、Yule-walker方程式を効率よく解け、計算量は\(2M^2\)程度になります。
最尤法
自己回帰モデルでは、各時点での値がその前のいくつかの時点の値に依存しています。
つまり、各時点\(t\)での値\(y_t\)は過去の値\(y_{t-1},…y_{t-m}\)に条件付けられています。
$$L(a)=p(y_1,…,y_N|θ)$$
$$L(a_1,…,a_j,σ^2)=p(y_1,…,y_N|θ)$$
$$L(a_1,…,a_N,σ^2)=\prod_{n=1}^{N}p(y_n|y_1,…,y_{N-1},a_1,…,a_N)$$
同時分布を各時点での条件付き分布に分解していきます。
この尤度関数を最大化することで、モデルのパラメータを推定します。
尤度関数については以下のコンテンツで詳しく解説しています。
【AICで使う】KL divergence(カルバック-ライブラー情報量)をわかりやすく解説|python
この条件付き分布の誤差項\(v_n\)が正規分布に従うと仮定することで、次のように導出されます。
$$y_{n}|y_{n-1},…,y_{n-m}〜N(a_1y_{n-1}+…+a_my_{n-m},σ^2)$$
$$L(a_1,…,a_N)=p(y_1,…,y_n|a_1,…,a_N)=(2π)^{-\frac{n}{2}}|C|^{-\frac{1}{2}}exp(-\frac{1}{2}(x-μ)^TC^{-1}(x-μ))$$
\(l(a_1,…,a_N)=log L(a_1,…,a_N)\)の数値的最適化により、最尤推定値を求めます。
これはNが大きくなればなるほど、計算量がかなり増えてしまいます。
ただ、線形の数値最適を行うので精度は良いですね。
これまで紹介した方法は、pythonのstatsmodelで引数として指定できます。
メリットデメリットがあるので、状況に応じて使い分けてください。
import statsmodels.api as sm
model = sm.tsa.AR(data)
result = model.fit(method='yw') # Yule-Walker法を使用
#method='yw'ならYule-Walker法
#method='mle'なら最尤法
#最小二乗法はARMAクラスで使用され、ARMA.fit(method='ls')で指定可能事前知識|定常性とは
ARモデルを作るにあたり、定常性が前提となります。
定常性とは、時系列データの統計的な性質が時間依存しないことです。
定常は弱定常と強定常に分類されます。
弱定常
・平均と分散(2次のモーメントまで)が変化しない
・自己共分散がラグ\(h\)にのみ依存する。
→2時点間の自己強分散関数はそれらの時間差ラグ\(h\)にのみ決定し、\(t-1\)期と\(t+1\)期間の自己共分散は、\(t+1\)期と\(t+3\)期間の自己共分散と変わらないということです。
どちらも\(h=2\)で一緒だからですね。
強定常
・時間によって確率分布が変化しない
補足：グレンジャー因果(Granger causality)
今回解説したARモデル(自己回帰モデル)に関して、グレンジャー因果というものをご紹介します。
$$y_{t} = a_1y_{t-1}+a_2y_{t-2}+..+a_py_{t-p}$$
この自己回帰モデルを基準モデルとしたときに、外生変数x_t-τを加えた以下のモデルを考えます。
$$y_{t} = a_1y_{t-1}+a_2y_{t-2}+..+a_py_{t-p}+b_1y_{t-1}+b_2y_{t-2}+..+a_py_{t-p}$$
外生変数を加えたモデルの方が、元のモデルよりも\(y_t\)の予測誤差が減少する時に、「xはyにグレンジャー因果」を持つと言います。
このグレンジャー因果は、因果関係の中でもかなり弱い因果であり、交絡因子の問題を解決できないなどのデメリットがあります。
必ずしも正確な因果ではないが、予測には有用な関係があるということを示しています。
補足：ARモデルの最適な次数pの決定にAICを使う。
時系列モデルの次数の決定方法として、AICというものがあります。
AIC(赤池情報量基準)を詳しく知りたい方は、以下をご覧ください。
自己回帰モデルは、以下のようにAR(p)とします。
$$AR(p)：Y_t = β_0+β_1Y_{t-1}+β_2X_{t-2}+…+β_pX_{t-p}+ε_t$$
AR(p)の次数pは、ラグ変数と呼ばれ、ARモデルにおける説明変数の個数です。
例えば、AR(1)なら一期前のデータまでモデルに入っているモデルであり、その係数\(φ=1\)であれば単位根とかランダムウォークと呼ばれます。
さて、このモデルを推定しようと考える時、pを決定するためにAICを使います。
$$AIC(p)=log(\frac{RSS(p)}{T})+(p+1){\frac{2}{T}}$$
ここでRSS(p)とは、AR(p)モデルにおける残差平方和です。Tはサンプルサイズですね。
残差平方和は、以下のコンテンツで復習できます。
そして、AICを最小にするpが最適なpになります。
最大のラグ次数をp<5などと設定しておき、一つずつAICをチェックするということです。
以下のように、第1項は残差平方和の単調変換なので、ラグ次数pが大きくなるほど単調に減少します。
→データの説明力は上がる
一方、第2項は説明変数の個数へのペナルティを表すので、ラグ次数pが大きくなるほど単調に増加します。
→データの簡潔性は失われる
| ラグ次数p | 0 | 1 | 2 | 3 |
| 第1項 | 403.7 | 118.6 | 118.4 | 117.6 |
| 第2項 | 1.031 | -0.194 | -0.195 | -0.202 |
| AIC | 0.014 | 0.028 | 0.042 | 0.056 |
AICは時系列モデル以外にも、多くのモデルの予測精度の指標の一つとして使われます。
時系列モデル以外では、以下のように最大対数尤度と-2の積を第1項にし、パラメーター数と2の積を第2項にします。
$$AIC=-2logL + 2k$$
基本的に意味合いは同じで、データの説明力が上がると第1項は小さくなりますが、同時にパラメーターも増えるので第2項も増えてしまい、そこをうまくバランスさせているということです。
