Forecast performance — mortality (expanding-window backtest)¶
A pseudo-real-time backtest of the mortality forecasts:
- fit to 1970, forecast 1971-1980; fit to 1980, forecast 1981-1990; ... through fit to 2020 — an expanding window with non-overlapping 10-year forecast blocks, so every year 1971-2024 is forecast once at a horizon of 1-10 years.
- Three models: Lee-Carter, Minnesota BVAR, steady-state BVAR (SSA Alt-2 prior), scored on $e_0$ and $e_{65}$ by sex, with 90% intervals.
Output: the realized series with the forecasts made along the way (and their bands), error by horizon, and accuracy + calibration scores. (The steady-state model uses a fixed SSA Alt-2 prior, so at early origins it injects a modern long-run assumption.)
Data: US series from the Human Mortality Database (mortality.org) and the Human Fertility Database (humanfertility.org).
1. Setup¶
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_LO,
lee_carter, forecast_lee_carter, project_rates)
import bvar
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
DIR = "."
ORIGINS, H, AGES = [1970, 1980, 1990, 2000, 2010, 2020], 10, [0, 65]
MODELS = ["Lee-Carter", "Minnesota BVAR", "Steady-state BVAR"]
MCOL = {"Lee-Carter": "#c1121f", "Minnesota BVAR": "#118ab2", "Steady-state BVAR": "#8338ec"}
NB, BB = 1200, 300
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"]}
ex_actual = {s: {a: ex_series(g[s], s, a) for a in AGES} for s in g}
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
IDX = {a: list(GROUP_LO).index(a) for a in AGES}
def forecast_ex(model, sex, origin, H=H):
"""Return {age: (lo, mid, hi)} 5/50/95 percentile e_x forecast paths (horizons 1..H),
for all AGES, from a single model fit."""
sample = np.log(g[sex].loc[1933:origin]); obs = np.log(g[sex].loc[origin].values)
if model == "Lee-Carter":
fit = lee_carter(sample); fc = forecast_lee_carter(fit, H)
def ex_paths(kv):
LT = [life_table_abridged(project_rates(fit, kv, obs)[h], sex)["ex"] for h in range(H)]
return {a: np.array([LT[h][IDX[a]] for h in range(H)]) for a in AGES}
hi, mid, lo = ex_paths(fc["k"] - 1.645*fc["se"]), ex_paths(fc["k"]), ex_paths(fc["k"] + 1.645*fc["se"])
return {a: (lo[a], mid[a], hi[a]) for a in AGES}
Yd = np.diff(sample.values, axis=0)
if model == "Minnesota BVAR":
fit = bvar.gibbs_minnesota(Yd, p=1, lam1=0.2, ndraw=NB, burn=BB, seed=1)
else:
fit = bvar.gibbs_steady_state(Yd, p=1, lam1=0.2, psi_prior_mean=PSI0,
psi_prior_sd=0.005, ndraw=NB, burn=BB, seed=1)
fit = bvar.filter_stationary(fit)
cum = np.cumsum(bvar.forecast(fit, Yd, H), axis=1)[::4] # thinned draws
# EX[draw, horizon, age]
EX = np.empty((len(cum), H, len(AGES)))
for d in range(len(cum)):
for h in range(H):
ex = life_table_abridged(np.exp(obs + cum[d, h]), sex)["ex"]
EX[d, h] = [ex[IDX[a]] for a in AGES]
out = {}
for ai, a in enumerate(AGES):
lo, mid, hi = np.percentile(EX[:, :, ai], [5, 50, 95], axis=0)
out[a] = (lo, mid, hi)
return out
print(f"origins {ORIGINS}; ages {AGES}; {len(ORIGINS)*2*len(MODELS)} model fits (each gives e0 & e65)")
origins [1970, 1980, 1990, 2000, 2010, 2020]; ages [0, 65]; 36 model fits (each gives e0 & e65)
2. Run the backtest¶
t0 = time.time(); rows = []
for origin in ORIGINS:
for sex in ["female", "male"]:
for model in MODELS:
res = forecast_ex(model, sex, origin)
for a in AGES:
lo, mid, hi = res[a]
for h in range(H):
yr = origin + h + 1
act = ex_actual[sex][a].get(yr, np.nan)
rows.append(dict(model=model, sex=sex, origin=origin, age=a, h=h + 1, year=yr,
ex_lo=lo[h], ex_fc=mid[h], ex_hi=hi[h], ex_act=act))
R = pd.DataFrame(rows)
R["err"] = R["ex_fc"] - R["ex_act"]
R["covered"] = (R["ex_act"] >= R["ex_lo"]) & (R["ex_act"] <= R["ex_hi"])
print(f"done in {time.time()-t0:.0f}s; {R['ex_act'].notna().sum()} scored forecasts "
f"({len(AGES)} ages x {len(MODELS)} models x 2 sexes)")
done in 373s; 648 scored forecasts (2 ages x 3 models x 2 sexes)
3. Realized series and the forecasts made along the way¶
The full realized series (black) with each model's 10-year forecasts drawn from the origin where they'd have been made — each carrying its 90% band (BVAR = posterior predictive; Lee-Carter = plug-in from $k_t$ uncertainty). The gap between the median and the black line is the forecast error; the band should mostly contain the black line if it is well-calibrated. One figure for $e_0$, one for $e_{65}$.
def plot_reconstruction(age):
fig, axes = plt.subplots(len(MODELS), 2, figsize=(14, 10), sharex=True)
for r, model in enumerate(MODELS):
for c, sex in enumerate(["female", "male"]):
ax = axes[r, c]; ea = ex_actual[sex][age]
ax.plot(ea.loc[1965:2024].index, ea.loc[1965:2024], "k-", lw=1.8, zorder=5, label="realized")
for origin in ORIGINS:
d = R[(R.model == model) & (R.sex == sex) & (R.origin == origin) & (R.age == age)]
ax.fill_between(d["year"], d["ex_lo"], d["ex_hi"], color=MCOL[model], alpha=.13,
label="90% band" if origin == ORIGINS[0] else None)
ax.plot(d["year"], d["ex_fc"], "-", color=MCOL[model], lw=1.4, alpha=.9,
label="forecast median" if origin == ORIGINS[0] else None)
ax.plot(d["year"].iloc[0], d["ex_fc"].iloc[0], "o", color=MCOL[model], ms=3.5)
for o in ORIGINS: ax.axvline(o, color="gray", lw=.3, ls=":")
ax.set_title(f"{model} — {sex}", fontsize=10); ax.grid(alpha=.25)
if c == 0: ax.set_ylabel(f"e{age} (years)")
axes[-1, 0].set_xlabel("Year"); axes[-1, 1].set_xlabel("Year")
axes[0, 0].legend(fontsize=8, loc="lower right")
fig.suptitle(f"Realized $e_{{{age}}}$ and the 10-year forecasts made along the way", y=1.0)
fig.tight_layout(); plt.show()
plot_reconstruction(0)
plot_reconstruction(65)
4. Error by forecast horizon¶
RMSE at each horizon 1-10 (pooled over origins and sexes), for $e_0$ and $e_{65}$.
fig, axes = plt.subplots(1, 2, figsize=(14, 4.5))
for ax, age in zip(axes, AGES):
by_h = R[R.age == age].dropna(subset=["err"]).groupby(["model", "h"])["err"].apply(
lambda e: np.sqrt((e**2).mean()))
for model in MODELS:
ax.plot(by_h[model].index, by_h[model].values, "o-", color=MCOL[model], lw=1.8, ms=4, label=model)
ax.set(title=f"e{age} RMSE by horizon", xlabel="years ahead", ylabel="RMSE (years)")
ax.grid(alpha=.25)
axes[0].legend(fontsize=8); fig.tight_layout(); plt.show()
5. Accuracy and calibration scores¶
D = R.dropna(subset=["err"])
tab = D.groupby(["age", "model"]).apply(lambda x: pd.Series({
"MAE": x["err"].abs().mean(), "RMSE": np.sqrt((x["err"]**2).mean()),
"bias": x["err"].mean(), "cover90": x["covered"].mean()})).round(3)
print("e0 and e65 forecast performance, 1971-2024 (all horizons, both sexes):")
print(tab.to_string())
print("\n bias>0 = over-predicting (too optimistic); cover90 = share of actuals in the "
"90% band (~0.90 = calibrated). COVID 2020-24 inflates the last block.")
e0 and e65 forecast performance, 1971-2024 (all horizons, both sexes):
MAE RMSE bias cover90
age model
0 Lee-Carter 0.559 0.772 -0.011 0.713
Minnesota BVAR 0.547 0.773 0.029 0.889
Steady-state BVAR 0.584 0.814 -0.263 0.824
65 Lee-Carter 0.357 0.490 -0.148 0.537
Minnesota BVAR 0.343 0.472 -0.105 0.870
Steady-state BVAR 0.337 0.472 -0.197 0.861
bias>0 = over-predicting (too optimistic); cover90 = share of actuals in the 90% band (~0.90 = calibrated). COVID 2020-24 inflates the last block.
6. Age profiles — 1- vs 10-year-ahead projections¶
A cross-sectional view: the realized log-mortality age profile in a target year (2019, pre-COVID) against the model's forecast of that same profile made 1 year ahead (fit to 2018) and 10 years ahead (fit to 2009), each with a 90% band. The 1-year forecast should sit almost on the realized profile with a tight band; the 10-year forecast is wider and can drift — most visibly at the ages that were moving fastest.
TARGET = 2019
def logm_profile(model, sex, origin, horizon):
sample = np.log(g[sex].loc[1933:origin]); obs = np.log(g[sex].loc[origin].values)
if model == "Lee-Carter":
fit = lee_carter(sample); fc = forecast_lee_carter(fit, horizon)
prof = lambda kv: np.log(project_rates(fit, kv, obs)[horizon - 1])
return prof(fc["k"] + 1.645*fc["se"]), prof(fc["k"]), prof(fc["k"] - 1.645*fc["se"]) # lo,mid,hi
Yd = np.diff(sample.values, axis=0)
if model == "Minnesota BVAR":
fit = bvar.gibbs_minnesota(Yd, p=1, lam1=0.2, ndraw=NB, burn=BB, seed=1)
else:
fit = bvar.gibbs_steady_state(Yd, p=1, lam1=0.2, psi_prior_mean=PSI0, psi_prior_sd=0.005,
ndraw=NB, burn=BB, seed=1)
fit = bvar.filter_stationary(fit)
logm = obs + np.cumsum(bvar.forecast(fit, Yd, horizon), axis=1)[:, horizon - 1, :]
return np.percentile(logm, [5, 50, 95], axis=0)
H1C, H10C = "#2a9d8f", "#e76f51"
fig, axes = plt.subplots(len(MODELS), 2, figsize=(13, 10), sharex=True)
for r, model in enumerate(MODELS):
for c, sex in enumerate(["female", "male"]):
ax = axes[r, c]
realized = np.log(g[sex].loc[TARGET].values)
l1, m1, h1 = logm_profile(model, sex, TARGET - 1, 1)
l10, m10, h10 = logm_profile(model, sex, TARGET - 10, 10)
ax.plot(GROUP_LO, realized, "k-", lw=2, label=f"realized {TARGET}", zorder=5)
ax.fill_between(GROUP_LO, l1, h1, color=H1C, alpha=.2)
ax.plot(GROUP_LO, m1, "-", color=H1C, lw=1.4, label="1-yr ahead")
ax.fill_between(GROUP_LO, l10, h10, color=H10C, alpha=.2)
ax.plot(GROUP_LO, m10, "--", color=H10C, lw=1.4, label="10-yr ahead")
ax.set_title(f"{model} — {sex}", fontsize=10); ax.grid(alpha=.25)
if c == 0: ax.set_ylabel("log m")
axes[-1, 0].set_xlabel("age"); axes[-1, 1].set_xlabel("age")
axes[0, 0].legend(fontsize=8, loc="upper left")
fig.suptitle(f"Log-mortality age profile {TARGET}: realized vs 1- and 10-year-ahead forecasts (90% bands)", y=1.0)
fig.tight_layout(); plt.show()
7. Notes and analysis¶
Point accuracy is similar; calibration is not. Over 1971-2024, Lee-Carter and the Minnesota BVAR are effectively tied on both $e_0$ (MAE ~0.55) and $e_{65}$ (~0.34), and the steady-state model is close (marginally worse on $e_0$, marginally best on $e_{65}$). The decisive difference is interval calibration: Lee-Carter's plug-in bands cover only ~71% of realized $e_0$ and just ~54% of $e_{65}$ against a nominal 90%, while the BVAR posterior bands are near nominal (~0.87). For anything needing a credible range — pension funding, solvency stress tests — Lee-Carter understates uncertainty by roughly half at age 65, even though its point forecast is just as good.
The errors are structural, not random. Every model under-forecasts old-age improvement ($e_{65}$ bias -0.11 to -0.20), and in the age-profile view (Section 6) the 10-year forecasts drift below realized mortality precisely at the young-adult/middle ages whose trend broke after ~2010 (the improvement slowdown; the male "deaths of despair"). These are coherent, age-patterned biases from regime changes that no extrapolative model anticipates — the same lesson as the sex-gap and fertility episodes.
The informative prior helps fertility but hurts mortality. The steady-state SSA-Alt-2 prior makes $e_0$ systematically pessimistic (bias -0.26): it caps improvement below the fast gains of the 1970s-90s. That is the mirror image of fertility, where anchoring the long run to $F^*$ reduces error — because there the diffuse model was badly biased (reverting to the baby-boom mean) and here it was not.
Caveats. The steady-state model uses a fixed (2024) SSA prior at every origin (a mild look-ahead at early origins); non-overlapping blocks give only ~6 origins; and the 2021-2024 block is four years carrying the COVID shock, penalising every model at the long horizons. Natural extensions: rate-level scores by age group and finer/rolling origins.