Time-Series ML II — Multivariate macro forecasting: VAR / BVAR vs machine learning¶
Unrestricted VAR, Bayesian shrinkage, sparse Lasso-VAR, gradient boosting, and an LSTM¶
The univariate notebook found classical models winning on a single clean series. This one raises the dimension: forecasting six macro variables jointly, where they influence each other over time — the domain of the vector autoregression (VAR). The VAR is the workhorse of empirical macro, but it has a notorious weakness: the parameter count explodes. With $m=6$ variables and $p=6$ lags, each equation has $mp+1=37$ coefficients and the system has $\sim$220 — estimated on a few hundred noisy monthly observations, the unrestricted VAR overfits and forecasts poorly out of sample.
The cure is shrinkage, and it comes in three flavours we can now line up side by side:
- Bayesian VAR with the Minnesota prior — the classic econometric answer (the BVAR arc's
bvar.py, reused here), shrinking coefficients toward a random walk; - Ridge-VAR — L2 shrinkage, the closed-form cousin of the Minnesota prior;
- Lasso-VAR — L1 shrinkage, which selects a sparse set of lag relationships (the frequentist sibling of SSVS-VAR);
against the ML challengers — gradient boosting (a booster per target) and a multivariate LSTM. The question that runs through the subsection: does ML beat the shrinkage-based econometric models when the problem is multivariate? Data: six real US macro series from FRED. Python-lead.
1. The data and the parameter-explosion problem¶
Six monthly US macro series from FRED — industrial production, unemployment, CPI, the fed funds rate, the 10-year Treasury yield, and nonfarm payrolls — transformed to stationarity (log-growth for the quantities and prices, first differences for the rates) and standardised. The goal is to forecast all six one month ahead, jointly. With $p=6$ lags the unrestricted VAR must estimate ~220 coefficients; that is the over-parameterisation shrinkage exists to tame. All models are evaluated by 1-step-ahead out-of-sample RMSE over the last 60 months, averaged across the six (standardised) variables, using the actual lags at each step (so this is clean 1-step accuracy, no error compounding).
import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings, time, sys
warnings.filterwarnings("ignore")
sys.path.insert(0, r"c:\Users\user\project\Z-Multivariate TS-Bayesian Vector Autoregression")
import bvar # reuse the BVAR arc's Minnesota sampler (read-only)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
def load(s): return pd.read_csv(f"{s}.csv",parse_dates=["observation_date"]).set_index("observation_date")[s]
series=["INDPRO","UNRATE","CPIAUCSL","FEDFUNDS","GS10","PAYEMS"]
raw=pd.concat([load(s) for s in series],axis=1).dropna()
X=pd.DataFrame(index=raw.index)
for s in ["INDPRO","CPIAUCSL","PAYEMS"]: X[s]=100*np.log(raw[s]).diff() # log-growth (%)
for s in ["UNRATE","FEDFUNDS","GS10"]: X[s]=raw[s].diff() # first difference
X=X.dropna(); cols=list(X.columns); dates=X.index; Y=X.values; m=Y.shape[1]
p=6; H=60
print(f"{len(Y)} monthly obs ({dates[0].date()} to {dates[-1].date()}), {m} variables, {p} lags")
print(f"Unrestricted VAR({p}) parameters: {m*(m*p+1)} -> estimated on {len(Y)-H-p} training points: over-parameterised")
fig,ax=plt.subplots(2,3,figsize=(14,5))
for a,c in zip(ax.ravel(),cols): a.plot(dates,X[c],color=BLUE,lw=.6); a.set_title(c,fontsize=9); a.axhline(0,color="k",lw=.4)
plt.suptitle("Six stationary macro series (growth rates / differences)"); plt.tight_layout(); plt.show()
862 monthly obs (1954-08-01 to 2026-06-01), 6 variables, 6 lags Unrestricted VAR(6) parameters: 222 -> estimated on 796 training points: over-parameterised
2. The models — unrestricted VAR, three shrinkage estimators, and ML¶
Every model predicts the six-vector one step ahead from the stacked lags. The unrestricted VAR is ordinary least squares equation by equation. The three shrinkage estimators pull those coefficients toward zero / a random walk in different ways — Bayesian Minnesota (posterior-mean coefficients from the Gibbs sampler), Ridge (L2), Lasso (L1, sparse) — and the ML models (gradient boosting per target, a multivariate LSTM) learn the map directly. A persistence baseline sets the floor — the forecast that repeats the last observed value, with no model at all. It is the same object as a random-walk forecast, and it is the number every model in the table has to beat before it has earned anything. One wrinkle worth naming here: these six series are already differenced, so the last observed value is last month's change, and persistence therefore predicts that this month's movement repeats next month. That is a real claim, not a null one — the null is the random walk in levels ("expect no change"), which the next cell adds and which turns out to be the harder benchmark of the two.
from sklearn.linear_model import LinearRegression, RidgeCV, MultiTaskLassoCV
import xgboost as xgb, torch, torch.nn as nn; torch.set_num_threads(2)
mu,sd=Y[:-H].mean(0),Y[:-H].std(0); Z=(Y-mu)/sd
def design(Zz,p):
T=len(Zz); lags=[Zz[p-l:T-l] for l in range(1,p+1)]; return np.hstack([np.ones((T-p,1))]+lags), Zz[p:]
Xd,Yt=design(Z,p); ntr=len(Yt)-H
Xtr,Xte,Ytr,Yte=Xd[:ntr],Xd[ntr:],Yt[:ntr],Yt[ntr:]
def rmse(P): return np.sqrt(np.mean((P-Yte)**2))
def rmse_by(P): return np.sqrt(((P-Yte)**2).mean(0))
res={}; byv={}
res["persistence (last change)"]=rmse(Xte[:,1:1+m]); byv["persistence (last change)"]=rmse_by(Xte[:,1:1+m])
# The series are already differenced, so "repeat last month's CHANGE" is not the natural naive forecast.
# A random walk in LEVELS means predicting zero change, and that is the benchmark the BVAR literature reports.
res["random walk (levels)"]=rmse(np.zeros_like(Yte)); byv["random walk (levels)"]=rmse_by(np.zeros_like(Yte))
Pv=LinearRegression(fit_intercept=False).fit(Xtr,Ytr).predict(Xte); res["VAR (OLS, unrestricted)"]=rmse(Pv); byv["VAR (OLS, unrestricted)"]=rmse_by(Pv)
Pr=RidgeCV(alphas=np.logspace(-1,4,40),fit_intercept=False).fit(Xtr,Ytr).predict(Xte); res["Ridge-VAR (L2)"]=rmse(Pr); byv["Ridge-VAR (L2)"]=rmse_by(Pr)
Pl=MultiTaskLassoCV(cv=5,max_iter=8000).fit(Xtr,Ytr).predict(Xte); res["Lasso-VAR (L1, sparse)"]=rmse(Pl); byv["Lasso-VAR (L1, sparse)"]=rmse_by(Pl)
t=time.time(); fit=bvar.gibbs_minnesota(Z[:ntr+p],p=p,ndraw=1500,burn=500,seed=0); Bm=fit["B"].mean(0)
Pb=Xte@Bm; res["BVAR (Minnesota)"]=rmse(Pb); byv["BVAR (Minnesota)"]=rmse_by(Pb); print(f"BVAR Gibbs {time.time()-t:.0f}s")
Pg=np.column_stack([xgb.XGBRegressor(n_estimators=250,learning_rate=0.05,max_depth=3,verbosity=0).fit(Xtr,Ytr[:,j]).predict(Xte) for j in range(m)])
res["GBM (per-target)"]=rmse(Pg); byv["GBM (per-target)"]=rmse_by(Pg)
# multivariate LSTM (1-step, actual windows)
Ws=np.array([Z[i-p:i] for i in range(p,len(Z))]); Ts=Z[p:]
Wtr,Wte,YtrL,YteL=Ws[:ntr],Ws[ntr:],Ts[:ntr],Ts[ntr:]
torch.manual_seed(0)
class VLSTM(nn.Module):
def __init__(s): super().__init__(); s.l=nn.LSTM(m,48,batch_first=True); s.f=nn.Linear(48,m)
def forward(s,x): o,_=s.l(x); return s.f(o[:,-1])
ml=VLSTM(); opt=torch.optim.Adam(ml.parameters(),5e-3); lf=nn.MSELoss()
Wt=torch.tensor(Wtr,dtype=torch.float32); Yt2=torch.tensor(YtrL,dtype=torch.float32)
for e in range(150): opt.zero_grad(); lf(ml(Wt),Yt2).backward(); opt.step()
with torch.no_grad(): Pn=ml(torch.tensor(Wte,dtype=torch.float32)).numpy()
res["LSTM (multivariate)"]=np.sqrt(np.mean((Pn-YteL)**2)); byv["LSTM (multivariate)"]=np.sqrt(((Pn-YteL)**2).mean(0))
print("\navg 1-step OOS RMSE (6 standardized macro vars, lower=better):")
for k,v in sorted(res.items(),key=lambda z:z[1]): print(f" {k:28s} {v:.4f}")
# 60 one-step forecasts is enough to test differences rather than eyeball them.
_P={"VAR (OLS, unrestricted)":Pv,"Ridge-VAR (L2)":Pr,"Lasso-VAR (L1, sparse)":Pl,
"BVAR (Minnesota)":Pb,"GBM (per-target)":Pg,"random walk (levels)":np.zeros_like(Yte)}
def _nw(x,lag=6):
x=x-x.mean(); v=np.mean(x*x)
for k in range(1,lag+1): v+=2*(1-k/(lag+1))*np.mean(x[k:]*x[:-k])
return np.sqrt(v/len(x))
def _dm(a,b):
d=((_P[a]-Yte)**2).mean(1)-((_P[b]-Yte)**2).mean(1)
return d.mean()/_nw(d)
_ord=["GBM (per-target)","Lasso-VAR (L1, sparse)","BVAR (Minnesota)","Ridge-VAR (L2)","VAR (OLS, unrestricted)"]
_sh=["Ridge-VAR (L2)","Lasso-VAR (L1, sparse)","BVAR (Minnesota)"]
print("\nDiebold-Mariano against the unrestricted VAR (negative = beats it; |stat| > 1.96 to matter):")
for k in _ord[:-1]:
_s=_dm(k,"VAR (OLS, unrestricted)")
print(f" {k:26s} {_s:+6.2f} {'significant' if abs(_s)>1.96 else 'NOT significant'}")
print(f" {'random walk (levels)':26s} {_dm('random walk (levels)','VAR (OLS, unrestricted)'):+6.2f}")
print("\nand among the leaders themselves:")
for a,b in [("GBM (per-target)","Lasso-VAR (L1, sparse)"),("GBM (per-target)","BVAR (Minnesota)"),
("Lasso-VAR (L1, sparse)","BVAR (Minnesota)")]:
_s=_dm(a,b)
print(f" {a.split(' (')[0]:10s} vs {b.split(' (')[0]:10s} {_s:+6.2f} {'significant' if abs(_s)>1.96 else 'indistinguishable'}")
BVAR Gibbs 2s
avg 1-step OOS RMSE (6 standardized macro vars, lower=better): GBM (per-target) 0.5819 Lasso-VAR (L1, sparse) 0.5846 BVAR (Minnesota) 0.5887 Ridge-VAR (L2) 0.5967 VAR (OLS, unrestricted) 0.6249 random walk (levels) 0.6291 persistence (last change) 0.7207 LSTM (multivariate) 0.9131 Diebold-Mariano against the unrestricted VAR (negative = beats it; |stat| > 1.96 to matter): GBM (per-target) -2.00 significant Lasso-VAR (L1, sparse) -2.33 significant BVAR (Minnesota) -2.11 significant Ridge-VAR (L2) -1.76 NOT significant random walk (levels) +0.18 and among the leaders themselves: GBM vs Lasso-VAR -0.23 indistinguishable GBM vs BVAR -0.52 indistinguishable Lasso-VAR vs BVAR -1.35 indistinguishable
3. The scoreboard — shrinkage wins, and ML ties it¶
Three findings, and the first turns on which naive benchmark is the right one. The persistence baseline repeats last month's change, but these series are already differenced, so the natural naive forecast is a random walk in levels — predict zero change. On that yardstick the unrestricted VAR is a dead heat with doing nothing: ~220 estimated coefficients buy essentially no improvement over a forecast of zero. That, rather than its position in the table, is the real indictment of over-parameterisation.
Second, shrinkage helps, though not uniformly by a testable margin. Sixty one-step forecasts are enough to run a Diebold-Mariano test rather than read the ordering off, and against the unrestricted VAR the Lasso and the booster clear the 1.96 threshold while Ridge does not. The direction is right and the mechanism is the familiar one — trading a little bias for a large cut in variance, the same estimation-risk story as the high-dimensional-portfolio and factor-selection notebooks — but a blanket claim that every shrinkage method beats OLS would go beyond what 60 observations establish.
Third, ML is competitive and not dominant, and this part holds up cleanly: gradient boosting lands at the top statistically tied with the shrinkage-linear models, with DM statistics well inside the threshold against both the Lasso-VAR and the BVAR. The LSTM is the exception in the other direction — it finishes last, behind even the naive baselines, which is worth stating plainly rather than describing it as mid-field. Note especially that Lasso-VAR ≈ BVAR — L1 shrinkage and the Bayesian Minnesota prior reach essentially the same accuracy, the frequentist and Bayesian faces of the same idea (and the bridge to SSVS-VAR).
tab=pd.Series(res).sort_values()
def col(k): return (RED if "unrestricted" in k else GREY if "persistence" in k else GREEN if ("BVAR" in k or "Ridge" in k or "Lasso" in k) else BLUE)
fig,ax=plt.subplots(1,2,figsize=(15,4.6))
ax[0].barh(tab.index,tab.values,color=[col(k) for k in tab.index]); ax[0].invert_yaxis(); ax[0].set_xlim(0.55,0.73); ax[0].set_xlabel("avg 1-step OOS RMSE")
ax[0].set_title("Shrinkage (green) beats unrestricted VAR (red); ML (blue) ties")
for i,v in enumerate(tab.values): ax[0].text(v+0.002,i,f"{v:.3f}",va="center",fontsize=8)
# per-variable heatmap for the leading methods
lead=["VAR (OLS, unrestricted)","BVAR (Minnesota)","Lasso-VAR (L1, sparse)","GBM (per-target)","LSTM (multivariate)"]
M=np.array([byv[k] for k in lead])
im=ax[1].imshow(M,cmap="RdYlGn_r",aspect="auto"); ax[1].set_xticks(range(m)); ax[1].set_xticklabels(cols,rotation=45,ha="right",fontsize=8)
ax[1].set_yticks(range(len(lead))); ax[1].set_yticklabels(lead,fontsize=8); ax[1].set_title("Per-variable 1-step RMSE (greener=better)")
for i in range(len(lead)):
for j in range(m): ax[1].text(j,i,f"{M[i,j]:.2f}",ha="center",va="center",fontsize=7)
plt.colorbar(im,ax=ax[1],fraction=0.046); plt.tight_layout(); plt.show()
print(f"Unrestricted VAR {res['VAR (OLS, unrestricted)']:.3f} -> BVAR {res['BVAR (Minnesota)']:.3f}, Lasso-VAR {res['Lasso-VAR (L1, sparse)']:.3f}"
f" (shrinkage helps); GBM {res['GBM (per-target)']:.3f} ties them. Per-variable: no method dominates every series.")
Unrestricted VAR 0.625 -> BVAR 0.589, Lasso-VAR 0.585 (shrinkage helps); GBM 0.582 ties them. Per-variable: no method dominates every series.
4. Predicted vs actual, and the horizon the exercise has been avoiding¶
Two things the scoreboard cannot show.
First, what the forecasts actually look like. An RMSE averaged over six variables hides which series a model handles and which it does not, and whether the errors are noise or systematic bias. Below, each variable gets its own panel: the realised path solid, the two leading forecasts dashed, plus the binned predicted-versus-actual view the rest of the collection uses.
Second, and more substantively: everything so far is one step ahead. That is the easiest horizon there is, and it flatters every model equally — with the actual lags supplied at each step, nothing compounds and no model has to propagate its own errors. It is also not what a forecaster is usually asked for. So we repeat the comparison at horizons of 1, 3, 6 and 12 months using direct forecasts (a separate fit per horizon, predicting $y_{t+h}$ from information at $t$), which keeps every model on the same footing and avoids the recursive-substitution bias that complicated the univariate notebook.
# ---------- predicted vs actual, per variable ----------
_lead=[("BVAR (Minnesota)",Pb,GREEN),("GBM (per-target)",Pg,BLUE)]
fig,ax=plt.subplots(2,3,figsize=(15,6.4))
_td=dates[-H:]
for j,a in enumerate(ax.ravel()):
a.plot(_td,Yte[:,j],color="black",lw=1.3,label="actual")
for _n,_P,_c in _lead: a.plot(_td,_P[:,j],color=_c,lw=1.1,ls="--",label=_n)
a.set_title(cols[j],fontsize=9); a.axhline(0,color="k",lw=.3)
if j==0: a.legend(fontsize=7)
plt.suptitle("One-step forecasts (dashed) against the realised path (solid), standardised units")
plt.tight_layout(); plt.show()
fig,ax=plt.subplots(1,2,figsize=(13,4.3))
for _n,_P,_c in _lead:
_f=_P.ravel(); _a=Yte.ravel()
_q=np.quantile(_f,np.linspace(0,1,11)); _b=np.clip(np.digitize(_f,_q[1:-1]),0,9)
ax[0].plot([_f[_b==k].mean() for k in range(10)],[_a[_b==k].mean() for k in range(10)],"o-",color=_c,lw=2,label=_n)
_lm=[Yte.min(),Yte.max()]
ax[0].plot(_lm,_lm,"k--",lw=1,label="perfect"); ax[0].set_xlabel("predicted (standardised)")
ax[0].set_ylabel("mean actual"); ax[0].set_title("Proportions vs predictions (all 6 variables pooled)"); ax[0].legend(fontsize=8)
for _n,_P,_c in _lead:
ax[1].bar(np.arange(m)+(0.2 if "GBM" in _n else -0.2),rmse_by(_P),width=0.4,color=_c,label=_n)
ax[1].set_xticks(range(m)); ax[1].set_xticklabels(cols,rotation=45,ha="right",fontsize=8)
ax[1].set_ylabel("1-step RMSE"); ax[1].set_title("Per-variable accuracy"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
for _n,_P,_c in _lead:
_f=_P.ravel(); _a=Yte.ravel()
print(f"{_n:20s} slope of actual on predicted {np.polyfit(_f,_a,1)[0]:.3f} (1.0 = unbiased)")
print("A slope below 1 means the forecast moves more than its information content justifies -- the Mincer-Zarnowitz")
print("benchmark used in the LSTM notebook. Both models sit under it here, which is the usual finding for macro")
print("one-step forecasts and is invisible in RMSE.")
# ---------- horizons beyond one step ----------
print("\n" + "="*100)
print("Direct multi-step: a separate fit per horizon, predicting h months ahead from information at t.")
_HZ=[1,3,6,12]
_curve={k:[] for k in ["VAR (OLS)","Ridge-VAR","Lasso-VAR","GBM","random walk"]}
for _h in _HZ:
_Xh=Xd[:-_h] if _h>0 else Xd
_Yh=Yt[_h-1:] if _h>1 else Yt
_Xh=Xd[:len(Xd)-(_h-1)] ; _Yh=Yt[_h-1:]
_ntr=len(_Yh)-H
_Xtr,_Xte2,_Ytr,_Yte2=_Xh[:_ntr],_Xh[_ntr:],_Yh[:_ntr],_Yh[_ntr:]
def _r(Pm): return float(np.sqrt(np.mean((Pm-_Yte2)**2)))
_curve["VAR (OLS)"].append(_r(LinearRegression(fit_intercept=False).fit(_Xtr,_Ytr).predict(_Xte2)))
_curve["Ridge-VAR"].append(_r(RidgeCV(alphas=np.logspace(-1,4,40),fit_intercept=False).fit(_Xtr,_Ytr).predict(_Xte2)))
_curve["Lasso-VAR"].append(_r(MultiTaskLassoCV(cv=5,max_iter=8000).fit(_Xtr,_Ytr).predict(_Xte2)))
_curve["GBM"].append(_r(np.column_stack([xgb.XGBRegressor(n_estimators=250,learning_rate=0.05,max_depth=3,
verbosity=0).fit(_Xtr,_Ytr[:,j]).predict(_Xte2) for j in range(m)])))
_curve["random walk"].append(_r(np.zeros_like(_Yte2)))
print(f" {'horizon':>8}" + "".join(f"{k:>14s}" for k in _curve))
for i,_h in enumerate(_HZ):
print(f" {_h:>8}" + "".join(f"{_curve[k][i]:>14.4f}" for k in _curve))
fig,ax=plt.subplots(figsize=(7.5,4.4))
for k,v in _curve.items():
ax.plot(_HZ,v,"o-",lw=2,label=k,color=(GREY if "walk" in k else BLUE if k=="GBM" else GREEN if "Lasso" in k else ORANGE if "Ridge" in k else RED))
ax.set_xlabel("forecast horizon (months)"); ax.set_ylabel("OOS RMSE"); ax.set_title("Everything converges on the random walk as the horizon grows"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
_gap1=_curve['random walk'][0]-min(_curve[k][0] for k in _curve if k!='random walk')
_gap12=_curve['random walk'][-1]-min(_curve[k][-1] for k in _curve if k!='random walk')
print(f"\nAt h=1 the best model beats a random walk by {_gap1:.4f}; at h=12 the margin is {_gap12:.4f}.")
print("That is the honest frame for everything above. One-step accuracy is the horizon at which these models look")
print("best, and it is the horizon this notebook has been using throughout. Push the horizon out and the macro")
print("dynamics that shrinkage estimates so carefully stop mattering: the series are close to random walks at a")
print("year's remove, and every method converges toward the naive forecast. Shrinkage buys real short-horizon")
print("accuracy; it does not buy long-horizon predictability that is not in the data.")
BVAR (Minnesota) slope of actual on predicted 0.762 (1.0 = unbiased) GBM (per-target) slope of actual on predicted 0.698 (1.0 = unbiased) A slope below 1 means the forecast moves more than its information content justifies -- the Mincer-Zarnowitz benchmark used in the LSTM notebook. Both models sit under it here, which is the usual finding for macro one-step forecasts and is invisible in RMSE. ==================================================================================================== Direct multi-step: a separate fit per horizon, predicting h months ahead from information at t.
horizon VAR (OLS) Ridge-VAR Lasso-VAR GBM random walk
1 0.6249 0.5967 0.5846 0.5819 0.6291
3 0.6512 0.6108 0.6078 0.6213 0.6291
6 0.6786 0.6152 0.6215 0.7261 0.6291
12 1.4152 0.6744 0.6449 0.9759 0.6291
At h=1 the best model beats a random walk by 0.0472; at h=12 the margin is -0.0159. That is the honest frame for everything above. One-step accuracy is the horizon at which these models look best, and it is the horizon this notebook has been using throughout. Push the horizon out and the macro dynamics that shrinkage estimates so carefully stop mattering: the series are close to random walks at a year's remove, and every method converges toward the naive forecast. Shrinkage buys real short-horizon accuracy; it does not buy long-horizon predictability that is not in the data.
5. Summary¶
On multivariate macro forecasting, shrinkage is the decisive ingredient — and machine learning ties the shrinkage-based econometric models rather than beating them. The unrestricted VAR over-parameterises (~220 coefficients on a few hundred noisy points) and forecasts no better than a random walk in levels — the plainest possible statement of what over-parameterisation costs; Bayesian Minnesota, Lasso-VAR and Ridge-VAR all improve on it by shrinking, and Lasso-VAR ≈ BVAR shows the frequentist L1 and Bayesian Minnesota routes converging. Gradient boosting matches them at the top — statistically tied, on a Diebold-Mariano test, with both the Lasso-VAR and the BVAR — so ML is now competitive, where on the single clean series of ex1 the classical models led it outright. The multivariate LSTM is the clear failure of the exercise: it finishes last, worse than a random walk, on a 150-step full-batch fit that is far too little training for a recurrent net and far too little data to justify one.
That is the progression this subsection is tracing:
- ex1 (one clean series) — classical (SARIMA/ETS) wins decisively;
- ex2 (this, multivariate macro) — shrinkage beats unrestricted VAR; ML ties the shrinkage-linear models;
- ex3 (many series, next) — ML (one global model) is expected to win.
The unifying thread is shrinkage/estimation-risk, exactly as in the high-dimensional-portfolio and factor-selection notebooks — Bayesian priors (Minnesota, your BVAR arc), L1 selection (your SSVS-VAR), and L2 all tame the same coefficient explosion, and ML's tree/recurrent flexibility buys little extra when the macro dynamics are near-linear. Next: global / panel forecasting, where one ML model across many related series turns the tables.