このコードは、2変量の時系列データに対してVARモデルを適用し、変数間の相互依存関係を推定しています。

import numpy as np
from statsmodels.tsa.api import VAR

np.random.seed(42)
n, e = 300, np.random.randn(2, 300)*0.5
y = np.zeros((n, 2))
for t in range(2, n):
    y[t] = [[.7,.2],[.3,.6]] @ y[t-1] + e[:, t]

model = VAR(y).fit(maxlags=5, ic='aic')
print(
    f"ラグ次数: {model.k_ar}\n"
    f"係数行列:\n{model.coefs[0]}\n"
    f"予測:\n{model.forecast(y[-model.k_ar:], 3)}")

簡単に説明します。

ライブラリのインポート

• numpy：数値計算を行うためのライブラリ。ここでは乱数生成や行列演算に使用しています。
• statsmodels.tsa.api.VAR：VARモデルを実装したクラスです。

データの生成

• np.random.randn(2, 300)*0.5：2変数分のホワイトノイズ（標準偏差0.5）を生成しています。
• y[t] = [[.7,.2],[.3,.6]] @ y[t-1] + e[:, t]：VAR(1)過程に従うデータを生成しています。@は行列積の演算子です。この式を展開すると以下のようになります。
  • y_1(t) = 0.7 \cdot y_1(t-1) + 0.2 \cdot y_2(t-1) + e_1(t)
  • y_2(t) = 0.3 \cdot y_1(t-1) + 0.6 \cdot y_2(t-1) + e_2(t)
• つまり、各変数は自分自身の過去の値だけでなく、他方の変数の過去の値にも依存する構造を持っています。

VARモデルの構築

• VAR(y).fit(maxlags=5, ic='aic')：VARモデルを推定します。
  • maxlags=5：最大5ラグまでを探索対象とします。
  • ic='aic'：AICを基準に最適なラグ次数を自動選択します。
• model.k_ar：選択されたラグ次数です。今回はラグ1が選択されており、データ生成のVAR(1)過程と一致しています。
• model.coefs[0]：ラグ1の係数行列です。真の係数[[0.7, 0.2], [0.3, 0.6]]に近い値が推定されています。
• model.forecast()：学習データの末尾から3期先までの予測を行います。

(A) 1つの変数の過去の値のみから将来を予測している → ×

1つの変数の過去の値のみを使うのはAR（自己回帰）モデルです。VARモデルは複数の変数を同時にモデル化し、各変数の予測に他の変数の過去の値も活用します。

(B) 複数の変数の相互依存関係を同時にモデル化している → ○

正解です。VARモデルでは、すべての変数がすべての変数のラグ値に依存する構造を仮定します。係数行列の非対角要素が変数間の影響度を表しており、今回の結果ではy_1からy_2への影響（0.363）とy_2からy_1への影響（0.2162）が推定されています。

(C) 外生変数を用いて目的変数を回帰分析している → ×

VARモデルでは、すべての変数が内生変数として扱われます。外生変数を含むモデルはVARXモデルと呼ばれ、VARの拡張版です。通常の回帰分析のような「説明変数→目的変数」の一方向の関係ではなく、変数間の双方向の関係をモデル化しています。

(D) 変数間の因果関係の方向を直接推定している → ×

VARモデルの係数は変数間の予測に役立つ関係を示しますが、因果関係を直接意味するものではありません。因果関係の方向を検定するには、VARモデルを前提としたグレンジャー因果検定を追加で実施する必要があります。

VAR(p)モデルは、以下の式で定義されます。

$$
\mathbf{y}_t = \mathbf{c} + A_1 \mathbf{y}_{t-1} + A_2 \mathbf{y}_{t-2} + \cdots + A_p \mathbf{y}_{t-p} + \mathbf{e}_t
$$

ここで、\mathbf{y}_tはK次元の変数ベクトル、\mathbf{c}は定数項ベクトル、A_1, \ldots, A_pはK×Kの係数行列、\mathbf{e}_tはホワイトノイズベクトルです。

今回のような2変量VAR(1)モデルの場合、K=2、p=1で、係数行列は1つ（A_1）です。このA_1の各要素が、変数間の1期前の影響度を表します。

$$
A_1 = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}
$$

• a_{11}, a_{22}（対角要素）：各変数の自己依存の強さ
• a_{12}, a_{21}（非対角要素）：変数間の交差的影響の強さ

VARモデルの係数行列から「変数間に予測的な関係があるか」はわかりますが、それが因果関係かどうかは別の問題です。グレンジャー因果検定は、VARモデルを前提として「変数y_2の過去の値を使うことでy_1の予測精度が有意に改善するか」を検定する手法です。

statsmodelsでは、VARモデルの推定後に以下のように実行できます。

model.test_causality('y1', 'y2', kind='f')

この検定では、係数行列の非対角要素がまとめて0であるかどうかをF検定で判定します。ただし、グレンジャー因果はあくまで「予測に役立つかどうか」の検定であり、経済学や自然科学における「真の因果関係」を証明するものではない点に注意が必要です。

VARモデルが安定（定常）であるためには、係数行列の固有値がすべて単位円の内部にある必要があります。今回のデータ生成に使った係数行列[[0.7, 0.2], [0.3, 0.6]]の固有値は約0.87と0.43で、いずれも絶対値が1未満のため安定条件を満たしています。

statsmodelsではmodel.is_stable()で安定性を確認できます。安定でないVARモデルは予測が発散するため、差分処理などで事前に定常化する必要があります。