Time Series Econometrics

Multivariate Time Series Econometrics

TL;DR — Fit VAR models, test Granger causality, test for cointegration with Johansen, estimate VECMs, plot impulse response functions with bootstrap confidence intervals, and decompose forecast error variance.

When to Use

This Skill is the starting point for any empirical work involving multiple time series in economics and finance — GDP, inflation, interest rates, exchange rates, commodity prices, and asset prices.

Background

VAR(p) Model

A VAR(p) model for k-dimensional vector Y_t is:

Y_t = c + A_1 Y_{t-1} + A_2 Y_{t-2} + ... + A_p Y_{t-p} + u_t

where u_t ~ N(0, Σ) is the innovation vector. Lag p is selected by information criteria: AIC, BIC (Schwarz), or HQIC.

VAR stability condition: All eigenvalues of the companion matrix must lie strictly inside the unit circle. If any eigenvalue >= 1, the system is non-stationary and you should consider VECM.

Granger Causality

Variable X Granger-causes Y if knowing past X improves forecasts of Y beyond past Y alone. Tested via an F-test on the restriction that all lagged X coefficients in the Y equation are jointly zero.

This is a predictive concept, not a structural causal one. Granger causality ≠ structural causality.

Johansen Cointegration

Two or more I(1) series are cointegrated if a linear combination is I(0). Johansen (1988) provides two likelihood-ratio tests:

• Trace test: H0 = at most r cointegrating vectors
• Maximum eigenvalue test: H0 = r vs H1 = r+1 cointegrating vectors

Critical values differ by deterministic specification (no constant, constant in CE, trend in CE). Use statsmodels.tsa.vector_ar.vecm.coint_johansen.

VECM

A VECM re-parameterizes the VAR for I(1) cointegrated variables:

ΔY_t = c + Π Y_{t-1} + Γ_1 ΔY_{t-1} + ... + Γ_{p-1} ΔY_{t-p+1} + u_t

where Π = αβ' is the error correction matrix, β is the cointegrating vector, and α is the speed-of-adjustment vector.

Impulse Response Functions

The IRF measures the dynamic response of Y_{t+h} to a one-unit shock in u_{j,t}. Cholesky decomposition of Σ identifies structural shocks via a lower-triangular ordering (Cholesky ordering matters for interpretation).

FEVD

FEVD at horizon h decomposes the forecast variance of variable i into contributions from each structural shock j. FEVD_{ij}(h) ∈ [0, 1] and Σ_j FEVD_{ij}(h) = 1.

Environment Setup

# Create environment
conda create -n tsecono python=3.11 -y
conda activate tsecono
pip install statsmodels>=0.14 arch>=6.0 numpy>=1.23 pandas>=1.5 matplotlib>=3.6

# Verify
python -c "import statsmodels; print('statsmodels', statsmodels.__version__)"
python -c "import arch; print('arch', arch.__version__)"

Core Workflow

Step 1 — Stationarity Tests and Data Preparation

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.vector_ar.var_model import VAR
from statsmodels.tsa.vector_ar.vecm import coint_johansen, VECM
import warnings

np.random.seed(42)


def adf_summary(series: pd.Series, name: str = "", max_diff: int = 2) -> pd.DataFrame:
    """
    Run ADF test on a series and its differences until stationarity is achieved.

    Args:
        series:   Time series to test.
        name:     Label for display.
        max_diff: Maximum order of differencing to attempt.

    Returns:
        DataFrame with columns: variable, adf_stat, p_value, lags, conclusion.
    """
    records = []
    for d in range(max_diff + 1):
        s = series.diff(d).dropna() if d > 0 else series.dropna()
        label = f"Δ^{d} {name}" if d > 0 else name
        adf_result = adfuller(s, autolag="AIC")
        stat, pval, lags_used = adf_result[0], adf_result[1], adf_result[2]
        conclusion = "stationary (I(0))" if pval < 0.05 else "unit root (non-stationary)"
        records.append({
            "variable": label,
            "adf_stat": round(stat, 4),
            "p_value": round(pval, 4),
            "lags": lags_used,
            "conclusion": conclusion,
        })
        if pval < 0.05:
            break
    return pd.DataFrame(records)


