Simulating Future Stock Prices using Monte Carlo Methods in Python

Using Monte Carlo methods, we'll write a quick simulation to predict future stock price outcomes for Apple ($AAPL) using Python. You can read more about Monte Carlo simulation (in a finance context) here.

1) Pull the data

First, we can import the libraries, and pull the historical stock data for Apple. For this example, I picked the last ~10 years, although it would be valuable to test sensitivities of different ranges as this alone is subjective:

# Import required libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from pandas_datareader import data
apple = data.DataReader('AAPL', 'yahoo',start='1/1/2009')
apple.head()

	High	Low	Open	Close	Volume	Adj Close
Date
2009-01-02	13.005714	12.165714	12.268572	12.964286	186503800.0	11.357091
2009-01-05	13.740000	13.244286	13.310000	13.511429	295402100.0	11.836406
2009-01-06	13.881429	13.198571	13.707143	13.288571	322327600.0	11.641174
2009-01-07	13.214286	12.894286	13.115714	13.001429	188262200.0	11.389630
2009-01-08	13.307143	12.862857	12.918571	13.242857	168375200.0	11.601127

#Next, we calculate the number of days that have elapsed in our chosen time window
time_elapsed = (apple.index[-1] - apple.index[0]).days

2) Calculate the compounded annualized growth rate over the length of the dataset + standard deviation (to feed into simulation)

#Current price / first record (e.g. price at beginning of 2009)
#provides us with the total growth %
total_growth = (apple['Adj Close'][-1] / apple['Adj Close'][1])

#Next, we want to annualize this percentage
#First, we convert our time elapsed to the # of years elapsed
number_of_years = time_elapsed / 365.0
#Second, we can raise the total growth to the inverse of the # of years
#(e.g. ~1/10 at time of writing) to annualize our growth rate
cagr = total_growth ** (1/number_of_years) - 1

#Now that we have the mean annual growth rate above,
#we'll also need to calculate the standard deviation of the
#daily price changes
std_dev = apple['Adj Close'].pct_change().std()

#Next, because there are roughy ~252 trading days in a year,
#we'll need to scale this by an annualization factor
#reference: https://www.fool.com/knowledge-center/how-to-calculate-annualized-volatility.aspx

number_of_trading_days = 252
std_dev = std_dev * math.sqrt(number_of_trading_days)

#From here, we have our two inputs needed to generate random
#values in our simulation
print ("cagr (mean returns) : ", str(round(cagr,4)))
print ("std_dev (standard deviation of return : )", str(round(std_dev,4)))

cagr (mean returns) : 0.3746 std_dev (standard deviation of return : ) 0.2878

3) Generate random values, run Monte Carlo simulation

#Generate random values for 1 year's worth of trading (252 days),
#using numpy and assuming a normal distribution
daily_return_percentages = np.random.normal(cagr/number_of_trading_days, std_dev/math.sqrt(number_of_trading_days)number_of_trading_days)+1

#Now that we have created a random series of future
#daily return %s, we can simply apply these forward-looking
#to our last stock price in the window, effectively carrying forward
#a price prediction for the next year

#This distribution is known as a 'random walk'

price_series = [apple['Adj Close'][-1]]

for j in daily_return_percentages:
    price_series.append(price_series[-1] * j)

#Great, now we can plot of single 'random walk' of stock prices
plt.plot(price_series)
plt.show()

#Now that we've created a single random walk above,
#we can simulate this process over a large sample size to
#get a better sense of the true expected distribution
number_of_trials = 3000

#set up an additional array to collect all possible
#closing prices in last day of window.
#We can toss this into a histogram
#to get a clearer sense of possible outcomes
closing_prices = []

for i in range(number_of_trials):
    #calculate randomized return percentages following our normal distribution
    #and using the mean / std dev we calculated above
    daily_return_percentages = np.random.normal(cagr/number_of_trading_days, std_dev/math.sqrt(number_of_trading_days),number_of_trading_days)+1
    price_series = [apple['Adj Close'][-1]]

    for j in daily_return_percentages:
        #extrapolate price out for next year
        price_series.append(price_series[-1] * j)

    #append closing prices in last day of window for histogram
    closing_prices.append(price_series[-1])

    #plot all random walks
    plt.plot(price_series)

plt.show()

#plot histogram
plt.hist(closing_prices,bins=40)

plt.show()

4) Analyze results

#from here, we can check the mean of all ending prices
#allowing us to arrive at the most probable ending point
mean_end_price = round(np.mean(closing_prices),2)
print("Expected price: ", str(mean_end_price))

Expected price: 193.85

#lastly, we can split the distribution into percentiles
#to help us gauge risk vs. reward

#Pull top 10% of possible outcomes
top_ten = np.percentile(closing_prices,100-10)

#Pull bottom 10% of possible outcomes
bottom_ten = np.percentile(closing_prices,10);

#create histogram again
plt.hist(closing_prices,bins=40)
#append w/ top 10% line
plt.axvline(top_ten,color='r',linestyle='dashed',linewidth=2)
#append w/ bottom 10% line
plt.axvline(bottom_ten,color='r',linestyle='dashed',linewidth=2)
#append with current price
plt.axvline(apple['Adj Close'][-1],color='g', linestyle='dashed',linewidth=2)

plt.show()

Here, we can see (based solely on using Monte Carlo simulation, of course) there looks to be more upside than downside for the next year, with the expected price running about $193 and only a 10% chance of the price landing below $128.

There are of course many, many things that drive a stock price beyond historical percent changes (meaning I would not make investment decisions solely off Monte Carlo simulation!), but this offers a thorough example of using Monto Carlo simulation to better understand the distribution of possible outcomes.