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')
High Low Open Close Volume Adj Close
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

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

q184 p1

#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

    #plot all random walks

#plot histogram

q184 p2

q184 p3

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
#append w/ top 10% line
#append w/ bottom 10% line
#append with current price
plt.axvline(apple['Adj Close'][-1],color='g', linestyle='dashed',linewidth=2)

q184 p4

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.