def generate_var_data(
    n: int = 300,
    k: int = 3,
    p: int = 2,
    cointegrated: bool = False,
) -> pd.DataFrame:
    """
    Simulate a stationary VAR(p) or cointegrated I(1) system.

    Args:
        n:             Number of observations.
        k:             Number of variables.
        p:             Lag order of the true VAR.
        cointegrated:  If True, generate cointegrated I(1) system.

    Returns:
        DataFrame with k columns and n rows.
    """
    if not cointegrated:
        # Stationary VAR: companion matrix has spectral radius < 1
        A1 = np.array([[0.5, 0.1, 0.0],
                       [0.0, 0.4, 0.2],
                       [0.1, 0.0, 0.3]])
        A2 = np.array([[0.1, 0.0, 0.0],
                       [0.0, 0.1, 0.0],
                       [0.0, 0.0, 0.1]])
        Sigma = np.array([[1.0, 0.3, 0.1],
                          [0.3, 1.0, 0.2],
                          [0.1, 0.2, 1.0]])
        L = np.linalg.cholesky(Sigma)

        Y = np.zeros((n + 50, k))
        for t in range(2, n + 50):
            eps = L @ np.random.randn(k)
            Y[t] = A1 @ Y[t - 1] + A2 @ Y[t - 2] + eps
        Y = Y[50:]
        cols = [f"y{i+1}" for i in range(k)]
        return pd.DataFrame(Y, columns=cols)
    else:
        # I(1) cointegrated: three series share one common stochastic trend
        common_trend = np.cumsum(np.random.randn(n))
        Y = np.column_stack([
            common_trend + 0.5 * np.random.randn(n),
            2 * common_trend + np.random.randn(n),
            -0.5 * common_trend + 1.5 * np.random.randn(n),
        ])
        cols = ["x1", "x2", "x3"]
        return pd.DataFrame(Y, columns=cols)

Step 2 — VAR Estimation and Granger Causality

def fit_var_and_granger(
    df: pd.DataFrame,
    maxlags: int = 8,
    ic: str = "aic",
    verbose: bool = True,
) -> dict:
    """
    Fit a VAR(p) model, run stability check, and test pairwise Granger causality.

    Args:
        df:       DataFrame of stationary time series (k columns).
        maxlags:  Maximum lag order to consider.
        ic:       Information criterion for lag selection: 'aic', 'bic', or 'hqic'.
        verbose:  Print results to console.

    Returns:
        Dictionary with keys: model, results, lag_order, granger_table, stable.
    """
    model = VAR(df)

    # Lag selection
    lag_selection = model.select_order(maxlags=maxlags)
    if verbose:
        print(lag_selection.summary())

    # Fit at IC-optimal lag
    p_opt = getattr(lag_selection, ic)
    p_opt = max(1, p_opt)  # ensure at least 1 lag
    results = model.fit(p_opt)

    if verbose:
        print(results.summary())

    # Stability check
    eigenvalues = np.abs(results.roots)
    stable = bool(np.all(eigenvalues > 1))  # roots of char polynomial outside unit circle
    if verbose:
        print(f"\nVAR stability: {'STABLE' if stable else 'UNSTABLE'}")
        print(f"  Min |root| = {eigenvalues.min():.4f} (must be > 1)")

    # Granger causality: all pairwise tests
    variables = df.columns.tolist()
    granger_rows = []
    for caused in variables:
        for causing in variables:
            if caused == causing:
                continue
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                gc_test = results.test_causality(caused, causing, kind="f")
            granger_rows.append({
                "H0": f"{causing} does NOT Granger-cause {caused}",
                "F_stat": round(gc_test.test_statistic, 4),
                "p_value": round(gc_test.pvalue, 4),
                "reject_H0": gc_test.pvalue < 0.05,
            })

    granger_table = pd.DataFrame(granger_rows)
    if verbose:
        print("\nGranger Causality Tests:")
        print(granger_table.to_string(index=False))

    return {
        "model": model,
        "results": results,
        "lag_order": p_opt,
        "granger_table": granger_table,
        "stable": stable,
    }


