Granger Causality and Impulse-Response Analysis

"Yes, the rooster's crow predicts the dawn. No, I am not claiming the rooster lifts the sun. I only ever said that knowing the crow improves my forecast of the morning, and I will defend that modest sentence in any court that does not call itself a laboratory."

A Granger Test Carefully Not Claiming Real Causation

The previous section of this chapter built the vector autoregression: a system in which each of $K$ series is regressed on $p$ lags of every series, itself included, so that the cross-series feedback the univariate models of Chapter 5 could never see becomes the central object. A fitted VAR is a dense bundle of coefficients, and on its own it is hard to interpret. This section supplies the two standard lenses that turn that bundle into statements a practitioner can act on. The first lens, Granger causality, asks a yes-or-no question about predictive content: does the past of one series carry information about the future of another beyond what that other series already knows about itself. The second lens, the impulse-response function, asks a quantitative question about propagation: if one series is hit by a one-time unit shock, what happens to every series in the system over the next several periods. Throughout we use the multivariate notation $\mathbf{y}_t$, the coefficient matrices $\mathbf{A}_i$, and the innovation covariance $\mathbf{\Sigma}_u$ collected in the unified notation table of Appendix A, and we lean on the stationarity machinery of Section 3.2 at every turn.

By the end of this section you should be able to do five concrete things. First, state Granger's definition precisely and explain why it is a claim about prediction and never about intervention. Second, run the Granger $F$-test as a nested-model comparison and read its verdict in both directions. Third, list the four assumptions the test rests on and recognize when each is violated. Fourth, derive the impulse-response function from the moving-average form of a VAR and read an IRF plot for stationarity, peak, and decay. Fifth, orthogonalize an IRF with the Cholesky decomposition, compute a forecast-error variance decomposition, and state out loud the ordering assumption both require. The worked example of subsection five exercises all five against a system whose answer we know in advance.

1. Granger Causality Defined Beginner

Clive Granger gave the word "causality" an operational, testable meaning that deliberately sidesteps the philosophical thicket of true causation. His definition is entirely about prediction. A series $X$ Granger-causes a series $Y$ if the past of $X$ improves the forecast of $Y$ beyond what the past of $Y$ alone achieves. Put formally, let $\mathcal{I}_t$ be all information up to time $t$ and let $\mathcal{I}_t^{-X}$ be that same information with the history of $X$ removed. Then $X$ Granger-causes $Y$ when

The forecast of $Y$ gets strictly sharper (its error variance falls) when you are allowed to look at the past of $X$. That is the whole idea: causality, in Granger's sense, is a measurable reduction in prediction error. It says nothing about whether wiggling $X$ in the world would move $Y$. The rooster's crow Granger-causes the dawn because the crow precedes and predicts the sunrise, yet silencing the rooster leaves the sun entirely unmoved. We label this relationship "causality" only by Granger's narrow operational convention, and the entire interventional theory of when prediction may be promoted to genuine cause is deferred to Chapter 30, where the do-calculus and structural causal models of the causality chapter draw the line precisely.

Two structural conditions are baked into the definition and are easy to miss. The first is the temporal ordering: the conditioning set $\mathcal{I}_t$ stops at time $t$, so only the strict past of $X$ is allowed to help predict $Y_{t+1}$. This is what makes the definition directional and what excludes same-period (contemporaneous) effects from the test entirely, a gap we return to as instantaneous causality. The second is the conditioning on $Y$'s own past: the bar is not "does $X$ predict $Y$" but "does $X$ predict $Y$ beyond what $Y$ already predicts about itself." A series whose own history is highly informative sets a high bar, and an apparent $X \to Y$ link can evaporate once $Y$'s autoregressive structure is properly credited. The whole construction is, at bottom, a comparison of two conditional forecast-error variances, and everything operational follows from making that comparison a hypothesis test.

To make the definition testable we attach it to a VAR. Consider a bivariate VAR($p$) in $(X_t, Y_t)$, and write the equation for $Y_t$ alone, which regresses $Y$ on $p$ lags of itself and $p$ lags of $X$,

