Safe Retirement Spending Using Certainty Equivalent Cash Flow and Tensorflow

TensorFlow code to find a retirement spending strategy that would have optimized certainty equivalent cash flow over all 30-year historical retirement cohorts, 1928-1986.

To calculate CE cash flow, first calculate CRRA utility for a retiree's cash flows and risk aversion parameter γ (gamma):

$$U = \frac{1}{n}\sum_{i}\frac{C_i^{1-\gamma}-1}{1-\gamma}$$

Then convert average CRRA utility back to a cash flow equivalent using the inverse function.

$$CE = [U(1-\gamma) + 1] ^ {\frac{1}{1-\gamma}}$$
  • γ = 0 means you’re risk neutral. There is no discount, however variable or uncertain the cash flows. The CE value equals the average of the cash flows.
  • γ = 8 means you’re fairly risk averse. The CE value reflects a large discount vs. the average cash flow.
  • The higher the variability of the cash flows, the greater the discount. And the higher the γ parameter, the greater the discount.

We model retirement cash flows as a function of:

  • Constant spending (a single value): A constant inflation-adjusted amount you withdraw each year in retirement. This is like the 4% in Bengen’s 4% rule. The inflation-adjusted value of this annual withdrawal never changes.
  • Variable spending (30 values, one for each year, i.e. a list or vector): A variable percentage of your portfolio value you withdraw each year. In contrast to the Bengen 4% rule, we’re, saying, if the portfolio appreciates, you can safely withdraw an additional amount based on the current value of the portfolio. Your total spending is the sum of 1) constant spending and 2) variable spending * portfolio value.
  • Stock allocation (30 values, one for each year): We study a portfolio with 2 assets: S&P 500 stocks and 10-year Treasurys. Bond allocation = 1 - stock allocation.

Using these variables, we construct a TensorFlow graph.

Constants:

  • γ = 8.
  • A portfolio starting value: 100.
  • Inflation-adjusted stock returns 1928-2015 (all numbers we use are inflation-adjusted, and we maximize inflation-adjusted cash flow).
  • Inflation-adjusted bond returns 1928-2015.

Operations:

  • Calculate 59 30-vectors, each one representing the cash flow of one 30-year retirement cohort 1928-1986, using the given constant spending, variable spending, and stock allocation.
  • Calculate the certainty equivalent cash flow of each cohort using γ.
  • Calculate the certainty equivalent cash flow over all cohorts using γ.

Use TensorFlow GradientDescentOptimizer to find the variables that resulted in the highest CE spending over all cohorts.


In [1]:
from __future__ import print_function

import argparse
import pickle
from time import strftime
import sys
import six
import random

import tensorflow as tf
import numpy as np
import pandas as pd
import lifetable

import matplotlib.pyplot as plt
%matplotlib inline

import matplotlib.mlab as mlab
import plotly.plotly as py


# TensorFlow numeric type to use for floating point variables
# tf.float32 is 2x faster but less accurate
# tf.float64 will run out of accuracy for high gamma (> 8)
float_type = tf.float64

In [2]:
############################################################
# returns 1928-2015
############################################################

first_year = 1928
last_year = 2015
years = range(first_year, last_year+1) # pythonically yields [1928, 1929...2015]
years_history = len(years)
years_retired = 30
num_cohorts = years_history - years_retired + 1
num_assets = 2

bestfile = "jupyter"
#gamma = 1.0

sp500 = pd.Series([
    0.4381,-0.083,-0.2512,-0.4384,-0.0864,0.4998,-0.0119,0.4674,0.3194,-0.3534,0.2928,-0.011,
    -0.1067,-0.1277,0.1917,0.2506,0.1903,0.3582,-0.0843,0.052,0.057,0.183,0.3081,0.2368,0.1815,
    -0.0121,0.5256,0.326,0.0744,-0.1046,0.4372,0.1206,0.0034,0.2664,-0.0881,0.2261,0.1642,0.124,
    -0.0997,0.238,0.1081,-0.0824,0.0356,0.1422,0.1876,-0.1431,-0.259,0.37,0.2383,-0.0698,0.0651,
    0.1852,0.3174,-0.047,0.2042,0.2234,0.0615,0.3124,0.1849,0.0581,0.1654,0.3148,-0.0306,0.3023,
    0.0749,0.0997,0.0133,0.372,0.2268,0.331,0.2834,0.2089,-0.0903,-0.1185,-0.2197,0.2836,0.1074,
    0.0483,0.1561,0.0548,-0.3655,0.2594,0.1482,0.021,0.1589,0.3215,0.1352,0.0136],
    index = years)

bonds=pd.Series([
    0.0084,0.042,0.0454,-0.0256,0.0879,0.0186,0.0796,0.0447,0.0502,0.0138,0.0421,0.0441,
    0.054,-0.0202,0.0229,0.0249,0.0258,0.038,0.0313,0.0092,0.0195,0.0466,0.0043,-0.003,
    0.0227,0.0414,0.0329,-0.0134,-0.0226,0.068,-0.021,-0.0265,0.1164,0.0206,0.0569,0.0168,
    0.0373,0.0072,0.0291,-0.0158,0.0327,-0.0501,0.1675,0.0979,0.0282,0.0366,0.0199,0.0361,
    0.1598,0.0129,-0.0078,0.0067,-0.0299,0.082,0.3281,0.032,0.1373,0.2571,0.2428,-0.0496,
    0.0822,0.1769,0.0624,0.15,0.0936,0.1421,-0.0804,0.2348,0.0143,0.0994,0.1492,-0.0825,
    0.1666,0.0557,0.1512,0.0038,0.0449,0.0287,0.0196,0.1021,0.201,-0.1112,0.0846,0.1604,
    0.0297,-0.091,0.1075,0.0128],
                index=years)

cpi=pd.Series([
    -0.0115607,0.005848,-0.0639535,-0.0931677,-0.1027397,0.0076336,0.0151515,0.0298507,
    0.0144928,0.0285714,-0.0277778,0,0.0071429,0.0992908,0.0903226,0.0295858,0.0229885,
    0.0224719,0.1813187,0.0883721,0.0299145,-0.0207469,0.059322,0.06,0.0075472,0.0074906,
    -0.0074349,0.0037453,0.0298507,0.0289855,0.0176056,0.017301,0.0136054,0.0067114,0.0133333,
    0.0164474,0.0097087,0.0192308,0.0345912,0.0303951,0.0471976,0.0619718,0.0557029,0.0326633,
    0.0340633,0.0870588,0.1233766,0.0693642,0.0486486,0.0670103,0.0901771,0.1329394,0.125163,
    0.0892236,0.0382979,0.0379098,0.0394867,0.0379867,0.010979,0.0443439,0.0441941,0.046473,
    0.0610626,0.0306428,0.0290065,0.0274841,0.026749,0.0253841,0.0332248,0.017024,0.016119,
    0.0268456,0.0338681,0.0155172,0.0237691,0.0187949,0.0325556,0.0341566,0.0254065,0.0408127,
    0.0009141,0.0272133,0.0149572,0.0296,0.0174,0.015,0.0076,0.0073],
              index=years)

real_stocks = sp500 - cpi
real_bonds = bonds - cpi

startval = 100
years_retired = 30
# 1% constant spending
const_spend_pct = .0225
const_spend = startval * const_spend_pct

# var spending a function of years left
var_spend_pcts = pd.Series([ 0.5/(30-ix) for ix in range(30)])
var_spend_pcts[29] = 1.0

# stocks starting at 82%, decreasing 0.5% per year
stock_allocations = pd.Series([0.82 - 0.005* ix for ix in range(30)])
bond_allocations = 1 - stock_allocations

pickle_list = [const_spend, var_spend_pcts, stock_allocations, bond_allocations]
pickle.dump( pickle_list, open( bestfile, "wb" ) )

