xjtulyc/demographic-analysis

Demographic Analysis

When to Use

Use this skill when you need to:

• Construct period or cohort life tables from age-specific death rates (mx) following standard actuarial conventions
• Compute life expectancy at birth (e0) with Chiang confidence intervals
• Load and process HMD (Human Mortality Database) data files for multiple countries and years
• Build a Leslie matrix for a female population and project population structure 50 years forward
• Derive stable population parameters from the Leslie matrix dominant eigenvalue
• Estimate the Lee-Carter mortality model (log m_{x,t} = a_x + b_x × k_t) via SVD and project future mortality rates with ARIMA forecasts of the time index k_t
• Decompose the life expectancy gap between two populations using Arriaga's method (contribution by age group to the e0 difference)

Background

Life Table Notation:

| Symbol | Formula | Meaning | |---|---|---| | m_x | observed | Death rate in age interval [x, x+5) | | q_x | 2m_x / (2 + m_x) | Probability of dying in interval | | l_x | l_{x-1} × (1 - q_{x-1}) | Survivors to exact age x (radix l_0=100,000) | | d_x | l_x × q_x | Deaths in interval | | L_x | (l_x + l_{x+1}) / 2 × n | Person-years lived in interval | | T_x | Σ_{a≥x} L_a | Total person-years above age x | | e_x | T_x / l_x | Remaining life expectancy at age x |

For the open-ended age interval (e.g., 85+): L_{85+} = l_{85} / m_{85}.

Chiang (1984) CI on e0: Based on variance of d_x:

Var(e0) ≈ Σ_x (l_x/l_0)² × (q_x(1-q_x) / d_x) × (e_x + (1-a_x) × n)²

where a_x is the average fraction of interval lived by those dying in it (typically 0.5, except infants where a_0 ≈ 0.1).

Leslie Matrix: A structured population projection matrix where:

• Top row = age-specific fertility rates (F_x = l_{x+1}/l_x × m_x, where m_x is maternity rate)
• Sub-diagonal = survival proportions (P_x = L_{x+1} / L_x)

Population projection: N(t+5) = A × N(t). After many iterations, the population grows at rate λ (dominant eigenvalue) and converges to the stable age distribution (dominant eigenvector).

Lee-Carter Model: The dominant model for mortality forecasting since 1992:

log(m_{x,t}) = a_x + b_x × k_t + ε_{x,t}

• a_x: average age profile of log mortality (row means)
• b_x: sensitivity of age x to the time trend k_t
• k_t: time index capturing overall mortality improvement

Estimated via SVD of the residual matrix after subtracting a_x. k_t is then modeled as a random walk with drift and projected forward using ARIMA(0,1,0) with drift.

Arriaga Decomposition: Decomposes the difference e0(pop1) - e0(pop2) into contributions from each age group. Useful for understanding which ages contribute most to a mortality gap.

Environment Setup

pip install pandas>=1.5 numpy>=1.23 scipy>=1.9 matplotlib>=3.6

HMD data: Register at https://www.mortality.org/ and download country life table files. Typical HMD file structure: columns Year, Age, mx, qx, ax, lx, dx, Lx, Tx, ex.

export HMD_DIR="/data/hmd/"

Core Workflow

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings("ignore")


# ---------------------------------------------------------------------------
# 1. Life Table Construction
# ---------------------------------------------------------------------------