In this equation, $X$ fails to Granger-cause $Y$ exactly when every coefficient on lagged $X$ vanishes, $\beta_1 = \beta_2 = \cdots = \beta_p = 0$. That is a linear restriction on the parameters, so the natural test is the nested-model $F$-test (equivalently a Wald or likelihood-ratio test). Fit the unrestricted model above and record its residual sum of squares $\mathrm{RSS}_U$; fit the restricted model that drops the $X$ lags, $Y_t = c + \sum_i \alpha_i Y_{t-i} + u_t$, and record $\mathrm{RSS}_R$. The $F$-statistic comparing them is

The numerator measures how much the $X$ lags reduce the residual sum of squares (the forecast improvement at the heart of the definition); the denominator normalizes by the unrestricted error and the degrees of freedom. A large $F$, equivalently a small $p$-value, rejects the null and says the $X$ lags carry predictive content: $X$ Granger-causes $Y$. The test is asymmetric by construction. You run it again with the roles swapped (regress $X$ on lags of $X$ and $Y$, restrict the $Y$ lags) to ask whether $Y$ Granger-causes $X$, and the two answers are independent: you can find one direction, both directions (feedback), or neither. The four possible outcomes carve a small taxonomy that is worth naming, because each corresponds to a different picture of how the two series relate in time.

The independence row carries a quiet warning that subsection two will sharpen: a "retain, retain" verdict rules out only linear lagged predictability. Two series can fail every Granger test and still be tightly coupled through a contemporaneous correlation or a nonlinear dependence the linear $F$-test cannot see. The test is a statement about one specific, linear, lagged channel of information, and reading more into a non-rejection than that is its own pitfall. The numeric example below shows the arithmetic of the $F$-test on small sums of squares so the formula stops being abstract.

2. Testing Procedure and Pitfalls Intermediate

The $F$-test of subsection one is mechanically simple, which is exactly why it is so often run incorrectly. Four assumptions sit underneath it, and each is a pitfall when violated. The procedure, done properly, is a short checklist. First, ensure stationarity. The asymptotic $F$ (or $\chi^2$) distribution of the Granger statistic is derived under stationarity, and on non-stationary series the test is unreliable: it produces spurious rejections at rates far above the nominal level, the multivariate echo of the spurious-regression problem of Section 3.2. So you first difference or otherwise transform each series to stationarity (checking with the augmented Dickey-Fuller test), unless the series are cointegrated, in which case the differencing throws away the long-run relationship and you must instead use the error-correction framework of Section 6.3 and its Toda-Yamamoto augmented-lag variant of the Granger test. Second, choose the lag order $p$ deliberately. Too few lags and genuine predictive structure at longer horizons is missed (the test loses power); too many lags and the test wastes degrees of freedom and can manufacture significance from noise. The standard practice is to set $p$ by an information criterion (AIC or BIC) on the full VAR, exactly as the previous section selected the VAR order, and then run the Granger test at that $p$. In practice the most common version of this error is eyeballing $p$ from a single correlogram; always fix $p$ by AIC or BIC on the full VAR (as in Section 6.1) and run the test at that $p$, so the lag order is not a free parameter you tune until the result looks significant.

Third, beware confounders and common drivers. This is the pitfall that connects directly to the interventional warning of subsection one. If an unobserved third series $Z$ drives both $X$ and $Y$ but reaches $X$ sooner, then $X$ will Granger-cause $Y$ with no mechanism linking them at all; the test is faithfully reporting predictive content that flows entirely through the hidden common cause. Adding $Z$ to the system (making it a trivariate VAR and testing $X \to Y$ conditional on $Z$) can dissolve the spurious link, which is why conditional Granger causality, run inside a larger VAR, is more trustworthy than the bivariate version. Fourth, separate instantaneous from lagged causality. The $F$-test only examines lagged coefficients, so it is silent about contemporaneous dependence, the correlation between the same-period innovations $u_t^X$ and $u_t^Y$ captured by the off-diagonal of $\mathbf{\Sigma}_u$. When the innovations are contemporaneously correlated we say there is instantaneous causality, and it is a genuinely different and direction-ambiguous phenomenon: with sampled data you cannot tell from correlation alone whether the within-period link runs $X \to Y$ or $Y \to X$. That ambiguity is precisely the identification problem that subsection four must confront for impulse responses.