In [3]:
class SafeWithdrawalModel:
    """initialize graph and parameters shared by all retirement cohorts"""
    def __init__(self,
                 returns_list, # series with returns for assets
                 names_list, # names of assets
                 allocations_list, # list of % allocated to each asset class
                 start_val, # starting portfolio value e.g. 100
                 const_spend_pct,
                 var_spend_pcts,
                 gamma,
                 survival,
                 verbose=False):

        # read params, initialize Tensorflow graph and session
        # set up ops specific to model
        self.verbose=verbose
        self.startval=startval
        self.returns_list = returns_list
        self.names_list = names_list
        self.num_assets = len(self.names_list)
        self.start_val = start_val
        self.ret_years = len(allocations_list[0])
        self.const_spend_pct = const_spend_pct
        self.var_spend_pcts = var_spend_pcts
        self.survival=survival
        self.gamma = gamma

        # model will have a cohort_history object, optimizer object
        # initialize with placeholder, needs rest of model initialized first
        self.cohort_history = None
        self.optimizer = None

        self.first_year = returns_list[0].index[0]
        self.last_year = returns_list[0].index[-1]
        self.total_cohorts = len(returns_list[0])
        self.ret_cohorts = self.total_cohorts - self.ret_years + 1

        print('%s Create TensorFlow graph and session' % strftime("%H:%M:%S"))
        self.graph = tf.Graph()
        self.sess = tf.Session(graph = self.graph)
        self.return_ops = []
        self.allocation_ops = []

        with self.graph.as_default():
            with tf.device("/cpu:0"):
                # some constants
                self.zero = tf.constant(0.0, dtype=float_type, name="zero")
                self.one = tf.constant(1.0, dtype=float_type, name="one")
                self.one_hundred = tf.constant(100.0, dtype=float_type, name="one_hundred")
                self.ten_thousand = tf.constant(10000.0, dtype=float_type, name="ten_thousand")
                self.one_hundred_thousand = tf.constant(100000.0, dtype=float_type, name="one_million")
                self.one_million = tf.constant(1000000.0, dtype=float_type, name="one_million")
                self.very_small_amts = tf.constant(np.array([0.000001] * self.ret_years),
                                                   dtype=float_type, name="very_small_amts")
                self.zero_years = tf.constant(np.zeros(self.ret_years),
                                              dtype=float_type, name = "zero_years")
                self.one_years = tf.constant(np.ones(self.ret_years), dtype=float_type, name="one_years")
                self.ret_years_op = tf.constant(self.ret_years, dtype=float_type, name="ret_years")
                #gamma
                self.gamma_op = tf.constant(gamma, dtype=float_type, name="gamma")
                self.one_minus_gamma = tf.sub(self.one, self.gamma, name="one_minus_gamma")
                self.inv_one_minus_gamma = tf.div(self.one, self.one_minus_gamma,
                                                  name="inv_one_minus_gamma")

                self.cost_multiplier = self.ten_thousand

                # generate op for start_val
                self.start_val_op = tf.constant(100.0, dtype=float_type, name ="port_start_val")

                # generate ops for returns
                for prefix, return_series in zip(names_list, returns_list):
                    self.return_ops.append(self.gen_tf_const_list(return_series, "%s_return" % prefix,
                                                                  verbose=self.verbose))

                # only implemented for n=2 assets
                # generate ops for allocations for first n-1 assets
                prefix = names_list[0]
                alloc_series = allocations_list[0]
                stock_alloc_ops = self.gen_tf_var_list(alloc_series, "%s_alloc" % prefix,
                                                       verbose=self.verbose)
                self.allocation_ops.append(stock_alloc_ops)

                # ops for soft constraints: 0 < stock allocation < 1
                self.alloc_min_0_ops = self.gen_zero_min_list(stock_alloc_ops, "alloc_min_0",
                                                              verbose=self.verbose)
                self.cost_alloc_min_0_op = tf.mul(self.cost_multiplier,
                                                  tf.add_n(self.alloc_min_0_ops,
                                                           name="cost_alloc_min_0"),
                                                  name="cost_alloc_min_0_mult")

                self.alloc_max_1_ops = self.gen_one_max_list(stock_alloc_ops, "alloc_max_1",
                                                             verbose=self.verbose)
                self.cost_alloc_max_1_op = tf.mul(self.cost_multiplier,
                                                  tf.add_n(self.alloc_max_1_ops, name = "cost_alloc_max_1"))

                # ops for soft constraints: declining stock allocation
                # why? for example, 1966 is the worst cohort, and 1974 is its worst stock return (-40%)
                # to maximize CE, optimization sets stock allocation at a minimum to not run out of money
                # in worst year. It will go e.g. 80% stock alloc in year 8 and 56% in year 9, return to
                # 80% in year 10.To avoid artifacts like that, knowing stock allocation should decline
                # over time, we add a large penalty to objective when stock allocation increases
                # from one year to next.
                self.alloc_decrease_ops = self.gen_diff_list(stock_alloc_ops, "alloc_decrease",
                                                             verbose=self.verbose)
                self.cost_alloc_decrease_op = tf.mul(self.cost_multiplier,
                                                     tf.add_n(self.alloc_decrease_ops,
                                                              name="alloc_decrease_cost_op"))
                # last asset is 1-previous assets
                bond_alloc_ops = []

                var_prefix = "%s_alloc" % names_list[1]
                print ('%s Create ops for %s' % (strftime("%H:%M:%S"), var_prefix))
                for ix, op in enumerate(stock_alloc_ops):
                    var_name = "%s_%d" % (var_prefix, ix)
                    if self.verbose:
                        print('Create %s' % var_name)
                    var_op = tf.sub(self.one, stock_alloc_ops[ix], name=var_name)
                    bond_alloc_ops.append(var_op)
                self.allocation_ops.append(bond_alloc_ops)

                # generate ops for const, var spending
                self.const_spend_pct_op = tf.Variable(const_spend_pct, dtype=float_type, name="const_spend_pct")
                self.sess.run(self.const_spend_pct_op.initializer)
                self.const_spending_op = tf.mul(self.const_spend_pct_op, self.one_hundred, name="const_spend")

                self.var_spending_ops = self.gen_tf_var_list(self.var_spend_pcts, "var_spend",
                                                             verbose=self.verbose)

                # all ops to be trained
                self.all_var_ops = [self.const_spend_pct_op] + \
                                   self.var_spending_ops + \
                                   self.allocation_ops[0]

                # op for soft constraint: const spending > 0
                self.cspend_min_0_op = tf.maximum(self.zero, tf.neg(self.const_spend_pct_op,
                                                                    name="neg_cspend_min_0_op"),
                                                  name="cspend_min_0_op")
                self.cost_cspend_min_0_op = tf.mul(self.cost_multiplier,
                                                   self.cspend_min_0_op,
                                                   name="cost_cspend_min_0")

                # op for soft constraint: var spending > 0
                self.vspend_min_0_ops = self.gen_zero_min_list(self.var_spending_ops, "vspend_min_0",
                                                               verbose=self.verbose)
                self.cost_vspend_min_0_op = tf.mul(self.cost_multiplier,
                                                   tf.add_n(self.vspend_min_0_ops,
                                                            name="cost_vspend_min_0"))

                if survival is not None:
                    survival_array=np.array(survival)
                    self.survival_tensor = tf.constant(survival_array, dtype=float_type,
                                                       name="survival_tensor")

                # global step counter
                self.step_count = tf.Variable(0, dtype=float_type, name="step_count", trainable=False)
                self.increment_step = self.step_count.assign_add(1)

                #init op
                self.init_op = tf.initialize_all_variables()

    def __del__(self):
        """When deleting model, close session, clear default graph"""
        print("Destructor reset graph")
        try:
            with self.graph.as_default():
                tf.reset_default_graph()
        except Exception, e:
            print ("Destructor couldn't reset graph: %s" % str(e))

        try:
            print ("Destructor close Tensorflow session")
            self.sess.close()
        except Exception, e:
            print ("Destructor couldn't close session: %s" % str(e))


    def gen_tf_const_list(self, const_iter, const_prefix, start_index=0, verbose=False):
        """take a list or iterator of values, generate and return tensorflow constant ops for each"""
        print ('%s Create constants %s' % (strftime("%H:%M:%S"), const_prefix))
        with self.graph.as_default():
            with tf.device("/cpu:0"):

                const_list = []
                for ix, const in enumerate(const_iter):
                    const_name = "%s_%d" % (const_prefix, start_index + ix)
                    if verbose:
                        print("Set constant %s to %f" % (const_name, const))
                    const_list.append(tf.constant(const, dtype=float_type, name=const_name))

        return const_list

    def gen_tf_var_list(self, var_iter, var_prefix, start_index=0, verbose=False):
        """take a list or iterator of values, generate and return tensorflow Variable ops for each"""
        print ('%s Create variables %s' % (strftime("%H:%M:%S"), var_prefix))
        with self.graph.as_default():
            with tf.device("/cpu:0"):
                var_list = []
                for ix, var in enumerate(var_iter):
                    var_name = "%s_%d" % (var_prefix, start_index + ix)
                    if verbose:
                        print("Create variable %s to %f" % (var_name, var))
                    var_op = tf.Variable(var, dtype=float_type, name=var_name)
                    self.sess.run(var_op.initializer)
                    var_list.append(var_op)

        return var_list

    def get_op_from_list(self, op_list, op_index):
        """take a list of ops, return value of op specified by op_index"""
        op = op_list[op_index]
        retval = self.sess.run([op])
        return retval

    def gen_zero_min_list(self, op_iter, op_prefix, start_index=0, verbose=False):
        """take a list or iterator of ops, generate and return an op which is max(-op, 0)
        for soft constraints > 0"""
        print ('%s Create ops for soft constraint %s > 0' % (strftime("%H:%M:%S"), op_prefix))
        with self.graph.as_default():
            with tf.device("/cpu:0"):
                op_list = []
                for ix, op in enumerate(op_iter):
                    op_name = "%s_%d" % (op_prefix, start_index + ix)
                    if verbose:
                        print("Zero_min op %s" % op_name)
                    new_op = tf.maximum(self.zero, tf.neg(op, name="neg_%s" % op_name), name=op_name)
                    op_list.append(new_op)

        return op_list

    def gen_one_max_list(self, op_iter, op_prefix, start_index=0, verbose=False):
        """take a list or iterator of ops, generate and return an op with is max(op-1, 0)
        for soft constraints > 0"""
        print ('%s Create ops for soft constraint %s < 1' % (strftime("%H:%M:%S"), op_prefix))
        with self.graph.as_default():
            with tf.device("/cpu:0"):
                op_list = []
                for ix, op in enumerate(op_iter):
                    op_name = "%s_%d" % (op_prefix, start_index + ix)
                    if verbose:
                        print('One_max op %s' % op_name)
                    new_op = tf.maximum(self.zero, tf.sub(op, self.one, name="one_minus_%s" % op_name),
                        name=op_name)
                    op_list.append(new_op)

        return op_list

    def gen_diff_list(self, op_iter, op_prefix, start_index=0, verbose=False):
        """generate and return an op for declining stock alloc constraint over time, max of 0 and decrease"""
        print ('%s Create ops for soft constraint, declining stock alloc %s' % (strftime("%H:%M:%S"),
                                                                                op_prefix))
        with self.graph.as_default():
            with tf.device("/cpu:0"):
                op_list = []
                for ix, op in enumerate(op_iter):
                    if ix == 0:
                        continue
                    op_name = "%s_%d" % (op_prefix, start_index + ix)
                    if verbose:
                        print("diff op %s" % op_name)

                    new_op = tf.maximum(self.zero, tf.sub(op_iter[ix], op_iter[ix-1]))
                    op_list.append(new_op)

        return op_list

    def gen_ce(self, input_tensor, prefix, survival_tensor=None, verbose=False):
      with tf.device("/cpu:0"):
        with self.graph.as_default():
            input_length = np.float64(input_tensor.get_shape().as_list()[0])

            if verbose:
                print("%s Create ce op with gamma: %f" % (strftime("%H:%M:%S"), self.gamma))
            if self.gamma == 1.0:
                u = tf.reduce_mean(tf.log(input_tensor), name="%s_u" % prefix)
                #print(self.sess.run(u))
                if survival_tensor is not None:
                    u0 = u
                    u = tf.reduce_mean(tf.mul(u0, survival_tensor, name="%s_u_surv" % prefix),
                                       name="%s_u" % prefix)
                ce = tf.exp(u, name="%s_ce" % prefix)
                if verbose:
                    print ('%s Create CE op %f' % (strftime("%H:%M:%S"), self.sess.run(ce)))
            else:
                # for high gamma numerical error is significant, calculation is most accurate near 1
                # so divide by mean
                input_mean = tf.reduce_mean(input_tensor, name="%s_mean" % prefix)
                input_conditioned = tf.div(input_tensor, input_mean, name="%s_conditioned" % prefix)

                u1 = tf.pow(input_conditioned, self.one_minus_gamma, name="%s_u1" % prefix)
                u2 = tf.sub(u1, self.one, name="%s_u2" % prefix)
                u3 = tf.mul(u2, self.inv_one_minus_gamma, name="%s_u3" % prefix)
                u = tf.reduce_mean(u3, name="%s_u" % prefix)

                if survival_tensor is not None:
                    u4 = u
                    u = tf.reduce_mean(tf.mul(u4, survival_tensor, name="%s_u_surv" % prefix),
                        name="%s_u" % prefix)

                ce1 = tf.mul(self.one_minus_gamma, u, name="%s_ce1" % prefix)
                ce2 = tf.add(ce1, self.one, name="%s_ce2" % prefix)
                ce3 = tf.pow(ce2, self.inv_one_minus_gamma, name="%s_ce3" % prefix)
                ce = tf.mul(input_mean, ce3, name="%s_ce" % prefix)

                if verbose:
                    print ('%s Create CE op %f' % (strftime("%H:%M:%S"), self.sess.run(ce)))
            return ce

In [4]:
class Cohort:
    """Cohort represents experience of an individual
       - retiring in a given year
       - using the specified SafeWithdrawal model"""
    def __init__(self, model, cohort_start_year):
        self.model = model
        self.cohort_start_year = cohort_start_year
        self.name = "cohort_%d" % cohort_start_year
        self.gen_tf_ops()

    def gen_tf_ops(self, verbose=False):
        
        if verbose:
            print("%s Instantiating cohort %s" % (strftime("%H:%M:%S"), self.name))

        stock_returns = self.model.return_ops[0]
        bond_returns = self.model.return_ops[1]

        stock_allocs = self.model.allocation_ops[0]
        bond_allocs = self.model.allocation_ops[1]

        self.port_returns_list = []
        self.port_prespend_list = []
        self.port_end_vals_list = []
        self.spend_amts_list = []
        self.spend_amts_nonzero_list = []

        with self.model.graph.as_default():
          with tf.device("/cpu:0"):

            if verbose:
                print ("%s Generating %d years from %d" % (strftime("%H:%M:%S"),
                                                           self.model.ret_years,
                                                           self.cohort_start_year))

            start_year_ix = self.cohort_start_year - self.model.first_year
            for ix in range(self.model.ret_years):
                op_stock_return = stock_returns[start_year_ix + ix]
                op_stock_alloc = stock_allocs[ix]
                op_bond_return = bond_returns[start_year_ix + ix]
                op_bond_alloc = bond_allocs[ix]

                op_const_spend = self.model.const_spending_op
                op_var_spend = self.model.var_spending_ops[ix]

                op_total_real_return = tf.add(tf.mul(op_stock_alloc, op_stock_return, name="%s_stock_%d"
                                                     % (self.name, ix)),
                                              tf.mul(op_bond_alloc, op_bond_return, name="%s_bond_%d"
                                                     % (self.name, ix)),
                                              name="%s_total_return_%d" % (self.name, ix))
                self.port_returns_list.append(op_total_real_return)

                if ix == 0:
                    prev_val = self.model.start_val_op
                else:
                    prev_val = self.port_end_vals_list[ix-1]

                op_port_end_val_prespend = tf.add(prev_val,
                                                  tf.mul(prev_val, self.port_returns_list[ix],
                                                         name="%s_dolreturn_%d" % (self.name, ix)),
                                                  name="%s_prespend_%d" % (self.name, ix))
                self.port_prespend_list.append(op_port_end_val_prespend)

                desired_spend_amt = tf.add(tf.mul(op_var_spend, op_port_end_val_prespend,
                    name="%s_des_vspend_%d" % (self.name, ix)),
                                           op_const_spend,
                                           name="%s_desired_spend_amt_%d" % (self.name, ix))
                #spend minimum of tmp_spend_amt, port value
                spend_amt = tf.minimum(desired_spend_amt, op_port_end_val_prespend,
                                       name="%s_actual_spend_amt_%d" % (self.name, ix))
                self.spend_amts_list.append(spend_amt)

                op_port_end_val = tf.sub(op_port_end_val_prespend, spend_amt, name="%s_endval_%d" %
                                                                                   (self.name, ix))
                self.port_end_vals_list.append(op_port_end_val)

            #now that we've computed cohort paths we pack results into 1D Tensors to calc objective
            self.spend_amts = tf.pack(self.spend_amts_list, name="%s_spend_amts" % self.name)
            self.port_end_vals = tf.pack(self.port_end_vals_list, name="%s_port_end_vals" % self.name)

            self.mean_spending = tf.reduce_mean(self.spend_amts, name="%s_mean_spending" % self.name)
            self.sd_spending = tf.sqrt(tf.reduce_mean(tf.pow(tf.sub(self.spend_amts,
                                                                    self.mean_spending), 2)),
                                       name="%s_sd_spending" % self.name)
            self.min_spending = tf.reduce_min(self.spend_amts, name="%s_min_spending" % self.name)
            self.max_spending = tf.reduce_max(self.spend_amts, name="%s_max_spending" % self.name)

            if self.model.gamma == 1.0:
                #spend a tiny amount even if spend is 0 so log is not NaN
                #doesn't really seem like best practice but...
                #0 spend years can't be in final solution
                #and don't want divide by zero errors if optimizer attempts one

                #chain new op off old op but keep a reference to old op around just in case
                self.spend_amts_maybe_zero = self.spend_amts
                self.spend_amts = tf.maximum(self.spend_amts_maybe_zero,
                                             self.model.very_small_amts,
                                             name="%s_actual_spend_nonzero" % self.name)
                self.total_spending = tf.reduce_sum(self.spend_amts, name="%s_total_spending_nonzero" %
                                                    self.name)
            else:
                self.total_spending = tf.reduce_sum(self.spend_amts, name="%s_total_spending" %
                                                    self.name)

            if self.model.survival is not None:
                self.ce = self.model.gen_ce_survival(self.spend_amts,
                                                     self.model.survival_tensor,
                                                     "%s_ce" % self.name)
            else:
                self.ce = self.model.gen_ce(self.spend_amts,
                                            "%s_ce" % self.name)

            #print (self.as_dataframe())

    def get_tf_ops(self):
        return self.model.start_val, self.port_returns_list, self.port_prespend_list, \
            self.spend_amts_list, self.port_end_vals_list, self.total_spending

    def as_dataframe(self):
        port_returns_ops = self.port_returns_list
        port_prespend_ops = self.port_prespend_list
        spend_amts_ops = self.spend_amts_list
        port_end_vals_ops = self.port_end_vals_list

        port_returns = self.model.sess.run(port_returns_ops)
        port_prespend = self.model.sess.run(port_prespend_ops)
        spend_amts = self.model.sess.run(spend_amts_ops)
        port_end_vals = self.model.sess.run(port_end_vals_ops)

        retlist = []
        for ix in range(self.model.ret_years):
            retlist.append([port_returns[ix],
                            port_prespend[ix],
                            spend_amts[ix],
                            port_end_vals[ix]
                        ])

        years = range(self.cohort_start_year, self.cohort_start_year+self.model.ret_years)
        return pd.DataFrame(retlist,
                            index = years,
                            columns=['portreturn', 'prespend', 'spend_amt', 'end_val'])

