Multivariate Time Series in Python — Core Concepts

When to go multivariate

Use multivariate analysis when:

• Multiple variables evolve together and influence each other (GDP, unemployment, inflation).
• You want to understand lead-lag relationships between variables.
• Forecasting one variable improves when you include others.
• You need to model the entire system, not just one output.

VAR: Vector AutoRegression

VAR is the multivariate extension of AR models. Instead of one equation, you have a system where each variable depends on past values of all variables:

y₁ₜ = c₁ + a₁₁y₁,ₜ₋₁ + a₁₂y₂,ₜ₋₁ + ε₁ₜ y₂ₜ = c₂ + a₂₁y₁,ₜ₋₁ + a₂₂y₂,ₜ₋₁ + ε₂ₜ

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

# DataFrame with multiple columns, DatetimeIndex
data = pd.DataFrame({
    "gdp_growth": gdp_series,
    "unemployment": unemp_series,
    "inflation": inflation_series,
})

model = VAR(data)

# Select optimal lag order
lag_order = model.select_order(maxlags=12)
print(lag_order.summary())

# Fit with selected lags
fitted = model.fit(lag_order.aic)
forecast = fitted.forecast(data.values[-fitted.k_ar:], steps=8)

Choosing the lag order

VAR models are sensitive to the number of lags. Too few miss important dynamics; too many waste parameters. Use information criteria:

lag_selection = model.select_order(maxlags=12)
print(f"AIC suggests: {lag_selection.aic} lags")
print(f"BIC suggests: {lag_selection.bic} lags")
# BIC usually suggests fewer lags — better for forecasting

Granger causality

Granger causality tests whether past values of one variable help predict another, beyond what the other variable’s own past provides. It is a statistical test, not true causality.

from statsmodels.tsa.stattools import grangercausalitytests

# Does 'temperature' Granger-cause 'ice_cream_sales'?
test_result = grangercausalitytests(
    data[["ice_cream_sales", "temperature"]],  # target first, predictor second
    maxlag=7,
)

# Check p-values at each lag
for lag, results in test_result.items():
    f_test_p = results[0]["ssr_ftest"][1]
    print(f"Lag {lag}: p-value = {f_test_p:.4f}")

If the p-value is below 0.05, temperature Granger-causes ice cream sales at that lag — meaning temperature history adds predictive value beyond sales history alone.

Cointegration

Two non-stationary series are cointegrated if a linear combination of them is stationary. Classic example: a dog on a leash in a park. Both the dog and the owner wander randomly (non-stationary), but the distance between them (the leash) stays bounded (stationary).

from statsmodels.tsa.stattools import coint

# Test if two series are cointegrated
score, p_value, _ = coint(series_a, series_b)
print(f"Cointegration p-value: {p_value:.4f}")
# p < 0.05 → evidence of cointegration

When variables are cointegrated, use a VECM (Vector Error Correction Model) instead of VAR on differences:

from statsmodels.tsa.vector_ar.vecm import VECM

model = VECM(data, k_ar_diff=2, coint_rank=1)
fitted = model.fit()
forecast = fitted.predict(steps=10)

The coint_rank is the number of cointegrating relationships — use the Johansen test to determine it:

from statsmodels.tsa.vector_ar.vecm import coint_johansen

result = coint_johansen(data, det_order=0, k_ar_diff=2)
print("Trace statistics:", result.lr1)
print("Critical values (5%):", result.cvt[:, 1])

Impulse Response Functions

IRFs show how a shock to one variable propagates through the system over time:

irf = fitted.irf(periods=20)

# Plot how a shock to GDP affects unemployment
irf.plot(impulse="gdp_growth", response="unemployment")

This is essential for policy analysis — “if GDP drops by 1%, how does unemployment respond over the next 8 quarters?”

Forecast Error Variance Decomposition

FEVD shows what percentage of each variable’s forecast error is explained by shocks to each other variable:

fevd = fitted.fevd(periods=10)
fevd.summary()
fevd.plot()

If 60% of the forecast error for sales comes from shocks to advertising (not sales itself), then advertising is a dominant driver.

Key assumptions and preparation

All variables in a VAR model must be stationary. The standard workflow:

1. Test each variable for stationarity (ADF test).
2. Difference non-stationary variables — unless they are cointegrated (use VECM instead).
3. Fit the VAR model on the stationary data.
4. Check residuals for remaining autocorrelation.

Common misconception

Granger causality is not real causality. If ice cream sales Granger-cause drowning deaths, it does not mean ice cream causes drowning — both are caused by hot weather. Granger causality only means one variable has predictive power for another, which can be due to confounding variables.

The one thing to remember: Multivariate time series analysis captures the dynamic interactions between variables — VAR models the system, Granger causality identifies predictive relationships, cointegration reveals long-run equilibria, and impulse response functions trace how shocks propagate.