Monte Carlo Simulations for Predicting Stock Prices [Python] (2024)

Monte Carlo Simulations are an incredibly powerful tool in numerous contexts, including operations research, game theory, physics, business and finance, among others. It is a technique used to understand the impact of risk and uncertainty when making a decision. Simply put, a Monte Carlo simulation runs an enourmous amount of trials with different random numbers generated from an underlying distribution for the uncertain variables.

Here, we will dive into how to predict stock prices using a Monte Carlo simulation!

What do we need to understand before we start?

We want to predict the price of the stock today. We know the price of the stock yesterday. So, what’s missing? We obviously do not know the daily return that the stock is going to yield today. This is where Monte Carlo comes in! But first… how can we estimate the return?

How do we predict the daily return of the stock? Brownian Motion.

Brownian motion will be the main driver for estimating the return. It is a stochastic process used for modeling random behavior over time. Brownian motion has two main components:

1. Drift — the direction that rates of returns have had in the past. That is, the expected return of the stock. You may ask yourself: why is the variance multiplied by 0.5? Becasue historical values are eroded in the future.
2. Volatility — the historical volatility multiplied by a random, standard normal variable.

Together, these compute the brownian motion — ie the daily return of a stock!

This technique will be used for every day into the future you want to predict, and for however many trials the monte carlo simulation will run! Let’s get to coding!

import numpy as np
import pandas as pd
from pandas_datareader import data as wb
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm

Step 1: Import the stock data. We will use Google (GOOG) as the ongoing example.

ticker = 'GOOG'
data = pd.DataFrame()
data[ticker] = wb.DataReader(ticker, data_source = 'yahoo', start = '2010-1-1')['Adj Close]#Plot
data.plot(figsize=(15,6))

This method is straightforward and simple. Using pandas_datareader, you can import stock data for free and without access keys. Highly convenient for our purposes!

Step 2: Compute the logarithmic returns of Google stock

log_return = np.log(1 + data.pct_change())#Plot
sns.distplot(log_returns.iloc[1:])
plt.xlabel("Daily Return")
plt.ylabel("Frequency")

Step 3: Compute the Drift

u = log_returns.mean()
var = log_returns.var()
drift = u - (0.5*var)

Step 4: Compute the Variance and Daily Returns

In this step we have to generate random variables for every day forecasted and for every simulation trial we will run.

stdev = log_returns.std()
days = 50
trials = 10000
Z = norm.ppf(np.random.rand(days, trials)) #days, trials
daily_returns = np.exp(drift.values + stdev.values * Z)

So close! Now that we have randomly generated 50 random variables for each of the days in every one of the ten thousand trials, all we need is to calculate the price path for each of the trials!

Step 5: Calculating the stock price for every trial

price_paths = np.zeros_like(daily_returns)
price_paths[0] = data.iloc[-1]
for t in range(1, days):
 price_paths[t] = price_paths[t-1]*daily_returns[t]

Does the above loop look familiar? Exaclty! It’s simply our stock price equation, the first one we saw in this article!

There it goes! We’ve ran a Monte Carlo simulation that predicts Google’s stock price 50 days into the future.

With the price_paths matrix, now you can calculate the probability of profitability, or the expected annualized returm, for example. We do both in the next section!

In this section, we briefly go over the functions created to complete the steps from the above process so that we can automate the simulations with as a many stocks as you want! We will also display some metrics from the Capital Asset Pricing Model.

Step 1: Importing stock data

def import_stock_data(tickers, start = '2010-1-1'):
 data = pd.DataFrame()
 if len([tickers]) ==1:
 data[tickers] = wb.DataReader(tickers, data_source='yahoo', start = start)['Adj Close']
 data = pd.DataFrame(data)
 else:
 for t in tickers:
 data[t] = wb.DataReader(t, data_source='yahoo', start = start)['Adj Close']
 return(data)#Example use
data = import_stock_data(["FB","GOOG","AAPL"], start = '2005-1-1')

Step 2: Compute Logarithmic Daily Returns

def log_returns(data):
 return (np.log(1+data.pct_change()))#Example use
log_return = log_returns(data)

Step 3: Calculate Drift

def drift_calc(data):
 lr = log_returns(data)
 u = lr.mean()
 var = lr.var()
 drift = u-(0.5*var)
 try:
 return drift.values
 except:
 return drift#Example use
drift_calc(data)

Step 4: Compute Daily Returns

def daily_returns(data, days, iterations):
 ft = drift_calc(data, return_type)
 try:
 stv = log_returns(data).std().values
 except:
 stv = log_returns(data).std() dr = np.exp(ft + stv * norm.ppf(np.random.rand(days, iterations))) return dr#Example use
daily_returns(data, 50, 1000)

Step 5 (CAPM): Compute the Sharpe Ratio, the risk-adjusted returns and the Beta of every stock.

As we can observe from the equations, we must compare each stock against the market. So we will import the market data (S&P 500). We also need the risk-free rate. We will approximate this rate by the yield of a 10-year US bond (2.5%).

def beta_sharpe(data, mark_ticker = "^GSPC", start='2010-1-1', riskfree = 0.025):
 # Beta
 dd, mark_ret = market_data_combination(data, mark_ticker, start)
 log_ret = log_returns(dd)
 covar = log_ret.cov()*252 # Annualized
 covar = pd.DataFrame(covar.iloc[:-1,-1])
 mrk_var = log_ret.iloc[:,-1].var()*252 #Annualized
 beta = covar/mrk_var
 stdev_ret = pd.DataFrame(((log_ret.std()*250**0.5)[:-1]), columns=['STD'])
 beta = beta.merge(stdev_ret, left_index=True, right_index=True)
 # CAPM
 for i, row in beta.iterrows():
 beta.at[i,'CAPM'] = riskfree + (row[mark_ticker] * (mark_ret-riskfree))
 # Sharpe
 for i, row in beta.iterrows():
 beta.at[i,'Sharpe'] = ((row['CAPM']-riskfree)/(row['STD']))
 beta.rename(columns={"^GSPC":"Beta"}, inplace=True)
 return beta#Example use
#Make sure the start date here is the same as the start day of the original data.
beta_sharpe(data, '2005-1-1') 

Step 7 (Probability): Calculate the probability of a certain outcome

Once we run the monte carlo simulation for several stocks, we may want to calculate the probability of our investment having a positive return, or 25% return, or simply the probability that the stock reach a specific price. So, we created this equation that does it for you!

def probs_find(predicted, higherthan, on = 'value'):
 if on == 'return':
 predicted0 = predicted.iloc[0,0]
 predicted = predicted.iloc[-1]
 predList = list(predicted)
 over = [(i*100)/predicted0 for i in predList if ((i-predicted0)*100)/predicted0 >= higherthan]
 less = [(i*100)/predicted0 for i in predList if ((i-predicted0)*100)/predicted0 < higherthan]
 elif on == 'value':
 predicted = predicted.iloc[-1]
 predList = list(predicted)
 over = [i for i in predList if i >= higherthan]
 less = [i for i in predList if i < higherthan]
 else:
 print("'on' must be either value or return")
 return (len(over)/(len(over)+len(less)))#Example use (probability our investment will return at least 20% over the days specified in our prediction
probs_find(predicted, 0.2, on = 'return')

Step 8: Run the Monte Carlo simulation for a single stock

This step is passive. We define the function so we can loop through it on the last step.

def simulate_mc(data, days, iterations, plot=True): # Generate daily returns
 returns = daily_returns(data, days, iterations)
 # Create empty matrix
 price_list = np.zeros_like(returns)
 # Put the last actual price in the first row of matrix. 
 price_list[0] = data.iloc[-1]
 # Calculate the price of each day
 for t in range(1,days):
 price_list[t] = price_list[t-1]*returns[t]
 # Plot Option
 if plot == True:
 x = pd.DataFrame(price_list).iloc[-1]
 fig, ax = plt.subplots(1,2, figsize=(14,4))
 sns.distplot(x, ax=ax[0])
 sns.distplot(x, hist_kws={'cumulative':True},kde_kws={'cumulative':True},ax=ax[1])
 plt.xlabel("Stock Price")
 plt.show()
 #CAPM and Sharpe Ratio
 # Printing information about stock
 try:
 [print(nam) for nam in data.columns]
 except:
 print(data.name)
 print(f"Days: {days-1}")
 print(f"Expected Value: ${round(pd.DataFrame(price_list).iloc[-1].mean(),2)}")
 print(f"Return: {round(100*(pd.DataFrame(price_list).iloc[-1].mean()-price_list[0,1])/pd.DataFrame(price_list).iloc[-1].mean(),2)}%")
 print(f"Probability of Breakeven: {probs_find(pd.DataFrame(price_list),0, on='return')}")
 return pd.DataFrame(price_list)#Example use
simulate_mc(data, 252, 1000)

Step 9: Run the complete Monte Carlo simulation for as many stocks as needed!

def monte_carlo(tickers, days_forecast, iterations, start_date = '2000-1-1', plotten=False): data = import_stock_data(tickers, start=start_date) inform = beta_sharpe(data, mark_ticker="^GSPC", start=start_date) simulatedDF = []
 for t in range(len(tickers)):
 y = simulate_mc(data.iloc[:,t], (days_forecast+1), iterations)
 if plotten == True:
 forplot = y.iloc[:,0:10]
 forplot.plot(figsize=(15,4))
 print(f"Beta: {round(inform.iloc[t,inform.columns.get_loc('Beta')],2)}")
 print(f"Sharpe: {round(inform.iloc[t,inform.columns.get_loc('Sharpe')],2)}") 
 print(f"CAPM Return: {round(100*inform.iloc[t,inform.columns.get_loc('CAPM')],2)}%")
 y['ticker'] = tickers[t]
 cols = y.columns.tolist()
 cols = cols[-1:] + cols[:-1]
 y = y[cols]
 simulatedDF.append(y)
 simulatedDF = pd.concat(simulatedDF)
 return simulatedDF#Example use
ret_sim_df = monte_carlo(['GOOG','FB','AAPL'], 252, 10000, start_date='2015-1-1')

Once we run the simulation, we can observe the distribution of values at the end date, which depends on the amount of days forecasted. The metrics for each stock, including the average return, the risk-adjusted return, the probability of breaking even, the Beta, and the Sharpe ratio are also printed. Hopefully, this will help you as an investor!

Voila! Just by loading/running the functions created in steps 1–9, you can run your own Monte Carlo simulation for the stocks you want! You’re only required to have the appropiate libraries installed.

It’s important to note that the Monte Carlo simulations here assume a normal distribution of returns. However, history has shown us that this is not precisely true. While similar, the most precise distribution is the Cauchy distribution. Nonetheless, for most purposes, a normal distribution is accurate enough over thousands of trials.

Hope you enjoyed this introduction to Monte Carlo simulations for predicting stock prices! Feel free to download and use the Jupyter notebook from Github and follow me on LinkedIn and Medium for other projects!

Github: https://github.com/eliasmelul/finance_portfolio

LinkedIn: https://www.linkedin.com/in/melulelias/

Medium: https://medium.com/@melulfrescoelias