In [5]:
class CohortHistory:
    """represents a set of cohorts retiring in different years using a strategy,
    to enabling aggregating and summarizing their experiences"""
    def __init__(self, model, cohort_years = None):
        self.model = model
        if cohort_years is None:
            cohort_years = [year for year in range(self.model.first_year,
                                                   self.model.first_year + self.model.ret_cohorts)]

        print('%s Create cohort history, years %d to %d' % (strftime("%H:%M:%S"),
                                                            cohort_years[0], cohort_years[-1]))
        self.cohort_list = [Cohort(model, year) for year in cohort_years]
        self.total_spending_ops = [cohort.total_spending for cohort in self.cohort_list]

    def as_dataframe(self):
        """report on on each cohort by year, e.g. 1928"""
        total_spending_ops = [cohort.total_spending for cohort in self.model.cohort_history.cohort_list]
        mean_spending_ops = [cohort.mean_spending for cohort in self.model.cohort_history.cohort_list]
        sd_spending_ops = [cohort.sd_spending for cohort in self.model.cohort_history.cohort_list]
        min_spending_ops = [cohort.min_spending for cohort in self.model.cohort_history.cohort_list]
        max_spending_ops = [cohort.max_spending for cohort in self.model.cohort_history.cohort_list]
        ce_ops = [cohort.ce for cohort in self.model.cohort_history.cohort_list]

        retlist = []
        years = range(self.model.first_year, self.model.first_year + self.model.ret_cohorts)
        for year, \
            meanspend, \
            sdspend, \
            minspend, \
            maxspend, \
            totalspend, \
            ce in zip(years, self.model.sess.run(mean_spending_ops),
                      self.model.sess.run(sd_spending_ops),
                      self.model.sess.run(min_spending_ops),
                      self.model.sess.run(max_spending_ops),
                      self.model.sess.run(total_spending_ops),
                      self.model.sess.run(ce_ops)):

            retlist.append([meanspend, sdspend, minspend, maxspend, totalspend, ce])

        return pd.DataFrame(retlist, index = years,
                            columns=['mean_spend', 'sd_spend', 'min_spend', 'max_spend', 
                                     'total_spend', 'ce'])

    def spend_by_year(self):
        """report spending by year for each cohort (ret_years rows x num_cohorts)"""
        dataframes = [cohort.as_dataframe() for cohort in self.model.cohort_history.cohort_list]
        years = range(self.model.ret_years)
        cohorts = range(len(dataframes))

        retlist = []
        for ix in years:
            spendlist = [df.spend_amt.iloc[ix] for df in dataframes]
            retlist.append(spendlist)

        colnames = ["%d" % (cohort+self.model.first_year) for cohort in cohorts]
        return pd.DataFrame(retlist, index = years, columns=colnames)

    def returns_by_year(self):
        """report returns by year for each cohort (ret_years rows x num_cohorts)"""
        dataframes = [cohort.as_dataframe() for cohort in self.model.cohort_history.cohort_list]
        years = range(self.model.ret_years)
        cohorts = range(len(dataframes))

        retlist = []
        for ix in years:
            returnlist = [df.portreturn.iloc[ix] for df in dataframes]
            retlist.append(returnlist)

        colnames = ["%d" % (cohort+self.model.first_year) for cohort in cohorts]
        return pd.DataFrame(retlist, index = years, columns=colnames)

    def summarize_by_year(self):
        """report on outcomes by retirement year, e.g. retirement year 1, 2...30"""
        dataframes = [cohort.as_dataframe() for cohort in self.model.cohort_history.cohort_list]
        years = range(self.model.ret_years)
        retlist = []
        for ix in years:
            spendlist = np.array([df.spend_amt.iloc[ix] for df in dataframes])
            spend_mean = np.mean(spendlist)
            spend_sd = np.std(spendlist)
            spend_min = np.min(spendlist)
            spend_max = np.max(spendlist)
            retlist.append([spend_mean, spend_sd, spend_min, spend_max])

        return pd.DataFrame(retlist, index = years,
                            columns=['spend_mean', 'spend_sd', 'spend_min', 'spend_max'])

In [6]:
# Optimizer
# Create an op which is the sum of spending in all years
#  - negate it so it will be minimized
#  - add large penalty when a stock allocation is < 0 as a soft constraint
#  - add large penalty when a stock allocation is > 1 as a soft constraint
#  - add large penalty when const or var spencint is < 0 as a soft constraint
#  - result is an op which can be minimized by gradient descent

class CohortHistoryOptimize():
    def __init__(self, model):
        self.model = model
        self.best_objective = 0.0
        self.best_step = 0

        graph = self.model.graph

        with graph.as_default():
          with tf.device("/cpu:0"):

            print ('%s Create optimizer class' % strftime("%H:%M:%S"))
            print ('%s Run variable initializers' % strftime("%H:%M:%S"))
            self.model.sess.run(model.init_op)

            print('%s Create cost ops' % strftime("%H:%M:%S"))
            print('%s Sum %d ce ops' % (strftime("%H:%M:%S"), len(self.model.cohort_history.cohort_list)))
            ce_ops = [cohort.ce for cohort in self.model.cohort_history.cohort_list]
            ce_tensor = tf.pack(ce_ops, name="all_cohorts_ce_tensor")

            # ce over ret_cohorts years
            self.total_ce_op = self.model.gen_ce(ce_tensor, "all_cohorts_ce")

            print("%s Total CE spend, all cohorts: %f" % (strftime("%H:%M:%S"),
                                                          self.model.sess.run(self.total_ce_op)))
            # basic cost
            cost_op_1 = tf.neg(self.total_ce_op, name="basic_cost")
            print("%s Raw cost objective: %f" % (strftime("%H:%M:%S"), self.model.sess.run(cost_op_1)))

            cost_op_2 = tf.add(cost_op_1, model.cost_alloc_min_0_op, name="cost_add_alloc_min_0")
            print("%s Add soft constraint penalty if stock alloc < 0: %f" %
                  (strftime("%H:%M:%S"), self.model.sess.run(cost_op_2)))

            cost_op_3 = tf.add(cost_op_2, model.cost_alloc_max_1_op, name="cost_add_alloc_max_1")
            print("%s Add soft constraint penalty if stock alloc > 1: %f" % 
                  (strftime("%H:%M:%S"), self.model.sess.run(cost_op_3)))

            cost_op_4 = tf.add(cost_op_3, model.cost_vspend_min_0_op, name="cost_vspend_min_0")
            print("%s Add soft constraint penalty if var spending < 0: %f" % 
                  (strftime("%H:%M:%S"), self.model.sess.run(cost_op_4)))

            cost_op_5 = tf.add(cost_op_4, model.cost_cspend_min_0_op, name="cost_cspend_min_0")
            print("%s Add soft constraint if const spending < 0: %f" % 
                  (strftime("%H:%M:%S"), self.model.sess.run(cost_op_5)))

            self.cost_op = tf.add(cost_op_5, model.cost_alloc_decrease_op, name="cost_alloc_decrease")
            print("%s Add soft constraint if stock alloc increases in any year: %f" %
                  (strftime("%H:%M:%S"), self.model.sess.run(self.cost_op)))

            self.best_objective = -self.model.sess.run(self.cost_op)
            print("%s All inclusive objective to be minimized: %f" % (strftime("%H:%M:%S"),
                                                                      -self.best_objective))
            self.best_const_spend = self.model.sess.run(model.const_spending_op)
            self.best_var_spend = self.model.sess.run(model.var_spending_ops)
            self.best_stock_alloc = self.model.sess.run(model.allocation_ops[0])

    def run_step(self, report_steps=1):
        """run one step of optimizer
           calc gradients
           apply gradients * learning rate to each variable to descend gradient and improve objective
           increment global step to remember how many steps we've run
           if (hopefully) new objective is best to date, save params and objective"""

        _, step = self.model.sess.run([self.optimize_step,
                                       self.model.increment_step])
        self.steps_ago +=1
        cost = self.model.sess.run(self.cost_op)
        assert not(np.isnan(cost)), "Objective is nan"
        objective = - cost

        #print objective each step
        #print("objective %f best %f" %(objective, self.best_objective))
        if np.isnan(cost):
            sys.stdout.write('X')
            sys.stdout.flush()
        elif objective > self.best_objective:
            self.best_objective = objective
            self.best_const_spend = self.model.sess.run(model.const_spending_op)
            self.best_var_spend = self.model.sess.run(model.var_spending_ops)
            self.best_stock_alloc = self.model.sess.run(model.allocation_ops[0])
            self.best_step = step
            self.steps_ago = 0
            sys.stdout.write('!')
            sys.stdout.flush()
        else:
            sys.stdout.write('.')
            sys.stdout.flush()

        if step % report_steps == 0:
            sys.stdout.write("\n%s step %d objective %f best %f (%d steps ago)\n" %
                             (strftime("%H:%M:%S"),
                              step,
                              objective,
                              self.best_objective,
                              self.steps_ago))
            # print variables optimized and gradients for debugging
            # sys.stdout.write("\n")
            # var_vals = self.model.sess.run(self.model.all_var_ops)
            # print("%s Variables" % strftime("%H:%M:%S"))
            # print(var_vals)

            # grad_vals = self.model.sess.run([grad[0] for grad in self.grads])
            # print("%s Gradients" % strftime("%H:%M:%S"))
            # print(grad_vals)

            sys.stdout.flush()

            self.best_bond_alloc = pd.Series([1 - bsa for bsa in self.best_stock_alloc])
            pickle_list = [self.best_const_spend, self.best_var_spend, self.best_stock_alloc,
                           self.best_bond_alloc]
            pickle.dump( pickle_list, open( picklefile, "wb" ) )

            # every 10 report_steps show current best
            if step % (report_steps * 10) == 0:
                print ("\n#Objective: %f\n" % self.best_objective)
                print ("const_spend = %f" % self.best_const_spend)
                print ("var_spend_pcts = pd.Series(%s)" % str(self.best_var_spend))
                print ("stock_allocations = pd.Series(%s)\n" %str(self.best_stock_alloc))

    def optimize(self, learning_rate, steps):
        """create the op for the optimizer using specified learning_rate, run for specified steps"""
        self.learning_rate = learning_rate
        self.steps = steps
        self.steps_ago = 0 # how many steps since objective improved

        print("%s Objective: %f" % (strftime("%H:%M:%S"), self.best_objective))
        print("%s Constant spending: %f" % (strftime("%H:%M:%S"), self.best_const_spend))
        print("%s Variable spending by year" % strftime("%H:%M:%S"))
        print(self.best_var_spend)
        print("%s Stock allocation by year" % strftime("%H:%M:%S"))
        print(self.best_stock_alloc)

        with self.model.graph.as_default():
            with tf.device("/cpu:0"):

                # minimize op
                print('%s Create optimizer (learning rate %.12f)' % (strftime("%H:%M:%S"),
                                                                     self.learning_rate))
                self.optimizer = tf.train.GradientDescentOptimizer(self.learning_rate)
                self.grads = self.optimizer.compute_gradients(self.cost_op)
                self.optimize_step = self.optimizer.apply_gradients(self.grads)
                # following line is equivalent to previous 2 lines
                # self.optimize_step = self.optimizer.minimize(self.cost_op)

        print('%s Create optimizer op and run %d steps' % (strftime("%H:%M:%S"), self.steps))
        for i in range(self.steps):
            self.run_step()

