Population projection — cohort-component method

The capstone: project the US population forward by feeding the forecast mortality (survival) and forecast fertility (births) through the cohort-component identity (Meseguer 2010, eq. 22), with net migration treated as an exogenous assumption.

$$P_{t+1}(0)=\text{surviving births}+M(0),\qquad P_{t+1}(a{+}1)=P_t(a)\,e^{-m_x(a)}+M(a{+}1).$$

• Mortality — Minnesota BVAR (best point accuracy in the backtest), by sex, disaggregated from the 22 groups to single years of age.
• Fertility — steady-state BVAR ($F^*=2$, best in the backtest) → births, split by a sex ratio at birth of 1.05.
• Net migration — exogenous: a fixed age/sex schedule (recent 2010-19 average, ~1.1M/yr) held constant, not modelled. Its level is a scenario dial (Section 5).
• Jump-off: Jan-1 2025; projected to 2060, single year of age.
• Uncertainty bands propagate the BVAR mortality + fertility posterior draws through the projection (migration fixed).

Data: US series from the Human Mortality Database (mortality.org) and the Human Fertility Database (humanfertility.org).

1. Setup and jump-off

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, _gid as fert_gid_fn
from population import net_migration
import bvar, projection as pj
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
CF, CM = "#d1495b", "#00798c"

DIR = "."
JUMP, TARGET = 2025, 2060
H = TARGET - JUMP                          # annual steps (rate-years); H+1 population snapshots
NDRAW, NFAN = 1200, 200
ages = np.arange(111)
mort_gid = pj.group_of_age(ages, GROUP_LO)
fert_gid = np.array([fert_gid_fn(a) for a in ages])

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(JUMP)].values.astype(float) for s in ["female", "male"]}
RATE_YEARS = np.arange(JUMP, TARGET)          # annual rate/flow years (length H)
POP_YEARS = np.arange(JUMP, TARGET + 1)       # Jan-1 population snapshots (length H+1)
print(f"jump-off {JUMP}: {sum(P0[s].sum() for s in P0)/1e6:.1f}M; projecting to {TARGET}")

jump-off 2025: 340.4M; projecting to 2060

2. Forecast inputs

2.1 Mortality (Minnesota BVAR) — survival

mort_cum = {}                              # posterior draws of cumulative Delta log m, by sex
mx_base = {s: mx[s].loc[2024].values for s in ["female", "male"]}
for s in ["female", "male"]:
    gm = aggregate_to_groups(mx, deaths, pop, s)
    Y = np.diff(np.log(gm.loc[1933:2024].values), axis=0)
    fit = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=1, lam1=0.2, ndraw=NDRAW, burn=300, seed=1))
    mort_cum[s] = np.cumsum(bvar.forecast(fit, Y, H), axis=1)      # (nd, H, 22)

# implied e0 forecast (median + 90% band) as a check
fig, ax = plt.subplots(figsize=(8, 3.6))
for s, c in [("female", CF), ("male", CM)]:
    gm = aggregate_to_groups(mx, deaths, pop, s).loc[2024].values
    E = np.array([[life_table_abridged(gm * np.exp(mort_cum[s][d, h]), s)["e0"] for h in range(H)]
                  for d in range(0, len(mort_cum[s]), 6)])
    lo, mid, hi = np.percentile(E, [5, 50, 95], axis=0)
    ax.fill_between(RATE_YEARS, lo, hi, color=c, alpha=.15); ax.plot(RATE_YEARS, mid, color=c, lw=1.8, label=s)
ax.set(title="Forecast e0 feeding the projection", xlabel="Year", ylabel="e0"); ax.grid(alpha=.25); ax.legend()
plt.show()

2.2 Fertility (steady-state BVAR) — births

Yf = np.log(G.loc[1933:2024].values); psi0 = np.log(G.loc[2024].values * 2.0 / tfr[2024])
ffit = bvar.filter_stationary(bvar.gibbs_steady_state(Yf, p=2, lam1=0.05, lam2=0.5, lam3=1.0, own_mean=0.8,
                              psi_prior_mean=psi0, psi_prior_sd=0.05, ndraw=NDRAW, burn=400, seed=1))
asfr_g_draws = np.exp(bvar.forecast(ffit, Yf, H))                  # (nd, H, 7) group ASFR
asfr_g_obs = G.loc[2024].values
asfr_base = np.zeros(111); ab = asfr_by_age(f"{DIR}/USAasfrTR.txt").loc[2024]
for a in ab.index: asfr_base[a] = ab[a]

