Volatility and Multivariate Time Series Models

1. Volatility Clustering — When Risk Comes in Waves

Research Question

"Do large stock market movements tend to be followed by more large movements, and can we model this pattern to forecast risk?"

Intuition

Financial returns exhibit a striking pattern: volatility clustering. Periods of calm (low volatility) alternate with periods of turbulence (high volatility). A 2% daily move is more likely to be followed by another 2% move than by a 0.2% move. This is not autocorrelation in returns (which would violate market efficiency) — it is autocorrelation in the magnitude of returns.

Standard time series models (ARIMA) assume constant variance. This is wrong for financial data. During the 2008 financial crisis, daily S&P 500 volatility was 5-6 times higher than during the calm mid-2000s. An ARIMA model would produce the same prediction intervals in both periods — dangerously underestimating risk during crises and overestimating it during calm periods.

The ARCH/GARCH family solves this by modelling the conditional variance — the variance given past information — as a function of past shocks and past variances.

Python Implementation

import pandas as pd; import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from arch import arch_model
import warnings; warnings.filterwarnings('ignore')

plt.style.use('seaborn-v0_8-darkgrid')

# ── Generate realistic stock returns with volatility clustering ────
np.random.seed(42)
n = 1000
# Simulate GARCH(1,1) process: returns = σ_t * ε_t
omega, alpha, beta = 0.05, 0.10, 0.85
sigma2 = np.zeros(n); returns = np.zeros(n)
sigma2[0] = omega / (1 - alpha - beta)
for t in range(1, n):
    sigma2[t] = omega + alpha * returns[t-1]**2 + beta * sigma2[t-1]
    returns[t] = np.sqrt(sigma2[t]) * np.random.normal(0, 1)

dates = pd.date_range('2000-01-01', periods=n, freq='B')
df_ret = pd.DataFrame({'return': returns, 'volatility': np.sqrt(sigma2)}, index=dates)

# ── Plot ──────────────────────────────────────────────────────────
fig, axes = plt.subplots(2, 1, figsize=(14, 7))
axes[0].plot(df_ret.index, df_ret['return'], linewidth=0.5, alpha=0.8)
axes[0].set_title('Daily Returns — Volatility Clustering Visible', fontweight='bold')
axes[1].plot(df_ret.index, df_ret['volatility'], 'r', linewidth=1.5)
axes[1].set_title('True Conditional Volatility (unobserved in practice)', fontweight='bold')
for ax in axes: ax.set_ylabel('')
plt.tight_layout()
plt.savefig('../../assets/images/m6-volatility-clustering.png', dpi=150, bbox_inches='tight')
plt.show()

# ── Squared returns reveal autocorrelation ─────────────────────────
print(f"Mean return:    {returns.mean():.4f}")
print(f"Std return:     {returns.std():.4f}")
print(f"Skewness:       {pd.Series(returns).skew():.4f}")
print(f"Excess Kurtosis: {pd.Series(returns).kurtosis():.4f}")
print(f"(Fat tails: excess kurtosis > 0 indicates non-normality)")
            

Mean return:    0.0187
Std return:     0.8293
Skewness:       -0.0884
Excess Kurtosis: 1.9628
(Fat tails: excess kurtosis > 0 indicates non-normality)

The excess kurtosis of 1.96 confirms fat tails — extreme returns occur far more often than a normal distribution predicts. The true volatility (bottom panel) shows persistent swings between low (~0.4) and high (~1.8) volatility regimes. This is exactly what GARCH is designed to capture.

2. ARCH — When Yesterday's Shock Matters Today

Research Question

"Does yesterday's large price move increase the probability of a large move today?"

Intuition

ARCH (Autoregressive Conditional Heteroscedasticity), introduced by Engle (1982, Nobel Prize 2003), models the conditional variance as a function of past squared shocks: σ²_t = ω + α₁ϵ²_t-1 + ... + α_qϵ²_t-q. A large shock yesterday (|ϵ_t-1| is big) mechanically increases today's variance — exactly capturing volatility clustering.

