Monte Carlo Simulation of net worth in Python..

[prodigalson8523]

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.

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)

 

[rhino521]

@prodigalson8523 Yeah don't use a normal distribution. It's the wrong one to use.

The accepted theoretical way to do this is by using Geometric Brownian Motion

 

[nayakara]

@rhino521 Geometric brown Ian motion? Is that where the market goes down the toilet?

 

[trishn40]

@prodigalson8523 Have you considered incorporating drift? Markets over the long term 'drift' upwards, your model doesn't account for this - as such, your model will tend to understate realised outcomes.

Consider looking into geometric brownian motion as a way to model markets, it's not right (indeed, nothing is really 'right') but it's perhaps more accurate than your current approach.

 

[der]

@trishn40 This model very obviously incorporates drift. That's what the mean of the normal distribution represents.

 

[prodigalson8523]

@trishn40 I'll explore brownian motion. Hopefully there is a Python library for that!

 

[prodigalson8523]

I've found a formula for Brownian motion... I've found something that suggests using random.gauss as an approximation to Brownian motion. I'll have a look at that.

 

[trishn40]

@prodigalson8523 That works as a decent approximation. The application of geometric brownian motion is actually pretty simple (or at least when i've applied it) - you're just adding some constant (which is determined by your return series) to your model outcomes (which are determined by the standard deviation of returns) - the constant ensures that over time, the return series 'drifts' upwards.

Truthfully, if you've incorporated some form of drift (in whatever way you see fit) - i'd consider that accurate enough. Outcomes will never realistically be reliable and i'd caution against point estimates - having some sort of range of acceptable outcomes/outputs would be better (which is what would be generated from your 500,000 trials - i'm certain they plot into a nice and normally distributed cone (as per your design)). I work in finance and have spent some time thinking/tinkering with these things - i have come to the conclusion that trying to forecast markets is a waste of time - you're not actually modelling the market, you're modelling how humans as a collective think and as everyone knows, humans are anything but consistent.

 

[ecsid]

@prodigalson8523 How did you check if the distribution of per unit time growth over the sampled time period was normal? Can you show this data?

 

[prodigalson8523]

@ecsid I don’t believe it is perfectly normal but is a reasonable approximation to normal. This is the beauty of Monte Carlo; it smooths out extremes

I’m out at the moment but happy to post the source of the raw performance data when i get back

 

[ecsid]

@prodigalson8523

This is the beauty of Monte Carlo; it smooths out extremes

I suppose this is precisely my concern.

 

[der]

@ecsid Luckily for you, OP has no idea what he's talking about and that's neither a property of or the beauty of Monte Carlo methods. Monte Carlo is useful for things that are hard to model analytically. Normal distributions are, relatively in advanced math, extremely easy to model analytically. There's no need to ever use Monte Carlo on a normal distribution.

 

[prodigalson8523]

Data was from here

https://people.stern.nyu.edu/adamodar/pc/datasets/histretSP.xls

 

[elfblitz]

@prodigalson8523 I'm no data scientist but market returns are typically considered to have excess kurtosis. I'll bet 20 bucks that your model significantly underestimates tail risk.