Life expectancy: published vs. derived from mortality rates¶
This notebook rebuilds the HMD period life table from the age-specific death
rates m_x (US Mx_1x1.txt) and compares the derived life expectancy e_x
against the published HMD values (the ex column of the *ltper_1x1 files),
for females, males, and both sexes combined, at the ages Meseguer (2006, 2010)
reports: birth, 19, 35, 65, 80.
If the life-table code is correct, the derived curves should sit right on top of the published ones.
Data: US series from the Human Mortality Database (mortality.org) and the Human Fertility Database (humanfertility.org).
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from hmd_io import read_mx, read_deaths, read_ltper
from lifetable import life_table
%matplotlib inline
plt.rcParams["figure.dpi"] = 110
1. Load the data¶
DIR = "."
mx = read_mx(f"{DIR}/US Mx_1x1.txt") # dict: 'female'/'male'/'total' -> [Year x Age]
deaths = read_deaths(f"{DIR}/US Deaths_1x1.txt") # for the both-sexes a0 weighting
# Published life tables (their 'ex' column is the published life expectancy).
pub_lt = {
"female": read_ltper(f"{DIR}/US fltper_1x1.txt"),
"male": read_ltper(f"{DIR}/US mltper_1x1.txt"),
"total": read_ltper(f"{DIR}/US bltper_1x1.txt"),
}
years = mx["total"].index.to_numpy()
print(f"Years {years.min()}-{years.max()} ({len(years)} years), ages 0-110+")
Years 1933-2024 (92 years), ages 0-110+
Methodology — how a death rate becomes life expectancy¶
Life expectancy is not measured directly; it is built up from age-specific
death rates through a life table. The input is the central death rate
$m_x$ (deaths at age $x$ per person-year lived at age $x$) from
US Mx_1x1.txt; the output is $e_x$, expected years of life remaining at age
$x$. The chain (implemented in lifetable.py):
1. Rate → probability, $m_x \to q_x$. $m_x$ is a rate; $q_x$ is the probability that someone alive at exact age $x$ dies before $x+1$. People who die were only at risk for part of the year, so
$$q_x = \frac{m_x}{1 + (1-a_x)\,m_x}$$
where $a_x$ is the average fraction of the year lived by those dying in it. Deaths are spread roughly evenly, so $a_x = 0.5$ for most ages. Age 0 is special — infant deaths cluster in the first days of life, so $a_0 \approx 0.05\text{–}0.15$; we set it with the Andreev–Kingkade rule, a function of $m_0$.
2. Survivorship, $l_x$. Start a synthetic cohort of $l_0 = 100{,}000$ newborns and thin it year by year: $\;l_{x+1} = l_x\,(1-q_x)$, with deaths $d_x = l_x - l_{x+1}$.
3. Person-years lived in the interval, $L_x$. Survivors live the full year, those who die live a fraction $a_x$: $\;L_x = l_{x+1} + a_x\,d_x$. Open interval (110+): everyone dies, $L = l/m$.
4. Person-years above age $x$: $\;T_x = \sum_{k\ge x} L_k$ (reverse cumulative sum).
5. Life expectancy: $\;e_x = T_x / l_x$. At $x=0$ this is life expectancy at birth: lay out every life-year the 100,000 cohort will ever live, divide by 100,000.
Unlike fertility rates (Part 2), death rates cannot simply be summed — a rate must first be turned into survivorship before it means anything, which is what the life table does.
2. Derive e_x from the mortality rates¶
For every year and sex we rebuild the life table from m_x and keep the whole
e_x vector (life expectancy at each age).
SEXES = ["female", "male", "total"]
AGES = [0, 19, 35, 65, 80] # Meseguer's reported ages
# derived[sex] -> DataFrame indexed by Year, columns = age, values = derived e_x
derived, published = {}, {}
for sex in SEXES:
d_rows, p_rows = {}, {}
for yr in years:
rates = mx[sex].loc[yr].to_numpy()
if np.isnan(rates).any():
continue
infant = (deaths["female"].loc[yr, 0], deaths["male"].loc[yr, 0]) if sex == "total" else None
ex = life_table(rates, sex, infant_deaths_by_sex=infant)["ex"]
d_rows[yr] = {a: ex[a] for a in AGES}
pub = pub_lt[sex]
p_rows[yr] = {a: pub.loc[(pub.Year == yr) & (pub.Age == a), "ex"].iloc[0] for a in AGES}
derived[sex] = pd.DataFrame(d_rows).T
published[sex] = pd.DataFrame(p_rows).T
derived["total"].head()
| 0 | 19 | 35 | 65 | 80 | |
|---|---|---|---|---|---|
| 1933 | 60.901119 | 48.169212 | 34.609541 | 12.729459 | 5.945596 |
| 1934 | 60.191303 | 47.991807 | 34.427193 | 12.623508 | 5.885641 |
| 1935 | 60.910652 | 48.154327 | 34.557711 | 12.692838 | 5.955585 |
| 1936 | 60.353946 | 47.588278 | 34.008346 | 12.282669 | 5.606858 |
| 1937 | 61.038208 | 48.088404 | 34.398964 | 12.548633 | 5.803086 |
3. Compare in graphs¶
One panel per age. Within each panel: published = solid line, derived = dashed line with markers, one colour per sex. The dashed curves land on top of the solid ones — that's the visual check.
COLORS = {"female": "#d1495b", "male": "#00798c", "total": "#4b4e6d"}
LABELS = {"female": "Female", "male": "Male", "total": "Both combined"}
fig, axes = plt.subplots(len(AGES), 1, figsize=(9, 3.1 * len(AGES)), sharex=True)
for ax, age in zip(axes, AGES):
for sex in SEXES:
c = COLORS[sex]
ax.plot(published[sex].index, published[sex][age], "-", color=c, lw=2.4, alpha=.55,
label=f"{LABELS[sex]} — published")
ax.plot(derived[sex].index, derived[sex][age], "--", color=c, lw=1.1,
marker="o", ms=2.2, markevery=6, label=f"{LABELS[sex]} — derived")
ax.set_title(f"Life expectancy at age {age}" + (" (at birth)" if age == 0 else ""),
fontsize=11, loc="left")
ax.set_ylabel("years")
ax.grid(alpha=.25)
axes[0].legend(ncol=3, fontsize=8, loc="lower right")
axes[-1].set_xlabel("Year")
fig.suptitle("Published vs. derived period life expectancy, US 1933-2024", y=1.001, fontsize=13)
fig.tight_layout()
plt.show()
4. How close is the match?¶
Absolute difference (derived − published), in years, summarised over all
years for each age and sex. Values are ~0.01 yr — the residual is just the
6-decimal rounding of the m_x file.
summ = {}
for sex in SEXES:
err = (derived[sex] - published[sex]).abs()
summ[(LABELS[sex], "mean")] = err.mean()
summ[(LABELS[sex], "max")] = err.max()
tbl = pd.DataFrame(summ)
tbl.index.name = "age"
tbl.round(4)
| Female | Male | Both combined | ||||
|---|---|---|---|---|---|---|
| mean | max | mean | max | mean | max | |
| age | ||||||
| 0 | 0.0060 | 0.0214 | 0.0037 | 0.0177 | 0.0048 | 0.0166 |
| 19 | 0.0062 | 0.0203 | 0.0036 | 0.0142 | 0.0049 | 0.0151 |
| 35 | 0.0064 | 0.0210 | 0.0040 | 0.0145 | 0.0050 | 0.0184 |
| 65 | 0.0080 | 0.0275 | 0.0050 | 0.0182 | 0.0067 | 0.0274 |
| 80 | 0.0143 | 0.0526 | 0.0102 | 0.0525 | 0.0128 | 0.0464 |
# Error over time, at birth — shows there is no systematic drift.
fig, ax = plt.subplots(figsize=(9, 3))
for sex in SEXES:
ax.plot(years[:len(derived[sex])], (derived[sex][0] - published[sex][0]),
color=COLORS[sex], lw=1.3, label=LABELS[sex])
ax.axhline(0, color="k", lw=.6)
ax.set(title="Derived − published e0 (years)", xlabel="Year", ylabel="years")
ax.grid(alpha=.25); ax.legend(fontsize=8)
plt.show()
Part 2 — Fertility (Human Fertility Database)¶
The Human Mortality Database has no fertility rates (US Births.txt is only
births by the baby's sex). The Total Fertility Rate comes from HMD's sister
database, the Human Fertility Database (HFD), whose US series also runs
1933-2024.
We check the same thing we did for mortality: derive the TFR from the age-specific rates and compare it to the published value.
USAasfrTR.txt- period ASFR by Lexis triangle (Year x Age x cohort).USAtfrRR.txt- HFD's published period TFR (our validation target).
Triangle catch: each age-year square is split into two triangles, each on ~half a year of exposure, so the two rows sum to ~2x the annual rate -> average them (x0.5). TFR = 0.5 * sum of all triangle rates in the year.
Methodology — ASFR and TFR¶
ASFR (age-specific fertility rate) at age $x$ = births to women aged $x$ divided by the number of women aged $x$, in a given year. Units: births per woman per year. Plotted against age it is the fertility schedule (section 2.2): near zero in the early teens, peaking in the 20s–30s, back to zero by ~50.
TFR (total fertility rate) is the area under that schedule — the sum of the age-specific rates over all ages:
$$\text{TFR} = \sum_x \text{ASFR}(x)$$
Because each ASFR is births per woman per year, summing across all childbearing ages gives births per woman over a lifetime. It is a period synthetic-cohort measure: the number of children a woman would have if she lived through all her childbearing years experiencing this year's rates at each age (a snapshot of one year, not a real cohort's completed family size). Replacement level is ≈ 2.1.
Lexis-triangle detail of our file. USAasfrTR.txt splits each age-year cell
into two triangles by the mother's birth cohort, and each triangle's rate is on
~half a year of exposure — so the two rows sum to ≈ 2× the annual rate. We
average them (the $0.5\times$ in tfr_from_triangles) to recover the
conventional one-year ASFR, then sum over ages for the TFR.
from hmd_io import asfr_by_age, tfr_from_triangles, read_tfr
asfr = asfr_by_age(f"{DIR}/USAasfrTR.txt") # annual ASFR, Year x Age
tfr_derived = tfr_from_triangles(f"{DIR}/USAasfrTR.txt")
tfr_published = read_tfr(f"{DIR}/USAtfrRR.txt")["TFR"]
fert = pd.DataFrame({"derived": tfr_derived, "published": tfr_published}).dropna()
fert["abs_err"] = (fert["derived"] - fert["published"]).abs()
print(f"TFR years {fert.index.min()}-{fert.index.max()} (n={len(fert)})")
print(f"mean abs err = {fert['abs_err'].mean():.5f} max = {fert['abs_err'].max():.5f}")
fert.head()
TFR years 1933-2024 (n=92) mean abs err = 0.00039 max = 0.00124
| derived | published | abs_err | |
|---|---|---|---|
| Year | |||
| 1933 | 2.011440 | 2.012 | 0.000560 |
| 1934 | 2.071745 | 2.072 | 0.000255 |
| 1935 | 2.037715 | 2.038 | 0.000285 |
| 1936 | 2.006750 | 2.007 | 0.000250 |
| 1937 | 2.038320 | 2.038 | 0.000320 |
2.1 Total Fertility Rate: derived vs. published¶
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 5.4), sharex=True,
gridspec_kw={"height_ratios": [3, 1]})
ax1.plot(fert.index, fert["published"], "-", color="#00798c", lw=2.4, alpha=.55,
label="published (USAtfrRR.txt)")
ax1.plot(fert.index, fert["derived"], "--", color="#00798c", lw=1.1,
marker="o", ms=2.4, markevery=4, label="derived (triangles x0.5)")
ax1.axhline(2.1, color="gray", ls=":", lw=1, label="replacement ~2.1")
ax1.annotate("baby boom", (1957, fert.loc[1957, "published"]),
xytext=(1962, 3.4), fontsize=8, color="#444")
ax1.set(ylabel="TFR (births per woman)",
title="US Total Fertility Rate, 1933-2024 - derived vs published")
ax1.grid(alpha=.25); ax1.legend(fontsize=8)
ax2.plot(fert.index, (fert["derived"] - fert["published"]) * 1000, color="#d1495b", lw=1)
ax2.axhline(0, color="k", lw=.6)
ax2.set(xlabel="Year", ylabel="err (x10^-3)", title="derived - published")
ax2.grid(alpha=.25)
fig.tight_layout(); plt.show()
2.2 The fertility schedule (ASFR age profile)¶
Age-specific fertility rates for selected years. The schedule falls sharply after the baby boom and shifts to older ages over time (childbearing postponement).
fig, ax = plt.subplots(figsize=(9, 4))
cmap = plt.cm.viridis
show_years = [1933, 1957, 1975, 1990, 2005, 2024]
for i, yr in enumerate(show_years):
ax.plot(asfr.columns, asfr.loc[yr], color=cmap(i / (len(show_years) - 1)),
lw=1.8, label=str(yr))
ax.set(xlabel="Age of mother", ylabel="ASFR (births per woman per year)",
title="US age-specific fertility schedule, selected years", xlim=(12, 50))
ax.grid(alpha=.25); ax.legend(title="Year", fontsize=8, ncol=2)
plt.show()
# TFR is the area under each schedule:
print("TFR (area under schedule) check:")
for yr in show_years:
print(f" {yr}: {asfr.loc[yr].sum():.3f}")
TFR (area under schedule) check: 1933: 2.011 1957: 3.739 1975: 1.768 1990: 2.067 2005: 2.061 2024: 1.593
Part 3 — Population by age & net migration¶
Population by age is built with the cohort-component identity (Meseguer 2010, eq. 22). Using Jan-1 population $P$, calendar-year deaths $D$, births $B$ and net migration $M$:
$$P_{t+1}(0)=B_t-D_t(0)+M_t(0),\qquad P_{t+1}(a{+}1)=P_t(a)-D_t(a)+M_t(a).$$
Net migration is not published and not modelled by the BVAR. Following
Meseguer we treat it as exogenous; here we derive the historical schedule as the
residual that closes the identity (net_migration in population.py):
$$M_t(a)=P_{t+1}(a)-\underbrace{[\text{survivors if }M=0]}_{\text{closed}}.$$
The residual = true net migration + census/estimation error + a small Lexis timing
approximation (HMD deaths are calendar-year squares, so the cohort that ages one
year loses ~half of two adjacent death squares — the 0.5 split in the code).
Territorial break: Jan-1 1959 appears twice in US Population.txt — 1959-
(contiguous US) and 1959+ (incl. Alaska & Hawaii). The reconstruction bridges
this so the ~0.8M reclassification is not mistaken for migration.
from hmd_io import read_population, read_births
from population import net_migration
pop = read_population(f"{DIR}/US Population.txt") # Jan-1, by Year x Age x sex
births = read_births(f"{DIR}/US Births.txt")
# deaths already loaded in Part 1 as `deaths`
mig = {s: net_migration(pop, deaths, births, s)[0] for s in ["female", "male", "total"]}
net_tot = mig["total"].sum(axis=1) # total net migration per year
print(f"Net migration derived for {net_tot.index.min()}-{net_tot.index.max()}")
print(f"mean 2000-2019: {net_tot.loc[2000:2019].mean():,.0f} / yr "
f"(known US net international migration ~1.0-1.1M/yr)")
for y in [1935, 1970, 2000, 2019, 2020, 2023]:
print(f" {y}: {net_tot[y]:>12,.0f}")
Net migration derived for 1933-2024
mean 2000-2019: 1,061,376 / yr (known US net international migration ~1.0-1.1M/yr) 1935: -136,094 1970: 874,740 2000: 1,309,820 2019: 827,715 2020: 770,234 2023: 2,524,732
3.1 Total net migration over time¶
The derived residual tracks the known US history: near-zero/negative in the 1930s Depression, rising after the 1965 Immigration Act, ~1M/yr in the 2000s, the 2020 COVID dip, and the post-2021 surge.
fig, ax = plt.subplots(figsize=(9, 4))
ax.fill_between(net_tot.index, net_tot / 1e6, 0, color="#00798c", alpha=.25)
ax.plot(net_tot.index, net_tot / 1e6, color="#00798c", lw=1.6)
ax.axhline(0, color="k", lw=.6)
for yr, txt, dy in [(1933, "Depression", .3), (1970, "post-1965 Act", .35),
(2020, "COVID dip", -.5), (2023, "2021+ surge", .2)]:
ax.annotate(txt, (yr, net_tot.get(yr, 0) / 1e6), fontsize=8, color="#444",
xytext=(yr, net_tot.get(yr, 0) / 1e6 + dy))
ax.set(title="US net migration (derived residual), 1933-2024",
xlabel="Year", ylabel="net migrants (millions/yr)")
ax.grid(alpha=.25)
plt.show()
3.2 Net migration age profile¶
Averaged over recent decades and grouped into 5-year bands to suppress the single-year census-interpolation noise. Migration concentrates at children and young adults and turns negative at the oldest ages — the classic immigration shape. (Meseguer's exogenous assumption schedules are smooth versions of this.)
def grouped_profile(m, years, width=5, max_age=84):
prof = m.loc[years[0]:years[1], [a for a in m.columns if a <= max_age]].mean(axis=0)
g = prof.groupby((prof.index // width) * width).sum()
return g
fig, ax = plt.subplots(figsize=(9, 4))
for (lo, hi), c in [((1950, 1969), "#9aa0b0"), ((2000, 2019), "#00798c")]:
g = grouped_profile(mig["total"], (lo, hi)) / 1e3
ax.bar(g.index + (2.2 if lo == 2000 else -0.2), g.values, width=2.2,
color=c, label=f"{lo}-{hi} avg")
ax.axhline(0, color="k", lw=.6)
ax.set(title="Net migration by age (5-yr groups, period averages)",
xlabel="Age", ylabel="net migrants (thousands/yr)")
ax.grid(alpha=.25, axis="y"); ax.legend(fontsize=8)
plt.show()
3.3 The two engines of population growth¶
Total population change decomposes exactly into natural increase (births − deaths, driven by Part 1 & 2) and net migration (Part 3). This closes the loop across all three vital processes.
B = births["Total"]
Dt = deaths["total"].sum(axis=1) # total deaths per year
nat = (B - Dt) # natural increase
comp = pd.DataFrame({"natural increase": nat, "net migration": net_tot}).dropna()
fig, ax = plt.subplots(figsize=(9, 4))
ax.stackplot(comp.index, comp["natural increase"] / 1e6, comp["net migration"] / 1e6,
labels=["natural increase (births - deaths)", "net migration"],
colors=["#d1495b", "#00798c"], alpha=.8)
ax.plot(comp.index, comp.sum(axis=1) / 1e6, color="k", lw=1, label="total change")
ax.axhline(0, color="k", lw=.5)
ax.set(title="Components of US population change, 1933-2024",
xlabel="Year", ylabel="millions / yr")
ax.grid(alpha=.25); ax.legend(fontsize=8, loc="upper left")
plt.show()
print("Share of population growth from net migration, by decade:")
for d in range(1960, 2020, 10):
sl = comp.loc[d:d+9]
print(f" {d}s: {100*sl['net migration'].sum()/sl.sum(axis=1).sum():4.0f}%")
Share of population growth from net migration, by decade: 1960s: 14% 1970s: 39% 1980s: 24% 1990s: 46% 2000s: 38% 2010s: 47%
Part 4 — Population pyramids: derived vs. published¶
Can we reproduce the population by age? With net migration taken as the residual
(Part 3) the cohort-component identity closes by construction, so that check is
trivial. The honest test is a closed reconstruction (project_closed): take the
observed pyramid in a base year and roll it forward with only survival (mortality)
and observed births — no migration. The gap between this migration-free pyramid
and the published one is exactly the accumulated net migration.
We use base year 1933 — the first year of data, so every later pyramid can be derived. (Rolling across the 1959 Alaska/Hawaii territorial break adds a small one-off ~0.8M to the post-1959 gap, on top of true migration.)
from population import project_closed
BASE = 1933
closed_e = {s: project_closed(pop, deaths, births, s, base_year=BASE) for s in ["male", "female"]}
YEARS = [1940, 1960, 1980, 2000, 2020]
print("year published derived(no mig.) gap = cumulative net migration since 1933")
for yr in YEARS:
pub = sum(pop[s].loc[str(yr)].sum() for s in ["male", "female"])
der = sum(closed_e[s].loc[yr].sum() for s in ["male", "female"])
print(f"{yr} {pub/1e6:7.1f}M {der/1e6:7.1f}M {(pub-der)/1e6:6.1f}M")
year published derived(no mig.) gap = cumulative net migration since 1933 1940 131.5M 132.2M -0.8M 1960 178.6M 175.0M 3.6M 1980 225.8M 210.5M 15.3M 2000 280.1M 246.2M 33.9M 2020 330.9M 278.8M 52.1M
4.1 Evolution 1940 → 2020: published (bars) vs. migration-free (outline)¶
Watch the age structure change — the Depression-pinched base in 1940, the baby-boom-wide base in 1960, and that bulge climbing the pyramid through 1980, 2000, 2020 as it ages. The outline (survival + births, no migration) pulls progressively inside the published bars: in 1940 the two nearly coincide (little migration yet), and by 2020 the gap is the ~52M of accumulated net migration.
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
CM, CF = "#00798c", "#d1495b"
ages = np.arange(0, 101)
maxv = max(pop[s].loc[str(y)].reindex(ages).max() for s in ["male", "female"] for y in YEARS) / 1e6
fig, axes = plt.subplots(1, len(YEARS), figsize=(15, 5), sharey=True)
for ax, yr in zip(axes, YEARS):
pm = pop["male"].loc[str(yr)].reindex(ages) / 1e6
pf = pop["female"].loc[str(yr)].reindex(ages) / 1e6
dm = closed_e["male"].loc[yr].reindex(ages) / 1e6
df_ = closed_e["female"].loc[yr].reindex(ages) / 1e6
ax.barh(ages, -pm, height=1.0, color=CM, alpha=.30)
ax.barh(ages, pf, height=1.0, color=CF, alpha=.30)
ax.step(-dm, ages, where="mid", color=CM, lw=1.0)
ax.step(df_, ages, where="mid", color=CF, lw=1.0)
ax.axvline(0, color="k", lw=.5)
tot = (pop["male"].loc[str(yr)].sum() + pop["female"].loc[str(yr)].sum()) / 1e6
ax.set_title(f"{yr} ({tot:.0f}M)", fontsize=10)
ax.set_xlim(-maxv * 1.05, maxv * 1.05)
ax.set_xticks([-2, -1, 0, 1, 2]); ax.set_xticklabels(["2", "1", "0", "1", "2"], fontsize=7)
ax.grid(alpha=.2, axis="x")
axes[0].set_ylabel("Age")
handles = [Patch(facecolor=CM, alpha=.3, label="published (male)"),
Patch(facecolor=CF, alpha=.3, label="published (female)"),
Line2D([0], [0], color=CM, label="derived, no migration (male)"),
Line2D([0], [0], color=CF, label="derived, no migration (female)")]
fig.legend(handles=handles, loc="upper center", ncol=4, fontsize=8, framealpha=.9,
bbox_to_anchor=(0.5, 1.06))
fig.supxlabel("Population (millions) — male left / female right", y=-0.02)
fig.suptitle("US population pyramids, base 1933 rolled forward with no migration", y=1.12)
fig.tight_layout()
plt.show()
4.2 Survival validated at old ages (short horizon)¶
Over a long roll, migration accumulates at every age, so a big gap is expected everywhere. To isolate the survival machinery, roll only five years (2015 → 2020): almost no one migrates into ages 75+ over five years, so the migration-free reconstruction should land on the published pyramid — and it does, within ~2%.
YEAR, SBASE = 2020, 2015
closed_s = {s: project_closed(pop, deaths, births, s, base_year=SBASE) for s in ["male", "female"]}
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
hi = np.arange(70, 101)
for ax, s, c in zip(axes, ["male", "female"], [CM, CF]):
ax.plot(hi, pop[s].loc[str(YEAR)].reindex(hi) / 1e3, "-", color=c, lw=2.4, alpha=.55, label="published")
ax.plot(hi, closed_s[s].loc[YEAR].reindex(hi) / 1e3, "--", color="k", lw=1.1,
label=f"derived, no mig. ({SBASE}->{YEAR})")
ax.set(title=f"{s.capitalize()}, ages 70+", xlabel="Age", ylabel="thousands")
ax.grid(alpha=.25); ax.legend(fontsize=8)
fig.suptitle(f"Survival check, {SBASE}->{YEAR} (5-year horizon)", y=1.03)
fig.tight_layout(); plt.show()
h = np.arange(75, 101)
for s in ["male", "female"]:
d = np.abs(pop[s].loc[str(YEAR)].reindex(h) - closed_s[s].loc[YEAR].reindex(h)).sum()
p = pop[s].loc[str(YEAR)].reindex(h).sum()
print(f"{s}: ages 75-100 abs diff {d/1e6:.2f}M = {100*d/p:.1f}% of {p/1e6:.1f}M")
male: ages 75-100 abs diff 0.20M = 2.2% of 9.2M female: ages 75-100 abs diff 0.27M = 2.1% of 12.8M