In [7]:
learning_rate = 0.0000001
fileprefix = "gamma_8"
picklefile = '%s.pickle' % fileprefix
csvfile = "summary_%s.csv" % fileprefix
yearsfile = "years_%s.csv" % fileprefix
returnsfile = "retyears_%s.csv" % fileprefix


max_steps_unimproved = 500
gamma = 8

print('%s Start optimization session' % strftime("%H:%M:%S"))
print('%s learning_rate: %.12f steps %d picklefile %s' % (strftime("%H:%M:%S"), learning_rate, max_steps_unimproved, picklefile))

#print("opening picklefile %s" % picklefile)
#const_spend, var_spend_pcts, stock_allocations, bond_allocations  = pickle.load( open(picklefile, "rb" ) )

print ("const spend_pct: %f" % const_spend_pct)
print ("variable spend:")
print (var_spend_pcts)
print ("stock allocation:" )
print (stock_allocations)

model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = gamma,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

print('%s Summary by cohort' % strftime("%H:%M:%S"))
print(model.cohort_history.as_dataframe())
all_years = model.cohort_history.spend_by_year()
all_years.to_csv(yearsfile, format="%.18f")
ret_years = model.cohort_history.returns_by_year()
ret_years.to_csv(returnsfile, format="%.18f")

summary = model.cohort_history.summarize_by_year()
print(summary)
summary.to_csv(csvfile, format="%.18f")

# run optimizer
# set an initial learning rate that improves objective by a reasonable amount each step

model.optimizer = CohortHistoryOptimize(model)
model.optimizer.optimize(learning_rate, steps=1)

# continue optimizing without re-initializing vars or optimizer
# reduce learning rate if no improvement for a while
# end when learning rate is too small to make significant improvement

max_steps = 1001 # don't run for hours because notebook will eventually crash
report_steps = 50
learning_rate = model.optimizer.learning_rate

for i in range(max_steps):
    model.optimizer.run_step(report_steps=report_steps)
    if model.optimizer.steps_ago >= max_steps_unimproved: # no improvement for too long
        break

const_spend = model.optimizer.best_const_spend
var_spend_pcts = pd.Series(model.optimizer.best_var_spend)
stock_allocations = pd.Series(model.optimizer.best_stock_alloc)
bond_allocations   = 1 - stock_allocations
pickle_list = [const_spend, var_spend_pcts, stock_allocations, bond_allocations]
pickle.dump( pickle_list, open( picklefile, "wb" ) )


16:18:21 Start optimization session
16:18:21 learning_rate: 0.000000100000 steps 500 picklefile gamma_8.pickle
const spend_pct: 0.022500
variable spend:
0     0.016667
1     0.017241
2     0.017857
3     0.018519
4     0.019231
5     0.020000
6     0.020833
7     0.021739
8     0.022727
9     0.023810
10    0.025000
11    0.026316
12    0.027778
13    0.029412
14    0.031250
15    0.033333
16    0.035714
17    0.038462
18    0.041667
19    0.045455
20    0.050000
21    0.055556
22    0.062500
23    0.071429
24    0.083333
25    0.100000
26    0.125000
27    0.166667
28    0.250000
29    1.000000
dtype: float64
stock allocation:
0     0.820
1     0.815
2     0.810
3     0.805
4     0.800
5     0.795
6     0.790
7     0.785
8     0.780
9     0.775
10    0.770
11    0.765
12    0.760
13    0.755
14    0.750
15    0.745
16    0.740
17    0.735
18    0.730
19    0.725
20    0.720
21    0.715
22    0.710
23    0.705
24    0.700
25    0.695
26    0.690
27    0.685
28    0.680
29    0.675
dtype: float64
16:18:21 Create TensorFlow graph and session
16:18:21 Create constants stocks_return
16:18:21 Create constants bonds_return
16:18:21 Create variables stocks_alloc
16:18:21 Create ops for soft constraint alloc_min_0 > 0
16:18:21 Create ops for soft constraint alloc_max_1 < 1
16:18:21 Create ops for soft constraint, declining stock alloc alloc_decrease
16:18:21 Create ops for bonds_alloc
16:18:21 Create variables var_spend
16:18:22 Create ops for soft constraint vspend_min_0 > 0
16:18:22 Create cohort history, years 1928 to 1986
16:18:30 Summary by cohort
      mean_spend   sd_spend  min_spend   max_spend  total_spend        ce
1928    7.103008   6.346463   3.609243   37.559135   213.090225  4.643645
1929    4.931037   2.850893   3.220363   18.891738   147.931102  3.927663
1930    5.925849   4.854840   3.295279   29.700473   177.775481  4.117642
1931    7.924386   8.582596   3.465120   49.803788   237.731583  4.503485
1932   13.610733  20.455730   3.996189  115.859799   408.321990  5.587242
1933   13.847575  19.429963   4.570448  108.115107   415.427255  5.790120
1934    9.919830  12.572192   3.899031   71.789912   297.594889  4.851973
1935   11.235387  15.423238   4.081657   87.596596   337.061613  5.020730
1936    8.024943   9.134418   3.547464   52.765278   240.748289  4.284851
1937    5.975075   4.725102   3.253967   27.792250   179.252244  3.845672
1938   11.590917  15.091473   3.788194   84.922410   347.727496  4.761143
1939    8.999374  10.204098   3.442112   58.120192   269.981217  4.256801
1940    9.476911   9.838279   3.449368   54.558966   284.307321  4.314835
1941   11.468200  12.702490   3.570599   70.703614   344.046006  4.581398
1942   16.791024  21.436332   4.034989  119.916818   503.730721  5.360601
1943   16.966374  22.265720   4.217314  125.724308   508.991217  5.394454
1944   14.330796  15.489738   3.987815   85.468515   429.923870  5.105430
1945   12.222867   9.913253   3.775813   51.797070   366.686011  4.911167
1946    9.603656   7.283004   3.436568   41.698092   288.109680  4.456034
1947   14.384591  13.997187   3.843206   81.213358   431.537743  5.234968
1948   15.646485  14.076276   3.950559   80.770604   469.394562  5.645485
1949   15.852519  13.379560   4.215325   77.623040   475.575580  6.040702
1950   13.773480  10.645199   4.240157   63.358408   413.204408  5.975190
1951   11.935441   8.866539   4.139393   54.428916   358.063222  5.817038
1952   10.697537   6.745577   4.134621   42.166271   320.926105  5.819612
1953    9.793289   6.527503   3.900066   42.029831   293.798655  5.778506
1954   10.492669   7.875785   4.657248   50.182212   314.780058  6.491792
1955    7.422228   4.141696   4.186486   28.060199   222.666846  5.528992
1956    6.114897   3.000229   3.779466   21.364709   183.446907  4.970847
1957    6.335115   3.810415   3.745804   25.857817   190.053451  5.029950
1958    7.563460   5.607652   4.478531   35.657479   226.903813  5.661326
1959    5.707103   2.884216   4.044702   20.069798   171.213076  4.797619
1960    5.587263   3.113681   3.933558   21.367040   167.617891  4.655229
1961    5.789478   3.284632   4.021656   21.814771   173.684353  4.721407
1962    4.805286   1.943842   3.627647   14.397159   144.158577  4.198353
1963    5.698170   3.630369   3.787890   23.631695   170.945089  4.499244
1964    4.886873   2.291664   3.516002   16.010685   146.606204  4.112266
1965    4.288050   1.201389   3.339557    9.565156   128.641501  3.842409
1966    4.047540   1.009610   3.237003    8.655287   121.426207  3.690863
1967    5.308811   3.653913   3.403824   23.385678   159.264336  4.017470
1968    4.565352   2.489663   3.220999   16.893296   136.960556  3.719474
1969    4.678672   3.063094   3.174779   19.978332   140.360170  3.659766
1970    6.910801   8.185465   3.350516   47.763073   207.324042  4.009532
1971    7.532961   8.958852   3.351297   50.877442   225.988815  4.044391
1972    6.879558   6.971568   3.247201   39.305564   206.386744  3.889541
1973    5.921728   4.622463   3.134198   26.011757   177.651843  3.706021
1974    9.354891  10.972260   3.363042   61.334286   280.646717  4.153282
1975   17.872324  26.009169   4.299000  143.101514   536.169727  5.262786
1976   14.590081  19.457642   3.888588  107.437121   437.702437  4.808366
1977   12.856045  16.620105   3.627015   93.068867   385.681349  4.517707
1978   16.532428  22.528160   3.853002  125.330292   495.972849  4.928438
1979   17.773045  21.145701   3.938639  113.758951   533.191350  5.196985
1980   18.634246  22.891358   3.915696  127.109058   559.027393  5.356443
1981   17.419744  21.721699   3.728327  122.761078   522.592316  5.221420
1982   21.479058  27.699562   4.230340  156.191132   644.371746  5.940978
1983   19.112005  24.656700   4.168397  140.470198   573.360148  5.776562
1984   17.957604  24.395430   3.976096  140.142303   538.728127  5.682128
1985   19.044310  27.139533   4.357432  155.328991   571.329286  6.020979
1986   15.596985  20.404056   4.206998  116.379295   467.909551  5.593521
    spend_mean   spend_sd  spend_min   spend_max