def plot_irf_fevd(
    var_results,
    periods: int = 12,
    signif: float = 0.05,
    output_prefix: str = "var",
) -> None:
    """
    Plot impulse response functions and FEVD for a fitted VAR.

    Args:
        var_results: Fitted VAR results object (statsmodels VARResults).
        periods:     IRF horizon in periods.
        signif:      Significance level for bootstrap CI (default 0.05 = 95%).
        output_prefix: Prefix for saved figure filenames.
    """
    irf = var_results.irf(periods)
    k = var_results.neqs
    variables = var_results.names

    # IRF plot
    fig_irf = irf.plot(
        orth=True,  # Cholesky orthogonalization
        signif=signif,
        figsize=(12, 9),
    )
    fig_irf.suptitle(f"Orthogonalized IRFs (Cholesky ordering: {', '.join(variables)})")
    fig_irf.tight_layout()
    fig_irf.savefig(f"{output_prefix}_irf.png", dpi=150)
    print(f"Saved IRF plot to {output_prefix}_irf.png")

    # FEVD at horizons 1, 4, 8, 12
    fevd = var_results.fevd(periods)
    fig_fevd = fevd.plot(figsize=(12, 6))
    fig_fevd.suptitle("Forecast Error Variance Decomposition")
    fig_fevd.tight_layout()
    fig_fevd.savefig(f"{output_prefix}_fevd.png", dpi=150)
    print(f"Saved FEVD plot to {output_prefix}_fevd.png")

    # Print FEVD table at selected horizons
    print("\nFEVD at horizons 1 / 4 / 8 / 12:")
    for h in [1, 4, 8, 12]:
        if h <= periods:
            print(f"\n  Horizon {h}:")
            fevd_df = pd.DataFrame(fevd.decomp[h - 1],
                                   index=variables, columns=variables)
            print(fevd_df.round(3).to_string())

Step 3 — Johansen Test and VECM Estimation

