Forecast performance — fertility (expanding-window backtest)

The fertility counterpart of Forecast_Performance_Mortality.ipynb: a pseudo-real-time backtest with an expanding window and non-overlapping 10-year forecast blocks (fit to 1970 -> forecast 1971-1980; ... fit to 2020), so every year 1971-2024 is forecast once at a horizon of 1-10 years.

• Three models: Lee-Tuljapurkar (constrained ARMA, $F^*=2$), Minnesota BVAR (diffuse long run), steady-state BVAR ($F^*=2$ prior), scored on the TFR with 90% intervals.
• Output: the realized TFR with the forecasts made along the way, error by horizon, accuracy + calibration scores, and a 1- vs 10-year-ahead age-profile view.

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

1. Setup

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 = "."
ORIGINS, H, FSTAR = [1970, 1980, 1990, 2000, 2010, 2020], 10, 2.0
MODELS = ["Lee-Tuljapurkar", "Minnesota BVAR", "Steady-state BVAR"]
MCOL = {"Lee-Tuljapurkar": "#c1121f", "Minnesota BVAR": "#118ab2", "Steady-state BVAR": "#8338ec"}
NB, BB = 2500, 500

G   = group_asfr(f"{DIR}/USAasfrTR.txt")
tfr = tfr_from_groups(G)

def _bvar_fit(model, origin):
    Y = np.log(G.loc[1933:origin].values)
    if model == "Minnesota BVAR":
        fit = bvar.gibbs_minnesota(Y, p=2, lam1=0.05, lam2=0.5, lam3=1.0, own_mean=0.8,
                                   ndraw=NB, burn=BB, seed=1)
    else:
        psi0 = np.log(G.loc[origin].values * FSTAR / tfr[origin])
        fit = 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=0.05, ndraw=NB, burn=BB, seed=1)
    return bvar.filter_stationary(fit), Y

def forecast_tfr(model, origin, H=H):
    if model == "Lee-Tuljapurkar":
        mean, lo, hi, _ = constrained_arma_forecast(tfr.loc[1933:origin], FSTAR, H)
        return lo, mean, hi
    fit, Y = _bvar_fit(model, origin)
    T = np.exp(bvar.forecast(fit, Y, H)).sum(axis=2)          # (nd, H) TFR draws
    return tuple(np.percentile(T, [5, 50, 95], axis=0))

def forecast_group_profile(model, origin, horizon):
    if model == "Lee-Tuljapurkar":
        b = lee_carter(np.log(G.loc[1933:origin]))["b"].values
        mean, lo, hi, _ = constrained_arma_forecast(tfr.loc[1933:origin], FSTAR, horizon)
        base = G.loc[origin].values
        f = lambda path: reconstruct_group_rates(base, b, path)[horizon - 1]
        return f(lo), f(mean), f(hi)
    fit, Y = _bvar_fit(model, origin)
    R = np.exp(bvar.forecast(fit, Y, horizon))[:, horizon - 1, :]
    return tuple(np.percentile(R, [5, 50, 95], axis=0))

print(f"origins {ORIGINS}; models {MODELS}; horizon {H}")

origins [1970, 1980, 1990, 2000, 2010, 2020]; models ['Lee-Tuljapurkar', 'Minnesota BVAR', 'Steady-state BVAR']; horizon 10

2. Run the backtest

t0 = time.time(); rows = []
for origin in ORIGINS:
    for model in MODELS:
        lo, mid, hi = forecast_tfr(model, origin)
        for h in range(H):
            yr = origin + h + 1
            rows.append(dict(model=model, origin=origin, h=h + 1, year=yr,
                             tfr_lo=lo[h], tfr_fc=mid[h], tfr_hi=hi[h], tfr_act=tfr.get(yr, np.nan)))
R = pd.DataFrame(rows)
R["err"] = R["tfr_fc"] - R["tfr_act"]
R["covered"] = (R["tfr_act"] >= R["tfr_lo"]) & (R["tfr_act"] <= R["tfr_hi"])
print(f"done in {time.time()-t0:.0f}s; {R['tfr_act'].notna().sum()} scored forecasts")

done in 13s; 162 scored forecasts

3. Realized TFR and the forecasts made along the way

The realized TFR (black) with each model's 10-year forecast blocks (± 90% band) drawn from the origin where they'd have been made.

fig, axes = plt.subplots(1, 3, figsize=(15, 4.4), sharey=True)
for ax, model in zip(axes, MODELS):
    ax.plot(tfr.loc[1965:2024].index, tfr.loc[1965:2024], "k-", lw=1.8, zorder=5, label="realized")
    for origin in ORIGINS:
        d = R[(R.model == model) & (R.origin == origin)]
        ax.fill_between(d["year"], d["tfr_lo"], d["tfr_hi"], color=MCOL[model], alpha=.13,
                        label="90% band" if origin == ORIGINS[0] else None)
        ax.plot(d["year"], d["tfr_fc"], "-", color=MCOL[model], lw=1.4,
                label="forecast median" if origin == ORIGINS[0] else None)
        ax.plot(d["year"].iloc[0], d["tfr_fc"].iloc[0], "o", color=MCOL[model], ms=3.5)
    ax.axhline(FSTAR, color="gray", lw=.5, ls=":");
    for o in ORIGINS: ax.axvline(o, color="gray", lw=.3, ls=":")
    ax.set_title(model, fontsize=10); ax.set_xlabel("Year"); ax.grid(alpha=.25)
