BVAR fertility forecasting II — the steady-state model¶

Companion to BVAR_Fertility.ipynb (Minnesota), structured like BVAR_Mortality_SteadyState.ipynb. Same data, sample, from-scratch Gibbs, and Lee-Tuljapurkar benchmark — the difference is the steady-state (Villani mean-adjusted) parametrization with an informative prior on the ultimate fertility level.

  • Model log fertility levels for 7 groups, 1933-2024; sample 1933-2000, forecast 2001-2024. Fertility settings: own-lag prior mean 0.8, $p=2$, tight $\lambda_1=0.05$, scalar $\lambda_2=0.5$.
  • Steady-state prior centred on the SSA Alternative-2 ultimate TFR $F^*=2$ (Meseguer centres $\Psi$ on the log Alt-2 ultimate rates). This is the key case for fertility, because the diffuse model mean-reverts to the baby-boom-inflated sample average.

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

1. Data — log fertility by group¶

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

DIR = "."
SAMPLE, FYEARS, FSTAR, PSI_SD = (1933, 2000), np.arange(2001, 2025), 2.0, 0.05
SS, MN, LC = "#8338ec", "#118ab2", "#c1121f"
G   = group_asfr(f"{DIR}/USAasfrTR.txt")
tfr = tfr_from_groups(G)
Gs  = G.loc[SAMPLE[0]:SAMPLE[1]]; tfrs = tfr.loc[SAMPLE[0]:SAMPLE[1]]
Y   = np.log(Gs.values)
print(f"7 groups; sample {Y.shape[0]} yrs; TFR 2000 jump-off {tfr[2000]:.2f}, actual 2024 {tfr[2024]:.2f}")
print(f"sample-mean TFR 1933-2000 = {tfrs.mean():.2f} (baby-boom lifted; the diffuse model reverts here)")

7 groups; sample 68 yrs; TFR 2000 jump-off 2.05, actual 2024 1.59
sample-mean TFR 1933-2000 = 2.44 (baby-boom lifted; the diffuse model reverts here)

2. The steady-state (mean-adjusted) BVAR and its Gibbs sampler¶

The Villani mean-adjusted VAR (same 3-block Gibbs as the mortality steady-state notebook): $$y_t = \Psi + \Pi_1(y_{t-1}-\Psi)+\dots+\Pi_p(y_{t-p}-\Psi)+\varepsilon_t,$$ with $y_t=\log f_{g,t}$. Here $\Psi$ is the steady-state vector of log fertility rates — its exponentiated sum is the ultimate TFR the forecast reverts to.

Priors. Minnesota prior on the AR coefficients $\Pi_l$ with the fertility settings (own-lag mean 0.8, $\lambda_1=0.05$, $\lambda_2=0.5$, $\lambda_3=1$); a Normal prior $\Psi\sim N(\psi_0,\Lambda_0)$ on the steady state (Section 2.1); $\Sigma\sim\mathcal{IW}$. Non-stationary draws are discarded.

3-block Gibbs. (1) $\Pi\mid\Psi,\Sigma$ — Normal on the demeaned data; (2) $\Psi\mid\Pi,\Sigma$ — Normal (the informative prior enters here); (3) $\Sigma\mid\Pi,\Psi$ — inverse-Wishart. bvar.gibbs_steady_state. Forecasting simulates around $\Psi$; rates $=\exp(\cdot)$, no cumulation.

2.1 The informative prior (ultimate TFR $F^*=2$)¶

We centre $\Psi$ on the log Alt-2 ultimate group rates: the 2000 age shape scaled to an ultimate TFR of $F^*=2$ (SSA Alternative-2, ~replacement). The prior sd $\lambda_\Psi$ sets how firmly the long run is tied to $F^*$; we use a tight value (probed in Section 6). The prior centre lies below the baby-boom-inflated sample average, so it pulls the long run down toward replacement.

