Geometric Brownian Motion (Random Walk) Process with Drift in Python; Simulate the Future Distribution of Stock Prices in order to Forecast the Future Value of a Stock

11 min readMar 5, 2023

The Random Walk Geometric Brownian Motion process can be used to forecast stock prices, prices of commodities, and other stochastic time-series data given a drift or growth rate and a volatility around the drift path.

Motivation

In order to make a prediction for an unknown future value of a stock, I will show in this article a path-dependent Monte Carlo simulation in Python to simulate future distribution of stock prices using a stochastic process called geometric Brownian motion (random walk) process with drift.

Forecasting is the act of predicting the future, whether it is based on historical data or speculation about the future when no history exists. When historical data exist, a quantitative or statistical approach is best, but if no historical data exist, then a qualitative or judgmental approach is usually the only recourse.

The Basics of Forecasting with Stochastic Processes

A stochastic process is nothing but a mathematically defined equation that can create a series of outcomes over time, outcomes that are not deterministic in nature. That is, it does not follow any simple discernible rule such as price will increase X percent every year or revenues will increase by this factor of X plus Y percent.

A stochastic process is by definition nondeterministic, and one can plug numbers into a stochastic process equation and obtain different results every time. For instance, the path of a stock price is stochastic in nature, and one cannot reliably predict the exact stock price path with any certainty.

However, the price evolution over time is enveloped in a process that generates these prices. The process is fixed and predetermined, but the outcomes are not.

Hence, by stochastic simulation, we create multiple pathways of prices, obtain a statistical sampling of these simulations, and make inferences on the potential pathways that the actual price may undertake given the nature and parameters of the stochastic process used to generate the time series.

Random Walk: Brownian Motion

Assume a process X, where

Figure 2

if and only if

Figure 3

is continuous, where the starting point is

Figure 4

where X is normally distributed with mean zero and variance one or

Figure 5

and were each increment in time is dependent of each other previous increment and is itself normally distributed with mean zero and variance t, such that

Figure 6

Then, the process

Figure 7

follows a Brownian Motion, where α is the drift parameter, σ the volatility measure, and

Figure 8

such that

or X and dX are normally distributed.

If at time zero,

Figure 10

then the expected value of the process X at any time t is such that

Figure 11

and the variance of the process X at any time t is

Figure 12

In the continuous case where there is a drift parameter α, the expected value becomes

Geometric Brownian Motion (Random Walk) Process

The path-dependent geometric Brownian motion (random walk) process for an unknown future value takes the form of

Figure 14

where

where

S0 = the variable’s today value (a.k.a starting value).

T = the forecast horizon (years).

N = the number of steps.

dt = the small interval of time.

St = the variable’s future value from one step to the next

μ = the annualized growth or drift rate.

σ = the annualized volatility.

ε = a random drawing from a standardized normal distribution, Φ(0, 1).

To estimate the parameters from a set of time-series data, the drift rate μ to be the average of the relative returns

while the σ is the standard deviation of all the

values.

Monte Carlo Simulation

The Monte Carlo simulation is usually used in situations where there is no analytical solution and other approaches (such as binomial tree) are less appropriate.

When used to predict the future value of a stock, Monte Carlo simulation uses the real-world valuation result. Real-world valuation means that stocks must grow on average at a higher rate than a risk-free asset. Therefore, the drift rate under a real-world valuation is the expected rate of return from the stock, while the drift rate under a risk-neutral valuation/world is the risk-free rate.

We sample paths to obtain the expected future value in a real-world. Consider a stock S with a value of S0 today, which we want to predict its future value T years from today. Assuming that the expected return and volatility are constant, we can predict the future value of the stock T years from today as follows:

Sample a random path (a.k.a iteration) for S in a real-world.
Calculate the future value of the stock.
Repeat steps 1 and 2 to get many samples values of the future value of the stock in a real-world.
Calculate the mean of the sample future values to get an estimate of the expected future value for the stock.

Suppose that the process for an unknown future value followed by the stock in a real-world is

Figure 18

Geometric Brownian Motion (Random Walk) Process with Drift in Python

Consider a stock with a starting value of 100, drift rate of 5%, annualized volatility of 25% and a forecast horizon of 10 years. Let’s predict the stock future value 10 years from today, based on those parameters and a stochastic process called Brownian motion (random walk) with drift.

