How do I plot the average of a random walk/monte carlo sim Python

Question

I have a problem where I am plotting a random walk and I am looking to get the average line for all the simulations. My code is below:

import numpy as np
from math import sqrt
import matplotlib.pyplot as plt

# starting stock price
S = 100

# number of trading days
T = 252

# random data inputs
mu = 0.061
stdev = 0.165

if __name__ == '__main__':

    average= []

    # running the simulation multiple times (# specified in range) times
    for i in range(100):

        daily_returns = np.random.normal((mu/T), stdev/sqrt(T), T) + 1

        # set starting price and create price series generated by above random daily returns
        price_list = [S]

        for x in daily_returns:
            price_list.append(price_list[-1]*x)

        # Generate Plots - price series and histogram of daily returns
        plt.plot(price_list, color='gray')
        plt.hist(daily_returns-1, 100) # Note that we run the line plot and histogram separately, not simultaneously.
        average.append(np.mean(price_list))
    plt.plot(average, color='red')
    plt.show()

The issue I am having is the average line I've been able to come up (probably incorrectly) seems to stop halfway through the plot. I am sure this is an easy fix but it is going to drive me crazy!

Thanks!

Answer 1

The reason why it's failing because you run the simulation 100 times, so len(avarage) will be 100, but len(price_list) is always 252 + 1. The easiest way to fix is to make those two identical. But that won't fix another huge problem: you calculate the avarage price of 252 + 1 days each simulation, so that's why your avarage is wrong at the starting day. You should be taking avarage by days. A better solution would be:

import numpy as np
import matplotlib.pyplot as plt

S = 100
T = 10
mu = 0.061
stdev = 0.165
SIMULATIONS = 100

if __name__ == '__main__':
    # the array to store simulation results
    full_array = np.empty(shape=(SIMULATIONS, T + 1))

    for i in range(SIMULATIONS):
        daily_returns = np.random.normal((mu/T), stdev/np.sqrt(T), T) + 1
        
        # A more efficient way of calculating the same as you have
        # It's a simple geometric series.
        price_list = np.empty(shape=len(daily_returns) + 1)
        price_list[0] = S
        price_list[1:] = daily_returns
        price_list = np.cumprod(price_list)
        plt.plot(price_list, color="gray")
        
        # save that simulation run
        full_array[i, :] = price_list

# plot the mean by days
plt.plot(np.mean(full_array, axis=0), color="red")