0     4.030408   0.298318   3.363042    4.657248
1     4.129680   0.437146   3.134198    5.285261
2     4.229364   0.522057   3.220363    5.401103
3     4.354657   0.657713   3.244310    6.217383
4     4.491963   0.756426   3.253967    6.378949
5     4.649857   0.865255   3.174779    6.886115
6     4.818928   1.036802   3.220999    7.744866
7     5.012877   1.211305   3.316392    7.947275
8     5.188269   1.328391   3.177389    8.098496
9     5.398312   1.527354   3.273882    9.219797
10    5.670773   1.777047   3.254359    9.830153
11    5.976920   2.071741   3.308724   10.264965
12    6.327206   2.408398   3.186369   11.858104
13    6.728654   2.831714   3.239022   13.537225
14    7.144893   3.224169   3.386575   15.083094
15    7.567687   3.609645   3.237003   16.112155
16    8.030148   4.051538   3.364358   19.321692
17    8.531063   4.419612   3.481640   21.741202
18    9.046367   4.776763   3.489725   21.351997
19    9.605672   5.098866   3.680113   21.965542
20   10.253281   5.447986   3.994624   22.581419
21   10.993393   5.870603   3.937472   22.499508
22   11.788984   6.358721   4.062717   24.836824
23   12.703260   6.957020   4.424327   27.235292
24   13.721394   7.505990   4.292667   31.673663
25   14.984045   8.078475   4.727285   34.160823
26   16.695166   8.775084   4.885216   34.730011
27   19.310819  10.301901   5.180413   40.852131
28   23.972998  13.573035   5.215851   53.220372
29   65.383673  42.556422   8.655287  156.191132
16:18:51 Create optimizer class
16:18:51 Run variable initializers
16:18:51 Create cost ops
16:18:51 Sum 59 ce ops
16:18:51 Total CE spend, all cohorts: 4.458856
16:18:52 Raw cost objective: -4.458856
16:18:53 Add soft constraint penalty if stock alloc < 0: -4.458856
16:18:54 Add soft constraint penalty if stock alloc > 1: -4.458856
16:18:55 Add soft constraint penalty if var spending < 0: -4.458856
16:18:56 Add soft constraint if const spending < 0: -4.458856
16:18:56 Add soft constraint if stock alloc increases in any year: -4.458856
16:18:57 All inclusive objective to be minimized: -4.458856
16:18:57 Objective: 4.458856
16:18:57 Constant spending: 2.250000
16:18:57 Variable spending by year
[0.016666666666666666, 0.017241379310344827, 0.017857142857142856, 0.018518518518518517, 0.019230769230769232, 0.02, 0.020833333333333332, 0.021739130434782608, 0.022727272727272728, 0.023809523809523808, 0.025000000000000001, 0.026315789473684209, 0.027777777777777776, 0.029411764705882353, 0.03125, 0.033333333333333333, 0.035714285714285712, 0.038461538461538464, 0.041666666666666664, 0.045454545454545456, 0.050000000000000003, 0.055555555555555552, 0.0625, 0.071428571428571425, 0.083333333333333329, 0.10000000000000001, 0.125, 0.16666666666666666, 0.25, 1.0]
16:18:57 Stock allocation by year
[0.81999999999999995, 0.81499999999999995, 0.80999999999999994, 0.80499999999999994, 0.79999999999999993, 0.79499999999999993, 0.78999999999999992, 0.78499999999999992, 0.77999999999999992, 0.77499999999999991, 0.76999999999999991, 0.7649999999999999, 0.76000000000000001, 0.75499999999999989, 0.75, 0.745, 0.73999999999999999, 0.73499999999999999, 0.72999999999999998, 0.72499999999999998, 0.71999999999999997, 0.71499999999999997, 0.70999999999999996, 0.70499999999999996, 0.69999999999999996, 0.69499999999999995, 0.68999999999999995, 0.68499999999999994, 0.67999999999999994, 0.67499999999999993]
16:18:57 Create optimizer (learning rate 0.000000100000)
16:20:32 Create optimizer op and run 1 steps
!
16:21:23 step 1 objective 4.459686 best 4.459686 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:21:38 step 50 objective 4.499532 best 4.499532 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:21:54 step 100 objective 4.538456 best 4.538456 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:22:11 step 150 objective 4.575187 best 4.575187 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:22:27 step 200 objective 4.607514 best 4.607514 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:22:42 step 250 objective 4.623554 best 4.623554 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:22:57 step 300 objective 4.624527 best 4.624527 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:23:13 step 350 objective 4.624715 best 4.624715 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:23:28 step 400 objective 4.624901 best 4.624901 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:23:43 step 450 objective 4.625085 best 4.625085 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:23:58 step 500 objective 4.625268 best 4.625268 (0 steps ago)

#Objective: 4.625268

const_spend = 2.451341
var_spend_pcts = pd.Series([0.016885938821778582, 0.01744496319734206, 0.01803798051766901, 0.018676430579853142, 0.019364949417501696, 0.020113948308229799, 0.020932209454706991, 0.021829421748703906, 0.022808914589830288, 0.023886301870811314, 0.025072995942117225, 0.026384768422203237, 0.027841642267812355, 0.029466383421367068, 0.031294196989241158, 0.033368199547562351, 0.035742289444631065, 0.038484683060579776, 0.041684408572901441, 0.04546682421643148, 0.050005601872476546, 0.055556622020647067, 0.062498066119545316, 0.071423684280120026, 0.08332717502673688, 0.099993384970494051, 0.1249932477807101, 0.16666038757204479, 0.24999485538097657, 1.0])
stock_allocations = pd.Series([0.81998746963760849, 0.81500684602996243, 0.81000281757623882, 0.80499918148079042, 0.79999333428579267, 0.795000800941847, 0.79000908348439425, 0.78499215545818135, 0.77997990794511596, 0.7750104069161895, 0.77000440513098645, 0.76499787527758523, 0.76000301231652967, 0.75500717506908532, 0.75001152412563798, 0.7449975407879833, 0.73999795360963416, 0.7350039898494819, 0.72999913738754352, 0.72500117456291902, 0.7199994839302033, 0.71500196482479683, 0.7100012246455113, 0.70500154990072972, 0.69999897666341693, 0.69500115614223379, 0.68999991284851314, 0.68499977138979828, 0.6800004321246198, 0.67500022181664676])

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:24:14 step 550 objective 4.625450 best 4.625450 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:24:29 step 600 objective 4.625631 best 4.625631 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:24:44 step 650 objective 4.625810 best 4.625810 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:24:59 step 700 objective 4.625988 best 4.625988 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:25:14 step 750 objective 4.626166 best 4.626166 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:25:30 step 800 objective 4.626342 best 4.626342 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:25:46 step 850 objective 4.626517 best 4.626517 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:26:02 step 900 objective 4.626691 best 4.626691 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:26:18 step 950 objective 4.626864 best 4.626864 (0 steps ago)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16:26:35 step 1000 objective 4.627035 best 4.627035 (0 steps ago)

#Objective: 4.627035

const_spend = 2.449279
var_spend_pcts = pd.Series([0.017031358899133917, 0.017578907100408648, 0.018153913176439, 0.018773929350320188, 0.019442976867359504, 0.02017525195237083, 0.020981426219935866, 0.021872245874610469, 0.022845390511961917, 0.02392083410972845, 0.025106346990624251, 0.026416451908836715, 0.027870825712112704, 0.029489354583179927, 0.031309749019701627, 0.033376993023987174, 0.035747084176808781, 0.03848709198815143, 0.041683566864681783, 0.045462594615727062, 0.049996426607145449, 0.055544269217236403, 0.062484086065450327, 0.07140802044810142, 0.08331130267543374, 0.099978316063519923, 0.12497911719401604, 0.16664787964068284, 0.24998475062647796, 1.0])
stock_allocations = pd.Series([0.81996814522394967, 0.81502364606846311, 0.81000899621286016, 0.80499687079078441, 0.79998105161942878, 0.79500299530143004, 0.79002391116876625, 0.78497588772189675, 0.779944216924226, 0.77503013324606707, 0.77001127432726146, 0.7649929216969743, 0.76000826860455695, 0.75501934400792037, 0.7500313158798908, 0.74499186110935423, 0.73999366868905647, 0.7350112742117203, 0.72999707180587459, 0.72500307966269051, 0.7199982976167042, 0.71500541865015044, 0.71000339566623971, 0.70500432732090546, 0.69999703658508805, 0.69500330560163537, 0.68999972551609978, 0.68499933506319477, 0.68000124374357573, 0.67500063694877588])

!!

In [8]:
# This is a good solution, objective ~4.7 after running many times at decreasing learning rates

const_spend_pct = 0.018033736915
var_spend_pcts = pd.Series([0.027420789200650691, 0.028952687799914788, 0.029987640366612587, 0.030563608815231274, 0.0314680099541\
00477, 0.0327892165458714, 0.033812194391987072, 0.035195378137199251, 0.037327252304951188, 0.039532656754857148, 0.04171133948016\
5768, 0.044205216561892961, 0.047114695564658048, 0.049375968207551808, 0.051345728157015109, 0.05421706340461499, 0.05717491967092\
8206, 0.060803826336315064, 0.065058176584335367, 0.069765865219855477, 0.073504852805054741, 0.079201426822341936, 0.0872191175663\
84868, 0.095745358220279422, 0.10919882479478223, 0.12781818289102989, 0.15488310518171028, 0.20138642112957703, 0.2991516975921870\
2, 1.0])
stock_allocations = pd.Series([0.85995594522555152, 0.85919296367537157, 0.85473424640534235, 0.85205353144872398, 0.84754303090697\
891, 0.84731613048078414, 0.84674123511449917, 0.83207329884146453, 0.81886882987838872, 0.81886882879077094, 0.81787428397721995, \
0.81082294011238876, 0.80935846467721007, 0.80915247713021987, 0.80869082405588244, 0.79235039525216477, 0.78697518924653442, 0.785\
66021347854842, 0.77606352982164872, 0.77226108408967586, 0.76418701504255893, 0.75957408240400748, 0.7506635662237392, 0.742109947\
15180123, 0.68810212447786767, 0.68212008875931995, 0.66794616058167322, 0.66270963157142881, 0.65979617905489518, 0.61499306500466\
799])

model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 8.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

print('%s Summary by cohort' % strftime("%H:%M:%S"))
print(model.cohort_history.as_dataframe())


16:26:35 Create TensorFlow graph and session
16:26:35 Create constants stocks_return
16:26:35 Create constants bonds_return
16:26:35 Create variables stocks_alloc
16:26:36 Create ops for soft constraint alloc_min_0 > 0
16:26:36 Create ops for soft constraint alloc_max_1 < 1
16:26:36 Create ops for soft constraint, declining stock alloc alloc_decrease
16:26:36 Create ops for bonds_alloc
16:26:36 Create variables var_spend
16:26:36 Create ops for soft constraint vspend_min_0 > 0
16:26:36 Create cohort history, years 1928 to 1986
16:26:54 Summary by cohort
      mean_spend   sd_spend  min_spend  max_spend  total_spend        ce
1928    6.231747   2.586673   3.915274  17.821277   186.952419  5.084862
1929    4.349565   0.782169   3.328030   7.501703   130.486944  4.030072
1930    5.168108   1.878688   3.481071  13.678188   155.043250  4.372770
1931    6.942199   4.156061   3.757318  25.990977   208.265967  4.988632
1932   11.907819  10.995819   4.657191  64.203340   357.234583  6.717394
1933   12.409915  10.975671   5.707059  62.326652   372.297460  7.122779
1934    8.770417   6.634870   4.493893  39.331985   263.112496  5.597494
1935    9.994297   8.516653   4.582492  49.637098   299.828902  5.872640
1936    7.244735   5.007556   3.797345  30.000082   217.342042  4.749586
1937    5.573977   2.641160   3.355214  16.233739   167.219303  4.108138
1938   10.701712   8.883767   4.258534  50.196110   321.051347  5.540202
1939    8.391608   5.971544   3.746160  34.362074   251.748240  4.759495
1940    9.036516   6.006617   3.778838  32.865000   271.095473  4.849111
1941   11.086558   7.879285   3.964307  43.227815   332.596728  5.263955
1942   16.164608  13.065494   4.758616  72.484401   484.938254  6.510528
1943   16.113163  13.026037   5.064820  73.730542   483.394904  6.552541
1944   14.097082   9.526212   4.634658  51.871333   422.912468  6.064975
1945   12.434738   6.464632   4.310231  31.720537   373.042153  5.734365
1946    9.730501   4.468377   3.776355  24.239083   291.915035  4.983624
1947   14.496389   8.142194   4.429282  47.040961   434.891662  6.196933
1948   16.054581   8.269578   4.605323  47.037054   481.637441  6.816020
1949   16.512693   7.756138   5.051763  44.796038   495.380795  7.488450
1950   14.339591   5.889525   5.110959  35.315021   430.187732  7.416719
1951   12.263212   4.507644   4.938166  28.548033   367.896358  7.126248
1952   11.099920   3.458155   4.955773  22.376416   332.997593  7.077327
1953   10.079154   3.117718   4.512278  21.909049   302.374628  6.916262
1954   10.719456   3.528773   5.817874  25.765312   321.583693  8.216409
1955    7.500280   1.732817   5.024713  13.625401   225.008394  6.480559
1956    6.041473   1.166329   4.342685   9.323010   181.244201  5.402160
1957    6.188017   1.430991   4.245431  11.928845   185.640500  5.469803
1958    7.400427   2.256392   5.481788  17.755333   222.012800  6.468176
1959    5.492784   0.962050   4.153474   8.606843   164.783518  5.017501
1960    5.266278   0.969556   3.957972   8.831636   157.988334  4.805155
1961    5.504569   1.105926   4.127639   9.953767   165.137075  5.002865
1962    4.545435   0.611007   3.599994   5.653342   136.363045  4.259504
1963    5.304891   1.239745   3.995944  10.804729   159.146741  4.785339
1964    4.536707   0.666107   3.561182   6.649972   136.101216  4.230332
1965    4.034187   0.464631   3.317734   4.933661   121.025608  3.853209
1966    3.798965   0.436429   3.056128   4.761497   113.968938  3.631914
1967    4.780540   1.299114   3.532308  10.670690   143.416193  4.242993
1968    4.080396   0.698098   3.262251   6.921534   122.411883  3.800844
1969    4.154946   1.052271   3.220717   8.984043   124.648373  3.763955
1970    5.851634   3.778231   3.500764  23.740204   175.549007  4.308289
1971    6.439688   4.524693   3.540727  26.899853   193.190649  4.382039
1972    6.026275   3.644721   3.416846  21.269286   180.788262  4.169551
1973    5.387303   2.587010   3.239628  14.836384   161.619096  3.909019
1974    8.366438   6.280137   3.604047  34.799294   250.993126  4.589014
1975   15.811972  15.157204   5.095055  81.969083   474.359153  6.370086
1976   13.054721  11.353804   4.432720  61.566076   391.641627  5.612846
1977   11.563246   9.477082   4.030354  52.439163   346.897375  5.147522
1978   14.809805  12.776853   4.443772  70.262519   444.294141  5.792671
1979   16.133861  12.211305   4.620209  64.034778   484.015841  6.202600
1980   16.807785  12.534120   4.597383  68.041751   504.233557  6.417250
1981   15.636378  11.645574   4.221454  65.610214   469.091349  6.105409
1982   19.516503  15.047263   5.047948  85.260856   585.495079  7.288942
1983   17.418456  13.042985   4.980582  75.605946   522.553678  7.058821
1984   16.115919  12.403931   4.634923  73.494326   483.477558  6.852123
1985   16.799990  13.359735   5.276680  79.663902   503.999714  7.478969
1986   13.877007  10.017654   5.044592  59.901777   416.310211  6.808008