In [2]:
PSI0 = np.log(G.loc[2000].values * FSTAR / tfr[2000])     # ultimate log rates, sum(exp)=F*
samp_mean_rate = np.exp(Y.mean(axis=0))                   # sample-mean group rate
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(range(7), np.exp(PSI0), "s-", color=SS, lw=1.6, ms=5, label=f"prior centre (ultimate, TFR={FSTAR})")
ax.plot(range(7), samp_mean_rate, "o--", color="gray", lw=1, ms=4, label=f"1933-2000 mean of log-rates (TFR={np.exp(Y.mean(0)).sum():.2f})")
ax.plot(range(7), G.loc[2000].values, "^:", color="k", lw=1, ms=4, label=f"observed 2000 (TFR={tfr[2000]:.2f})")
ax.set_xticks(range(7)); ax.set_xticklabels(FERT_LABELS, fontsize=8)
ax.set(title="Steady-state prior centre vs sample mean vs 2000", xlabel="Age group", ylabel="group fertility rate")
ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()
print(f"implied ultimate TFR: prior {np.exp(PSI0).sum():.2f} vs mean-of-log-rates {np.exp(Y.mean(0)).sum():.2f} (arithmetic sample-mean TFR is 2.44)")
No description has been provided for this image 
implied ultimate TFR: prior 2.00 vs mean-of-log-rates 2.34 (arithmetic sample-mean TFR is 2.44)

3. Estimation, convergence, stationarity¶

In [3]:
NDRAW, BURN = 4000, 800
t = time.time()
fit = bvar.filter_stationary(bvar.gibbs_steady_state(
    Y, p=2, lam1=0.05, lam2=0.5, lam3=1.0, own_mean=0.8,
    psi_prior_mean=PSI0, psi_prior_sd=PSI_SD, ndraw=NDRAW, burn=BURN, seed=1))
print(f"{NDRAW} draws in {time.time()-t:.1f}s; stationary {fit['stationary_frac']:.3f} -> {len(fit['Sigma'])} kept")

fig, axes = plt.subplots(1, 3, figsize=(12, 3))
axes[0].plot(fit["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$")
axes[1].plot(fit["Pi"][:, 1, 1], lw=.5); axes[1].set_title("$\Pi$ own lag, 20-24")
axes[2].plot(np.exp(fit["mu"]).sum(1), lw=.5); axes[2].set_title("implied ultimate TFR $\sum e^{\Psi}$")
for a in axes: a.set_xlabel("kept draw"); a.grid(alpha=.2)
fig.suptitle("Trace plots (steady-state fertility BVAR, p=2)", y=1.03); fig.tight_layout(); plt.show()

<>:9: SyntaxWarning: invalid escape sequence '\S'
<>:10: SyntaxWarning: invalid escape sequence '\P'
<>:11: SyntaxWarning: invalid escape sequence '\s'
<>:9: SyntaxWarning: invalid escape sequence '\S'
<>:10: SyntaxWarning: invalid escape sequence '\P'
<>:11: SyntaxWarning: invalid escape sequence '\s'
C:\Users\user\AppData\Local\Temp\ipykernel_42088\2426170431.py:9: SyntaxWarning: invalid escape sequence '\S'
  axes[0].plot(fit["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$")
C:\Users\user\AppData\Local\Temp\ipykernel_42088\2426170431.py:10: SyntaxWarning: invalid escape sequence '\P'
  axes[1].plot(fit["Pi"][:, 1, 1], lw=.5); axes[1].set_title("$\Pi$ own lag, 20-24")
C:\Users\user\AppData\Local\Temp\ipykernel_42088\2426170431.py:11: SyntaxWarning: invalid escape sequence '\s'
  axes[2].plot(np.exp(fit["mu"]).sum(1), lw=.5); axes[2].set_title("implied ultimate TFR $\sum e^{\Psi}$")

4000 draws in 1.4s; stationary 0.997 -> 3989 kept
No description has been provided for this image

3.1 The estimated steady state $\Psi$¶

The posterior of the implied ultimate TFR $\sum_g e^{\Psi_g}$ blends the prior ($F^*=2$) with the data (which would pull it to the ~2.4 sample mean). With the tight prior it settles near $F^*$.

In [4]:
ult = np.exp(fit["mu"]).sum(1)                    # implied ultimate TFR per draw
fig, ax = plt.subplots(figsize=(8, 3.4))
ax.hist(ult, bins=40, color=SS, alpha=.5, density=True)
ax.axvline(FSTAR, color=LC, lw=1.5, ls="--", label=f"prior centre F*={FSTAR}")
ax.axvline(tfrs.mean(), color="gray", lw=1.5, ls=":", label=f"sample mean {tfrs.mean():.2f}")
ax.axvline(np.median(ult), color=SS, lw=1.5, label=f"posterior median {np.median(ult):.2f}")
ax.set(title="Posterior of the implied ultimate TFR", xlabel="ultimate TFR"); ax.legend(fontsize=8); ax.grid(alpha=.2)
plt.show()
No description has been provided for this image

4. Forecasts vs. actual and the benchmark¶

TFR and all 7 group rates (median + 90% band) vs actual, the Minnesota BVAR, and the Lee-Tuljapurkar benchmark.

In [5]:
def bands(fit):
    R = np.exp(bvar.forecast(fit, Y, len(FYEARS)))
    return np.percentile(R, [5, 50, 95], axis=0), np.percentile(R.sum(axis=2), [5, 50, 95], axis=0)

gb, tb = bands(fit)                               # steady-state
mn = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=2, lam1=0.05, lam2=0.5, lam3=1.0,
                                                 own_mean=0.8, ndraw=NDRAW, burn=BURN, seed=1))