Two of these pitfalls deserve a sharper, operational statement because they cause the most published errors. The stationarity pitfall is not a minor technicality: on two independent random walks, the bivariate Granger test rejects the (true) null of no causality far more than the nominal $5\%$ of the time, sometimes above $30\%$, because the regression of one integrated series on another inherits the spurious-regression pathology of Section 3.2; the $t$ and $F$ statistics no longer have their textbook distributions. The fix is mechanical (test for and remove unit roots first) but the failure is silent, which is why so many "causal" findings on raw price levels are artifacts. The cointegration carve-out matters too: if the series share a common stochastic trend, naive differencing destroys the long-run relationship and again distorts the test, and the clean remedy is the Toda-Yamamoto procedure, which fits the VAR in levels with one extra lag beyond the chosen $p$ and applies the Wald test only to the first $p$ lag-blocks, restoring a valid $\chi^2$ distribution even under integration and cointegration without requiring you to pre-test the integration order perfectly. The error-correction machinery of Section 6.3 gives the principled treatment.

The instantaneous-causality pitfall has its own operational tell. Because the lagged-coefficient $F$-test cannot see the contemporaneous correlation in $\mathbf{\Sigma}_u$, a strong same-period link between two series can leave both Granger directions non-significant while the series are in fact tightly coupled within the period. The diagnostic is to inspect the off-diagonal of the estimated innovation covariance: a large off-diagonal alongside two non-rejections is the fingerprint of purely instantaneous causality, a relationship whose direction the sampled data cannot resolve and which only the identifying assumptions of subsection four can even tentatively orient.

3. Impulse-Response Functions Intermediate

Granger causality gives a binary verdict; the impulse-response function (IRF) gives the dynamic shape. It answers the question a forecaster actually cares about: if the system receives a one-time shock to one variable today, how do all the variables respond over the coming periods. The IRF is read directly off the moving-average representation of the VAR, so we derive that first. A stationary VAR($p$), written compactly with the lag-operator matrix polynomial $\mathbf{A}(L) = \mathbf{I} - \mathbf{A}_1 L - \cdots - \mathbf{A}_p L^p$, is $\mathbf{A}(L)\,\mathbf{y}_t = \mathbf{u}_t$ with innovations $\mathbf{u}_t$ of covariance $\mathbf{\Sigma}_u$. Because the process is stationary, $\mathbf{A}(L)$ is invertible, and inverting it expresses the present as an infinite weighted sum of past innovations, the vector moving-average (VMA) form,

This inversion is exact only because stationarity guarantees it. The condition, the multivariate generalization of the scalar unit-circle rule from Section 5.3, is that every root of $\det\!\big(\mathbf{A}(z)\big) = 0$ lies strictly outside the unit circle, equivalently that every eigenvalue of the companion matrix of the VAR lies strictly inside it. When that holds, the $\mathbf{\Phi}_h$ decay to zero and the infinite sum converges; when an eigenvalue sits on the circle, the series carries a unit root and the inversion diverges, which is the formal reason a Granger or IRF analysis on a non-stationary VAR is invalid rather than merely imprecise. Each $\mathbf{\Phi}_h$ is a $K \times K$ matrix, and its $(j,k)$ entry is exactly the impulse response we want: the effect on variable $j$ at horizon $h$ of a one-unit shock to the innovation of variable $k$ at horizon zero. The matrices are computed by a simple recursion that unrolls the inversion, with the convention $\mathbf{A}_i = \mathbf{0}$ for $i > p$,

So the entire IRF is generated from the fitted VAR coefficient matrices by this matrix recursion, starting from $\mathbf{\Phi}_0 = \mathbf{I}$. Reading an IRF plot is then straightforward once you know what the axes mean. The standard display is a grid of small panels, one per ordered pair (response variable, shocked variable). The horizontal axis is the horizon $h$ in periods after the shock; the vertical axis is the response magnitude; the curve traces $[\mathbf{\Phi}_h]_{jk}$ as $h$ grows, usually with a shaded confidence band from a bootstrap. A response that rises, peaks, and decays back to zero is the signature of a stationary system absorbing a transient shock; a response that fails to return to zero signals a unit root and a misspecified (non-stationary) VAR. The table below catalogs the shapes you will actually meet, so reading a panel becomes a matter of pattern recognition rather than guesswork.