In [9]:
df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.show()



In [10]:
df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


Out[10]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 1.803374 2.742079 85.995595 14.004405 4.740133 3.604047 5.817874
1 1.803374 2.895269 85.919296 14.080704 4.954235 3.239628 6.990869
2 1.803374 2.998764 85.473425 14.526575 5.102155 3.338954 7.131929
3 1.803374 3.056361 85.205353 14.794647 5.231042 3.328030 8.398137
4 1.803374 3.146801 84.754303 15.245697 5.401859 3.355214 8.518801
5 1.803374 3.278922 84.731613 15.268387 5.640160 3.220717 9.473607
6 1.803374 3.381219 84.674124 15.325876 5.843391 3.262251 10.780489
7 1.803374 3.519538 83.207330 16.792670 6.106355 3.482782 11.021653
8 1.803374 3.732725 81.886883 18.113117 6.406250 3.268162 11.260031
9 1.803374 3.953266 81.886883 18.113117 6.743933 3.372960 13.164527
10 1.803374 4.171134 81.787428 18.212572 7.143623 3.348702 14.116656
11 1.803374 4.420522 81.082294 18.917706 7.594771 3.409277 14.611927
12 1.803374 4.711470 80.935846 19.064154 8.120312 3.246629 17.223627
13 1.803374 4.937597 80.915248 19.084752 8.577191 3.273653 18.662278
14 1.803374 5.134573 80.869082 19.130918 8.939804 3.422242 20.469916
15 1.803374 5.421706 79.235040 20.764960 9.347002 3.205772 21.213132
16 1.803374 5.717492 78.697519 21.302481 9.727362 3.317734 24.702617
17 1.803374 6.080383 78.566021 21.433979 10.143674 3.450519 27.386865
18 1.803374 6.505818 77.606353 22.393647 10.531063 3.401243 25.689757
19 1.803374 6.976587 77.226108 22.773892 10.881431 3.542791 25.921611
20 1.803374 7.350485 76.418702 23.581298 11.038810 3.749375 25.284224
21 1.803374 7.920143 75.957408 24.042592 11.355122 3.645121 23.647565
22 1.803374 8.721912 75.066357 24.933643 11.747497 3.685455 24.225980
23 1.803374 9.574536 74.210995 25.789005 11.987134 3.848685 26.258617
24 1.803374 10.919882 68.810212 31.189788 12.405500 3.660422 29.108868
25 1.803374 12.781818 68.212009 31.787991 12.913144 3.878824 28.843689
26 1.803374 15.488311 66.794616 33.205384 13.578350 3.829700 27.946354
27 1.803374 20.138642 66.270963 33.729037 14.794993 3.828073 30.237365
28 1.803374 29.915170 65.979618 34.020382 17.323086 3.581394 37.396578
29 1.803374 100.000000 61.499307 38.500693 35.762553 3.056128 85.260856

In [11]:
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)

ax.plot(df_summary['const_spend'], label="constant spending")
ax.plot(df_summary['var_spend'], label="variable spending")
ax.plot(df_summary['stocks'], label="stocks")
ax.plot(df_summary['bonds'], label="bonds")

plt.ylabel("% of portfolio (starting portfolio for const spending)")
plt.xlabel("Retirement year")

#legend causes probs for plotly
#leg = ax.legend(loc="upper left", fancybox=True)
#leg.get_frame().set_alpha(0.3)

plt.show()



In [12]:
plot_url = py.iplot_mpl(fig, filename='chart')
plot_url


Out[12]:

In [13]:
from decimal import *
getcontext().prec = 32

def u (cashflow, gamma):
    d_cashflow = Decimal(cashflow)
    d_gamma = Decimal(gamma)
    d_one = Decimal('1.00000000000000000000000000000000')
    
    if gamma != 1:
        one_minus_gamma = d_one - d_gamma
        u = (d_cashflow ** one_minus_gamma - d_one) / one_minus_gamma
        return u
    else:
        return d_cashflow.ln()


x = np.linspace(0,10,1001)
y0 = np.array([u(z, 0.0) for z in x])
y1 = np.array([u(z, 1.0) for z in x])
y2 = np.array([u(z, 2.0) for z in x])
y4 = np.array([u(z, 4.0) for z in x])
y8 = np.array([u(z, 8.0) for z in x])
y16 = np.array([u(z, 16.0) for z in x])
y32 = np.array([u(z, 32.0) for z in x])

# note, at high gamma, can barely distinguish 9.9 from 10
# even with extra decimal places vs. float64
y32


Out[13]:
array([Decimal('-Infinity'),
       Decimal('-3.2258064516129011441382804407354E+60'),
       Decimal('-1.5021331848636740558493325676468E+51'), ...,
       Decimal('0.032258064516129032258064516129029'),
       Decimal('0.032258064516129032258064516129029'),
       Decimal('0.032258064516129032258064516129029')], dtype=object)

In [14]:
plt.axis([ -0.2, 6.2, -3.0, 3,])

plt.plot(x, y0, color='k', label = r'$\gamma = 0$')
plt.plot(x, y1, color='b', label = r'$\gamma = 1$')
plt.plot(x, y2, color='r', label = r'$\gamma = 2$')
plt.plot(x, y4, color='g', label = r'$\gamma = 4$')
plt.plot(x, y8, color='purple', label = r'$\gamma = 8$')
plt.plot(x, y16, color='cyan', label = r'$\gamma = 16$')
plt.plot(x, y32, color='orange', label = r'$\gamma = 32$')

plt.ylabel("CRRA utility")
plt.xlabel("Cash flow ")
plt.legend(loc="lower right", ncol=2, bbox_to_anchor=[1, 0])


Out[14]:
<matplotlib.legend.Legend at 0x7f972ffc4090>

In [15]:
# gamma = 0

