irfCompute#

Purpose#

Compute orthogonalized impulse response functions using Cholesky identification.

Format#

irf = irfCompute(result, n_ahead)#
Parameters:
  • result (struct) – an instance of a varResult or bvarResult structure.

  • n_ahead (scalar) – number of horizons to compute (e.g., 20).

  • var_names (Mx1 string array) – Optional keyword, override variable names from the estimation result.

  • quiet (scalar) – Optional keyword, set to 1 to suppress printed output. Default = 0.

Returns:

irf (struct) –

An instance of an irfResult structure containing:

irf.irf

Array of (n_ahead+1) mxm matrices. irf.irf[h+1] is the mxm impulse response matrix at horizon h. Element [i, j] is the response of variable i to a one-standard-deviation shock to variable j. Index 1 is the impact response (h=0).

irf.n_ahead

Scalar, number of horizons computed.

irf.m

Scalar, number of variables.

irf.var_names

Mx1 string array, variable names.

irf.ident

String, identification method: "cholesky" or "generalized".

Model#

An impulse response function (IRF) traces the dynamic effect of a one-standard-deviation structural shock to variable \(j\) on variable \(i\) over \(h\) periods.

For a VAR(p) in companion form \(Y_t = F Y_{t-1} + G \varepsilon_t\), the reduced-form IRF at horizon \(h\) is:

\[\Phi_h = J \, F^h \, J'\]

where \(F\) is the \(mp \times mp\) companion matrix and \(J = [I_m \; 0 \; \cdots \; 0]\) selects the first \(m\) rows.

Cholesky identification: To give shocks a structural interpretation, the reduced-form innovations are orthogonalized via the Cholesky factorization \(\Sigma = P P'\) where \(P\) is lower triangular. The structural IRF is:

\[\Theta_h = \Phi_h \, P\]

Element \(\Theta_h[i, j]\) is the response of variable \(i\) at horizon \(h\) to a one-standard-deviation shock to variable \(j\).

Identification assumption: Cholesky identification imposes a recursive causal ordering. Variable 1 can affect all others contemporaneously; variable \(m\) is affected by all others but affects none contemporaneously. This assumption is appropriate when there is a natural fast-to-slow ordering (e.g., financial variables respond faster than real activity). Algorithm ———

  1. Extract companion matrix \(F\) and Cholesky factor \(P = \text{chol}(\Sigma)'\) from the VAR estimates.

  2. Iterate: For \(h = 0, 1, \ldots, n\_ahead\):

    \[\Theta_h = J \, F^h \, J' \, P\]

    The companion power \(F^h\) is computed iteratively (matrix multiplication, not matrix exponentiation) for numerical stability.

  3. Store \(\Theta_0, \Theta_1, \ldots, \Theta_{n\_ahead}\) as an array of \(m \times m\) matrices.

Complexity: \(O(n\_ahead \cdot m^2 p^2)\) — dominated by the \(mp \times mp\) matrix multiplications. Sub-millisecond for typical systems. Examples ——–

Monetary Policy Shock#

Trace the effect of a federal funds rate shock on GDP and CPI:

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

// Variable ordering: GDP (slow), CPI (medium), FFR (fast policy instrument)
// This ordering means: FFR shocks can affect GDP and CPI contemporaneously,
// but GDP shocks take one period to reach FFR.
result = varFit(data, 4);

struct irfResult irf;
irf = irfCompute(result, 20);

Output:

================================================================================
Impulse Response Functions (cholesky)
Horizons: 0-20
================================================================================

Shock to: GDP
  h          GDP       CPI       FFR
--------------------------------------------------------------------------------
  0     0.5280     0.0456     0.0919
  1     0.1859     0.0810     0.2753
  2     0.1600     0.0612     0.4442
 ...
================================================================================

The impact response (h=0) shows that a 1-SD GDP shock raises GDP by 0.528, CPI by 0.046, and FFR by 0.092 — consistent with the central bank responding to output movements within the quarter.

IRF from BVAR with Shrinkage#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

ctl = bvarControlCreate();
ctl.p = 4;

struct bvarResult br;
br = bvarFit(data, ctl, quiet=1);

// IRF at the posterior mean of B and Sigma
irf = irfCompute(br, 20);

For posterior IRF bands (credible intervals), use irfSvCompute() with an SV-BVAR result.

Plotting IRFs#

Reshape IRF results into a plot-ready dataframe:

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));
result = varFit(data, 4, quiet=1);
irf = irfCompute(result, 20, quiet=1);

// Get plot data: (n_ahead+1) x (m*m) matrix with column names
struct DataFrame plot_data;
plot_data = irfPlotData(irf);

// Plot GDP response to FFR shock
plotXY(seqa(0, 1, 21), plot_data[., "GDP<-FFR"]);

Troubleshooting#

IRFs don’t decay to zero: If the VAR is near-nonstationary (max eigenvalue close to 1), IRFs can be very persistent. This is not a numerical error — it reflects the model’s dynamics. Check result.max_eigenvalue. For near-unit-root systems, consider:

  • Using longer horizons (40-60 periods instead of 20).

  • Differencing the data.

  • Adding sum-of-coefficients priors (bvarFit() with lambda6 > 0).

IRFs are sensitive to variable ordering: This is inherent to Cholesky identification — different orderings produce different structural shocks. If the ordering is uncertain, use girfCompute() (generalized IRF, ordering-invariant) or svarIdentify() (sign restrictions).

Impact response has wrong sign: Check the variable ordering. In Cholesky identification, the first variable’s shock is unrestricted; later variables’ shocks are residualized. A monetary policy variable (FFR) should typically be ordered last so its shock is “purged” of contemporaneous output and price movements. Remarks ——-

Identification: The Cholesky decomposition of \(\Sigma\) is used to orthogonalize the innovations. The ordering of variables in the data determines the recursive causal structure: variable 1 can affect all others contemporaneously, variable 2 can affect variables 3, …, m but not 1, and so on.

To change the identification ordering, reorder the columns of the data before calling varFit() or bvarFit().

For ordering-invariant responses, use girfCompute() (generalized IRF, Pesaran & Shin 1998). For theory-based identification with sign/zero restrictions, see svarIdentify().

For BVAR, the IRF is computed at the posterior mean of B and \(\Sigma\). For posterior IRF bands, use irfSvCompute() with an bvarSvResult.

Indexing convention: irf.irf[1, ., .] is the impact response (h=0). irf.irf[h+1, ., .] is the response at horizon h. Element irf.irf[h+1, i, j] is the response of variable i to a shock to variable j. Verification ————

Verified against R vars::irf() with boot=FALSE at \(10^{-6}\) tolerance on a 2-variable VAR(1) with known DGP. Tests cover impact values, decay patterns at h=1 and h=2, and the Cholesky lower-triangularity constraint (zero upper-off-diagonal at h=0). See gausslib-var/tests/r_benchmark.rs.

Additionally verified against ECB BEAR Cholesky IRFs on matched-prior BVAR(4), covering all 9 shock-response pairs at horizons 0, 10, and 20 (17 tests). See crossval/bear_matched_irf.e. References ———-

  • Lutkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer. Chapter 2.3 (IRF computation), Chapter 9 (structural identification).

  • Pesaran, M.H. and Y. Shin (1998). “Generalized impulse response analysis in linear multivariate models.” Economics Letters, 58(1), 17-29.

  • Sims, C.A. (1980). “Macroeconomics and reality.” Econometrica, 48(1), 1-48.

Library#

timeseries

Source#

irf.src