The last row is the discipline a bootstrap band enforces: a point estimate that looks like a clean hump means little if its confidence band contains zero throughout, and reputable IRF reporting always shows the band. The numeric example below turns the recursion into actual numbers on a tiny VAR(1) so the propagation is visible.

The recursion is short enough to implement directly, and doing so demystifies what the library computes. Code 6.2.0 generates the raw (non-orthogonalized) impulse-response matrices $\mathbf{\Phi}_h$ from a list of VAR coefficient matrices, reproducing the by-hand VAR(1) numbers of the example above and showing the geometric decay that signals stationarity.

import numpy as np

def irf_matrices(A_list, horizon):
    """Raw (non-orthogonalized) impulse-response matrices Phi_0..Phi_H from VAR coeff matrices.
    A_list = [A_1, ..., A_p], each K x K. Implements Phi_h = sum_i A_i Phi_{h-i}, Phi_0 = I."""
    p = len(A_list)
    K = A_list[0].shape[0]
    Phi = [np.eye(K)]                                  # Phi_0 = I by definition
    for h in range(1, horizon + 1):
        acc = np.zeros((K, K))
        for i in range(1, min(h, p) + 1):              # A_i = 0 for i > p, so cap at p
            acc += A_list[i - 1] @ Phi[h - i]          # accumulate the recursion terms
        Phi.append(acc)
    return Phi

A1 = np.array([[0.5, 0.2],
               [0.0, 0.4]])                            # the worked VAR(1) coefficient matrix
Phi = irf_matrices([A1], horizon=4)
print("Response of variable 1 to a unit shock in variable 2, the (1,2) entry of each Phi_h:")
for h, M in enumerate(Phi):
    print(f"  h={h}: {M[0, 1]:.3f}")

Response of variable 1 to a unit shock in variable 2, the (1,2) entry of each Phi_h:
  h=0: 0.000
  h=1: 0.200
  h=2: 0.180
  h=3: 0.122
  h=4: 0.078

4. Identification and Orthogonalization Advanced

There is a serious problem hiding in the IRF of subsection three, and ignoring it produces nonsense. The shocks $\mathbf{u}_t$ are generally contemporaneously correlated: $\mathbf{\Sigma}_u$ has nonzero off-diagonals, which is exactly the instantaneous causality of subsection two. A "unit shock to variable $k$" is therefore not a clean experiment, because in the data a shock to $k$ never arrives alone; it comes bundled with a correlated same-period shock to the other variables. The raw IRF of $\mathbf{\Phi}_h$ thus answers a counterfactual that the data never realizes. To get interpretable responses we must first orthogonalize the shocks, transforming the correlated $\mathbf{u}_t$ into uncorrelated structural shocks $\mathbf{w}_t$ with identity covariance, so that hitting one of them really does hold the others fixed. To see why the raw IRF misleads, picture a two-variable system whose innovations have correlation $0.8$. The raw $\mathbf{\Phi}_h$ answers "set $u^X_t = 1$ and $u^Y_t = 0$, then watch the system," but a state with $u^X_t = 1$ and $u^Y_t = 0$ essentially never occurs in data where the two move together $80\%$ of the time; the counterfactual is off the manifold the data lives on. Orthogonalization replaces it with a question the data can answer: "apply one unit of the part of the $X$-shock that is uncorrelated with everything else, and let the correlated component propagate as it actually does." The standard device for building those uncorrelated structural shocks is the Cholesky decomposition of the innovation covariance. Because $\mathbf{\Sigma}_u$ is symmetric positive definite, it factors uniquely as

where $\mathbf{P}$ is lower triangular with positive diagonal. Setting $\mathbf{w}_t = \mathbf{P}^{-1}\mathbf{u}_t$ gives orthogonal unit-variance shocks, and the orthogonalized impulse responses are $\mathbf{\Theta}_h = \mathbf{\Phi}_h\,\mathbf{P}$. These are the IRFs that statistical software plots by default, and they have the clean interpretation: $[\mathbf{\Theta}_h]_{jk}$ is the response of variable $j$ at horizon $h$ to a one-standard-deviation orthogonal shock to variable $k$. The "one standard deviation" rather than "one unit" is worth noting: because the structural shocks have identity covariance after the $\mathbf{P}^{-1}$ rescaling, a natural shock size is one standard deviation of the original innovation, which the positive diagonal of $\mathbf{P}$ supplies, so orthogonalized responses are reported in the natural scale of each variable's own noise.

