Compute probability that one person survives another

age-difference-mortality.pl

#! /usr/local/bin/perl -w

use strict;
use warnings;

# Get data from <URL: https://www.insee.fr/fr/statistiques/3311422?sommaire=3311425 >

# Typical plotting command (assuming output of this program is in "morta-hh.dat"):
# (echo 'set terminal pngcairo size 1000,800' ; echo 'set output "morta-hh.png"' ; echo 'set xrange [0:100]' ; echo 'set yrange [0:100]' ; echo 'set size square' ; echo 'set cbrange [0:1]' ; echo 'set palette defined ( 0 "black", 0.2 "#00FFFF", 0.4 "red", 0.5 "blue", 0.6 "green", 0.8 "#FF00FF", 1 "white" )' ; echo 'set title "Probabilité de survie ♂/♂"' ; echo 'set xlabel "Âge personne 1"' ; echo 'set ylabel "Âge personne 2"' ; echo 'set key title "Proba"' ; echo 'plot "morta-hh.dat" using ($1+0.5):($2+0.5):3 with image, [x=0:100] x lt rgb "black" notitle, [x=5:100] x-5 lt rgb "black" dashtype 2 notitle, [x=0:95] x+5 lt rgb "black" dashtype 2 notitle') | gnuplot

# Survival rate for 100000 men at a given age:
my @surva_h = ( 100000.00, 99646.38, 99584.48, 99560.65, 99542.91, 99529.11, 99517.40, 99506.68, 99497.36, 99489.02, 99481.17, 99473.28, 99465.23, 99456.63, 99446.25, 99433.04, 99416.08, 99394.21, 99364.61, 99325.52, 99277.60, 99223.36, 99165.14, 99104.38, 99041.02, 98975.56, 98906.62, 98835.96, 98763.15, 98688.44, 98611.89, 98532.92, 98451.44, 98366.94, 98277.47, 98183.71, 98085.12, 97981.69, 97870.62, 97748.41, 97615.30, 97471.03, 97314.06, 97144.09, 96952.84, 96743.99, 96515.56, 96262.12, 95980.04, 95671.84, 95333.22, 94962.76, 94560.57, 94122.31, 93639.94, 93106.27, 92510.59, 91868.98, 91164.55, 90404.11, 89603.44, 88726.32, 87811.08, 86843.01, 85817.07, 84749.15, 83631.21, 82468.54, 81255.23, 79963.47, 78609.85, 77171.48, 75664.91, 74074.52, 72372.36, 70591.31, 68690.65, 66661.32, 64472.47, 62143.84, 59602.88, 56899.27, 53995.91, 50879.45, 47590.06, 44103.65, 40456.55, 36669.42, 32798.20, 28877.75, 25031.96, 21323.65, 17793.68, 14522.03, 11572.22, 8965.70, 6746.70, 4945.35, 3519.28, 2414.38, 1597.26, 1046.63, 670.51, 414.35, 246.26 );

# Survival rate for 100000 women at a given age:
my @surva_f = ( 100000.00, 99709.56, 99653.44, 99632.06, 99618.85, 99608.19, 99599.28, 99591.19, 99583.82, 99577.03, 99570.56, 99564.03, 99557.25, 99550.27, 99542.39, 99532.98, 99521.71, 99509.19, 99495.12, 99478.57, 99459.49, 99438.31, 99417.14, 99395.89, 99374.43, 99352.55, 99329.59, 99306.28, 99281.48, 99255.03, 99225.95, 99195.28, 99162.72, 99128.74, 99092.06, 99053.52, 99011.43, 98963.89, 98909.98, 98850.36, 98785.28, 98713.71, 98635.64, 98550.09, 98452.29, 98341.74, 98218.83, 98083.20, 97935.39, 97766.63, 97588.96, 97394.55, 97180.49, 96952.72, 96707.35, 96439.52, 96153.77, 95842.86, 95508.46, 95151.03, 94766.53, 94365.48, 93935.19, 93486.05, 93003.14, 92500.80, 91955.40, 91387.85, 90785.65, 90138.55, 89440.14, 88671.27, 87854.59, 86952.35, 85965.66, 84912.19, 83742.14, 82441.89, 81017.82, 79454.64, 77705.78, 75724.62, 73497.10, 71021.92, 68256.28, 65174.13, 61752.51, 57992.80, 53893.42, 49548.25, 44936.82, 40138.32, 35235.66, 30337.45, 25567.35, 21029.21, 16874.09, 13147.01, 9937.02, 7272.47, 5171.31, 3556.91, 2337.24, 1491.06, 916.02 );

