actuary.py ( http://blog.xkcd.com/2012/07/12/a-morbid-python-script/ )
#!/usr/bin/python | |
import sys | |
import datetime | |
# Calculates death probabilities based on Social Security | |
# actuarial tables for a given group of people. | |
# Run with a list of ages/genders and an optional timespan (or year in the future): | |
# python actuary.py 63m 80m 75f 73m 10 | |
# or: | |
# python actuary.py 63m 80m 75f 73m 2022 | |
# This will give statistics for that group, including | |
# various probabilities over 10 years. Years can be | |
# ommitted and it will still give some statistics. | |
# If "Years" exceeds the current calendar year, it will be interpreted as a date. | |
# replaced with Japanese one | |
# http://www.mhlw.go.jp/toukei/saikin/hw/life/21th/index.html | |
bothtables=[[0.00037, 0.00026, 0.00018, 0.00013, 0.00011, 0.00010, 0.00009, 0.00008, 0.00008, 0.00008, 0.00010, 0.00011, 0.00013, 0.00015, 0.00019, 0.00024, 0.00030, 0.00038, 0.00045, 0.00051, 0.00057, 0.00061, 0.00064, 0.00064, 0.00064, 0.00065, 0.00066, 0.00067, 0.00068, 0.00069, 0.00071, 0.00074, 0.00077, 0.00081, 0.00085, 0.00090, 0.00098, 0.00108, 0.00118, 0.00128, 0.00140, 0.00152, 0.00166, 0.00181, 0.00198, 0.00216, 0.00238, 0.00263, 0.00289, 0.00317, 0.00347, 0.00381, 0.00419, 0.00461, 0.00507, 0.00558, 0.00612, 0.00669, 0.00732, 0.00810, 0.00888, 0.00961, 0.01037, 0.01121, 0.01214, 0.01319, 0.01434, 0.01553, 0.01685, 0.01842, 0.02023, 0.02227, 0.02466, 0.02753, 0.03087, 0.03478, 0.03919, 0.04420, 0.04974, 0.05568, 0.06208, 0.06937, 0.07793, 0.08752, 0.09785, 0.10827, 0.11926, 0.13135, 0.14503, 0.16041, 0.17569, 0.19195, 0.20922, 0.22755, 0.24695, 0.26744, 0.28905, 0.31177, 0.33560, 0.36051, 0.38649, 0.41348, 0.44142, 0.47023, 0.49980, 0.53002, 0.56075, 0.59182, 0.62304, 0.65422], [0.00033, 0.00023, 0.00015, 0.00011, 0.00009, 0.00008, 0.00008, 0.00007, 0.00006, 0.00006, 0.00006, 0.00007, 0.00008, 0.00010, 0.00012, 0.00014, 0.00016, 0.00019, 0.00021, 0.00024, 0.00025, 0.00026, 0.00026, 0.00026, 0.00026, 0.00027, 0.00028, 0.00031, 0.00034, 0.00036, 0.00038, 0.00040, 0.00042, 0.00045, 0.00048, 0.00052, 0.00056, 0.00061, 0.00066, 0.00071, 0.00076, 0.00082, 0.00091, 0.00100, 0.00108, 0.00115, 0.00124, 0.00138, 0.00153, 0.00167, 0.00179, 0.00191, 0.00204, 0.00219, 0.00236, 0.00254, 0.00273, 0.00292, 0.00313, 0.00340, 0.00370, 0.00401, 0.00434, 0.00465, 0.00498, 0.00537, 0.00584, 0.00634, 0.00692, 0.00767, 0.00856, 0.00961, 0.01083, 0.01224, 0.01381, 0.01558, 0.01763, 0.02006, 0.02284, 0.02600, 0.02956, 0.03371, 0.03867, 0.04458, 0.05155, 0.05937, 0.06837, 0.07843, 0.08949, 0.10160, 0.11490, 0.12964, 0.14628, 0.16458, 0.18367, 0.20330, 0.22398, 0.24573, 0.26854, 0.29242, 0.31734, 0.34328, 0.37021, 0.39806, 0.42678, 0.45627, 0.48643, 0.51715, 0.54827, 0.57965, 0.61112, 0.64247, 0.67352, 0.70404]] | |
def deathprob(age, years): | |
#negative ages = female | |
act=[] | |
if age<0: | |
act=bothtables[1] | |
age=-1*age | |
else: | |
act=bothtables[0] | |
while(len(act)<int(age+years+2)): # slower/bloaiter but keeps things clean | |
act.append(act[-1]**0.5) | |
liveprob=1 | |
i=0 | |
iage=int(age) | |
fage=age%1 | |
while i<=years-1: | |
thisyear=(1-fage)*act[iage+i]+fage*act[iage+i+1] | |
liveprob*=1-thisyear | |
i+=1 | |
if years%1: # Amortizes risk of dying over a partial year, which is | |
# 1-P(living last full year)^(year fraction) | |
lastyear=(1-fage)*act[iage+i]+fage*act[iage+i+1] | |
lastyearlive=1-lastyear | |
lastyearlive=lastyearlive**((years%1)) | |
liveprob*=lastyearlive | |
return 1-liveprob | |
def proballdie(ages, years): | |
probsliving=[] | |
for i in ages: | |
probsliving.append(1-deathprob(i, years)) | |
prod=1 | |
for i in probsliving: | |
prod*=(1-i) | |
return prod | |
def probanydie(ages, years): | |
probsliving=[] | |
for i in ages: | |
probsliving.append(1-deathprob(i, years)) | |
prod=1 | |
for i in probsliving: | |
prod*=i | |
return 1-prod | |
def calcexp(ages, prob, flag): | |
i=0 | |
for interval in (10, 1, 0.1, 0.01): | |
probs=0 | |
while(probs<prob): | |
i+=interval | |
if flag==0: | |
probs=proballdie(ages, i) | |
else: | |
probs=probanydie(ages, i) | |
i-=interval | |
return i | |
ages=[] | |
# print sys.argv[1:] | |
for arg in sys.argv[1:]: | |
gender=1 | |
years=1.0 | |
if arg[-1]=='m' or arg[-1]=='M': | |
try: | |
ages.append(1*float(arg[:-1])) | |
except: | |
print "Error parsing argument", arg | |
elif arg[-1]=='f' or arg[-1]=='F': | |
try: | |
ages.append(-1*float(arg[:-1])) | |
except: | |
print "Error parsing argument", arg | |
else: | |
try: | |
years=float(arg) | |
break | |
except: | |
print "Error parsing argument", arg | |
if not sys.argv[1:]: | |
print "The format is 'actuary.py 15m 80f 23', with a list of ages and a number of years to run the projections." | |
raise SystemExit | |
if not ages: | |
print "No ages specified. Format is 12m, 17f, etc." | |
raise SystemExit | |
# print "Ages:", ages | |
# print "Years:", years | |
(datetime.date.today()+datetime.timedelta(days=365.242191*1)).year | |
someone_years=[calcexp(ages, 0.05, 1), | |
calcexp(ages, 0.5, 1), | |
calcexp(ages, 0.95, 1)] | |
someone_dates=[(datetime.date.today()+datetime.timedelta(days=365.242191*someone_years[0])).year, | |
(datetime.date.today()+datetime.timedelta(days=365.242191*someone_years[1])).year, | |
(datetime.date.today()+datetime.timedelta(days=365.242191*someone_years[2])).year] | |
print "There is a 5% chance of someone dying within", someone_years[0], "years (by", str(someone_dates[0])+")." | |
print "There is a 50% chance of someone dying within", someone_years[1], "years (by", str(someone_dates[1])+")." | |
print "There is a 95% chance of someone dying within", someone_years[2], "years (by", str(someone_dates[2])+")." | |
print "" | |
if len(ages)>1: | |
everyone_years=[calcexp(ages, 0.05, 0), | |
calcexp(ages, 0.5, 0), | |
calcexp(ages, 0.95, 0)] | |
everyone_dates=[(datetime.date.today()+datetime.timedelta(days=365.242191*everyone_years[0])).year, | |
(datetime.date.today()+datetime.timedelta(days=365.242191*everyone_years[1])).year, | |
(datetime.date.today()+datetime.timedelta(days=365.242191*everyone_years[2])).year] | |
print "There is a 5% chance of everyone dying within", everyone_years[0], "years (by", str(everyone_dates[0])+")." | |
print "There is a 50% chance of everyone dying within", everyone_years[1], "years (by", str(everyone_dates[1])+")." | |
print "There is a 95% chance of everyone dying within", everyone_years[2], "years (by", str(everyone_dates[2])+")." | |
if years: | |
yearword="years" | |
if years==1: | |
yearword="year" | |
print "" | |
if years>datetime.date.today().year: | |
years=years-datetime.date.today().year | |
if len(ages)>1: | |
p=100*proballdie(ages, years) | |
printable="" | |
if p<0.001: | |
printable="<0.001" | |
elif p>99.99: | |
printable=">99.99" | |
else: | |
printable=str(p)[:5] | |
print "Probability of all dying in", years, yearword+": ", printable+"%" | |
p=100*probanydie(ages, years) | |
printable="" | |
if p<0.001: | |
printable="<0.001" | |
elif p>99.99: | |
printable=">99.99" | |
print p | |
else: | |
printable=str(p)[:5] | |
print "Probability of a death within", years, yearword+":", printable+"%" | |
raise SystemExit |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment