Monte Carlo simulations are a technique where a user randomly samples paths for variables to then estimate potential scenarios and likelihood of events. They are particularly useful when we need to value something that does not follow a smooth distribution.

Black-Scholes can be represented as a set of discrete movements by using a difference equation. In this way, we are trying to predict going from time t to time t+1. In this example we will simulate the discrete path as days but in other cases it can be wise to use a much smaller increment. The equation to do this is:

First, to pull random numbers from the standard distribution we can use numpy’s random capability and give it 0 as the first argument (meaning a mean of 0), and 1 for the second argument (meaning standard deviation of 1 because the standard deviation is taken care of in what is multiplied into it). The third argument is 252 meaning to pull 252 samples. We also set the random seed to 1 so that any time someone runs this they get the same “random” numbers. The random seed is useful when you want to get the same random numbers pulled every time. I wouldn’t do this if I wanted truly random numbers, I am just doing it so that everyone will get the same results on their computer.

import pandas as pd
import numpy as np
from scipy.stats import norm
#Set the random seed
np.random.seed(1)
#Sample 252 random number that are uniform from 0 to 1
print(np.random.normal(0,1, 252))

If we use numpy’s reshape function we are able to change a one dimensional array to a two dimensional one. If, for example, we want to get 10 sets of 252, we can run for 252 * 10 values then reshape like so:

#If we multiply 252 by 10 and then re-shape it into an array with dimensions of 252 x 10
np.random.seed(1)
print(np.random.normal(0,1, 252 * 10).reshape(252,10))

Now that we know how to do this, our first step is to get the z-scores from the standard normal distribution (W in the equation).

#Now, we are going to using this for random returns, this first step gets us z-scores
np.random.seed(1)
returns = pd.DataFrame(np.random.normal(0,1, 252 * 10).reshape(252, 10))
print(returns)

The next step is to convert these z-scores into the part of the returns that we will raise e to.

#Set the variables
rf = .02
timestep = 1/252
sigma = .08

#Turn the z-scores to returns
returns = returns * sigma * timestep ** .5 + (rf-.5 * sigma ** 2) * timestep
print(returns)

Now we raise e to this.

#Second part of the conversion to returns
returns = np.exp(returns)
print(returns)

The function cumprod() will find the cumulative product at each point. This will let us see at each point what the total return should be.

#Using cumprod will give us the cummulative product which can be used to track the total return at each point
returns = returns.cumprod()
print(returns)

One final step is to make it so that at time 0 the value is 1 meaning whatever the starting price is. First, move up the index by 1.

#Shift the index forward by 1
returns.index = returns.index + 1
print(returns)

Now set the index of 0 to 1.

#Set the value at 0, the beginning to be 1
returns.loc[0] = 1
returns = returns.sort_index()
print(returns)

And plot simulated paths if we use $100 as the starting price.

import matplotlib.pyplot as plt

#Now we can plot simulated stock paths, multiply in the starting stock price of $100
returns = returns * 100
returns.plot(kind='line')
plt.xlabel("Timestep")
plt.ylabel("Stock Price")
plt.title("Simulated Stock Price Paths")
plt.show()