prodigalson8523
New member
I'll start by saying I'm not a Python expert - just a dabbler so there is probably a far better way to do this...
Background
I'm 100% in equities at the moment, but approaching that final 1/3 of my investing life. I've started to think about my 'transition to bonds' plan and what this might mean for my retirement FIRE number so I built a python app that uses Monte Carlo simulation to model my portfolio.
I've taken historic performance data for the S&P500, Treasury & Corporate Bonds, and inflation since 1929 and found the distribution of this data in terms of average and standard deviation. I''ve then created weight growths for 100% equities, 90/10, 80/20, 70/30 etc. I've used all of this data to run 500,000 sims of my lifetime for each of these allocations. I get some simple stats and a distribution curve of my balance for each of the simulations. I understand this US centric view might not give a 100% accurate picture, but given the US is ~70% of the global market it's a fair guess to say that over time what ever the US market does the world is likely to follow.
Anyway I wanted to share the code for review, feedback, or just to use yourself in your own models.
Background
I'm 100% in equities at the moment, but approaching that final 1/3 of my investing life. I've started to think about my 'transition to bonds' plan and what this might mean for my retirement FIRE number so I built a python app that uses Monte Carlo simulation to model my portfolio.
I've taken historic performance data for the S&P500, Treasury & Corporate Bonds, and inflation since 1929 and found the distribution of this data in terms of average and standard deviation. I''ve then created weight growths for 100% equities, 90/10, 80/20, 70/30 etc. I've used all of this data to run 500,000 sims of my lifetime for each of these allocations. I get some simple stats and a distribution curve of my balance for each of the simulations. I understand this US centric view might not give a 100% accurate picture, but given the US is ~70% of the global market it's a fair guess to say that over time what ever the US market does the world is likely to follow.
Anyway I wanted to share the code for review, feedback, or just to use yourself in your own models.
Code:
import random
import matplotlib.pyplot as plt
import numpy as np
import statistics as st
from scipy.stats import norm
def monte_carlo_simulation(starting_balance, average_growth, std_dev, years, periodic_payment):
# list to store the results of each iteration
results = []
return_per = []
# number of iterations to run
iterations = 500000
for i in range(iterations):
# current balance set to starting balance
balance = starting_balance
for j in range(years):
# add periodic payment to balance
balance += periodic_payment
# increase balance by growth prediction percentage
growth = random.normalvariate(average_growth, std_dev)
return_per.append(growth)
balance *= (1 + growth)
# add final balance to results list
results.append(balance)
# calculate min, max, and average of results
min_result = min(results)
max_result = max(results)
avg_result = sum(results)/len(results)
twenty = np.percentile(results,80)
eighty = np.percentile(results,20)
stdevAmount = st.stdev(results)
lessThanStart = norm.cdf(starting_balance,avg_result,stdevAmount)*100
# plot results as normal distribution curve
# plt.hist(results, bins=100, color='g')
plt.hist(results,bins=100,range=(100000,5000000))
plt.show()
print(f"Min :${min_result:,.0f}")
print(f"Avg :${avg_result:,.0f}")
print(f"Max :${max_result:,.0f}")
print(f"20% chance of at least :${twenty:,.0f}")
print(f"80% chance of at least :${eighty:,.0f}")
print(f"68% (1sd) chance of between : ${avg_result-stdevAmount:,.0f} to ${avg_result+stdevAmount:,.0f}")
print(f"There is a {lessThanStart:,.2f}% probability of having less than ${starting_balance}")
# print(f"Most likely Range 2sd : ${avg_result-(2*stdevAmount):,.0f} to ${avg_result+(2*stdevAmount):,.0f}")
# print(f"20% return :{np.percentile(return_per,.2):,.4f}")
# print(f"95% return :{np.percentile(return_per,.95):,.4f}")
# print(f"Best Return :{max(return_per):,.4f}")
# print(f"Worst Return :{min(return_per):,.4f}")
# print(f"Av Return :{sum(return_per)/len(return_per):,.4f}")
# print(f"Std Dev Return :{st.stdev(return_per):,.4f}")
# put your data here...
myBalance = 25000
myTimescale = 30
myPeriodicPayments = 8000
# 100% Equities - Growths and StdDev are historic and adjusted for historic inflation
print("100% Equities")
monte_carlo_simulation(myBalance, 0.086, 0.1923, myTimescale, myPeriodicPayments)
# 80%/20% Equities/Bonds
print("-----------------------------------------------")
print("80% Equities")
monte_carlo_simulation(myBalance, 0.0733, 0.1689, myTimescale, myPeriodicPayments)
# 70%/30%
print("-----------------------------------------------")
print("70% Equities")
monte_carlo_simulation(myBalance, 0.0669, 0.1572, myTimescale, myPeriodicPayments)
# 60%/40%
print("-----------------------------------------------")
print("60% Equities")
monte_carlo_simulation(myBalance, 0.0605, 0.1455, myTimescale, myPeriodicPayments)
# 50%/50%
print("-----------------------------------------------")
print("50% Equities")
monte_carlo_simulation(myBalance, 0.0541, 0.1338, myTimescale, myPeriodicPayments)
#100% Bonds
print("-----------------------------------------------")
print("100% Bonds")
monte_carlo_simulation(myBalance, 0.0222, 0.0753, myTimescale, myPeriodicPayments)
# For Comparison - non inflation adjusted 100% equities
print("-----------------------------------------------")
print("100% Equities - not inflation adjusted")
monte_carlo_simulation(myBalance, 0.1182, 0.1936, myTimescale, myPeriodicPayments)