const_spend_pct = 0.000000
var_spend_pcts = pd.Series([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
stock_allocations = pd.Series([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])

model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 0.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:27:28 Create TensorFlow graph and session
16:27:28 Create constants stocks_return
16:27:28 Create constants bonds_return
16:27:28 Create variables stocks_alloc
16:27:28 Create ops for soft constraint alloc_min_0 > 0
16:27:28 Create ops for soft constraint alloc_max_1 < 1
16:27:28 Create ops for soft constraint, declining stock alloc alloc_decrease
16:27:28 Create ops for bonds_alloc
16:27:29 Create variables var_spend
16:27:29 Create ops for soft constraint vspend_min_0 > 0
16:27:29 Create cohort history, years 1928 to 1986
Out[15]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
1 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
2 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
3 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
4 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
5 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
6 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
7 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
8 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
9 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
10 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
11 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
12 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
13 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
14 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
15 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
16 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
17 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
18 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
19 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
20 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
21 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
22 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
23 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
24 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
25 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
26 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
27 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
28 0.0 0.0 100.0 0.0 0.000000 0.000000 0.000000
29 0.0 100.0 100.0 0.0 824.468943 343.454154 1701.299092

In [16]:
# gamma = 1
#Objective: 9.308091

const_spend_pct = 0.000210991438
var_spend_pcts = pd.Series([0.033091538348404496, 0.034238168285506382, 0.035469702440122583, 0.036794378152860914, 0.038221756778048857, 0.03976288047503998, 0.041430947662087859, 0.043242686669740876, 0.045216839687742649, 0.047378561748440171, 0.04975812222825144, 0.052387429460561782, 0.055309651974463532, 0.058579229736857016, 0.062260304236144169, 0.066431027576827867, 0.071195216736901532, 0.076691471776952069, 0.083104783014253336, 0.090686442187799754, 0.099783243431139951, 0.11089638349654848, 0.12478101422016996, 0.14262272988447125, 0.16639735490582469, 0.19966642400185314, 0.24955549906491153, 0.33270239606716523, 0.49854063104140045, 1.0])
stock_allocations = pd.Series([0.99999999932534844, 0.99997650919314807, 0.99995061726681489, 0.99992473029357398, 0.99933647980151297, 0.99915824051192981, 0.99801222176142057, 0.9979130315525, 0.99767125200215423, 0.99762089448679248, 0.99748512848933513, 0.99744284142382889, 0.99734447386668101, 0.99709686697171851, 0.99700226658406843, 0.99670993806194219, 0.99638045411265985, 0.99628135924967531, 0.99610184313223671, 0.99588814803982129, 0.99574876269478818, 0.99557475808121887, 0.99508658081662371, 0.99501534744474973, 0.99468563684794653, 0.98122175151927205, 0.97403705686938635, 0.9563148825731429, 0.93175496886826248, 0.91435187362375314])


model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 1.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:28:03 Create TensorFlow graph and session
16:28:03 Create constants stocks_return
16:28:03 Create constants bonds_return
16:28:03 Create variables stocks_alloc
16:28:03 Create ops for soft constraint alloc_min_0 > 0
16:28:03 Create ops for soft constraint alloc_max_1 < 1
16:28:03 Create ops for soft constraint, declining stock alloc alloc_decrease
16:28:03 Create ops for bonds_alloc
16:28:03 Create variables var_spend
16:28:04 Create ops for soft constraint vspend_min_0 > 0
16:28:04 Create cohort history, years 1928 to 1986
Out[16]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 0.021099 3.309154 100.000000 6.746515e-08 3.597118 2.064910 5.094147
1 0.021099 3.423817 99.997651 2.349081e-03 3.853408 1.594728 6.730730
2 0.021099 3.546970 99.995062 4.938273e-03 4.108692 1.626304 7.032329
3 0.021099 3.679438 99.992473 7.526971e-03 4.432245 1.652666 9.042422
4 0.021099 3.822176 99.933648 6.635202e-02 4.781124 1.852413 9.191793
5 0.021099 3.976288 99.915824 8.417595e-02 5.182618 1.712576 11.383236
6 0.021099 4.143095 99.801222 1.987778e-01 5.612291 1.817066 13.699620
7 0.021099 4.324269 99.791303 2.086968e-01 6.102872 2.191610 14.312215
8 0.021099 4.521684 99.767125 2.328748e-01 6.552272 1.899995 14.695080
9 0.021099 4.737856 99.762089 2.379106e-01 7.090885 2.097695 17.605828
10 0.021099 4.975812 99.748513 2.514872e-01 7.779066 2.336610 19.423694
11 0.021099 5.238743 99.744284 2.557159e-01 8.553012 2.479428 20.127833
12 0.021099 5.530965 99.734447 2.655526e-01 9.428448 2.404902 24.226963
13 0.021099 5.857923 99.709687 2.903133e-01 10.432121 2.552166 24.875795
14 0.021099 6.226030 99.700227 2.997733e-01 11.429301 2.863576 28.956837
15 0.021099 6.643103 99.670994 3.290062e-01 12.382714 2.668600 30.429766
16 0.021099 7.119522 99.638045 3.619546e-01 13.395461 2.947546 35.091477
17 0.021099 7.669147 99.628136 3.718641e-01 14.438561 3.402328 40.921069
18 0.021099 8.310478 99.610184 3.898157e-01 15.441457 3.669829 35.900669
19 0.021099 9.068644 99.588815 4.111852e-01 16.483097 3.697111 42.143027
20 0.021099 9.978324 99.574876 4.251237e-01 17.636579 4.308710 44.343573
21 0.021099 11.089638 99.557476 4.425242e-01 18.824673 4.608361 43.230767
22 0.021099 12.478101 99.508658 4.913419e-01 19.891423 5.082879 46.466793
23 0.021099 14.262273 99.501535 4.984653e-01 20.979373 5.924389 51.172062
24 0.021099 16.639735 99.468564 5.314363e-01 21.901249 6.247370 57.736473
25 0.021099 19.966642 98.122175 1.877825e+00 22.775050 6.536648 53.348776
26 0.021099 24.955550 97.403706 2.596294e+00 23.727347 7.634730 56.503904
27 0.021099 33.270240 95.631488 4.368512e+00 24.794141 8.815078 48.446182
28 0.021099 49.854063 93.175497 6.824503e+00 26.041660 9.211319 54.611926
29 0.021099 100.000000 91.435187 8.564813e+00 27.371094 11.333554 55.578990

In [17]:
# gamma = 2

#Objective: 7.232013

const_spend_pct = 0.006633196051
var_spend_pcts = pd.Series([0.041651782115649724, 0.04298351383672569, 0.044307405688720158, 0.045644563126512307, 0.046935371105069679, 0.048256594099694364, 0.049876320240762828, 0.051675640492043345, 0.053725939553432397, 0.056085607010629972, 0.058531950342040277, 0.061306089393194735, 0.064303573280149864, 0.067511781568926968, 0.070864795122752594, 0.074644853152054047, 0.079354357195997025, 0.084874368957208082, 0.091405892559806134, 0.098722229027103148, 0.10732707056538419, 0.11779806370242611, 0.13130833549322588, 0.1483069190268195, 0.17052212615693424, 0.20140944216392345, 0.24724327540518604, 0.32492166343973394, 0.48312975250382012, 1.0])
stock_allocations = pd.Series([0.99128140981577162, 0.9890186164406034, 0.9889931746406847, 0.98899317444042156, 0.98897294365892874, 0.98892702392459531, 0.98017623383520169, 0.97434793556657606, 0.96296843154230494, 0.96175969914252735, 0.96164944382766815, 0.95833816100296942, 0.95788049761186367, 0.95748409443075122, 0.94935349502063981, 0.94491171918437289, 0.93545099514061691, 0.9305036848920506, 0.92943528021152122, 0.92302678056280996, 0.92027686390605101, 0.91845762356824345, 0.91278569545635735, 0.91231424523945759, 0.90537710239072444, 0.89686206541538249, 0.88878472270546383, 0.87921904341544299, 0.86565496875257597, 0.84597976461898239])



model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 2.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.axis([ 0, 30, 0, 60,])

plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:28:40 Create TensorFlow graph and session
16:28:40 Create constants stocks_return
16:28:40 Create constants bonds_return
16:28:41 Create variables stocks_alloc
16:28:41 Create ops for soft constraint alloc_min_0 > 0
16:28:41 Create ops for soft constraint alloc_max_1 < 1
16:28:41 Create ops for soft constraint, declining stock alloc alloc_decrease
16:28:41 Create ops for bonds_alloc
16:28:41 Create variables var_spend
16:28:42 Create ops for soft constraint vspend_min_0 > 0
16:28:42 Create cohort history, years 1928 to 1986
Out[17]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 0.66332 4.165178 99.128141 0.871859 5.161895 3.245959 7.030791
1 0.66332 4.298351 98.901862 1.098138 5.396760 2.617961 8.928972
2 0.66332 4.430741 98.899317 1.100683 5.607994 2.623224 9.141206
3 0.66332 4.564456 98.899317 1.100683 5.878749 2.597315 11.336183
4 0.66332 4.693537 98.897294 1.102706 6.146992 2.761720 11.238780
5 0.66332 4.825659 98.892702 1.107298 6.449318 2.557306 13.403070
6 0.66332 4.987632 98.017623 1.982377 6.784354 2.634121 15.647823
7 0.66332 5.167564 97.434794 2.565206 7.169577 3.013669 15.998777
8 0.66332 5.372594 96.296843 3.703157 7.497493 2.643221 16.021299
9 0.66332 5.608561 96.175970 3.824030 7.914637 2.820001 18.845901
10 0.66332 5.853195 96.164944 3.835056 8.439878 2.858989 20.281296
11 0.66332 6.130609 95.833816 4.166184 9.027518 2.947235 20.608590
12 0.66332 6.430357 95.788050 4.211950 9.667616 2.767598 24.225521
13 0.66332 6.751178 95.748409 4.251591 10.368093 2.842841 24.097100
14 0.66332 7.086480 94.935350 5.064650 10.978069 3.107622 27.070527
15 0.66332 7.464485 94.491172 5.508828 11.502707 2.789158 28.092316
16 0.66332 7.935436 93.545100 6.454900 12.093391 2.975199 30.969127
17 0.66332 8.487437 93.050368 6.949632 12.679081 3.270948 35.268139
18 0.66332 9.140589 92.943528 7.056472 13.205917 3.322945 31.097404
19 0.66332 9.872223 92.302678 7.697322 13.659439 3.269352 33.961858
20 0.66332 10.732707 92.027686 7.972314 14.133382 3.624376 34.482303
21 0.66332 11.779806 91.845762 8.154238 14.587435 3.865872 32.308234
22 0.66332 13.130834 91.278570 8.721430 14.962935 3.980520 34.330859
23 0.66332 14.830692 91.231425 8.768575 15.269924 4.532727 36.473540
24 0.66332 17.052213 90.537710 9.462290 15.381560 4.595348 39.668080
25 0.66332 20.140944 89.686207 10.313793 15.426583 4.664751 36.024495
26 0.66332 24.724328 88.878472 11.121528 15.468103 5.097378 36.976562
27 0.66332 32.492166 87.921904 12.078096 15.660284 5.389315 31.355033
28 0.66332 48.312975 86.565497 13.434503 16.097469 4.971545 32.414358
29 0.66332 100.000000 84.597976 15.402024 16.501998 4.614898 36.390796

In [18]:
# gamma = 4
#Objective: 5.653304

const_spend_pct = 0.013135043394
var_spend_pcts = pd.Series([0.035822978028495614, 0.037142880037028099, 0.038232214787301988, 0.039230415975672818, 0.040201921370750057, 0.04128306632974512, 0.042562738008764228, 0.044203187483426658, 0.04612997210464486, 0.048403252082914246, 0.050913960441316269, 0.053804872403093118, 0.056961354837341847, 0.059937273709488398, 0.062775232619198637, 0.066025078674551443, 0.070108404481856915, 0.075376643015682129, 0.081397919560765669, 0.087812983760457097, 0.094573649661557693, 0.10292659618285707, 0.11353980019663899, 0.12569730816023159, 0.14198879320676885, 0.16414995001568367, 0.19648563495632851, 0.25250659931593217, 0.3706052129208936, 1.0])
stock_allocations = pd.Series([0.90099032199651863, 0.89924731857030904, 0.89924731828348214, 0.89876830969341182, 0.89770857617754174, 0.89503554351652825, 0.89232180232291247, 0.88505129512940295, 0.8648266951006971, 0.86372639925213479, 0.85995907966929552, 0.85302430311396027, 0.85289466074670783, 0.85270233833392206, 0.85157005570710653, 0.84601302356094343, 0.83930541580936513, 0.83354308642897557, 0.82595372913293841, 0.82174793960063786, 0.80903946044432207, 0.80583306771707164, 0.79764354604204957, 0.79580225813480721, 0.75787634532529125, 0.75092813511924261, 0.73890315821297692, 0.73198221026156862, 0.72778224605194297, 0.69545344502127571])


model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 4.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.axis([ 0, 30, 0, 60,])

plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:29:15 Create TensorFlow graph and session
16:29:15 Create constants stocks_return
16:29:15 Create constants bonds_return
16:29:15 Create variables stocks_alloc
16:29:15 Create ops for soft constraint alloc_min_0 > 0
16:29:15 Create ops for soft constraint alloc_max_1 < 1
16:29:15 Create ops for soft constraint, declining stock alloc alloc_decrease
16:29:16 Create ops for bonds_alloc
16:29:16 Create variables var_spend
16:29:16 Create ops for soft constraint vspend_min_0 > 0
16:29:16 Create cohort history, years 1928 to 1986
Out[18]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 1.313504 3.582298 90.099032 9.900968 5.160262 3.624936 6.630540
1 1.313504 3.714288 89.924732 10.075268 5.359165 3.103391 8.097281
2 1.313504 3.823221 89.924732 10.075268 5.515571 3.173570 8.218728
3 1.313504 3.923042 89.876831 10.123169 5.704384 3.158267 9.933032
4 1.313504 4.020192 89.770858 10.229142 5.893790 3.215999 9.956513
5 1.313504 4.128307 89.503554 10.496446 6.117539 3.020760 11.226966
6 1.313504 4.256274 89.232180 10.767820 6.361940 3.071008 12.884046
7 1.313504 4.420319 88.505130 11.494870 6.670104 3.385726 13.142801
8 1.313504 4.612997 86.482670 13.517330 6.937730 3.080038 13.185069
9 1.313504 4.840325 86.372640 13.627360 7.282293 3.183739 15.414317
10 1.313504 5.091396 85.995908 14.004092 7.730806 3.167866 16.522639
11 1.313504 5.380487 85.302430 14.697570 8.235829 3.238491 16.958931
12 1.313504 5.696135 85.289466 14.710534 8.790754 3.052960 19.950505
13 1.313504 5.993727 85.270234 14.729766 9.338580 3.096768 20.912103
14 1.313504 6.277523 85.157006 14.842994 9.789087 3.295179 23.194224
15 1.313504 6.602508 84.601302 15.398698 10.187303 3.009584 23.653342
16 1.313504 7.010840 83.930542 16.069458 10.641891 3.156542 27.415610
17 1.313504 7.537664 83.354309 16.645691 11.161797 3.373851 30.782408
18 1.313504 8.139792 82.595373 17.404627 11.613556 3.332475 28.301408
19 1.313504 8.781298 82.174794 17.825206 11.967621 3.476046 29.244046
20 1.313504 9.457365 80.903946 19.096054 12.243826 3.692976 28.891891
21 1.313504 10.292660 80.583307 19.416693 12.534684 3.750771 26.515161
22 1.313504 11.353980 79.764355 20.235645 12.771957 3.763388 26.411370
23 1.313504 12.569731 79.580226 20.419774 12.869307 3.993518 28.844498
24 1.313504 14.198879 75.787635 24.212365 12.895273 3.747646 31.328059
25 1.313504 16.414995 75.092814 24.907186 12.908486 3.920047 29.091372
26 1.313504 19.648563 73.890316 26.109684 12.989382 3.894303 28.021177
27 1.313504 25.250660 73.198221 26.801779 13.446328 3.841805 26.212153
28 1.313504 37.060521 72.778225 27.221775 14.707136 3.511419 29.290128
29 1.313504 100.000000 69.545345 30.454655 22.213173 3.156638 48.777261

In [19]:
# gamma = 6

#Objective: 5.056278

const_spend_pct = 0.015637933571
var_spend_pcts = pd.Series([0.031606123286676821, 0.033096808496424615, 0.034176874786077449, 0.034949059621799367, 0.035890789552150887, 0.037116351576046536, 0.038256812396706585, 0.039784074922038948, 0.041850514856008901, 0.044113104719702513, 0.046479531291527841, 0.049188553520785189, 0.052237088536900282, 0.054826487679491116, 0.05719659838516862, 0.060236417433033075, 0.06372908174792162, 0.068170731957061817, 0.073291375378079365, 0.078832696364871821, 0.084046153747339689, 0.091046547868110006, 0.10034482939483094, 0.11066289255761518, 0.1255268715113966, 0.14590561544900946, 0.17559280570819452, 0.22684587500024678, 0.33478357325498115, 1.0])
stock_allocations = pd.Series([0.88033446277449279, 0.87918369248956052, 0.87699906807696382, 0.87541489848777843, 0.87254543385561067, 0.87112127695141239, 0.8696038278716699, 0.85852400031522713, 0.84145611984667457, 0.84145611982772339, 0.83898260822198845, 0.83191026515233402, 0.83112108486165093, 0.83100258791320347, 0.83025217396697193, 0.81919450523261073, 0.81311531977557161, 0.809624259366025, 0.80101985702945411, 0.79700856489108074, 0.78660821862681019, 0.7827212663418508, 0.77416817507231639, 0.76897655138636212, 0.72297797819724985, 0.7165312149934886, 0.70342618929327727, 0.69734141791783077, 0.69379260335071868, 0.65522743833716124])

model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 6.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.axis([ 0, 30, 0, 60,])

plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:29:49 Create TensorFlow graph and session
16:29:49 Create constants stocks_return
16:29:49 Create constants bonds_return
16:29:49 Create variables stocks_alloc
16:29:50 Create ops for soft constraint alloc_min_0 > 0
16:29:50 Create ops for soft constraint alloc_max_1 < 1
16:29:50 Create ops for soft constraint, declining stock alloc alloc_decrease
16:29:50 Create ops for bonds_alloc
16:29:50 Create variables var_spend
16:29:50 Create ops for soft constraint vspend_min_0 > 0
16:29:50 Create cohort history, years 1928 to 1986
Out[19]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 1.563793 3.160612 88.033446 11.966554 4.953237 3.621346 6.222775
1 1.563793 3.309681 87.918369 12.081631 5.167058 3.182037 7.550773
2 1.563793 3.417687 87.699907 12.300093 5.321259 3.269539 7.681328
3 1.563793 3.494906 87.541490 12.458510 5.478367 3.255841 9.172201
4 1.563793 3.589079 87.254543 12.745457 5.658951 3.296806 9.249442
5 1.563793 3.711635 87.112128 12.887872 5.892822 3.132238 10.357694
6 1.563793 3.825681 86.960383 13.039617 6.115328 3.177317 11.837691
7 1.563793 3.978407 85.852400 14.147600 6.402631 3.470715 12.091717
8 1.563793 4.185051 84.145612 15.854388 6.690486 3.204995 12.242438
9 1.563793 4.411310 84.145612 15.854388 7.033949 3.288447 14.320895
10 1.563793 4.647953 83.898261 16.101739 7.460332 3.267263 15.356741
11 1.563793 4.918855 83.191027 16.808973 7.939837 3.332161 15.824846
12 1.563793 5.223709 83.112108 16.887892 8.481095 3.156292 18.639529
13 1.563793 5.482649 83.100259 16.899741 8.978700 3.189385 19.849604
14 1.563793 5.719660 83.025217 16.974783 9.380353 3.360150 21.883548
15 1.563793 6.023642 81.919451 18.080549 9.779419 3.109486 22.478892
16 1.563793 6.372908 81.311532 18.688468 10.193235 3.237130 26.113748
17 1.563793 6.817073 80.962426 19.037574 10.661888 3.407857 29.143336
18 1.563793 7.329138 80.101986 19.898014 11.081542 3.360888 27.069028
19 1.563793 7.883270 79.700856 20.299144 11.434131 3.518694 27.624942
20 1.563793 8.404615 78.660822 21.339178 11.653879 3.717353 27.133932
21 1.563793 9.104655 78.272127 21.727873 11.962620 3.681786 25.068861
22 1.563793 10.034483 77.416818 22.583182 12.283094 3.729875 24.834411
23 1.563793 11.066289 76.897655 23.102345 12.459354 3.901329 27.646007
24 1.563793 12.552687 72.297798 27.702202 12.683537 3.682808 30.325249
25 1.563793 14.590562 71.653121 28.346879 12.940746 3.877009 29.080604
26 1.563793 17.559281 70.342619 29.657381 13.304662 3.831778 27.998996
27 1.563793 22.684588 69.734142 30.265858 14.123992 3.799812 28.257736
28 1.563793 33.478357 69.379260 30.620740 15.977814 3.506184 32.839585
29 1.563793 100.000000 65.522744 34.477256 28.039147 2.982893 63.160835

In [20]:
#gamma = 8

#Objective: 4.344292

const_spend_pct = 0.018033736915
var_spend_pcts = pd.Series([0.027420789200650691, 0.028952687799914788, 0.029987640366612587, 0.030563608815231274, 0.0314680099541\
00477, 0.0327892165458714, 0.033812194391987072, 0.035195378137199251, 0.037327252304951188, 0.039532656754857148, 0.04171133948016\
5768, 0.044205216561892961, 0.047114695564658048, 0.049375968207551808, 0.051345728157015109, 0.05421706340461499, 0.05717491967092\
8206, 0.060803826336315064, 0.065058176584335367, 0.069765865219855477, 0.073504852805054741, 0.079201426822341936, 0.0872191175663\
84868, 0.095745358220279422, 0.10919882479478223, 0.12781818289102989, 0.15488310518171028, 0.20138642112957703, 0.2991516975921870\
2, 1.0])
stock_allocations = pd.Series([0.85995594522555152, 0.85919296367537157, 0.85473424640534235, 0.85205353144872398, 0.84754303090697\
891, 0.84731613048078414, 0.84674123511449917, 0.83207329884146453, 0.81886882987838872, 0.81886882879077094, 0.81787428397721995, \
0.81082294011238876, 0.80935846467721007, 0.80915247713021987, 0.80869082405588244, 0.79235039525216477, 0.78697518924653442, 0.785\
66021347854842, 0.77606352982164872, 0.77226108408967586, 0.76418701504255893, 0.75957408240400748, 0.7506635662237392, 0.742109947\
15180123, 0.68810212447786767, 0.68212008875931995, 0.66794616058167322, 0.66270963157142881, 0.65979617905489518, 0.61499306500466\
799])



model = SafeWithdrawalModel(returns_list = [real_stocks, real_bonds],
                            names_list = ["stocks","bonds"],
                            allocations_list = [stock_allocations, bond_allocations],
                            start_val = 100.0,
                            const_spend_pct = const_spend_pct,
                            var_spend_pcts = var_spend_pcts,
                            gamma = 8.0,
                            survival=None
)

# generate cohorts
model.cohort_history = CohortHistory(model)

df_summary = model.cohort_history.summarize_by_year()
df_years = model.cohort_history.spend_by_year()

from pylab import rcParams
rcParams['figure.figsize'] = 8, 6

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    plt.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label='_nolegend_')


plt.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
plt.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
plt.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
plt.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=1, label='+1 SD')
plt.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=1, label='-1 SD')
plt.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
plt.legend(loc="upper left", bbox_to_anchor=[0, 1])
plt.axis([ 0, 30, 0, 60,])