axes[0].set_ylabel("TFR"); axes[0].legend(fontsize=8)
fig.suptitle("Realized TFR and the 10-year forecasts made along the way", y=1.02)
fig.tight_layout(); plt.show()

4. Error by forecast horizon

fig, ax = plt.subplots(figsize=(9, 4.4))
by_h = R.dropna(subset=["err"]).groupby(["model", "h"])["err"].apply(lambda e: np.sqrt((e**2).mean()))
for model in MODELS:
    ax.plot(by_h[model].index, by_h[model].values, "o-", color=MCOL[model], lw=1.8, ms=4, label=model)
ax.set(title="TFR forecast RMSE by horizon (pooled origins)",
       xlabel="years ahead", ylabel="RMSE of TFR"); ax.grid(alpha=.25); ax.legend(fontsize=9)
plt.show()

5. Accuracy and calibration scores

D = R.dropna(subset=["err"])
tab = D.groupby("model").apply(lambda x: pd.Series({
    "MAE": x["err"].abs().mean(), "RMSE": np.sqrt((x["err"]**2).mean()),
    "bias": x["err"].mean(), "cover90": x["covered"].mean()})).round(3)
print("TFR forecast performance, 1971-2024 (all horizons):")
print(tab.to_string())
print("\n bias>0 = over-predicting TFR; cover90 = share of actuals in the 90% band (~0.90 = calibrated).")

TFR forecast performance, 1971-2024 (all horizons):
                     MAE   RMSE   bias  cover90
model                                          
Lee-Tuljapurkar    0.154  0.250  0.130    0.944
Minnesota BVAR     0.254  0.409  0.246    0.759
Steady-state BVAR  0.138  0.213  0.116    0.815

 bias>0 = over-predicting TFR; cover90 = share of actuals in the 90% band (~0.90 = calibrated).

6. Age profile — 1- vs 10-year-ahead

The realized 7-group fertility schedule in 2019 vs the model's forecast of it made 1 year ahead (fit 2018) and 10 years ahead (fit 2009), with 90% bands. The 10-year forecast misses the postponement — younger groups over-predicted, older groups under-predicted — because the age shape is anchored at the origin.

TARGET = 2019
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4), sharey=True)
H1C, H10C = "#2a9d8f", "#e76f51"
for ax, model in zip(axes, MODELS):
    realized = G.loc[TARGET].values
    l1, m1, h1 = forecast_group_profile(model, TARGET - 1, 1)
    l10, m10, h10 = forecast_group_profile(model, TARGET - 10, 10)
    x = range(7)
    ax.plot(x, realized, "k-", lw=2, label=f"realized {TARGET}", zorder=5)
    ax.fill_between(x, l1, h1, color=H1C, alpha=.2); ax.plot(x, m1, "-", color=H1C, lw=1.4, label="1-yr ahead")
    ax.fill_between(x, l10, h10, color=H10C, alpha=.2); ax.plot(x, m10, "--", color=H10C, lw=1.4, label="10-yr ahead")
    ax.set_xticks(x); ax.set_xticklabels(FERT_LABELS, rotation=45, fontsize=7)
    ax.set_title(model, fontsize=10); ax.grid(alpha=.25)
axes[0].set_ylabel("group fertility rate"); axes[0].legend(fontsize=8)
fig.suptitle(f"Fertility age schedule {TARGET}: realized vs 1- and 10-year-ahead (90% bands)", y=1.02)
fig.tight_layout(); plt.show()

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

7. Notes and analysis

Anchoring the long run to $F^*=2$ clearly helps. Both models that pin the ultimate TFR near replacement — the steady-state BVAR (MAE 0.138, best overall) and the Lee-Tuljapurkar benchmark (0.154) — beat the diffuse Minnesota BVAR (0.254) by a wide margin. The Minnesota model reverts toward the 1933-origin sample mean, which the baby boom inflates (see the 1970-origin block shooting up to ~2.9 while the actual TFR was collapsing to 1.7), so it over-predicts most (bias +0.25). This is the mirror image of mortality, where the informative prior hurt: fertility's diffuse long run is badly biased, so constraining it pays off.

Everything still over-predicts. All three biases are positive (+0.12 to +0.25): the $F^*=2$ anchor (and the sample mean) sit above the actual post-2007 decline to ~1.6, so even the best model is systematically too high. $F^*$ is the binding assumption — a lower ultimate TFR would track the realized decline better (as the single-origin $F^*$ sweep showed).

Calibration is reversed from mortality. Here the benchmark is best-calibrated (cover90 0.94): its constrained ARMA(1,1) is highly persistent and throws very wide TFR bands (visible in Section 3) that contain the actual most of the time. The Minnesota BVAR is worst (0.76) — its diffuse-mean over-prediction pushes the actual toward the band edge — with the steady-state between (0.82). Contrast mortality, where it was Lee-Carter whose plug-in bands were far too narrow. The lesson is the same: point accuracy and interval calibration are separate questions, and the ranking on one need not match the other.

Age profile (Section 6). As with mortality, the 1-year-ahead schedule hugs the realized one, while the 10-year-ahead misses the postponement — over-predicting the young groups and under-predicting the older ones — because the age shape is frozen at the origin. This age-pattern error is orthogonal to the TFR-level error and would need genuine per-group dynamics, not just a TFR constraint, to fix.