Fertility forecasting — the Lee-Tuljapurkar benchmark

The fertility counterpart of LC_Mortality.ipynb: the benchmark Meseguer (2010) compares his fertility BVAR against. Following Lee (1993) / Lee-Tuljapurkar (1994):

• Model fertility for 7 age groups (<20, 20-24, ..., 45+), 1933-2024, sampling 1933-2000 and forecasting 2001-2024.
• A Lee-Carter age decomposition of the log rates gives the age pattern, but the total fertility rate (TFR) is not a random walk (fertility is mean-reverting, not trending). Instead it follows a constrained ARMA(1,1) with a fixed ultimate mean $F^*$ — Meseguer sets $F^*=2$, the SSA Alternative-2 assumption (roughly replacement).

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

1. Data — age-specific fertility and the TFR

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from fertility import (group_asfr, tfr_from_groups, constrained_arma_forecast,
                       rw_drift_forecast, reconstruct_group_rates, FERT_LABELS)
from mortality import lee_carter
%matplotlib inline
plt.rcParams["figure.dpi"] = 110

DIR = "."
SAMPLE, FYEARS, FSTAR = (1933, 2000), np.arange(2001, 2025), 2.0
G   = group_asfr(f"{DIR}/USAasfrTR.txt")      # Year x 7 groups
tfr = tfr_from_groups(G)
print(f"7 fertility groups: {FERT_LABELS}")
print(f"years {G.index.min()}-{G.index.max()}; sample {SAMPLE[0]}-{SAMPLE[1]}")
print(f"TFR: 1960 {tfr[1960]:.2f}, 2000 {tfr[2000]:.2f}, 2010 {tfr[2010]:.2f}, 2024 {tfr[2024]:.2f}")

7 fertility groups: ['<20', '20-24', '25-29', '30-34', '35-39', '40-44', '45+']
years 1933-2024; sample 1933-2000
TFR: 1960 3.67, 2000 2.05, 2010 1.92, 2024 1.59

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.2))
a1.plot(tfr.index, tfr, color="#6a4c93", lw=1.8)
a1.axvspan(SAMPLE[0], SAMPLE[1], color="gray", alpha=.08, label="sample")
a1.axhline(2.1, color="gray", ls=":", lw=1, label="replacement ~2.1")
a1.axhline(FSTAR, color="#c1121f", ls="--", lw=1, label="F* = 2 (SSA Alt-2)")
a1.set(title="US Total Fertility Rate, 1933-2024", xlabel="Year", ylabel="TFR")
a1.grid(alpha=.25); a1.legend(fontsize=8)

for y, c in [(1960, "#9aa0b0"), (1980, "#f4a259"), (2000, "#5b8e7d"), (2020, "#6a4c93")]:
    a2.plot(range(7), G.loc[y].values, "o-", color=c, ms=4, label=str(y))
a2.set_xticks(range(7)); a2.set_xticklabels(FERT_LABELS, fontsize=8)
a2.set(title="Age pattern of fertility (postponement to older ages)",
       xlabel="Age group", ylabel="group fertility rate")
a2.grid(alpha=.25); a2.legend(fontsize=8, title="Year")
plt.show()

2. The Lee-Tuljapurkar benchmark

Age decomposition. As in mortality, a Lee-Carter SVD of the log group rates, $\log f_g(t) = a_g + b_g\,k_t$, gives the average age profile $a_g$, the age pattern of change $b_g$ ($\sum_g b_g = 1$), and a single index $k_t$.

Constrained TFR (Meseguer eq. 8). A random-walk $k_t$ would let the TFR drift without bound — implausible for fertility. Following Lee (1993), the TFR is instead modelled as an ARMA(1,1) whose unconditional mean is fixed at an ultimate value $F^*$: $$\text{TFR}_t = F^*(1-\phi) + \phi\,\text{TFR}_{t-1} + e_t + \theta e_{t-1},$$ so the forecast reverts to $F^*$ (we use $F^*=2$). We fit this to the sample and forecast the TFR, then distribute it across the 7 groups using the age pattern $b_g$, anchored to the observed 2000 rates so the age-specific forecasts join the data (a Lee-Miller jump-off adjustment): $f_g(t)=f_g^{obs}(2000)\,e^{b_g\delta_t}$ with $\delta_t$ chosen so $\sum_g f_g(t)=\text{TFR}_t$. For contrast we also show an unconstrained random-walk-with-drift TFR forecast.

3. Fit — age pattern and index, 1933-2000

Gs = G.loc[SAMPLE[0]:SAMPLE[1]]; tfrs = tfr.loc[SAMPLE[0]:SAMPLE[1]]
fit = lee_carter(np.log(Gs))
fig, axes = plt.subplots(1, 3, figsize=(14, 3.6))
axes[0].plot(range(7), fit["a"].values, "o-", color="#6a4c93"); axes[0].set_title("$a_g$ (avg log rate)")
axes[1].plot(range(7), fit["b"].values, "o-", color="#6a4c93"); axes[1].set_title("$b_g$ (age pattern, sums to 1)")
for ax in axes[:2]:
    ax.set_xticks(range(7)); ax.set_xticklabels(FERT_LABELS, rotation=45, fontsize=7); ax.grid(alpha=.25)
axes[2].plot(fit["k"].index, fit["k"].values, color="#6a4c93"); axes[2].set_title("$k_t$ (index)")
axes[2].set_xlabel("Year"); axes[2].grid(alpha=.25)
fig.tight_layout(); plt.show()
print("Note: on log fertility the SVD's b_g is pulled toward the small, volatile "
      "tail groups (40-44, 45+) -- a known limitation of naive Lee-Carter for fertility.")