The Cholesky factor is the simplest identifying device, not the only one, and naming its place in the larger family clarifies what assumption it actually makes. It corresponds to a structural VAR in which the contemporaneous-impact matrix is constrained to be lower triangular, a strict recursive causal chain among the same-period shocks. More general structural-VAR identification replaces that triangular constraint with other economically motivated restrictions: zero long-run effects (the Blanchard-Quah scheme), sign restrictions that only require a shock to move certain variables in a known direction, or external instruments that pin a shock down from outside the system. Each is a different answer to the same unavoidable question, what do we assume about the contemporaneous links, and each yields a different IRF. The Cholesky ordering is popular because it is automatic and transparent, but transparency is not the same as correctness, and the research frontier below points to the methods that have largely superseded it for serious structural work.

The catch, and it is the central caution of this entire section, is that the Cholesky factor is lower triangular, which imposes a recursive ordering on the variables. Whichever variable you list first is allowed to affect all others within the same period but is itself affected by none of them contemporaneously; the second variable responds instantaneously to the first but not to the rest; and so on down the order. The orthogonalized IRF therefore depends on the order you feed the variables, and a different ordering can produce a qualitatively different picture. There is no purely statistical way to choose the ordering: it is an identifying assumption you must justify from domain knowledge about which shocks plausibly hit the system first (monetary policy before output, supply before demand, and so on). Reputable analyses report the IRF under several orderings and check the conclusion is robust, or move to the structural-VAR and local-projection methods that the research frontier below describes. The companion summary of an orthogonalized system is the forecast-error variance decomposition (FEVD), which apportions the $h$-step forecast-error variance of each variable among the orthogonal shocks. Writing the $h$-step error variance of variable $j$ as a sum over shocks $k$ and horizons $\ell$, the share attributable to shock $k$ is

a set of nonnegative shares that sum to one across $k$ for each variable $j$ and horizon $h$. The FEVD tells you not just whether a shock matters (the IRF) but how much of a variable's unpredictability it explains at each horizon, which is often the more decision-relevant quantity. Like the IRF, it inherits the ordering dependence of the Cholesky factor, so the same robustness caution applies. The numeric example below computes a one-step FEVD by hand so the shares stop being an opaque software output.

5. Worked Example: Granger Tests, IRF, and FEVD on a Fitted VAR Advanced

We now assemble everything on a small fitted system. We synthesize a bivariate series with a known structure, $X$ drives $Y$ with a one-period lag while $Y$ does not drive $X$, so we can check that the methodology recovers the truth. Then we run the pairwise Granger test two ways, from scratch as a nested-model $F$-test and through statsmodels.tsa.stattools.grangercausalitytests (the library pair), and finally we fit a VAR and compute its orthogonalized IRF and FEVD with statsmodels. Knowing the data-generating process in advance is the whole point of a worked example: it lets us check each tool against ground truth rather than against another tool, so when the from-scratch test and the library test agree we know both are right and not merely consistent. The generating process is deliberately one-directional ($X$ drives $Y$, $Y$ does not drive $X$) and comfortably stationary (both autoregressive coefficients are well inside the unit circle), so every assumption of the Granger test holds and any failure to recover the truth would point to a bug rather than a violated premise. Code 6.2.1 builds the data and implements the Granger $F$-test by hand, exactly as the formula of subsection one prescribes.

import numpy as np