plt.show()

df_summary['const_spend'] = const_spend_pct * 100
df_summary['var_spend'] = var_spend_pcts * 100
df_summary['stocks'] = stock_allocations * 100
df_summary['bonds'] = 100 - df_summary['stocks']
cols = ['const_spend', 'var_spend', 'stocks', 'bonds', 'spend_mean', 'spend_min', 'spend_max']
df_summary[cols]


16:30:30 Create TensorFlow graph and session
16:30:30 Create constants stocks_return
16:30:30 Create constants bonds_return
16:30:30 Create variables stocks_alloc
16:30:31 Create ops for soft constraint alloc_min_0 > 0
16:30:31 Create ops for soft constraint alloc_max_1 < 1
16:30:31 Create ops for soft constraint, declining stock alloc alloc_decrease
16:30:31 Create ops for bonds_alloc
16:30:31 Create variables var_spend
16:30:31 Create ops for soft constraint vspend_min_0 > 0
16:30:31 Create cohort history, years 1928 to 1986
Out[20]:
const_spend var_spend stocks bonds spend_mean spend_min spend_max
0 1.803374 2.742079 85.995595 14.004405 4.740133 3.604047 5.817874
1 1.803374 2.895269 85.919296 14.080704 4.954235 3.239628 6.990869
2 1.803374 2.998764 85.473425 14.526575 5.102155 3.338954 7.131929
3 1.803374 3.056361 85.205353 14.794647 5.231042 3.328030 8.398137
4 1.803374 3.146801 84.754303 15.245697 5.401859 3.355214 8.518801
5 1.803374 3.278922 84.731613 15.268387 5.640160 3.220717 9.473607
6 1.803374 3.381219 84.674124 15.325876 5.843391 3.262251 10.780489
7 1.803374 3.519538 83.207330 16.792670 6.106355 3.482782 11.021653
8 1.803374 3.732725 81.886883 18.113117 6.406250 3.268162 11.260031
9 1.803374 3.953266 81.886883 18.113117 6.743933 3.372960 13.164527
10 1.803374 4.171134 81.787428 18.212572 7.143623 3.348702 14.116656
11 1.803374 4.420522 81.082294 18.917706 7.594771 3.409277 14.611927
12 1.803374 4.711470 80.935846 19.064154 8.120312 3.246629 17.223627
13 1.803374 4.937597 80.915248 19.084752 8.577191 3.273653 18.662278
14 1.803374 5.134573 80.869082 19.130918 8.939804 3.422242 20.469916
15 1.803374 5.421706 79.235040 20.764960 9.347002 3.205772 21.213132
16 1.803374 5.717492 78.697519 21.302481 9.727362 3.317734 24.702617
17 1.803374 6.080383 78.566021 21.433979 10.143674 3.450519 27.386865
18 1.803374 6.505818 77.606353 22.393647 10.531063 3.401243 25.689757
19 1.803374 6.976587 77.226108 22.773892 10.881431 3.542791 25.921611
20 1.803374 7.350485 76.418702 23.581298 11.038810 3.749375 25.284224
21 1.803374 7.920143 75.957408 24.042592 11.355122 3.645121 23.647565
22 1.803374 8.721912 75.066357 24.933643 11.747497 3.685455 24.225980
23 1.803374 9.574536 74.210995 25.789005 11.987134 3.848685 26.258617
24 1.803374 10.919882 68.810212 31.189788 12.405500 3.660422 29.108868
25 1.803374 12.781818 68.212009 31.787991 12.913144 3.878824 28.843689
26 1.803374 15.488311 66.794616 33.205384 13.578350 3.829700 27.946354
27 1.803374 20.138642 66.270963 33.729037 14.794993 3.828073 30.237365
28 1.803374 29.915170 65.979618 34.020382 17.323086 3.581394 37.396578
29 1.803374 100.000000 61.499307 38.500693 35.762553 3.056128 85.260856

In [21]:
# make a plotly chart

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)

from matplotlib import colors
mycolors =  list(colors.cnames)
random.shuffle(mycolors)

for ix in range(model.first_year, model.first_year + model.ret_cohorts):
    ax.plot(df_years[str(ix)], linewidth=1, color=mycolors[ix-model.first_year], alpha=0.5, label=str(ix))


ax.plot(df_summary.spend_mean, color='k', linewidth=2, label='Mean')
ax.plot(df_summary.spend_max, color='g', linewidth=2, label='Best')
ax.plot(df_summary.spend_min, color='r', linewidth=2, label='Worst')
ax.plot(df_summary.spend_mean + df_summary.spend_sd, color='b', linewidth=2, label='+1 SD')
ax.plot(df_summary.spend_mean - df_summary.spend_sd, color='b', linewidth=2, label='-1 SD')
ax.fill_between(df_summary.index, df_summary.spend_mean + df_summary.spend_sd, 
                 df_summary.spend_mean - df_summary.spend_sd, color='blue', alpha='0.1')

plt.ylabel("Annual spending as % of starting portfolio")
plt.xlabel("Retirement year")
#plotly pukes on legend
#ax.legend(loc="upper left", bbox_to_anchor=[0, 1])
ax.axis([ 0, 30, 0, 60,])

plt.show()



In [22]:
plot_url = py.iplot_mpl(fig, filename='chart')
plot_url


/home/ubuntu/anaconda2/envs/tensorflow/lib/python2.7/site-packages/plotly/matplotlylib/renderer.py:445: UserWarning:

Dang! That path collection is out of this world. I totally don't know what to do with it yet! Plotly can only import path collections linked to 'data' coordinates

Out[22]:

In [ ]: