BVAR fertility forecasting I — the Minnesota-prior model

The fertility counterpart of BVAR_Mortality.ipynb, benchmarked against LC_Fertility.ipynb (Lee-Tuljapurkar). Estimated with the same from-scratch Gibbs sampler (bvar.py), but with fertility-specific settings from Meseguer (2010, Section 5.1) — fertility is mean-reverting, not trending, so the model differs from mortality in three ways:

• modelled on log fertility levels $\log f_{g,t}$ (not differences), 7 age groups, 1933-2024, sample 1933-2000, forecast 2001-2024;
• the Minnesota prior centres the own first lag on a stationary AR(0.8) (not a random walk), with $p=2$ lags;
• a tight $\lambda_1=0.05$ (log fertility can otherwise induce explosive draws), scalar $\lambda_2=0.5$, $\lambda_3=1$; non-stationary draws are discarded.

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

1. Data — log fertility by group

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 = (1933, 2000), np.arange(2001, 2025), 2.0
BV, LC = "#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)                      # log fertility levels, (T x 7)
print(f"7 groups {FERT_LABELS}; sample {Y.shape[0]} yrs x 7; forecast {FYEARS[0]}-{FYEARS[-1]}")

7 groups ['<20', '20-24', '25-29', '30-34', '35-39', '40-44', '45+']; sample 68 yrs x 7; forecast 2001-2024

fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.2))
for k, gname in enumerate(FERT_LABELS):
    a1.plot(G.index, np.log(G[gname]), lw=1, label=gname)
a1.axvspan(*SAMPLE, color="gray", alpha=.08)
a1.set(title="log fertility rate by age group", xlabel="Year", ylabel="log f")
a1.grid(alpha=.25); a1.legend(fontsize=7, ncol=2)
# correlation of the 7 groups (levels)
Cf = np.corrcoef(Y.T)
im = a2.imshow(Cf, cmap="RdBu_r", vmin=-1, vmax=1)
a2.set_xticks(range(7)); a2.set_yticks(range(7))
a2.set_xticklabels(FERT_LABELS, rotation=90, fontsize=7); a2.set_yticklabels(FERT_LABELS, fontsize=7)
a2.set_title("group correlation (log f)")
fig.colorbar(im, ax=a2, fraction=.046, pad=.04)
fig.tight_layout(); plt.show()
print("The 7 groups show the same contiguity band as mortality (adjacent ~0.9, "
      "most-distant <20/45+ ~0.27) and are all POSITIVELY correlated, so the "
      "'spatial trick' would be well-defined. But with only 7 groups -- and even "
      "'distant' pairs still moderately correlated -- its benefit is marginal, so "
      "Meseguer (and we) use a scalar lam2=0.5 instead.")

The 7 groups show the same contiguity band as mortality (adjacent ~0.9, most-distant <20/45+ ~0.27) and are all POSITIVELY correlated, so the 'spatial trick' would be well-defined. But with only 7 groups -- and even 'distant' pairs still moderately correlated -- its benefit is marginal, so Meseguer (and we) use a scalar lam2=0.5 instead.

2. The model and its Gibbs sampler

Same reduced-form Minnesota BVAR and 2-block Gibbs sampler as the mortality notebook — $y_t = c + \Pi_1 y_{t-1}+\dots+\Pi_p y_{t-p}+\varepsilon_t$, with $\mathrm{vec}(B)\mid\Sigma \sim N$ and $\Sigma\mid B \sim \mathcal{IW}$ — but here $y_t=\log f_{g,t}$ (levels) and the Minnesota prior is set for a mean-reverting series:

Prior means. The own first lag is centred on 0.8 (a stationary AR(0.8)), not on 1 (random walk) or 0 (differences). Fertility does not trend, so a random-walk prior is inappropriate; AR(0.8) encodes gradual mean reversion.

Prior variances (Minnesota). $\sigma_{ijl}=\lambda_1/l^{\lambda_3}$ (own) and $(\lambda_1\lambda_2/l^{\lambda_3})(s_i/s_j)$ (cross), with $\lambda_1=0.05$ (tight — loose priors on log fertility produce explosive parameter draws), $\lambda_2=0.5$ (scalar), $\lambda_3=1$, diffuse intercept. $\Sigma\sim\mathcal{IW}$.

Stationarity. The unconditional mean of a stationary VAR is a well-defined long-run fertility level, but explosive draws would blow up over a 24-year forecast. Following Meseguer we discard non-stationary draws (companion-matrix eigenvalues outside the unit circle) — bvar.filter_stationary.

Forecast. Because we model levels, the simulated paths are $\log f$ directly (no cumulation): rates $=\exp(\cdot)$ and $\mathrm{TFR}=\sum_g f_g$.

3. Estimation, convergence, stationarity

