BVAR mortality forecasting II — the steady-state model¶
Companion to BVAR_Mortality.ipynb (the Minnesota model). Same data, sample,
from-scratch Gibbs sampler, and Lee-Carter benchmark — the difference is the
parametrization and an informative long-run prior:
- Model the first differences of log mortality $\Delta\log m_t$, 22 abridged age groups, separately by sex; sample 1933-2000, forecast 2001-2024.
- Steady-state (Villani mean-adjusted) VAR with the same Minnesota prior on the AR coefficients, plus a Normal prior on the steady state $\Psi$ centred on the SSA Alternative-2 ultimate mortality-improvement assumptions.
- Estimated with a 3-block Gibbs sampler.
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, time
from hmd_io import read_mx, read_deaths, read_population
from mortality import (aggregate_to_groups, life_table_abridged, 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}
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)
The Minnesota prior on the AR coefficients is unchanged, so the spatial trick still applies: $\lambda_2(i,j)=0.8\,\mathrm{corr}(\Delta\log y_i,\Delta\log y_j)$, tightening cross-coefficients between weakly-correlated (distant) age groups. The correlation band (adjacent ages move together) is the same as in the Minnesota notebook; what is new here is the treatment of the long run.
2. The steady-state (mean-adjusted) BVAR and its Gibbs sampler¶
Model. Villani's mean-adjusted form writes the VAR around its unconditional mean (steady state) $\Psi$: $$y_t = \Psi + \Pi_1(y_{t-1}-\Psi)+\dots+\Pi_p(y_{t-p}-\Psi)+\varepsilon_t,\qquad \varepsilon_t\sim N(0,\Sigma).$$ Since $y_t=\Delta\log m_t$, $\Psi$ is the steady-state vector of mortality changes — the ultimate annual rate of mortality improvement by age. A stationary VAR's long-horizon forecasts converge to $\Psi$, so a prior placed directly on $\Psi$ is an intuitive way to encode a view about the long run — much easier to reason about than the intercept $c$ of the reduced form (which has no simple meaning). That is the whole point of the mean-adjusted parametrization.
Prior. The same Minnesota prior on the AR coefficients $\Pi_l$ as in the Minnesota notebook (own-lag mean 0, overall tightness $\lambda_1$, spatial cross-shrinkage $\lambda_2(i,j)$, lag decay $\lambda_3$). The new component is a Normal prior on the steady state, $\Psi\sim N(\psi_0,\Lambda_0)$ with $\Lambda_0$ diagonal (Section 2.1). $\Sigma\sim\mathcal{IW}(S_0,\nu_0)$ as before.
Gibbs sampler (3 blocks) — Villani (2009). Alternate the three closed-form conditionals:
- $\Pi\mid\Psi,\Sigma,Y\sim N$: regress the demeaned data $z_t=y_t-\Psi$ on its own lags with the Minnesota prior (precision $V_0^{-1}+\Sigma^{-1}\otimes Z'Z$).
- $\Psi\mid\Pi,\Sigma,Y\sim N$: with $u_t=y_t-\sum_l\Pi_l y_{t-l}$ and $D=I-\sum_l\Pi_l$, the model is $u_t=D\Psi+\varepsilon_t$ — a Gaussian regression for $\Psi$ combined with the prior $N(\psi_0,\Lambda_0)$.
- $\Sigma\mid\Pi,\Psi,Y\sim\mathcal{IW}(S_0+E'E,\ \nu_0+T-p)$ on the demeaned-VAR residuals $E$.
Discard the burn-in; keep $\{\Pi^{(g)},\Psi^{(g)},\Sigma^{(g)}\}$. Forecasting
simulates each draw forward around its own $\Psi$. Implemented in bvar.py
(gibbs_steady_state, forecast).
2.1 The informative steady-state prior (SSA Alternative 2)¶
We centre $\Psi$ on the SSA Trustees' Alternative-2 ultimate mortality-improvement assumptions (2024 Report): an age-graded annual rate of decline that reaches its ultimate value by ~2048 — roughly 1.6%/yr under age 15, ~0.74%/yr at 65+, and ~0.82%/yr overall on a death-weighted basis (~0.9%/yr as an unweighted group mean), declining with age. In the model's units the prior mean is $\psi_0(\text{age})=-\text{rate}/100$ (a negative number — mortality falling). The prior sd $\lambda_\Psi$ sets how firmly the long run is tied to SSA; we use a moderate value here and probe its effect in Section 6. SSA publishes no clean 22-group table, so this is an age-graded profile anchored on the published bands.
# SSA Alt-2 ultimate annual % reduction in central death rates, age-graded (22 groups)
SSA_PCT = np.array([1.7, 1.6, 1.5, 1.4, 1.0, 0.95, 0.9, 0.9, 0.88, 0.85, 0.83,
0.8, 0.78, 0.76, 0.75, 0.74, 0.72, 0.68, 0.62, 0.55, 0.5, 0.45])
PSI0 = -SSA_PCT / 100.0 # steady-state prior mean (Delta log m)
PSI_SD = 0.005 # moderate prior sd
samp_mean = {s: Y[s].mean(axis=0) for s in Y} # historical mean Delta log m
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(GROUP_LO, 100 * -PSI0, "s-", color="k", lw=1.6, ms=3, label="SSA Alt-2 prior centre")
for s in ["female", "male"]:
ax.plot(GROUP_LO, 100 * -samp_mean[s], "o--", color=COLOR[s], ms=3, lw=1, label=f"{s} sample mean (1933-2000)")
ax.set(title="Ultimate annual mortality improvement (% / yr): SSA prior vs historical",
xlabel="Age (group start)", ylabel="% decline / yr")
ax.grid(alpha=.25); ax.legend(fontsize=8)
plt.show()
print("The SSA ultimate rates are slower than the historical pace at most ages -> "
"the informative prior tempers long-run improvement.")
The SSA ultimate rates are slower than the historical pace at most ages -> the informative prior tempers long-run improvement.
3. Estimation and convergence¶
NDRAW, BURN = 2500, 600
fits = {}
for s in ["female", "male"]:
t = time.time()
fits[s] = bvar.gibbs_steady_state(Y[s], p=1, lam1=0.2, lam3=1.0,
psi_prior_mean=PSI0, psi_prior_sd=PSI_SD,
ndraw=NDRAW, burn=BURN, seed=1)
print(f"{s}: {NDRAW} draws in {time.time()-t:.1f}s "
f"(Pi {fits[s]['Pi'].shape}, mu {fits[s]['mu'].shape}, Sigma {fits[s]['Sigma'].shape})")
female: 2500 draws in 26.3s (Pi (2500, 22, 22), mu (2500, 22), Sigma (2500, 22, 22))
male: 2500 draws in 27.2s (Pi (2500, 22, 22), mu (2500, 22), Sigma (2500, 22, 22))
# convergence traces: Sigma, an AR coefficient, and a steady-state component
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}$")
axes[1].plot(fs["Pi"][:, 6, 6], lw=.5); axes[1].set_title("$\Pi$ own lag, 25-29")
axes[2].plot(fs["mu"][:, 15], lw=.5); axes[2].set_title("$\Psi$ (steady state), 70-74")
for a in axes: a.set_xlabel("kept draw"); a.grid(alpha=.2)
fig.suptitle("Trace plots (female, steady-state BVAR p=1)", y=1.03); fig.tight_layout(); plt.show()
<>:4: SyntaxWarning: invalid escape sequence '\S'
<>:5: SyntaxWarning: invalid escape sequence '\P'
<>:6: SyntaxWarning: invalid escape sequence '\P'
<>:4: SyntaxWarning: invalid escape sequence '\S'
<>:5: SyntaxWarning: invalid escape sequence '\P'
<>:6: SyntaxWarning: invalid escape sequence '\P'
C:\Users\user\AppData\Local\Temp\ipykernel_12700\3783769293.py:4: SyntaxWarning: invalid escape sequence '\S'
axes[0].plot(fs["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$")
C:\Users\user\AppData\Local\Temp\ipykernel_12700\3783769293.py:5: SyntaxWarning: invalid escape sequence '\P'
axes[1].plot(fs["Pi"][:, 6, 6], lw=.5); axes[1].set_title("$\Pi$ own lag, 25-29")
C:\Users\user\AppData\Local\Temp\ipykernel_12700\3783769293.py:6: SyntaxWarning: invalid escape sequence '\P'
axes[2].plot(fs["mu"][:, 15], lw=.5); axes[2].set_title("$\Psi$ (steady state), 65-69")
3.1 The estimated steady state $\Psi$¶
The posterior of $\Psi$ (the ultimate improvement rate by age) blends the SSA prior with the historical data. Below, the posterior mean and 90% band against the SSA prior centre — where the data are informative the posterior pulls away from the prior; where the very oldest ages are noisy the posterior is pulled below the SSA prior (into faster implied improvement) by the neighbouring age groups.
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for ax, s in zip(axes, ["female", "male"]):
mu = fits[s]["mu"] # (ndraw, 22), units Delta log m
lo, mid, hi = np.percentile(100 * -mu, [5, 50, 95], axis=0) # as % decline / yr
ax.fill_between(GROUP_LO, lo, hi, color=COLOR[s], alpha=.15)
ax.plot(GROUP_LO, mid, "-", color=COLOR[s], lw=1.8, label="posterior")
ax.plot(GROUP_LO, 100 * -PSI0, "s--", color="k", ms=3, lw=1, label="SSA Alt-2 prior")
ax.set(title=f"{s}: steady-state improvement $\Psi$", xlabel="Age", ylabel="% decline / yr")
ax.grid(alpha=.25); ax.legend(fontsize=8)
fig.tight_layout(); plt.show()
<>:8: SyntaxWarning: invalid escape sequence '\P'
<>:8: SyntaxWarning: invalid escape sequence '\P'
C:\Users\user\AppData\Local\Temp\ipykernel_12700\850272627.py:8: SyntaxWarning: invalid escape sequence '\P'
ax.set(title=f"{s}: steady-state improvement $\Psi$", xlabel="Age", ylabel="% decline / yr")
4. Forecasts vs. actual and the Lee-Carter benchmark¶
Posterior predictive $e_0$ and $e_{65}$ (median and 90% band) vs actual and Lee-Carter (also with a 90% plug-in interval). Anchored at the observed 2000 log-rates.
def ex_bands(fit, s, age, thin=2):
draws = bvar.forecast(fit, Y[s], len(FYEARS))[::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.array([[life_table_abridged(np.exp(base + cum[d, h]), s)["ex"][idx]
for h in range(len(FYEARS))] for d in range(nd)])
return np.percentile(E, [5, 50, 95], axis=0)
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)
e = lambda kv: np.array([life_table_abridged(project_rates(fit, kv, jo)[i], s)["ex"][idx]
for i in range(len(FYEARS))])
return e(fc["k"] + 1.645*fc["se"]), e(fc["k"]), e(fc["k"] - 1.645*fc["se"])
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="LC 90%")
ax.plot(FYEARS, lmid, ":", color=COLOR[s], lw=1.6, label="LC 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("Steady-state BVAR forecasts vs actual vs Lee-Carter", y=1.0); fig.tight_layout(); plt.show()
4.1 Log-mortality forecasts at selected ages¶
The rate-level view (female), with BVAR and Lee-Carter 90% bands. Compared with the Minnesota model, the steady-state forecasts are tempered (pulled toward the slower SSA ultimate rates) and their bands are narrower — the informative prior buys confidence. This helps where improvement genuinely slowed (age 0, and the 25-29 reversal from the opioid epidemic) but can over-temper the oldest ages.
def logm_bands(fit, s):
dr = bvar.forecast(fit, Y[s], len(FYEARS))
lm = logm[s].loc[2000].values + np.cumsum(dr, axis=1)
return np.percentile(lm, [5, 50, 95], axis=0)
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"])
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="LC 90%")
ax.plot(FYEARS, lcmid[:, j], ":", color="#555", lw=1.5, label="LC 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("Steady-state BVAR log-mortality forecasts at selected ages", y=1.0)
fig.tight_layout(); plt.show()
5. Lag length — $p=1,2,3$¶
As with the Minnesota model, the shrinkage keeps higher-lag forecasts in check.
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, 1100, "#e07b39"), (3, 700, "#3a7d44")]:
fp = bvar.gibbs_steady_state(Y["female"], p=p, lam1=0.2, lam3=1.0,
psi_prior_mean=PSI0, psi_prior_sd=PSI_SD, ndraw=nd, burn=300, seed=2)
ax.plot(FYEARS, ex_bands(fp, "female", 0)[1], "--", color=c, lw=1.4, label=f"BVAR p={p}")
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="Female e0: steady-state BVAR at p=1,2,3", xlabel="Year", ylabel="e0 (years)")
ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()
6. Prior tightness — how firmly to trust SSA¶
The steady-state prior sd $\lambda_\Psi$ controls the informativeness. Tight = the long run is pinned to SSA; diffuse = the data (historical pace) dominate and the model approaches an uninformative steady-state VAR. Female $e_0$ forecast under a range of $\lambda_\Psi$.
fig, ax = plt.subplots(figsize=(9, 4))
act = ex_series(g["female"], "female", 0)
ax.plot(FYEARS, act.reindex(FYEARS), color=COLOR["female"], lw=2, label="actual")
for sd, c in [(0.002, "#08519c"), (0.005, "#3182bd"), (0.02, "#6baed6"), (10.0, "#bdbdbd")]:
fsd = bvar.gibbs_steady_state(Y["female"], p=1, lam1=0.2, lam3=1.0,
psi_prior_mean=PSI0, psi_prior_sd=sd, ndraw=1500, burn=300, seed=3)
lab = f"$\lambda_\Psi$={sd}" + (" (diffuse)" if sd >= 10 else "")
ax.plot(FYEARS, ex_bands(fsd, "female", 0)[1], "--", color=c, lw=1.4, label=lab)
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="Female e0: steady-state forecast vs prior tightness $\lambda_\Psi$",
xlabel="Year", ylabel="e0 (years)")
ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()
print("Tighter prior -> lower e0 (tied to SSA's slower improvement); diffuse -> higher (data-driven).")
<>:7: SyntaxWarning: invalid escape sequence '\l'
<>:10: SyntaxWarning: invalid escape sequence '\l'
<>:7: SyntaxWarning: invalid escape sequence '\l'
<>:10: SyntaxWarning: invalid escape sequence '\l'
C:\Users\user\AppData\Local\Temp\ipykernel_12700\709369282.py:7: SyntaxWarning: invalid escape sequence '\l'
lab = f"$\lambda_\Psi$={sd}" + (" (diffuse)" if sd >= 10 else "")
C:\Users\user\AppData\Local\Temp\ipykernel_12700\709369282.py:10: SyntaxWarning: invalid escape sequence '\l'
ax.set(title="Female e0: steady-state forecast vs prior tightness $\lambda_\Psi$",
Tighter prior -> lower e0 (tied to SSA's slower improvement); diffuse -> higher (data-driven).
7. Steady-state vs. Minnesota — the rate / life-expectancy trade-off¶
Out-of-sample forecast error (pre-COVID 2001-2019): mean absolute error in log m (all 22 groups) and in $e_0$ / $e_{65}$, for both BVAR models and Lee-Carter.
def err_table(s):
base = logm[s].loc[2000].values; actual_logm = np.log(g[s].loc[2001:2024]).values
a0 = ex_series(g[s], s, 0).reindex(FYEARS).values; a65 = ex_series(g[s], s, 65).reindex(FYEARS).values
n = 19
out = {}
# steady-state (reuse fits) and Minnesota
models = {"Steady-state": fits[s],
"Minnesota": bvar.gibbs_minnesota(Y[s], p=1, lam1=0.2, lam3=1.0, ndraw=2000, burn=500, seed=1)}
for name, fit in models.items():
lm = np.median(base + np.cumsum(bvar.forecast(fit, Y[s], len(FYEARS)), axis=1), axis=0)
e0 = np.array([life_table_abridged(np.exp(lm[h]), s)["ex"][0] for h in range(len(FYEARS))])
e65 = np.array([life_table_abridged(np.exp(lm[h]), s)["ex"][list(GROUP_LO).index(65)] for h in range(len(FYEARS))])
out[name] = (np.abs(lm[:n]-actual_logm[:n]).mean(), np.abs(e0[:n]-a0[:n]).mean(), np.abs(e65[:n]-a65[:n]).mean())
# Lee-Carter
lcfit = lee_carter(np.log(g[s].loc[SAMPLE[0]:SAMPLE[1]])); fc = forecast_lee_carter(lcfit, len(FYEARS))
lclm = np.log(project_rates(lcfit, fc["k"], base))
e0 = np.array([life_table_abridged(np.exp(lclm[h]), s)["ex"][0] for h in range(len(FYEARS))])
e65 = np.array([life_table_abridged(np.exp(lclm[h]), s)["ex"][list(GROUP_LO).index(65)] for h in range(len(FYEARS))])
out["Lee-Carter"] = (np.abs(lclm[:n]-actual_logm[:n]).mean(), np.abs(e0[:n]-a0[:n]).mean(), np.abs(e65[:n]-a65[:n]).mean())
return out
print(f"{'sex':>7} {'model':>13} | {'logm MAE':>9} {'e0 MAE':>7} {'e65 MAE':>8}")
for s in ["female", "male"]:
for name, (r, e0, e65) in err_table(s).items():
print(f"{s:>7} {name:>13} | {r:9.4f} {e0:7.2f} {e65:8.2f}")
sex model | logm MAE e0 MAE e65 MAE
female Steady-state | 0.0866 0.27 0.51 female Minnesota | 0.1288 0.24 0.32 female Lee-Carter | 0.1259 0.32 0.22
male Steady-state | 0.0995 0.55 0.97 male Minnesota | 0.1204 0.31 0.90 male Lee-Carter | 0.1178 0.32 0.92
8. Summary¶
- The steady-state BVAR reparametrizes the VAR around its long-run mean $\Psi$ (the ultimate mortality-improvement rate), letting us place an informative SSA Alt-2 prior directly on the long run — estimated with a 3-block Gibbs sampler.
- The informative prior tempers the forecasts and narrows the bands relative to the Minnesota model.
- Trade-off (Section 7): the steady-state model tends to forecast the age-specific rates more accurately, while the Minnesota model does better on the $e_0$ summary — because $e_0$ is dominated by old-age mortality, where the SSA prior over-tempers improvement. Which model is "better" depends on the target.
- The result is sensitive to the prior tightness $\lambda_\Psi$ (Section 6): loosening it moves the steady-state model toward the data-driven Minnesota result.