def construct_life_table(
    mx: pd.Series,
    age_groups: list[int],
    radix: int = 100_000,
    last_open: bool = True,
) -> pd.DataFrame:
    """
    Construct a standard period life table from age-specific death rates.

    Parameters
    ----------
    mx : pd.Series
        Age-specific mortality rates (deaths per person-year) indexed by age_groups.
    age_groups : list of int
        Age group starting values (e.g., [0,1,5,10,...,85]).
    radix : int
        Starting cohort size (default 100,000).
    last_open : bool
        Whether the last age group is an open-ended interval.

    Returns
    -------
    pd.DataFrame with columns: age, mx, qx, ax, lx, dx, Lx, Tx, ex.
    """
    n_ages = len(age_groups)
    # Width of each age interval
    widths = [age_groups[i+1] - age_groups[i] for i in range(n_ages - 1)] + [np.inf]

    mx_vals = mx.values.astype(float)
    ax_vals = np.full(n_ages, 0.5)  # average fraction of interval lived
    ax_vals[0] = 0.1  # infant (WHO convention: ~0.1 for infants)

    qx = np.zeros(n_ages)
    for i, (m, a, w) in enumerate(zip(mx_vals, ax_vals, widths)):
        if np.isinf(w):
            qx[i] = 1.0  # open-ended interval
        else:
            qx[i] = (w * m) / (1 + (1 - a) * w * m)

    lx = np.zeros(n_ages)
    lx[0] = radix
    for i in range(1, n_ages):
        lx[i] = lx[i-1] * (1 - qx[i-1])

    dx = lx * qx

    Lx = np.zeros(n_ages)
    for i, (l, d, a, w) in enumerate(zip(lx, dx, ax_vals, widths)):
        if np.isinf(w):
            Lx[i] = lx[i] / mx_vals[i] if mx_vals[i] > 0 else lx[i] * 5
        else:
            Lx[i] = w * (lx[i] - d) + a * w * d

    Tx = np.cumsum(Lx[::-1])[::-1]
    ex = Tx / np.where(lx > 0, lx, np.nan)

    return pd.DataFrame({
        "age": age_groups,
        "width": [w if not np.isinf(w) else None for w in widths],
        "mx": mx_vals.round(6),
        "qx": qx.round(6),
        "ax": ax_vals,
        "lx": lx.round(1),
        "dx": dx.round(1),
        "Lx": Lx.round(1),
        "Tx": Tx.round(1),
        "ex": ex.round(3),
    })


def life_expectancy_ci(lt: pd.DataFrame, alpha: float = 0.05) -> dict:
    """
    Compute Chiang (1984) confidence interval for life expectancy at birth (e0).

    Parameters
    ----------
    lt : pd.DataFrame
        Life table output from construct_life_table().
    alpha : float
        Significance level.

    Returns
    -------
    dict with e0, ci_lower, ci_upper, se.
    """
    e0 = lt.loc[0, "ex"]
    radix = lt.loc[0, "lx"]

    # Variance contribution from each age group (Chiang formula)
    var_contributions = []
    for _, row in lt.iterrows():
        if row["dx"] > 0 and row["qx"] < 1.0:
            term = (row["lx"] / radix) ** 2 * row["qx"] * (1 - row["qx"]) / row["dx"]
            var_contributions.append(term * row["ex"] ** 2)
    var_e0 = sum(var_contributions)
    se = np.sqrt(var_e0)

    z = stats.norm.ppf(1 - alpha / 2)
    return {
        "e0": round(e0, 3),
        "ci_lower": round(e0 - z * se, 3),
        "ci_upper": round(e0 + z * se, 3),
        "se": round(se, 4),
    }


# ---------------------------------------------------------------------------
# 2. Leslie Matrix Population Projection
# ---------------------------------------------------------------------------

def build_leslie_matrix(
    Lx: np.ndarray,
    lx: np.ndarray,
    Fx: np.ndarray,
) -> np.ndarray:
    """
    Construct a Leslie projection matrix for female population.

    Parameters
    ----------
    Lx : np.ndarray
        Person-years lived in each age interval (from life table).
    lx : np.ndarray
        Survivors at start of each age interval.
    Fx : np.ndarray
        Age-specific fertility rates (births per woman per 5-year period).
        Must be same length as Lx.

    Returns
    -------
    np.ndarray — Leslie matrix (n_ages × n_ages).
    """
    n = len(Lx)
    A = np.zeros((n, n))

    # Top row: fertility (accounting for infant survival to next census)
    A[0, :] = Fx * Lx / lx

    # Sub-diagonal: survival proportions
    for i in range(n - 1):
        if Lx[i] > 0:
            A[i + 1, i] = Lx[i + 1] / Lx[i]

    return A


def project_population(
    A: np.ndarray,
    N0: np.ndarray,
    n_steps: int = 10,
    step_years: int = 5,
) -> pd.DataFrame:
    """
    Project population n_steps periods into the future using the Leslie matrix.

    Parameters
    ----------
    A : np.ndarray
        Leslie matrix (n_ages × n_ages).
    N0 : np.ndarray
        Initial population vector by age group.
    n_steps : int
        Number of projection steps.
    step_years : int
        Years per step (typically 5).

    Returns
    -------
    pd.DataFrame with columns: step, year_offset, total_pop, growth_rate,
    plus age-group columns.
    """
    N = N0.copy().astype(float)
    rows = [{"step": 0, "year_offset": 0, "total_pop": N.sum(), **{f"age_{i}": N[i] for i in range(len(N))}}]

    prev_total = N.sum()
    for step in range(1, n_steps + 1):
        N = A @ N
        total = N.sum()
        row = {
            "step": step,
            "year_offset": step * step_years,
            "total_pop": total,
            "growth_rate": (total / prev_total) ** (1 / step_years) - 1,
        }
        row.update({f"age_{i}": N[i] for i in range(len(N))})
        rows.append(row)
        prev_total = total

    return pd.DataFrame(rows)


