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

# はじめに

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

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

* [京都大学講義資料　Levinson-Durbin のアルゴリズム](http://www.kumst.kyoto-u.ac.jp/kougi/time_series/Appendix1.pdf)
* [東京大学講義資料　時系列解析（６）](http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E8%A7%A3%E6%9E%90%EF%BC%88%EF%BC%96%EF%BC%89_2019.pdf)

# 本記事の特徴

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

![var.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/b757e3e1-969f-b321-c38b-c4773b7c9680.png)

# ARモデルの推定

自己回帰モデル（ARモデル）

```math
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$を推定します。

:::note info
一般的にはARモデルの誤差項$\epsilon$の分散（1次元正規分布の分散）は$\sigma^2$と記述されますが、ここでは、表記を単純化するためと、VARモデルと記載を統一するために、$\sigma$とします。
:::

## 自己共分散から連立方程式を作る

まず、表記を簡単にするために平均が0になるように調整します。$x_t$の平均$\mu$は、観測データから推定値を求めることができるので、これを既知とします。新しい時系列データ$y_t$を次の式で定義します。

```math
y_t = x_t - \mu
```


こうすると、$E(y_t)=0$となるので、式[A]は次のように変形することができます。

```math
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$を求めていきます。

:::note info

$a_0$は、$a_1, a_2, \cdots, a_k$を求めた後に次の方法で求めることができます。

```math
a_0 = (1-a_1-a_2-\cdots -a_k)\mu
```

なぜならば、
```math
\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)$をかけて、期待値をとります。

```math
\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})$とすると、

```math
\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}
```

となります。これをベクトル・行列化すると

```math
\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モデルで解くべき式


    ```math
    C_0 = \sigma_{(0)}
    ```

* [step1] 1次ARモデルで解くべき式

    ```math
    \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モデルで解くべき式

    ```math
    \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モデルで解くべき式

    ```math
    \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モデルのパラメータを求める

解くべき式は

> ```math
> C_0 = \sigma_{(0)}
> ```

でした。

未知数$\sigma_{(0)}$は次のようにして求めることができます（左辺と右辺を入れ替えるだけです）。
```math
\sigma_{(0)} = C_0  \ \ \ \ [0.1]
```

### [step1] 0次の結果を利用して、1次のARモデルのパラメータを求める

解くべき式は、

> ```math
> \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行目に着目すると

```math
C_0 - a_1^{(1)} C_1 = \sigma_{(1)} \ \ \ \ \ [1.1]
```

2行目に着目すると
```math
a_1^{(1)} C_0 = C_1  \ \ \ \ [1.2]
```

式[0.1]の結果を使うと、次のようにまとめることができます。
```math
\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モデルのパラメータを求める

解くべき式は、

> ```math
> \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行目に着目すると、
```math
C_0 - a_1^{(2)}C_1 - a_2^{(2)} C_2 = \sigma_{(2)} \ \ \ \ \ [2.1]
```

解くべき式の2,3行目に着目すると、

```math
\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行目を変形します。

```math
\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)}$も求めます。

```math
\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番下の行に着目すると

```math
a_2^{(2)} = 
(C_2 - a_1^{(1)} C_1) \sigma_{(1)}^{-1}\\
```

1番上の行に着目すると

```math
\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モデルのパラメータを求める

解くべき式は

> ```math
> \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行目に着目すると、
```math
C_0 - a_1^{(3)}C_1 - a_2^{(3)} C_2 - a_3^{(3)} C_3 = \sigma_{(3)} \ \ \ \ \ [3.1]
```
解くべき式の2,3,4行目に着目すると、

```math
\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行目を変形します。

```math
\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)}$も求めます。
```math
\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番下の行に着目すると

```math
a_3^{(3)} = (C_3 - a_1^{(2)} C_2 - a_2^{(2)} C_1) \sigma_{(2)}^{-1}
```

1番上の行に着目すると
```math
\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次を求める
```math
\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$次まで求める

```math
\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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/58f0f97e-0dc4-4bc7-0bae-f4f9d9f33311.png)

Pythonのコードの例は次の通りです。うまくパラメータを推定できていることが分かります。
（ただし、わかりやすさ重視のコードのため、遅いと思います。for文をベクトル演算に置き換えると高速になります。）

```Python
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モデル）

```math
\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$を次の式で定義します。

```math
\boldsymbol{y_t} = \boldsymbol{x}_t - \boldsymbol{\mu}
```

こうすると、$E(\boldsymbol{y}_t)=\boldsymbol{0}$となるので、式[B]は次のように変形することができます。


```math
\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$を求めていきます。

:::note info
$\boldsymbol{a}_0$は、$A_1, A_2, \cdots, A_k$を求めた後に次の方法で求めることができます。

```math
\boldsymbol{a}_0 = (E-A_1-A_2-\cdots -A_k)\boldsymbol{\mu}
```

なぜならば、

```math
\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モデルは次の式で与えられます。

```math
\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)$をかけて、期待値をとります。


```math
\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)$とすると、

```math
\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}
```

これをベクトル・行列化すると

