BVAR mortality forecasting I — the Minnesota-prior model¶
A Bayesian VAR on US mortality, following Meseguer (2010), estimated with a from-scratch Gibbs sampler (no PyMC). This first notebook covers the Minnesota-prior model; the steady-state (Villani mean-adjusted) model follows separately.
Design (matching Meseguer and our Lee-Carter benchmark):
- Model the first differences of log mortality $\Delta\log m_t$ (changes), 22 abridged age groups, separately by sex.
- Sample 1933-2000, out-of-sample forecast 2001-2024.
- Lag order $p=1$ (paper), with a lag-length comparison $p=1,2,3$ at the end.
- Benchmark: the Lee-Carter forecasts from
LC_Mortality.ipynb.
Data: US series from the Human Mortality Database (mortality.org) and the Human Fertility Database (humanfertility.org).
1. Data — mortality changes and their correlation structure¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from hmd_io import read_mx, read_deaths, read_population
from mortality import (aggregate_to_groups, life_table_abridged, e0_series, ex_series,
GROUP_LABEL, GROUP_LO, lee_carter, forecast_lee_carter, project_rates)
import bvar
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
DIR = "."
SAMPLE, FYEARS = (1933, 2000), np.arange(2001, 2025)
COLOR = {"female": "#d1495b", "male": "#00798c"}
mx = read_mx(f"{DIR}/US Mx_1x1.txt")
deaths = read_deaths(f"{DIR}/US Deaths_1x1.txt")
pop = read_population(f"{DIR}/US Population.txt")
g = {s: aggregate_to_groups(mx, deaths, pop, s) for s in ["female", "male"]}
logm = {s: np.log(g[s]) for s in g}
# modelled series: Delta log m over the sample
Y = {s: np.diff(logm[s].loc[SAMPLE[0]:SAMPLE[1]].values, axis=0) for s in g}
print(f"Delta log m sample matrix: {Y['female'].shape} (years x 22 age groups)")
Delta log m sample matrix: (67, 22) (years x 22 age groups)
1.1 The correlation structure — why the "spatial trick"¶
The 22 mortality-change series do not move independently: adjacent age groups are strongly correlated, distant ones weakly. The heatmap below (female $\Delta\log m$) shows this temporal-contiguity band. The Minnesota prior exploits it via $\lambda_2(i,j)=0.8\,\mathrm{corr}(\Delta\log y_i,\Delta\log y_j)$: the lower the correlation, the smaller $\lambda_2$, the harder that cross-coefficient is shrunk toward zero. The prior thus learns the age-adjacency structure from the data.
Cf = np.corrcoef(Y["female"].T)
fig, ax = plt.subplots(figsize=(6.4, 5.4))
im = ax.imshow(Cf, cmap="RdBu_r", vmin=-1, vmax=1)
ax.set_xticks(range(22)); ax.set_yticks(range(22))
ax.set_xticklabels(GROUP_LABEL, rotation=90, fontsize=6)
ax.set_yticklabels(GROUP_LABEL, fontsize=6)
ax.set_title("Correlation of $\Delta\log m$ across age groups (female)")
fig.colorbar(im, fraction=0.046, pad=0.04)
plt.show()
print(f"e.g. corr(25-29, 30-34) = {Cf[6,7]:.2f} vs corr(25-29, 90-94) = {Cf[6,20]:.2f}")
<>:7: SyntaxWarning: invalid escape sequence '\D'
<>:7: SyntaxWarning: invalid escape sequence '\D'
C:\Users\user\AppData\Local\Temp\ipykernel_37288\2231276805.py:7: SyntaxWarning: invalid escape sequence '\D'
ax.set_title("Correlation of $\Delta\log m$ across age groups (female)")
e.g. corr(25-29, 30-34) = 0.82 vs corr(25-29, 90-94) = 0.31
2. The Minnesota-prior BVAR and its Gibbs sampler¶
Model. For one sex, the 22 age-group mortality changes $y_t=\Delta\log m_t$ follow a VAR($p$): $$y_t = c + \Pi_1 y_{t-1}+\dots+\Pi_p y_{t-p}+\varepsilon_t,\qquad \varepsilon_t\sim N(0,\Sigma).$$ Stacking $t=p{+}1,\dots,T$ gives $Y = XB + E$ with $X=[\,1,\ y_{t-1},\dots,y_{t-p}\,]$ ($(T{-}p)\times k$, $k=1+mp$), $B$ ($k\times m$) collecting the intercept $c$ and the $\Pi_l$, and rows of $E\sim N(0,\Sigma)$.
Prior — coefficients $\mathrm{vec}(B)\sim N(\mathrm{vec}(B_0),V_0)$, $V_0$ diagonal:
- Means $B_0$: the own first lag is centred on 0, everything else on 0. On differenced data a zero own-lag coefficient is a random walk in log-mortality levels; the drift is carried by the (diffuse) intercept.
- Variances $V_0$ (Minnesota, Litterman 1986 / Meseguer eq. 16): prior sd of the coefficient on variable $j$, lag $l$, equation $i$, $$\sigma_{ijl}=\begin{cases}\lambda_1/l^{\lambda_3}, & j=i,\\[1mm]\big(\lambda_1\,\lambda_2(i,j)/l^{\lambda_3}\big)\,(s_i/s_j), & j\ne i,\end{cases}$$ with $s_i$ the residual sd of an AR($p$) fit (units scaling); the intercept is diffuse. $\lambda_1$ = overall tightness (0.2); $\lambda_3$ = lag decay (1 = harmonic, inert at $p{=}1$).
Spatial trick. $\lambda_2(i,j)=0.8\,\mathrm{corr}(\Delta\log y_i,\Delta\log y_j)$ — the correlation structure of Section 1.1 baked into the prior, so contiguous ages interact and distant ages decouple.
Prior — covariance $\Sigma\sim\mathcal{IW}(S_0,\nu_0)$, diffuse.
Gibbs sampler (2 blocks). Alternate the two closed-form full conditionals:
- $\mathrm{vec}(B)\mid\Sigma,Y\sim N(\bar b,\bar V)$, with $\bar V^{-1}=V_0^{-1}+(\Sigma^{-1}\otimes X'X)$ and $\bar b=\bar V[V_0^{-1}\mathrm{vec}(B_0)+\mathrm{vec}(X'Y\Sigma^{-1})]$ (drawn via the Cholesky factor of the precision).
- $\Sigma\mid B,Y\sim\mathcal{IW}(S_0+E'E,\ \nu_0+T-p)$, $E=Y-XB$.
Discard the burn-in; keep $\{B^{(g)},\Sigma^{(g)}\}$.
Forecast. Simulate each draw forward, propagating parameter and shock
uncertainty, then cumulate $\Delta\log m\to\log m$, exponentiate, and run the
abridged life table for $e_0/e_x$ credible bands. Implemented in bvar.py
(gibbs_minnesota, forecast).
2.1 The four Minnesota hyperparameters (Meseguer 2010, eqs. 16-17)¶
Every prior standard deviation is governed by four hyperparameters $\lambda_1,\dots,\lambda_4$. For the coefficient on variable $j$ at lag $l$ in equation $i$, and for the intercept (or any exogenous regressor) $\phi_i$:
$$\sigma_{ijl}=\begin{cases}\dfrac{\lambda_1}{l^{\lambda_3}}, & j=i\\[2mm]\dfrac{\lambda_1\,\lambda_2(i,j)}{l^{\lambda_3}}\cdot\dfrac{s_i}{s_j}, & j\ne i\end{cases}\qquad\qquad \sigma_{\phi_i}=\lambda_4\,s_i .$$
$\lambda_1$ — overall tightness (we use 0.2). The prior sd on the first own lag ($i{=}j,\,l{=}1$). Decreasing $\lambda_1$ shrinks the diagonal of $\Pi_1$ toward its prior mean (1 in levels, 0 in our differenced data) and every other coefficient toward zero — it sets how tightly the random-walk specification is imposed.
$\lambda_2(i,j)$ — lag weight (cross-variable) (spatial trick). Scales the prior sd of other variables' lags relative to own lags. With $0<\lambda_2<1$ the off-diagonal $\Pi_l$ elements shrink toward zero, so a series is explained more by its own lags than by others; $\lambda_2=1$ treats all lags equally. Writing it as $\lambda_2(i,j)$ lets each pair carry its own weight — which we set to $0.8\,\mathrm{corr}(\Delta\log y_i,\Delta\log y_j)$: the lower the correlation between age groups $i$ and $j$, the smaller $\lambda_2(i,j)$, the harder that cross-coefficient is shrunk toward zero.
$\lambda_3$ — lag decay (we use 1). Shrinks coefficients toward their prior means faster at higher lags: $\lambda_3=1$ is harmonic decay, larger $\lambda_3$ downweights distant lags more, "giving greater importance to more recent lags." Inert at $p{=}1$; active in the $p=2,3$ comparison (Section 5).
$\lambda_4$ — intercept / exogenous tightness (large = diffuse). Sets the prior sd on the intercept as $\lambda_4 s_i$. Following standard practice the Minnesota prior is made non-informative on the intercept (zero mean, large $\lambda_4$), so the drift — the average pace of mortality decline — is left to the data rather than the prior.
3. Estimation and convergence¶
NDRAW, BURN = 2500, 600
fits = {}
import time
for s in ["female", "male"]:
t = time.time()
fits[s] = bvar.gibbs_minnesota(Y[s], p=1, lam1=0.2, lam3=1.0, ndraw=NDRAW, burn=BURN, seed=1)
print(f"{s}: {NDRAW} draws in {time.time()-t:.1f}s (B {fits[s]['B'].shape}, Sigma {fits[s]['Sigma'].shape})")
female: 2500 draws in 38.4s (B (2500, 23, 22), Sigma (2500, 22, 22))
male: 2500 draws in 42.5s (B (2500, 23, 22), Sigma (2500, 22, 22))
# convergence: trace of a few parameters (kept draws should look stationary)
fs = fits["female"]
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
axes[0].plot(fs["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$ (age 0 var)")
axes[1].plot(fs["B"][:, 1, 6], lw=.5); axes[1].set_title("cross lag age 0 -> 25-29") # row0=intercept, so row1 = lag-1 coef on var0 (age 0)
axes[2].plot(fs["B"][:, 7, 6], lw=.5); axes[2].set_title("own lag, eq 25-29") # row7 = lag-1 coef on var6 (25-29)
for a in axes: a.set_xlabel("kept draw"); a.grid(alpha=.2)
fig.suptitle("Trace plots (female, Minnesota BVAR p=1)", y=1.03); fig.tight_layout(); plt.show()
<>:4: SyntaxWarning: invalid escape sequence '\S'
<>:4: SyntaxWarning: invalid escape sequence '\S'
C:\Users\user\AppData\Local\Temp\ipykernel_37288\3929709419.py:4: SyntaxWarning: invalid escape sequence '\S'
axes[0].plot(fs["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$ (age 0 var)")
4. Forecasts vs. actual and the Lee-Carter benchmark¶
Posterior predictive $e_0$ and $e_{65}$ (median and 90% band) against the actual 2001-2024 values and the Lee-Carter forecast. Both carry 90% intervals, but they mean different things: the BVAR band propagates parameter + shock uncertainty from the full posterior, whereas the Lee-Carter band is a plug-in interval reflecting only the random-walk uncertainty in the index $k_t$ (the $\alpha_x,\beta_x$ are treated as fixed). The LC band is therefore expected to be narrower — it understates uncertainty. Both are anchored at the observed 2000 log-rates.
def ex_bands(fit, s, age, H=len(FYEARS), thin=2):
draws = bvar.forecast(fit, Y[s], H)[::thin]
idx = list(GROUP_LO).index(age)
base = logm[s].loc[2000].values
cum = np.cumsum(draws, axis=1)
nd = draws.shape[0]
E = np.empty((nd, H))
for d in range(nd):
for h in range(H):
E[d, h] = life_table_abridged(np.exp(base + cum[d, h]), s)["ex"][idx]
return np.percentile(E, [5, 50, 95], axis=0)
# Lee-Carter benchmark (jump-off adjusted) with 90% plug-in interval from k_t uncertainty
def lc_ex_bands(s, age):
fit = lee_carter(np.log(g[s].loc[SAMPLE[0]:SAMPLE[1]])); fc = forecast_lee_carter(fit, len(FYEARS))
jo = np.log(g[s].loc[2000].values); idx = list(GROUP_LO).index(age)
def e_at(kv):
R = project_rates(fit, kv, jo)
return np.array([life_table_abridged(R[i], s)["ex"][idx] for i in range(len(FYEARS))])
mid = e_at(fc["k"])
hi = e_at(fc["k"] - 1.645 * fc["se"]) # low k -> low mortality -> high e
lo = e_at(fc["k"] + 1.645 * fc["se"])
return lo, mid, hi
fig, axes = plt.subplots(2, 2, figsize=(12, 7.5), sharex=True)
for col, s in enumerate(["female", "male"]):
for row, age in enumerate([0, 65]):
ax = axes[row, col]
lo, mid, hi = ex_bands(fits[s], s, age)
llo, lmid, lhi = lc_ex_bands(s, age)
act = ex_series(g[s], s, age)
ax.plot(act.loc[1980:2000].index, act.loc[1980:2000], color="gray", lw=1)
ax.plot(FYEARS, act.reindex(FYEARS), color=COLOR[s], lw=2, label="actual")
ax.fill_between(FYEARS, lo, hi, color="k", alpha=.12, label="BVAR 90%")
ax.plot(FYEARS, mid, "--", color="k", lw=1.3, label="BVAR median")
ax.fill_between(FYEARS, llo, lhi, color=COLOR[s], alpha=.10, label="Lee-Carter 90%")
ax.plot(FYEARS, lmid, ":", color=COLOR[s], lw=1.6, label="Lee-Carter median")
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set_title(f"{s} e{age}"); ax.grid(alpha=.25)
if row == 0 and col == 0: ax.legend(fontsize=7.5)
axes[1, 0].set_xlabel("Year"); axes[1, 1].set_xlabel("Year")
fig.suptitle("Minnesota BVAR forecasts vs actual vs Lee-Carter", y=1.0); fig.tight_layout(); plt.show()
# out-of-sample e0 MAE (pre-COVID 2001-2019), BVAR median vs Lee-Carter
print("Out-of-sample e0 MAE, 2001-2019 (years):")
for s in ["female", "male"]:
bv = ex_bands(fits[s], s, 0)[1]
lc = lc_ex_bands(s, 0)[1]
act = ex_series(g[s], s, 0).reindex(FYEARS).to_numpy()
n = 19
print(f" {s}: BVAR {np.abs(bv[:n]-act[:n]).mean():.2f} Lee-Carter {np.abs(lc[:n]-act[:n]).mean():.2f}")
Out-of-sample e0 MAE, 2001-2019 (years):
female: BVAR 0.22 Lee-Carter 0.32
male: BVAR 0.32 Lee-Carter 0.32
4.1 Log-mortality forecasts at selected ages¶
Beyond the $e_0/e_{65}$ summaries, the model forecasts each age group's rate. Four representative female series below (BVAR and Lee-Carter median + 90% band, vs actual). Each panel tells a different demographic story:
- Age 0 (infants). Mortality fell steeply for decades but the decline decelerated after ~2000 — the actual line flattens above both model medians, which extrapolate the historical pace. Much of the easy progress (perinatal care, the post-1990s drop in SIDS) had already been banked.
- Age 25-29 (young adults). The cautionary panel: mortality rose after ~2010, driven by the drug-overdose / opioid epidemic and other "deaths of despair" (suicide, alcohol-related). Every model, trained on decades of decline, forecast the opposite, and the actual line walks clean out of both 90% bands. No purely extrapolative model can anticipate a trend reversal like this.
- Age 65-69 (early retirement). A steady decline, largely from cardiovascular gains, tracked well by both models — apart from the sharp 2020-21 COVID spike, which sits outside the bands.
- Age 85-89 (oldest-old). Mortality kept falling, if noisily; the forecasts track the level but tend to under-predict the improvement at the very oldest ages.
Throughout, the wide BVAR band (full posterior: parameter + shock uncertainty) versus the narrower Lee-Carter band (plug-in: only $k_t$ uncertainty, with $\alpha_x,\beta_x$ fixed) shows LC understating uncertainty — most starkly at 25-29.
def logm_bands(fit, s):
dr = bvar.forecast(fit, Y[s], len(FYEARS))
lm = logm[s].loc[2000].values + np.cumsum(dr, axis=1) # cumulate changes onto observed 2000
return np.percentile(lm, [5, 50, 95], axis=0) # (3, H, 22)
def lc_logm_bands(s):
fit = lee_carter(np.log(g[s].loc[SAMPLE[0]:SAMPLE[1]])); fc = forecast_lee_carter(fit, len(FYEARS))
jo = np.log(g[s].loc[2000].values)
lm = lambda kv: np.log(project_rates(fit, kv, jo))
return lm(fc["k"] - 1.645 * fc["se"]), lm(fc["k"]), lm(fc["k"] + 1.645 * fc["se"]) # low/mid/high log m
s = "female"
b = logm_bands(fits[s], s)
lclo, lcmid, lchi = lc_logm_bands(s)
act = np.log(g[s].loc[2001:2024]).values
hist = np.log(g[s].loc[1970:2000])
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for ax, a in zip(axes.ravel(), ["0", "25-29", "65-69", "85-89"]):
j = GROUP_LABEL.index(a)
ax.plot(hist.index, hist.values[:, j], color="gray", lw=1, label="history")
ax.plot(FYEARS, act[:, j], "k-", lw=2, label="actual")
ax.fill_between(FYEARS, b[0, :, j], b[2, :, j], color=COLOR[s], alpha=.15, label="BVAR 90%")
ax.plot(FYEARS, b[1, :, j], "--", color=COLOR[s], lw=1.6, label="BVAR median")
ax.fill_between(FYEARS, lclo[:, j], lchi[:, j], color="#555", alpha=.12, label="Lee-Carter 90%")
ax.plot(FYEARS, lcmid[:, j], ":", color="#555", lw=1.5, label="Lee-Carter median")
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set_title(f"log m, age {a} (female)"); ax.grid(alpha=.25)
axes[0, 0].legend(fontsize=7.5)
fig.suptitle("Minnesota BVAR log-mortality forecasts at selected ages", y=1.0)
fig.tight_layout(); plt.show()
5. Lag length — $p=1,2,3$¶
The Minnesota shrinkage lets us add lags without over-fitting. Below, the female $e_0$ forecast under $p=1,2,3$ (fewer draws at higher $p$). The lag-decay $\lambda_3=1$ shrinks higher lags harder, so the forecasts stay close.
fig, ax = plt.subplots(figsize=(9, 4))
act = ex_series(g["female"], "female", 0)
ax.plot(act.loc[1990:2000].index, act.loc[1990:2000], color="gray", lw=1)
ax.plot(FYEARS, act.reindex(FYEARS), color=COLOR["female"], lw=2, label="actual")
for p, nd, c in [(1, 2000, "#222"), (2, 1200, "#e07b39"), (3, 700, "#3a7d44")]:
fp = bvar.gibbs_minnesota(Y["female"], p=p, lam1=0.2, lam3=1.0, ndraw=nd, burn=300, seed=2)
lo, mid, hi = ex_bands(fp, "female", 0)
ax.plot(FYEARS, mid, "--", color=c, lw=1.4, label=f"BVAR p={p}")
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="Female e0: Minnesota BVAR at p=1,2,3", xlabel="Year", ylabel="e0 (years)")
ax.grid(alpha=.25); ax.legend(fontsize=8)
plt.show()
6. Summary¶
- The Minnesota-prior BVAR is estimated with a hand-written 2-block Gibbs sampler on $\Delta\log m$ (22 groups, per sex), with the spatial trick encoding age-adjacency via $\lambda_2(i,j)=0.8\,\mathrm{corr}$.
- Out-of-sample $e_0$/$e_{65}$ forecasts come with full posterior credible bands (parameter + shock uncertainty), unlike the plug-in Lee-Carter intervals.
- Forecasts are robust to lag length thanks to the Minnesota shrinkage.
Next: the steady-state (Villani mean-adjusted) BVAR with the informative SSA Alt-2 prior, then a head-to-head of both BVARs against Lee-Carter.