def stable_population_params(A: np.ndarray) -> dict:
    """
    Derive stable population parameters from the Leslie matrix dominant eigenvalue.

    Returns
    -------
    dict with lambda_1 (growth multiplier per step), intrinsic_r (per year),
    stable_age_structure (proportions), NRR (Net Reproduction Rate, trace-based approximation).
    """
    eigenvalues, eigenvectors = np.linalg.eig(A)
    dominant_idx = np.argmax(np.real(eigenvalues))
    lambda_1 = float(np.real(eigenvalues[dominant_idx]))
    r = np.log(lambda_1) / 5  # assuming 5-year step

    stable_vec = np.real(eigenvectors[:, dominant_idx])
    stable_vec = np.abs(stable_vec) / np.abs(stable_vec).sum()

    return {
        "lambda_1": round(lambda_1, 6),
        "intrinsic_r": round(r, 6),
        "doubling_time_years": round(np.log(2) / r, 1) if r > 0 else np.inf,
        "stable_age_structure": stable_vec.round(5),
    }


# ---------------------------------------------------------------------------
# 3. Lee-Carter Mortality Model
# ---------------------------------------------------------------------------

def fit_lee_carter(
    mx_matrix: pd.DataFrame,
) -> dict:
    """
    Fit the Lee-Carter model via SVD.

    Parameters
    ----------
    mx_matrix : pd.DataFrame
        Matrix of log(m_{x,t}) with age groups as index, years as columns.
        NaN values are replaced with column means.

    Returns
    -------
    dict with a_x (pd.Series), b_x (pd.Series), k_t (pd.Series), svd_singular_values.
    """
    log_mx = np.log(mx_matrix.replace(0, np.nan))
    # a_x: mean over time
    a_x = log_mx.mean(axis=1)
    # Residual matrix
    residual = log_mx.subtract(a_x, axis=0).fillna(0)

    # SVD: keep first singular vector
    U, s, Vt = np.linalg.svd(residual.values, full_matrices=False)
    b_x = pd.Series(U[:, 0], index=mx_matrix.index)
    k_t_raw = s[0] * Vt[0, :]

    # Normalize: b_x sums to 1, k_t adjusted accordingly
    b_norm = b_x / b_x.sum()
    k_t = pd.Series(k_t_raw * b_x.sum(), index=mx_matrix.columns)

    # Re-center: shift k_t so that Σ k_t = 0 and adjust a_x
    k_mean = k_t.mean()
    k_t_centered = k_t - k_mean
    a_x_adjusted = a_x + b_norm * k_mean

    return {
        "a_x": a_x_adjusted,
        "b_x": b_norm,
        "k_t": k_t_centered,
        "singular_values": s[:3],
    }


def forecast_lee_carter(
    lc: dict,
    n_years: int = 20,
    age_groups: list[int] | None = None,
) -> pd.DataFrame:
    """
    Project Lee-Carter mortality rates using a random walk with drift for k_t.

    Parameters
    ----------
    lc : dict
        Output from fit_lee_carter().
    n_years : int
        Number of years to forecast.
    age_groups : list of int, optional
        Age groups for the result index.

    Returns
    -------
    pd.DataFrame — projected m_{x,t} for each forecast year, age groups as index.
    """
    k_t = lc["k_t"]
    last_year = k_t.index[-1]

    # Fit random walk with drift: k_{t+1} = k_t + c + e
    diffs = k_t.diff().dropna()
    drift = diffs.mean()
    se_drift = diffs.std()

    forecast_years = [last_year + i for i in range(1, n_years + 1)]
    k_last = float(k_t.iloc[-1])
    k_forecasts = [k_last + drift * (i + 1) for i in range(n_years)]

    a_x = lc["a_x"]
    b_x = lc["b_x"]
    projected = {}
    for yr, k_proj in zip(forecast_years, k_forecasts):
        log_mx_proj = a_x + b_x * k_proj
        projected[yr] = np.exp(log_mx_proj)

    return pd.DataFrame(projected, index=a_x.index)

Advanced Usage

Arriaga Decomposition of e0 Gap

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


