Monte Carlo Option Pricing with GBM

1. Geometric Brownian Motion — Stock Price Dynamics

Under the risk-neutral measure, the stock price follows a Geometric Brownian Motion (GBM):

$$dS = (r - q)\, S \, dt + \sigma \, S \, dW_t$$

Applying Itô's Lemma to $\ln S$ and integrating over $[0, T]$ gives the exact solution:

$$S_T = S_0 \; \exp\!\left[\left(r - q - \tfrac{\sigma^2}{2}\right)T + \sigma\sqrt{T}\; Z\right]$$

where $Z \sim \mathcal{N}(0, 1)$.

Log-Return Decomposition

The continuously compounded return is:

$$\ln\!\left(\frac{S_T}{S_0}\right) = \underbrace{\left(r - q - \tfrac{\sigma^2}{2}\right)T}_{\text{drift}} + \underbrace{\sigma \sqrt{T}\; Z}_{\text{diffusion}}$$

The drift term includes a $-\tfrac{\sigma^2}{2}$ correction — this is the Itô adjustment that ensures $\mathbb{E}[S_T] = S_0 \, e^{(r-q)T}$ under the risk-neutral measure. The drift is the deterministic component that governs the expected direction of the stock price over time.
The diffusion term scales the random shock $Z$ by $\sigma\sqrt{T}$. The diffusion is the stochastic component that introduces randomness into each simulated path.

2. Terminal Price Distribution — Volatility Comparison

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# ── Parameters ──
S0    = 50            # Initial stock price
K     = 50            # Strike price
T     = 9 / 12        # Time to maturity (years)
r     = 0.04          # Risk-free interest rate
q     = 0.00          # Continuous dividend yield
sigma = 0.30          # Volatility (base)
N     = 10_000_000    # Number of simulations

np.random.seed(42)
z = np.random.normal(0, 1, N)

st = S0 * np.exp((r - q - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)

# Terminal prices at two different implied volatilities
sigma_low  = 0.20
sigma_high = 0.30

S_T_low  = S0 * np.exp((r - q - 0.5 * sigma_low**2)  * T + sigma_low  * np.sqrt(T) * z)
S_T_high = S0 * np.exp((r - q - 0.5 * sigma_high**2) * T + sigma_high * np.sqrt(T) * z)

# Both share the same risk-neutral mean
mean_both   = np.mean(S_T_low)
median_low  = np.median(S_T_low)
median_high = np.median(S_T_high)

fig, ax = plt.subplots(figsize=(10, 5))

ax.hist(S_T_low,  bins=100, density=True, alpha=0.55, color="steelblue",  label=f"sigma = {sigma_low:.0%}")
ax.hist(S_T_high, bins=100, density=True, alpha=0.45, color="indianred", label=f"sigma = {sigma_high:.0%}")

ax.axvline(mean_both, color="black", linestyle="-", linewidth=1.3, label=f"Mean = {mean_both:.2f}")
ax.axvline(median_low,  color="steelblue", linestyle="--", linewidth=1.2, label=f"Median (20%) = {median_low:.2f}")
ax.axvline(median_high, color="indianred", linestyle="--", linewidth=1.2, label=f"Median (30%) = {median_high:.2f}")

ax.set_title("Terminal Price Distribution -- sigma = 20% vs 30%", fontweight="bold")
ax.set_xlabel("$S_T$")
ax.set_ylabel("Density")
ax.legend()
plt.tight_layout()
plt.show()

3. Option Pricing — Monte Carlo vs Black-Scholes

# Monte Carlo simulation (single-step terminal value)

# Call
payoff_call  = np.maximum(st - K, 0)
premium_call = np.mean(payoff_call) * np.exp(-r * T)
se_call      = np.std(payoff_call, ddof=1) / np.sqrt(N)

# Put
payoff_put  = np.maximum(K - st, 0)
premium_put = np.mean(payoff_put) * np.exp(-r * T)
se_put      = np.std(payoff_put, ddof=1) / np.sqrt(N)

# Black-Scholes analytical prices (for comparison)
d1 = (np.log(S0 / K) + (r - q + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
bs_call = S0 * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
bs_put  = K * np.exp(-r * T) * norm.cdf(-d2) - S0 * np.exp(-q * T) * norm.cdf(-d1)

print(premium_call, bs_call)
print(premium_call - bs_call)

4. Monte Carlo Simulation Paths

## Paths matrix, now we go with two dimensions -> (n_paths, n_steps)
n_paths = 175
n_steps = int(T * 252)
dt = T / n_steps

np.random.seed(7)
Z = np.random.normal(0, 1, (n_paths, n_steps))

# Build price matrix: each row is one path
paths = np.zeros((n_paths, n_steps + 1))
paths[:, 0] = S0

# Paths composition and terminal value
for t in range(n_steps):
    paths[:, t + 1] = paths[:, t] * np.exp(
        (r - q - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z[:, t]
    )

S_terminal = paths[:, -1]

# Classification of terminal values, call option
atm_band = 0.5
is_itm = S_terminal > (K + atm_band)
is_otm = S_terminal < (K - atm_band)
is_atm = ~is_itm & ~is_otm

days = np.arange(n_steps + 1)

plt.figure(figsize=(10, 5))

for i in np.where(is_otm)[0]:
    plt.plot(days, paths[i], color="red", alpha=0.15, linewidth=0.5)

for i in np.where(is_atm)[0]:
    plt.plot(days, paths[i], color="gold", alpha=0.6, linewidth=0.8)

for i in np.where(is_itm)[0]:
    plt.plot(days, paths[i], color="green", alpha=0.15, linewidth=0.5)

plt.axhline(K, color="black", linestyle="--", linewidth=1.2, label=f"Strike K = {K}")

plt.ylabel("Stock Price ($)")
plt.xlabel("Days")
plt.title(f"MC simulation -- {n_paths} paths (sigma = {sigma:.0%})")
plt.legend([
    plt.Line2D([], [], color="green",  linewidth=2),
    plt.Line2D([], [], color="gold",   linewidth=2),
    plt.Line2D([], [], color="red",    linewidth=2),
    plt.Line2D([], [], color="black",  linewidth=1.2, linestyle="--"),
], [
    f"ITM  ({np.sum(is_itm)})",
    f"ATM  ({np.sum(is_atm)})",
    f"OTM  ({np.sum(is_otm)})",
    f"Strike = {K}",
], loc="upper left")
plt.tight_layout()
plt.show()

Monte Carlo Option Pricing — GBM

1. Geometric Brownian Motion — Stock Price Dynamics

Log-Return Decomposition

2. Terminal Price Distribution — Volatility Comparison

3. Option Pricing — Monte Carlo vs Black-Scholes

4. Monte Carlo Simulation Paths

David Arias, CFA