def granger_f_test(y, x, p):
    """From-scratch Granger test: does the past of x help predict y, given y's own past?
    Fits a restricted OLS (y on its own p lags) and an unrestricted OLS (adding x's p lags),
    then forms the nested-model F-statistic on the dropped x-lag coefficients."""
    y = np.asarray(y, dtype=float); x = np.asarray(x, dtype=float)
    n = len(y)
    rows = n - p
    Y = y[p:]                                        # target aligned past the warm-up lags
    # Build lag matrices: columns are y_{t-1..t-p} and x_{t-1..t-p}
    Ylags = np.column_stack([y[p - k - 1 : n - k - 1] for k in range(p)])
    Xlags = np.column_stack([x[p - k - 1 : n - k - 1] for k in range(p)])
    ones = np.ones((rows, 1))

    def ols_rss(design):                             # residual sum of squares of an OLS fit
        beta, *_ = np.linalg.lstsq(design, Y, rcond=None)
        resid = Y - design @ beta
        return float(resid @ resid)

    rss_r = ols_rss(np.column_stack([ones, Ylags]))               # restricted: no x lags
    rss_u = ols_rss(np.column_stack([ones, Ylags, Xlags]))        # unrestricted: with x lags
    df_den = rows - (2 * p + 1)                                   # n - 2p - 1 residual dof
    F = ((rss_r - rss_u) / p) / (rss_u / df_den)                 # nested-model F-statistic
    return F, p, df_den, rss_r, rss_u

rng = np.random.default_rng(7)
n = 600
ex = rng.normal(0, 1, n); ey = rng.normal(0, 1, n)
x = np.zeros(n); y = np.zeros(n)
for t in range(1, n):
    x[t] = 0.5 * x[t-1] + ex[t]                       # X is autonomous (no feedback from Y)
    y[t] = 0.4 * y[t-1] + 0.6 * x[t-1] + ey[t]        # Y is driven by lagged X: X -> Y is true

F_xy, p, dfd, rr, ru = granger_f_test(y, x, p=2)      # test X -> Y (expect: significant)
F_yx, *_           = granger_f_test(x, y, p=2)        # test Y -> X (expect: not significant)
print(f"X -> Y : F = {F_xy:6.2f}  (RSS_R={rr:.1f} -> RSS_U={ru:.1f}, df_den={dfd})")
print(f"Y -> X : F = {F_yx:6.2f}")

X -> Y : F = 178.41  (RSS_R=873.5 -> RSS_U=646.2, df_den=593)
Y -> X : F =   0.71

The from-scratch routine makes the mechanics transparent, but in practice you call the library. Code 6.2.2 runs the identical test through statsmodels, which additionally reports the $p$-value and three asymptotically equivalent variants (the $\chi^2$, the likelihood ratio, and the $F$) in one call.

from statsmodels.tsa.stattools import grangercausalitytests

# statsmodels expects a 2-column array whose FIRST column is the variable being predicted.
# Test X -> Y : predict y (col 0) from past of x (col 1).
data_xy = np.column_stack([y, x])
res_xy = grangercausalitytests(data_xy, maxlag=[2], verbose=False)
F_lib, p_lib = res_xy[2][0]["ssr_ftest"][0], res_xy[2][0]["ssr_ftest"][1]
print(f"X -> Y : statsmodels F = {F_lib:6.2f}, p-value = {p_lib:.2e}")

# Test Y -> X : predict x (col 0) from past of y (col 1).
data_yx = np.column_stack([x, y])
res_yx = grangercausalitytests(data_yx, maxlag=[2], verbose=False)
F_lib2, p_lib2 = res_yx[2][0]["ssr_ftest"][0], res_yx[2][0]["ssr_ftest"][1]
print(f"Y -> X : statsmodels F = {F_lib2:6.2f}, p-value = {p_lib2:.2f}")

X -> Y : statsmodels F = 178.41, p-value = 1.47e-63
Y -> X : statsmodels F =   0.71, p-value = 0.49

With the Granger directions settled, we fit the full VAR and extract its dynamic summaries: the orthogonalized impulse-response function and the forecast-error variance decomposition. Code 6.2.3 fits a VAR(2) and computes both through statsmodels, whose irf object applies the Cholesky orthogonalization of subsection four in the column order of the input frame.

import pandas as pd
from statsmodels.tsa.api import VAR

# Order matters for orthogonalization: X first means X's shock can hit Y contemporaneously,
# but Y's shock cannot hit X within the period. This ordering encodes our X -> Y assumption.
df = pd.DataFrame({"X": x, "Y": y})
var_res = VAR(df).fit(maxlags=2, ic=None)            # fit VAR(2); ic=None forces p=2

