Retirement Simulator #10: Test and Compare Different Retirement dates (And make it pretty!)
Published Dec. 31, 2021, 9:39 p.m.
Hello there! Welcome to the 10th and final tutorial in the 'Retirement Simulator' series! Congrats on making it this far!
Python is a beautiful thing. We should make better looking plots then, don't you think? The bulk of this tutorial will be spent on fixing some plot formatting features, however, the real hero of this tutorial is the ease with which we can adjust our code so that we can pick 3 (or more) retirement dates and have python display the results for our enjoyment and fiscal edification.
First let us manually reset our figure size. Go to the code where we first initialized our 'fig' and 'axs' variables and add the following:
fig, axs = plt.subplots(2,2,figsize = (9, 4.5))At the top of our code we now need to import mdates:
import matplotlib.dates as mdatesNow scroll down to where we make the plot (in my code it is line 137) where our code has "if sim%10 == 0:" Inside the if statement we need to add the following:
ax = axs[0,1]
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
for label in ax.get_ticklabels(which="major"):
label.set(rotation=30, horizontalalignment = "right")
plt.subplots_adjust(wspace = 0.32, hspace = 0.42)
# axs[0,1].plot(dates[sim], money[sim])
ax.plot(dates[sim], money[sim])
plt.pause(0.0001)If all of that ran without errors then your plot should have looked better (format-wise). Now that we have the format looking good enough we can do some programming!
At this point we have several lines of code that make up the active portion of our program. It starts somewhere in the middle just before the for loop that uses the nsims variable and it ends nearly at the end of the file. What we would like to do is choose some retirement dates and see how an extra 5 or 10 years of employment affects our retirement account. To do this most simply we need to define some retirement dates in a list. I propose:
retirement_dates = [
date(2040, 3, 31),
date(2035, 3, 31),
date(2030, 3, 31)
]Also, we might like to color-code our results so we will make the following list as well:
colors = [
'red',
'green',
'blue',
]Now we will zip these two lists together; this will associate 'red' with the first date in our list, 'green' with the second date... etc.
Now for one last gasp at cleverness. We will start a new for loop as follows:
for col, retirement_date in zip(colors, retirement_dates):
...
...Then if you are using Notepad++ or a similar text editor, you can highlight all of the text that makes up our code, hit 'tab' and now the code is inside our new for loop! All that is left to do is to take advantage of the col (color) variable that we created.
Now inside our code, where we make a plot we can add:
ax.plot(... ... , color = col)
For the plots at the end where we make the thick black dashed line we can do the following:
axs[0,1].plot(dates[0], average_amount, linewidth=3, linestyle ='--', color = col, label=f"Retire {todays_date.year}")
axs[0,1].plot(dates[0], average_amount, linewidth=3, linestyle ='--', color = 'k')What we have done is plot the average _amount twice: first colored with our scheme, and with a label, then again in black without a label. Python will plot them both but the black colored plot will be on top while the color scheme will show up in our legend with the color that we want. Yay for trickery!
We can add the "color=col" statement for the histogram of final amounts also.
The final change we need to make is to establish a binning convention for our histograms. Python chooses bins that are very different for each simulation and that reduces the clarity of the histogram. The fix is simple enough. Before the histogram is plotted we need to add the following lines:
my_bins = np.linspace(5e4, 8e5, 200)Then inside our call to plot the histogram we must add:
axs[1,0].hist(finalamounts, bins=my_bins, ...Further we can add:
axs[1,0].hist(finalamounts, bins=my_bins, range = (5e4,8e5), ...And now our simulation gives us some plots that we can look at and analyze! Now for the rest of the video, I go on to explain how you can change the simulation to targe a retirement amount rather than a retirement date. (It's actually not that difficult).
In case you made it this far but can't get your code to work, here is a complete code example:
import sys, os
from datetime import date, timedelta
import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy.optimize import curve_fit
todays_date = date.today()
two_weeks = timedelta(days = 14)
one_year = timedelta(days = 365.25)
one_month = one_year/12
semi_monthly = one_month/2
pay_frequency = two_weeks
amount_already_in_retirement = 10000
paychecks_per_year = int(round(one_year/pay_frequency, 0))
Current_Salary = 50000
one_paycheck = Current_Salary/paychecks_per_year
add_to_retirement = 0.10*one_paycheck
# add_to_retirement = 100
retirement_date = date(2030, 3, 31)
retirement = amount_already_in_retirement
annual_interst_rate = 0.10
"""
want to convert annual_interst_rate into a bi-weekly interest rate equivalent.
bwir should be such that (1+bwir)**26 = air (annual_interst_rate)
26*log((1+bwir)) = log(1+air)
log((1+bwir)) = (1/26)*log(1+air)
1 + bwir = exp((1/26)*log(1+air))
bwir = exp((1/26)*log(1+air)) - 1
"""
def get_paycheck_interest_rate(annual_interest_rate, paychecks_per_year):
return np.exp((1/paychecks_per_year)*np.log(1+annual_interest_rate))-1
paycheck_interest_rate = get_paycheck_interest_rate(annual_interst_rate, paychecks_per_year)
mean = 0.07
sigma = 0.07
annual_interst_rate = random.gauss(mean, sigma)
paycheck_interest_rate = get_paycheck_interest_rate(annual_interst_rate, paychecks_per_year)
filename = "dow_market_returns.txt"
f = open(filename, 'r')
lines = f.readlines()
f.close()
years=[]
rates = []
for line in lines:
year, rate = line.split(',')
years.append(int(year))
rates.append(float(rate)/100)
annual_interst_rate = random.choice(rates)
year_check = todays_date.year
fig, axs = plt.subplots(2,2, figsize=(9, 4.5)) # axs --> axs[0,0] 1st plot, axs[0,1], 2nd plot ...
rate_bins = np.arange(-.55, .65, 0.1)
bin_centers = np.arange(-.50, .65, 0.1)
the_weights, thebins, zzz = axs[0,0].hist(rates, bins=rate_bins)
def gaus(x, a0, x0, sigma):
return a0*np.exp(-(x-x0)**2 / (2*sigma**2))
popt, pcov = curve_fit(gaus, bin_centers, the_weights)
print(popt)
print(*popt)
x_vals = np.arange(-.55, .75, 0.01)
y_vals = gaus(x_vals, *popt)
axs[0,0].plot(x_vals, y_vals)
def makesampledistribution(xvals, yvals):
ymax = max(yvals)
norm = 100/ymax
yvals = [norm*y for y in yvals]
x_dist = []
for i,y in enumerate(yvals):
inty = int(y)
if inty > 0:
for j in range(inty):
x_dist.append(xvals[i])
return x_dist
x_dist = makesampledistribution(x_vals, y_vals)
axs[0,0].hist(x_dist, bins =x_vals)
retirement_date = date(2030, 3, 31)
retirement_dates = [
date(2045, 3, 31),
date(2040, 3, 31),
date(2035, 3, 31),
]
retirement_amounts = [
100000,
400000,
800000,
]
colors = ['red', 'green', 'blue']
## 0.09+0.045+0.03
for col, retirement_date in zip(colors, retirement_dates):
i=0
nsims = 10000
money = [[] for x in range(nsims)]
dates = [[] for x in range(nsims)]
for sim in range(nsims):
todays_date = date.today()
# retirement_date = date(2030, 3, 31)
retirement = amount_already_in_retirement
Current_Salary = 50000
one_paycheck = Current_Salary/paychecks_per_year
add_to_retirement = 0.10*one_paycheck
while todays_date < retirement_date:
i+=1
retirement += add_to_retirement
retirement *= (1+paycheck_interest_rate)
todays_date += pay_frequency
this_year = todays_date.year
if this_year > year_check:
# annual_interst_rate = random.gauss(mean, sigma)
annual_interst_rate = random.choice(x_dist)
paycheck_interest_rate = get_paycheck_interest_rate(annual_interst_rate, paychecks_per_year)
# print(this_year, annual_interst_rate)
year_check = todays_date.year
money[sim].append(retirement)
dates[sim].append(todays_date)
if sim%100 == 0:
ax = axs[0,1]
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
for label in ax.get_xticklabels(which="major"):
label.set(rotation=30, horizontalalignment = 'right')
plt.subplots_adjust(wspace = 0.32, hspace = 0.42)
ax.plot(dates[sim], money[sim], color=col)
# print(f"I have {retirement:.2f} as of {todays_date}.")
finalamounts = []
average_amount = np.zeros(shape = len(dates[0]))
for sim in range(nsims):
finalamounts.append(money[sim][-1])
money[sim] = np.array(money[sim])
average_amount += money[sim]/nsims
axs[0,1].plot(dates[0], average_amount, linewidth = 3, linestyle = "--", color = col, label = f"Retire {todays_date.year}")
axs[0,1].plot(dates[0], average_amount, linewidth = 3, linestyle = "--", color = "k")
axs[0,1].legend()
my_range = (4e4, 1e6)
my_bins = np.linspace(*my_range, 100)
ax = axs[1,0]
# for label in ax.get_xticklabels(which="major"):
# label.set(rotation=30, horizontalalignment = 'right')
ax.hist(finalamounts, bins = my_bins, range = my_range, histtype="step", color = col, label = f"Retire {todays_date.year}")
ax.legend()
plt.pause(0.0001)
plt.show() Keep watching and let me know how it went! See you in the next series!
-
Retirement Simulator Introduction and Overview
-
Retirement Simulator #1: March through time with datetime and put money under the mattress!
-
Retirement Simulator #2: Apply an interest rate and verify it!
-
Retirement Simulator #3: Use a Gaussian distribution of rates!
-
Retirement Simulator #3+: Verify our Gaussian distribution!
-
Retirement Simulator #4: Make a US market rates distribution
-
Retirement Simulator #5: Fit the Market Data with scipy curve_fit
-
Retirement Simulator #6: Sampling from our custom market distribution
-
Retirement Simulator #7: Organize multiple plots with subplots! and learn about plt.pause()
-
Retirement Simulator #8: Multiverse Investing: Simulate 1000 random investors
-
Retirement Simulator #9: Track averages and final outcomes
-
Retirement Simulator #10: Test and Compare Different Retirement dates (And make it pretty!)
(currently viewing)