← 時系列分析 一覧

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：AR・MAモデル（単変量AR）・定常性と自己相関　|　数理：確率過程（マルコフ連鎖・ポアソン過程）（統計）

要点（BLUF）

• VAR($p$) ＝ 単変量 AR の多変量版。各変数を「全変数の過去ラグ」で回帰します。係数はスカラーではなく行列 $A_1,\dots,A_p$ になり、変数どうしの相互作用（フィードバック）を表せます。
• 安定性条件：係数行列を積み上げたコンパニオン行列の固有値がすべて単位円内（絶対値 < 1）。VAR(1) なら $A_1$ の固有値が単位円内。満たさないと非定常（ランダムウォークと単位根）。
• 推定は各式を OLS（系列ごとに同じ右辺）。予測は点予測だけでなく予測区間を添えます。安定 VAR の区間は有限値（無条件分散）へ収束し、非定常のように青天井には広がりません。

1. VAR($p$)：各変数を全変数の過去で説明する

$K$ 本の系列 $\mathbf{y}_t=(y_{1t},\dots,y_{Kt})^\top$ を、互いの過去ラグで同時に説明します。

• $\mathbf{c}$ は $K$ 次元の定数ベクトル、$A_i$ は $K\times K$ の係数行列、$\boldsymbol\varepsilon_t$ は平均 0・共分散 $\Sigma$ の多変量ホワイトノイズ。
• 行 $j$ を取り出すと「$y_{jt}$ ＝ 定数 ＋（全変数の過去ラグの線形結合）＋ 誤差」。$A_i$ の $(j,k)$ 成分は「$k$ 番目の変数のラグ $i$ が $j$ 番目の変数に与える係数」です。対角成分が自分自身の慣性（単変量 AR 相当）、非対角成分が他変数からの波及。

2 変数 VAR(1) を成分で書くと相互依存が見えます：

$a_{12}\ne0$ なら「$y_2$ の過去が $y_1$ に効く」、$a_{21}\ne0$ なら逆向き。両方非ゼロなら相互フィードバックです。単変量 AR（AR・MAモデル）は $K=1$ の特殊ケース。

安定性条件（コンパニオン行列の固有値）

VAR($p$) は、ラグを状態に畳み込んだ VAR(1) 形（コンパニオン形） に書けます。$Kp\times Kp$ のコンパニオン行列

の固有値がすべて単位円内（絶対値 < 1）なら定常（安定）。VAR(1) のときはこれが $A_1$ そのものなので、$A_1$ の固有値を見ればよい。単変量 AR の「特性根が単位円外」条件（AR・MAモデル）の行列版です。

コード①：真の係数行列を仕込んで VAR で復元

真の $A_1$（固有値 0.6, 0.3 で安定）と定数 $\mathbf{c}$ を仕込んだ 2 変数 VAR(1) を生成し、statsmodels の VAR で係数行列を復元します。次数 $p$ は select_order（AIC）で選ばせます。

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

# 真の 2変数 VAR(1):  y_t = c + A1 y_{t-1} + ε_t
# A1 の固有値が単位円内 → 安定（固有値 0.6, 0.3）
c_true = np.array([1.0, 0.5])
A1_true = np.array([[0.5, 0.2],
                    [0.1, 0.4]])
Sigma = np.array([[1.0, 0.3],
                  [0.3, 1.0]])   # 残差の同時相関

np.random.seed(3)
n = 2000
eps = np.random.multivariate_normal([0, 0], Sigma, size=n)
Y = np.zeros((n, 2))
for t in range(1, n):
    Y[t] = c_true + A1_true @ Y[t-1] + eps[t]

# 安定性：A1 の固有値（絶対値）がすべて 1 未満か
eig = np.abs(np.linalg.eigvals(A1_true))
print("A1 固有値の絶対値 =", np.round(eig, 3), "→", "安定" if np.all(eig < 1) else "不安定")

# VAR を当てて次数選択（AIC）
model = VAR(Y)
sel = model.select_order(maxlags=5)
print("AIC が選ぶ次数 p =", sel.aic)

res = model.fit(1)            # VAR(1) を推定
A1_hat = res.coefs[0]         # 推定された係数行列（ラグ1）
print("\n真の A1:\n", np.round(A1_true, 3))
print("推定 A1:\n", np.round(A1_hat, 3))
print("\n真の c =", np.round(c_true, 3), " 推定 c =", np.round(res.params[0], 3))

出力：

A1 固有値の絶対値 = [0.6 0.3] → 安定
AIC が選ぶ次数 p = 1

真の A1:
 [[0.5 0.2]
 [0.1 0.4]]
推定 A1:
 [[0.494 0.247]
 [0.081 0.451]]

真の c = [1.  0.5]  推定 c = [0.999 0.514]

