Widowhood simulation

import csv
import math
from pathlib import Path

SEED = 20260822
SIMULATIONS = 200_000
WEIBULL_SHAPE = 4
WOMAN_MEAN_REMAINING_YEARS = 25
MAN_MEAN_AT_SAME_AGE = 22
MALE_YEARS_LOST_PER_POSITIVE_AGE_GAP = 0.55
AGE_GAPS = range(-5, 11)

def create_rng(initial_seed):
    state = initial_seed & 0xFFFFFFFF

    def rng():
        nonlocal state
        state = (1664525 * state + 1013904223) & 0xFFFFFFFF
        return (state + 1) / 4294967297

    return rng

WEIBULL_MEAN_FACTOR = math.gamma(1 + 1 / WEIBULL_SHAPE)

def draw_weibull(mean, rng):
    scale = mean / WEIBULL_MEAN_FACTOR
    return scale * (-math.log(rng())) ** (1 / WEIBULL_SHAPE)

def rounded(value, digits=4):
    return round(value, digits)

def simulate_age_gap(age_gap, rng):
    male_mean = MAN_MEAN_AT_SAME_AGE - MALE_YEARS_LOST_PER_POSITIVE_AGE_GAP * age_gap
    woman_outlives = 0
    survivor_years = 0
    conditional_survivor_years = 0
    conditional_count = 0

    for _ in range(SIMULATIONS):
        woman_lifetime = draw_weibull(WOMAN_MEAN_REMAINING_YEARS, rng)
        man_lifetime = draw_weibull(male_mean, rng)
        difference = woman_lifetime - man_lifetime
        if difference > 0:
            woman_outlives += 1
            conditional_survivor_years += difference
            conditional_count += 1
        survivor_years += max(difference, 0)

    return {
        'age_gap': age_gap,
        'male_mean_remaining_years': rounded(male_mean, 3),
        'probability_woman_outlives': rounded(woman_outlives / SIMULATIONS),
        'expected_survivor_years': rounded(survivor_years / SIMULATIONS),
        'conditional_survivor_years': rounded(conditional_survivor_years / conditional_count),
        'simulations': SIMULATIONS,
    }

results = [simulate_age_gap(age_gap, create_rng(SEED + age_gap + 10)) for age_gap in AGE_GAPS]
results[5], results[-1]

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

output_directory = Path("public")
(output_directory / "analysis").mkdir(parents=True, exist_ok=True)
(output_directory / "figures").mkdir(parents=True, exist_ok=True)

gaps = [row["age_gap"] for row in results]
probabilities = [row["probability_woman_outlives"] for row in results]
survivor_years = [row["expected_survivor_years"] for row in results]

plt.rcParams.update({"font.family": "DejaVu Sans Mono", "font.size": 9})
fig, axes = plt.subplots(2, 1, figsize=(9.5, 6.25), sharex=True, constrained_layout=True)
fig.patch.set_facecolor("#f7f5f0")

axes[0].plot(gaps, probabilities, color="#38538d", marker="o", markersize=3, linewidth=2)
axes[0].set_title("Woman outlives man", loc="left", color="#171717", fontweight="bold")
axes[0].set_ylabel("probability", color="#6d6a63")
axes[0].set_ylim(0.35, 0.9)
axes[0].grid(axis="y", color="#c9c5bc", linestyle=(0, (3, 5)), linewidth=0.8)

axes[1].plot(gaps, survivor_years, color="#9b8060", marker="o", markersize=3, linewidth=2)
axes[1].set_title("Expected survivor years", loc="left", color="#171717", fontweight="bold")
axes[1].set_ylabel("years", color="#6d6a63")
axes[1].set_xlabel("husband age relative to a 60-year-old woman (years)", color="#6d6a63")
axes[1].set_ylim(0, 12)
axes[1].grid(axis="y", color="#c9c5bc", linestyle=(0, (3, 5)), linewidth=0.8)
axes[1].set_xticks([-5, 0, 5, 10])

for axis in axes:
    axis.set_facecolor("#f7f5f0")
    axis.spines["top"].set_visible(False)
    axis.spines["right"].set_visible(False)
    axis.spines["left"].set_color("#171717")
    axis.spines["bottom"].set_color("#171717")
    axis.tick_params(colors="#6d6a63")
    axis.grid(axis="x", color="#c9c5bc", linestyle=(0, (1, 4)), linewidth=0.7)

fig.savefig(output_directory / "figures" / "widowhood-simulation.svg", format="svg", metadata={"Date": None})
svg_path = output_directory / "figures" / "widowhood-simulation.svg"
svg_path.write_text("\n".join(line.rstrip() for line in svg_path.read_text().splitlines()) + "\n")
fig.savefig(output_directory / "figures" / "widowhood-simulation.png", format="png", dpi=180, facecolor=fig.get_facecolor())
plt.close(fig)

with (output_directory / "analysis" / "widowhood-simulation.csv").open("w", newline="") as handle:
    writer = csv.DictWriter(handle, fieldnames=results[0].keys())
    writer.writeheader()
    writer.writerows(results)

print(f"Generated {len(results)} age-gap scenarios with {SIMULATIONS} simulations each.")
print(f"Baseline age gap 0: P(woman outlives)={results[5]["probability_woman_outlives"]}, E[survivor years]={results[5]["expected_survivor_years"]}")
print(f"Age gap +10: P(woman outlives)={results[-1]["probability_woman_outlives"]}, E[survivor years]={results[-1]["expected_survivor_years"]}")