Note: on log fertility the SVD's b_g is pulled toward the small, volatile tail groups (40-44, 45+) -- a known limitation of naive Lee-Carter for fertility.

4. Forecast — TFR and age-specific rates vs. actual

mean, lo, hi, params = constrained_arma_forecast(tfrs, FSTAR, len(FYEARS))
rw = rw_drift_forecast(tfrs, len(FYEARS))
act = tfr.reindex(FYEARS)

fig, ax = plt.subplots(figsize=(9, 4.2))
ax.plot(tfrs.index, tfrs, color="gray", lw=1.2, label="sample")
ax.plot(FYEARS, act, color="#6a4c93", lw=2.2, label="actual")
ax.fill_between(FYEARS, lo, hi, color="#c1121f", alpha=.12, label="F*=2 ARMA 90%")
ax.plot(FYEARS, mean, "--", color="#c1121f", lw=1.6, label="F*=2 ARMA median")
ax.plot(FYEARS, rw, ":", color="#457b9d", lw=1.6, label="random walk (drift)")
ax.axhline(FSTAR, color="#c1121f", lw=.6, ls=":")
ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="TFR forecast: constrained ARMA (F*=2) vs random walk vs actual",
       xlabel="Year", ylabel="TFR"); ax.grid(alpha=.25); ax.legend(fontsize=8)
plt.show()

mae_c = np.abs(mean - act.values).mean(); mae_rw = np.abs(rw - act.values).mean()
print(f"TFR forecast MAE 2001-2024: constrained F*=2 {mae_c:.3f} | random walk {mae_rw:.3f}")
print(f"Both revert to ~2 but actual fell to {act.iloc[-1]:.2f} -> the F*=2 assumption "
      f"over-predicts (post-2007 fertility decline, a regime change absent from the sample).")

TFR forecast MAE 2001-2024: constrained F*=2 0.202 | random walk 0.203
Both revert to ~2 but actual fell to 1.59 -> the F*=2 assumption over-predicts (post-2007 fertility decline, a regime change absent from the sample).

# age-specific: distribute the constrained TFR path (and its 90% bounds) across the
# 7 groups, anchored to the observed 2000 rates (jump-off adjusted)
base, bpat = G.loc[2000].values, fit["b"].values
gr_mid = reconstruct_group_rates(base, bpat, mean)
gr_lo  = reconstruct_group_rates(base, bpat, lo)
gr_hi  = reconstruct_group_rates(base, bpat, hi)

fig, axes = plt.subplots(2, 4, figsize=(15, 6.5)); axf = axes.ravel()
for k, gname in enumerate(FERT_LABELS):
    ax = axf[k]
    ax.plot(G.loc[1980:2000].index, G.loc[1980:2000, gname], color="gray", lw=1)
    ax.plot(FYEARS, G.loc[2001:2024, gname], color="#6a4c93", lw=1.8, label="actual")
    ax.fill_between(FYEARS, gr_lo[:, k], gr_hi[:, k], color="#c1121f", alpha=.15, label="benchmark 90%")
    ax.plot(FYEARS, gr_mid[:, k], "--", color="#c1121f", lw=1.4, label="benchmark median")
    ax.axvline(2000, color="k", lw=.4, ls=":"); ax.set_title(f"group {gname}"); ax.grid(alpha=.25)
ax = axf[7]                                        # 8th panel: the TFR total
ax.plot(tfr.loc[1980:2000].index, tfr.loc[1980:2000], color="gray", lw=1)
ax.plot(FYEARS, act, color="#6a4c93", lw=1.8, label="actual")
ax.fill_between(FYEARS, lo, hi, color="#c1121f", alpha=.15)
ax.plot(FYEARS, mean, "--", color="#c1121f", lw=1.4)
ax.axhline(FSTAR, color="#c1121f", lw=.5, ls=":"); ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set_title("TFR (total)"); ax.grid(alpha=.25)
axf[0].set_ylabel("group fertility rate"); axf[4].set_ylabel("group fertility rate")
axf[0].legend(fontsize=7.5)
fig.suptitle("Fertility benchmark (F*=2): all 7 groups + TFR, with 90% intervals", y=1.0)
fig.tight_layout(); plt.show()
print("The benchmark misses both the level (TFR too high) and the shift toward older "
      "ages: younger groups (<20, 20-24) are badly over-predicted while the actual "
      "fell; older groups (30-34, 35-39) rose above the fixed-age-pattern forecast.")

The benchmark misses both the level (TFR too high) and the shift toward older ages: younger groups (<20, 20-24) are badly over-predicted while the actual fell; older groups (30-34, 35-39) rose above the fixed-age-pattern forecast.

5. Summary

• The fertility benchmark is Lee-Tuljapurkar: a Lee-Carter age pattern plus a constrained ARMA(1,1) that reverts the TFR to an ultimate $F^*=2$.
• With the 1933-2000 sample (TFR near replacement at the jump-off), it forecasts the TFR staying near 2 -- but actual fertility fell to ~1.6 after 2007, so the benchmark over-predicts. The revert-to-$F^*$ assumption is doing the work, and it was set too high for what happened.
• Naive Lee-Carter on log fertility also handles the age pattern poorly (the $b_g$ is dominated by small, volatile tail groups), motivating a joint model.

Next: the fertility BVAR (Minnesota and steady-state) on the same 7 groups, benchmarked against this -- the steady-state's informative prior centres on the ultimate TFR, directly comparable to $F^*$.