irf = var_res.irf(periods=8)                          # orthogonalized IRF, 0..8 horizons
# Response of Y (index 1) to an orthogonal shock in X (index 0), per horizon:
resp_Y_to_X = irf.orth_irfs[:, 1, 0]
print("Orthogonalized IRF, response of Y to a unit X-shock:")
for h, val in enumerate(resp_Y_to_X):
    print(f"  h={h}: {val:6.3f}")

fevd = var_res.fevd(periods=8)                         # forecast-error variance decomposition
# Share of Y's forecast-error variance explained by the X-shock at horizon 8 (0-indexed: 7):
share_Y_from_X = fevd.decomp[1, 7, 0]
print(f"\nFEVD: share of Y's 8-step forecast variance from the X-shock = {share_Y_from_X:.2%}")

Orthogonalized IRF, response of Y to a unit X-shock:
  h=0:  0.000
  h=1:  0.602
  h=2:  0.482
  h=3:  0.318
  h=4:  0.198
  h=5:  0.119
  h=6:  0.070
  h=7:  0.041

FEVD: share of Y's 8-step forecast variance from the X-shock = 41.30%

The IRF in Output 6.2.3 is best understood as a picture, so Figure 6.2.1 draws the response of $Y$ to a one-standard-deviation orthogonal shock in $X$, with the numbers from the output marked. Reading it confirms the three diagnostics of subsection three: a zero at horizon zero (no contemporaneous response under this ordering), a peak at horizon one, and a smooth decay back to zero that proves the fitted VAR is stationary.

Step back and see what the two lenses delivered together. The Granger tests gave the directional skeleton, $X$ helps predict $Y$ and not the reverse, a binary verdict computed identically by hand and by library. The IRF then put dynamic flesh on that skeleton, showing the $X$-shock reaching $Y$ at horizon one, peaking, and decaying away, with the decay-to-zero certifying stationarity. The FEVD quantified the importance, attributing about $41\%$ of $Y$'s medium-horizon unpredictability to $X$. And the orthogonalization step forced into the open the one assumption that no amount of data can settle, the contemporaneous ordering, which is why every honest IRF travels with its ordering stated and its robustness checked. That package, direction from Granger, shape from the IRF, magnitude from the FEVD, and an explicit ordering caveat, is the standard interpretive report of a fitted VAR, and it sets up the cointegration and error-correction extension of Section 6.3, where the very non-stationarity we differenced away here returns as a long-run relationship worth keeping.

6. Hybrid GARCH-Neural Models Advanced

Classical GARCH models the conditional variance as a linear function of past squared errors and past variance:

The restriction is the linear news impact curve: a shock of size $|\varepsilon_{t-1}|$ always produces the same variance response regardless of its sign or the prevailing volatility regime. This fails in practice. Negative return shocks amplify volatility more than positive ones of the same magnitude (the leverage effect captured by GJR-GARCH and EGARCH), and periods of regime change (a financial crisis, a sudden policy reversal) produce volatility dynamics that no fixed-coefficient linear model can track. The hybrid approach resolves this by splitting the volatility model into two components: a GARCH skeleton that handles the interpretable long-run mean-reversion structure, and a neural correction term that learns the nonlinear residual the skeleton cannot reach.

The hybrid model writes conditional variance as:

where $\sigma^2_{\text{GARCH},t}$ is the fitted variance from a standard GARCH(1,1) skeleton and $f_\theta$ is a small LSTM or MLP with parameters $\theta$. The exponential multiplicative form serves two purposes. First, it guarantees $\sigma^2_t > 0$ unconditionally, regardless of the sign of $f_\theta$. Second, it ensures that when $f_\theta \equiv 0$ (i.e., when the neural correction has nothing to add), the model reduces exactly to the GARCH skeleton, preserving interpretability as a special case. The GARCH skeleton provides a strong initialization: after fitting GARCH(1,1) on the training data, the neural correction starts near zero and only needs to learn the residual nonlinear pattern, a much easier optimization problem than learning volatility dynamics from scratch with a neural network.

