株価が上がると消費が増えるのでしょうか、それとも消費が増えるから株価が上がるのでしょうか？金利を上げるとインフレ率はどのくらいの期間で低下するのでしょうか？こうした「複数の変数が互いに影響し合う」システムを分析するために開発されたのがVARモデル（Vector AutoRegression: ベクトル自己回帰モデル）です。
通常のARモデルは1つの変数の過去の値だけを使って予測しますが、VARモデルは複数の変数が互いの過去の値から影響を受ける構造を明示的にモデル化します。計量経済学者クリストファー・シムズが1980年に提案し、2011年のノーベル経済学賞の受賞理由の一つとなったこの手法は、現在でもマクロ経済分析の標準ツールです。
VARモデルを理解することで、以下のような場面で活用できます。
- マクロ経済分析 — GDP、失業率、金利、インフレ率など複数のマクロ変数の相互影響を分析する
- グレンジャー因果の検定 — 「変数Aが変数Bの予測に役立つか」を統計的に検定する
- インパルス応答分析 — あるショック（例えば金利の突然の上昇）が他の変数にどのように波及するかを追跡する
本記事の内容
- 単変量ARの限界と多変量への拡張の動機
- VAR(p)モデルの定式化（行列表現）
- 安定性条件
- OLSによるパラメータ推定
- グレンジャー因果テスト
- インパルス応答関数（IRF）
- 予測誤差の分散分解（FEVD）
- Python（statsmodels）での実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
単変量ARモデルの限界
単変量のAR(p)モデルは、1つの変数の過去の値のみを使って将来を予測します。
$$ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + \varepsilon_t $$
しかし、現実の多くのシステムでは、複数の変数が互いに影響し合っています。例えば、気温と電力需要の関係を考えてみましょう。
- 気温が上がると → エアコンの使用が増えて → 電力需要が増える
- 電力需要が増えると → 発電所がフル稼働して → 排熱で気温がわずかに上がる（ヒートアイランド効果）
このような双方向の影響関係は、単変量モデルでは捉えられません。変数 $x_t$ と $y_t$ がある場合、$x_t$ の予測に $y$ の過去の値も使い、同時に $y_t$ の予測にも $x$ の過去の値を使うモデルが必要です。これがVARモデルの出発点です。
VARモデルは、このような多変量間の相互依存を「全ての変数が全ての変数の過去の値に依存する」という対称的な定式化で表現します。次のセクションでは、この定式化を数学的に見ていきましょう。
VAR(p)モデルの定式化
2変量VAR(1)モデル
まず、最も単純な2変量VAR(1)モデルから始めましょう。$K = 2$ 個の変数 $y_{1,t}$ と $y_{2,t}$ の場合、
$$ y_{1,t} = c_1 + \phi_{11} y_{1,t-1} + \phi_{12} y_{2,t-1} + \varepsilon_{1,t} $$
$$ y_{2,t} = c_2 + \phi_{21} y_{1,t-1} + \phi_{22} y_{2,t-1} + \varepsilon_{2,t} $$
ここで $\phi_{12}$ は「$y_2$ の過去の値が $y_1$ に与える影響」を表す係数です。もし $\phi_{12} \neq 0$ ならば、$y_2$ の過去の情報が $y_1$ の予測に役立つことを意味します。
行列表現
上の2式をベクトルと行列を使ってまとめると、
$$ \begin{pmatrix} y_{1,t} \\ y_{2,t} \end{pmatrix} = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} + \begin{pmatrix} \phi_{11} & \phi_{12} \\ \phi_{21} & \phi_{22} \end{pmatrix} \begin{pmatrix} y_{1,t-1} \\ y_{2,t-1} \end{pmatrix} + \begin{pmatrix} \varepsilon_{1,t} \\ \varepsilon_{2,t} \end{pmatrix} $$
ベクトル $\bm{y}_t = (y_{1,t}, y_{2,t})^\top$、$\bm{c} = (c_1, c_2)^\top$、係数行列 $\bm{A}_1$、ノイズベクトル $\bm{\varepsilon}_t$ を定義すると、
$$ \bm{y}_t = \bm{c} + \bm{A}_1 \bm{y}_{t-1} + \bm{\varepsilon}_t $$
一般のVAR(p)モデル
$K$ 個の変数、$p$ 次のラグを持つ一般的なVAR(p)モデルは、次のように定式化されます。
$$ \begin{equation} \bm{y}_t = \bm{c} + \bm{A}_1 \bm{y}_{t-1} + \bm{A}_2 \bm{y}_{t-2} + \cdots + \bm{A}_p \bm{y}_{t-p} + \bm{\varepsilon}_t \end{equation} $$
ここで、
- $\bm{y}_t = (y_{1,t}, y_{2,t}, \dots, y_{K,t})^\top$ は $K \times 1$ のベクトル
- $\bm{c}$ は $K \times 1$ の定数ベクトル
- $\bm{A}_i$ は $K \times K$ の係数行列（$i = 1, 2, \dots, p$）
- $\bm{\varepsilon}_t$ は $K \times 1$ のホワイトノイズベクトル
ノイズベクトル $\bm{\varepsilon}_t$ は以下の性質を満たすと仮定します。
$$ \mathbb{E}[\bm{\varepsilon}_t] = \bm{0}, \quad \mathbb{E}[\bm{\varepsilon}_t \bm{\varepsilon}_t^\top] = \bm{\Sigma}, \quad \mathbb{E}[\bm{\varepsilon}_t \bm{\varepsilon}_s^\top] = \bm{O} \quad (t \neq s) $$
$\bm{\Sigma}$ は $K \times K$ の分散共分散行列で、一般に非対角要素がゼロでないことに注意してください。つまり、同時刻の異なる変数のノイズは相関を持ちうるのです。
パラメータの数
VAR(p)モデルのパラメータ数を数えてみましょう。各 $\bm{A}_i$ は $K \times K$ の行列で $p$ 個あるので $K^2 p$ 個、定数ベクトル $\bm{c}$ が $K$ 個、分散共分散行列 $\bm{\Sigma}$ が $K(K+1)/2$ 個（対称行列なので）。合計で、
$$ K^2 p + K + \frac{K(K+1)}{2} $$
例えば、$K = 5$ 変数、$p = 4$ ラグの場合、$5^2 \times 4 + 5 + 15 = 120$ 個のパラメータが必要です。変数の数やラグの数が増えると急速にパラメータ数が増大するため、サンプルサイズに対してモデルが過大にならないよう注意が必要です。
VARモデルの定式化を理解したところで、次にこのモデルが定常過程を生成するための条件を見ていきましょう。
安定性条件
VAR(1)の場合
VAR(1)モデル $\bm{y}_t = \bm{c} + \bm{A}_1 \bm{y}_{t-1} + \bm{\varepsilon}_t$ が定常であるための条件は、係数行列 $\bm{A}_1$ の全ての固有値の絶対値が1未満であることです。
$$ \begin{equation} |\lambda_i(\bm{A}_1)| < 1, \quad i = 1, 2, \dots, K \end{equation} $$
この条件は、単変量AR(1)の定常性条件 $|\phi| < 1$ の自然な一般化です。固有値の絶対値が1以上だと、ショックの影響が減衰せず、過程が爆発的に振る舞います。
コンパニオン形式
VAR(p)モデルの安定性条件は、コンパニオン行列を用いてVAR(1)の形に帰着させることで統一的に扱えます。
VAR(p)モデルを次のようにVAR(1)形式に変換します。新しい状態ベクトルを
$$ \bm{Y}_t = \begin{pmatrix} \bm{y}_t \\ \bm{y}_{t-1} \\ \vdots \\ \bm{y}_{t-p+1} \end{pmatrix} $$
と定義すると、$Kp \times 1$ のベクトルになります。コンパニオン行列は
$$ \bm{A} = \begin{pmatrix} \bm{A}_1 & \bm{A}_2 & \cdots & \bm{A}_{p-1} & \bm{A}_p \\ \bm{I}_K & \bm{O} & \cdots & \bm{O} & \bm{O} \\ \bm{O} & \bm{I}_K & \cdots & \bm{O} & \bm{O} \\ \vdots & & \ddots & & \vdots \\ \bm{O} & \bm{O} & \cdots & \bm{I}_K & \bm{O} \end{pmatrix} $$
この $Kp \times Kp$ のコンパニオン行列の全ての固有値の絶対値が1未満であれば、VAR(p)過程は定常です。
$$ \begin{equation} |\lambda_i(\bm{A})| < 1, \quad \forall i \end{equation} $$
安定性条件を満たすVARモデルに対しては、パラメータの推定が可能です。次に、最も基本的な推定方法であるOLSについて見ていきましょう。
OLSによるパラメータ推定
各方程式を個別にOLS推定
VARモデルの重要な性質として、各方程式を個別に最小二乗法（OLS）で推定しても一致推定量が得られることが知られています。
VAR(p)モデルの第 $k$ 番目の方程式は、
$$ y_{k,t} = c_k + \sum_{i=1}^{p} \sum_{j=1}^{K} a_{kj}^{(i)} y_{j,t-i} + \varepsilon_{k,t} $$
これは通常の線形回帰モデルの形をしています。説明変数は $\{y_{j,t-i}\}$ で、全ての方程式に共通です。このような構造をSUR（Seemingly Unrelated Regressions）モデルと呼びます。全ての方程式の説明変数が同じ場合、SURのGLS推定量はOLS推定量と一致します。
行列表現での推定
データが $T$ 期分あるとき、各方程式を行列形式で書くと、
$$ \bm{Y} = \bm{X} \bm{B} + \bm{E} $$
ここで $\bm{Y}$ は $(T-p) \times K$ の観測行列、$\bm{X}$ は説明変数の行列（定数項と $p$ 期分のラグ変数）、$\bm{B}$ はパラメータ行列です。OLS推定量は、
$$ \begin{equation} \hat{\bm{B}} = (\bm{X}^\top \bm{X})^{-1} \bm{X}^\top \bm{Y} $$
誤差の分散共分散行列は、
$$ \hat{\bm{\Sigma}} = \frac{1}{T-p-Kp-1} (\bm{Y} – \bm{X}\hat{\bm{B}})^\top (\bm{Y} – \bm{X}\hat{\bm{B}}) $$
で推定されます。
ラグ次数の選択
VAR(p)のラグ次数 $p$ の選択は、情報量基準を用いて行います。主要な基準は以下の通りです。
$$ \text{AIC}(p) = \ln |\hat{\bm{\Sigma}}_p| + \frac{2K^2 p}{T} $$
$$ \text{BIC}(p) = \ln |\hat{\bm{\Sigma}}_p| + \frac{K^2 p \ln T}{T} $$
BICはAICよりも大きなペナルティを課すため、よりパーシモニアスな（少ないラグの）モデルを選択する傾向があります。
パラメータの推定方法を理解したところで、VARモデルの最も重要な応用であるグレンジャー因果テストに進みましょう。
グレンジャー因果テスト
グレンジャー因果の定義
クライブ・グレンジャー（2003年ノーベル経済学賞）は、因果関係の操作的な定義として次の概念を提案しました。
「変数 $y_2$ の過去の情報が、$y_1$ の予測を有意に改善するとき、$y_2$ は $y_1$ のグレンジャー原因（Granger cause）である」
注意: グレンジャー因果は「予測における有用性」を測るものであり、真の因果関係を示すものではありません。
検定の枠組み
2変量VAR(p)モデルにおいて、「$y_2$ は $y_1$ のグレンジャー原因ではない」という帰無仮説は、
$$ H_0: \phi_{12}^{(1)} = \phi_{12}^{(2)} = \cdots = \phi_{12}^{(p)} = 0 $$
つまり、$y_1$ の方程式における $y_2$ のラグ変数の係数が全てゼロであることです。これは、通常のF検定で検定できます。
制約なしモデル（全てのラグを含む）と制約ありモデル（$y_2$ のラグを除外）の残差二乗和を比較し、F統計量を計算します。
$$ \begin{equation} F = \frac{(\text{RSS}_{\text{restricted}} – \text{RSS}_{\text{unrestricted}}) / p}{\text{RSS}_{\text{unrestricted}} / (T – 2p – 1)} \end{equation} $$
F統計量が有意であれば、$y_2$ のラグが $y_1$ の予測に統計的に有意な貢献をしている、すなわちグレンジャー因果が存在すると結論します。
グレンジャー因果は変数間の「予測的な方向性」を示しますが、ある変数へのショックが他の変数にどのように波及するかの動的な経路を示すものではありません。そこで登場するのがインパルス応答関数です。
インパルス応答関数（IRF）
MA($\infty$)表現
安定なVAR(p)過程は、単変量ARモデルと同様に、MA($\infty$)表現に書き換えることができます。
$$ \begin{equation} \bm{y}_t = \bm{\mu} + \sum_{i=0}^{\infty} \bm{\Psi}_i \bm{\varepsilon}_{t-i} $$
ここで $\bm{\Psi}_0 = \bm{I}_K$、$\bm{\Psi}_i$ は $K \times K$ の行列で、
$$ \bm{\Psi}_i = \sum_{j=1}^{i} \bm{\Psi}_{i-j} \bm{A}_j \quad (i = 1, 2, \dots) $$
により再帰的に計算できます（ただし $j > p$ のとき $\bm{A}_j = \bm{O}$）。
IRFの定義
$\bm{\Psi}_i$ の $(j, k)$ 成分 $\psi_{jk}^{(i)}$ は、
$$ \psi_{jk}^{(i)} = \frac{\partial y_{j,t+i}}{\partial \varepsilon_{k,t}} $$
を意味します。つまり、「時刻 $t$ に変数 $k$ のノイズに1単位のショックが加わったとき、$i$ 期後の変数 $j$ がどれだけ変化するか」を表します。
$\psi_{jk}^{(i)}$ を $i$ の関数としてプロットしたものがインパルス応答関数（Impulse Response Function, IRF）です。これにより、ショックの動的な波及過程を視覚化できます。
直交化インパルス応答
しかし、$\bm{\varepsilon}_t$ の成分間は一般に同時相関を持つため（$\bm{\Sigma}$ の非対角要素が非ゼロ）、「変数 $k$ のみにショックを与える」ということが厳密にはできません。
この問題を解決するために、$\bm{\Sigma}$ をコレスキー分解 $\bm{\Sigma} = \bm{P}\bm{P}^\top$（$\bm{P}$ は下三角行列）して、直交化されたノイズ $\bm{u}_t = \bm{P}^{-1} \bm{\varepsilon}_t$ を定義します。$\bm{u}_t$ の成分は互いに無相関で分散1です。
直交化インパルス応答は、
$$ \begin{equation} \bm{\Theta}_i = \bm{\Psi}_i \bm{P} \end{equation} $$
で定義されます。ただし、コレスキー分解の結果は変数の順序に依存するため、変数の順序付けには経済理論等に基づく正当化が必要です。
予測誤差の分散分解（FEVD）
定義
$h$ 期先の予測誤差は、
$$ \bm{y}_{t+h} – \hat{\bm{y}}_{t+h|t} = \sum_{i=0}^{h-1} \bm{\Psi}_i \bm{\varepsilon}_{t+h-i} = \sum_{i=0}^{h-1} \bm{\Theta}_i \bm{u}_{t+h-i} $$
変数 $j$ の $h$ 期先予測誤差の分散は、
$$ \text{MSE}_j(h) = \sum_{i=0}^{h-1} \sum_{k=1}^{K} (\theta_{jk}^{(i)})^2 $$
そのうち、変数 $k$ のショックに起因する部分の割合は、
$$ \begin{equation} \text{FEVD}_{jk}(h) = \frac{\sum_{i=0}^{h-1} (\theta_{jk}^{(i)})^2}{\sum_{i=0}^{h-1} \sum_{m=1}^{K} (\theta_{jm}^{(i)})^2} \end{equation} $$
分散分解により、「変数 $j$ の予測誤差のうち、何%が変数 $k$ のショックによって引き起こされたか」がわかります。自分自身のショックによる寄与が大きい変数は外生的（他の変数からの影響が少ない）、他の変数のショックによる寄与が大きい変数は内生的（他の変数からの影響を強く受ける）と解釈できます。
理論を一通り理解したところで、Pythonで実際にVARモデルを推定し、グレンジャー因果、インパルス応答、分散分解を実行してみましょう。
PythonによるVARモデルの実装
シミュレーションデータの生成
まず、真のパラメータが既知のVAR(2)データを生成し、モデルの推定とツールの使い方を確認します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# VAR(2)モデルのパラメータ
K = 2 # 変数の数
T = 500 # サンプルサイズ
# 係数行列
A1 = np.array([[0.5, 0.1],
[0.2, 0.3]])
A2 = np.array([[-0.2, 0.0],
[0.1, -0.1]])
# ノイズの分散共分散行列
Sigma = np.array([[1.0, 0.3],
[0.3, 0.8]])
# データ生成
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
for k in range(K):
axes[k].plot(Y[:, k], linewidth=0.8)
axes[k].set_title(f'Variable y_{k+1}')
axes[k].set_xlabel('Time')
plt.tight_layout()
plt.show()
生成された2つの時系列は、互いの過去の値から影響を受け合う相互依存構造を持っています。$y_1$ から $y_2$ への影響（$A_1$ の $(2,1)$ 要素 $= 0.2$）と $y_2$ から $y_1$ への影響（$A_1$ の $(1,2)$ 要素 $= 0.1$）が設定されているため、両方向のグレンジャー因果が検出されることが期待されます。
statsmodelsによるVAR推定
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
np.random.seed(42)
# データ生成（前のセクションと同じ）
K, T = 2, 500
A1 = np.array([[0.5, 0.1], [0.2, 0.3]])
A2 = np.array([[-0.2, 0.0], [0.1, -0.1]])
Sigma = np.array([[1.0, 0.3], [0.3, 0.8]])
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
# DataFrameに変換
df = pd.DataFrame(Y, columns=['y1', 'y2'])
# VARモデルの推定
model = VAR(df)
# ラグ次数の選択
lag_order = model.select_order(maxlags=10)
print("Lag Order Selection:")
print(lag_order.summary())
# 最適ラグで推定
result = model.fit(maxlags=2)
print("\n" + "=" * 50)
print("VAR(2) Estimation Results")
print("=" * 50)
print(result.summary())
ラグ次数選択の結果と推定されたパラメータを確認しましょう。
- ラグ次数選択: AICやBICの値が最小になるラグが $p = 2$ であれば、真のモデル次数を正しく同定しています。
- 推定されたパラメータ: $\bm{A}_1$ と $\bm{A}_2$ の推定値が真の値に近いことが確認できます。
グレンジャー因果テスト
import numpy as np
import pandas as pd
from statsmodels.tsa.api import VAR
np.random.seed(42)
# データ生成
K, T = 2, 500
A1 = np.array([[0.5, 0.1], [0.2, 0.3]])
A2 = np.array([[-0.2, 0.0], [0.1, -0.1]])
Sigma = np.array([[1.0, 0.3], [0.3, 0.8]])
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
df = pd.DataFrame(Y, columns=['y1', 'y2'])
model = VAR(df)
result = model.fit(maxlags=2)
# グレンジャー因果テスト
print("Granger Causality Tests:")
print("=" * 50)
# y2 → y1 のグレンジャー因果
print("\n[Test: y2 Granger-causes y1]")
gc_result1 = result.test_causality('y1', 'y2', kind='f')
print(gc_result1.summary())
# y1 → y2 のグレンジャー因果
print("\n[Test: y1 Granger-causes y2]")
gc_result2 = result.test_causality('y2', 'y1', kind='f')
print(gc_result2.summary())
グレンジャー因果テストの結果から、因果の方向性を読み取れます。
- $y_2 \to y_1$ のグレンジャー因果: $A_1$ の $(1,2)$ 要素が $0.1$ と比較的小さいため、検出されるかどうかはサンプルサイズに依存します。$T = 500$ であれば検出されることが多いですが、p値がやや大きい場合もあります。
- $y_1 \to y_2$ のグレンジャー因果: $A_1$ の $(2,1)$ 要素が $0.2$ とやや大きいため、こちらの方が有意になりやすいです。
インパルス応答関数の可視化
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
np.random.seed(42)
# データ生成
K, T = 2, 500
A1 = np.array([[0.5, 0.1], [0.2, 0.3]])
A2 = np.array([[-0.2, 0.0], [0.1, -0.1]])
Sigma = np.array([[1.0, 0.3], [0.3, 0.8]])
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
df = pd.DataFrame(Y, columns=['y1', 'y2'])
model = VAR(df)
result = model.fit(maxlags=2)
# インパルス応答関数
irf = result.irf(periods=20)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 直交化IRF（コレスキー分解ベース）
for i in range(2):
for j in range(2):
irf_vals = irf.orth_irfs[:, i, j]
lower = irf.orth_ci[:, i, j, 0] # 信頼区間下限
upper = irf.orth_ci[:, i, j, 1] # 信頼区間上限
axes[i, j].plot(irf_vals, linewidth=2, color='steelblue')
axes[i, j].fill_between(range(len(irf_vals)), lower, upper,
color='steelblue', alpha=0.15)
axes[i, j].axhline(y=0, color='gray', linestyle='--', alpha=0.5)
axes[i, j].set_title(f'Response of y{i+1} to y{j+1} shock')
axes[i, j].set_xlabel('Periods')
plt.suptitle('Orthogonalized Impulse Response Functions', fontsize=13)
plt.tight_layout()
plt.show()
4つのパネルから、ショックの波及過程を読み取れます。
- $y_1 \to y_1$（左上）: 自分自身へのショックは初期に最も大きく、その後指数的に減衰しています。約10期で概ねゼロに近づいています。
- $y_2 \to y_1$（右上）: $y_2$ のショックが $y_1$ に波及するクロス効果です。初期は小さいですが、1〜2期のラグを持って影響が現れています。
- $y_1 \to y_2$（左下）: $y_1$ のショックが $y_2$ に波及する効果です。$A_1(2,1) = 0.2$ が比較的大きいため、より明確な波及効果が見られます。
- $y_2 \to y_2$（右下）: $y_2$ 自身のショックへの応答です。
全てのIRFがゼロに収束していることは、VARモデルの安定性条件が満たされていることの反映です。
分散分解
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
np.random.seed(42)
# データ生成
K, T = 2, 500
A1 = np.array([[0.5, 0.1], [0.2, 0.3]])
A2 = np.array([[-0.2, 0.0], [0.1, -0.1]])
Sigma = np.array([[1.0, 0.3], [0.3, 0.8]])
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
df = pd.DataFrame(Y, columns=['y1', 'y2'])
model = VAR(df)
result = model.fit(maxlags=2)
# 分散分解
fevd = result.fevd(periods=20)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for idx in range(2):
decomp = fevd.decomp[:, idx, :] # (periods, K)
axes[idx].stackplot(range(20), decomp[:, 0], decomp[:, 1],
labels=['y1 shock', 'y2 shock'],
colors=['steelblue', 'coral'], alpha=0.7)
axes[idx].set_title(f'FEVD of y{idx+1}')
axes[idx].set_xlabel('Forecast Horizon')
axes[idx].set_ylabel('Proportion')
axes[idx].set_ylim(0, 1)
axes[idx].legend(loc='center right')
plt.suptitle('Forecast Error Variance Decomposition', fontsize=13)
plt.tight_layout()
plt.show()
# 数値での確認
print("FEVD at horizon 1, 5, 10, 20:")
for h in [0, 4, 9, 19]:
print(f"\n Horizon {h+1}:")
for k in range(K):
contrib = fevd.decomp[h, k, :]
print(f" y{k+1}: y1_shock={contrib[0]:.3f}, y2_shock={contrib[1]:.3f}")
分散分解の結果から、各変数の変動の源泉がわかります。
- $y_1$ の予測誤差: 短期的には $y_1$ 自身のショックがほぼ全てを説明しますが、予測ホライズンが長くなるにつれて $y_2$ のショックの寄与が増加します。これは $y_2 \to y_1$ の影響が時間をかけて蓄積されることを反映しています。
- $y_2$ の予測誤差: 初期から $y_1$ のショックの寄与がある程度あります。これは $A_1(2,1) = 0.2$ が比較的大きいため、$y_1$ のショックが $y_2$ に速やかに波及することを示しています。
予測
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
np.random.seed(42)
# データ生成
K, T = 2, 500
A1 = np.array([[0.5, 0.1], [0.2, 0.3]])
A2 = np.array([[-0.2, 0.0], [0.1, -0.1]])
Sigma = np.array([[1.0, 0.3], [0.3, 0.8]])
Y = np.zeros((T, K))
eps = np.random.multivariate_normal(np.zeros(K), Sigma, T)
for t in range(2, T):
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps[t]
df = pd.DataFrame(Y, columns=['y1', 'y2'])
# 訓練/テスト分割
train = df[:450]
test = df[450:]
# 推定と予測
model = VAR(train)
result = model.fit(maxlags=2)
forecast = result.forecast(train.values[-2:], steps=50)
forecast_df = pd.DataFrame(forecast, columns=['y1', 'y2'],
index=range(450, 500))
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
for k in range(K):
axes[k].plot(train.index[-100:], train.iloc[-100:, k],
linewidth=1, label='Train')
axes[k].plot(test.index, test.iloc[:, k],
linewidth=1.5, color='black', label='Test')
axes[k].plot(forecast_df.index, forecast_df.iloc[:, k],
linewidth=2, linestyle='--', color='red', label='Forecast')
axes[k].set_title(f'y{k+1} Forecast')
axes[k].legend()
axes[k].axvline(x=450, color='gray', linestyle=':', alpha=0.5)
plt.tight_layout()
plt.show()
予測結果を見ると、VARモデルの予測はホライズンが長くなるにつれて急速にゼロ（平均）に収束しています。これは定常VARモデルの一般的な特性で、長期予測は無条件平均に近づきます。短期（数ステップ先）の予測ではデータのパターンを捉えていますが、長期予測の精度は限定的です。
まとめ
本記事では、VARモデルの理論と実装について詳しく解説しました。
- VARモデルの構造: 複数の変数が互いの過去の値から影響を受け合う多変量自己回帰モデル。行列表現でコンパクトに定式化できる
- 安定性条件: コンパニオン行列の全ての固有値の絶対値が1未満であることが、定常性の必要十分条件
- 推定: 各方程式を個別にOLSで推定しても一致推定量が得られる。ラグ次数はAIC/BICで選択
- グレンジャー因果: ある変数の過去の情報が他の変数の予測に統計的に有意な貢献をするかをF検定で検定
- インパルス応答関数: ショックの動的な波及過程を可視化。コレスキー分解による直交化が標準
- 分散分解: 予測誤差のうち各変数のショックに起因する割合を定量化
VARモデルは多変量時系列分析の基礎であり、構造VARモデル（SVAR）、ベクトル誤差修正モデル（VECM）、動的因子モデルなど、多くの発展的手法の出発点になっています。
次のステップとして、以下の記事も参考にしてください。