Population backtest — full-uncertainty model comparison (2001-2024)¶
Each model family forecasts both its own mortality and fertility (fit to 1933-2000), and all posterior uncertainty propagates through the cohort-component projection to the population. We jump off at Jan-1 2000, project to 2024, and compare the predicted male/female population to actual.
- Benchmark = Lee-Carter mortality + Lee-Tuljapurkar fertility ($F^*=2$).
- Minnesota BVAR = Minnesota mortality + Minnesota fertility.
- Steady-state BVAR = steady-state mortality + steady-state fertility ($F^*=2$).
Net migration is exogenous (2000-19 average, held constant, common to all three). Every joint mortality+fertility draw is projected, so the bands reflect the combined mortality and fertility forecast uncertainty.
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, read_births, asfr_by_age
from mortality import (aggregate_to_groups, GROUP_LO, life_table_abridged,
lee_carter, forecast_lee_carter, project_rates)
from fertility import (group_asfr, tfr_from_groups, constrained_arma_draws, reconstruct_group_rates)
from population import net_migration
import bvar, projection as pj
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
DIR = "."
JB, H, ND = 2000, 24, 100
BYEARS = np.arange(JB, JB + H + 1)
FAMILIES = ["Benchmark", "Minnesota BVAR", "Steady-state BVAR"]
FCOL = {"Benchmark": "#c1121f", "Minnesota BVAR": "#118ab2", "Steady-state BVAR": "#8338ec"}
BANDS = [("0-19", 0, 20), ("20-64", 20, 65), ("65-84", 65, 85), ("85+", 85, 111)]
rng = np.random.default_rng(0)
ages = np.arange(111); mort_gid = pj.group_of_age(ages, GROUP_LO)
from fertility import _gid as _fg
fert_gid = np.array([_fg(a) for a in ages])
PSI0_MORT = -np.array([1.7,1.6,1.5,1.4,1.0,.95,.9,.9,.88,.85,.83,.8,.78,.76,.75,.74,.72,.68,.62,.55,.5,.45]) / 100
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"); births = read_births(f"{DIR}/US Births.txt")
G = group_asfr(f"{DIR}/USAasfrTR.txt"); tfr = tfr_from_groups(G)
P0 = {s: pop[s].loc[str(JB)].values.astype(float) for s in ["female", "male"]}
mig = {s: net_migration(pop, deaths, births, s)[0].loc[2000:2019].mean().values for s in P0}
mx_base = {s: mx[s].loc[JB].values for s in P0}
G0 = G.loc[JB].values
asfr_base = np.zeros(111); ab = asfr_by_age(f"{DIR}/USAasfrTR.txt").loc[JB]
for a in ab.index: asfr_base[a] = ab[a]
actb = {s: np.array([pop[s].loc[str(y)].sum() for y in BYEARS]) for s in P0}
print(f"jump-off {JB}: {sum(P0[s].sum() for s in P0)/1e6:.1f}M; comparing to actual through {BYEARS[-1]}")
jump-off 2000: 280.1M; comparing to actual through 2024
2. Forecast inputs — each family's mortality and fertility (with uncertainty)¶
For every family we draw the mortality group-changes (by sex) and the fertility group rates, fit to 1933-2000. Below, the implied forecast $e_0$ and TFR that feed the projection.
def mort_change(fam, sex):
gm = aggregate_to_groups(mx, deaths, pop, sex); Ys = np.log(gm.loc[1933:JB])
if fam == "Benchmark":
fit = lee_carter(Ys); fc = forecast_lee_carter(fit, H); jo = np.log(gm.loc[JB].values); g0 = np.log(gm.loc[JB].values)
return np.array([np.log(project_rates(fit, fc["k"] + np.cumsum(rng.normal(0, fc["sigma"], H)), jo)) - g0[None, :]
for _ in range(ND)])
Y = np.diff(Ys.values, axis=0)
if fam == "Minnesota BVAR":
fit = bvar.gibbs_minnesota(Y, p=1, lam1=.2, ndraw=1000, burn=250, seed=1)
else:
fit = bvar.gibbs_steady_state(Y, p=1, lam1=.2, psi_prior_mean=PSI0_MORT, psi_prior_sd=.005, ndraw=1000, burn=250, seed=1)
cum = np.cumsum(bvar.forecast(bvar.filter_stationary(fit), Y, H), axis=1)
return cum[np.linspace(0, len(cum) - 1, ND).astype(int)]
def fert_asfr(fam):
if fam == "Benchmark":
b = lee_carter(np.log(G.loc[1933:JB]))["b"].values
td = constrained_arma_draws(tfr.loc[1933:JB], 2.0, H, ND, seed=1)
return np.array([reconstruct_group_rates(G0, b, td[d]) for d in range(ND)])
Yf = np.log(G.loc[1933:JB].values)
if fam == "Minnesota BVAR":
fit = bvar.gibbs_minnesota(Yf, p=2, lam1=.05, lam2=.5, lam3=1.0, own_mean=.8, ndraw=1500, burn=400, seed=1)
else:
psi = np.log(G0 * 2.0 / tfr[JB])
fit = bvar.gibbs_steady_state(Yf, p=2, lam1=.05, lam2=.5, lam3=1.0, own_mean=.8,
psi_prior_mean=psi, psi_prior_sd=.05, ndraw=1500, burn=400, seed=1)
d = np.exp(bvar.forecast(bvar.filter_stationary(fit), Yf, H))
return d[np.linspace(0, len(d) - 1, ND).astype(int)]
t0 = time.time()
MC = {fam: {s: mort_change(fam, s) for s in P0} for fam in FAMILIES} # (ND, H, 22)
FA = {fam: fert_asfr(fam) for fam in FAMILIES} # (ND, H, 7)
print(f"forecast draws in {time.time()-t0:.0f}s")
forecast draws in 60s
fig, axes = plt.subplots(1, 3, figsize=(15, 3.8))
for fam in FAMILIES:
for s, ls in [("female", "-"), ("male", "--")]:
gm0 = aggregate_to_groups(mx, deaths, pop, s).loc[JB].values
E = np.array([[life_table_abridged(gm0 * np.exp(MC[fam][s][d, h]), s)["e0"] for h in range(H)]
for d in range(0, ND, 4)])
axes[0 if s == "female" else 1].plot(BYEARS[1:], np.median(E, 0), ls, color=FCOL[fam], lw=1.5,
label=fam if s == "female" else None)
T = FA[fam].sum(axis=2)
lo, mid, hi = np.percentile(T, [5, 50, 95], axis=0)
axes[2].fill_between(BYEARS[1:], lo, hi, color=FCOL[fam], alpha=.12)
axes[2].plot(BYEARS[1:], mid, color=FCOL[fam], lw=1.6, label=fam)
axes[0].set_title("forecast e0 — female"); axes[1].set_title("forecast e0 — male")
axes[2].set_title("forecast TFR"); axes[2].axhline(2, color="gray", ls=":", lw=1)
for a in axes: a.grid(alpha=.25); a.set_xlabel("Year")
axes[0].legend(fontsize=7.5); axes[2].legend(fontsize=7.5)
fig.tight_layout(); plt.show()
3. Project every joint draw¶
def project_family(fam):
tot = {s: np.zeros((ND, H + 1)) for s in P0}
pyr = {s: np.zeros((ND, 111)) for s in P0}
band = {name: np.zeros((ND, H + 1)) for name, _, _ in BANDS}
for j in range(ND):
mxfc = {s: pj.disaggregate(mx_base[s], np.exp(MC[fam][s][j]), mort_gid) for s in P0}
asfrfc = pj.disaggregate(asfr_base, FA[fam][j] / G0, fert_gid)
tr = pj.project(P0, mxfc, asfrfc, mig)
both = tr["female"] + tr["male"] # total by age, (H+1, 111)
for s in P0:
tot[s][j] = tr[s].sum(axis=1); pyr[s][j] = tr[s][-1]
for name, lo, hi in BANDS:
band[name][j] = both[:, lo:hi].sum(axis=1)
return tot, pyr, band
t0 = time.time()
RES = {}; PYR = {}; BAND = {}
for fam in FAMILIES:
RES[fam], PYR[fam], BAND[fam] = project_family(fam)
print(f"projected {ND} draws x {len(FAMILIES)} families in {time.time()-t0:.0f}s")
projected 100 draws x 3 families in 0s
4. Predicted vs actual population, 2001-2024¶
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4))
panels = [("female", "Female"), ("male", "Male"), ("total", "Total")]
for ax, (key, lab) in zip(axes, panels):
act = actb["female"] + actb["male"] if key == "total" else actb[key]
ax.plot(BYEARS, act / 1e6, "k-", lw=2.4, label="actual", zorder=6)
for fam in FAMILIES:
R = RES[fam]["female"] + RES[fam]["male"] if key == "total" else RES[fam][key]
q = np.percentile(R, [5, 50, 95], axis=0) / 1e6
ax.fill_between(BYEARS, q[0], q[2], color=FCOL[fam], alpha=.13)
ax.plot(BYEARS, q[1], "--", color=FCOL[fam], lw=1.4, label=fam)
ax.axvline(2020, color="gray", lw=.5, ls=":")
ax.set_title(f"{lab} population"); ax.set_xlabel("Year"); ax.grid(alpha=.25)
axes[0].set_ylabel("millions"); axes[0].legend(fontsize=8)
fig.suptitle("Full-uncertainty projection (mortality + fertility) vs actual, from Jan-1 2000", y=1.02)
fig.tight_layout(); plt.show()
print("2024 total population (millions): actual "
f"{(actb['female'][-1]+actb['male'][-1])/1e6:.1f}")
for fam in FAMILIES:
R = RES[fam]["female"] + RES[fam]["male"]; q = np.percentile(R, [5, 50, 95], axis=0) / 1e6
inside = (act[-1] >= np.percentile(R, 5, 0)[-1]) & (act[-1] <= np.percentile(R, 95, 0)[-1])
print(f" {fam:18}: {q[1,-1]:.1f} (90% {q[0,-1]:.1f}-{q[2,-1]:.1f}) actual in band: {bool(inside)}")
2024 total population (millions): actual 338.5 Benchmark : 343.6 (90% 323.1-380.0) actual in band: True Minnesota BVAR : 357.6 (90% 333.3-377.1) actual in band: True Steady-state BVAR : 345.2 (90% 333.1-362.3) actual in band: True
5. Population by age group¶
The total hides the age structure — and splitting it out localizes exactly where the models agree and where they don't. Below, actual vs. the three families (full- uncertainty 90% bands) for four broad bands.
What to look for, and why:
- 20-64 (working age): near-perfect for all three families. Everyone aged 20-64 in 2024 was born before 2005 — mostly before the 2000 jump-off — so their count is just the year-2000 cohorts surviving forward plus migration, with no fertility involved. The tight bands and near-exact fit (196.9 / 197.0 / 196.0M vs actual 196.1M) confirm the survival + migration machinery is sound; this is the projection's most reliable band.
- 65-84 and 85+ (old): tight and accurate, with the families indistinguishable — these are already-born cohorts, so purely mortality-driven, and the bands comfortably contain the actual (the COVID dip is visible at 85+).
- 0-19 (children): the fertility battleground. This is the only band where the families separate, and where the bands are wide. Minnesota over-projects by ~16M (98 vs 82M actual) — its diffuse fertility reverts to the baby-boom-inflated mean and produces far too many births — while the $F^*=2$ families (benchmark, steady-state) sit within a few million.
Bottom line: the entire ~19M total-population overshoot of the Minnesota family is essentially all in the 0-19 band — a pure fertility-assumption artifact. At the population level the mortality models are interchangeable; fertility is what separates the projections, and it does so entirely through the young ages.
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
for ax, (name, lo, hi) in zip(axes.ravel(), BANDS):
act = np.array([sum(pop[s].loc[str(y)].values[lo:hi].sum() for s in P0) for y in BYEARS])
ax.plot(BYEARS, act / 1e6, "k-", lw=2.2, label="actual", zorder=6)
for fam in FAMILIES:
q = np.percentile(BAND[fam][name], [5, 50, 95], axis=0) / 1e6
ax.fill_between(BYEARS, q[0], q[2], color=FCOL[fam], alpha=.12)
ax.plot(BYEARS, q[1], "--", color=FCOL[fam], lw=1.3, label=fam)
ax.axvline(2020, color="gray", lw=.5, ls=":")
ax.set_title(f"Age {name}"); ax.set_ylabel("millions"); ax.grid(alpha=.25)
axes[0, 0].legend(fontsize=8)
axes[1, 0].set_xlabel("Year"); axes[1, 1].set_xlabel("Year")
fig.suptitle("Population by age group: full-uncertainty projection vs actual", y=1.0)
fig.tight_layout(); plt.show()
print("2024 by age group (millions): actual | benchmark / Minnesota / steady-state (medians)")
for name, lo, hi in BANDS:
a = sum(pop[s].loc['2024'].values[lo:hi].sum() for s in P0) / 1e6
med = [np.median(BAND[f][name], 0)[-1] / 1e6 for f in FAMILIES]
print(f" {name:>6}: {a:5.1f} | " + " / ".join(f"{m:.1f}" for m in med))
2024 by age group (millions): actual | benchmark / Minnesota / steady-state (medians)
0-19: 82.1 | 86.5 / 98.3 / 88.2
20-64: 196.1 | 196.9 / 197.0 / 196.0
65-84: 53.9 | 55.4 / 55.4 / 54.8
85+: 6.4 | 6.4 / 6.2 / 6.2
6. Notes¶
- Every family now forecasts both vital rates with full posterior uncertainty, propagated jointly through the projection — so the bands are wider and more honest than the mortality-only backtest.
- All families track actual closely to ~2017, then over-project: COVID (2020-22 excess deaths, unforecastable), the exogenous migration held at its pre-pandemic average, and fertility forecasts that revert above the realised post-2007 decline all push the prediction high.
- The families differ mainly through fertility (mortality is near-identical at the population level): the Minnesota family, whose diffuse fertility reverts to the baby-boom-inflated mean, predicts the most births and the highest population; the $F^*=2$ families (benchmark, steady-state) are lower and closer to actual.
Final note — why the steady-state BVAR¶
Across this whole project the steady-state BVAR stands out as the most complete methodology, for three reasons the analysis makes concrete:
Shrinkage lets a large system be modelled simultaneously. The Minnesota prior makes it feasible to estimate one joint VAR over all age groups at once — 22 mortality groups (484 autoregressive coefficients) or 7 fertility groups — on only ~67 years of data, where an unrestricted VAR would be hopelessly over-parametrised. This captures the full cross-age dynamic and correlation structure (the "spatial trick"), rather than collapsing it to a single Lee-Carter factor.
It accommodates prior information in a principled, non-ad-hoc way. The mean-adjusted (Villani) parametrisation places a genuine prior directly on the long-run mean — the ultimate mortality-improvement rate or the ultimate TFR (e.g. the SSA Alternative-2 assumptions). Its influence is governed by an explicit prior variance, and the data can override it. This replaces the ad-hoc constraints the benchmark needs (forcing $k_t$ / the TFR to a target via a hand-imposed ARMA) with a coherent probability statement whose strength you can dial.
It delivers full uncertainty. The Gibbs sampler yields the entire posterior — parameter and shock uncertainty — which propagates coherently all the way through the cohort-component projection to the population (as the bands above show). The plug-in Lee-Carter / Lee-Tuljapurkar intervals capture only part of this and were badly under-calibrated in the backtests (e.g. Lee-Carter's 90% band covered just ~54% of realised $e_{65}$, versus ~87% for the BVAR). Honest, calibrated intervals are exactly what long-range demographic and actuarial applications require.