NDRAW, BURN = 4000, 800
t = time.time()
raw = 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)
fit = bvar.filter_stationary(raw)
print(f"{NDRAW} draws in {time.time()-t:.1f}s; stationary fraction {fit['stationary_frac']:.3f} "
      f"-> {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["B"][:, 1, 1], lw=.5); axes[1].set_title("own lag-1, group 20-24")
axes[2].plot(fit["B"][:, 2, 1], lw=.5); axes[2].set_title("cross 25-29 -> 20-24")
for a in axes: a.set_xlabel("kept draw"); a.grid(alpha=.2)
fig.suptitle("Trace plots (Minnesota fertility BVAR, p=2)", y=1.03); fig.tight_layout(); plt.show()

<>:10: SyntaxWarning: invalid escape sequence '\S'
<>:10: SyntaxWarning: invalid escape sequence '\S'
C:\Users\user\AppData\Local\Temp\ipykernel_39296\1643195348.py:10: SyntaxWarning: invalid escape sequence '\S'
  axes[0].plot(fit["Sigma"][:, 0, 0], lw=.5); axes[0].set_title("$\Sigma_{0,0}$")

4000 draws in 1.5s; stationary fraction 0.999 -> 3995 kept

4. Forecasts vs. actual and the Lee-Tuljapurkar benchmark

Posterior predictive TFR and all 7 group rates (median + 90% band), against the actual values and the benchmark (constrained ARMA $F^*=2$ for the TFR, distributed by the LC age pattern) — both carry 90% bands (BVAR = full posterior; benchmark = the ARMA TFR interval propagated through the fixed age pattern). The BVAR mean-reverts toward the sample average. The baby-boom years (1946-64, TFR ~3.4) within the 1933-2000 sample lift that average to ~2.44 (it would be ~2.07 without them, near the 2000 jump-off of 2.05), so the AR(0.8) forecasts fertility rising back toward ~2.3 — over-shooting the actual decline even more than the benchmark.

def group_bands(fit):
    R = np.exp(bvar.forecast(fit, Y, len(FYEARS)))          # (nd, H, 7)
    return np.percentile(R, [5, 50, 95], axis=0)            # (3, H, 7)

bnd = group_bands(fit)                                       # BVAR group bands
tfr_bnd = np.percentile(np.exp(bvar.forecast(fit, Y, len(FYEARS))).sum(axis=2), [5, 50, 95], axis=0)

# Lee-Tuljapurkar benchmark
lcfit = lee_carter(np.log(Gs))
bmean, blo, bhi, _ = constrained_arma_forecast(tfrs, FSTAR, len(FYEARS))
base_g = G.loc[2000].values
bgr    = reconstruct_group_rates(base_g, lcfit["b"].values, bmean)
bgr_lo = reconstruct_group_rates(base_g, lcfit["b"].values, blo)
bgr_hi = reconstruct_group_rates(base_g, lcfit["b"].values, bhi)

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, bnd[0, :, k], bnd[2, :, k], color=BV, alpha=.15, label="BVAR 90%")
    ax.plot(FYEARS, bnd[1, :, k], "--", color=BV, lw=1.4, label="BVAR median")
    ax.fill_between(FYEARS, bgr_lo[:, k], bgr_hi[:, k], color=LC, alpha=.10, label="benchmark 90%")
    ax.plot(FYEARS, bgr[:, k], ":", color=LC, 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]
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, tfr_bnd[0], tfr_bnd[2], color=BV, alpha=.15)
ax.plot(FYEARS, tfr_bnd[1], "--", color=BV, lw=1.4)
ax.fill_between(FYEARS, blo, bhi, color=LC, alpha=.10)
ax.plot(FYEARS, bmean, ":", color=LC, lw=1.4)
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=7)
fig.suptitle("Minnesota fertility BVAR vs benchmark vs actual (90% bands)", y=1.0)
fig.tight_layout(); plt.show()

act = tfr.reindex(FYEARS).values
print(f"TFR forecast MAE 2001-2024: BVAR {np.abs(tfr_bnd[1]-act).mean():.3f} | "
      f"benchmark(F*=2) {np.abs(bmean-act).mean():.3f} | actual 2024 = {act[-1]:.2f}")

TFR forecast MAE 2001-2024: BVAR 0.347 | benchmark(F*=2) 0.202 | actual 2024 = 1.59

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

Meseguer uses $p=2$. With the tight $\lambda_1$ and stationarity filter the TFR median is stable across lag order.

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, "#118ab2"), (3, "#3a7d44")]:
    fp = bvar.filter_stationary(bvar.gibbs_minnesota(Y, p=p, lam1=0.05, lam2=0.5, lam3=1.0,
                                                     own_mean=0.8, 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"BVAR p={p}")
ax.axhline(FSTAR, color=LC, lw=.5, ls=":"); ax.axvline(2000, color="k", lw=.4, ls=":")
ax.set(title="TFR: Minnesota fertility BVAR at p=1,2,3", xlabel="Year", ylabel="TFR")
ax.grid(alpha=.25); ax.legend(fontsize=8); plt.show()

6. Summary

• The Minnesota fertility BVAR models log fertility levels with an AR(0.8)-centred prior ($p=2$, tight $\lambda_1=0.05$, scalar $\lambda_2$), estimated by the same Gibbs sampler with non-stationary draws discarded (~99.9% kept here).
• It gives full posterior credible bands for every group and the TFR.
• But with a diffuse long run it mean-reverts toward the 1933-2000 sample average (~2.44) — lifted by the 1946-64 baby boom, versus ~2.07 without it — so it forecasts the TFR rising above the 2000 jump-off and over-predicts even more than the benchmark. It also inherits the age-pattern miss (postponement).
• This is exactly the case for an informative steady-state prior on the ultimate fertility level (the fertility analogue of $F^*$) — the next notebook.