Monte Carlo Simulation With Brownian Motion: Modelling Price Risk the Right Way

Treasury teams, investors and risk managers rarely need to know the price of an asset next year. They need to know the range: how high it could go, how low it could fall and how likely each outcome is. Monte Carlo simulation answers that question by generating thousands of possible price paths and reading the distribution of results.

The most common model behind those paths is Geometric Brownian Motion (GBM). It assumes that prices move randomly around a trend, described by two parameters: the drift, the average rate at which the price changes, and the volatility, how much it moves around that trend. It also assumes that logarithmic returns are normally distributed.

To show how it works, we simulate Apple's share price from 2021 to 2024 and compare the simulated paths with what actually happened.

Method one: the Black-Scholes formula

The Black-Scholes framework assumes log-normal returns, which gives the standard expression for the price at time t:

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

Here μ is the drift, σ the volatility and W(t) a standard Brownian motion. The term -0.5 * σ² corrects for volatility drag, and it is part of the formula.

The code below estimates the drift and volatility from Apple's daily returns and simulates 1,000 paths.

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)

Method two: model the log-returns directly

If returns are log-normal, their logarithms are normal, and they can be simulated with a standard Brownian motion and then turned back into prices. It is the same model seen through a different lens, and it is often easier to explain.

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)

What the comparison shows

Both methods produce practically identical distributions, which confirms that they are two valid routes to the same model. For decision makers, the value is not in the individual paths but in the range they describe: the band within which the price is likely to move, and the probability of reaching a level that would hurt the business.