Since we wish to run a path-dependent simulation, we will use 100 time steps. Finally, we will use 50,000 iterations in the calculation.

Thus, S0 = 100, T = 10.0, μ = 0.05, σ = 0.25, N=100, M = 50,000

First, let’s load the relevant python libraries for data science.

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from numpy import random as rn
from scipy import stats
import scipy.stats as si
import seaborn as sns
G = 1629562571
np.random.seed(G) 
from datetime import datetime
# datetime object containing current date and time
now = datetime.now()
dt_string = now.strftime("%d/%m/%Y %H:%M:%S")

Second, let’s set the model’s parameters.

S0 = 100
T = 10.0
μ = 0.05
σ = 0.25
M = 50000
N = 100
N

100

Third, let’s build metrics for the future values and the random drawings

ε = rn.randn(M,N)
S = S0*np.ones((M,N+1))
dt = T/N
dt

0.1

Fourth, let’s build the model

start_time = datetime.now() 
for i in range(0,N):
    S[:,i+1] = S[:,i]*(1 + μ*dt + σ*ε[:,i]*np.sqrt(dt))

Fifth, let’s visualize the model results for 8 paths

plt.figure(figsize=(13,7))
fontsize=15
plt.title('Path-Dependent Monte Carlo Simulation - Geometric Brownian Motion Process with Drift',fontsize=fontsize)
plt.xlabel('Years',fontsize=fontsize)
plt.ylabel('Stock prices (USD)',fontsize=fontsize)
plt.grid(axis='y')
a = [ rn.randint(0,M) for j in range(1,8)]
for runer in a:
    plt.plot(np.arange(0,T+dt,dt), S[runer], 'red')

Sixth, let’s present the main results

V = (S[:,-1])
print("\033[1m The estimate of the future value of the stock is ${:.2f}".format(np.mean(V)))
print("\033[1m The accuracy of the estimate of the future value of the stock is ${:.2f}".format((np.std(V)/np.sqrt(M))))

The estimate of the future value of the stock is $164.35
The accuracy of the estimate is $0.67

Seventh, let’s print Intrinsic Value property outputs

time_elapsed = datetime.now() - start_time
def NORMSINV(x):
    x = si.norm.ppf(x)
    return (x)
Workbook_Name = "Geometric Brownian Motion.ipynb"
Number_of_Steps = "{:,.0f}".format(N)
Number_of_Iterations = "{:,.0f}".format(M)
Number_of_Inputs = "{:,.0f}".format(6)
Number_of_Outputs = 1
Sampling_Type = "Latin Hypercube"
Simulation_Start_Time = dt_string
Simulation_Duration = '{}'.format(time_elapsed)
Random_N_Generator = 'Mersenne Twister'
e = ['Workbook Name','Number of Steps','Number of Iterations','Number of Inputs','Number of Outputs','Sampling Type',\
 'Simulation Start Time','Simulation Duration','Random # Generator']
f = [Workbook_Name, Number_of_Steps, Number_of_Iterations, Number_of_Inputs, Number_of_Outputs, Sampling_Type,\
 Simulation_Start_Time, Simulation_Duration, Random_N_Generator]