def johansen_vecm_workflow(
    df_levels: pd.DataFrame,
    det_order: int = 0,
    k_ar_diff: int = 1,
    verbose: bool = True,
) -> dict:
    """
    Run Johansen cointegration test and estimate VECM.

    Args:
        df_levels:  DataFrame of I(1) series in levels (not differenced).
        det_order:  Deterministic terms: -1 (none), 0 (constant), 1 (linear trend).
        k_ar_diff:  Lag order in the VECM (lags of differences).
        verbose:    Print results.

    Returns:
        Dictionary with keys: johansen_result, r_selected, vecm_results, beta, alpha.
    """
    # Step 1: Johansen test
    johansen = coint_johansen(df_levels.values, det_order, k_ar_diff)

    variables = df_levels.columns.tolist()
    k = len(variables)

    if verbose:
        print("=" * 60)
        print("Johansen Cointegration Test")
        print("=" * 60)
        trace_cv = johansen.cvt   # (k, 3): 90%, 95%, 99% critical values
        maxeig_cv = johansen.cvm

        print("\nTrace Test:")
        print(f"{'H0: r<=':>12}  {'Trace Stat':>12}  {'CV 5%':>10}  {'Reject?':>8}")
        for r in range(k):
            stat = johansen.lr1[r]
            cv5 = trace_cv[r, 1]
            reject = stat > cv5
            print(f"  r <= {r:>4}  {stat:>12.4f}  {cv5:>10.4f}  {'YES' if reject else 'NO':>8}")

        print("\nMax-Eigenvalue Test:")
        print(f"{'H0: r=':>12}  {'Max-Eig Stat':>12}  {'CV 5%':>10}  {'Reject?':>8}")
        for r in range(k):
            stat = johansen.lr2[r]
            cv5 = maxeig_cv[r, 1]
            reject = stat > cv5
            print(f"  r = {r:>5}  {stat:>12.4f}  {cv5:>10.4f}  {'YES' if reject else 'NO':>8}")

    # Select cointegrating rank r
    r_selected = 0
    for r in range(k):
        if johansen.lr1[r] > johansen.cvt[r, 1]:
            r_selected = r + 1

    if verbose:
        print(f"\nSelected cointegrating rank: r = {r_selected}")

    if r_selected == 0:
        print("No cointegration found — estimate VAR in differences.")
        return {"johansen_result": johansen, "r_selected": 0,
                "vecm_results": None, "beta": None, "alpha": None}

    # Step 2: VECM estimation
    vecm_model = VECM(df_levels, k_ar_diff=k_ar_diff, coint_rank=r_selected,
                      deterministic="co")  # constant in cointegrating equation
    vecm_result = vecm_model.fit()

    if verbose:
        print("\nVECM Results:")
        print(vecm_result.summary())

    beta = vecm_result.beta           # (k, r) cointegrating vectors
    alpha = vecm_result.alpha         # (k, r) speed of adjustment

    if verbose:
        print("\nCointegrating Vector(s) [normalized]:")
        beta_df = pd.DataFrame(beta, index=variables,
                               columns=[f"CE{i+1}" for i in range(r_selected)])
        print(beta_df.round(4))
        print("\nSpeed of Adjustment (α):")
        alpha_df = pd.DataFrame(alpha, index=variables,
                                columns=[f"CE{i+1}" for i in range(r_selected)])
        print(alpha_df.round(4))

    return {
        "johansen_result": johansen,
        "r_selected": r_selected,
        "vecm_results": vecm_result,
        "beta": beta,
        "alpha": alpha,
    }

Advanced Usage

Bootstrap Confidence Intervals for IRF

def bootstrap_irf(
    var_results,
    periods: int = 12,
    n_boot: int = 500,
    signif: float = 0.05,
    seed: int = 0,
) -> dict:
    """
    Bootstrap percentile confidence intervals for Cholesky IRF.

    Uses residual bootstrap: resample VAR residuals with replacement,
    simulate new data, refit VAR, compute IRF.

    Args:
        var_results: Fitted VARResults object.
        periods:     IRF horizon.
        n_boot:      Number of bootstrap replications.
        signif:      Two-sided CI level (0.05 => 95% CI).
        seed:        Random seed for reproducibility.

    Returns:
        Dictionary with keys: irf_point, irf_lower, irf_upper.
            Each has shape (periods+1, k, k).
    """
    rng = np.random.default_rng(seed)
    k = var_results.neqs
    p = var_results.k_ar
    residuals = var_results.resid  # (T - p, k)
    T_eff = residuals.shape[0]

    # Centered residuals
    resid_centered = residuals - residuals.mean(axis=0)
    Y_orig = var_results.endog       # (T, k)
    coefs = var_results.coefs        # (p, k, k)
    intercept = var_results.intercept  # (k,)

    irf_boots = np.zeros((n_boot, periods + 1, k, k))

    for b in range(n_boot):
        # Resample residuals
        idx = rng.integers(0, T_eff, size=T_eff)
        boot_resid = resid_centered[idx]

        # Simulate new series
        Y_boot = np.zeros((p + T_eff, k))
        Y_boot[:p] = Y_orig[:p]
        for t in range(p, p + T_eff):
            Y_boot[t] = intercept.copy()
            for lag in range(p):
                Y_boot[t] += coefs[lag] @ Y_boot[t - lag - 1]
            Y_boot[t] += boot_resid[t - p]

        df_boot = pd.DataFrame(Y_boot[p:], columns=var_results.names)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            try:
                var_boot = VAR(df_boot).fit(p)
                irf_boot = var_boot.irf(periods)
                irf_boots[b] = irf_boot.orth_irfs
            except Exception:
                irf_boots[b] = np.nan

    alpha_lo = signif / 2
    alpha_hi = 1 - signif / 2
    irf_lower = np.nanquantile(irf_boots, alpha_lo, axis=0)
    irf_upper = np.nanquantile(irf_boots, alpha_hi, axis=0)

    irf_point = var_results.irf(periods).orth_irfs

    return {"irf_point": irf_point, "irf_lower": irf_lower, "irf_upper": irf_upper}


