はじめに
この記事では、自己回帰モデル（ARモデル）や、ベクトル自己回帰モデル（VARモデル）のパラメータをレビンソン・ダービン（Levinson-Durbin）法という方法で、推定する方法を説明します。
レビンソン・ダービン法は、次の資料が参考になります。
本記事では、上記の参考資料には含まれていない次の内容を含めて、かつ、本記事だけで完結する内容にしました。
- 参考資料では、「$i$次のパラメータを使った$i+1$次のパラメータの導出」の方法が書かれていましたが、本記事では、次のように具体的な例で示しました。
- 0次のパラメータを導出
- 0次のパラメータを使って1次のパラメータを導出
- 1次のパラメータを使って2次のパラメータを導出
- 2次のパラメータを使って3次のパラメータを導出
- 参考資料では、VARモデルのパラメータ推定方法（すなわち、レビンソン・ダービン法の多次元化）の導出について触れられていなかったので、そこを詳しく記述しました。
- カラーの図を使ってアルゴリズムの流れがわかるようにしました。VARモデルのパラメータを推定する流れを図で表現したものがこちらです。美しい構造をしています。
ARモデルの推定
自己回帰モデル（ARモデル）
x_t = a_0 + a_1 x_{t-1} + a_2 x_{t-2} + \cdots + a_k x_{t-k} + \epsilon_t, \ \ \ \epsilon_t \sim N(0, \sigma)\ \ \ \ \ [A]
の係数$a$と分散$\sigma$を推定します。
一般的にはARモデルの誤差項$\epsilon$の分散（1次元正規分布の分散）は$\sigma^2$と記述されますが、ここでは、表記を単純化するためと、VARモデルと記載を統一するために、$\sigma$とします。
自己共分散から連立方程式を作る
まず、表記を簡単にするために平均が0になるように調整します。$x_t$の平均$\mu$は、観測データから推定値を求めることができるので、これを既知とします。新しい時系列データ$y_t$を次の式で定義します。
y_t = x_t - \mu
こうすると、$E(y_t)=0$となるので、式[A]は次のように変形することができます。
y_t = a_1 y_{t-1} + a_2 y_{t-2} + \cdots + a_k y_{t-k} + \epsilon_t, \ \ \ \epsilon_t \sim N(0, \sigma)
となります。$a_0$を消去して問題を単純化できました。 ここから$a_1, a_2, \cdots, a_k$と$\sigma$を求めていきます。
$a_0$は、$a_1, a_2, \cdots, a_k$を求めた後に次の方法で求めることができます。
a_0 = (1-a_1-a_2-\cdots -a_k)\mu
なぜならば、
\begin{align}
&y_t = a_1 y_{t-1} + a_2 y_{t-2} + \cdots + a_k y_{t-k} + \epsilon_t\\
\Leftrightarrow\ &x_t-\mu = a_1 (x_{t-1} - \mu) + a_2 (x_{t-2}-\mu) + \cdots + (a_k x_{t-k}-\mu) + \epsilon_t\\
\Leftrightarrow\ &
x_t = (1-a_1-a_2-\cdots- a_k)\mu + a_1 x_{t-1} + a_2 x_{t-2} + \cdots + a_k x_{t-k} + \epsilon_t
\end{align}
だからです。
では、$a_1, a_2, \cdots, a_k$と$\sigma$を求めていきます。
両辺に$y_{t-i}(i=0,1,\cdots,k)$をかけて、期待値をとります。
\begin{align}
E(y_t y_t)\ \ \ \ &= a_1 E(y_{t-1} y_t) \ \ \ + a_2 E(y_{t-2} y_t) \ \ \ \ + \cdots + a_k E(y_{t-k} y_t)\ \ \ \ + E(\epsilon_t y_t)\\
E(y_t y_{t-1}) \ &= a_1 E(y_{t-1} y_{t-1}) + a_2 E(y_{t-2} y_{t-1}) + \cdots + a_k E(y_{t-k} y_{t-1}) + E(\epsilon_t y_{t-1})\\
E(y_t y_{t-2}) \ &= a_1 E(y_{t-1} y_{t-2}) + a_2 E(y_{t-2} y_{t-2}) + \cdots + a_k E(y_{t-k} y_{t-2}) + E(\epsilon_t y_{t-2})\\
\vdots\\
E(y_t y_{t-k}) &= a_1 E(y_{t-1} y_{t-k}) + a_2 E(y_{t-2} y_{t-k}) + \cdots + a_k E(y_{t-k} y_{t-k})\ + E(\epsilon_t y_{t-k})
\end{align}
$E(y_t) = 0$なので、$E(y_t y_{t-i})$ は $y_t$ と $y_{t－i}$ の共分散、すなわち、$i$次の自己共分散を示します。$C_i = E(y_t y_{t-i})$とすると、
\begin{align}
C_0\ &= a_1 C_{-1} \ \ + a_2 C_{-2} + \cdots + a_k C_{-k}\ \ \ \ + \sigma \\
C_1\ &= a_1 C_0 \ \ \ \ + a_2 C_{-1} + \cdots + a_k C_{-(k-1)}\\
C_2\ &= a_1 C_1 \ \ \ \ + a_2 C_0 \ \ + \cdots + a_k C_{-(k-2)}\\
\vdots\\
C_k\ &= a_1 C_{k-1} + a_2 C_{k-2} + \cdots + a_k C_0
\end{align}
となります。これをベクトル・行列化すると
\begin{align}
&
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} & & C_{-k} \\
C_1 & C_0 & C_{-1} & \cdots & C_{-(k-1)} \\
C_2 & C_1 & C_0 & & C_{-(k-2)} \\
& \vdots & &\ddots& \\
C_k & C_{k-1} & C_{k-2} & & C_0
\end{bmatrix}
\begin{bmatrix}
1 \\ - a_1\\ \vdots\\ - a_k
\end{bmatrix}
=
\begin{bmatrix}
\sigma \\ 0\\ \vdots\\ 0
\end{bmatrix}\\
\Leftrightarrow&
\begin{bmatrix}
C_0 & C_1 & C_2 & & C_k \\
C_1 & C_0 & C_1 & \cdots & C_{k-1} \\
C_2 & C_1 & C_0 & & C_{k-2} \\
& \vdots & &\ddots& \\
C_k & C_{k-1} & C_{k-2} & & C_0
\end{bmatrix}
\begin{bmatrix}
1 \\ - a_1\\ \vdots\\ - a_k
\end{bmatrix}
=
\begin{bmatrix}
\sigma \\
0\\ \vdots\\ 0
\end{bmatrix}
\end{align}
となります。変形には、$C_i = E(y_t y_{t-i}) = E(y_{t-i} y_t ) = C_{-i}$ という性質を使いました。
自己共分散$C_0, C_1, \cdots, C_k$は、データから推定値を求めることができるので、これらを既知とすると、上式は未知数が$k+1$個（$\sigma, a_1, \cdots, a_k$）で式が$k+1$個の連立方程式なので、解くことができます。
レビンソン・ダービン法で連立方程式を解く
レビンソン・ダービン法では、$i$次のARモデルのパラメータを求めるのに、$i-1$次のARモデルのパラメータを使用します。これにより、高速に解くことができます。
ここでは、例として3次のARモデルのパラメータを求めてみます。レビンソン・ダービン法では、いきなり3次のARモデルを求めようとしません。次の手順で行います。
- [step0] 0次のARモデルのパラメータを求める
- [step1] 0次の結果を利用して、1次のARモデルのパラメータを求める
- [step2] 1次の結果を利用して、2次のARモデルのパラメータを求める
- [step3] 2次の結果を利用して、3次のARモデルのパラメータを求める
それぞれのstepで解くべき式は次の通りです。
-
[step0] 0次ARモデルで解くべき式
C_0 = \sigma_{(0)} -
[step1] 1次ARモデルで解くべき式
\begin{bmatrix} C_0 & C_1 \\ C_1 & C_0 \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(1)} \end{bmatrix} = \begin{bmatrix} \sigma_{(1)} \\ 0 \end{bmatrix} -
[step2] 2次ARモデルで解くべき式
\begin{bmatrix} C_0 & C_1 & C_2 \\ C_1 & C_0 & C_1 \\ C_2 & C_1 & C_0 \\ \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(2)}\\ - a_2^{(2)} \end{bmatrix} = \begin{bmatrix} \sigma_{(2)} \\ 0\\ 0 \end{bmatrix} -
[step3] 3次ARモデルで解くべき式
\begin{bmatrix} C_0 & C_1 & C_2 & C_3 \\ C_1 & C_0 & C_1 & C_2 \\ C_2 & C_1 & C_0 & C_1 \\ C_3 & C_2 & C_1 & C_0 \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(3)}\\ - a_2^{(3)}\\ - a_3^{(3)} \end{bmatrix} = \begin{bmatrix} \sigma_{(3)} \\ 0\\ 0\\ 0 \end{bmatrix}
$a$や$\sigma$の解は、ARモデルの次数によって異なるので、$^{()}$を付けて区別しています。
[step0] 0次のARモデルのパラメータを求める
解くべき式は
C_0 = \sigma_{(0)}
でした。
未知数$\sigma_{(0)}$は次のようにして求めることができます（左辺と右辺を入れ替えるだけです）。
\sigma_{(0)} = C_0 \ \ \ \ [0.1]
[step1] 0次の結果を利用して、1次のARモデルのパラメータを求める
解くべき式は、
\begin{bmatrix} C_0 & C_1 \\ C_1 & C_0 \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(1)} \end{bmatrix} = \begin{bmatrix} \sigma_{(1)} \\ 0 \end{bmatrix}
でした。
1行目に着目すると
C_0 - a_1^{(1)} C_1 = \sigma_{(1)} \ \ \ \ \ [1.1]
2行目に着目すると
a_1^{(1)} C_0 = C_1 \ \ \ \ [1.2]
式[0.1]の結果を使うと、次のようにまとめることができます。
\begin{align}
a_1^{(1)} &= C_1 \sigma_{(0)}^{-1} \\
\sigma_{(1)} &= \sigma_{(0)} - a_1^{(1)}a_1^{(1)} \sigma_{(0)}
\end{align}
[step2] 1次の結果を利用して、2次のARモデルのパラメータを求める
解くべき式は、
\begin{bmatrix} C_0 & C_1 & C_2 \\ C_1 & C_0 & C_1 \\ C_2 & C_1 & C_0 \\ \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(2)}\\ - a_2^{(2)} \end{bmatrix} = \begin{bmatrix} \sigma_{(2)} \\ 0\\ 0 \end{bmatrix}
でした。
解くべき式の1行目に着目すると、
C_0 - a_1^{(2)}C_1 - a_2^{(2)} C_2 = \sigma_{(2)} \ \ \ \ \ [2.1]
解くべき式の2,3行目に着目すると、
\begin{align}
&\begin{bmatrix}
C_0 & C_1 \\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_1^{(2)}\\a_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2
\end{bmatrix}\ \ \ \ \ [2.2]\\
\Leftrightarrow&
\begin{bmatrix}
C_0 & C_1 \\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_2^{(2)}\\a_1^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
C_2\\C_1
\end{bmatrix}\ \ \ \ \ [2.3]
\end{align}
式[1.2]の結果を用いて、解くべき式の2行目を変形します。
\begin{align}
& a_1^{(2)} C_0 = C_1 - a_2^{(2)} C_1\\
\Leftrightarrow\ & a_1^{(2)} C_0 = a_1^{(1)} C_0 - a_2^{(2)} a_1^{(1)} C_0 \\
\Leftrightarrow\ & a_1^{(2)} = a_1^{(1)} - a_2^{(2)} a_1^{(1)}
\end{align}
つまり、$a_1^{(2)}$ は $a_2^{(2)}$ から求めることができます。
この結果と式[1.1][1.2]を使って、解くべき式を変形していくことで、$a_2^{(2)}$を求めます。ついでに$\sigma_{(2)}$も求めます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_1 & C_0 & C_1 \\
C_2 & C_1 & C_0 \\
\end{bmatrix}
\begin{bmatrix}
1 \\ - a_1^{(2)}\\- a_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(2)} \\0\\0
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_1 & C_0 & C_1 \\
C_2 & C_1 & C_0 \\
\end{bmatrix}
\begin{bmatrix}
1 \\ - a_1^{(1)} + a_2^{(2)} a_1^{(1)}\\ -a_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(2)} \\ 0\\ 0
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_1 & C_0 & C_1 \\
C_2 & C_1 & C_0 \\
\end{bmatrix}
\begin{bmatrix}
\begin{bmatrix}
1 \\ - a_1^{(1)}\\ 0
\end{bmatrix}
- a_2^{(2)}
\begin{bmatrix}
0 \\ - a_1^{(1)}\\ 1
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(2)} \\ 0\\ 0
\end{bmatrix}\\
\Leftrightarrow\ & \begin{bmatrix}
\sigma_{(1)} \\ 0\\ C_2 - a_1^{(1)} C_1
\end{bmatrix}
- a_2^{(2)}
\begin{bmatrix}
C_2 - a_1^{(1)} C_1 \\ 0\\ \sigma_{(1)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(2)} \\
0\\
0
\end{bmatrix}
\end{align}
1番下の行に着目すると
a_2^{(2)} =
(C_2 - a_1^{(1)} C_1) \sigma_{(1)}^{-1}\\
1番上の行に着目すると
\begin{align}
\sigma_{(2)} &= \sigma_{(1)} - a_2^{(2)} (C_2 - a_1^{(1)} C_1 )\\
&= \sigma_{(1)} - {a_2^{(2)}} a_2^{(2)} \sigma_{(1)}
\end{align}
となります。
[step3] 2次の結果を利用して、3次のARモデルのパラメータを求める
解くべき式は
\begin{bmatrix} C_0 & C_1 & C_2 & C_3 \\ C_1 & C_0 & C_1 & C_2 \\ C_2 & C_1 & C_0 & C_1 \\ C_3 & C_2 & C_1 & C_0 \end{bmatrix} \begin{bmatrix} 1 \\ - a_1^{(3)}\\ - a_2^{(3)}\\ - a_3^{(3)} \end{bmatrix} = \begin{bmatrix} \sigma_{(3)} \\ 0\\ 0\\ 0 \end{bmatrix}
でした。
解くべき式の1行目に着目すると、
C_0 - a_1^{(3)}C_1 - a_2^{(3)} C_2 - a_3^{(3)} C_3 = \sigma_{(3)} \ \ \ \ \ [3.1]
解くべき式の2,3,4行目に着目すると、
\begin{align}
&
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_1 & C_0 & C_1 \\
C_2 & C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_1^{(3)}\\a_2^{(3)}\\a_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2\\C_3
\end{bmatrix}\ \ \ \ \ [3.2]\\
\Leftrightarrow&\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_1 & C_0 & C_1 \\
C_2 & C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_3^{(3)}\\a_2^{(3)}\\a_1^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_3\\C_2\\C_1
\end{bmatrix}\ \ \ \ \ [3.3]
\end{align}
式[2.2][2.3]の結果を用いて、解くべき式の2,3行目を変形します。
\begin{align}
&\begin{bmatrix}
C_0 & C_1\\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_1^{(3)}\\a_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2\\
\end{bmatrix}
-a_3^{(3)}
\begin{bmatrix}
C_2\\C_1\\
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1\\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_1^{(3)}\\a_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_0 & C_1\\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_1^{(2)}\\a_2^{(2)}
\end{bmatrix}
-a_3^{(3)}
\begin{bmatrix}
C_0 & C_1\\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
a_2^{(2)}\\a_1^{(2)}
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
a_1^{(3)}\\a_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
a_1^{(2)}\\a_2^{(2)}
\end{bmatrix}
-a_3^{(3)}
\begin{bmatrix}
a_2^{(2)}\\a_1^{(2)}
\end{bmatrix}
\end{align}
つまり、$a_1^{(3)}, a_2^{(3)}$ は $a_3^{(3)}$ から求めることができます。
この結果と式[2.1][2.2]を使って、解くべき式を変形していくことで、$a_3^{(3)}$を求めます。ついでに$\sigma_{(3)}$も求めます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_1 & C_0 & C_1 & C_2 \\
C_2 & C_1 & C_0 & C_1 \\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
1 \\
- a_1^{(3)}\\
- a_2^{(3)}\\
- a_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(3)} \\
0\\
0\\
0
\end{bmatrix}\\
\Leftrightarrow
&
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_1 & C_0 & C_1 & C_2 \\
C_2 & C_1 & C_0 & C_1 \\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
1 \\
-a_1^{(2)}+a_3^{(3)}a_2^{(2)}\\
-a_2^{(2)}+a_3^{(3)}a_1^{(2)}\\
- a_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(3)} \\
0\\
0\\
0
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_1 & C_0 & C_1 & C_2 \\
C_2 & C_1 & C_0 & C_1 \\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
\begin{bmatrix}
1 \\ - a_1^{(2)}\\ - a_2^{(2)}\\ 0
\end{bmatrix}
- a_3^{(3)}
\begin{bmatrix}
0 \\ - a_2^{(2)}\\ - a_1^{(2)}\\ 1
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(3)} \\ 0\\ 0\\ 0
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
\sigma_{(2)} \\
0\\ 0\\ C_3 - a_1^{(2)} C_2 - a_2^{(2)} C_1
\end{bmatrix}
- a_3^{(3)}
\begin{bmatrix}
C_3 - a_1^{(2)} C_2 - a_2^{(2)} C_1 \\ 0\\ 0\\ \sigma_{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\sigma_{(3)} \\ 0\\ 0\\ 0
\end{bmatrix}
\end{align}
1番下の行に着目すると
a_3^{(3)} = (C_3 - a_1^{(2)} C_2 - a_2^{(2)} C_1) \sigma_{(2)}^{-1}
1番上の行に着目すると
\sigma_{(3)} = \sigma_{(2)} - a_3^{(3)} (C_3 - a_1^{(2)} C_2 - a_2^{(2)} C_1) = \sigma_{(2)} - {a_3^{(3)}} a_3^{(3)} \sigma_{(2)}
まとめ）K次のARモデルのパラメータを求めるには
4次以上も3次と同じ方法で求めることができます。$K$次のARモデルのパラメータを求める方法は次の通りです。
- 0,1次を求める
\begin{align}
\sigma_{(0)} &= C_0\\
a_1^{(1)} &= C_1 \sigma_{(0)}^{-1} \\
\sigma_{(1)} &= \sigma_{(0)} - a_1^{(1)}a_1^{(1)} \sigma_{(0)}
\end{align}
- $k=2,\cdots,K$を次の式に順に代入し、逐次的に$K$次まで求める
\begin{align}
{a_k^{(k)}} &= \Bigl(C_k - \sum_{i=1}^{k-1} a_i^{(k-1)}C_{k-i}\Bigr)\sigma_{(k-1)}^{-1} \\
{a_i^{(k)}} &= {a_i^{(k-1)}} - {a_k^{(k)}}{a_{k-i}^{(k-1)}} \ \ \ (i=1,\cdots,k-1)\\ \\
\sigma_{(k)} &= \sigma_{(k-1)} - {a_k^{(k)}}{a_k^{(k)}}\sigma_{(k-1)}
\end{align}
図にすると、次のようなイメージです。
Pythonのコードの例は次の通りです。うまくパラメータを推定できていることが分かります。
（ただし、わかりやすさ重視のコードのため、遅いと思います。for文をベクトル演算に置き換えると高速になります。）
import numpy as np
import scipy
np.random.seed(42)
#######################
# データ作成
#######################
n = 100000 # データ数
a = [1/2, 1/3, - 1/4, 1/3] # 係数
sigma = 1 # 誤差の標準偏差
k = 3 # 次数
x = np.zeros([n])
epsilons = scipy.stats.norm.rvs(loc=0, scale=sigma, size=n)
for i in range(k):
x[i] = epsilons[i] # x[0] = 乱数, x[1] = 乱数, x[2] = 乱数
for i in range(k, n):
# x[i] = a[0] + a[1] * x[i - 1] + a[2] * x[i - 2] + + a[3] * x[i - 3] + 乱数
x[i] = a[0]
for j in range(1, k + 1):
x[i] += a[j] * x[i - j]
x[i] += epsilons[i]
##############################################
# レビンソン・ダービン法でa, sigmaを推定
##############################################
# 平均を0に
mu = np.mean(x)
y = x - mu
# C[0], C[1], C[2], C[3]を導出
C = []
for i in range(k + 1):
s = 0
for j in range(n - i):
s += y[j + i] * y[j]
s /= (n - i)
C.append(s)
# 3次のARモデルの係数をレビンソン・ダービン法で導出
for _k in range(k + 1):
if _k == 0: # 0次
sigma = C[0]
a = [mu]
else: # 1次、2次、3次
a_old = a.copy()
a = [None] * (_k + 1)
s = 0
for i in range(1, _k):
s += a_old[i] * C[_k - i]
a[_k] = (C[_k] - s) / sigma
for i in range(1, _k):
a[i] = a_old[i] - a[_k] * a_old[_k - i]
sigma = sigma - a[_k] ** 2 * sigma
a[0] = (1 - np.sum(a[1:])) * mu
print(a)
# [0.4970878153369187, 0.3339784918894054, -0.2492625989152757, 0.3364405532924277] (真の値 [0.5 , 0.33, -0.25, 0.33])
print(sigma) # 1.0017802899102914 (真の値 1.0)
VARモデルの推定
ベクトル自己回帰モデル（VARモデル）
\boldsymbol{x}_t = \boldsymbol{a}_0 + A_1 \boldsymbol{x}_{t-1} + A_2 \boldsymbol{x}_{t-2} + \cdots + A_k \boldsymbol{x}_{t-k} + \boldsymbol{\epsilon}_t, \ \ \ \boldsymbol{\epsilon}_t \sim N(\boldsymbol{0}, \Sigma)\ \ \ \ \ [B]
の係数$\boldsymbol{a},A$と分散共分散行列$\Sigma$を推定する方法を説明します。
式は基本的に、行列・ベクトルで表現しています。$\boldsymbol{x}$の次元数（予測対象の特徴量数）を$m$とすると、次元数は下のようになります。
- $\boldsymbol{a}$: $m$次元ベクトル
- $A$: $m \times m$次元行列
- $\boldsymbol{\epsilon}$: $m$次元ベクトル
- $\boldsymbol{0}$: $m$次元ベクトル
- $\Sigma$: $m \times m$次元行列
相互共分散行列から連立方程式を作る
まず、表記を簡単にするために平均が$\boldsymbol{0}$になるように調整します。$\boldsymbol{x}_t$の平均$\boldsymbol{\mu}$は、観測データから推定値を求めることができるので、これを既知とします。新しい時系列データ$\boldsymbol{y}_t$を次の式で定義します。
\boldsymbol{y_t} = \boldsymbol{x}_t - \boldsymbol{\mu}
こうすると、$E(\boldsymbol{y}_t)=\boldsymbol{0}$となるので、式[B]は次のように変形することができます。
\boldsymbol{y}_t = A_1 \boldsymbol{y}_{t-1} + A_2 \boldsymbol{y}_{t-2} + \cdots + A_k \boldsymbol{y}_{t-k} + \boldsymbol{\epsilon}_t, \ \ \ \boldsymbol{\epsilon}_t \sim N(\boldsymbol{0}, \Sigma)
となります。$\boldsymbol{a}_0$を消去して問題を単純化できました。 ここから$A_1, A_2, \cdots, A_k$と$\Sigma$を求めていきます。
$\boldsymbol{a}_0$は、$A_1, A_2, \cdots, A_k$を求めた後に次の方法で求めることができます。
\boldsymbol{a}_0 = (E-A_1-A_2-\cdots -A_k)\boldsymbol{\mu}
なぜならば、
\begin{align}
&\boldsymbol{y}_t = A_1 \boldsymbol{y}_{t-1} + A_2 \boldsymbol{y}_{t-2} + \cdots + A_k \boldsymbol{y}_{t-k} + \boldsymbol{\epsilon}_t\\
\Leftrightarrow\ &\boldsymbol{x}_t-\boldsymbol{\mu} = A_1 (\boldsymbol{x}_{t-1} - \boldsymbol{\mu}) + A_2 (\boldsymbol{x}_{t-2}-\boldsymbol{\mu}) + \cdots + (A_k \boldsymbol{x}_{t-k}-\boldsymbol{\mu}) + \boldsymbol{\epsilon}_t\\
\Leftrightarrow\ &
\boldsymbol{x}_t = (E-A_1-A_2-\cdots- A_k)\boldsymbol{\mu} + A_1 \boldsymbol{x}_{t-1} + A_2 \boldsymbol{x}_{t-2} + \cdots + A_k \boldsymbol{x}_{t-k} + \boldsymbol{\epsilon}_t
\end{align}
だからです。
ARモデルの時と同様に、平均を$\boldsymbol{0}$に調整済みの次の時系列データで考えることにします。平均を$\boldsymbol{0}$に調整済みのVARモデルは次の式で与えられます。
\boldsymbol{y}_t = A_1 \boldsymbol{y}_{t-1} + A_2 \boldsymbol{y}_{t-2} + \cdots + A_k \boldsymbol{y}_{t-k} + \boldsymbol{\epsilon}_t, \ \ \ \boldsymbol{\epsilon}_t \sim N(\boldsymbol{0}, \Sigma)
両辺に$\boldsymbol{y}_{t-i}^T(i=0,1,\cdots,k)$をかけて、期待値をとります。
\begin{align}
E(\boldsymbol{y}_t \boldsymbol{y}_t^T)\ \ \ \ &= A_1 E(\boldsymbol{y}_{t-1} \boldsymbol{y}_t^T) \ \ \ \ + A_2 E(\boldsymbol{y}_{t-2} \boldsymbol{y}_t^T) \ \ \ + \cdots + A_k E(\boldsymbol{y}_{t-k} \boldsymbol{y}_t^T)\ \ \ \ + E(\epsilon_t \boldsymbol{y}_t^T)\\
E(\boldsymbol{y}_t \boldsymbol{y}_{t-1}^T) \ &= A_1 E(\boldsymbol{y}_{t-1} \boldsymbol{y}_{t-1}^T)\ + A_2 E(\boldsymbol{y}_{t-2} \boldsymbol{y}_{t-1}^T)\ + \cdots + A_k E(\boldsymbol{y}_{t-k} \boldsymbol{y}_{t-1}^T) \ + E(\epsilon_t \boldsymbol{y}_{t-1}^T)\\
\vdots\\
E(\boldsymbol{y}_t \boldsymbol{y}_{t-k}^T) &= A_1 E(\boldsymbol{y}_{t-1} \boldsymbol{y}_{t-k}^T) + A_2 E(\boldsymbol{y}_{t-2} \boldsymbol{y}_{t-k}^T) + \cdots + A_k E(\boldsymbol{y}_{t-k} \boldsymbol{y}_{t-k}^T)\ \ \ + E(\epsilon_t \boldsymbol{y}_{t-k}^T)
\end{align}
$ E(\boldsymbol{y}_t)=\boldsymbol{0} $ なので、$E(\boldsymbol{y}_t \boldsymbol{y}_{t-i}^T)$ は $\boldsymbol{y}_t$ と $\boldsymbol{y}_{t－i}$ の共分散、すなわち、$i$次の相互共分散を示します。$C_i = E(\boldsymbol{y}_t \boldsymbol{y}_{t-i}^T)$とすると、
\begin{align}
C_0\ &= A_1 C_{-1} \ \ + A_2 C_{-2} \ \ + \cdots + A_k C_{-k}\ \ \ \ + \Sigma \\
C_1\ &= A_1 C_0 \ \ \ \ + A_2 C_{-1} \ \ + \cdots + A_k C_{-(k-1)}\\
C_2\ &= A_1 C_1 \ \ \ \ + A_2 C_{0} \ \ \ \ + \cdots + A_k C_{-(k-2)}\\
\vdots\\
C_k &= A_1 C_{k-1} + A_2 C_{k-2} + \cdots + A_{-k} C_0
\end{align}
これをベクトル・行列化すると
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} & & C_{-k} \\
C_1 & C_0 & C_{-1} & \cdots & C_{-(k-1)} \\
C_2 & C_1 & C_0 & & C_{-(k-2)} \\
& \vdots & &\ddots& \\
C_k & C_{k-1} & C_{k-2} & & C_0
\end{bmatrix}
*
\begin{bmatrix}
E \\- A_1\\- A_2\\\vdots\\- A_k
\end{bmatrix}
=
\begin{bmatrix}
\Sigma \\O\\\vdots\\O
\end{bmatrix}
となります。$C,A,\Sigma,E,O$は$m\times m$次元の行列であることに注意してください。
ここで、独自の記号である$*$について説明します。
通常の行列の積を使った場合には下のようになります。
\begin{align}
\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}
\begin{bmatrix}
E \\- A_1
\end{bmatrix}
&=
\begin{bmatrix}
C_0 E - C_{-1}A_1\\
C_1 E - C_0 A_1
\end{bmatrix}
\\
&=
\begin{bmatrix}
C_0 - C_{-1}A_1\\
C_1 - C_0 A_1
\end{bmatrix}
\end{align}
$*$を使った場合は下のようになるとします。
\begin{align}
\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}
*
\begin{bmatrix}
E \\- A_1
\end{bmatrix}
&=
\begin{bmatrix}
E C_0 - A_1 C_{-1}\\
E C_1 -A_1 C_0
\end{bmatrix}
\\
&=
\begin{bmatrix}
C_0 - A_1 C_{-1}\\
C_1 -A_1 C_0
\end{bmatrix}
\end{align}
1次元の時とは違い$C_i=C_{-i}$は成立しないので、この形にすることはできません。
\begin{bmatrix}
C_0 & C_1 & C_2 & & C_k \\
C_1 & C_0 & C_1 & \cdots & C_{k-1} \\
C_2 & C_1 & C_0 & & C_{k-2} \\
& \vdots & &\ddots& \\
C_k & C_{k-1} & C_{k-2} & & C_0
\end{bmatrix}
*
\begin{bmatrix}
E \\- A_1\\- A_2\\\vdots\\- A_k
\end{bmatrix}
=
\begin{bmatrix}
\Sigma \\O\\\vdots\\O
\end{bmatrix}
これが原因で、1次元の時と同様の方法で解くことができず、「逆向きのVARモデル」というものを導入する必要があります。
逆向きのVARモデルは次の式で与えられます。
\boldsymbol{y}_t = B_1 \boldsymbol{y}_{t+1} + B_2 \boldsymbol{y}_{t+2} + \cdots + B_k \boldsymbol{y}_{t+k} + \boldsymbol{\xi}_t, \ \ \ \boldsymbol{\xi}_t \sim N(\boldsymbol{0}, \Phi)
両辺に$\boldsymbol{y}_{t+i}^T(i=0,1,\cdots,k)$をかけて、期待値をとります。
\begin{align}
E(\boldsymbol{y}_t \boldsymbol{y}_t^T)\ \ \ \ &= B_1 E(\boldsymbol{y}_{t+1} \boldsymbol{y}_t^T) \ \ \ + B_2 E(\boldsymbol{y}_{t+2} \boldsymbol{y}_t^T) \ \ + \cdots + B_k E(\boldsymbol{y}_{t+k} \boldsymbol{y}_t^T)\ \ \ + E(\epsilon_t \boldsymbol{y}_t^T)\\
E(\boldsymbol{y}_t \boldsymbol{y}_{t+1}^T) \ &= B_1 E(\boldsymbol{y}_{t+1} \boldsymbol{y}_{t+1}^T) + B_2 E(\boldsymbol{y}_{t+2} \boldsymbol{y}_{t+1}^T) + \cdots + B_k E(\boldsymbol{y}_{t+k} \boldsymbol{y}_{t+1}^T) \ + E(\epsilon_t \boldsymbol{y}_{t+1}^T)\\
\vdots\\
E(\boldsymbol{y}_t \boldsymbol{y}_{t+k}^T) &= B_1 E(\boldsymbol{y}_{t+1} \boldsymbol{y}_{t+k}^T) + B_2 E(\boldsymbol{y}_{t+2} \boldsymbol{y}_{t+k}^T) + \cdots + B_k E(\boldsymbol{y}_{t+k} \boldsymbol{y}_{t+k}^T) + E(\epsilon_t \boldsymbol{y}_{t+k}^T)
\end{align}
となります。これをベクトル・行列化すると
\begin{bmatrix}
C_0 & C_1 & C_2 & & C_k \\
C_{-1} & C_0 & C_1 & \cdots & C_{k-1} \\
C_{-2} & C_{-1} & C_0 & & C_{k-2} \\
& \vdots & &\ddots& \\
C_{-k} & C_{-(k-1)} & C_{-(k-2)} & & C_0
\end{bmatrix}
*
\begin{bmatrix}
E \\- B_1\\- B_2\\\vdots\\- B_k
\end{bmatrix}
=
\begin{bmatrix}
\Phi \\O\\\vdots\\O
\end{bmatrix}
となります。
レビンソン・ダービン法で連立方程式を解く
ここでは、例として3次のVARモデルのパラメータを求めることとします。次の手順で行います。
- [step0] 0次のVARモデルのパラメータを求める
- [step1] 0次の結果を利用して、1次のVARモデルのパラメータを求める
- [step2] 1次の結果を利用して、2次のVARモデルのパラメータを求める
- [step3] 2次の結果を利用して、3次のVARモデルのパラメータを求める
後で分かりますが、逆向きのVARモデルのパラメータを求めることで、通常のVARモデルのパラメータも面白いように求めることができます。それぞれのstepで解くべき式は次の通りです。
-
[step0] 0次VARモデルで解くべき式
- 通常の向き
C_0 = \Sigma_{(0)}\\ - 逆向き
C_0 = \Phi_{(0)}
- 通常の向き
-
[step1] 1次VARモデルで解くべき式
-
通常の向き
\begin{bmatrix} C_0 & C_{-1} \\ C_1 & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(1)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(1)} \\ O \end{bmatrix} -
逆向き
\begin{bmatrix} C_0 & C_1 \\ C_{-1} & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(1)} \end{bmatrix} = \begin{bmatrix} \Phi_{(1)} \\ O \end{bmatrix}
-
-
[step2] 2次VARモデルで解くべき式
-
通常の向き
\begin{bmatrix} C_0 & C_{-1} & C_{-2} \\ C_1 & C_0 & C_{-1} \\ C_2 & C_1 & C_0 \\ \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(2)}\\ - A_2^{(2)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(2)} \\ O\\ O \end{bmatrix}\\ -
逆向き
\begin{bmatrix} C_0 & C_1 & C_2 \\ C_{-1} & C_0 & C_1 \\ C_{-2} & C_{-1} & C_0 \\ \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(2)}\\ - B_2^{(2)} \end{bmatrix} = \begin{bmatrix} \Phi_{(2)} \\ O\\ O \end{bmatrix}
-
-
[step3] 3次VARモデルで解くべき式
-
通常の向き
\begin{bmatrix} C_0 & C_{-1} & C_{-2} & C_{-3} \\ C_1 & C_0 & C_{-1} & C_{-2} \\ C_2 & C_1 & C_0 & C_{-1} \\ C_3 & C_2 & C_1 & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(3)}\\ - A_2^{(3)}\\ - A_3^{(3)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(3)} \\ O\\ O\\ O \end{bmatrix} -
逆向き
\begin{bmatrix} C_0 & C_1 & C_2 & C_3 \\ C_{-1} & C_0 & C_1 & C_2 \\ C_{-2} & C_{-1} & C_0 & C_1 \\ C_{-3} & C_{-2} & C_{-1} & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(3)}\\ - B_2^{(3)}\\ - B_3^{(3)} \end{bmatrix} = \begin{bmatrix} \Phi_{(3)} \\ O\\ O\\ O \end{bmatrix}
-
[step0] 0次のARモデルのパラメータを求める
解くべき式は
C_0 = \Sigma_{(0)}\\ C_0 = \Phi_{(0)}
でした。
未知数$\Sigma_{(0)}, \Phi_{(0)}$は次のようにして求めることができます（左辺右辺を入れ替えるだけ）。
\Sigma_{(0)} = C_0\ \ \ \ \ [0.1]\\
\Phi_{(0)} = C_0\ \ \ \ \ [0.2]
[step1] 0次の結果を利用して、1次のVARモデルのパラメータを求める
解くべき式は、
\begin{bmatrix} C_0 & C_{-1} \\ C_1 & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(1)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(1)} \\ O \end{bmatrix}\\ \begin{bmatrix} C_0 & C_1 \\ C_{-1} & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(1)} \end{bmatrix} = \begin{bmatrix} \Phi_{(1)} \\ O \end{bmatrix}
でした。
1行目に着目すると
C_0 - A_1^{(1)} C_{-1} = \Sigma_{(1)} \ \ \ \ \ [1.1]\\
C_0 - B_1^{(1)} C_1 = \Phi_{(1)} \ \ \ \ \ [1.2]
2行目に着目すると
A_1^{(1)} C_0 = C_1 \ \ \ \ \ [1.3]\\
B_1^{(1)} C_0 = C_{-1} \ \ \ \ \ [1.4]
式[0.1][0.2]の結果を使うと、次のようにまとめることができます。
\begin{align}
A_1^{(1)} &= C_1 \Phi_{(0)}^{-1}\\
B_1^{(1)} &= C_{-1} \Sigma_{(0)}^{-1}\\
\Sigma_{(1)} &= \Sigma_{(0)} - A_1^{(1)}B_1^{(1)} \Sigma_{(0)}\\
\Phi_{(1)} &= \Phi_{(0)} - B_1^{(1)}A_1^{(1)} \Phi_{(0)}
\end{align}
[step2] 1次の結果を利用して、2次のVARモデルのパラメータを求める
解くべき式は、
\begin{bmatrix} C_0 & C_{-1} & C_{-2} \\ C_1 & C_0 & C_{-1} \\ C_2 & C_1 & C_0 \\ \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(2)}\\ - A_2^{(2)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(2)} \\ O\\ O \end{bmatrix}\\ \begin{bmatrix} C_0 & C_1 & C_2 \\ C_{-1} & C_0 & C_1 \\ C_{-2} & C_{-1} & C_0 \\ \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(2)}\\ - B_2^{(2)} \end{bmatrix} = \begin{bmatrix} \Phi_{(2)} \\ O\\ O \end{bmatrix}
でした。
解くべき式の1行目に着目すると
C_0 - A_1^{(2)} C_{-1} - A_2^{(2)} C_{-2} = \Sigma_{(2)} \ \ \ \ \ [2.1]\\
C_0 - B_1^{(2)} C_1 - B_2^{(2)} C_2 = \Phi_{(2)} \ \ \ \ \ [2.2]
解くべき式の2,3行目に着目すると、
\begin{bmatrix}
C_0 & C_{-1} \\
C_1 & C_0
\end{bmatrix}
*
\begin{bmatrix}
A_1^{(2)}\\A_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2
\end{bmatrix}\ \ \ \ \ [2.3]\\
\begin{bmatrix}
C_0 & C_1 \\
C_{-1} & C_0
\end{bmatrix}
*
\begin{bmatrix}
B_1^{(2)}\\
B_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
C_{-1}\\
C_{-2}
\end{bmatrix}\ \ \ \ \ [2.4]
式[1.3][1.4]の結果を用いて、解くべき式の2行目を変形します。
\begin{align}
& A_1^{(2)} C_0 = C_1 - A_2^{(2)} C_{-1}\\
\Leftrightarrow\ & A_1^{(2)} C_0 = A_1^{(1)} C_0 - A_2^{(2)} B_1^{(1)} C_0 \\
\Leftrightarrow\ & A_1^{(2)} = A_1^{(1)} - A_2^{(2)} B_1^{(1)}
\\
\\
& B_1^{(2)} C_0 = C_{-1} - B_2^{(2)} C_1\\
\Leftrightarrow\ & B_1^{(2)} C_0 = B_1^{(1)} C_0 - B_2^{(2)} A_1^{(1)} C_0 \\
\Leftrightarrow\ & B_1^{(2)} = B_1^{(1)} - B_2^{(2)} A_1^{(1)}
\end{align}
つまり、$A_1^{(2)}, B_1^{(2)}$ は $A_2^{(2)}, B_2^{(2)}$ から求めることができます。
この結果と式[1.1]~[1.4]を使って、解くべき式を変形していくことで、$A_2^{(2)}, B_2^{(2)}$を求めます。ついでに$\Sigma_{(2)}, \Phi_{(2)}$も求めます。
まず、通常の向きに関する解くべき式を変形していきます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} \\
C_1 & C_0 & C_{-1} \\
C_2 & C_1 & C_0 \\
\end{bmatrix}*
\begin{bmatrix}
E \\- A_1^{(2)}\\- A_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(2)} \\O\\O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} \\
C_1 & C_0 & C_{-1} \\
C_2 & C_1 & C_0 \\
\end{bmatrix}*
\begin{bmatrix}
E \\ - A_1^{(1)} + A_2^{(2)} B_1^{(1)} \\ - A_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(2)} \\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} \\
C_1 & C_0 & C_{-1} \\
C_2 & C_1 & C_0 \\
\end{bmatrix}
*
\begin{bmatrix}
\begin{bmatrix}
E \\ - A_1^{(1)}\\ O
\end{bmatrix}
- A_2^{(2)}
\begin{bmatrix}
O\\ - B_1^{(1)}\\ E
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(2)} \\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ & \begin{bmatrix}
\Sigma_{(1)} \\ O\\C_2 - A_1^{(1)} C_1
\end{bmatrix}
- A_2^{(2)}
\begin{bmatrix}
C_{-2} - B_1^{(1)} C_{-1} \\O\\\Phi_{(1)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(2)} \\O\\O
\end{bmatrix}
\end{align}
1番下に着目すると
A_2^{(2)} = (C_2 - A_1^{(1)} C_1 )\Phi_{(1)}^{-1}\ \ \ \ \ [2.5]
1番上に着目すると
\Sigma_{(2)} = \Sigma_{(1)} - A_2^{(2)} (C_{-2} - B_1^{(1)} C_{-1} )\ \ \ \ \ [2.6]
となります。
逆向きのVARモデルに関する解くべき式も同じように展開していきます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_{-1} & C_0 & C_1 \\
C_{-2} & C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\- B_1^{(2)}\\- B_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(2)} \\O\\O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_{-1} & C_0 & C_1 \\
C_{-2} & C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\ - B_1^{(1)} + B_2^{(2)} A_1^{(1)} \\ - B_2^{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(2)} \\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_{-1} & C_0 & C_1 \\
C_{-2} & C_{-1} & C_0
\end{bmatrix}
*
\begin{bmatrix}
\begin{bmatrix}
E \\ - B_1^{(1)}\\ O
\end{bmatrix}
- B_2^{(2)}
\begin{bmatrix}
0 \\ - A_1^{(1)}\\ E
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(2)} \\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ & \begin{bmatrix}
\Phi_{(1)} \\ O\\C_{-2} - B_1^{(1)} C_{-1}
\end{bmatrix}
- B_2^{(2)}
\begin{bmatrix}
C_2 - A_1^{(1)} C_1\\O\\\Sigma_{(1)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(2)} \\O\\O
\end{bmatrix}
\end{align}
1番下に着目すると
B_2^{(2)} = (C_{-2} - B^{(1)} C_{-1} )\Sigma_{(1)}^{-1}\ \ \ \ \ [2.7]
1番上に着目すると
\Phi_{(2)} = \Phi_{(1)} - B_2^{(2)} (C_2 - A_1^{(1)} C_1 )\ \ \ \ \ [2.8]
となります。
式[2.5]～[2.8]より、$\Sigma_{(2)},\Phi_{(2)}$は次のように表現できます。
\begin{align}
\Sigma_{(2)} &= \Sigma_{(1)} - {A_2^{(2)}} B_2^{(2)} \Sigma_{(1)} \\
\Phi_{(2)}&= \Phi_{(1)} - {B_2^{(2)}} A_2^{(2)} \Phi_{(1)}
\end{align}
[step3] 2次の結果を利用して、3次のVARモデルのパラメータを求める
解くべき式は、
\begin{bmatrix} C_0 & C_{-1} & C_{-2} & C_{-3} \\ C_1 & C_0 & C_{-1} & C_{-2} \\ C_2 & C_1 & C_0 & C_{-1} \\ C_3 & C_2 & C_1 & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - A_1^{(3)}\\ - A_2^{(3)}\\ - A_3^{(3)} \end{bmatrix} = \begin{bmatrix} \Sigma_{(3)} \\ O\\ O\\ O \end{bmatrix}\begin{bmatrix} C_0 & C_1 & C_2 & C_3 \\ C_{-1} & C_0 & C_1 & C_2 \\ C_{-2} & C_{-1} & C_0 & C_1 \\ C_{-3} & C_{-2} & C_{-1} & C_0 \end{bmatrix} * \begin{bmatrix} E \\ - B_1^{(3)}\\ - B_2^{(3)}\\ - B_3^{(3)} \end{bmatrix} = \begin{bmatrix} \Phi_{(3)} \\ O\\ O\\ O \end{bmatrix}
でした。
解くべき式の1行目に着目すると
C_0 - A_1^{(3)} C_{-1} - A_2^{(3)} C_{-2} - A_3^{(3)} C_{-3} = \Sigma_{(3)} \ \ \ \ \ [3.1]\\
C_0 - B_1^{(3)} C_1 - B_2^{(3)} C_2 - B_3^{(3)} C_3 = \Phi_{(3)} \ \ \ \ \ [3.2]
解くべき式の2,3行目に着目すると、
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} \\
C_1 & C_0 & C_{-1} \\
C_2 & C_1 & C_0
\end{bmatrix}
*
\begin{bmatrix}
A_1^{(3)}\\A_2^{(3)}\\A_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2\\C_3
\end{bmatrix}\ \ \ \ \ [3.3]\\
\begin{bmatrix}
C_0 & C_1 & C_2 \\
C_{-1} & C_0 & C_1 \\
C_{-2} & C_{-1} & C_0
\end{bmatrix}
*
\begin{bmatrix}
B_1^{(3)}\\B_2^{(3)}\\B_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_{-1}\\
C_{-2}\\
C_{-3}
\end{bmatrix}\ \ \ \ \ [3.4]
式[2.3][2.4]の結果を用いて、解くべき式の2,3行目を変形します。
\begin{align}
&\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
A_1^{(3)}\\A_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_1\\C_2\\
\end{bmatrix}
-A_3^{(3)}
\begin{bmatrix}
C_{-2}\\C_{-1}\\
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
A_1^{(3)}\\A_2^{(3)}
\end{bmatrix}
=\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
A_1^{(2)}\\A_2^{(2)}
\end{bmatrix}
-A_3^{(3)}
\begin{bmatrix}
C_0 & C_{-1}\\
C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
B_2^{(2)}\\B_1^{(2)}
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
A_1^{(3)}\\A_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
A_1^{(2)}\\A_2^{(2)}
\end{bmatrix}
-A_3^{(3)}
\begin{bmatrix}
B_2^{(2)}\\B_1^{(2)}
\end{bmatrix}
\end{align}
\begin{align}
&\begin{bmatrix}
C_0 & C_1\\
C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
B_1^{(3)}\\B_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
C_{-1}\\C_{-2}\\
\end{bmatrix}
-B_3^{(3)}
\begin{bmatrix}
C_2\\C_1\\
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1\\
C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
B_1^{(3)}\\B_2^{(3)}
\end{bmatrix}
=\begin{bmatrix}
C_0 & C_1\\
C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
B_1^{(2)}\\B_2^{(2)}
\end{bmatrix}
-B_3^{(3)}
\begin{bmatrix}
C_0 & C_1\\
C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
A_2^{(2)}\\A_1^{(2)}
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
B_1^{(3)}\\B_2^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
B_1^{(2)}\\B_2^{(2)}
\end{bmatrix}
-B_3^{(3)}
\begin{bmatrix}
A_2^{(2)}\\A_1^{(2)}
\end{bmatrix}
\end{align}
つまり、$A_1^{(3)}, A_2^{(3)},B_1^{(3)},B_2^{(3)}$ は $A_3^{(3)}, B_3^{(3)}$ から求めることができます。
この結果と式[2.1]~[2.4]を使って、解くべき式を変形していくことで、$A_2^{(2)}, B_2^{(2)}$を求めます。ついでに$\Sigma_{(2)}, \Phi_{(2)}$も求めます。
まず、通常の向きに関する解くべき式を変形していきます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} & C_{-3} \\
C_1 & C_0 & C_{-1} & C_{-2} \\
C_2 & C_1 & C_0 & C_{-1}\\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\- A_1^{(3)}\\- A_2^{(3)}\\- A_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(3)} \\O\\O\\O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} & C_{-3} \\
C_1 & C_0 & C_{-1} & C_{-2} \\
C_2 & C_1 & C_0 & C_{-1}\\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\ - A_1^{(2)} + A_3^{(3)} B_2^{(2)} \\- A_2^{(2)} + A_3^{(3)} B_1^{(2)} \\ - A_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(3)} \\ O\\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_{-1} & C_{-2} & C_{-3} \\
C_1 & C_0 & C_{-1} & C_{-2} \\
C_2 & C_1 & C_0 & C_{-1}\\
C_3 & C_2 & C_1 & C_0
\end{bmatrix}
*
\begin{bmatrix}
\begin{bmatrix}
E \\ - A_1^{(2)}\\ - A_2^{(2)}\\ O
\end{bmatrix}
- A_3^{(3)}
\begin{bmatrix}
O \\ - B_2^{(2)}\\ - B_1^{(2)}\\ E
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(3)} \\ O\\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ & \begin{bmatrix}
\Sigma_{(2)} \\ O\\O\\C_3 - A_1^{(2)} C_2 - A_2^{(2)} C_1
\end{bmatrix}
- A_3^{(3)}
\begin{bmatrix}
C_{-3} - B_1^{(2)} C_{-2} - B_2^{(2)} C_{-1} \\O\\O\\\Phi_{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Sigma_{(3)} \\O\\O\\O
\end{bmatrix}
\end{align}
1番下に着目すると
A_3^{(3)} = (C_3 - A_1^{(2)} C_2 - A_2^{(2)} C_1)\Phi_{(2)}^{-1}\ \ \ \ \ [3.5]
1番上に着目すると
\Sigma_{(3)} = \Sigma_{(2)} - A_3^{(3)} (C_{-3} - B_1^{(2)}C_{-2} - B_2^{(2)} C_{-1} ) \ \ \ \ \ [3.6]
となります。
逆向きのVARモデルに関する解くべき式も同じように展開していきます。
\begin{align}
&
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_{-1} & C_0 & C_1 & C_2 \\
C_{-2} & C_{-1} & C_0 & C_1\\
C_{-3} & C_{-2} & C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\- B_1^{(3)}\\- B_2^{(3)}\\- B_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(3)} \\O\\O\\O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_{-1} & C_0 & C_1 & C_2 \\
C_{-2} & C_{-1} & C_0 & C_1\\
C_{-3} & C_{-2} & C_{-1} & C_0
\end{bmatrix}*
\begin{bmatrix}
E \\ - B_1^{(2)} + B_3^{(3)} A_2^{(2)} \\- B_2^{(2)} + B_3^{(3)} A_1^{(2)} \\ - B_3^{(3)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(3)} \\ O\\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ &
\begin{bmatrix}
C_0 & C_1 & C_2 & C_3 \\
C_{-1} & C_0 & C_1 & C_2 \\
C_{-2} & C_{-1} & C_0 & C_1\\
C_{-3} & C_{-2} & C_{-1} & C_0
\end{bmatrix}
*
\begin{bmatrix}
\begin{bmatrix}
E \\ - B_1^{(2)}\\ - B_2^{(2)}\\ O
\end{bmatrix}
- B_3^{(3)}
\begin{bmatrix}
O\\ - A_2^{(2)}\\ - A_1^{(2)}\\ E
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(3)} \\ O\\ O\\ O
\end{bmatrix}\\
\Leftrightarrow\ & \begin{bmatrix}
\Phi_{(2)} \\ O\\O\\C_{-3} - B_1^{(2)} C_{-2} - B_2^{(2)} C_{-1}
\end{bmatrix}
- B_3^{(3)}
\begin{bmatrix}
C_3 - A_1^{(2)} C_2 - A_2^{(2)} C_1 \\O\\O\\\Sigma_{(2)}
\end{bmatrix}
=
\begin{bmatrix}
\Phi_{(3)} \\O\\O\\O
\end{bmatrix}
\end{align}
1番下に着目すると
B_3^{(3)} = (C_{-3} - B_1^{(2)} C_{-2} - B_2^{(2)} C_{-1})\Sigma_{(2)}^{-1}\ \ \ \ \ [3.7]
1番上に着目すると
\Phi_{(3)} = \Phi_{(2)} - B_3^{(3)} (C_3 - A_1^{(2)} C_2 - A_2^{(2)} C_1) \ \ \ \ \ [3.8]
となります。
式[3.5]～[3.8]より、$\Sigma_{(3)},\Phi_{(3)}$は次のように表現できます。
\begin{align}
\Sigma_{(3)} &= \Sigma_{(2)} - {A_3^{(3)}} B_3^{(3)} \Sigma_{(2)}\\
\Phi_{(3)} &= \Phi_{(2)} - {B_3^{(3)}} A_3^{(3)} \Phi_{(2)}
\end{align}
まとめ）K次のARモデルのパラメータを求めるには
4次以上も3次と同じ方法で求めることができます。$K$次のVARモデルのパラメータを求める方法は次の通りです。
- 0,1次を求める
\begin{align}
\Sigma_{(0)} &= C_0\\
\Phi_{(0)} &= C_0\\
A_1^{(1)} &= C_1 \Phi_{(0)}^{-1}\\
B_1^{(1)} &= C_{-1} \Sigma_{(0)}^{-1}\\
\Sigma_{(1)} &= \Sigma_{(0)} - A_1^{(1)}B_1^{(1)} \Sigma_{(0)}\\
\Phi_{(1)} &= \Phi_{(0)} - B_1^{(1)}A_1^{(1)} \Phi_{(0)}
\end{align}
- $k=2,\cdots,K$を次の式に順に代入し、逐次的に$K$次まで求める
\begin{align}
A_k^{(k)} &= \Bigl(C_k - \sum_{i=1}^{k-1} A_i^{(k-1)}C_{k-i}\Bigr)\Phi_{(k-1)}^{-1}
\\
B_k^{(k)} &=\Bigl(C_{-k} - \sum_{i=1}^{k-1} B_i^{(k-1)}C_{-(k-i)}\Bigr)\Sigma_{(k-1)}^{-1}
\\
{A_i^{(k)}} &= {A_i^{(k-1)}} - {A_k^{(k)}}{B_{k-i}^{(k-1)}} \ \ \ (i=1,\cdots,k-1)\\
{B_i^{(k)}} &= {B_i^{(k-1)}} - {B_k^{(k)}}{A_{k-i}^{(k-1)}} \ \ \ (i=1,\cdots,k-1)\\
\Sigma_{(k)} &= \Sigma_{(k-1)} - {A_k^{(k)}} B_k^{(k)} \Sigma_{(k-1)}\\
\Phi_{(k)} &= \Phi_{(k-1)} - {B_k^{(k)}} A_k^{(k)} \Phi_{(k-1)}
\end{align}
図にすると次のようなイメージです。
実は、下式が成り立つので、これを利用して$B_k^{(k)}$の計算することも可能なのですが、図の美しさが崩れるので、今回は不採用としています。
\Bigl(C_{-k} - \sum_{i=1}^{k-1} B_i^{(k-1)}C_{-(k-i)}\Bigr) = \Bigl(C_k - \sum_{i=1}^{k-1} A_i^{(k-1)}C_{k-i}\Bigr)^T
Pythonのコードの例は次の通りです。うまくパラメータを推定できていることが分かります。
（ただし、わかりやすさ重視のコードのため、遅いと思います。for文をベクトル演算に置き換えると高速になります。）
import numpy as np
import scipy
np.random.seed(42)
#######################
# データ作成
#######################
n = 100000 # データ数
m = 2 # 特徴量数
a = np.array([-1, 1]) # バイアス項
A = [
None,
np.array([[1/2, -1/3], [1/4, 1/5]]), # 1次の係数はA[0]ではなくA[1]に入れることとする
np.array([[-1/4, -1/5], [1/8, 1/6]]),
np.array([[-1/3, 1/3], [-1/5, 1/3]])
]
k = 3 # 次数
Sigma = np.eye(m) # 誤差項の分散共分散行列
epsilons = scipy.stats.multivariate_normal.rvs(mean=np.zeros([m]), cov=Sigma, size=n)
X = np.zeros([n, m])
for i in range(k):
X[i, :] = epsilons[i, :] # X[0] = 乱数, X[1] = 乱数, X[2] = 乱数
for i in range(k, n):
X[i] = a
for j in range(1, k + 1):
X[i] += np.dot(A[j], X[i - j])
X[i] += epsilons[i]
#############################################
# レビンソン・ダービン法でa, A, Sigmaを推定
#############################################
# 平均を0に
mu = X.mean(axis=0)
X = X - mu
# C[0], C[1], C[2], C[3]を導出
C = []
for i in range(k + 1):
s = 0
for j in range(n - i):
s += np.dot(X[[j + i], :].T, X[[j], :])
s /= (n - i)
C.append(s)
# 3次のVARモデルの係数をレビンソン・ダービン法で導出
for _k in range(k + 1):
if _k == 0: # 0次
Sigma = C[0].copy()
Phi = C[0].copy()
a = mu
A = [None]
B = [None]
else: # 1次、2次、3次
A_old = A.copy()
B_old = B.copy()
A = [None] * (_k + 1)
B = [None] * (_k + 1)
s = 0
for i in range(1, _k):
s += np.dot(A_old[i], C[_k - i])
A[_k] = np.dot(C[_k] - s, scipy.linalg.inv(Phi))
s = 0
for i in range(1, _k):
s += np.dot(B_old[i], C[_k - i].T)
B[_k] = np.dot(C[_k].T - s, scipy.linalg.inv(Sigma))
#B[_k] = np.dot((C[_k] - s).T, scipy.linalg.inv(Sigma))
for i in range(1, _k):
A[i] = A_old[i] - np.dot(A[_k], B_old[_k - i])
B[i] = B_old[i] - np.dot(B[_k], A_old[_k - i])
s = np.eye(m)
for i in range(1, _k + 1):
s -= A[i]
a = np.dot(s, mu)
Sigma -= np.dot(np.dot(A[_k], B[_k]), Sigma)
Phi -= np.dot(np.dot(B[_k], A[_k]), Phi)
print(a)
# [-1.00380545 0.98742596] 真の値 [-1.0, 1.0]
print(A)
# [None,
# array([[ 0.49859278, -0.32981643], 真の値 [[ 0.500, -0.333],
# [ 0.24812199, 0.20107688]]) [ 0.250, 0.200]]
# array([[-0.24964478, -0.19875606], [[-0.250, -0.200],
# [ 0.12833023, 0.16483817]]) [-0.125, 0.166]]
# array([[-0.33257546, 0.3322076 ], [[-0.333, 0.333],
# [-0.20459696, 0.33583133]]) [-0.200, 0.333]]
# ]
print(Sigma)
# [[1.00348475 0.00197257] 真の値 [[1.0, 0.0],
# [0.00197257 0.99619334]] [0.0, 1.0]]