How to Use Monte Carlo Simulation with Brownian Motion

In this article we implement a Monte Carlo simulation of stock prices using Geometric Brownian Motion (GBM) in a simple way, explaining the underlying logic so you can run it correctly.

In mathematical finance, an asset is often modelled as following a Geometric Brownian Motion, which means its price obeys a stochastic differential equation (SDE). GBM assumes the price moves up and down randomly, but with a general trend of rising or falling over time. It captures this with two ingredients: how fast the price changes on average (the drift) and how much it randomly jumps around (the volatility). It also assumes that the logarithmic returns of the asset follow a normal distribution.

For this example we use Apple's stock price and run a GBM Monte Carlo simulation over the same period. We then plot the simulations alongside the stock's actual path.

Two ways to model GBM

1. The Black-Scholes approach

The Black-Scholes framework assumes that asset returns follow a log-normal distribution. Applying the maths gives the standard formula for the asset price:

S(t) = S(0) * exp((μ - 0.5 * σ²) * t + σ * W(t))

Here μ is the drift (the annualised average return of the stock), σ is the volatility (the standard deviation of returns) and W(t) is a standard Brownian motion. The term -0.5 * σ² is a correction for the “drag” effect on returns, and it is built into the formula.

GBM Monte Carlo simulation for Apple stock (2021-2024)

Now let's put it into practice. The code below simulates Apple's price from 2021 to 2024 in Python.

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt

apple = yf.download("AAPL", start="2021-01-01", end="2024-01-01")
apple.columns = apple.columns.levels[0]
price = apple.Close
returns = price.pct_change().dropna()

t = len(returns)           # number of daily steps
dt = 1                     # one trading day per step
mu = returns.mean()        # drift
sigma = returns.std()      # volatility
simulations = 1000
S0 = price.iloc[0]

def gbm_paths(mu, sigma, W):
    return np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * W)

Ws = np.random.standard_normal(size=(simulations, t))
gbm = gbm_paths(mu, sigma, Ws)
gbm = np.insert(gbm, 0, S0, axis=1).cumprod(axis=1)

x_axis = range(0, t + 1)
fig, ax = plt.subplots(figsize=(12, 5.5))
for path in gbm:
    ax.plot(x_axis, path, color="grey", linewidth=0.1)
ax.plot(x_axis, price, color="red", linewidth=1)

2. The log-transformation approach

If returns are log-normal, we can simply take the logarithm of the returns. After that transformation the data can be modelled with a standard Brownian motion. It leads to the same underlying model, seen through a different mathematical lens.

log_returns = np.log(returns + 1)
mu_log = log_returns.mean()
sigma_log = log_returns.std()

def bm_paths(mu, sigma, W):
    return mu + sigma * W

bm = bm_paths(mu_log, sigma_log, Ws)
bm = S0 * np.exp(bm.cumsum(axis=1))
bm = np.insert(bm, 0, S0, axis=1)

fig, ax = plt.subplots(figsize=(12, 5.5))
for path in bm:
    ax.plot(x_axis, path, color="grey", linewidth=0.1)
ax.plot(x_axis, price, color="red", linewidth=1)

Comparing both results

The two charts are practically identical. With perfectly normal log-returns they converge to the same result, which confirms that both methods are valid ways to run a Geometric Brownian Motion Monte Carlo simulation.