def plot_bootstrap_irf(boot_dict: dict, variable_names: list, output_path: str = None) -> None:
    """
    Plot point-estimate IRF with bootstrap CI bands.

    Args:
        boot_dict:      Output of bootstrap_irf().
        variable_names: List of variable names.
        output_path:    If provided, save figure here.
    """
    k = len(variable_names)
    irf_point = boot_dict["irf_point"]
    irf_lower = boot_dict["irf_lower"]
    irf_upper = boot_dict["irf_upper"]
    H = irf_point.shape[0]
    horizons = np.arange(H)

    fig, axes = plt.subplots(k, k, figsize=(4 * k, 3 * k), sharex=True)
    for i in range(k):
        for j in range(k):
            ax = axes[i, j]
            ax.plot(horizons, irf_point[:, i, j], color="#2980B9", linewidth=2)
            ax.fill_between(horizons, irf_lower[:, i, j], irf_upper[:, i, j],
                            alpha=0.2, color="#2980B9")
            ax.axhline(0, color="black", linewidth=0.8, linestyle="--")
            if i == 0:
                ax.set_title(f"Shock: {variable_names[j]}", fontsize=9)
            if j == 0:
                ax.set_ylabel(f"Response: {variable_names[i]}", fontsize=9)
    fig.suptitle("Bootstrap IRF (95% CI)")
    fig.tight_layout()
    if output_path:
        fig.savefig(output_path, dpi=150)
        print(f"Saved bootstrap IRF to {output_path}")
    plt.show()

Troubleshooting

External Resources

• Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.
• Sims, C.A. (1980). "Macroeconomics and Reality." Econometrica, 48(1), 1–48.
• Johansen, S. (1988). "Statistical Analysis of Cointegration Vectors." JEDS, 12(2–3), 231–254.
• statsmodels VAR documentation: https://www.statsmodels.org/stable/vector_ar.html
• statsmodels VECM documentation: https://www.statsmodels.org/stable/vector_ar.html#vector-error-correction-models-vecm

Examples

Example 1 — Stationary VAR: Full Pipeline

# Generate stationary data
df_stat = generate_var_data(n=400, k=3, cointegrated=False)

# Check stationarity
for col in df_stat.columns:
    print(adf_summary(df_stat[col], name=col).to_string(index=False))

# Fit VAR and test Granger causality
var_out = fit_var_and_granger(df_stat, maxlags=8, ic="aic")

# IRF and FEVD plots
plot_irf_fevd(var_out["results"], periods=12, output_prefix="var_example")

Example 2 — Cointegrated System: Johansen + VECM

# Generate cointegrated I(1) data
df_coint = generate_var_data(n=400, k=3, cointegrated=True)

# Confirm unit roots in levels, stationarity in differences
for col in df_coint.columns:
    print(adf_summary(df_coint[col], name=col).to_string(index=False))

# Johansen test and VECM
vecm_out = johansen_vecm_workflow(df_coint, det_order=0, k_ar_diff=2)

# Confirm error-correction spreads are stationary
if vecm_out["beta"] is not None:
    beta = vecm_out["beta"]
    spread = df_coint.values @ beta[:, 0]
    adf_spread = adfuller(spread, autolag="AIC")
    print(f"\nSpread ADF p-value: {adf_spread[1]:.4f} (should be < 0.05)")

Changelog