# Survival rate for 100000 men of lowest income vigintile:
my @surva_q01 = ( 100000.00, 99521.81, 99438.28, 99406.01, 99382.10, 99363.49, 99347.79, 99333.30, 99320.71, 99309.42, 99298.80, 99288.18, 99277.34, 99265.72, 99251.67, 99233.81, 99210.89, 99181.46, 99141.74, 99089.21, 99024.44, 98949.56, 98865.00, 98774.52, 98677.67, 98573.96, 98462.82, 98343.63, 98215.71, 98078.28, 97930.48, 97771.33, 97599.75, 97414.53, 97214.27, 96997.44, 96762.28, 96506.80, 96228.81, 95925.90, 95595.44, 95234.56, 94840.15, 94408.87, 93937.08, 93420.92, 92856.33, 92239.16, 91565.21, 90830.34, 90030.56, 89162.16, 88221.82, 87206.83, 86115.10, 84945.26, 83696.80, 82370.18, 80966.94, 79489.80, 77942.71, 76330.89, 74660.24, 72936.70, 71166.05, 69353.78, 67504.93, 65624.01, 63714.93, 61780.91, 59824.18, 57845.58, 55844.60, 53819.32, 51766.47, 49681.43, 47558.25, 45389.73, 43168.57, 40889.12, 38547.85, 36143.66, 33678.45, 31157.55, 28590.33, 25990.60, 23377.15, 20774.04, 18210.31, 15719.21, 13336.69, 11071.65, 8993.53, 7127.94, 5493.83, 4101.89, 2964.74, 2065.44, 1379.92, 878.80, 529.60, 314.91, 182.04, 100.48, 52.68 );

# Survival rate for 100000 men of highest income vigintile:
my @surva_q20 = ( 100000.00, 99740.18, 99694.69, 99677.11, 99664.08, 99653.94, 99645.39, 99637.49, 99630.63, 99624.48, 99618.69, 99612.90, 99606.99, 99600.66, 99593.00, 99583.26, 99570.76, 99554.71, 99533.05, 99504.39, 99469.05, 99432.26, 99395.45, 99358.60, 99321.65, 99284.56, 99247.26, 99209.67, 99171.69, 99133.21, 99094.11, 99054.25, 99013.47, 98971.58, 98928.40, 98883.68, 98837.16, 98788.54, 98737.49, 98683.65, 98626.60, 98565.93, 98501.15, 98431.71, 98357.03, 98276.46, 98189.25, 98094.61, 97991.65, 97879.40, 97756.79, 97622.66, 97475.73, 97314.62, 97137.83, 96943.76, 96730.65, 96496.66, 96239.82, 95958.03, 95649.14, 95310.88, 94940.84, 94536.25, 94093.96, 93610.35, 93081.21, 92501.66, 91866.00, 91167.55, 90398.59, 89550.26, 88612.45, 87573.59, 86420.46, 85137.99, 83708.99, 82114.01, 80331.64, 78339.50, 76114.94, 73636.00, 70882.82, 67839.28, 64495.21, 60848.81, 56908.74, 52695.54, 48243.56, 43602.38, 38837.26, 34028.30, 29267.87, 24656.12, 20293.79, 16273.40, 12672.20, 9546.15, 6925.46, 4812.42, 3231.07, 2147.91, 1398.47, 880.35, 534.02 );

# Complete the data to avoid boundary effects:
my @alltabs = ( \@surva_h, \@surva_f, \@surva_q01, \@surva_q20 );
foreach my $sr ( @alltabs ) {
    my $len = scalar(@{$sr});
    # Assume the mortality rate for very high ages stays constant:
    my $rat = $sr->[$len-1] / $sr->[$len-2];
    my $v = $sr->[$len-1];
    while ( $v>0.1 ) {
	$v *= $rat;
	push @{$sr}, $v;
    }
    # Kill off the one-in-a-millionth survivor:
    push @{$sr}, 0;
}

my @survax = @surva_h;
my @survay = @surva_h;

for ( my $agex=0 ; $agex<100 ; $agex++ ) {
    for ( my $agey=0 ; $agey<100 ; $agey++ ) {
	my $sum = 0;
	for ( my $u=0 ; $agex+$u<scalar(@survax) && $agey+$u<scalar(@survay) ; $u++ ) {
	    # Basically $diff is the number dying in year $u (aged
	    # $agey+$u), but make it symmetric:
	    my $diff = (($agey+$u+1<scalar(@survay) ? $survay[$agey+$u]-$survay[$agey+$u+1] : 0) + ($u>0 ? $survay[$agey+$u-1]-$survay[$agey+$u] : 0))/2;
	    $sum += $survax[$agex+$u] * $diff;
	}
	my $prob = $sum / $survax[$agex] / $survay[$agey];
	printf "%d\t%d\t%.4f\n", $agex, $agey, $prob;
    }
}