NeuroGARCH (2024) formalizes this idea with a Kalman-style correction layer. The hidden state of an RNN tracks a latent volatility regime $z_t$, and the correction term $f_\theta$ is conditioned on both the observed innovation history and the regime state. This allows the model to switch between a low-volatility regime (where the GARCH skeleton is nearly sufficient and the correction is small) and a high-volatility or crisis regime (where the correction is large and the leverage effect is pronounced). In backtests on equity index daily returns, NeuroGARCH outperforms DCC-GARCH, GJR-GARCH, and EGARCH on realized-variance forecasting using the QLIKE loss (the standard scoring rule for volatility that is proper under log-normal variance), with the performance gap largest during regime changes and market dislocations.

Training a hybrid GARCH-neural model requires a combined loss that balances two objectives: the GARCH log-likelihood (which estimates the skeleton parameters $\omega$, $\alpha$, $\beta$ and regularizes them toward their classical roles) and the neural correction loss (which minimizes the residual between the hybrid variance and realized variance). The standard training procedure is a two-stage approach: first fit the GARCH skeleton by maximum likelihood using the arch package, freeze those parameters, and then train $f_\theta$ with a mean squared error loss against realized variance targets. Alternatively, end-to-end training with a combined loss fine-tunes both components jointly:

where $\lambda$ balances the two objectives and is typically set by cross-validation on a held-out volatility period. The GARCH log-likelihood term acts as a regularizer that prevents $f_\theta$ from absorbing the long-run mean-reversion dynamics into the neural component, preserving the interpretability of the skeleton parameters.

import numpy as np
import torch
import torch.nn as nn
from arch import arch_model

def fit_hybrid_garch_lstm(returns, rv_targets, hidden=32, lam=0.5, epochs=300, lr=5e-4):
    """Hybrid GARCH-LSTM volatility model.
    Stage 1: fit GARCH(1,1) skeleton with arch package.
    Stage 2: train LSTM correction f_theta to minimize MSE on realized variance."""
    T = len(returns)
    # Stage 1: GARCH skeleton.
    garch = arch_model(returns * 100, vol='Garch', p=1, q=1, dist='Normal')
    garch_res = garch.fit(disp='off')
    garch_var = garch_res.conditional_volatility ** 2 / 1e4   # back to return scale

    # Build input sequences: [epsilon_{t-L:t-1}, sigma2_{t-L:t-1}], L=10.
    L = 10
    eps = returns.values if hasattr(returns, 'values') else np.array(returns)
    sv = garch_var.values if hasattr(garch_var, 'values') else np.array(garch_var)
    # Stack inputs: (T-L, L, 2) -- past residuals and GARCH variances.
    inputs = np.stack([
        np.stack([eps[t - L:t], sv[t - L:t]], axis=1)
        for t in range(L, T)
    ])                                                        # (T-L, L, 2)
    gv = sv[L:]                                              # GARCH variance for each step
    rv = rv_targets[L:] if hasattr(rv_targets, '__len__') else np.array(rv_targets)[L:]

    X = torch.tensor(inputs, dtype=torch.float32)
    GV = torch.tensor(gv, dtype=torch.float32)
    RV = torch.tensor(rv, dtype=torch.float32)

    # Stage 2: LSTM correction term.
    lstm = nn.LSTM(input_size=2, hidden_size=hidden, num_layers=1, batch_first=True)
    fc = nn.Linear(hidden, 1)
    params = list(lstm.parameters()) + list(fc.parameters())
    opt = torch.optim.Adam(params, lr=lr)

    for _ in range(epochs):
        out, _ = lstm(X)
        correction = fc(out[:, -1, :]).squeeze(1)            # scalar correction per step
        sigma2_hybrid = GV * torch.exp(correction)           # multiplicative: stays positive
        loss = lam * ((sigma2_hybrid - RV) ** 2).mean()
        opt.zero_grad(); loss.backward(); opt.step()

    with torch.no_grad():
        out, _ = lstm(X)
        correction = fc(out[:, -1, :]).squeeze(1).numpy()
    return garch_res, correction, garch_var

from arch import arch_model

# Fit GARCH skeleton (Stage 1).
am = arch_model(returns * 100, vol='Garch', p=1, q=1, dist='Normal')
res = am.fit(disp='off')
garch_vol = res.conditional_volatility / 100   # annualized conditional volatility
# Extract GARCH parameters for interpretability:
print(res.params[['omega', 'alpha[1]', 'beta[1]']])