出力の意味：$A_1$ の固有値の絶対値は $0.6,\ 0.3$ でともに 1 未満——安定（定常）。AIC は正しく $p=1$ を選びました。推定行列 $\hat A_1$ は真値 $A_1$ を成分ごとに復元（$0.494,\,0.247,\,0.081,\,0.451$ vs 真 $0.5,\,0.2,\,0.1,\,0.4$）、定数も $0.999,\,0.514$ と当たっています。**非対角成分 $a_{12}\approx0.25,\ a_{21}\approx0.08$ が「変数間の波及」**で、これがあるからこそ単変量 AR では捉えられない相互作用を VAR が表せます。

2. 多変量予測と予測区間

VAR の $h$ 期先予測は、推定した係数行列で漸化的に回します（$\hat{\mathbf{y}}_{t+1}=\hat{\mathbf{c}}+\hat A_1\mathbf{y}_t+\cdots$、次は予測値を代入）。予測誤差の分散も漸化式で積み上がり、各変数に予測区間が付きます。安定 VAR では区間幅が有限値（無条件分散）へ収束する点が、非定常（ARIMA の和分）で青天井に広がるのと対照的です（ARMA・ARIMAモデル）。

flowchart LR
    A["多変量系列 y_t（K本）"] --> B["VAR(p) を各式 OLS で推定"]
    B --> C["次数 p を AIC/BIC で選択"]
    C --> D["forecast_interval で h 期先を予測"]
    D --> E["点予測 + 95%予測区間（各変数）"]
    B --> F["IRF: ショックの波及を追跡（任意）"]

コード②：予測区間つき多変量予測をホールドアウトで評価

末尾 20 期を検証に回し（時間順分割・シャッフルなし）、forecast_interval で 点予測＋95% 予測区間を出して評価・図示します。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from statsmodels.tsa.api import VAR

# 同じ真の VAR(1) を生成（安定）
c_true = np.array([1.0, 0.5])
A1_true = np.array([[0.5, 0.2],
                    [0.1, 0.4]])
Sigma = np.array([[1.0, 0.3], [0.3, 1.0]])
np.random.seed(3)
n = 600
eps = np.random.multivariate_normal([0, 0], Sigma, size=n)
Y = np.zeros((n, 2))
for t in range(1, n):
    Y[t] = c_true + A1_true @ Y[t-1] + eps[t]

# ホールドアウト：末尾 20 期を検証に回し、過去だけで推定
h = 20
train, test = Y[:-h], Y[-h:]
res = VAR(train).fit(1)

# forecast_interval：点予測 + 95% 予測区間（直近 p 期を渡す）
point, lower, upper = res.forecast_interval(train[-res.k_ar:], steps=h, alpha=0.05)

# 検証：各変数の RMSE と区間カバレッジ
for j, name in enumerate(["y1", "y2"]):
    rmse = np.sqrt(np.mean((point[:, j] - test[:, j])**2))
    cover = np.mean((lower[:, j] <= test[:, j]) & (test[:, j] <= upper[:, j]))
    w1 = upper[0, j] - lower[0, j]
    wH = upper[-1, j] - lower[-1, j]
    print(f"{name}: RMSE={rmse:.3f}  95%被覆={cover*100:.0f}%  区間幅 1期={w1:.2f}→{h}期={wH:.2f}")

# 図示（y1 のみ：実測・点予測・予測区間）
t_axis = np.arange(n)
plt.figure(figsize=(9, 4.5))
plt.plot(t_axis[-60:-h], train[-60+h:, 0], color="gray", lw=1, label="訓練(実測)")
plt.plot(t_axis[-h:], test[:, 0], color="k", lw=1.5, marker="o", ms=3, label="検証(実測)")
plt.plot(t_axis[-h:], point[:, 0], color="C1", lw=2, label="点予測")
plt.fill_between(t_axis[-h:], lower[:, 0], upper[:, 0], color="C1", alpha=0.25, label="95%予測区間")
plt.axvline(n-h-1, ls=":", color="k")
plt.xlabel("時点 t"); plt.ylabel("y1"); plt.legend()
plt.title("VAR(1) による多変量予測（y1・予測区間つき）")
plt.tight_layout(); plt.show()

出力：

y1: RMSE=1.047  95%被覆=100%  区間幅 1期=3.97→20期=5.09
y2: RMSE=1.065  95%被覆=100%  区間幅 1期=3.86→20期=4.49

出力の意味：各変数の RMSE は約 1.0（残差の標準偏差 1.0 と整合）で、検証 20 点はすべて 95% 区間に収まりました。注目は区間幅が（y1 で）$3.97\to5.09$ と少し広がってから頭打ちになること——安定 VAR の予測分散は無条件分散へ収束するので、ARIMA の和分系列のように $\sqrt{h}$ で青天井に広がりません（ARMA・ARIMAモデル）。図では橙の点予測がやがて無条件平均へ落ち着き、帯の幅も一定に近づきます。検証 20 点全被覆は標本が 1 本のため目安で、較正は本来ウォークフォワードで複数回確かめます（予測の評価指標と時系列CV）。

3. インパルス応答（IRF）：ショックの波及を追う（任意）