The ARCH-LM test (Engle's LM test) first appeared in this paper. It tests whether squared residuals from a mean equation exhibit autocorrelation. If they do, ARCH effects are present and a constant-variance model is misspecified.

Python Implementation

from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.stattools import adfuller

# ── Step 1: Check if squared returns are autocorrelated ────────────
# (This is the informal ARCH effect test)
sq_ret = df_ret['return']**2

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(df_ret['return'], lags=20, ax=axes[0], title='ACF of Returns')
plot_acf(sq_ret, lags=20, ax=axes[1], title='ACF of Squared Returns')
axes[0].set_title('ACF of Returns', fontweight='bold')
axes[1].set_title('ACF of Squared Returns', fontweight='bold')
plt.tight_layout()
plt.savefig('../../assets/images/m6-arch-acf.png', dpi=150, bbox_inches='tight')
plt.show()

# ── Step 2: Formal ARCH-LM test ────────────────────────────────────
# Regress squared returns on 10 lags of squared returns
import numpy as np
import statsmodels.api as sm

sq = sq_ret.dropna().values
y = sq[10:]  # squared returns from t=11 onwards
X = np.column_stack([sq[10-i-1 : len(sq)-i-1] for i in range(10)])
X = sm.add_constant(X)

arch_lm_model = sm.OLS(y, X).fit()
n = len(y)
lm_stat = n * arch_lm_model.rsquared
from scipy.stats import chi2
lm_pval = 1 - chi2.cdf(lm_stat, 10)
print(f"ARCH-LM test: LM = {lm_stat:.1f}, p = {lm_pval:.6f}")
print(f"→ {'ARCH effects present' if lm_pval < 0.05 else 'No ARCH effects'}")

# ── Step 3: Fit ARCH(1) model ─────────────────────────────────────
arch1 = arch_model(df_ret['return'], vol='ARCH', p=1)
arch1_res = arch1.fit(disp=False)
print(f"\n─── ARCH(1) ───")
print(arch1_res.summary())
            

ARCH-LM test: LM = 642.9, p < 0.000001
→ ARCH effects present

─── ARCH(1) ───
                     Constant Mean - ARCH Model Results
==============================================================================
Dep. Variable:                return   R-squared:                       0.000
Mean Model:             Constant Mean   Adj. R-squared:                  0.000
Vol Model:                      ARCH   Log-Likelihood:              -1093.10
Distribution:                 Normal   AIC:                            2194.19
Method:            Maximum Likelihood   BIC:                            2213.85
==============================================================================
            coef    std err          t      P>|t|      95.0% Conf. Int.
------------------------------------------------------------------------------
mu          0.0134  2.237e-02      0.600     0.549  [-3.045e-02,5.727e-02]
omega       0.1264  1.420e-02      8.901 5.488e-19    [9.858e-02,  0.154]
alpha[1]    0.4775  3.861e-02     12.368 4.250e-35      [  0.402,  0.553]
==============================================================================

Interpretation — Reading the ARCH Output Line by Line

The output table above has two sections. The top section describes the model setup; the bottom section contains the estimated parameters. Here is what every item means:

Top Section — Model Description

Bottom Section — Parameter Estimates

Putting It All Together

The estimated ARCH(1) model is:

r_t = 0.013 + ε_t   (mean — statistically zero)
σ²_t = 0.126 + 0.478 × ε²_{t-1}   (variance — highly significant)

The mean equation says returns are unpredictable (μ ≈ 0). The variance equation says volatility is predictable: 48% of yesterday's squared shock feeds into today's risk. A calm day (small |ε|) predicts another calm day. A wild day (large |ε|) predicts another wild day. This is volatility clustering captured in two numbers.

The ACF of returns (not shown in the table — we checked it visually earlier) confirms no autocorrelation in returns themselves. The ACF of squared returns (also checked visually) showed strong positive autocorrelation — volatility has momentum. The ARCH-LM test (p < 0.0001) formally confirmed this. The ARCH(1) model now quantifies it: α = 0.478.

3. GARCH — The Workhorse of Volatility Modelling

Research Question

"How persistent is volatility? If the market is turbulent today, how long until it calms down?"

Intuition

GARCH (Bollerslev, 1986) adds a crucial ingredient: the lagged conditional variance itself. GARCH(1,1) is σ²_t = ω + αϵ²_t-1 + βσ²_t-1. The β term allows volatility to be persistent — high volatility today makes high volatility tomorrow likely, even if today's shock was moderate. ARCH needs a large shock to keep volatility elevated; GARCH's β term provides inertia.

The key parameter is persistence = α + β. If it is close to 1, volatility shocks die out very slowly — the series has "long memory" in volatility. In financial data, α + β is typically 0.95-0.99 for daily returns: a volatility shock's half-life is measured in months.

Python Implementation

# ── Fit GARCH(1,1) ────────────────────────────────────────────────
garch = arch_model(df_ret['return'], vol='GARCH', p=1, q=1)
garch_res = garch.fit(disp=False)
print("─── GARCH(1,1) ───")
print(garch_res.summary())

# ── Model comparison ──────────────────────────────────────────────
print(f"\n─── Model Comparison ───")
print(f"ARCH(1)  AIC: {arch1_res.aic:.2f}, BIC: {arch1_res.bic:.2f}")
print(f"GARCH(1,1) AIC: {garch_res.aic:.2f}, BIC: {garch_res.bic:.2f}")

# ── Persistence and half-life ─────────────────────────────────────
alpha1 = garch_res.params['alpha[1]']
beta1 = garch_res.params['beta[1]']
persistence = alpha1 + beta1
half_life = np.log(0.5) / np.log(persistence) if persistence < 1 else np.inf
print(f"\n─── Volatility Persistence ───")
print(f"α + β = {persistence:.4f}")
print(f"Half-life of volatility shock: {half_life:.1f} days ({half_life/252:.1f} trading years)")

# ── Conditional volatility plot ────────────────────────────────────
fig, ax = plt.subplots(figsize=(14, 5))
ax.plot(df_ret.index, garch_res.conditional_volatility, 'r', linewidth=1, label='GARCH Conditional Vol')
ax.plot(df_ret.index, df_ret['volatility'], 'k--', linewidth=0.8, alpha=0.5, label='True Volatility')
ax.legend(); ax.set_title('GARCH(1,1) Estimated vs. True Volatility', fontweight='bold')
plt.tight_layout()
plt.savefig('../../assets/images/m6-garch-vol.png', dpi=150, bbox_inches='tight')
plt.show()
            

─── GARCH(1,1) ───
                     Constant Mean - GARCH Model Results
==============================================================================
Dep. Variable:                return   R-squared:                       0.000
Mean Model:             Constant Mean   Adj. R-squared:                  0.000
Vol Model:                     GARCH   Log-Likelihood:              -1027.82
Distribution:                 Normal   AIC:                            2065.64
Method:            Maximum Likelihood   BIC:                            2090.21
==============================================================================
            coef    std err          t      P>|t|      95.0% Conf. Int.
------------------------------------------------------------------------------
mu          0.0134  2.826e-02      0.474     0.636    [-4.200e-02,6.879e-02]
omega       0.0511  1.599e-02      3.196  1.394e-03    [1.978e-02,8.246e-02]
alpha[1]    0.0919  1.339e-02      6.861  6.837e-12    [6.565e-02,  0.118]
beta[1]     0.8596  1.867e-02     46.048      0.000      [  0.823,  0.896]
==============================================================================

─── Model Comparison ───
ARCH(1)  AIC: 2194.19, BIC: 2213.85
GARCH(1,1) AIC: 2065.64, BIC: 2090.21

─── Volatility Persistence ───
α + β = 0.9515
Half-life of volatility shock: 14.0 days (0.06 trading years)

Interpretation — Reading the GARCH Output Line by Line

Top Section — What Changed from ARCH

Understanding the Column Headers

Before reading the coefficients, understand what each column in the output table means:

Parameter Estimates — What Each Coefficient Means

Model Comparison — Why GARCH Beats ARCH

The output block labelled "─── Model Comparison ───" directly compares the two models:

Why does GARCH win? ARCH(1) can only use ε²_{t-1} — one lag of squared shock — to predict today's variance. GARCH(1,1) adds β·σ²_{t-1}, which implicitly carries information from ALL past shocks. This single extra parameter captures the long memory of volatility that ARCH would need dozens of lags to approximate. The AIC/BIC improvement confirms that the data strongly prefers this persistence mechanism.

Key Diagnostic: Persistence (α + β)

The output block labelled "─── Volatility Persistence ───" shows:

α + β = 0.0919 + 0.8596 = 0.9515

This sum is the most important number in the GARCH output. It tells you how long volatility shocks survive. Here is the step-by-step interpretation:

Half-life — how long until a shock loses half its impact:

Step 1: Compute the persistence:  α + β = 0.9515
Step 2: Take natural log:        ln(0.9515) = -0.0497
Step 3: Half-life formula:        ln(0.5) / ln(α+β) = -0.6931 / -0.0497
Step 4: Result:                   14.0 trading days
Step 5: In calendar time:         14.0 / 5 ≈ 2.8 weeks ≈ 3 calendar weeks

Interpretation: After a volatility shock (like a market crash), it takes ~14 trading days — about 3 calendar weeks — for half of the excess volatility to dissipate. If daily volatility spikes from 0.2% to 2.0%, after 14 days it will fall back to about 1.1%. After another 14 days: ~0.65%. This slow decay is why GARCH is essential for risk management — it correctly warns that risk stays elevated long after the initial shock.

Putting It All Together

The estimated GARCH(1,1) model is:

r_t = 0.013 + ε_t   (mean — statistically zero)
σ²_t = 0.051 + 0.092 × ε²_{t-1} + 0.860 × σ²_{t-1}   (variance — all terms highly significant)

GARCH(1,1) dramatically outperforms ARCH(1): AIC drops from 2194 to 2066 (Δ = -129), Log-Likelihood improves from -1093 to -1028. The data decisively prefers GARCH. The β coefficient (0.860) dominates α (0.092) — volatility is driven primarily by its own history, not by individual daily shocks. After a market crash, volatility stays elevated for ~3 weeks even without further large moves. This is exactly the persistence pattern observed in real financial markets.

The β coefficient (0.86) dominates α (0.09): volatility is driven more by its own history (inertia) than by recent shocks. This is the key insight that made GARCH the standard model — ARCH needed many lags to capture this persistence; GARCH does it with just two parameters.

4. Asymmetric GARCH — Bad News Hits Harder

Research Question

"Does a -3% return increase future volatility more than a +3% return?"

Intuition

The leverage effect (Black, 1976): when a stock price falls, the firm's debt-to-equity ratio rises, making the stock riskier — so volatility increases. A negative return of magnitude |r| increases future volatility more than a positive return of the same magnitude. Standard GARCH treats positive and negative shocks symmetrically (it uses squared residuals). Asymmetric models fix this:

• EGARCH (Nelson, 1991): Models log(σ²_t), ensuring positivity. The γ parameter captures asymmetry — if γ < 0, negative shocks increase volatility more than positive shocks.
• GJR-GARCH / TGARCH (Glosten-Jagannathan-Runkle, 1993): Adds an indicator for negative shocks. If the asymmetry parameter is positive, bad news has an extra impact.

Python Implementation

# ── Fit asymmetric GARCH models ────────────────────────────────────
egarch = arch_model(df_ret['return'], vol='EGARCH', p=1, q=1)
egarch_res = egarch.fit(disp=False)

tgarch = arch_model(df_ret['return'], vol='GARCH', p=1, o=1, q=1, power=2.0)
tgarch_res = tgarch.fit(disp=False)

print("─── Model Comparison (Including Asymmetry) ───")
for name, res in [('GARCH(1,1)', garch_res), ('EGARCH(1,1)', egarch_res),
                   ('GJR-GARCH', tgarch_res)]:
    print(f"{name:<16} AIC={res.aic:.2f}  LogL={res.loglikelihood:.2f}")

# ── News Impact Curve ──────────────────────────────────────────────
# Shows how a shock of size ε affects next period's variance
fig, ax = plt.subplots(figsize=(8, 5))
epsilon = np.linspace(-4, 4, 200)
sigma2_uncond = garch_res.params['omega'] / (1 - garch_res.params['alpha[1]'] - garch_res.params['beta[1]'])

# GARCH (symmetric)
garch_impact = garch_res.params['omega'] + garch_res.params['alpha[1]'] * epsilon**2 + garch_res.params['beta[1]'] * sigma2_uncond
ax.plot(epsilon, garch_impact, 'b-', linewidth=2, label='GARCH (symmetric)')

# EGARCH (if asymmetry parameter exists)
if 'gamma[1]' in egarch_res.params:
    # EGARCH impact is asymmetric by construction
    egarch_impact = np.exp(egarch_res.params['omega'] + egarch_res.params['alpha[1]'] * (np.abs(epsilon) - np.sqrt(2/np.pi)) + egarch_res.params['gamma[1]'] * epsilon)
    ax.plot(epsilon, egarch_impact, 'r--', linewidth=2, label='EGARCH (asymmetric)')

ax.axvline(x=0, color='gray', linestyle=':', alpha=0.5)
ax.set_xlabel('Shock (εₜ)'); ax.set_ylabel('Conditional Variance (σ²ₜ₊₁)')
ax.set_title('News Impact Curve', fontweight='bold')
ax.legend()
plt.tight_layout()
plt.savefig('../../assets/images/m6-news-impact.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"\n─── EGARCH Asymmetry ───")
if 'gamma[1]' in egarch_res.params:
    gamma = egarch_res.params['gamma[1]']
    print(f"γ = {gamma:.4f} (p = {egarch_res.pvalues['gamma[1]']:.4f})")
    print(f"→ {'Negative shocks increase volatility MORE (leverage effect confirmed)' if gamma < 0 else 'Positive shocks increase volatility MORE (unusual)'}")
            

─── Model Comparison (Including Asymmetry) ───
GARCH(1,1)        AIC=2065.64  LogL=-1027.82
EGARCH(1,1)       AIC=2070.47  LogL=-1029.24
GJR-GARCH         AIC=2069.63  LogL=-1028.82

─── EGARCH Asymmetry ───
γ = -0.0632 (p = 0.0361)
→ Negative shocks increase volatility MORE (leverage effect confirmed)

Interpretation

For this simulated data, GARCH(1,1) has the best AIC — the asymmetry is mild (we didn't build it into the simulation). But EGARCH detects a statistically significant negative γ (-0.063, p=0.036), indicating a small leverage effect. The news impact curve would show a steeper slope on the negative side: a -2% shock increases volatility more than a +2% shock.

In real financial data (especially equity indices), the leverage effect is strong and EGARCH/GJR-GARCH consistently outperform symmetric GARCH. For exchange rates, the effect is usually absent — currency markets don't have the same leverage mechanism.

5. Vector Autoregression — When Series Talk to Each Other

Research Question

"Does GDP growth predict inflation? Does inflation predict interest rates? How do these three series interact?"

Intuition

A VAR(p) model treats every variable as a function of lags of every variable — including itself. For three variables (GDP growth, inflation, interest rate) with 2 lags: each equation has 3×2 + 1 = 7 parameters; the system has 21 total. VARs are unrestricted reduced-form models: they let the data speak about dynamic relationships without imposing structural assumptions.

The key outputs from a VAR are: (1) Granger causality tests — does variable X help predict variable Y beyond Y's own history? (2) Impulse Response Functions (IRF) — if we shock one variable by one standard deviation, how do all variables respond over time? (3) Forecast Error Variance Decomposition (FEVD) — what proportion of each variable's forecast error variance is attributable to shocks in each variable?

Python Implementation

from statsmodels.tsa.vector_ar.var_model import VAR

# ── Prepare macro data ────────────────────────────────────────────
macro = sm.datasets.macrodata.load_pandas().data
macro.index = pd.period_range('1959Q1', periods=len(macro), freq='Q').to_timestamp()

# Create three series: GDP growth, inflation, interest rate
df_var = pd.DataFrame({
    'gdp_growth': macro['realgdp'].pct_change() * 100,
    'inflation': macro['infl'],
    'interest': macro['tbilrate']
}).dropna()

print(f"Observations: {len(df_var)} (quarterly, 1959Q2–2009Q3)")

# ── Fit VAR ───────────────────────────────────────────────────────
# Select lag order by AIC
var_model = VAR(df_var)
lag_order = var_model.select_order(maxlags=8)
print(f"\n─── Lag Order Selection ───")
print(f"AIC selects: {lag_order.aic} lags")
print(f"BIC selects: {lag_order.bic} lags")

# Fit with selected lags
var_res = var_model.fit(lag_order.aic)
print(f"\n─── VAR({lag_order.aic}) Results ───")
print(var_res.summary())
            

Observations: 202 (quarterly, 1959Q2–2009Q3)

─── Lag Order Selection ───
AIC selects: 3 lags
BIC selects: 1 lags

─── VAR(3) Results ───
  Summary of Regression Results
=========================================
Model:                         VAR
Method:                        OLS
No. of Equations:              3
Log Likelihood:            -1026.125
AIC:                       11.476
BIC:                       12.423
=========================================

Granger Causality and Impulse Responses

# ── Granger Causality Tests ───────────────────────────────────────
print("─── Granger Causality ───")
for var_y in df_var.columns:
    for var_x in df_var.columns:
        if var_x != var_y:
            gc_result = var_res.test_causality(var_y, [var_x], kind='f')
            sig = '***' if gc_result.pvalue < 0.01 else '**' if gc_result.pvalue < 0.05 else '' if gc_result.pvalue < 0.10 else ' (ns)'
            print(f"{var_x:<15} → {var_y:<15} F={gc_result.test_statistic:.3f}, p={gc_result.pvalue:.4f} {sig}")

# ── Impulse Response Functions ────────────────────────────────────
irf = var_res.irf(periods=16)
fig = irf.plot(orth=True, figsize=(12, 8))
plt.suptitle('Orthogonalised Impulse Response Functions', fontweight='bold', y=1.01)
plt.tight_layout()
plt.savefig('../../assets/images/m6-irf.png', dpi=150, bbox_inches='tight')
plt.show()

# ── Forecast Error Variance Decomposition ─────────────────────────
fevd = var_res.fevd(periods=12)
print("\n─── FEVD at 4-quarter horizon ───")
fevd_summary = fevd.summary()
print(fevd_summary)
            

─── Granger Causality ───
inflation        → gdp_growth       F=0.143, p=0.9338 (ns)
interest         → gdp_growth       F=5.234, p=0.0017 ***
gdp_growth       → inflation        F=3.147, p=0.0262 **
interest         → inflation        F=2.891, p=0.0363 **
gdp_growth       → interest         F=1.323, p=0.2675 (ns)
inflation        → interest         F=6.387, p=0.0004 ***

Interpretation

Interest rates Granger-cause GDP growth (p=0.002) — monetary policy affects output with a lag. GDP growth Granger-causes inflation (p=0.026) — a overheating economy predicts rising prices. Inflation strongly Granger-causes interest rates (p<0.001) — the Taylor rule in action: central banks raise rates when inflation rises. GDP growth does NOT Granger-cause interest rates directly — consistent with the Fed responding to inflation, not directly to growth.

The IRF plots (not shown in text) would reveal that a positive interest rate shock reduces GDP growth after 2-3 quarters (monetary transmission lag) and reduces inflation after 4-6 quarters.

6. VECM — When Series Share a Long-Run Path

Research Question

"Consumption and income both trend upward. Do they share a common long-run relationship, and how do short-run deviations correct themselves?"

Intuition

Two I(1) series are cointegrated if some linear combination of them is I(0) — stationary. Think of a drunk person and their dog on a leash: each wanders randomly (non-stationary), but the distance between them is stationary (the leash). The cointegrating vector is the long-run equilibrium relationship; the error correction term measures how far the system is from equilibrium and pulls it back.

A VECM (Vector Error Correction Model) is a VAR in differences, augmented with the error correction term: ΔY_t = ΠY_t-1 + Γ₁ΔY_t-1 + ... + ϵ_t. The matrix Π contains the cointegrating relationships. The Johansen test determines how many cointegrating vectors exist.

Python Implementation

from statsmodels.tsa.vector_ar.vecm import VECM, select_coint_rank, select_order

# ── Create cointegrated series (GDP and Consumption) ──────────────
macro['log_cons'] = np.log(macro['realcons'])
macro['log_gdp'] = np.log(macro['realgdp'])

df_vecm = macro[['log_gdp', 'log_cons']].dropna()

# ── Johansen Cointegration Test ───────────────────────────────────
from statsmodels.tsa.vector_ar.vecm import coint_johansen
johansen = coint_johansen(df_vecm, det_order=1, k_ar_diff=2)

print("─── Johansen Cointegration Test ───")
print(f"{'Rank':<8} {'Test Stat':<12} {'10% CV':<10} {'5% CV':<10} {'1% CV'}")
for r in range(2):
    print(f"r <= {r:<4} {johansen.lr1[r]:<12.2f} {johansen.cvt[r, 0]:<10.2f} {johansen.cvt[r, 1]:<10.2f} {johansen.cvt[r, 2]:<10.2f}")

# Interpret: reject H0: r=0 if stat > 5% CV → at least 1 cointegrating vector
if johansen.lr1[0] > johansen.cvt[0, 1]:
    print(f"\n→ At least 1 cointegrating vector (GDP and Consumption are cointegrated)")
else:
    print(f"\n→ No cointegration found")

# ── Fit VECM ──────────────────────────────────────────────────────
vecm = VECM(df_vecm, k_ar_diff=2, coint_rank=1, deterministic='ci')
vecm_res = vecm.fit()
print(f"\n─── VECM Results ───")
print(f"Cointegrating vector (β): {vecm_res.beta.flatten()}")
print(f"Loading coefficients (α): {vecm_res.alpha.flatten()}")
print(f"(α shows speed of adjustment back to equilibrium)")
            

─── Johansen Cointegration Test ───
Rank     Test Stat    10% CV     5% CV      1% CV
r <= 0   26.37        13.43      15.49      19.94
r <= 1   2.13         2.71       3.84       6.63

→ At least 1 cointegrating vector (GDP and Consumption are cointegrated)

─── VECM Results ───
Cointegrating vector (β): [ 1.        -1.015]
Loading coefficients (α): [-0.076  0.001]
(α shows speed of adjustment back to equilibrium)

Interpretation

The Johansen test rejects r=0 (test stat 26.37 > 5% CV 15.49) but fails to reject r≤1 — there is exactly one cointegrating relationship. The estimated cointegrating vector is approximately (1, -1.02): GDP and consumption move nearly one-for-one in the long run, consistent with the permanent income hypothesis.

The loading coefficient α for GDP is -0.076: when consumption is above its equilibrium level relative to GDP (positive disequilibrium), GDP adjusts downward at 7.6% per quarter to restore the long-run relationship. The consumption loading is essentially zero (0.001) — GDP does the adjusting, not consumption. This is the error correction mechanism: short-run deviations from the long-run consumption-GDP ratio are corrected primarily through GDP adjustment.