mgb, mtb = bands(mn)                              # Minnesota (diffuse long run)
lcfit = lee_carter(np.log(Gs)); bmean, blo, bhi, _ = constrained_arma_forecast(tfrs, FSTAR, len(FYEARS))
bgr = reconstruct_group_rates(G.loc[2000].values, lcfit["b"].values, bmean)

fig, axes = plt.subplots(2, 4, figsize=(15, 6.8)); 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="#333", lw=1.8, label="actual")
    ax.fill_between(FYEARS, gb[0, :, k], gb[2, :, k], color=SS, alpha=.15, label="steady-state 90%")
    ax.plot(FYEARS, gb[1, :, k], "-", color=SS, lw=1.4, label="steady-state median")
    ax.plot(FYEARS, mgb[1, :, k], "--", color=MN, lw=1.1, label="Minnesota median")
    ax.plot(FYEARS, bgr[:, k], ":", color=LC, lw=1.2, label="benchmark median")
    ax.axvline(2000, color="k", lw=.4, ls=":"); ax.set_title(f"group {gname}"); ax.grid(alpha=.25)
ax = axf[7]
ax.plot(tfr.loc[1980:2000].index, tfr.loc[1980:2000], color="gray", lw=1)
ax.plot(FYEARS, tfr.reindex(FYEARS), color="#333", lw=1.8, label="actual")
ax.fill_between(FYEARS, tb[0], tb[2], color=SS, alpha=.15)
ax.plot(FYEARS, tb[1], "-", color=SS, lw=1.4)
ax.plot(FYEARS, mtb[1], "--", color=MN, lw=1.1)
ax.plot(FYEARS, bmean, ":", color=LC, lw=1.2)
ax.axhline(FSTAR, color=LC, lw=.5, ls=":"); ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set_title("TFR (total)"); ax.grid(alpha=.25)
axf[0].legend(fontsize=6.5)
fig.suptitle("Steady-state fertility BVAR vs Minnesota vs benchmark vs actual", y=1.0)
fig.tight_layout(); plt.show()
No description has been provided for this image

5. Lag length — $p=1,2,3$¶

In [6]:
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(tfr.loc[1990:2000].index, tfr.loc[1990:2000], color="gray", lw=1)
ax.plot(FYEARS, tfr.reindex(FYEARS), color="#333", lw=2, label="actual")
for p, c in [(1, "#222"), (2, SS), (3, "#3a7d44")]:
    fp = bvar.filter_stationary(bvar.gibbs_steady_state(Y, p=p, lam1=0.05, lam2=0.5, lam3=1.0, own_mean=0.8,
                                                        psi_prior_mean=PSI0, psi_prior_sd=PSI_SD, ndraw=2500, burn=500, seed=2))
    med = np.percentile(np.exp(bvar.forecast(fp, Y, len(FYEARS))).sum(axis=2), 50, axis=0)
    ax.plot(FYEARS, med, "--", color=c, lw=1.4, label=f"p={p}")
ax.axhline(FSTAR, color=LC, lw=.5, ls=":"); ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="TFR: steady-state fertility BVAR at p=1,2,3", xlabel="Year", ylabel="TFR")
ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()
No description has been provided for this image

6. Prior tightness — how firmly to trust $F^*$¶

