書記が数学やるだけ#938 VARモデル
多変量の自己回帰モデルについて。
問題
説明
VAR（ベクトル自己回帰）モデルは、単変量自己回帰モデルを多変量へと一般化した動学的統計モデルである。本モデルを用いた分析では、主に以下の3点を通じて変数の動態を明らかにする。
Granger因果性検定：変数間における情報の先行性と予測に対する有効性を確認する。
インパルス応答関数：特定のショックに対する各変数の反応から、影響の方向性と持続性を把握する。
予測誤差分散分解（FEVD）：各変数の変動に対する他変数の寄与度を算出し、相互依存関係の強さを定量的に評価する。
解答
定常性の判定には，係数行列の固有値を用いる。
実際に可視化してみる。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import grangercausalitytests
# ---------------------------------------------------------
# 1. パラメータ設定 (演習問題の数値)
# ---------------------------------------------------------
phi_1 = np.array([[0.5, 0.3], # y1_t = 0.5*y1_{t-1} + 0.3*y2_{t-1}
[0.2, 0.4]]) # y2_t = 0.2*y1_{t-1} + 0.4*y2_{t-1}
sigma_u = np.array([[1.0, 0.6], # 誤差項の分散共分散行列 (同時点の相関あり)
[0.6, 1.0]])
# ---------------------------------------------------------
# 2. データのシミュレーション生成
# ---------------------------------------------------------
n_obs = 1000
np.random.seed(42)
# 相関を持つホワイトノイズの生成
errors = np.random.multivariate_normal([0, 0], sigma_u, size=n_obs)
y = np.zeros((n_obs, 2))
# VAR(1) プロセスの生成
for t in range(1, n_obs):
y[t] = phi_1 @ y[t-1] + errors[t]
df = pd.DataFrame(y, columns=['y1', 'y2'])
# ---------------------------------------------------------
# 3. 生データの可視化 (EDA)
# ---------------------------------------------------------
# 3-1. 時系列プロット (最初の100期間)
plt.figure(figsize=(12, 5))
plt.plot(df['y1'], label='y1 (Target 1)', alpha=0.8, color='dodgerblue')
plt.plot(df['y2'], label='y2 (Target 2)', alpha=0.8, color='darkorange')
plt.title("Simulated VAR(1) Time Series Data (t=0 to 100)", fontsize=14)
plt.xlabel("Time (t)")
plt.ylabel("Value")
plt.legend(loc='upper right')
plt.grid(True, linestyle='--', alpha=0.6)
plt.xlim(0, 100)
plt.show()
# 3-2. 散布図 (同時点の連動性を確認)
plt.figure(figsize=(6, 6))
plt.scatter(df['y1'], df['y2'], alpha=0.3, color='purple')
plt.title("Correlation Check: y1 vs y2", fontsize=14)
plt.xlabel("y1")
plt.ylabel("y2")
plt.grid(True)
plt.show()
「y2 が y1 をグレンジャー因果しない」とは、y1 の予測に y2 の過去値が寄与しないことを指し，y2,t−1 の係数を 0 とした形式で表される。
実際に検定してみると，y2 → y1・y1 → y2ともに有意であり，グレンジャー因果することが示せる。
# ---------------------------------------------------------
# 4. VARモデルの推定
# ---------------------------------------------------------
model = VAR(df)
results = model.fit(maxlags=1)
print("\n" + "="*50)
print("--- 推定された係数行列 ---")
print(results.params)
print("="*50)
# ---------------------------------------------------------
# 5. グレンジャー因果性検定
# ---------------------------------------------------------
print("\n【検定1】y2 が y1 をグレンジャー因果するか (y2 -> y1)")
grangercausalitytests(df[['y1', 'y2']], maxlag=1)
print("\n" + "-"*50)
print("【検定2】y1 が y2 をグレンジャー因果するか (y1 -> y2)")
grangercausalitytests(df[['y2', 'y1']], maxlag=1)==================================================
--- 推定された係数行列 ---
y1 y2
const -0.057518 -0.003857
L1.y1 0.472356 0.200295
L1.y2 0.323139 0.396099
==================================================
【検定1】y2 が y1 をグレンジャー因果するか (y2 -> y1)
Granger Causality
number of lags (no zero) 1
ssr based F test: F=73.8761 , p=0.0000 , df_denom=996, df_num=1
ssr based chi2 test: chi2=74.0986 , p=0.0000 , df=1
likelihood ratio test: chi2=71.4793 , p=0.0000 , df=1
parameter F test: F=73.8761 , p=0.0000 , df_denom=996, df_num=1
--------------------------------------------------
【検定2】y1 が y2 をグレンジャー因果するか (y1 -> y2)
Granger Causality
number of lags (no zero) 1
ssr based F test: F=37.4069 , p=0.0000 , df_denom=996, df_num=1
ssr based chi2 test: chi2=37.5196 , p=0.0000 , df=1
likelihood ratio test: chi2=36.8322 , p=0.0000 , df=1
parameter F test: F=37.4069 , p=0.0000 , df_denom=996, df_num=1インパルス応答関数の計算により「ある変数のショックが、他の変数に何期後に、どの程度の大きさで伝わるか」というダイナミクスを記述できる。
通常のインパルス応答では攪乱項 ϵ1,t と ϵ2,t が互いに独立であると仮定して計算したが、現実のデータではこれらは相関していることが多いため、コレスキー分解を用いてショックを直交化する。
予測誤差分散分解（FEVD: Forecast Error Variance Decomposition）では，「各変数の予測誤差のばらつきが、どの変数のショックによって説明されるか」の割合が示される。
可視化は以下の通り。
# ---------------------------------------------------------
# 6. 波及効果と寄与度の分析 (IRF & FEVD)
# ---------------------------------------------------------
# 6-1. 直交化インパルス応答
irf = results.irf(10)
irf.plot(orth=True)
plt.suptitle("Orthogonalized Impulse Response", fontsize=15)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
# 6-2. 予測誤差分散分解
fevd = results.fevd(10)
fevd.plot()
plt.suptitle("Forecast Error Variance Decomposition", fontsize=15)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
ここまでの内容を，変数間の相関がないモデルで示すと以下のようになる。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import grangercausalitytests
# ---------------------------------------------------------
# 1. パラメータ設定 (因果関係なし)
# ---------------------------------------------------------
# 係数をすべてゼロにする (ラグ効果なし)
phi_1_no_causality = np.zeros((2, 2))
# 分散共分散行列を単位行列にする (同時点の相関なし)
sigma_u_no_causality = np.eye(2)
# ---------------------------------------------------------
# 2. データのシミュレーション生成 (N=1000)
# ---------------------------------------------------------
n_obs = 1000
np.random.seed(123) # 別のシードを使用
# 独立したホワイトノイズの生成
errors = np.random.multivariate_normal([0, 0], sigma_u_no_causality, size=n_obs)
y = np.zeros((n_obs, 2))
# VAR(1) プロセスの生成 (ただのホワイトノイズの蓄積)
for t in range(1, n_obs):
# y[t] = 0 @ y[t-1] + errors[t] => y[t] = errors[t]
y[t] = phi_1_no_causality @ y[t-1] + errors[t]
df_no_gc = pd.DataFrame(y, columns=['y1', 'y2'])
# ---------------------------------------------------------
# 3. 分析と可視化
# ---------------------------------------------------------
# 3-1. 時系列プロット
plt.figure(figsize=(12, 5))
plt.plot(df_no_gc['y1'], label='y1 (White Noise)', alpha=0.8, color='blue')
plt.plot(df_no_gc['y2'], label='y2 (White Noise)', alpha=0.8, color='orange')
plt.title("Simulated Time Series: No Causality (White Noise)", fontsize=14)
plt.xlabel("Time (t)")
plt.ylabel("Value")
plt.legend(loc='upper right')
plt.grid(True, linestyle='--', alpha=0.6)
plt.xlim(0, 100) # 拡大
plt.show()
# 3-2. 散布図 (同時点の相関なし)
plt.figure(figsize=(6, 6))
plt.scatter(df_no_gc['y1'], df_no_gc['y2'], alpha=0.3, color='purple')
plt.title("Scatter Plot: y1 vs y2 (No Correlation)", fontsize=14)
plt.xlabel("y1")
plt.ylabel("y2")
plt.grid(True)
plt.show()
# 3-3. インパルス応答 (IRF)
model_no_gc = VAR(df_no_gc)
results_no_gc = model_no_gc.fit(maxlags=1)
irf_no_gc = results_no_gc.irf(10)
irf_no_gc.plot(orth=True) # 直交化しても意味はない
plt.suptitle("Impulse Response: No Causality", fontsize=15)
plt.show()
# 3-4. グレンジャー因果性検定 (y2 -> y1)
print("\n" + "="*50)
print("【検定】y2 が y1 をグレンジャー因果するか (y2 -> y1) - 因果なし例")
grangercausalitytests(df_no_gc[['y1', 'y2']], maxlag=1)
==================================================
【検定】y2 が y1 をグレンジャー因果するか (y2 -> y1) - 因果なし例
Granger Causality
number of lags (no zero) 1
ssr based F test: F=0.2400 , p=0.6243 , df_denom=996, df_num=1
ssr based chi2 test: chi2=0.2407 , p=0.6237 , df=1
likelihood ratio test: chi2=0.2407 , p=0.6237 , df=1
parameter F test: F=0.2400 , p=0.6243 , df_denom=996, df_num=1本記事のもくじはこちら：
