top of page

Monte Carlo Simulation

Updated: Jul 28, 2022

  • A Monte Carlo simulation is a model used to predict the probability of different outcomes when the intervention of random variables is present.

  • Monte Carlo simulations help to explain the impact of risk and uncertainty in prediction and forecasting models.

  • A variety of fields utilize Monte Carlo simulations, including finance, engineering, supply chain, and science.

  • The basis of a Monte Carlo simulation involves assigning multiple values to an uncertain variable to achieve multiple results and then averaging the results to obtain an estimate.

  • Monte Carlo simulations assume perfectly efficient markets.

If you need to solve a problem that you cannot figure out, chances are Monte Carlo simulations can be used to get you pretty close to correct. Recently I was trying to find a way to optimize a complex strategy that utilizes spreads of spreads in Futures markets, and this optimization algorithm was killing me. I decided that I could set the spread ratio as a random variable, and run it as a Monte Carlo simulation, and at least get in the right direction — in 5 minutes & 100k iterations I had a simple 15-line solution to a problem that had taken me maybe 350 lines of Python when I initially tried to use a minimization function. This is the adjustable wrench of your toolbox.


Let’s dive in, and I’m going to over-comment this code so it could not be clearer what’s doing what. First off, we need your security returns. This can be pulled from system performance (Like my example) or pulled with quantrautil using simple price data. The DF should be Securities Tickers as columns, with date rows and daily return values. Technically, it should be log returns so I’ve included that calculation — but the difference is typically immaterial.

Begin by initializing the arrays for saving the performance values from the Monte Carlo runs, and set up the Monte Carlo simulations loop, defining the # of runs as high as your PC will allow (start with 1k or so, scale up).


The important thing to remember is the weights, this is the Monte Carlo simulations magic, as this will select a new value each time and provide the power to the optimization, making each run unique. Also notice how each value your saving (ret, vol, Sharpe) is really an array being saved with indexing per each run — so run 0 begins, randomizes weights, and calculates the return, saving it to index 0 in ret arr, and then vol saves in 0 index, and finally SR. Once you’ve looped through, you run an argmax() function in the SR arr (or whatever value you maximizing), which will give you weight, 477 in my example. This mean run 477 gave you the highest SR value, and that is your ideal weight! You can find the optimal values in all weights, indexing with your max run number.

(I also included a helper function to save the weights and tickers into a DF, and Pkl it for reference).


import numpy as np
import numpy.random as nrand
import matplotlib.pyplot as plt
import math

#Define vars
S = 1000000
T = 252
mu = 1.92246
vol = .86123  #Between .76 and .86

result = []
for i in range(1000):
    daily_returns = np.random.normal(mu/T,vol/math.sqrt(T),T)+1
    
    price_list = [S]
    
    for x in daily_returns:
        price_list.append(price_list[-1]*x)
    result.append(price_list[-1]) #appending each runs end value --to calculate the mean return
        
    plt.plot(price_list)#This is key, KEEP THIS IN LOOP,votherwise it will plot one iteration/return path.
plt.title('Daily MC')
plt.show()

from scipy import stats
from statistics import stdev

print('Mean:',np.mean(result))
print('Mean Ret:',np.mean(result)/S*100)
print('Median:',np.median(result))
print('Median Ret:',np.median(result)/S*100)
print('Min:',np.min(result))
print('Min Ret:',np.min(result)/S*100)
print('Max:',np.max(result))
print('Max Ret:',np.max(result)/S*100)
#print('Mode:',stats.mode(result))
print('Stdev', stdev(result))
mc_mu = np.mean(result)
med = np.median(result)
mc_min = np.min(result)
mx = np.max(result)
std = stdev(result)

metrics = [mc_mu,med,mc_min,mx]
print('sharpe:',mu/vol)

print('5% Quantile',np.percentile(result,5))
print('5% Quantile %',np.percentile(result,5)/S*100)
print('95% Quantile',np.percentile(result,95))
print('95% Quantile %',np.percentile(result,95)/S*100)