Python

HOME
Python
【時系列】ARモデルをわかりやすく解説|Yule-Walker法や最尤法も

05/03/2022 / 最終更新日時 : 01/16/2024 bluest Python

【時系列】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')で指定可能