```math
\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$次元の行列であることに注意してください。

:::note info
ここで、独自の記号である$*$について説明します。

通常の行列の積を使った場合には下のようになります。

```math
\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}
```

$*$を使った場合は下のようになるとします。

```math
\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}
```
:::

:::note warn
1次元の時とは違い$C_i=C_{-i}$は成立しないので、この形にすることは**できません**。

```math
\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モデルは次の式で与えられます。

```math
\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)$をかけて、期待値をとります。


```math
\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}
```

となります。これをベクトル・行列化すると
```math
\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モデルで解くべき式
    * 通常の向き
        ```math
        C_0 = \Sigma_{(0)}\\
        ```
    * 逆向き
        ```math
        C_0 = \Phi_{(0)}
        ```

* [step1] 1次VARモデルで解くべき式
    * 通常の向き

        ```math
        \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}
        ```
    * 逆向き
        ```math
        \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モデルで解くべき式
    * 通常の向き

        ```math
        \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}\\
        ```
    * 逆向き
        ```math
        \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モデルで解くべき式
    * 通常の向き

        ```math
        \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}
        ```
    * 逆向き
        ```math
        \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モデルのパラメータを求める

解くべき式は

> ```math
> C_0 = \Sigma_{(0)}\\
> C_0 = \Phi_{(0)}
> ```

でした。

未知数$\Sigma_{(0)}, \Phi_{(0)}$は次のようにして求めることができます（左辺右辺を入れ替えるだけ）。

```math
\Sigma_{(0)} = C_0\ \ \ \ \ [0.1]\\
\Phi_{(0)} = C_0\ \ \ \ \ [0.2]
```

### [step1] 0次の結果を利用して、1次のVARモデルのパラメータを求める

解くべき式は、

> ```math
> \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行目に着目すると

```math
C_0 - A_1^{(1)} C_{-1} = \Sigma_{(1)}  \ \ \ \ \ [1.1]\\
C_0 - B_1^{(1)} C_1 = \Phi_{(1)}  \ \ \ \ \ [1.2]
```

2行目に着目すると

```math
A_1^{(1)} C_0  = C_1  \ \ \ \ \ [1.3]\\
B_1^{(1)} C_0  = C_{-1}  \ \ \ \ \  [1.4]
```

式[0.1][0.2]の結果を使うと、次のようにまとめることができます。

```math
\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モデルのパラメータを求める

解くべき式は、

> ```math
> \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行目に着目すると
```math
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行目に着目すると、

```math
\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行目を変形します。

```math
\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)}$も求めます。

まず、通常の向きに関する解くべき式を変形していきます。

```math
\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番下に着目すると
```math
A_2^{(2)} = (C_2 - A_1^{(1)} C_1 )\Phi_{(1)}^{-1}\ \ \ \ \ [2.5]
```

1番上に着目すると
```math
\Sigma_{(2)} = \Sigma_{(1)} - A_2^{(2)} (C_{-2} - B_1^{(1)} C_{-1} )\ \ \ \ \ [2.6]
```

となります。

逆向きのVARモデルに関する解くべき式も同じように展開していきます。

```math
\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番下に着目すると
```math
B_2^{(2)} = (C_{-2} - B^{(1)} C_{-1}  )\Sigma_{(1)}^{-1}\ \ \ \ \ [2.7]
```

1番上に着目すると

```math
\Phi_{(2)} = \Phi_{(1)} - B_2^{(2)} (C_2 - A_1^{(1)} C_1  )\ \ \ \ \ [2.8]
```

となります。

式[2.5]～[2.8]より、$\Sigma_{(2)},\Phi_{(2)}$は次のように表現できます。

```math
\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モデルのパラメータを求める

解くべき式は、
> ```math
> \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}
> ```
> ```math
> \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行目に着目すると

```math
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行目に着目すると、

```math
\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行目を変形します。
```math
\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}
```

```math
\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)}$も求めます。

まず、通常の向きに関する解くべき式を変形していきます。
```math
\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番下に着目すると
```math
A_3^{(3)} = (C_3 - A_1^{(2)} C_2 -  A_2^{(2)} C_1)\Phi_{(2)}^{-1}\ \ \ \ \ [3.5]
```

1番上に着目すると
```math
\Sigma_{(3)} = \Sigma_{(2)} - A_3^{(3)} (C_{-3} -  B_1^{(2)}C_{-2}  - B_2^{(2)} C_{-1} ) \ \ \ \ \ [3.6]
```

となります。


逆向きのVARモデルに関する解くべき式も同じように展開していきます。

```math
\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番下に着目すると
```math
B_3^{(3)} = (C_{-3} -  B_1^{(2)} C_{-2} - B_2^{(2)} C_{-1})\Sigma_{(2)}^{-1}\ \ \ \ \ [3.7]
```

1番上に着目すると
```math
\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)}$は次のように表現できます。

```math
\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次を求める

```math
\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$次まで求める

```math
\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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/571330/b757e3e1-969f-b321-c38b-c4773b7c9680.png)

実は、下式が成り立つので、これを利用して$B_k^{(k)}$の計算することも可能なのですが、図の美しさが崩れるので、今回は不採用としています。

```math
\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文をベクトル演算に置き換えると高速になります。）


```Python
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]]
```