Per5 = "{:,.1f}".format(np.percentile(V, 5))
P5 = "{:.0%}".format(0.05)
Per10 = "{:,.1f}".format(np.percentile(V, 10))
P10 = "{:.0%}".format(0.10)
Per15 = "{:,.1f}".format(np.percentile(V, 15))
P15 = "{:.0%}".format(0.15)
Per20 = "{:,.1f}".format(np.percentile(V, 20))
P20 = "{:.0%}".format(0.20)
Per25 = "{:,.1f}".format(np.percentile(V, 25))
P25 = "{:.0%}".format(0.25)
Per30 = "{:,.1f}".format(np.percentile(V, 30))
P30 = "{:.0%}".format(0.30)
Per35 = "{:,.1f}".format(np.percentile(V, 35))
P35 = "{:.0%}".format(0.35)
Per40 = "{:,.1f}".format(np.percentile(V, 40))
P40 = "{:.0%}".format(0.40)
Per45 = "{:,.1f}".format(np.percentile(V, 45))
P45 = "{:.0%}".format(0.45)
Per50 = "{:,.1f}".format(np.percentile(V, 50))
P50 = "{:.0%}".format(0.50)
Per55 = "{:,.1f}".format(np.percentile(V, 55))
P55 = "{:.0%}".format(0.55)
Per60 = "{:,.1f}".format(np.percentile(V, 60))
P60 = "{:.0%}".format(0.60)
Per65 = "{:,.1f}".format(np.percentile(V, 65))
P65 = "{:.0%}".format(0.65)
Per70 = "{:,.1f}".format(np.percentile(V, 70))
P70 = "{:.0%}".format(0.70)
Per75 = "{:,.1f}".format(np.percentile(V, 75))
P75 = "{:.0%}".format(0.75)
Per80 = "{:,.1f}".format(np.percentile(V, 80))
P80 = "{:.0%}".format(0.80)
Per85 = "{:,.1f}".format(np.percentile(V, 85))
P85 = "{:.0%}".format(0.85)
Per90 = "{:,.1f}".format(np.percentile(V, 90))
P90 = "{:.0%}".format(0.90)
Per95 = "{:,.1f}".format(np.percentile(V, 95))
P95 = "{:.0%}".format(0.95)
Minimum = "{:,.1f}".format(np.min(V))
Maximum = "{:,.1f}".format(np.max(V))
Mean = "{:,.1f}".format(np.mean(V))
Std_Dev = "{:,.1f}".format(np.std(V))
Variance = int(np.var(V))
Std_Error = "{:,.1f}".format(np.std(V)/np.sqrt(M))
Skewness = round(stats.skew(V),9)
Kurtosis = round((stats.kurtosis(V)+3),9)
Median = "{:,.1f}".format(np.median(V))
Mode = "{:,.1f}".format(stats.mode(V)[0][0])
Left_X = Per5
Left_P = P5
Right_X = Per95
Right_P = P95
Diff_X = "{:,.1f}".format((np.percentile(V, 95) - np.percentile(V, 5)))
Diff_P = "{:.0%}".format(0.90)
Confidence_Level = P95
Lower_Bound = "{:,.1f}".format((np.mean(V) - (np.std(V)/np.sqrt(M))*NORMSINV(0.975)))
Upper_Bound = "{:,.1f}".format((np.mean(V) + (np.std(V)/np.sqrt(M))*NORMSINV(0.975)))
g = {'Information': e, 'Result': f}
st = pd.DataFrame(data=g)
a = ['Minimum','Maximum','Mean','Std Dev','Variance','Std Error', 'Skewness','Kurtosis','Median','Mode',\
 'Left X','Left P','Right X','Right P','Diff X','Diff P','Confidence Level','Lower 95.0%','Upper 95.0%']
b = [Minimum, Maximum, Mean, Std_Dev, Variance, Std_Error, Skewness, Kurtosis, Median, Mode, Left_X, Left_P,\
 Right_X, Right_P, Diff_X, Diff_P, Confidence_Level, Lower_Bound, Upper_Bound]
c = [P5,P10,P15,P20,P25,P30,P35,P40,P45,P50,P55,P60,P65,P70,P75,P80,P85,P90,P95]
d = [Per5, Per10, Per15, Per20, Per25, Per30, Per35, Per40, Per45, Per50, Per55, Per60, Per65,\
 Per70, Per75, Per80, Per85, Per90, Per95]
d = {'Statistics': a, 'Statistics Result': b, 'Percentile': c, 'Percentile Result': d}
st1 = pd.DataFrame(data=d)
from datetime import date
today = date.today()
now = datetime.now()
import calendar
curr_date = date.today()
print("\033[1m Simulation Summary Information")
print("\033[0m ================================================")
print("\033[1m Performed By:","\033[0mIntrinsic Value Team #2")
print("\033[1m Date:","\033[0m",calendar.day_name[curr_date.weekday()],",",today.strftime("%B %d, %Y"),",",now.strftime("%H:%M:%S AM"))
st

print("\033[1m Summary Statistics for Future Value of Stock")
print("\033[0m ======================================================")
print("\033[1m Performed By:","\033[0mIntrinsic Value Team #2")
print("\033[1m Date:","\033[0m",calendar.day_name[curr_date.weekday()],",",today.strftime("%B %d, %Y"),",",now.strftime("%H:%M:%S AM"))
st1