$\lambda_\Psi$ controls the informativeness: tight pins the long run to $F^*=2$; diffuse lets the data dominate and the model collapses to the (over-predicting) Minnesota result at the baby-boom-inflated mean.

In [7]:
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(FYEARS, tfr.reindex(FYEARS), color="#333", lw=2, label="actual")
for sd, c in [(0.02, "#3a0ca3"), (0.05, "#8338ec"), (0.1, "#b56cff"), (5.0, "#c9b6e4")]:
    fsd = bvar.filter_stationary(bvar.gibbs_steady_state(Y, p=2, lam1=0.05, lam2=0.5, lam3=1.0, own_mean=0.8,
                                                         psi_prior_mean=PSI0, psi_prior_sd=sd, ndraw=2000, burn=400, seed=3))
    med = np.percentile(np.exp(bvar.forecast(fsd, Y, len(FYEARS))).sum(axis=2), 50, axis=0)
    lab = f"$\lambda_\Psi$={sd}" + (" (diffuse)" if sd >= 5 else "")
    ax.plot(FYEARS, med, "--", color=c, lw=1.4, label=lab)
ax.axhline(FSTAR, color=LC, lw=.5, ls=":"); ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="TFR: steady-state forecast vs prior tightness $\lambda_\Psi$ (F*=2)",
       xlabel="Year", ylabel="TFR"); ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()
print("Tight -> reverts to F*=2; diffuse -> reverts to the ~2.4 sample mean (= Minnesota).")

<>:7: SyntaxWarning: invalid escape sequence '\l'
<>:10: SyntaxWarning: invalid escape sequence '\l'
<>:7: SyntaxWarning: invalid escape sequence '\l'
<>:10: SyntaxWarning: invalid escape sequence '\l'
C:\Users\user\AppData\Local\Temp\ipykernel_42088\3146925271.py:7: SyntaxWarning: invalid escape sequence '\l'
  lab = f"$\lambda_\Psi$={sd}" + (" (diffuse)" if sd >= 5 else "")
C:\Users\user\AppData\Local\Temp\ipykernel_42088\3146925271.py:10: SyntaxWarning: invalid escape sequence '\l'
  ax.set(title="TFR: steady-state forecast vs prior tightness $\lambda_\Psi$ (F*=2)",
No description has been provided for this image 
Tight -> reverts to F*=2; diffuse -> reverts to the ~2.4 sample mean (= Minnesota).

7. Model comparison — TFR forecast error¶

In [8]:
act = tfr.reindex(FYEARS).values
rows = {"Steady-state (F*=2, tight)": tb[1], "Minnesota (diffuse)": mtb[1], "Benchmark (Lee-Tuljapurkar)": bmean}
print(f"{'model':>28} | {'TFR MAE 2001-24':>15} | {'TFR 2024':>8}")
for name, med in rows.items():
    print(f"{name:>28} | {np.abs(med-act).mean():15.3f} | {med[-1]:8.2f}")
print(f"{'actual 2024':>28} | {'':>15} | {act[-1]:8.2f}")

                       model | TFR MAE 2001-24 | TFR 2024
  Steady-state (F*=2, tight) |           0.189 |     2.04
         Minnesota (diffuse) |           0.347 |     2.27
 Benchmark (Lee-Tuljapurkar) |           0.202 |     2.04
                 actual 2024 |                 |     1.59

8. Summary¶

  • The steady-state fertility BVAR reparametrizes around the long-run mean $\Psi$, letting us place an informative prior on the ultimate TFR ($F^*=2$).
  • Unlike mortality, the informative prior clearly helps here: the diffuse Minnesota model reverts to the baby-boom-inflated ~2.4 sample mean and badly over-predicts, whereas anchoring the long run near replacement (with a tight $\lambda_\Psi$) brings the TFR forecast down to ~2.0 and the lowest error of the three models — though the margin over the Lee–Tuljapurkar benchmark (0.189 vs 0.202) is slim; the clear gain is over the diffuse Minnesota model.
  • The result is dial-able: Section 6 shows loosening $\lambda_\Psi$ recovers the over-predicting Minnesota result, and $F^*$ itself is the key lever (a lower $F^*$ would track the actual post-2007 decline better).
  • Both BVARs still inherit the age-pattern (postponement) miss, which a fixed-shape long run cannot capture.