TFR = asfr_g_draws.sum(axis=2)
lo, mid, hi = np.percentile(TFR, [5, 50, 95], axis=0)
fig, ax = plt.subplots(figsize=(8, 3.4))
ax.fill_between(RATE_YEARS, lo, hi, color="#8338ec", alpha=.15); ax.plot(RATE_YEARS, mid, color="#8338ec", lw=1.8)
ax.axhline(2.0, color="gray", ls=":", lw=1)
ax.set(title="Forecast TFR feeding the projection", xlabel="Year", ylabel="TFR"); ax.grid(alpha=.25)
plt.show()

2.3 Net migration — exogenous assumption

Migration is not forecast; following Meseguer we supply it as an external assumption. The baseline is the age/sex schedule of net migration averaged over 2010-2019 (~1.1M/yr), held constant every projection year. Its age profile concentrates at young working ages.

mig = {}
for s in ["female", "male"]:
    m, _ = net_migration(pop, deaths, births, s)
    mig[s] = m.loc[2010:2019].mean().values          # (111,) constant assumption
fig, ax = plt.subplots(figsize=(8, 3.4))
for s, c in [("female", CF), ("male", CM)]:
    ax.plot(ages[:90], mig[s][:90] / 1e3, color=c, lw=1.5, label=s)
ax.set(title=f"Exogenous net-migration assumption (total {sum(mig[s].sum() for s in mig)/1e6:.2f}M/yr)",
       xlabel="Age", ylabel="net migrants (000s/yr)"); ax.grid(alpha=.25); ax.legend()
plt.show()

3. Projection (central + uncertainty)

def mx_single(cum_s, base_s):
    return pj.disaggregate(base_s, np.exp(cum_s), mort_gid)        # (H, 111)
def asfr_single(asfr_g):
    return pj.disaggregate(asfr_base, asfr_g / asfr_g_obs, fert_gid)

# central (median forecast)
mxfc_c = {s: mx_single(np.median(mort_cum[s], axis=0), mx_base[s]) for s in P0}
asfr_c = asfr_single(np.median(asfr_g_draws, axis=0))
central = pj.project(P0, mxfc_c, asfr_c, mig)

# uncertainty: propagate NFAN joint draws
nd = min(len(mort_cum["female"]), len(asfr_g_draws))
idx = np.linspace(0, nd - 1, NFAN).astype(int)
totals = np.empty((NFAN, H + 1)); old_share = np.empty((NFAN, H + 1)); pyr_end = {s: np.empty((NFAN, 111)) for s in P0}
t0 = time.time()
for j, i in enumerate(idx):
    mxfc = {s: mx_single(mort_cum[s][i], mx_base[s]) for s in P0}
    tr = pj.project(P0, mxfc, asfr_single(asfr_g_draws[i]), mig)
    tot = np.array([sum(tr[s][t].sum() for s in tr) for t in range(H + 1)])
    old = np.array([sum(tr[s][t][65:].sum() for s in tr) for t in range(H + 1)])
    totals[j] = tot; old_share[j] = old / tot
    for s in P0: pyr_end[s][j] = tr[s][-1]
print(f"{NFAN} projections in {time.time()-t0:.0f}s")

200 projections in 0s

4. Results

tot_c = np.array([sum(central[s][t].sum() for s in central) for t in range(H + 1)]) / 1e6
tlo, thi = np.percentile(totals, [5, 95], axis=0) / 1e6
os_c = np.array([sum(central[s][t][65:].sum() for s in central) for t in range(H + 1)]) /        np.array([sum(central[s][t].sum() for s in central) for t in range(H + 1)]) * 100
oslo, oshi = np.percentile(old_share, [5, 95], axis=0) * 100

fig, (a1, a2) = plt.subplots(1, 2, figsize=(13, 4.2))
a1.fill_between(POP_YEARS, tlo, thi, color="#4b4e6d", alpha=.15, label="90%")
a1.plot(POP_YEARS, tot_c, color="#4b4e6d", lw=2, label="central")
a1.set(title="US total population", xlabel="Year", ylabel="millions"); a1.grid(alpha=.25); a1.legend()
a2.fill_between(POP_YEARS, oslo, oshi, color="#e07b39", alpha=.15)
a2.plot(POP_YEARS, os_c, color="#e07b39", lw=2)
a2.set(title="Share aged 65+", xlabel="Year", ylabel="%"); a2.grid(alpha=.25)
plt.tight_layout(); plt.show()
print(f"Total pop {TARGET}: {tot_c[-1]:.0f}M (90% {tlo[-1]:.0f}-{thi[-1]:.0f}); "
      f"65+ share {os_c[-1]:.1f}% (from {os_c[0]:.1f}% in {JUMP})")

