WER of VPW


In [150]:
%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import itertools
import functools
import operator

sns.set(color_codes=True)

Monte Carlo for Mortality


In [151]:
# based on http://www.ssa.gov/OACT/STATS/table4c6.html Period Life Table, 2009
# https://drive.google.com/open?id=0Bx-R9jPONgMgX0ktTUVFbVU0TEU
lifetable_for_women = [
	0.010298, # 65
	0.011281,
	0.012370,
	0.013572,
	0.014908,
	0.016440, # 70
	0.018162,
	0.020019,
	0.022003,
	0.024173,
	0.026706, # 75
	0.029603,
	0.032718,
	0.036034,
	0.039683,
	0.043899, # 80
	0.048807,
	0.054374,
	0.060661,
	0.067751,
	0.075729, # 85
	0.084673,
	0.094645,
	0.105694,
	0.117853,
	0.131146, # 90
	0.145585,
	0.161175,
	0.177910,
	0.195774,
	0.213849, # 95
    0.231865,
    0.249525,
    0.266514,
    0.282504,
    0.299455, # 100
    0.317422, 
    0.336467,
    0.356655, 
    0.378055, 
    0.400738, # 105
    0.424782, 
    0.450269, 
    0.477285, 
    0.505922, 
    0.536278, # 110
    0.568454, 
    0.602561, 
    0.638715, 
    0.677038, 
    0.717660, # 115
    0.760720, 
    0.806363, 
    0.851378,
    0.893947 # 119 
]

In [152]:
import random
def die(age):
    if age >= len(lifetable_for_women): return age
    if random.random() <= lifetable_for_women[age]:
        return age
    else:
        return die(age + 1)

In [153]:
def gen_age():
    return die(0) + 65

In [204]:
# Just to show what the overall mortality curve looks like
ages = []
for _ in range(250000):
    ages.append(gen_age())
ages = pd.Series(ages, name='Age at death')
sns.distplot(ages)


Out[204]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f11caede590>

Monte Carlo for Market Returns


In [155]:
# import US based data from 1871 from the simba backtesting spreadsheet found on bogleheads.org
# https://www.bogleheads.org/forum/viewtopic.php?p=2753812#p2753812
dataframe = pd.read_csv('returns.csv')
dataframe.head()


Out[155]:
Year Bonds Stocks Inflation
0 1871 2.46 15.12 1.53
1 1872 3.46 11.13 2.26
2 1873 8.00 -2.51 -4.41
3 1874 9.58 4.33 -6.92
4 1875 5.53 4.67 -5.79

In [156]:
def shuffle(df, n=1, axis=0):
    index = list(df.index)
    random.shuffle(index)
    df = df.ix[index]
    df.reset_index()
    return df

Generate monte carlo returns based on shuffling historical data.

The actual WER calculated is extremely sensitive to the model used for simulating returns. The more generous the model, the higher the SSR will be (since it can calculate a maximal strategy) and the lower the WER of most strategies.

This means that it is hard to compare WER's across research unless are using the same monte carlo simulations, as far as I can tell.

The relative performance of strategies should be the same across research but the absolute WER will vary.


In [206]:
def personal_returns(dataframe, age=None, inflation=False, adjustment=0.0):
    df = shuffle(dataframe.copy())
    if age == None:
        age = gen_age()
    lifetime_returns = []
    for _ in range(age - 64):
        row = df.iloc[_]
        bonds = row['Bonds']
        stocks = row['Stocks']
        annual_returns = (bonds * .4) + (stocks * .6)
        # The Morningstar paper isn't clear whether they use real returns or not.
        # I'll assume not and leave out the inflation adjustment.
        if inflation:
            annual_returns -= row['Inflation']
        if adjustment:
            annual_returns += adjustment
        lifetime_returns.append(annual_returns / 100)
    return lifetime_returns

In [207]:
def prod(x):
    return functools.reduce(operator.mul, x, 1)

In [208]:
def ssr(r):
    def ssr_seq(r):
        if len(r) == 0: return 1
        else:
            denom = prod(map(lambda x: x + 1, r))
            return (1 / denom) + ssr_seq(r[1:])

    # skip the final year of returns (since you have withdrawn the entire portfolio at the beginning of the year,
    # knowing that you will die this year)
    r_r = list(reversed(r[:-1]))
    return 1 / ssr_seq(r_r)

In [209]:
# test that I've implemented SSR correctly.

pr = [0.067040000000000002, -0.049419999999999999, 0.20482000000000003, 0.28746000000000005]
print('SSR  %f' % ssr(pr))

a0 = pr[0] + 1
a1 = pr[1] + 1
a2 = pr[2] + 1
a3 = pr[3] + 1

hand_ssr = (1.0
            / (1
                + (1/ (a0))
                + (1/ (a0 * a1))
                + (1/ (a0 * a1 * a2))
            )
)
print('Hand %f' % hand_ssr)


SSR  0.267283
Hand 0.267283

In [210]:
def cew(cashflows):
    # the paper says that results aren't sensitive to changing this (which turns out to be true!)
    gamma = 4.0
    
    def sigma(c):
        return pow(c, -gamma) / gamma
    
    constant_factor = 1.0 / len(cashflows) * gamma
    base = constant_factor * sum(map(sigma, cashflows))
    return pow(base, -1/gamma)

In [211]:
# test CEW
print(cew([0.05, 0.055, 0.06]))
print('0.05424659419') # calculated from a google spreadsheet


0.0542465941947
0.05424659419

