Retirement simulation

bond_returns.csv

Retirement.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "\n",
    "from IPython.display import display\n",
    "import ipywidgets as widgets\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.lines as mlines\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "stocks = pd.read_csv('./stock_returns.csv', header=None, names=['year', 'yield'], index_col=0, squeeze=True)\n",
    "bonds = pd.read_csv('./bond_returns.csv', header=None, names=['year', 'yield'], index_col=0, squeeze=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "starting_year = datetime.date.today().year\n",
    "\n",
    "social_security = 2000\n",
    "full_soc_sec = 67\n",
    "\n",
    "inflation = 1.0322"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.19 s, sys: 175 ms, total: 2.37 s\n",
      "Wall time: 2.38 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Generate one million simulation chains of the next 100 years of stock and bond yields.\n",
    "sim_stocks = pd.DataFrame(\n",
    "    np.random.choice(stocks, size=(1000000, 100), replace=True),\n",
    "    columns=np.arange(starting_year + 1, starting_year + 101)\n",
    ")\n",
    "\n",
    "sim_bonds = pd.DataFrame(\n",
    "    np.random.choice(bonds, size=(1000000, 100), replace=True),\n",
    "    columns=np.arange(starting_year + 1, starting_year + 101)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bond_function(allocation, birth_year, year):\n",
    "    if allocation == 'Age as bond percentage':\n",
    "        return min((year - birth_year) / 100, 1.)\n",
    "    return 1. - int(allocation[:2]) / 100\n",
    "\n",
    "def run_simulation(birth_year, retirement_age, starting_value, contribution_rate, burn_rate, allocation):\n",
    "    bond_allocation = pd.Series({\n",
    "        year: bond_function(allocation, birth_year, year) for year in sim_stocks.columns\n",
    "    })\n",
    "    contributions = pd.Series({\n",
    "        year: contribution_rate * (year - birth_year < retirement_age) for year in sim_stocks.columns\n",
    "    })\n",
    "    spending = pd.Series({\n",
    "        year: (burn_rate * (year - birth_year >= retirement_age) -\n",
    "               12 * social_security * (year - birth_year >= full_soc_sec))\n",
    "        for year in sim_stocks.columns\n",
    "    })\n",
    "\n",
    "    results = pd.DataFrame({starting_year: pd.Series(starting_value for _ in range(1000000))})\n",
    "    for year in sim_stocks.columns:\n",
    "        yield_ = 1. + (1 - bond_allocation[year]) * sim_stocks[year] + bond_allocation[year] * sim_bonds[year]\n",
    "        results[year] = yield_ * (results[year - 1] + contributions[year] - spending[year]) / inflation\n",
    "\n",
    "    fig, axs = plt.subplots(2, 1, sharex=True, figsize=(12, 16))\n",
    "\n",
    "    ax = axs[0]\n",
    "    df = pd.DataFrame({\n",
    "        year: results[year].describe([.03, .5, .97])[['97%', '50%', '3%']]\n",
    "        for year in sim_stocks.columns\n",
    "    }).T\n",
    "    df.columns = ['97th percentile', '50th percentile', '3rd percentile']\n",
    "    df.plot(ax=ax, logy=True)\n",
    "    ax.set_title('Value of Retirement Fund')\n",
    "    ax.set_xlabel('Year')\n",
    "    ax.set_ylabel('Value (log USD)')\n",
    "\n",
    "    ax = axs[1]\n",
    "    s = results.lt(0).sum() / 1e6\n",
    "    s.plot(subplots=True, ax=ax)\n",
    "    ax.axhline(0.03, color='red', linestyle='-.')\n",
    "    ax.legend(handles=[\n",
    "        mlines.Line2D([], [], label='Chance of running out of money (%)'),\n",
    "        mlines.Line2D([], [], color='red', linestyle='-.', label='3% threshold')\n",
    "    ])\n",
    "    ax.set_title('Risk of Ruin')\n",
    "    ax.set_xlabel('Year')\n",
    "    ax.set_ylabel('Probability of negative value')\n",
    "    \n",
    "    for ax in axs:\n",
    "        ax.axvline(birth_year + full_soc_sec, color='blue', linestyle=':')\n",
    "        ax.axvspan(birth_year + retirement_age, birth_year + 100, facecolor='green', alpha=0.15)\n",
    "\n",
    "    handles = [\n",
    "        mpatches.Patch(color='green', alpha=0.15, label='Retirement up to age 100'),\n",
    "        mlines.Line2D([], [], color='blue', linestyle=':', label='Full Social Security benefits')\n",
    "    ]\n",
    "    fig.legend(handles=handles, bbox_to_anchor=(0.95, 0.92), borderaxespad=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "47906b1da6584db8a9231bd2009d9fcc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedIntText(value=1977, description='Birth Year:', layout=Layout(width='50%'), max=20…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.run_simulation(birth_year, retirement_age, starting_value, contribution_rate, burn_rate, allocation)>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "layout = widgets.Layout(width='50%')\n",
    "style = {'description_width': '150px'}\n",
    "\n",
    "birth_widget = widgets.BoundedIntText(\n",
    "    value=1977,\n",
    "    min=starting_year - 100,\n",
    "    max=starting_year,\n",
    "    step=1,\n",
    "    description='Birth Year:',\n",
    "    disabled=False,\n",
    "    layout=layout,\n",
    "    style=style\n",
    ")\n",
    "retirement_widget = widgets.BoundedIntText(\n",
    "    value=67,\n",
    "    min=40,\n",
    "    max=80,\n",
    "    step=1,\n",
    "    description='Retirement Age:',\n",
    "    disabled=False,\n",
    "    layout=layout,\n",
    "    style=style\n",
    ")\n",
    "starting_widget = widgets.IntText(\n",
    "    value=100000,\n",
    "    step=100,\n",
    "    description='Starting Value:', disabled=False, layout=layout, style=style\n",
    ")\n",
    "contribution_widget = widgets.IntText(\n",
    "    value=25000,\n",
    "    step=100,\n",
    "    description='Yearly Contribution:', disabled=False, layout=layout, style=style\n",
    ")\n",
    "burn_widget = widgets.IntText(\n",
    "    value=65000,\n",
    "    step=100,\n",
    "    description='Retirement Burn Rate:', disabled=False, layout=layout, style=style\n",
    ")\n",
    "allocation_widget = widgets.Dropdown(\n",
    "    options=[\n",
    "        'Age as bond percentage',\n",
    "        '80 / 20 stocks vs bonds',\n",
    "        '60 / 40 stocks vs bonds'\n",
    "    ],\n",
    "    value='Age as bond percentage',\n",
    "    description='Allocation Strategy:',\n",
    "    disabled=False,\n",
    "    layout=layout,\n",
    "    style=style\n",
    ")\n",
    "\n",
    "widgets.interact_manual(\n",
    "    run_simulation,\n",
    "    birth_year=birth_widget, retirement_age=retirement_widget, starting_value=starting_widget,\n",
    "    contribution_rate=contribution_widget, burn_rate=burn_widget, allocation=allocation_widget,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

stock_returns.csv