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.

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

#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

#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()

#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 percent 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.

Get better at data science interviews by solving a few questions per week

Learn more