Total pop 2060: 384M (90% 365-405); 65+ share 24.3% (from 18.2% in 2025)

# pyramids: jump-off vs target
fig, ax = plt.subplots(figsize=(8, 7))
A = 100
for s, sign, c in [("male", -1, CM), ("female", 1, CF)]:
    ax.barh(ages[:A+1], sign * P0[s][:A+1] / 1e6, height=1, color=c, alpha=.28,
            label=f"{s} {JUMP}")
    med = np.median(pyr_end[s], axis=0)
    ax.step(sign * med[:A+1] / 1e6, ages[:A+1], where="mid", color=c, lw=1.6,
            label=f"{s} {TARGET}")
ax.axvline(0, color="k", lw=.6); ax.set_ylim(0, A); ax.set_ylabel("Age")
xt = np.arange(-2.5, 2.6, 1.25); ax.set_xticks(xt); ax.set_xticklabels([f"{abs(v):.1f}" for v in xt])
ax.set_xlabel("Population (millions)   (male left / female right)")
ax.set_title(f"US population pyramid: {JUMP} (bars) vs {TARGET} (outline)")
ax.legend(fontsize=8); ax.grid(alpha=.2, axis="x"); plt.show()

5. Migration scenarios (exogenous)

Because migration is an exogenous input, its level is a policy/assumption dial. Below, the central projection under zero, baseline (~1.1M/yr), and high (1.5x) net migration — holding the forecast mortality and fertility fixed.

fig, ax = plt.subplots(figsize=(9, 4.2))
for factor, name, c in [(0.0, "zero", "#bbb"), (1.0, "baseline", "#4b4e6d"), (1.5, "high (1.5x)", "#2a9d8f")]:
    mg = {s: mig[s] * factor for s in P0}
    tr = pj.project(P0, mxfc_c, asfr_c, mg)
    tt = np.array([sum(tr[s][t].sum() for s in tr) for t in range(H + 1)]) / 1e6
    ax.plot(POP_YEARS, tt, lw=2, color=c, label=f"{name}: {TARGET} = {tt[-1]:.0f}M")
ax.set(title="Total population under exogenous migration scenarios",
       xlabel="Year", ylabel="millions"); ax.grid(alpha=.25); ax.legend()
plt.show()

6. Validation — 2000-2024 backtest

Does the projection method reproduce known population? We jump off at Jan-1 2000 and project to 2024 under the three mortality models (Lee-Carter, Minnesota, steady-state), each fit to 1933-2000, with fertility (steady-state forecast) and exogenous migration (2000-19 average, held constant) common to all three — so the spread is the mortality model. Intervals come from each mortality model's forecast draws; the projected male and female totals are compared to actual.

JB, HB, NDB = 2000, 24, 80
BYEARS = np.arange(JB, JB + HB + 1)
rng = np.random.default_rng(0)
P0b = {s: pop[s].loc[str(JB)].values.astype(float) for s in ["female", "male"]}
migb = {s: net_migration(pop, deaths, births, s)[0].loc[2000:2019].mean().values for s in P0b}
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

# common fertility (steady-state) fit to 1933-2000; base = observed 2000 single-year ASFR
Yfb = np.log(G.loc[1933:JB].values); psib = np.log(G.loc[JB].values * 2.0 / tfr[JB])
ffb = bvar.filter_stationary(bvar.gibbs_steady_state(Yfb, p=2, lam1=.05, lam2=.5, lam3=1.0, own_mean=.8,
                             psi_prior_mean=psib, psi_prior_sd=.05, ndraw=800, burn=200, seed=1))
asfr_base_b = np.zeros(111); abb = asfr_by_age(f"{DIR}/USAasfrTR.txt").loc[JB]
for a in abb.index: asfr_base_b[a] = abb[a]
asfr_cb = pj.disaggregate(asfr_base_b, np.median(np.exp(bvar.forecast(ffb, Yfb, HB)), axis=0) / G.loc[JB].values, fert_gid)

def mort_change_draws(model, sex):
    gm = aggregate_to_groups(mx, deaths, pop, sex); Y = np.diff(np.log(gm.loc[1933:JB].values), axis=0)
    if model == "Lee-Carter":
        fit = lee_carter(np.log(gm.loc[1933:JB])); fc = forecast_lee_carter(fit, HB)
        jo = np.log(gm.loc[JB].values); g0 = np.log(gm.loc[JB].values)
        return [np.log(project_rates(fit, fc["k"] + np.cumsum(rng.normal(0, fc["sigma"], HB)), jo)) - g0[None, :]
                for _ in range(NDB)]
    if model == "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, HB), axis=1)
    return [cum[j] for j in range(min(NDB, len(cum)))]

def family_totals(model):
    cg = {s: mort_change_draws(model, s) for s in P0b}
    n = min(len(cg["female"]), len(cg["male"]))
    tot = {s: np.zeros((n, HB + 1)) for s in P0b}
    for j in range(n):
        mxfc = {s: pj.disaggregate(mx[s].loc[JB].values, np.exp(cg[s][j]), mort_gid) for s in P0b}
        tr = pj.project(P0b, mxfc, asfr_cb, migb)
        for s in P0b: tot[s][j] = tr[s].sum(axis=1)
    return tot

bt = {m: family_totals(m) for m in ["Lee-Carter", "Minnesota BVAR", "Steady-state BVAR"]}
actb = {s: np.array([pop[s].loc[str(y)].sum() for y in BYEARS]) for s in P0b}
print("2024 projected vs actual (millions):")
for s in P0b:
    print(f"  {s}: actual {actb[s][-1]/1e6:.1f}  |  " +
          "  ".join(f"{m.split()[0]} {np.median(bt[m][s],0)[-1]/1e6:.1f}" for m in bt))

2024 projected vs actual (millions):
  female: actual 170.9  |  Lee-Carter 176.3  Minnesota 175.9  Steady-state 175.2
  male: actual 167.5  |  Lee-Carter 170.5  Minnesota 170.4  Steady-state 169.9

BCOL = {"Lee-Carter": "#c1121f", "Minnesota BVAR": "#118ab2", "Steady-state BVAR": "#8338ec"}
fig, axes = plt.subplots(1, 2, figsize=(13, 4.4))
for ax, s in zip(axes, ["female", "male"]):
    ax.plot(BYEARS, actb[s] / 1e6, "k-", lw=2.2, label="actual", zorder=5)
    for m, c in BCOL.items():
        q = np.percentile(bt[m][s], [5, 50, 95], axis=0) / 1e6
        ax.fill_between(BYEARS, q[0], q[2], color=c, alpha=.12)
        ax.plot(BYEARS, q[1], "--", color=c, lw=1.4, label=m)
    ax.axvline(2020, color="gray", lw=.5, ls=":")
    ax.set_title(f"{s.capitalize()} population, 2000-2024"); ax.set_xlabel("Year")
    ax.set_ylabel("millions"); ax.grid(alpha=.25)
axes[0].legend(fontsize=8)
fig.suptitle("Population projection backtest: actual vs 3 mortality models (from Jan-1 2000)", y=1.02)
fig.tight_layout(); plt.show()
print("Projections track actual closely through ~2019, then over-project: COVID (2020-22 "
      "excess deaths + migration collapse) is unforecastable, and migration is held at its "
      "pre-pandemic average. The three mortality models are nearly indistinguishable at the "
      "population-total level -- mortality-model choice barely moves the 24-year total.")

Projections track actual closely through ~2019, then over-project: COVID (2020-22 excess deaths + migration collapse) is unforecastable, and migration is held at its pre-pandemic average. The three mortality models are nearly indistinguishable at the population-total level -- mortality-model choice barely moves the 24-year total.

7. Notes

• This closes the loop: the BVAR mortality and fertility forecasts drive survival and births, and the cohort-component identity ages the population forward, with net migration exogenous (a fixed age/sex schedule, dialled in Section 5).
• The projection reproduces the expected US trajectory — total population still growing (driven largely by migration; Section 5 shows it flattens without it) and rapid population aging (65+ share rising sharply) as the boom cohorts reach old age.
• Approximations: group->single-year mortality/fertility via jump-off-scaled factors; annual constant-hazard survival; SRB fixed at 1.05. Bands reflect mortality+fertility forecast uncertainty only (migration is a fixed assumption, varied by scenario).