def arriaga_decomposition(
    lt1: pd.DataFrame,
    lt2: pd.DataFrame,
) -> pd.DataFrame:
    """
    Decompose the difference e0(lt1) - e0(lt2) into age-group contributions
    using Arriaga's (1984) method.

    Parameters
    ----------
    lt1 : pd.DataFrame
        Life table for population 1 (output of construct_life_table).
    lt2 : pd.DataFrame
        Life table for population 2.

    Returns
    -------
    pd.DataFrame with age, direct_effect, indirect_effect, total_contribution.
    """
    radix = lt1.loc[0, "lx"]

    def direct(x, lt_i, lt_j):
        l_i = lt_i.loc[x, "lx"] / radix
        L_i = lt_i.loc[x, "Lx"] / radix
        L_j = lt_j.loc[x, "Lx"] / radix
        T_j = lt_j.loc[x, "Tx"] / radix
        return l_i * (L_i / lt_j.loc[x, "lx"] * radix - L_j)

    rows = []
    for idx in lt1.index[:-1]:  # Exclude last open-ended group
        l_i = lt1.loc[idx, "lx"] / radix
        T_j_next = lt2.loc[idx + 1, "Tx"] / radix if idx + 1 in lt2.index else 0
        T_i = lt1.loc[idx, "Tx"] / radix
        T_j = lt2.loc[idx, "Tx"] / radix
        L_i = lt1.loc[idx, "Lx"] / radix
        L_j = lt2.loc[idx, "Lx"] / radix
        lj_x = lt2.loc[idx, "lx"]
        if lj_x <= 0:
            continue
        direct_e = l_i * (L_i / lj_x * radix - L_j)
        indirect_e = T_j_next / radix * (lt1.loc[idx, "lx"] / radix - lt2.loc[idx, "lx"] / radix
                                          ) if T_j_next > 0 else 0
        rows.append({
            "age": lt1.loc[idx, "age"],
            "direct_effect": round(direct_e, 5),
            "indirect_effect": round(indirect_e, 5),
            "total_contribution": round(direct_e + indirect_e, 5),
        })
    df = pd.DataFrame(rows)
    df["cumulative"] = df["total_contribution"].cumsum()
    return df

Troubleshooting

| Problem | Cause | Solution | |---|---|---| | Negative qx values | mx too high combined with ax near 0 | Clamp qx to [0,1] after calculation | | Leslie matrix eigenvalue < 1 | Declining population (below-replacement fertility) | Expected for many developed countries; report r < 0 | | Lee-Carter SVD gives wrong sign | SVD sign indeterminacy | Flip sign of b_x and k_t together if k_t is increasing (mortality rising) | | e0 Chiang CI too narrow | Large dx values (large sample) | This is correct behavior; CI shrinks with more exposure | | HMD file format varies | Different versions use different delimiters | Use pd.read_table(..., sep=r'\s+', skiprows=2) for whitespace-delimited HMD files |

External Resources

• Human Mortality Database: https://www.mortality.org/
• Preston, S., Heuveline, P. & Guillot, M. (2001). Demography: Measuring and Modeling Population Processes. Wiley-Blackwell.
• Lee, R. & Carter, L. (1992). Modeling and forecasting U.S. mortality. JASA, 87(419), 659-671.
• Arriaga, E. (1984). Measuring and explaining the change in life expectancies. Demography, 21(1), 83-96.
• HMD User Manual: https://www.mortality.org/File/GetDocument/Public/Docs/MethodsProtocolV6.pdf

Examples

Example 1: Life Table from mx + e0 with CI

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Illustrative mx values (approximate US female mortality, 2019)
age_groups = [0, 1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85]
mx_values = [
    0.00558, 0.00027, 0.00013, 0.00011, 0.00048, 0.00060, 0.00064, 0.00080,
    0.00111, 0.00161, 0.00253, 0.00388, 0.00609, 0.00940, 0.01408, 0.02214,
    0.03567, 0.05959, 0.15200,
]
mx_series = pd.Series(mx_values, index=age_groups)
lt = construct_life_table(mx_series, age_groups)

print("=== Life Table (selected ages) ===")
print(lt[["age", "mx", "qx", "lx", "Lx", "ex"]].to_string(index=False))

ci = life_expectancy_ci(lt)
print(f"\ne0 = {ci['e0']} years [{ci['ci_lower']}, {ci['ci_upper']}] (95% CI)")

