Vector Autoregression (VAR) Model

A Vector Autoregression (VAR) Model is a multivariate time-series model where each variable is predicted from its own past values and the past values of all other variables. It captures dynamic interdependencies, enabling forecasting and causal-style analysis via impulse responses and Granger causality. Example: jointly modeling inflation, unemployment, and interest rates to forecast and assess policy shocks.

Imagine a group of friends whose moods affect each other: how you feel today depends a bit on how you felt yesterday, and also on how your friends felt yesterday. A Vector Autoregression (VAR) model works like that for several things measured over time—like inflation, unemployment, and interest rates. Instead of predicting one timeline by itself, it predicts a set of timelines together. Each value is estimated using past values of itself and the past values of the other related series. In AI and statistics, VAR helps capture “feedback loops” between multiple time-based signals.

import numpy as np  # numerical arrays for simulation
import pandas as pd  # time-indexed tabular data
from statsmodels.tsa.api import VAR  # VAR estimation and forecasting tools

# --- Create a realistic multivariate time series (e.g., inflation, unemployment, interest rate) ---
np.random.seed(7)  # make results reproducible
n = 180  # monthly observations (~15 years)
dates = pd.date_range("2010-01-01", periods=n, freq="MS")  # month-start index

# Simulate three correlated macro series with cross-lag effects (a VAR-like data-generating process)
eps = np.random.multivariate_normal(mean=[0, 0, 0], cov=[[1.0, 0.4, 0.2],
                                                        [0.4, 1.0, 0.3],
                                                        [0.2, 0.3, 1.0]], size=n)
Y = np.zeros((n, 3))  # container for the 3 series
A1 = np.array([[0.6, 0.1, 0.0],   # lag-1 effects across all variables
               [0.2, 0.5, 0.1],
               [0.1, 0.2, 0.4]])
for t in range(1, n):  # generate a stable multivariate process
    Y[t] = A1 @ Y[t - 1] + eps[t]  # cross-variable lag dependence + correlated shocks

df = pd.DataFrame(Y, index=dates, columns=["inflation", "unemployment", "policy_rate"])  # name series

# --- Fit a VAR model and forecast jointly ---
train = df.iloc[:-12]  # keep last 12 months for out-of-sample evaluation
test = df.iloc[-12:]   # "future" values we will compare against

model = VAR(train)  # set up the multivariate time series model
res = model.fit(maxlags=6, ic="aic")  # *** KEY LINE *** fit VAR with lag order chosen by AIC

fc = res.forecast(y=train.values[-res.k_ar:], steps=len(test))  # multi-step joint forecast
fc = pd.DataFrame(fc, index=test.index, columns=df.columns)  # align forecast with dates/columns

rmse = np.sqrt(((fc - test) ** 2).mean())  # per-series forecast error summary
print("Selected lag order (k_ar):", res.k_ar)  # show chosen lag length
print("RMSE by series:\n", rmse.round(3))  # quick diagnostic for each variable

• The key output is the fitted statsmodels.tsa.api.VAR result (res) and the joint forecasts from VARResults.forecast, producing a multivariate prediction for all series at once.
• fit(maxlags=6, ic="aic") searches lag orders up to 6 and selects the one minimizing AIC; changing maxlags can materially affect stability and forecast quality.
• Common mistake: forecasting with the wrong “initial” lag window—VARResults.forecast requires the last k_ar observations (train.values[-res.k_ar:]), not the whole training set.
• In ML workflows, the forecast matrix can be used as multi-target features/targets (e.g., scenario simulation), and the residuals can be checked for remaining structure before downstream decision models.

library(vars)      # VAR estimation, lag selection, forecasting
library(ggplot2)   # quick plotting for forecast vs actual

# --- Create a realistic multivariate time series (e.g., inflation, unemployment, interest rate) ---
set.seed(7)  # make results reproducible
n <- 180  # monthly observations (~15 years)
dates <- seq(as.Date("2010-01-01"), by = "month", length.out = n)  # month index

# Simulate correlated shocks for 3 series
Sigma <- matrix(c(1.0, 0.4, 0.2,
                  0.4, 1.0, 0.3,
                  0.2, 0.3, 1.0), nrow = 3, byrow = TRUE)
eps <- MASS::mvrnorm(n = n, mu = c(0, 0, 0), Sigma = Sigma)  # correlated innovations

# Generate a stable VAR(1)-like process with cross-lag effects
A1 <- matrix(c(0.6, 0.1, 0.0,
               0.2, 0.5, 0.1,
               0.1, 0.2, 0.4), nrow = 3, byrow = TRUE)
Y <- matrix(0, nrow = n, ncol = 3)  # container for the 3 series
for (t in 2:n) {  # build time dependence
  Y[t, ] <- A1 %*% Y[t - 1, ] + eps[t, ]  # cross-variable lag dependence + shocks
}

df <- data.frame(date = dates,
                 inflation = Y[, 1],
                 unemployment = Y[, 2],
                 policy_rate = Y[, 3])

# --- Fit a VAR model and forecast jointly ---
train <- df[1:(n - 12), c("inflation", "unemployment", "policy_rate")]  # training window
test  <- df[(n - 11):n, c("inflation", "unemployment", "policy_rate")]  # holdout window

lag_sel <- VARselect(train, lag.max = 6, type = "const")  # compute IC-based lag suggestions
p <- lag_sel$selection[["AIC(n)"]]  # pick the AIC-selected lag order

fit <- VAR(train, p = p, type = "const")  # *** KEY LINE *** fit VAR(p) with chosen lag length
pred <- predict(fit, n.ahead = nrow(test))  # multi-step forecasts for each series

fc <- cbind(pred$fcst$inflation[, "fcst"],
            pred$fcst$unemployment[, "fcst"],
            pred$fcst$policy_rate[, "fcst"])
colnames(fc) <- colnames(test)

rmse <- sqrt(colMeans((fc - as.matrix(test))^2))  # per-series forecast error summary
print(paste("Selected lag order (p):", p))
print(round(rmse, 3))

• The key output is the fitted vars::VAR object (fit) and the multi-series forecasts from predict, which returns a forecast path for each variable conditional on the others.
• vars::VARselect provides multiple information criteria; this example uses AIC via selection[["AIC(n)"]], but BIC often chooses smaller (more regularized) lag orders.
• Common mistake: mixing date columns into the model matrix—vars::VAR expects only numeric series, so the code subsets to the three value columns.
• In AI/ML usage, the forecasted trajectories can feed planning/optimization, and the selected lag order can be treated as a hyperparameter tuned via rolling-origin evaluation.