VAR and vector error correction models (VECM) in Python

Output, inflation and interest rates respond to one another with lags, and it is rarely clear in advance which of them is exogenous. A vector autoregression (VAR) sidesteps that choice by regressing every variable on lags of every variable in the system. The fitted VAR shows whether one variable helps predict another, how a shock spreads through the system over time and which shocks drive each variable's forecast errors. When the variables wander without a fixed mean but are tied together in the long run, the vector error-correction model (VECM) keeps that long-run link instead of differencing it away.

What VAR and VECM models estimate

A VAR(p) for K variables stacks K regressions with the same right-hand side:

y_t = v + A₁ y_{t−1} + … + A_p y_{t−p} + u_t,        E[u_t u_t'] = Σ

Each A_j is a K × K matrix of lag coefficients, and the errors may be correlated across equations. Because every equation has the same regressors, ordinary least squares equation by equation is the Gaussian maximum-likelihood estimator. A VAR is stable when every eigenvalue of its companion matrix lies inside the unit circle; shocks then die out.

The fitted VAR is read through three tools:

• Granger causality. Variable x Granger-causes y if lags of x improve the prediction of y given the lags of y itself (Granger 1969). The test is a Wald test that the lag coefficients of x in the y equation are jointly zero. It is a statement about prediction, not about structural causation.
• Impulse responses. These trace the effect of a one-time shock on every variable over the following periods, the way Sims (1980) proposed to read a VAR. Because the errors are correlated, orthogonalized responses separate the shocks with a Cholesky factor of Σ, so they depend on the order of the variables: the first variable's shock moves everything on impact. Generalized responses (Pesaran and Shin 1998) do not depend on the order.
• Variance decomposition. This reports the share of each variable's h-step forecast-error variance that is due to each orthogonalized shock.

A VAR assumes the variables are stationary. If they are integrated of order one, I(1), you can difference them, but if some linear combination β'y_t is stationary the variables are cointegrated, and differencing throws away the long-run relation. The VECM writes the VAR in levels as

Δy_t = α β' y_{t−1} + Γ₁ Δy_{t−1} + … + Γ_{p−1} Δy_{t−p+1} + deterministic terms + e_t

The r columns of β are the cointegrating vectors and α holds the adjustment speeds: how strongly each variable responds when the system is away from its long-run relation. Johansen (1991) derives the maximum-likelihood estimator, a reduced-rank regression, and the tests of the number of cointegrating vectors. The cointegrating rank r is chosen first with Johansen's trace test; the cointegration page covers that step.

In OpenEconometrics

oe.var(data=df, y=["growth", "inflation"], time="quarter", lags=1) fits a VAR. With the data declared by tsset quarter, it corresponds to the Stata command var growth inflation, lags(1). Note the difference in convention: lags=p includes all lags 1 to p, which Stata writes lags(1/p).

• Order and time. The order of y is the Cholesky order of orthogonalized responses and the order of Johansen's normalization, so y must be a list. time= is optional: rows are sorted by it, integer periods must be consecutive and distinct, and datetime values are taken as consecutive periods with a recorded warning. Without it, the row order is the time order.
• Options. x= adds exogenous regressors (Stata's exog()), trend=True a linear trend and constant=False drops the constant. small=True switches to t and F statistics with T − m degrees of freedom and dfk=True estimates Σ with the divisor T − m, as Stata's small and dfk options do; covariance="robust" gives a heteroskedasticity-robust covariance.
• Lag selection. oe.varsoc(data=df, y=[...], maxlag=4) fits lags 0 to maxlag on a common sample and reports LR, FPE, AIC, HQIC and SBIC, like Stata's varsoc; attrs["selected"] holds each criterion's choice.
• Stored diagnostics. The result carries Granger causality tests (fit.extra["granger"]), lag-exclusion, normality and residual-autocorrelation LM tests, and the stability eigenvalues (fit.extra["stability"]).
• Responses and forecasts. oe.irf(fit, steps=8, kind="orthogonalized") returns simple, orthogonalized or generalized responses with delta-method standard errors, plus the variance decomposition for the orthogonalized kind. oe.forecast(fit, steps) gives dynamic forecasts with standard errors.

oe.vec(data=df, y=["long_rate", "short_rate"], lags=2, rank=1, trend="rconstant") fits a VECM by Johansen's method and corresponds to vec long_rate short_rate, lags(2) rank(1) trend(rconstant). lags counts lags of the VAR in levels, so the VECM has lags - 1 lagged differences. trend= takes Stata's five cases: "none", "rconstant" (constant only inside the cointegrating equations), "constant" (the default), "rtrend" and "trend". The cointegrating vectors are normalized so that the first r variables have coefficient one; order y accordingly. The result holds the cointegrating equations in extra["beta"], the adjustment coefficients as the terms D_<variable>:L._ce1 and Johansen's rank estimate in extra["selected_rank"]. The VAR and cointegration family guide gives the formulas and conventions.

Supported options

Options of oe.var

Options of oe.vec

oe.varsoc

Stata command: varsoc.

oe.irf

Stata command: irf create.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script runs two simulations. The first is 300 quarters of a stable VAR(1) with constants 1.0 and 0.5: growth depends on its own lag (0.5) but not on lagged inflation, and inflation on lagged growth (0.2) and its own lag (0.6). The shocks have unit variance and correlation 0.3, so with growth ordered first, the true orthogonalized response of inflation to a growth shock is 0.30, 0.38, 0.33 and 0.25 at steps 0 to 3. The second is 300 months of two I(1) interest rates tied by the long-run relation long_rate − short_rate − 1 = 0: each month the long rate closes 0.1 and the short rate 0.2 of the gap, so β = (1, −1) with constant −1 and α = (−0.1, 0.2).

import numpy as np
import pandas as pd
import openecon as oe


def simulate_var(n=300, burn=100, seed=21):
    """Stationary VAR(1): lagged growth moves inflation, lagged inflation does not move growth."""
    rng = np.random.default_rng(seed)
    c = np.array([1.0, 0.5])
    A = np.array([[0.5, 0.0],    # growth_t    = 1.0 + 0.5 growth_{t-1}                  + u1
                  [0.2, 0.6]])   # inflation_t = 0.5 + 0.2 growth_{t-1} + 0.6 infl_{t-1} + u2
    chol = np.linalg.cholesky(np.array([[1.0, 0.3], [0.3, 1.0]]))
    y = np.zeros((n + burn, 2))
    for t in range(1, n + burn):
        y[t] = c + A @ y[t - 1] + chol @ rng.normal(size=2)
    df = pd.DataFrame(y[burn:], columns=["growth", "inflation"])
    df["quarter"] = np.arange(1, n + 1)
    return df


def simulate_rates(n=300, seed=22):
    """Two I(1) interest rates whose spread minus 1 is stationary: one cointegrating vector."""
    rng = np.random.default_rng(seed)
    alpha = np.array([-0.1, 0.2])            # speed at which each rate closes the gap
    y = np.zeros((n, 2))
    y[0] = [3.0, 2.0]
    for t in range(1, n):
        gap = y[t - 1, 0] - y[t - 1, 1] - 1.0
        y[t] = y[t - 1] + alpha * gap + rng.normal(scale=0.2, size=2)
    return pd.DataFrame({"month": np.arange(1, n + 1), "long_rate": y[:, 0], "short_rate": y[:, 1]})


if __name__ == "__main__":
    macro = simulate_var()
    lag = oe.varsoc(data=macro, y=["growth", "inflation"], time="quarter", maxlag=4).attrs["selected"]["aic"]
    print("Lag order chosen by AIC:", lag)
    fit = oe.var(data=macro, y=["growth", "inflation"], time="quarter", lags=lag)
    print(fit.summary())
    irf = oe.irf(fit, steps=3, kind="orthogonalized")
    print(irf.query("impulse == 'growth' and response == 'inflation'")
             [["step", "irf", "std_error", "ci_low", "ci_high"]].round(3).to_string(index=False))

    rates = simulate_rates()
    vec = oe.vec(data=rates, y=["long_rate", "short_rate"], time="month", lags=2, rank=1, trend="rconstant")
    print("Johansen rank estimate:", vec.extra["selected_rank"])
    for c in vec.extra["beta"][0]["coefficients"]:
        se = "normalized" if c["std_error"] is None else f"{c['std_error']:.4f}"
        print(f"beta  {c['variable']:<10} {c['estimate']:8.4f}  {se}")
    for c in vec.coefficients:
        if c.term.endswith("L._ce1"):
            print(f"alpha {c.term:<20} {c.estimate:8.4f}  {c.std_error:.4f}")

With openecon 0.3.20a1 from PyPI the script prints:

Lag order chosen by AIC: 1
Vector autoregression — growth
Observations: 299  |  Covariance: nonrobust  |  Confidence: 95%

Term                      Estimate  Std. error         z     P>|stat|   CI lower   CI upper
----------------------  ----------  ----------  --------  -----------  ---------  ---------
[growth]
growth:L1.growth          0.498526    0.056088   8.88828  6.20606e-19   0.388595   0.608456
growth:L1.inflation     -0.0295344   0.0460748  -0.64101     0.521516  -0.119839  0.0607705
growth:Intercept           1.16531    0.139955   8.32632   8.3391e-17   0.891001    1.43961
[inflation]
inflation:L1.growth       0.204078   0.0543721   3.75336  0.000174483  0.0975104   0.310645
inflation:L1.inflation    0.629053   0.0446652   14.0837   4.7816e-45   0.541511   0.716595
inflation:Intercept       0.494182    0.135673   3.64245  0.000270054   0.228268   0.760096

Excluded observations: 1

log_likelihood: -821.15  |  aic: 1654.3  |  bic: 1676.5  |  hqic: 1663.19
aic_per_obs: 5.53277  |  hqic_per_obs: 5.56249  |  sbic_per_obs: 5.60703  |  fpe: 0.866772
det_sigma_ml: 0.832673  |  n_equations: 2  |  n_lags: 1  |  df_eq: 3
df_model: 6  |  df_resid: 296  |  T: 299
Granger causality: all other variables excluded from growth: chi2(1) = 0.410894, p = 0.5215
Granger causality: all other variables excluded from inflation: chi2(1) = 14.0877, p = 0.0001745
Wald test that all endogenous variables at lag 1 are jointly zero: chi2(4) = 394.914, p = 3.493e-84
Jarque-Bera test that the disturbances are jointly normal: chi2(4) = 1.58762, p = 0.811
LM test of no residual autocorrelation at lag 1: chi2(4) = 3.10209, p = 0.5409
LM test of no residual autocorrelation at lag 2: chi2(4) = 1.67604, p = 0.7951
Warning: The first 1 observation(s) supply the lagged values; the estimation sample has 299 observations.
 step   irf  std_error  ci_low  ci_high
    0 0.384      0.054   0.277    0.490
    1 0.448      0.063   0.325    0.570
    2 0.382      0.066   0.254    0.511
    3 0.288      0.059   0.172    0.404
Johansen rank estimate: 1
beta  long_rate    1.0000  normalized
beta  short_rate  -0.9348  0.0779
beta  _cons       -1.0831  0.1974
alpha D_long_rate:L._ce1    -0.0913  0.0329
alpha D_short_rate:L._ce1    0.1909  0.0352

Reading the output

• The lag order is found. AIC picks the true single lag. The own lags are 0.499 and 0.629 (true 0.5 and 0.6), the effect of growth on inflation 0.204 (true 0.2) and that of inflation on growth −0.030 with a standard error of 0.046 (true 0).
• Granger causality runs one way. Lagged inflation does not help predict growth (χ²(1) = 0.41, p = 0.52), while lagged growth helps predict inflation (χ²(1) = 14.09, p = 0.00017). With one lag, each test is the square of the coefficient's z statistic: (−0.641)² = 0.41.
• The residuals pass the checks. The joint Jarque–Bera test (p = 0.81) and the LM tests at lags 1 and 2 (p = 0.54 and 0.80) find no problem.
• The impulse responses track the truth. A one-standard-deviation growth shock raises inflation by 0.384 on impact, 0.448 after one quarter, then 0.382 and 0.288. The true values 0.30, 0.38, 0.33 and 0.25 all lie inside the intervals. The impact effect is non-zero only because the shocks are correlated and growth is ordered first.
• The VECM recovers the long-run relation. Johansen's trace test picks rank 1. The cointegrating vector is normalized on long_rate; the estimate for short_rate is −0.935 with a standard error of 0.078, about 0.8 standard errors from the true −1, and the restricted constant is −1.083 (0.197) against −1.
• Both rates adjust. When the spread is above its long-run value, the long rate falls (α = −0.091, standard error 0.033) and the short rate rises (0.191, 0.035), close to the true −0.1 and 0.2.

Limitations

• oe.var accepts no weights.
• The common oe.predict and oe.margins adapters do not cover oe.var; use its own result tables and family helpers.
• oe.vec accepts no weights.
• The common oe.predict and oe.margins adapters do not cover oe.vec; use its own result tables and family helpers.
• lags is a single integer: lags 1 to p all enter (Stata's lags(1 3) with holes is not available), and exogenous regressors enter only contemporaneously.
• There are no coefficient constraints, no structural VAR (svar), no bootstrap standard errors for responses or forecasts and no cumulative responses.
• oe.vec uses only Johansen's normalization and the conventional covariance, has no seasonal indicators and no constraints on α or β. Its forecast standard errors ignore parameter uncertainty and its impulse responses have no standard errors.
• VAR forecast standard errors include parameter uncertainty only when the model has no exogenous regressors and no trend.
• Each series must be gap free; panels of series are not supported, and a system may have at most 1,000 coefficients and 50 variables.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3), 424–438. doi:10.2307/1912791
• Johansen, S. (1991). Estimation and hypothesis testing of cointegration vectors in Gaussian vector autoregressive models. Econometrica, 59(6), 1551–1580. doi:10.2307/2938278
• Pesaran, H. H. and Shin, Y. (1998). Generalized impulse response analysis in linear multivariate models. Economics Letters, 58(1), 17–29. doi:10.1016/S0165-1765(97)00214-0
• Sims, C. A. (1980). Macroeconomics and reality. Econometrica, 48(1), 1–48. doi:10.2307/1912017