原油価格が上がると、数週間後にガソリン価格が上がる — これは直感的に当然のことに思えます。では、逆にガソリン価格の変動が原油価格に影響を与えることはあるでしょうか？ あるいは、株式市場のある銘柄が急落したとき、その衝撃波が他の銘柄にどのように伝播し、何日後にピークを迎え、いつ収束するのでしょうか？
こうした「複数の時系列がどう関係し合っているのか」という問いは、単純な相関係数だけでは答えられません。なぜなら、時系列データには 時間のずれ（ラグ） と 方向性（AがBに影響するのか、BがAに影響するのか） という2つの次元があるからです。静的な相関係数はこの両方を無視してしまいます。
多変量時系列の関係性分析は、以下のような幅広い分野で不可欠です。
- マクロ経済分析: GDP・金利・失業率・インフレ率の相互作用を分析し、金融政策の効果を定量化します。中央銀行が金利を引き上げたとき、その影響が失業率に現れるまでのタイムラグと波及パターンを予測できます
- 金融工学: 複数の資産間のリード・ラグ関係を発見し、ポートフォリオのリスク管理やアービトラージ戦略に活用します
- 脳科学・神経科学: 脳の異なる領域間の信号伝達パターンを解析し、認知機能のメカニズムを解明します
- 気象・環境科学: 海面温度とエルニーニョ現象、CO2濃度と気温変動など、地球規模の因果連鎖を追跡します
本記事の内容
- 相互相関関数（CCF）による時間ずれの検出
- グレンジャー因果検定による予測的因果関係の判定
- VARモデルを基盤としたインパルス応答関数（IRF）の理論と導出
- 分散分解（FEVD）による変動要因の定量的分離
- 共和分による長期均衡関係の検出
- 5つの手法を統一的に使い分けるための実践ガイド
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
多変量時系列分析の全体像 — 5つのツールの位置づけ
多変量時系列の関係性を分析するツールは複数ありますが、それぞれが異なる「問い」に答えます。まず全体の地図を頭に入れてから、個々の手法に進みましょう。
ある日、あなたが2つの時系列 $\{x_t\}$ と $\{y_t\}$ を持っているとします。これらの関係について、自然に浮かぶ疑問は次の5つです。
問い1: 「$x$ と $y$ は時間的にずれて連動しているか？」 — これに答えるのが 相互相関関数（CCF） です。CCFは時間のずれ（ラグ）を横軸に取り、各ラグでの相関を計算します。「$x$ が $y$ より3期先行して動いている」といったリード・ラグ構造を発見できます。
問い2: 「$x$ の過去の情報は $y$ の予測に役立つか？」 — これに答えるのが グレンジャー因果検定 です。統計的検定の枠組みで「予測的因果」の有無を判定します。CCFが単なる相関の大きさを見るのに対し、グレンジャー因果は「$y$ 自身の過去だけでは説明できない情報が $x$ の過去に含まれているか」を問います。
問い3: 「$x$ に突発的なショックが加わったとき、$y$ はどう応答するか？」 — これに答えるのが インパルス応答関数（IRF） です。VARモデルを推定した上で、ある変数への1単位のショックが他の変数にどのような時間的パターンで波及するかを追跡します。CCFやグレンジャー因果では見えない「動的な波及メカニズム」を可視化できます。
問い4: 「$y$ の予測誤差のうち、どれだけが $x$ のショックに起因するか？」 — これに答えるのが 分散分解（FEVD） です。IRFが「波及のパターン」を示すのに対し、FEVDは「波及の大きさ」を定量化します。各変数の予測誤差分散を、システム内の各ショックに帰属させます。
問い5: 「$x$ と $y$ は短期的にはずれるが、長期的に均衡関係を保っているか？」 — これに答えるのが 共和分（Cointegration） です。個々の系列が非定常（単位根を持つ）であっても、それらの線形結合が定常であれば「長期均衡」が存在します。
これら5つのツールは排他的ではなく、相補的です。実務では複数の手法を組み合わせて、多角的な視点から関係性を検証します。
それでは、最もシンプルな道具から始めましょう。まずは2つの時系列の「時間ずれ付きの相関」を測る相互相関関数です。
相互相関関数（CCF）の理論
「相関」に時間軸を入れる
1変量の自己相関関数（ACF）は、同じ時系列の異なる時点間の相関を測りました。これを2つの異なる時系列に拡張したのが 相互相関関数（Cross-Correlation Function, CCF） です。
日常的なアナロジーで考えてみましょう。あなたが川の上流と下流に温度計を置いたとします。上流で温泉水が流れ込んで水温が上がると、しばらく時間が経ってから下流でも水温が上がります。このとき、上流の温度変化と下流の温度変化は「完全に同時」には連動しませんが、「適切な時間ずれ」を入れれば強い相関を示すはずです。CCFは、まさにこの「最適な時間ずれ」を体系的に見つけるためのツールです。
数学的定義
2つの定常時系列 $\{x_t\}$ と $\{y_t\}$ に対して、相互共分散関数（Cross-Covariance Function） を次のように定義します。
$$ \gamma_{xy}(h) = \text{Cov}(x_{t+h}, y_t) = E[(x_{t+h} – \mu_x)(y_t – \mu_y)] $$
ここで $h$ はラグ（時間のずれ）、$\mu_x = E[x_t]$、$\mu_y = E[y_t]$ はそれぞれの平均です。定常性の仮定により、この共分散は $t$ に依存せず $h$ のみの関数になります。
これを標準化して 相互相関関数（CCF） を得ます。
$$ \rho_{xy}(h) = \frac{\gamma_{xy}(h)}{\sqrt{\gamma_{xx}(0) \cdot \gamma_{yy}(0)}} = \frac{\gamma_{xy}(h)}{\sigma_x \sigma_y} $$
$\gamma_{xx}(0) = \text{Var}(x_t) = \sigma_x^2$ および $\gamma_{yy}(0) = \text{Var}(y_t) = \sigma_y^2$ を分母に置くことで、$-1 \leq \rho_{xy}(h) \leq 1$ の範囲に正規化されます。
自己相関関数との重要な違い
自己相関関数では $\gamma_{xx}(h) = \gamma_{xx}(-h)$ — つまりラグの正負で対称でした。しかし相互相関関数では一般に次が成り立ちます。
$$ \gamma_{xy}(h) = \gamma_{yx}(-h) $$
これは「$x$ がラグ $h$ で $y$ に先行する」ことと「$y$ がラグ $-h$ で $x$ に先行する」ことが同じ情報を表していることを意味します。したがって、CCFのグラフでは ラグの正の側と負の側で異なる情報 が得られます。$h > 0$ のピークは $x$ が $y$ に先行していることを、$h < 0$ のピークは $y$ が $x$ に先行していることを示唆します。
標本相互相関関数
実際のデータ $x_1, x_2, \ldots, x_T$ と $y_1, y_2, \ldots, y_T$ からCCFを推定するには、標本版を使います。
$$ \hat{\gamma}_{xy}(h) = \frac{1}{T} \sum_{t=1}^{T-h} (x_{t+h} – \bar{x})(y_t – \bar{y}), \quad h \geq 0 $$
ここで除数を $T$（$T-h$ ではなく）とするのは、推定量のバイアスを減らす慣例です。そして標本CCFは次のように標準化します。
$$ \hat{\rho}_{xy}(h) = \frac{\hat{\gamma}_{xy}(h)}{\sqrt{\hat{\gamma}_{xx}(0) \cdot \hat{\gamma}_{yy}(0)}} $$
2つの系列が無相関であるとき、$\hat{\rho}_{xy}(h)$ は近似的に $N(0, 1/T)$ に従います。したがって、$\pm 2/\sqrt{T}$ を信頼区間として描くことで、統計的に有意なラグを視覚的に判断できます。
理論を理解したところで、実際にPythonで相互相関を計算し、可視化してみましょう。
相互相関関数のPython実装
まず、2つの合成時系列データを生成します。$x_t$ が $y_t$ に3期先行する構造を人工的に作り、CCFがこのリード・ラグ関係を正しく検出できるか確認します。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import ccf
np.random.seed(42)
# パラメータ
T = 500 # 時系列の長さ
lag_true = 3 # xがyに先行するラグ
# x_t: AR(1)プロセス
x = np.zeros(T)
for t in range(1, T):
x[t] = 0.7 * x[t-1] + np.random.randn()
# y_t: x_{t-3} に依存 + ノイズ
y = np.zeros(T)
for t in range(lag_true, T):
y[t] = 0.5 * x[t - lag_true] + 0.3 * y[t-1] + np.random.randn()
# 標本CCFを計算
ccf_values = ccf(x, y, nlags=20, alpha=None)
# 可視化
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# 時系列プロット
axes[0].plot(x[:200], label='$x_t$ (leading)', alpha=0.8)
axes[0].plot(y[:200], label='$y_t$ (lagging)', alpha=0.8)
axes[0].set_xlabel('Time')
axes[0].set_ylabel('Value')
axes[0].set_title('Two Time Series (first 200 points)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# CCFプロット
lags = np.arange(len(ccf_values))
conf_bound = 2 / np.sqrt(T)
axes[1].bar(lags, ccf_values, width=0.4, color='steelblue', alpha=0.8)
axes[1].axhline(y=conf_bound, color='red', linestyle='--', label=f'95% confidence (±{conf_bound:.3f})')
axes[1].axhline(y=-conf_bound, color='red', linestyle='--')
axes[1].axhline(y=0, color='black', linewidth=0.5)
axes[1].set_xlabel('Lag $h$')
axes[1].set_ylabel('$\\hat{\\rho}_{xy}(h)$')
axes[1].set_title('Cross-Correlation Function (CCF): $x_t$ leading $y_t$')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上のグラフから、2つの重要な特徴が読み取れます。
- CCFはラグ3で最大値を取ります — これはデータ生成時に $y_t = 0.5 x_{t-3} + \cdots$ と設定したリード・ラグ構造と完全に一致しています。CCFがラグ3でピークを示すということは、$x$ が $y$ に3期先行して影響を与えていることを意味します。
- ラグ3の前後（ラグ2, 4, 5あたり）にも有意な相関が残っています — これは $x_t$ 自体がAR(1)で自己相関を持ち、$y_t$ もAR成分を持つため、真のラグ3の影響が周辺のラグにも「にじむ」ためです。この現象は実データでも頻繁に起きるため、CCFのピーク位置だけでなく全体の形状を見ることが重要です。
- 赤い破線（95%信頼区間）を超えるラグが複数存在します — 信頼区間は「2系列が完全に無相関」という帰無仮説の下で構成されています。これを超えるラグでは、偶然ではなく実質的な相関があると判断できます。
CCFは「どのラグで相関が強いか」を教えてくれますが、それが本当に因果関係を意味するかは別問題です。例えば、$x$ と $y$ の両方に影響する第3の変数 $z$ が存在する場合、$x$ と $y$ の間にCCFのピークが現れても、$x \to y$ の直接的な因果関係は存在しないかもしれません。この「相関と因果のギャップ」を部分的に埋めるのが、次に紹介するグレンジャー因果検定です。
グレンジャー因果検定 — 予測的因果の判定
「因果」をどう定義するか
「$x$ が $y$ の原因である」と言いたいとき、日常的な因果の概念は哲学的に複雑です。クライヴ・グレンジャーは1969年、この難題をプラグマティックに回避する巧みな定義を提案しました。
身近な例で考えましょう。あなたが天気予報士で、明日の気温を予測しようとしています。これまでの気温データだけで予測するよりも、気圧のデータも加えたほうが予測精度が上がるなら、「気圧は気温のグレンジャー原因である」と言います。ポイントは「真の因果メカニズム」を問わないことです。あくまで 予測力の向上 という観測可能な基準に絞っています。
数学的定義
$y_t$ のAR($p$)モデル（制約モデル）を考えます。
$$ y_t = c + \sum_{i=1}^{p} \alpha_i y_{t-i} + \varepsilon_t $$
次に、$x$ の過去も追加した非制約モデルを考えます。
$$ y_t = c + \sum_{i=1}^{p} \alpha_i y_{t-i} + \sum_{i=1}^{p} \beta_i x_{t-i} + u_t $$
帰無仮説は $H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0$（$x$ の過去は $y$ の予測に貢献しない）です。この検定は標準的なF検定で実行できます。
$$ F = \frac{(\text{RSS}_R – \text{RSS}_U) / p}{\text{RSS}_U / (T – 2p – 1)} $$
ここで $\text{RSS}_R$ は制約モデルの残差平方和、$\text{RSS}_U$ は非制約モデルの残差平方和です。F統計量が有意に大きければ、$x$ の過去の情報が $y$ の予測に統計的に有意な貢献をしている — すなわち「$x$ は $y$ のグレンジャー原因である」と結論します。
重要な注意点
グレンジャー因果検定にはいくつかの落とし穴があります。
- 予測的因果であり、真の因果ではない: $x \to y$ のグレンジャー因果が検出されても、真に $x$ が $y$ の原因であるとは限りません。第3の共通原因 $z$ が存在し、$z$ が $x$ に先に影響し、遅れて $y$ に影響する場合でも、$x \to y$ のグレンジャー因果は検出されます。
- 定常性の仮定: 検定が正しく機能するには、両系列が定常である必要があります。非定常な系列に適用すると、見かけ上の因果関係（疑似因果）が検出されやすくなります。
- ラグ次数 $p$ の選択: 結果はラグ次数に敏感です。AICやBICなどの情報量基準で適切なラグを選ぶ必要があります。
グレンジャー因果検定の詳細な理論・導出・Pythonでの実装は、専用の記事で詳しく解説しています。
ここでは、先ほどの合成データに対してグレンジャー因果検定を実行し、結果を確認しましょう。
from statsmodels.tsa.stattools import grangercausalitytests
# x → y のグレンジャー因果検定（ラグ1〜6）
print("=== Granger Causality: x → y ===")
data_xy = np.column_stack([y, x]) # statsmodelsは [被説明変数, 説明変数] の順
result_xy = grangercausalitytests(data_xy, maxlag=6, verbose=True)
print("\n=== Granger Causality: y → x ===")
data_yx = np.column_stack([x, y])
result_yx = grangercausalitytests(data_yx, maxlag=6, verbose=True)
このコードの出力を見ると、$x \to y$ 方向ではラグ3以上で非常に小さなp値（0.05未満）が得られ、帰無仮説が棄却されます。一方、$y \to x$ 方向ではp値が大きく、帰無仮説は棄却されません。これは、データ生成の構造（$y_t$ は $x_{t-3}$ に依存するが、$x_t$ は $y$ の過去に依存しない）と完全に一致します。つまり、グレンジャー因果検定は因果の方向性を正しく識別しています。
CCFとグレンジャー因果は「関係の有無」と「因果の方向」を教えてくれましたが、ショックが波及する 動的なパターン — 最初はどのくらいの大きさで影響し、何期後にピークを迎え、いつ収束するのか — は教えてくれません。この動的な波及パターンを描くのがインパルス応答関数です。そしてインパルス応答を理解するには、まずその土台となるVARモデルを復習する必要があります。
VARモデルの復習
単変量ARから多変量VARへ
単変量のAR(1)モデルは $y_t = c + \phi y_{t-1} + \varepsilon_t$ でした。これを $k$ 個の変数からなるベクトル $\bm{y}_t = (y_{1t}, y_{2t}, \ldots, y_{kt})^\top$ に拡張したものが VARモデル です。
家計の経済を例にしましょう。今月の収入と支出は、先月の収入と支出の両方に影響を受けます。収入が増えれば支出も増え、支出が増えれば翌月の貯蓄（つまり可処分な収入の余裕）が減る、というように 変数同士が互いにフィードバック し合います。VARモデルはこの相互依存構造を自然に表現します。
VAR($p$)の定義
$k$ 変量のVAR($p$)モデルは次のように定義されます。
$$ \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{u}_t $$
ここで各記号は以下の意味を持ちます。
- $\bm{y}_t$: $k \times 1$ の内生変数ベクトル（時刻 $t$ での観測値）
- $\bm{c}$: $k \times 1$ の定数項ベクトル
- $\bm{A}_i$: $k \times k$ の係数行列（ラグ $i$ での変数間の影響を記述）
- $\bm{u}_t$: $k \times 1$ のホワイトノイズベクトル。$E[\bm{u}_t] = \bm{0}$、$E[\bm{u}_t \bm{u}_t^\top] = \bm{\Sigma}_u$
2変量VAR(1)の場合を書き下すと、構造が明確になります。
$$ \begin{pmatrix} y_{1t} \\ y_{2t} \end{pmatrix} = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} + \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} y_{1,t-1} \\ y_{2,t-1} \end{pmatrix} + \begin{pmatrix} u_{1t} \\ u_{2t} \end{pmatrix} $$
行列 $\bm{A}_1$ の非対角要素 $a_{12}$ と $a_{21}$ が、変数間のクロス効果を表します。$a_{12} \neq 0$ なら $y_{2,t-1}$（変数2の過去）が $y_{1t}$（変数1の現在）に影響し、$a_{21} \neq 0$ なら逆方向の影響が存在します。
VMA表現 — VARの無限移動平均変換
インパルス応答関数を理論的に導出するには、VARモデルを VMA（Vector Moving Average）表現 に変換する必要があります。これは、「$\bm{y}_t$ を過去と現在のショック（撹乱項）の線形結合として表す」ことを意味します。
安定なVAR(1)モデル $\bm{y}_t = \bm{A}_1 \bm{y}_{t-1} + \bm{u}_t$（定数項は簡単のため省略）を繰り返し代入してみましょう。
まず $\bm{y}_{t-1} = \bm{A}_1 \bm{y}_{t-2} + \bm{u}_{t-1}$ を代入すると、
$$ \bm{y}_t = \bm{A}_1(\bm{A}_1 \bm{y}_{t-2} + \bm{u}_{t-1}) + \bm{u}_t = \bm{A}_1^2 \bm{y}_{t-2} + \bm{A}_1 \bm{u}_{t-1} + \bm{u}_t $$
さらに $\bm{y}_{t-2}$ を代入して同じ操作を繰り返すと、
$$ \bm{y}_t = \bm{A}_1^3 \bm{y}_{t-3} + \bm{A}_1^2 \bm{u}_{t-2} + \bm{A}_1 \bm{u}_{t-1} + \bm{u}_t $$
この再帰を無限に繰り返し、安定性条件（$\bm{A}_1$ の固有値がすべて単位円内）の下で $\bm{A}_1^n \to \bm{0}$（$n \to \infty$）を利用すると、VMA($\infty$)表現が得られます。
$$ \bm{y}_t = \sum_{i=0}^{\infty} \bm{\Phi}_i \bm{u}_{t-i}, \quad \bm{\Phi}_i = \bm{A}_1^i, \quad \bm{\Phi}_0 = \bm{I}_k $$
一般のVAR($p$)の場合も、コンパニオン形式に変換すれば同様のVMA表現が得られます。行列 $\bm{\Phi}_i$ の $(j, l)$ 要素は「$i$ 期前の変数 $l$ へのショックが、現在の変数 $j$ にどれだけ影響しているか」を表します。
この $\bm{\Phi}_i$ こそがインパルス応答関数の核心です。VMA表現を手に入れたことで、インパルス応答の理論的な導出に進む準備が整いました。
インパルス応答関数（IRF）の理論と導出
インパルス応答とは何か
インパルス応答関数（Impulse Response Function, IRF）は、元々は電気工学・制御工学で使われていた概念です。電気回路にパルス信号（一瞬だけ電圧をかける）を入力したとき、出力がどう変化するかを追跡する — これがインパルス応答の原型です。時系列分析ではこれを借用し、「ある変数に1単位のショック（撹乱）を与えたとき、システム内の各変数がどのように反応するか」を時間軸に沿って追跡します。
具体的な場面で考えましょう。中央銀行が金利を突然0.25ポイント引き上げたとします。このショックは即座にGDPや失業率に影響するでしょうか？ おそらくすぐには影響せず、数四半期かけて徐々に効果が現れ、やがて元の水準に戻る（あるいは新しい水準に落ち着く）でしょう。IRFは、このプロセスの全体像を描きます。
数学的定義
VMA表現 $\bm{y}_t = \sum_{i=0}^{\infty} \bm{\Phi}_i \bm{u}_{t-i}$ から、インパルス応答関数は次のように定義されます。
$$ \text{IRF}_{j \leftarrow l}(h) = \frac{\partial y_{j,t+h}}{\partial u_{l,t}} = (\bm{\Phi}_h)_{jl} $$
これは「時刻 $t$ に変数 $l$ への撹乱項 $u_{l,t}$ が1単位増加したとき、$h$ 期後の変数 $j$ がどれだけ変化するか」を意味します。$\bm{\Phi}_0 = \bm{I}_k$ なので、ショック直後（$h=0$）には自分自身に1の影響を与え、他の変数への影響は0です（少なくとも非直交化IRFの場合）。
直交化インパルス応答 — コレスキー分解
ここで重要な問題が生じます。現実には、撹乱項 $\bm{u}_t$ の各成分は一般に相関を持っています。つまり $\bm{\Sigma}_u = E[\bm{u}_t \bm{u}_t^\top]$ は対角行列ではありません。
これは解釈上の困難を引き起こします。「変数1にだけショックを与える」つもりでも、$u_{1t}$ と $u_{2t}$ が相関していると、変数2にも同時にショックが波及してしまうのです。「変数1のショックの影響」を純粋に取り出すには、相関を除去する必要があります。
標準的な方法は コレスキー分解 を使うことです。正定値対称行列 $\bm{\Sigma}_u$ を下三角行列 $\bm{P}$ を用いて次のように分解します。
$$ \bm{\Sigma}_u = \bm{P}\bm{P}^\top $$
この $\bm{P}$ を使って新しい撹乱項 $\bm{\varepsilon}_t = \bm{P}^{-1}\bm{u}_t$ を定義すると、
$$ E[\bm{\varepsilon}_t \bm{\varepsilon}_t^\top] = \bm{P}^{-1}\bm{\Sigma}_u (\bm{P}^{-1})^\top = \bm{P}^{-1}\bm{P}\bm{P}^\top(\bm{P}^\top)^{-1} = \bm{I}_k $$
$\bm{\varepsilon}_t$ の各成分は互いに無相関で、分散が1に正規化されています。VMA表現を $\bm{\varepsilon}_t$ を用いて書き直すと、
$$ \bm{y}_t = \sum_{i=0}^{\infty} \bm{\Phi}_i \bm{P} \bm{\varepsilon}_{t-i} = \sum_{i=0}^{\infty} \bm{\Theta}_i \bm{\varepsilon}_{t-i} $$
ここで $\bm{\Theta}_i = \bm{\Phi}_i \bm{P}$ です。直交化インパルス応答関数 は次のように定義されます。
$$ \text{OIRF}_{j \leftarrow l}(h) = (\bm{\Theta}_h)_{jl} = (\bm{\Phi}_h \bm{P})_{jl} $$
この定義では、「変数 $l$ への1標準偏差のショック」が他の変数に純粋にどう伝わるかを追跡できます。
コレスキー分解の順序依存性
コレスキー分解を使う直交化IRFには、重要な注意点があります。それは 変数の順序に結果が依存する ということです。
$\bm{P}$ は下三角行列なので、順序が最初の変数のショックは他のすべての変数に同時に影響しますが、順序が最後の変数のショックは自分自身にしか同時に影響しません。これは「最初の変数が最も外生的（他から影響を受けにくい）」という暗黙の仮定を置いていることになります。
したがって、変数の順序は経済理論や物理的な因果関係に基づいて慎重に決めるべきです。例えばマクロ経済分析では「金利 → GDP → インフレ率」のように、政策変数を先に置く慣行があります。順序に対する頑健性を確認するために、複数の順序でIRFを計算して結果を比較することも推奨されます。
理論の準備が整いました。次は、VARモデルの推定からインパルス応答関数の計算・可視化まで、Pythonで一気通貫に実装しましょう。
インパルス応答関数のPython実装と可視化
ここでは、マクロ経済的な3変量システム（金利・GDP成長率・インフレ率）を模した合成データを生成し、VARモデルを推定した上でインパルス応答関数を計算・可視化します。
データ生成
まず、既知の構造を持つ3変量VAR(2)プロセスを生成します。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
np.random.seed(123)
T = 1000
k = 3 # 変数の数: interest_rate, gdp_growth, inflation
# VAR(2)の係数行列を設定
A1 = np.array([
[0.5, -0.1, 0.0], # 金利は自身の過去に強く依存
[0.1, 0.4, 0.05], # GDP成長率は金利にも影響される
[0.0, 0.15, 0.3] # インフレ率はGDP成長率に影響される
])
A2 = np.array([
[0.1, 0.0, 0.0],
[-0.05, 0.1, 0.0],
[0.0, 0.05, 0.1]
])
# 撹乱項の共分散行列（変数間にいくらかの同時相関を持たせる）
Sigma_u = np.array([
[1.0, 0.3, 0.1],
[0.3, 0.8, 0.2],
[0.1, 0.2, 0.6]
])
L_chol = np.linalg.cholesky(Sigma_u)
# VAR(2)プロセスをシミュレーション
Y = np.zeros((T, k))
for t in range(2, T):
eps = L_chol @ np.random.randn(k)
Y[t] = A1 @ Y[t-1] + A2 @ Y[t-2] + eps
# データに名前をつける
var_names = ['Interest Rate', 'GDP Growth', 'Inflation']
VARモデルの推定とIRFの計算
import pandas as pd
# pandasのDataFrameに変換
df = pd.DataFrame(Y, columns=var_names)
# VARモデルの推定（ラグ次数はAICで選択）
model = VAR(df)
lag_order = model.select_order(maxlags=8)
print("AIC, BIC, HQIC, FPEによるラグ次数の選択:")
print(lag_order.summary())
# AICで選択されたラグ次数で推定
results = model.fit(lag_order.aic)
print(f"\n選択されたラグ次数 (AIC): {lag_order.aic}")
print(f"推定された係数行列 A1:\n{results.coefs[0].round(3)}")
この出力では、AICが正しくラグ2を選択し、推定された係数行列 $\hat{\bm{A}}_1$ がデータ生成に使った真の $\bm{A}_1$ に近い値を示すことが確認できます。推定値と真値が完全に一致しないのは、有限サンプルの推定誤差によるものであり、$T = 1000$ のサンプルサイズでは十分に精度の高い推定が得られます。
インパルス応答関数の可視化
# インパルス応答関数を計算（直交化IRF、20期先まで）
irf = results.irf(20)
# 可視化
fig, axes = plt.subplots(3, 3, figsize=(14, 10))
for i in range(k):
for j in range(k):
ax = axes[i, j]
irf_values = irf.irfs[:, i, j]
lower = irf.ci[:, i, j, 0] # 信頼区間の下限
upper = irf.ci[:, i, j, 1] # 信頼区間の上限
periods = np.arange(len(irf_values))
ax.plot(periods, irf_values, 'b-', linewidth=2)
ax.fill_between(periods, lower, upper, alpha=0.2, color='blue')
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_title(f'{var_names[j]} → {var_names[i]}', fontsize=10)
ax.grid(True, alpha=0.3)
if i == 2:
ax.set_xlabel('Periods')
plt.suptitle('Orthogonalized Impulse Response Functions', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
3行3列のIRFグラフから、複数の重要な知見が得られます。
- 金利 → GDP成長率（1行目・2列目）: 金利への正のショックはGDP成長率を最初は押し下げ、2〜3期後に最大の負の影響を与えた後、徐々に0に収束します。これは金融引き締めの効果がタイムラグを持って実体経済に影響するという経済学的直感と整合します。
- GDP成長率 → インフレ率（2行目・3列目）: GDP成長率へのショックはインフレ率を一時的に上昇させます。需要の増加が物価上昇を引き起こすというディマンドプル・インフレのメカニズムが見えています。
- 対角要素（自己応答）: 各変数は自分自身のショックに対して時刻0で最大の応答を示し、その後指数的に減衰します。安定なVARプロセスなので、すべてのIRFは最終的に0に収束しています。
- 信頼区間: 青い帯で示された95%信頼区間を見ると、金利→GDP成長率の効果は統計的に有意（信頼区間が0をまたがない期間が存在する）ですが、一部の弱い効果では信頼区間が0をまたいでおり、統計的に有意とは言えません。
IRFは「ショックがどう伝わるか」というパターンを教えてくれましたが、次の自然な問いは「各変数の変動のうち、どの程度がどのショックに起因するか？」です。この定量的な帰属分析を行うのが分散分解です。
分散分解（FEVD）の理論
変動の源泉を分解する
インパルス応答関数は「パターン」を描きますが、分散分解（Forecast Error Variance Decomposition, FEVD）は「大きさ」を数値化します。日常の例で言えば、IRFが「台風が来ると翌日に気温が下がり、3日後に湿度が上がる」というパターンを示すのに対し、FEVDは「気温の変動のうち60%は台風に起因し、30%は季節変動、10%は測定誤差」のように割合を教えてくれます。
数学的導出
VMA表現 $\bm{y}_t = \sum_{i=0}^{\infty} \bm{\Theta}_i \bm{\varepsilon}_{t-i}$（直交化済み）を使います。$h$ 期先予測誤差は次のように書けます。
$$ \bm{y}_{t+h} – \hat{\bm{y}}_{t+h|t} = \sum_{i=0}^{h-1} \bm{\Theta}_i \bm{\varepsilon}_{t+h-i} $$
この予測誤差の $j$ 番目の成分を取り出すと、
$$ y_{j,t+h} – \hat{y}_{j,t+h|t} = \sum_{i=0}^{h-1} \sum_{l=1}^{k} (\bm{\Theta}_i)_{jl} \varepsilon_{l,t+h-i} $$
$\bm{\varepsilon}_t$ の各成分は互いに無相関で分散1なので、$j$ 番目の変数の $h$ 期先予測誤差分散は次のようになります。
$$ \text{MSE}_j(h) = \text{Var}(y_{j,t+h} – \hat{y}_{j,t+h|t}) = \sum_{i=0}^{h-1} \sum_{l=1}^{k} (\bm{\Theta}_i)_{jl}^2 $$
この合計を「変数 $l$ のショックに起因する部分」と「それ以外」に分けることで、分散分解が得られます。変数 $j$ の $h$ 期先予測誤差分散のうち、変数 $l$ のショックに帰属する割合は次のとおりです。
$$ \text{FEVD}_{j \leftarrow l}(h) = \frac{\sum_{i=0}^{h-1} (\bm{\Theta}_i)_{jl}^2}{\sum_{i=0}^{h-1} \sum_{m=1}^{k} (\bm{\Theta}_i)_{jm}^2} $$
分子はインパルス応答係数の二乗和なので、IRFの「大きさ」を累積的に集計していることが見て取れます。すべての変数 $l$ について足し上げると $\sum_{l=1}^{k} \text{FEVD}_{j \leftarrow l}(h) = 1$ となり、予測誤差分散が100%分配されます。
分散分解の時間的変化
FEVDは予測ホライズン $h$ に依存します。一般的な傾向として、
- 短期（$h$ が小さい）: 各変数の予測誤差分散は自分自身のショックに大きく帰属します。なぜなら、他の変数からの影響が伝搬するには時間がかかるためです
- 長期（$h$ が大きい）: 他の変数のショックの寄与が増えていきます。VARモデルの連立構造を通じて、ショックが徐々にシステム全体に浸透するためです
この時間的変化のパターン自体が、変数間の関係性について重要な情報を提供します。では、実際にPythonで分散分解を計算し、可視化しましょう。
分散分解のPython実装
# 分散分解を計算（20期先まで）
fevd = results.fevd(20)
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
colors = ['#2196F3', '#FF9800', '#4CAF50']
for i in range(k):
ax = axes[i]
decomp = fevd.decomp[i] # (h, k)の配列
periods = np.arange(1, decomp.shape[0] + 1)
# 積み上げ面グラフ
ax.stackplot(periods, decomp.T, labels=var_names, colors=colors, alpha=0.8)
ax.set_title(f'FEVD of {var_names[i]}', fontsize=12)
ax.set_xlabel('Forecast Horizon')
ax.set_ylabel('Proportion')
ax.set_ylim(0, 1)
ax.legend(loc='center right', fontsize=8)
ax.grid(True, alpha=0.3)
plt.suptitle('Forecast Error Variance Decomposition', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
# 数値でも確認（10期先と20期先）
print("=== 分散分解 (h=10) ===")
fevd.summary()
積み上げ面グラフから、いくつかの特徴的なパターンが読み取れます。
- 金利の予測誤差: ほとんどが自分自身のショックに帰属しており、GDP成長率やインフレ率からの影響は限定的です。これは金利が比較的「外生的」な変数であることを意味しています（データ生成時に $\bm{A}_1$ の1行目で $a_{12}$, $a_{13}$ を小さく設定したことと整合します）。
- GDP成長率の予測誤差: 短期（1〜3期先）では自身のショックが支配的ですが、ホライズンが伸びるにつれて金利のショックの寄与が徐々に増加しています。10期先では金利のショックが全体の15〜20%程度を占めるようになります。金融政策の波及効果がタイムラグを持って蓄積されていく様子が定量的に見えています。
- インフレ率の予測誤差: GDP成長率と自身のショックが主な源泉です。金利のショックの寄与は比較的小さく、金利がインフレに影響するのはGDP成長率を経由した間接的な経路が中心であることを示唆しています。
ここまでの分析ツール — CCF、グレンジャー因果、IRF、FEVD — はすべて 定常な 時系列を前提としています。しかし、経済データの多くは非定常（単位根を持つ）です。非定常な系列間に長期的な均衡関係が存在する場合、階差を取って定常化すると長期情報が失われてしまいます。この問題に対処するのが共和分の概念です。
共和分（Cointegration）の導入
「一緒にさまよう」2つの時系列
共和分の概念は、最初は直感に反するかもしれません。まず、身近な例で考えてみましょう。
酔っ払いが2人、同じ居酒屋を出て帰路につきます。それぞれランダムウォーク（千鳥足）で歩くので、個々の位置は予測不可能です。しかし、2人がロープで結ばれているとしたらどうでしょう？ それぞれの位置は依然としてランダムウォークですが、2人の 距離 は一定の範囲内に保たれます。この「個々には非定常だが、線形結合が定常」という関係が共和分です。
数学的定義
$k$ 個の時系列 $y_{1t}, y_{2t}, \ldots, y_{kt}$ がそれぞれ $I(1)$（1回の階差で定常になる）であるとします。あるゼロでないベクトル $\bm{\beta} = (\beta_1, \beta_2, \ldots, \beta_k)^\top$ が存在して、
$$ z_t = \bm{\beta}^\top \bm{y}_t = \beta_1 y_{1t} + \beta_2 y_{2t} + \cdots + \beta_k y_{kt} \sim I(0) $$
すなわち線形結合 $z_t$ が定常過程であるとき、$\bm{y}_t$ は 共和分（cointegrated） であると言い、$\bm{\beta}$ を 共和分ベクトル と呼びます。
経済学的な解釈
共和分関係は 長期均衡 を表します。有名な例として、消費と所得の関係があります。消費 $C_t$ と所得 $Y_t$ はそれぞれ非定常（成長トレンドを持つ）ですが、恒常所得仮説によれば両者には長期的な比例関係 $C_t = \beta Y_t + z_t$ があり、乖離 $z_t$ は定常（一時的な逸脱は修正される）です。
他にも、株価と配当、為替レートと金利差、長期金利と短期金利（金利の期間構造仮説）など、多くの経済変数のペアに共和分関係が見出されています。
ヨハンセン検定の概要
共和分関係の検定と推定には、ヨハンセン検定（Johansen test） が標準的に使われます。この方法はVECM（Vector Error Correction Model）を基盤とし、以下のステップで進みます。
まず、VAR($p$)モデルをVECM形式に書き換えます。
$$ \Delta \bm{y}_t = \bm{\Pi} \bm{y}_{t-1} + \sum_{i=1}^{p-1} \bm{\Gamma}_i \Delta \bm{y}_{t-i} + \bm{u}_t $$
ここで $\bm{\Pi} = \sum_{i=1}^{p} \bm{A}_i – \bm{I}_k$ は 長期影響行列 です。共和分の検定は、この $\bm{\Pi}$ のランク $r$ を検定することに帰着します。
- $\text{rank}(\bm{\Pi}) = 0$: 共和分関係なし（各変数は独立にランダムウォーク）
- $0 < \text{rank}(\bm{\Pi}) = r < k$: $r$ 本の共和分関係が存在
- $\text{rank}(\bm{\Pi}) = k$: すべての変数が定常（そもそも共和分の枠組みは不要）
$\bm{\Pi}$ のランクは $\bm{\Pi} = \bm{\alpha}\bm{\beta}^\top$ と分解されます。ここで $\bm{\beta}$ は $k \times r$ の共和分ベクトル行列、$\bm{\alpha}$ は $k \times r$ の調整速度行列です。$\bm{\beta}^\top \bm{y}_{t-1}$ は均衡からの乖離を表し、$\bm{\alpha}$ はその乖離がどの程度の速さで修正されるかを表します。
ヨハンセンはこのランクを、行列 $\bm{\Pi}$ の固有値に基づくトレース検定と最大固有値検定の2種類で検定する手法を開発しました。Pythonではstatsmodelsの coint_johansen 関数で実行できます。
from statsmodels.tsa.vector_ar.vecm import coint_johansen
# 共和分検定のデモ用データ生成
# 2つのI(1)系列が共和分関係を持つ例
np.random.seed(456)
T_coint = 500
# 共通の確率トレンド
trend = np.cumsum(np.random.randn(T_coint))
# 共和分関係を持つ2系列
x1 = trend + 0.5 * np.random.randn(T_coint) # x1 ≈ trend
x2 = 2.0 * trend + 0.8 * np.random.randn(T_coint) # x2 ≈ 2*trend
data_coint = np.column_stack([x1, x2])
# ヨハンセン検定（定数項あり、ラグ2）
result_joh = coint_johansen(data_coint, det_order=0, k_ar_diff=2)
print("=== ヨハンセン検定の結果 ===")
print(f"トレース統計量: {result_joh.lr1.round(3)}")
print(f"臨界値 (90%, 95%, 99%):\n{result_joh.cvt.round(3)}")
print(f"\n最大固有値統計量: {result_joh.lr2.round(3)}")
print(f"臨界値 (90%, 95%, 99%):\n{result_joh.cvm.round(3)}")
print(f"\n推定された共和分ベクトル: {result_joh.evec[:, 0].round(3)}")
出力を見ると、ランク0の帰無仮説（共和分関係なし）は棄却され、ランク1の帰無仮説は棄却されません。したがって、共和分関係は1本存在すると判定されます。推定された共和分ベクトルは $(1, -0.5)$ に近い値を示すはずです。なぜなら、データ生成時に $x_2 \approx 2 x_1$ としたので、$x_1 – 0.5 x_2 \approx 0$（定常）という関係が成り立つからです。
共和分が検出された場合は、VARモデルの代わりにVECM（誤差修正モデル）を使うことで、短期ダイナミクスと長期均衡の両方を適切にモデル化できます。VECMの詳細は別の記事で扱いますが、ここでは「非定常系列間の長期関係を検出する道具がある」ということを押さえておきましょう。
ここまでで5つの分析ツールを一通り見てきました。最後に、実務でどの手法をどの場面で使うべきかを整理しましょう。
手法の使い分けガイド
分析フローチャート
実際に多変量時系列データを分析する際には、以下のフローで手法を選択すると体系的に進められます。
ステップ1: 定常性の確認
まず、各系列にADF検定やKPSS検定を適用して定常性を調べます。
- すべて定常 → ステップ2へ
- 非定常系列が含まれる → ステップ1aへ
ステップ1a: 共和分の検定
非定常系列が複数ある場合、ヨハンセン検定で共和分関係を調べます。
- 共和分あり → VECM（誤差修正モデル）を使用。その上でIRFやFEVDを実行
- 共和分なし → 階差を取って定常化し、ステップ2へ
ステップ2: 関係性の探索的分析
CCFを計算してリード・ラグ構造を可視化します。どの変数がどの変数に先行しているかの大まかなアタリをつけます。
ステップ3: 因果関係の検定
グレンジャー因果検定を全変数のペアに対して実行し、予測的因果の方向を特定します。この結果はVARモデルの変数選択にも使えます（因果関係のない変数をモデルから除外する根拠になります）。
ステップ4: VARモデルの推定
適切なラグ次数を情報量基準で選択し、VARモデルを推定します。残差診断（自己相関の不在、正規性など）を確認します。
ステップ5: インパルス応答と分散分解
推定されたVARモデルに基づいて、IRFとFEVDを計算します。IRFでショックの動的な波及パターンを可視化し、FEVDで各変数の変動の源泉を定量化します。
各手法の比較表
それぞれの手法が「何に答えるか」と「何を仮定するか」を表にまとめます。
| 手法 | 答える問い | 必要な仮定 | モデル | 方向性 |
|---|---|---|---|---|
| CCF | 時間ずれ付きの相関は？ | 定常性 | なし（ノンパラメトリック） | 双方向を同時に表示 |
| グレンジャー因果 | 予測的因果は？ | 定常性、ラグ次数 | AR/VAR | 一方向ずつ検定 |
| IRF | ショックの動的波及は？ | 定常性、VARの正しい特定 | VAR | 直交化で方向を分離 |
| FEVD | 変動の源泉の割合は？ | 定常性、VARの正しい特定 | VAR | 各変数への帰属割合 |
| 共和分 | 長期均衡関係は？ | $I(1)$の変数 | VECM | 均衡からの調整速度 |
実践的な注意点
サンプルサイズ: VARモデルのパラメータ数は変数数とラグ次数に対して急速に増えます。$k$ 変量VAR($p$) では各方程式に $kp + 1$ 個のパラメータがあり、全体では $k(kp + 1)$ 個です。3変量VAR(4)でも39個のパラメータを推定する必要があり、信頼性のある推定には最低でも200〜300個の観測が欲しいところです。
変数の選択: VARモデルにあまりに多くの変数を含めると、パラメータ数が爆発して推定精度が低下します（次元の呪い）。理論的に重要な変数に絞り、グレンジャー因果検定の結果も参考にして変数を選択しましょう。
構造変化: 分析期間中に構造変化（政策レジームの転換、経済危機など）がある場合、VARモデルの係数が不安定になります。部分サンプル分析やChow検定で構造変化の有無を確認することが重要です。
頑健性チェック: 直交化IRFの結果が変数の順序に依存する問題を回避するために、一般化IRF（Pesaran and Shin, 1998）を計算する方法もあります。statsmodelsでは orth=False オプションで非直交化IRFを計算でき、これは順序に依存しません。
統合的な実装例 — 5つの手法を一気通貫で適用
最後に、これまで学んだ手法を1つのデータセットに対して体系的に適用する統合的な実装例を示します。ここでは、3変量の合成データに対して、CCF → グレンジャー因果 → VAR推定 → IRF → 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, adfuller, ccf
np.random.seed(2026)
# --- 1. データ生成 (3変量VAR(1)プロセス) ---
T = 600
A = np.array([
[0.6, 0.0, 0.0],
[0.3, 0.5, 0.0],
[0.0, 0.2, 0.4]
])
Y = np.zeros((T, 3))
for t in range(1, T):
Y[t] = A @ Y[t-1] + np.random.randn(3)
df = pd.DataFrame(Y, columns=['X1', 'X2', 'X3'])
# --- 2. 定常性の確認 (ADF検定) ---
print("=== ADF検定 ===")
for col in df.columns:
adf_stat, pval, _, _, _, _ = adfuller(df[col])
print(f"{col}: ADF統計量={adf_stat:.3f}, p値={pval:.4f} → {'定常' if pval < 0.05 else '非定常'}")
# --- 3. 相互相関関数 ---
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
pairs = [('X1', 'X2'), ('X1', 'X3'), ('X2', 'X3')]
conf = 2 / np.sqrt(T)
for idx, (a, b) in enumerate(pairs):
ccf_vals = ccf(df[a], df[b], nlags=15)
axes[idx].bar(range(len(ccf_vals)), ccf_vals, width=0.4, color='steelblue', alpha=0.8)
axes[idx].axhline(y=conf, color='red', linestyle='--', alpha=0.7)
axes[idx].axhline(y=-conf, color='red', linestyle='--', alpha=0.7)
axes[idx].axhline(y=0, color='black', linewidth=0.5)
axes[idx].set_title(f'CCF: {a} leading {b}')
axes[idx].set_xlabel('Lag')
axes[idx].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
CCFのグラフから、X1→X2の方向で正の相互相関（ラグ1付近でピーク）が見られ、X2→X3でも同様の構造が観察されます。一方、X1→X3の相互相関は弱く、間接的な影響（X1→X2→X3）が時間をかけて現れることがわかります。この結果は、データ生成の構造（$\bm{A}$ の非対角要素の配置）と整合しています。
# --- 4. グレンジャー因果検定 ---
print("\n=== グレンジャー因果検定 (ラグ=4) ===")
for i, col_y in enumerate(df.columns):
for j, col_x in enumerate(df.columns):
if i != j:
data_pair = df[[col_y, col_x]].values
result = grangercausalitytests(data_pair, maxlag=4, verbose=False)
min_pval = min([result[lag][0]['ssr_ftest'][1] for lag in range(1, 5)])
direction = "***" if min_pval < 0.01 else "**" if min_pval < 0.05 else "*" if min_pval < 0.1 else ""
print(f"{col_x} → {col_y}: 最小p値={min_pval:.4f} {direction}")
グレンジャー因果検定の結果は、X1→X2、X2→X3の方向で有意（p値 < 0.01）であり、逆方向は非有意となります。これはデータ生成の係数行列 $\bm{A}$ で $a_{21} = 0.3$（X1→X2の影響）と $a_{32} = 0.2$（X2→X3の影響）を設定し、逆方向はゼロとした構造と正確に一致しています。
# --- 5. VARモデル推定 → IRF → FEVD ---
model = VAR(df)
results = model.fit(ic='aic')
print(f"\n=== VARモデル (ラグ={results.k_ar}) ===")
# IRF
irf = results.irf(15)
fig = irf.plot(orth=True)
fig.suptitle('Orthogonalized IRF', fontsize=14)
plt.tight_layout()
plt.show()
# FEVD
fevd = results.fevd(15)
fig = fevd.plot()
fig.suptitle('Forecast Error Variance Decomposition', fontsize=14)
plt.tight_layout()
plt.show()
IRFとFEVDの結果は、ここまでのCCFやグレンジャー因果の結果と整合した全体像を描いています。X1のショックはX2に1期遅れで伝播し、X3にはX2を経由して2期以上遅れて伝播します。FEVDでは、X2の予測誤差分散のうちかなりの割合（長期的には30%前後）がX1のショックに帰属しており、X3にもX2を経由した間接的な影響が見えています。5つの手法が一貫したストーリーを語っていることが確認できます。
まとめ
本記事では、多変量時系列の関係性を分析する5つの手法を統一的に解説しました。
- 相互相関関数（CCF） は、2つの時系列間のリード・ラグ構造を探索的に検出するための最もシンプルなツールです。モデルを仮定せずに「どの時間ずれで相関が強いか」を可視化できます
- グレンジャー因果検定 は、「ある変数の過去の情報が別の変数の予測に寄与するか」をF検定で判定します。CCFでは捉えられない予測的因果の方向性を統計的に検定できます
- インパルス応答関数（IRF） は、VARモデルの上に構築され、ある変数への突発的ショックが他の変数にどのような時間的パターンで波及するかを追跡します。コレスキー分解による直交化で、純粋なショックの影響を分離できます
- 分散分解（FEVD） は、各変数の予測誤差分散をシステム内の各ショックに帰属させ、変動の源泉を定量的に分解します。IRFが「パターン」を示すのに対し、FEVDは「大きさ」を数値化します
- 共和分 は、非定常な時系列間に長期均衡関係が存在するかを検定し、VECMの枠組みで短期ダイナミクスと長期均衡を同時にモデル化する道を開きます
これらの手法は排他的ではなく相補的であり、多角的な視点から関係性を検証することで、より頑健な結論が得られます。分析のフローとしては、定常性の確認 → 共和分検定（必要に応じて）→ CCFで探索 → グレンジャー因果で方向性を検定 → VAR/VECMを推定 → IRFとFEVDで動的構造を分析、という順序が標準的です。
次のステップとして、以下の記事も参考にしてください。