VAR の魅力のひとつがインパルス応答関数（IRF）。「ある変数に大きさ 1 のショックを与えたら、各変数が何期にわたってどう反応するか」を係数行列の累乗で追います。安定系ならショックの影響は時間とともに 0 へ減衰します。

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

c_true = np.array([1.0, 0.5])
A1_true = np.array([[0.5, 0.2],
                    [0.1, 0.4]])
Sigma = np.array([[1.0, 0.3], [0.3, 1.0]])
np.random.seed(3)
n = 2000
eps = np.random.multivariate_normal([0, 0], Sigma, size=n)
Y = np.zeros((n, 2))
for t in range(1, n):
    Y[t] = c_true + A1_true @ Y[t-1] + eps[t]

res = VAR(Y).fit(1)
irf = res.irf(8)                 # 8 期先までのインパルス応答
resp = irf.irfs               # shape (steps+1, 2, 2): [期, 応答変数, ショック変数]

print("y1 に大きさ1のショック → 各変数の反応（直交化なし）")
print(f"{'h':>2}{'→y1':>8}{'→y2':>8}")
for hh in [0, 1, 2, 4, 8]:
    print(f"{hh:>2}{resp[hh,0,0]:>8.3f}{resp[hh,1,0]:>8.3f}")
print("反応は時間とともに 0 へ減衰 →", "安定系" if abs(resp[8,0,0]) < abs(resp[0,0,0]) else "発散")

出力：

y1 に大きさ1のショック → 各変数の反応（直交化なし）
 h     →y1     →y2
 0   1.000   0.000
 1   0.494   0.081
 2   0.264   0.076
 4   0.087   0.037
 8   0.012   0.006

出力の意味：$y_1$ への 1 単位ショックは、即時（$h=0$）に $y_1$ を 1、$y_2$ を 0 動かし、翌期（$h=1$）には $y_1\to0.494,\ y_2\to0.081$——これはちょうど $A_1$ の 1 列目です（係数行列がそのまま 1 期波及を決める）。以降は累乗で減衰し、$h=8$ では $0.012,\,0.006$ とほぼ消えます。安定系ではショックが永続しないことが固有値 < 1 から保証されます。なお実務の IRF は誤差の同時相関 $\Sigma$ を考慮した直交化 IRF（コレスキー分解）を使い、変数の順序付けに結果が依存する——構造的な解釈には注意が要ります（後述）。

4. 数式の直観

• 係数行列 ＝ 波及の地図：$A_i$ の非対角成分が「他変数の過去が自分に効く強さ」。VAR は『どの変数がどの変数を、どのラグで動かすか』を行列にまとめた装置です。
• 固有値 < 1 ＝ ショックが減衰する：コンパニオン行列の固有値が単位円内なら、過去のショックの寄与が幾何級数的に消え、系列は無条件平均へ回帰します。これが定常性の正体で、予測区間が有限に収束する理由でもあります。
• VAR は『関係の方向』までは決めない：相関と同時性を捉えますが、どちらが原因かは VAR 自体は語りません。次の グレンジャー因果 は「過去が予測を改善するか」で方向の手がかりを与えますが、それも予測的であって構造的因果ではありません。

⚠️ よくある誤解・落とし穴

• 「変数を増やすほど良い」ではない：VAR のパラメータ数は $K^2p$ で、変数や次数とともに爆発的に増えます。データに対して係数が多すぎると過学習し予測が悪化（次元の呪い）。AIC/BIC で倹約的に選び、多変量なら正則化 VAR やベイズ VAR を検討します。
• 「非定常系列にそのまま VAR」ではない：各系列が単位根（ランダムウォークと単位根）を持つなら、レベルの VAR は見せかけの関係を生みがち。共和分があるなら VECM（共和分と誤差修正モデル（VECM））、無いなら差分してから VAR が定石です。
• 「IRF や係数＝因果」ではない：直交化 IRF はコレスキー順序（変数の並べ方）に依存し、順序を変えると解釈が変わります。VAR の係数は予測的関係で、構造的因果は別の枠組み（グレンジャー因果 の⚠️、因果サイト 構造的因果モデルとdo演算子）。
• 「点予測だけ」ではない：多変量でも必ず予測区間を添えます。安定 VAR は有限へ収束、非定常は広がる——どちらかを区間の形で示すこと。

関連ノート

• AR・MAモデル（単変量 AR——VAR の $K=1$ 版）
• ARMA・ARIMAモデル（予測区間・和分での区間の広がり）
• グレンジャー因果（VAR の係数から予測的因果を検定）
• 共和分と誤差修正モデル（VECM）（非定常な多変量・長期均衡）
• 定常性と自己相関・ランダムウォークと単位根（定常性・単位根）
• 状態空間モデルの枠組み（多変量の状態空間表現）
• 確率過程（マルコフ連鎖・ポアソン過程）（統計・多変量定常過程の土台）
• 第5章 多変量時系列 目次
• 時系列分析・予測テキスト 全体目次