{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Widowhood simulation\n", "\n", "This notebook generates the article's illustrative Monte Carlo experiment. It is a pedagogical model, not an estimate fitted to a country's mortality table. The partners are simulated independently, with Weibull remaining-lifetime distributions, so shared behaviours, assortative health, and bereavement effects are intentionally excluded.\n", "\n", "The fixed seed, sample size, parameters, CSV output, and SVG figure make the result reproducible." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import csv\n", "import math\n", "from pathlib import Path\n", "\n", "SEED = 20260822\n", "SIMULATIONS = 200_000\n", "WEIBULL_SHAPE = 4\n", "WOMAN_MEAN_REMAINING_YEARS = 25\n", "MAN_MEAN_AT_SAME_AGE = 22\n", "MALE_YEARS_LOST_PER_POSITIVE_AGE_GAP = 0.55\n", "AGE_GAPS = range(-5, 11)\n", "\n", "def create_rng(initial_seed):\n", " state = initial_seed & 0xFFFFFFFF\n", "\n", " def rng():\n", " nonlocal state\n", " state = (1664525 * state + 1013904223) & 0xFFFFFFFF\n", " return (state + 1) / 4294967297\n", "\n", " return rng\n", "\n", "WEIBULL_MEAN_FACTOR = math.gamma(1 + 1 / WEIBULL_SHAPE)\n", "\n", "def draw_weibull(mean, rng):\n", " scale = mean / WEIBULL_MEAN_FACTOR\n", " return scale * (-math.log(rng())) ** (1 / WEIBULL_SHAPE)\n", "\n", "def rounded(value, digits=4):\n", " return round(value, digits)\n", "\n", "def simulate_age_gap(age_gap, rng):\n", " male_mean = MAN_MEAN_AT_SAME_AGE - MALE_YEARS_LOST_PER_POSITIVE_AGE_GAP * age_gap\n", " woman_outlives = 0\n", " survivor_years = 0\n", " conditional_survivor_years = 0\n", " conditional_count = 0\n", "\n", " for _ in range(SIMULATIONS):\n", " woman_lifetime = draw_weibull(WOMAN_MEAN_REMAINING_YEARS, rng)\n", " man_lifetime = draw_weibull(male_mean, rng)\n", " difference = woman_lifetime - man_lifetime\n", " if difference > 0:\n", " woman_outlives += 1\n", " conditional_survivor_years += difference\n", " conditional_count += 1\n", " survivor_years += max(difference, 0)\n", "\n", " return {\n", " 'age_gap': age_gap,\n", " 'male_mean_remaining_years': rounded(male_mean, 3),\n", " 'probability_woman_outlives': rounded(woman_outlives / SIMULATIONS),\n", " 'expected_survivor_years': rounded(survivor_years / SIMULATIONS),\n", " 'conditional_survivor_years': rounded(conditional_survivor_years / conditional_count),\n", " 'simulations': SIMULATIONS,\n", " }\n", "\n", "results = [simulate_age_gap(age_gap, create_rng(SEED + age_gap + 10)) for age_gap in AGE_GAPS]\n", "results[5], results[-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def line_path(values, x, y):\n", " commands = []\n", " for index, value in enumerate(values):\n", " command = 'M' if index == 0 else 'L'\n", " commands.append(f'{command}{x(index):.2f},{y(value):.2f}')\n", " return ' '.join(commands)\n", "\n", "def make_svg(rows):\n", " width, height = 760, 500\n", " left, right = 72, 28\n", " plot_width = width - left - right\n", " x = lambda index: left + plot_width * index / (len(rows) - 1)\n", " probability_y = lambda value: 210 - ((value - 0.35) / 0.35) * 150\n", " survivor_y = lambda value: 445 - (value / 12) * 150\n", " grid = ''.join(\n", " f''\n", " f'{value:.2f}'\n", " for value in (0.35, 0.45, 0.55, 0.65)\n", " )\n", " survivor_grid = ''.join(\n", " f'{value}'\n", " for value in (0, 4, 8, 12)\n", " )\n", " x_labels = ''.join(\n", " f'{row[\"age_gap\"]:+d}'\n", " for index, row in enumerate(rows)\n", " )\n", " probability_points = line_path([row['probability_woman_outlives'] for row in rows], x, probability_y)\n", " survivor_points = line_path([row['expected_survivor_years'] for row in rows], x, survivor_y)\n", " dots = ''.join(\n", " f''\n", " f''\n", " for index, row in enumerate(rows)\n", " )\n", " return f'''\n", " Illustrative simulation of age-gap effects on outsurvival probability and survivor years\n", " Two lines show an illustrative independent Weibull simulation for age gaps from five years younger to ten years older. The probability that the woman outlives the man and expected survivor years both increase as the man's age gap increases.\n", " \n", " Woman outlives man\n", " probability\n", " {grid}\n", " \n", " \n", " \n", " Expected survivor years\n", " years\n", " {survivor_grid}\n", " \n", " \n", " \n", " {dots}\n", " {x_labels}\n", " husband's age relative to a 60-year-old woman (years)\n", " outsurvival probability\n", " expected survivor years\n", "'''\n", "\n", "output_directory = Path('public')\n", "(output_directory / 'analysis').mkdir(parents=True, exist_ok=True)\n", "(output_directory / 'figures').mkdir(parents=True, exist_ok=True)\n", "\n", "with (output_directory / 'analysis' / 'widowhood-simulation.csv').open('w', newline='') as handle:\n", " writer = csv.DictWriter(handle, fieldnames=results[0].keys())\n", " writer.writeheader()\n", " writer.writerows(results)\n", "\n", "(output_directory / 'figures' / 'widowhood-simulation.svg').write_text(make_svg(results))\n", "print(f'Generated {len(results)} age-gap scenarios with {SIMULATIONS} simulations each.')\n", "print(f'Baseline age gap 0: P(woman outlives)={results[5][\"probability_woman_outlives\"]}, E[survivor years]={results[5][\"expected_survivor_years\"]}')\n", "print(f'Age gap +10: P(woman outlives)={results[-1][\"probability_woman_outlives\"]}, E[survivor years]={results[-1][\"expected_survivor_years\"]}')" ] } ], "metadata": { "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"name": "python", "version": "3.13"} }, "nbformat": 4, "nbformat_minor": 5 }