plt.figure(figsize = (4,4))
sns.set_style('white')
sns.set(font_scale = 1.2)
ax = sns.histplot(data=V,bins=50,color='red')
ax.set_xlabel('Values',fontsize=14)
ax.set_xlim( np.percentile(V, 0) , np.percentile(V, 99) )
print("\033[1m Probability Density Function for Future Value of Stock (Sim#1)")
print("\033[0m ======================================================")
print("\033[1m Performed By:","\033[0mIntrinsic Value Team #2")
print("\033[1m Date:","\033[0m",calendar.day_name[curr_date.weekday()],",",today.strftime("%B %d, %Y"),",",now.strftime("%H:%M:%S AM"))

plt.figure(figsize = (4,4))
kwargs = {'cumulative': True}
sns.set_style('white')
sns.set(font_scale = 1.2)
ax = sns.ecdfplot(V, color='red')
ax.set_xlabel('Values',fontsize=14)
ax.set_xlim( np.percentile(V, 0) , np.percentile(V, 99) )
print("\033[1m Cumulative Distribution Function for Future Value of Stock (Sim#1)")
print("\033[0m ======================================================")
print("\033[1m Performed By:","\033[0mIntrinsic Value Team #2")
print("\033[1m Date:","\033[0m",calendar.day_name[curr_date.weekday()],",",today.strftime("%B %d, %Y"),",",now.strftime("%H:%M:%S AM"))

Summary

The Monte Carlo model is widely used in situations where there is no analytical solution and other approaches (such as binomial tree) are less appropriate, such as complex path dependent options (our case). This model uses a simulation of risk factors according to an assumed probability distribution and then calculates the future values along each path separately.

Each future value is then averaged. The average value is referred as the predicted future value. Standard error of sample mean typically serves as an indication of the level of precision of the model.

Since this standard error declines very slow with the number of path typically one need to use dozens of thousands (or even more) of sample paths in order to arrive at an acceptable level of accuracy.

About the Author

Roi Polanitzer, CFV, QFV, FEM, F.IL.A.V.F.A., FRM, CRM, PDS, is a well-known authority in Israel the field of business valuation and has written hundreds of papers that articulate many of the concepts used in modern business valuation around the world. Mr. Polanitzer is the Owner and Chief Appraiser of Intrinsic Value — Independent Business Appraisers, a business valuation firm headquartered in Rishon LeZion, Israel. He is also the Owner and Chief Data Scientist of Prediction Consultants, a consulting firm that specializes in advanced analysis and model development.

Over more than 17 years, he has performed quantitative finance valuation consulting engagements in the areas of: consulting and conducting hedge effectiveness tests, fair value of embedded derivatives, valuations and risk analysis used for 2nd Galai report (VaR sensitivity), valuations and risk analysis for IFRS 7 & IFRS 9, convertible bond valuations, financial options and complex derivative positions, financial and inflation related models (IAS 39, AG7, AG8), credit risk analysis and independent expert opinions for litigation purposes.

Mr. Polanitzer holds an undergraduate degree in economics and a graduate degree in business administration, majoring in finance, both from the Ben-Gurion University of the Negev. He is a Full Actuary (Fellow), a Corporate Finance Valuator (CFV), a Quantitative Finance Valuator (QFV) and a Financial and Economic Modeler (FEM) from the Israel Association of Valuators and Financial Actuaries (IAVFA). Mr. Polanitzer is the Founder of the IAVFA and currently serves as its chairman.

Mr. Polanitzer’s professional recognitions include being designated a Financial Risk Manager (FRM) by the Global Association of Risk Professionals (GARP), a Certified Risk Manager (CRM) by the Israel Association of Risk Managers (IARM), as well as being designated a Python Data Analyst (PDA), a Machine Learning Specialist (MLS), an Accredited in Deep Learning (ADL) and a Professional Data Scientist (PDS) by the Professional Data Scientists’ Israel Association (PDSIA). Mr. Polanitzer is the Founder of the PDSIA and currently serves as its CEO.

He is the editor of IAVFA’s weekly newsletter since its inception (primarily for the professional appraisal community in Israel).

Mr. Polanitzer develops and teaches business valuation professional trainings and courses for the Israel Association of Valuators and Financial Actuaries, and frequently speaks on business valuation at professional meetings and conferences in Israel. He also developed IAVFA’s certification programs in the field of valuation and he is responsible for writing the IAVFA’s statement of financial valuation standards.