Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@yutera12(tera)

ARモデル、VARモデルのパラメータをレビンソン・ダービン法で推定【Pythonコード付き】

Last updated at Posted at 2023-03-27

はじめに

この記事では、自己回帰モデル（ARモデル）や、ベクトル自己回帰モデル（VARモデル）のパラメータをレビンソン・ダービン（Levinson-Durbin）法という方法で、推定する方法を説明します。

この2つの資料を参考にしました。

本記事の特徴

  • 上記の参考資料では、「$i$次のパラメータを使った$i+1$次のパラメータの導出」の方法が書かれていましたが、本記事では、次のように具体的な例で示しました。
    • 0次のパラメータを導出
    • 0次のパラメータを使って1次のパラメータを導出
    • 1次のパラメータを使って2次のパラメータを導出
    • 2次のパラメータを使って3次のパラメータを導出
  • 参考資料では、VARモデルのパラメータ推定方法（すなわち、レビンソン・ダービン法の多次元化）の導出について触れられていなかったので、そこを詳しく記述しました。
  • カラーの図を使ってアルゴリズムの流れがわかるようにしました。VARモデルのパラメータを推定する流れを図で表現したものがこちらです。美しい構造をしています。

var.png

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}

図にすると、次のようなイメージです。

ar.png

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}

図にすると次のようなイメージです。

var.png

実は、下式が成り立つので、これを利用して$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]]
2
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?