In [212]:
# starting at age 65, based on 60% domestic stocks, 40% domestic bonds
vpw_rates = map(lambda x: x/100.0, [5.0,
    5.1,
    5.1,
    5.2,
    5.3,
    5.4,
    5.5,
    5.6,
    5.7,
    5.9,
    6.0,
    6.2,
    6.3,
    6.5,
    6.7,
    6.9,
    7.2,
    7.5,
    7.8,
    8.1,
    8.5,
    9.0,
    9.5,
    10.1,
    10.9,
    11.7,
    12.8,
    14.2,
    15.9,
    18.2,
    21.5,
    # diverge from the actual VPW...never exhaust the porfolio, just take 25% forever
    25.0, #25.4,
    25.0, #34.6,
    25.0, #50.9,
    25.0, #100
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
    25.0,
])
cew(vpw_rates)


Out[212]:
0.07515574150628429

In [213]:
def vpw_wer(p, s):   
    lifespan = len(p)
    v = vpw_rates[0:lifespan]
    
    w = cew(v) / s
    return w

In [214]:
def endow_wer(p, s):
    ''' use a constant 5.5% withdrawal '''
    lifespan = len(p)
    r = list(itertools.repeat(.055, lifespan))
    w = cew(r) / s
    return w

In [215]:
life_expectancy = [
# from the Annuity 2000 table A-1
# starting age 65
    23.02,
    22.16,
    21.31,
    20.47,
    19.63,
    18.81, # 70
    17.99,
    17.17,
    16.40,
    15.62,
    14.86, # 75
    14.12,
    13.40,
    12.69,
    12.01,
    11.34, #80
    10.70,
    10.08,
    9.48,
    8.91,
    8.37, # 85
    7.85,
    7.37,
    6.91,
    6.48,
    6.08, # 90
    5.72,
    5.38,
    5.06,
    4.77,
    4.50, # 95
    4.24,
    4.0,
    3.76,
    3.52,
    3.29, # 100
    3.05,
    2.81,
    2.58,
    2.35,
    2.12, # 105
    1.91,
    1.70,
    1.51,
    1.33,
    1.16, # 110
    # now just pad it out
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
]
rmd_rates = map(lambda x: 1.0 / x, life_expectancy)

In [216]:
def rmd_wer(p, s):
    lifespan = len(p)
    v = rmd_rates[0:lifespan]
    
    w = cew(v) / s
    return w

Compare Three Strategies


In [223]:
vpw_wers = []
endow_wers = []
rmd_wers = []

for _ in range(25000):
    p = personal_returns(dataframe.copy())
    
    s = ssr(p)
    
    vpw_wers.append(vpw_wer(p, s))
    endow_wers.append(endow_wer(p, s))
    rmd_wers.append(rmd_wer(p, s))
    
series = pd.Series({
    'VPW' : np.average(vpw_wers),
    'Endowment' : np.average(endow_wers),
    'RMD': np.average(rmd_wers),
})

series.plot(kind='bar')
series.head()


Out[223]:
Endowment    0.594138
RMD          0.646805
VPW          0.650810
dtype: float64

Show impact of inflation-adjusting market monte carlo


In [222]:
vpw1_wers = []

for _ in range(25000):
    p = personal_returns(dataframe.copy())    
    s = ssr(p)
    vpw1_wers.append(vpw_wer(p, s))

vpw2_wers = []
for _ in range(25000):
    p = personal_returns(dataframe.copy(), inflation=True)
    s = ssr(p)
    vpw2_wers.append(vpw_wer(p, s))
    
series = pd.Series({
    'Base' : np.average(vpw1_wers),
    'Inflation-adjusted' : np.average(vpw2_wers),
})

series.plot(kind='bar')
series.head()


Out[222]:
Base                  0.653433
Inflation-adjusted    0.798863
dtype: float64

In [187]:
# chance of dying in N years of retirement
years_of_retirement = 10
1 - prod(map(lambda x: 1-x, lifetable_for_women[0:years_of_retirement]))


Out[187]:
0.1518308184017616

In [188]:
pr = personal_returns(dataframe, age=65 - 1 + years_of_retirement)
s = ssr(pr)
print('Years of retirement %s. SSR %s' % (len(pr), s))
PORTFOLIO_START = 1000000

ssr_port = PORTFOLIO_START
ssr_wd = ssr_port * s
ssr_vals = []

vpw_port = PORTFOLIO_START
vpw_vals = []

for _ in range(len(pr)):
    # make the withdrawal at the start of the year
    ssr_port -= ssr_wd
    # then add in all of that year's gains
    ssr_port += ssr_port * pr[_]
    ssr_vals.append(ssr_port)
    
    vpw_wd = (vpw_rates[_] * vpw_port)
    vpw_port -= vpw_wd
    vpw_port += vpw_port * pr[_]
    vpw_vals.append(vpw_port)

plt.plot(ssr_vals)
plt.plot(vpw_vals)


Years of retirement 10. SSR 0.139255610428
Out[188]:
[<matplotlib.lines.Line2D at 0x7f11cb44d090>]

In [199]:
sns.tsplot(vpw_rates[0:30], color='green')
sns.tsplot(rmd_rates[0:30], color='blue')


Out[199]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f11caf2ca10>

In [203]:
vpw_vs_rmd = []
ssr_vs_rmd = []

ssrate = ssr(personal_returns(dataframe.copy(), age=95))
print('Sustainable rate for 30 years: %f' % ssrate)

for _ in range(20):
    vpw_w = vpw_rates[_]
    rmd_w = rmd_rates[_]
    
    delta = vpw_w - rmd_w
    vpw_vs_rmd.append(delta)

    delta = ssrate - rmd_w
    ssr_vs_rmd.append(delta)

    
sns.tsplot(vpw_vs_rmd, color='blue')
sns.tsplot(ssr_vs_rmd, color='green')


Sustainable rate for 30 years: 0.088126
Out[203]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f11cb2a4e10>