Mortality forecasting — Lee-Carter benchmark¶
Forecast US all-cause mortality following Meseguer (2010), as groundwork for the Bayesian VAR. This notebook establishes the benchmark: the Lee-Carter model (and its coherent Li-Lee extension). Design choices (matching Meseguer):
- Age structure: 22 abridged groups — 0, 1-4, 5-9, ..., 95-99, 100+.
- Sexes modelled separately (Meseguer: "separate models for male and female mortality are estimated").
- Sample: 1933-2000. Out-of-sample forecast: 2001-2024 (24 held-out years).
Lee-Carter: $\log m_x(t) = a_x + b_x\,k_t + \varepsilon$, fit by SVD, with $k_t$ extrapolated as a random walk with drift. The Li-Lee variant adds a shared factor so the male/female gap stays coherent (bounded) in the long run.
Data: US series from the Human Mortality Database (mortality.org) and the Human Fertility Database (humanfertility.org).
1. Data — grouped mortality by sex¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from hmd_io import read_mx, read_deaths, read_population, read_e0, read_ltper
from mortality import (aggregate_to_groups, life_table_abridged, e0_series, ex_series,
GROUP_LABEL, GROUP_LO,
lee_carter, forecast_lee_carter, project_rates,
lee_carter_common, forecast_common, rates_from_common)
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
DIR = "."
SAMPLE = (1933, 2000)
FYEARS = np.arange(2001, 2025) # out-of-sample forecast years
SEXCOL = {"female": "Female", "male": "Male"}
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")
e0pub = read_e0(f"{DIR}/US E0per.txt")
# 22-group central death rates, by sex
g = {s: aggregate_to_groups(mx, deaths, pop, s) for s in ["female", "male"]}
print(f"{len(GROUP_LABEL)} age groups:", ", ".join(GROUP_LABEL))
print(f"years {g['female'].index.min()}-{g['female'].index.max()}; "
f"sample {SAMPLE[0]}-{SAMPLE[1]}, forecast {FYEARS[0]}-{FYEARS[-1]}")
g["female"].loc[[1933, 2000, 2024], ["0", "1-4", "25-29", "65-69", "100+"]].round(5)
22 age groups: 0, 1-4, 5-9, 10-14, 15-19, 20-24, 25-29, 30-34, 35-39, 40-44, 45-49, 50-54, 55-59, 60-64, 65-69, 70-74, 75-79, 80-84, 85-89, 90-94, 95-99, 100+ years 1933-2024; sample 1933-2000, forecast 2001-2024
| group | 0 | 1-4 | 25-29 | 65-69 | 100+ |
|---|---|---|---|---|---|
| Year | |||||
| 1933 | 0.05418 | 0.00439 | 0.00360 | 0.03556 | 0.35949 |
| 2000 | 0.00650 | 0.00029 | 0.00053 | 0.01515 | 0.45393 |
| 2024 | 0.00516 | 0.00023 | 0.00059 | 0.01167 | 0.41707 |
1.1 Grouping preserves life expectancy¶
Before modelling, confirm the abridged (22-group) life table reproduces the published single-year e0 — so the grouping choice doesn't distort our target.
fig, ax = plt.subplots(figsize=(9, 4))
for s in ["female", "male"]:
e0g = e0_series(g[s], s)
ax.plot(e0g.index, e0g, color=COLOR[s], lw=1.6, label=f"{s} — abridged (derived)")
ax.plot(e0pub.index, e0pub[SEXCOL[s]], ":", color="k", lw=1,
label="published" if s == "female" else None)
err = (e0g - e0pub[SEXCOL[s]]).abs()
print(f"{s}: abridged vs published e0 — mean abs err {err.mean():.3f} yr, max {err.max():.3f}")
ax.axvspan(SAMPLE[0], SAMPLE[1], color="gray", alpha=.08, label="sample")
ax.set(title="e0 from 22-group abridged life table vs published single-year",
xlabel="Year", ylabel="e0 (years)"); ax.grid(alpha=.25); ax.legend(fontsize=7.5)
plt.show()
female: abridged vs published e0 — mean abs err 0.017 yr, max 0.069 male: abridged vs published e0 — mean abs err 0.009 yr, max 0.063
2. The Lee-Carter model¶
Lee & Carter (1992) compress an entire age schedule of log death rates into a single time-varying index:
$$\log m_x(t) = a_x + b_x\,k_t + \varepsilon_{x,t}$$
- $a_x$ — average age profile: the mean of $\log m_x$ over time, i.e. the baseline shape of mortality across age.
- $k_t$ — mortality index: one number for the overall level of mortality in year $t$. As $k_t$ falls, mortality falls; it typically trends down over time.
- $b_x$ — age sensitivity: how strongly each age responds to a change in $k_t$. Large $b_x$ = fast improvement at that age; $b_x\approx0$ = little change.
Identification. The decomposition is only unique up to scale and shift, so two constraints are imposed: $\sum_x b_x = 1$ and $\sum_t k_t = 0$. With these, $a_x$ is exactly the time-mean and $k_t$ is centred.
Estimation (SVD). Subtract $a_x$ from the log-rates and take the singular value
decomposition of the centred $T\times A$ matrix; the leading singular vectors give
$b_x$ and $k_t$ (a rank-1 approximation). → lee_carter in mortality.py.
Forecasting. The age effects $a_x, b_x$ are held fixed and only the index $k_t$ is projected. $k_t$ is close to linear, so we model it as a random walk with drift, $k_{t+h}=k_t + h\,\delta + \text{errors}$, with drift $\delta=(k_T-k_1)/(T-1)$; the forecast standard error grows like $\sqrt{h}$ (the 95% fan). Projected rates are $\exp(a_x + b_x k_{t+h})$.
2.1 Jump-off adjustment (Lee-Miller 2001)¶
Because it is a rank-1 fit, the model's fitted rates in the last sample year do not exactly match the observed ones — for US males the fitted $e_0(2000)$ is 0.78 yr below actual (for females ~0). Launching the forecast from the fitted value makes it start disconnected from the data. We instead anchor to the observed last-year rates:
$$\log m_x(2000{+}h) = \log m_x^{\text{obs}}(2000) + b_x\,h\,\delta$$
so $h=0$ reproduces the data exactly and only the trend is projected. This removes
the discontinuity and, for males, cuts the out-of-sample $e_0$ error from ~1.0 to
~0.3 yr. → project_rates(..., logm_jumpoff=...).
2.2 Fitted parameters, sample 1933-2000¶
fits = {s: lee_carter(np.log(g[s].loc[SAMPLE[0]:SAMPLE[1]])) for s in ["female", "male"]}
fig, (axb, axk) = plt.subplots(1, 2, figsize=(12, 4))
for s in ["female", "male"]:
axb.plot(GROUP_LO, fits[s]["b"].to_numpy(), "o-", color=COLOR[s], ms=3, label=s)
axk.plot(fits[s]["k"].index, fits[s]["k"].to_numpy(), color=COLOR[s], lw=1.8, label=s)
axb.axhline(0, color="k", lw=.5)
axb.set(title="b_x — age pattern of mortality improvement", xlabel="Age (group start)",
ylabel="b_x (sums to 1)"); axb.grid(alpha=.25); axb.legend()
axk.set(title="k_t — mortality index (fit)", xlabel="Year", ylabel="k_t")
axk.grid(alpha=.25); axk.legend()
plt.show()
print("b_x>0 at nearly all ages means mortality falls almost everywhere as k_t declines; "
"it turns slightly negative only at the oldest ages; "
"larger b_x = faster improvement at that age.")
b_x>0 at nearly all ages means mortality falls almost everywhere as k_t declines; it turns slightly negative only at the oldest ages; larger b_x = faster improvement at that age.
3. Forecast and out-of-sample evaluation¶
$k_t$ is extrapolated as a random walk with drift; the 95% fan comes from the random-walk forecast standard error, and forecasts use the jump-off adjustment (Section 2.1) so they start from the observed 2000 rates. We convert each forecast rate vector to $e_0$ via the abridged life table and compare to the actual 2001-2024 values.
jumpoff = {s: np.log(g[s].loc[SAMPLE[1]].to_numpy()) for s in ["female", "male"]} # observed log-rates, 2000
def e0_from_k(fit, kvals, sex, jo):
R = project_rates(fit, np.atleast_1d(kvals), logm_jumpoff=jo)
return np.array([life_table_abridged(R[i], sex)["e0"] for i in range(R.shape[0])])
fig, axes = plt.subplots(1, 2, figsize=(12, 4.3), sharey=True)
oos = {}
for ax, s in zip(axes, ["female", "male"]):
fit = fits[s]; fc = forecast_lee_carter(fit, len(FYEARS)); jo = jumpoff[s]
e0_mid = e0_from_k(fit, fc["k"], s, jo)
e0_hi = e0_from_k(fit, fc["k"] - 1.96 * fc["se"], s, jo) # low k -> low mortality -> high e0
e0_lo = e0_from_k(fit, fc["k"] + 1.96 * fc["se"], s, jo)
actual = e0_series(g[s], s)
oos[s] = pd.DataFrame({"forecast": e0_mid, "actual": actual.reindex(FYEARS).to_numpy()}, index=FYEARS)
hist = actual.loc[SAMPLE[0]:SAMPLE[1]]
ax.plot(hist.index, hist, color="gray", lw=1.2, label="sample (fit)")
ax.plot(FYEARS, actual.reindex(FYEARS), color=COLOR[s], lw=2, label="actual (out-of-sample)")
ax.plot(FYEARS, e0_mid, "--", color="k", lw=1.3, label="LC forecast")
ax.fill_between(FYEARS, e0_lo, e0_hi, color="k", alpha=.12, label="95% interval")
ax.axvline(2000, color="k", lw=.5, ls=":")
ax.annotate("COVID", (2021, actual.get(2021, np.nan)), fontsize=8, color="#a00",
xytext=(2008, actual.get(2021, 76)))
ax.set(title=f"{s.capitalize()} e0: Lee-Carter forecast vs actual", xlabel="Year")
ax.grid(alpha=.25); ax.legend(fontsize=7.5)
axes[0].set_ylabel("e0 (years)")
plt.show()
# Out-of-sample e0 error, pre-COVID (2001-2019) vs full (2001-2024)
print("Lee-Carter out-of-sample e0 error (mean abs, years):")
for s in ["female", "male"]:
d = oos[s]; e = (d["forecast"] - d["actual"]).abs()
print(f" {s}: 2001-2019 {e.loc[2001:2019].mean():.2f} | 2001-2024 {e.mean():.2f} "
f"(inflated by the 2020-22 COVID shock no pre-2000 model can foresee)")
Lee-Carter out-of-sample e0 error (mean abs, years): female: 2001-2019 0.32 | 2001-2024 0.73 (inflated by the 2020-22 COVID shock no pre-2000 model can foresee) male: 2001-2019 0.32 | 2001-2024 0.65 (inflated by the 2020-22 COVID shock no pre-2000 model can foresee)
3.1 Life expectancy at other ages¶
e0 summarises the whole schedule, but e65 isolates old-age mortality and is the pension-relevant measure (Meseguer, at SSA, highlights it). Life expectancy is read off the abridged life table at any group-boundary age (0, 20, 65, 80, ...). Below: the abridged values match published, e65 is forecast vs actual, and out-of-sample errors are tabulated across ages.
pub_lt = {"female": read_ltper(f"{DIR}/US fltper_1x1.txt"),
"male": read_ltper(f"{DIR}/US mltper_1x1.txt")}
def pub_ex(s, age, year):
d = pub_lt[s]; return d[(d.Year == year) & (d.Age == age)]["ex"].iloc[0]
# abridged e65 vs published (validation)
for s in ["female", "male"]:
a = ex_series(g[s], s, 65)
print(f"{s} e65: abridged {a[2000]:.2f} vs published {pub_ex(s,65,2000):.2f} (2000); "
f"{a[2020]:.2f} vs {pub_ex(s,65,2020):.2f} (2020)")
def ex_forecast(s, age):
fit = fits[s]; fc = forecast_lee_carter(fit, len(FYEARS))
idx = list(GROUP_LO).index(age)
R = project_rates(fit, fc["k"], jumpoff[s])
return np.array([life_table_abridged(R[i], s)["ex"][idx] for i in range(len(FYEARS))])
# e65 forecast vs actual
fig, ax = plt.subplots(figsize=(9, 4))
for s in ["female", "male"]:
act = ex_series(g[s], s, 65)
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=f"{s} actual")
ax.plot(FYEARS, ex_forecast(s, 65), "--", color="k", lw=1.2,
label="LC forecast" if s == "female" else None)
ax.axvline(2000, color="k", lw=.5, ls=":")
ax.set(title="e65 (life expectancy at age 65): Lee-Carter forecast vs actual",
xlabel="Year", ylabel="e65 (years)"); ax.grid(alpha=.25); ax.legend(fontsize=8)
plt.show()
# out-of-sample MAE across ages
print("\nLee-Carter out-of-sample MAE by age (years), pre-COVID 2001-2019:")
print(f"{'age':>4} | {'female':>7} | {'male':>7}")
for age in [0, 20, 65, 80]:
row = {}
for s in ["female", "male"]:
f = ex_forecast(s, age); a = ex_series(g[s], s, age).reindex(FYEARS).to_numpy()
e = np.abs(f - a)[:19] # 2001-2019
row[s] = e.mean()
print(f"{age:>4} | {row['female']:>7.2f} | {row['male']:>7.2f}")
female e65: abridged 19.04 vs published 19.05 (2000); 19.78 vs 19.78 (2020) male e65: abridged 16.06 vs published 16.05 (2000); 17.14 vs 17.13 (2020)
Lee-Carter out-of-sample MAE by age (years), pre-COVID 2001-2019: age | female | male 0 | 0.32 | 0.32 20 | 0.24 | 0.37 65 | 0.22 | 0.92 80 | 0.14 | 0.54
4. Li-Lee coherent (common-factor) model¶
The Li-Lee model shares a common factor $B(x)K(t)$ across sexes (random walk with drift) and gives each sex a deviation $b_i(x)k_i(t)$ whose index is a stationary AR(1), so the sexes cannot drift apart without limit. Note that "coherence" is an assumption — mean-reversion of the gap toward its historical average — not a law. We test it against the actual female-minus-male $e_0$ gap, in the context of the gap's history.
logm_s = {s: np.log(g[s].loc[SAMPLE[0]:SAMPLE[1]]) for s in ["female", "male"]}
cf = lee_carter_common(logm_s)
def gap_paths(horizon, years):
fc_ind = {s: forecast_lee_carter(fits[s], horizon) for s in ["female", "male"]}
e_ind = {s: e0_from_k(fits[s], fc_ind[s]["k"], s, jumpoff[s]) for s in ["female", "male"]}
fcc = forecast_common(cf, horizon)
e_com = {s: np.array([life_table_abridged(rates_from_common(cf, fcc, s, i, jumpoff[s]), s)["e0"]
for i in range(horizon)]) for s in ["female", "male"]}
return (e_ind["female"] - e_ind["male"]), (e_com["female"] - e_com["male"]), fcc
gap_ind, gap_com, fcc = gap_paths(len(FYEARS), FYEARS)
gap_all = e0_series(g["female"], "female") - e0_series(g["male"], "male") # full history 1933-2024
gap_act = gap_all.reindex(FYEARS).to_numpy()
# long horizon to 2100 to show divergence vs coherence
H = 100; LY = np.arange(2001, 2101)
gap_ind_L, gap_com_L, _ = gap_paths(H, LY)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.2))
hist = gap_all.loc[SAMPLE[0]:SAMPLE[1]]
a1.plot(hist.index, hist, color="gray", lw=1.4, label="history (sample)")
a1.plot(FYEARS, gap_act, color="k", lw=2, label="actual (out-of-sample)")
a1.plot(FYEARS, gap_ind, "--", color="#00798c", lw=1.6, label="independent LC")
a1.plot(FYEARS, gap_com, "--", color="#d1495b", lw=1.6, label="Li-Lee common factor")
pk = hist.idxmax()
a1.annotate(f"peak {int(pk)}", (pk, hist[pk]), fontsize=8,
xytext=(pk - 2, hist[pk] + 0.25))
a1.axvline(2000, color="k", lw=.5, ls=":")
a1.set(title="Female - male e0 gap: history + forecasts", xlabel="Year", ylabel="gap (years)")
a1.grid(alpha=.25); a1.legend(fontsize=7.5)
a2.plot(LY, gap_ind_L, color="#00798c", lw=1.6, label="independent LC")
a2.plot(LY, gap_com_L, color="#d1495b", lw=1.6, label="Li-Lee common factor")
a2.set(title="Projected gap to 2100 (coherence)", xlabel="Year", ylabel="gap (years)")
a2.grid(alpha=.25); a2.legend(fontsize=8)
plt.show()
print(f"AR(1) phi: " + ", ".join(f"{s} {fcc['sex'][s]['phi']:.3f}" for s in ['female','male'])
+ " (<1 => mean-reverting => bounded gap)")
print(f"2024 gap -> actual {gap_act[-1]:.2f} | independent {gap_ind[-1]:.2f} | common {gap_com[-1]:.2f}")
print(f"2100 gap -> independent {gap_ind_L[-1]:.2f} | common {gap_com_L[-1]:.2f}")
AR(1) phi: female 0.944, male 0.989 (<1 => mean-reverting => bounded gap) 2024 gap -> actual 4.88 | independent 5.80 | common 4.53 2100 gap -> independent 6.33 | common 2.80
4.1 Is the widening "implausible"? The gap in context¶
The US gap peaked at ~7.7 yr in 1975 and has narrowed since (men's smoking fell earlier/faster; cardiovascular gains helped men more), reaching 5.3 at the sample end and ~4.9 out-of-sample. Independent LC extrapolates each sex's full-sample (1933-2000) linear drift — a window dominated by the 1933-75 gap expansion — so it projects continued widening, which the 2001-24 holdout falsifies.
So the widening is not "biologically implausible"; it is (a) an artifact of the long sample and (b) contradicted by the data. And it is highly sample-dependent — fit on a recent window, independent LC narrows the gap instead (and overshoots the other way). Li-Lee happens to be closer here because its mean-reversion lands near the leveled-off gap — a better assumption for this episode, not a proven law. The deeper lesson is about linear extrapolation through a regime change.
# Sensitivity of the independent-LC gap forecast to the fit window
print("Independent-LC projected 2024 F-M e0 gap, by fit start year:")
for start in [1933, 1950, 1970, 1980]:
e = {}
for s in ["female", "male"]:
fitw = lee_carter(np.log(g[s].loc[start:SAMPLE[1]]))
fcw = forecast_lee_carter(fitw, len(FYEARS))
Rw = project_rates(fitw, fcw["k"], jumpoff[s])
e[s] = life_table_abridged(Rw[-1], s)["e0"]
print(f" fit {start}-2000 -> 2024 gap = {e['female']-e['male']:.2f}")
print(f" {'actual 2024 gap':>14} = {gap_act[-1]:.2f} "
f"(long window widens, recent window narrows, truth leveled off)")
Independent-LC projected 2024 F-M e0 gap, by fit start year: fit 1933-2000 -> 2024 gap = 5.80 fit 1950-2000 -> 2024 gap = 5.40 fit 1970-2000 -> 2024 gap = 4.01 fit 1980-2000 -> 2024 gap = 3.17 actual 2024 gap = 4.88 (long window widens, recent window narrows, truth leveled off)
4.2 e0 forecasts: independent vs common factor, by sex¶
The coherent model ties both sexes to a shared trend, so relative to independent Lee-Carter it pulls the fast-improving females down and the slower males up — the two forecasts converge. The effect is small near term and grows with horizon (shown to 2100).
H = 100; HY = np.arange(2001, 2001 + H)
def e0_ind_path(s):
fc = forecast_lee_carter(fits[s], H); R = project_rates(fits[s], fc["k"], jumpoff[s])
return np.array([life_table_abridged(R[i], s)["e0"] for i in range(H)])
def e0_com_path(s):
fcc = forecast_common(cf, H)
return np.array([life_table_abridged(rates_from_common(cf, fcc, s, i, jumpoff[s]), s)["e0"] for i in range(H)])
Ei = {s: e0_ind_path(s) for s in ["female", "male"]}
Ec = {s: e0_com_path(s) for s in ["female", "male"]}
fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.3))
for s in ["female", "male"]:
a1.plot(HY, Ei[s], "--", color=COLOR[s], lw=1.3, label=f"{s} independent")
a1.plot(HY, Ec[s], "-", color=COLOR[s], lw=2.0, label=f"{s} common factor")
a1.plot(FYEARS, e0_series(g[s], s).reindex(FYEARS), ":", color="k", lw=1,
label="actual (to 2024)" if s == "female" else None)
a1.axvline(2024, color="k", lw=.4, ls=":")
a1.set(title="e0 forecasts: independent vs Li-Lee common factor", xlabel="Year", ylabel="e0 (years)")
a1.grid(alpha=.25); a1.legend(fontsize=7.5)
for s in ["female", "male"]:
a2.plot(HY, Ec[s] - Ei[s], color=COLOR[s], lw=1.8, label=s)
a2.axhline(0, color="k", lw=.5)
a2.set(title="common factor - independent (e0 difference)", xlabel="Year", ylabel="years")
a2.grid(alpha=.25); a2.legend(fontsize=8)
plt.show()
print("common-factor minus independent e0 (years):")
for y in [2024, 2050, 2100]:
i = y - 2001
print(f" {y}: female {Ec['female'][i]-Ei['female'][i]:+.2f}, male {Ec['male'][i]-Ei['male'][i]:+.2f}")
common-factor minus independent e0 (years): 2024: female -0.51, male +0.75 2050: female -0.94, male +1.37 2100: female -1.42, male +2.11
5. Summary¶
- The 22-group abridged life table reproduces published e0 to ~0.02 yr, so the Meseguer grouping is safe.
- Lee-Carter tracks e0 well through 2019; 2020-2022 errors are the COVID shock, not model failure. It is the benchmark the BVAR must beat.
- With the 1933-2000 sample, independent LC projects the female-male gap widening, but that is a sample artifact (the window is dominated by the 1933-75 gap expansion) and is falsified by the 2001-24 holdout, where the gap narrowed. Li-Lee mean-reverts the gap and fits this episode better — a different assumption, not a proven law. The lesson: beware linear extrapolation through a regime change.
Next: the Bayesian VAR (Meseguer's method) on these same 22 groups, compared head to head against this Lee-Carter benchmark.