# Plot survival curve
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(lt["age"], lt["lx"] / 100_000, "b-", linewidth=2)
ax.fill_between(lt["age"], 0, lt["lx"] / 100_000, alpha=0.15, color="blue")
ax.set_xlabel("Age")
ax.set_ylabel("Probability of Survival (l_x / l_0)")
ax.set_title(f"Survival Curve — e0 = {ci['e0']} years")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("survival_curve.png", dpi=150)
plt.show()

Example 2: Leslie Matrix 50-Year Population Projection

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simplified female population (5-year age groups, 0-85+)
n_groups = len(age_groups)
lx_vals = lt["lx"].values
Lx_vals = lt["Lx"].values

# Stylized fertility: ages 20-40 fertile (TFR ≈ 1.8)
Fx = np.zeros(n_groups)
fertile_ages = [4, 5, 6, 7, 8]  # indices for 20-24, 25-29, 30-34, 35-39, 40-44
for idx in fertile_ages:
    Fx[idx] = 0.18 / 5  # births per woman per 5-year interval
# Adjust for sex ratio at birth and infant survival
Fx *= 0.4878  # fraction female × sex ratio

A = build_leslie_matrix(Lx_vals, lx_vals, Fx)

# Initial population (approximately stable distribution)
N0 = np.array([6000, 5800, 5600, 5500, 5400, 5200, 5100, 4900,
               4700, 4400, 4100, 3700, 3200, 2700, 2100, 1600, 1100, 700, 400]).astype(float)

proj = project_population(A, N0, n_steps=10, step_years=5)

print("=== 50-Year Population Projection ===")
print(proj[["step", "year_offset", "total_pop", "growth_rate"]].round(4).to_string(index=False))

params = stable_population_params(A)
print(f"\nStable population parameters:")
print(f"  Lambda (5-year): {params['lambda_1']}")
print(f"  Intrinsic rate r: {params['intrinsic_r']:.5f} per year")
print(f"  Doubling time: {params['doubling_time_years']} years" if params['lambda_1'] > 1
      else f"  Halving time: {abs(round(np.log(2) / abs(params['intrinsic_r']), 1))} years")

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(proj["year_offset"], proj["total_pop"] / 1000, "o-", color="#2166ac", linewidth=2)
ax.set_xlabel("Years from baseline")
ax.set_ylabel("Female Population (thousands)")
ax.set_title("Leslie Matrix: 50-Year Population Projection")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("leslie_projection.png", dpi=150)
plt.show()

Example 3: Lee-Carter SVD Estimation + 20-Year Forecast

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
years_lc = list(range(1960, 2021))
age_idx = [0, 1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85]
base_mx = np.array(mx_values)

# Simulate mortality decline over time
mx_sim = {}
for t, yr in enumerate(years_lc):
    decline = 0.015 * t  # linear mortality improvement
    mx_sim[yr] = base_mx * np.exp(-decline + rng.normal(0, 0.03, len(age_idx)))
mx_df = pd.DataFrame(mx_sim, index=age_idx)

lc = fit_lee_carter(mx_df)
print("=== Lee-Carter Model Parameters ===")
print("k_t (time index, first 5 values):", lc["k_t"].values[:5].round(4))
print("b_x (age sensitivity, top 5 ages):")
print(lc["b_x"].head(5).round(4).to_string())

# Forecast 20 years
mx_forecast = forecast_lee_carter(lc, n_years=20)
print(f"\nForecasted mx for age 65 (2021-2040):")
print(mx_forecast.loc[65].round(6))

# Plot k_t trend and forecast
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(lc["k_t"].index, lc["k_t"].values, "b-", linewidth=2, label="k_t (observed)")
forecast_k = [float(lc["k_t"].iloc[-1]) + lc["k_t"].diff().dropna().mean() * (i + 1)
              for i in range(20)]
forecast_yrs = list(range(years_lc[-1] + 1, years_lc[-1] + 21))
axes[0].plot(forecast_yrs, forecast_k, "r--", linewidth=2, label="k_t (forecast)")
axes[0].axvline(years_lc[-1], color="gray", linestyle=":", linewidth=1)
axes[0].set_xlabel("Year")
axes[0].set_ylabel("k_t")
axes[0].set_title("Lee-Carter Time Index k_t")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(lc["b_x"].index, lc["b_x"].values, "g-o", linewidth=2, markersize=5)
axes[1].set_xlabel("Age Group")
axes[1].set_ylabel("b_x (sensitivity)")
axes[1].set_title("Lee-Carter Age Sensitivity b_x")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("lee_carter_model.png", dpi=150)
plt.show()