Monte Carlo Simulation for Derivative Pricing › Simulating Geometric Brownian Motion

To generate sample paths, we need to simulate the stochastic differential equation governing the asset price. The key mathematical tool that enables this simulation is Itô's lemma, which we applied in earlier chapters to derive the behavior of functions of stochastic processes. From Itô's lemma, we derived that if $S_t$ follows geometric Brownian motion, then the terminal price has an exact analytical form:

$$S_T = S_0 \exp\left[\left(r - \frac{\sigma^2}{2}\right)T + \sigma W_T\right]$$

where:

• $S_T$: asset price at time $T$
• $S_0$: initial asset price
• $r$: risk-free interest rate
• $\sigma$: volatility
• $T$: time to maturity
• $W_T$: value of the Brownian motion at time $T$

This equation reveals that the terminal stock price is determined by an exponential function of two terms: a deterministic drift and a random component. The drift term $(r - \sigma^2/2)T$ includes the Itô correction, or "volatility drag," of $-\sigma^2/2$. This adjustment arises from the convexity of the exponential function in stochastic calculus. Without it, we would systematically overestimate expected prices.

Since $W_T \sim \mathcal{N}(0, T)$, meaning the Brownian motion at time $T$ follows a normal distribution with mean zero and variance $T$, we can write $W_T = \sqrt{T} \cdot Z$ where $Z \sim \mathcal{N}(0, 1)$ is a standard normal random variable. This substitution gives us the exact simulation formula for the terminal price:

$$S_T = S_0 \exp\left[\left(r - \frac{\sigma^2}{2}\right)T + \sigma \sqrt{T} \cdot Z\right]$$

where:

• $S_T$: asset price at time $T$
• $S_0$: initial asset price
• $r$: risk-free interest rate
• $\sigma$: volatility
• $T$: time to maturity
• $Z$: standard normal random variable, $Z \sim \mathcal{N}(0, 1)$

The formula decomposes the price evolution into two interpretable parts. The first part, $(r - \sigma^2/2)T$, represents the deterministic drift of the log-price, combining the risk-free rate growth with the volatility drag adjustment. The second part, $\sigma \sqrt{T} \cdot Z$, captures the random diffusion term, scaled by volatility (which amplifies randomness) and by the square root of time (reflecting the well-known result that standard deviation of returns scales with the square root of time horizon).

For European options that depend only on $S_T$, this single-step simulation is all we need. We simply draw a standard normal random variable, plug it into the formula, and obtain one realization of the terminal price. Repeating this process thousands or millions of times gives us the sample distribution needed for Monte Carlo estimation. However, for path-dependent options where the payoff depends on the entire trajectory of prices, we must simulate the complete path from time zero to maturity.

To generate a complete price path, we divide the time interval $[0, T]$ into $M$ discrete steps, each of size $\Delta t = T/M$. At each step, we must advance the price from its current value to the next point in time. The Euler-Maruyama discretization provides a straightforward method for this task by applying a first-order Taylor approximation to the stochastic differential equation:

$$S_{t+\Delta t} = S_t + rS_t \Delta t + \sigma S_t \sqrt{\Delta t} \cdot Z$$

where:

• $S_{t+\Delta t}$: asset price at the next time step
• $S_t$: current asset price
• $r$: risk-free interest rate
• $\sigma$: volatility
• $\Delta t$: time step size ($T/M$)
• $Z$: independent standard normal random variable

The approximation adds two terms to the current price to obtain the next price. The first term, $rS_t \Delta t$, represents the expected growth over the time step, capturing the deterministic drift component. The second term, $\sigma S_t \sqrt{\Delta t} \cdot Z$, introduces the random fluctuation, with the square root of the time step appearing because variance scales linearly with time while standard deviation scales with the square root.

However, this additive form has a subtle but important drawback: it can produce negative prices when $\Delta t$ is large or volatility is high. Since we are adding a random term that can be arbitrarily negative to a positive price, there is always some probability of ending up with a negative value, which is nonsensical for asset prices. A more stable approach uses the exact solution over each interval, applying our closed-form expression for terminal prices but treating each small interval as a mini-simulation:

$$S_{t+\Delta t} = S_t \exp\left[\left(r - \frac{\sigma^2}{2}\right)\Delta t + \sigma \sqrt{\Delta t} \cdot Z\right]$$

where:

• $S_{t+\Delta t}$: asset price at the next time step
• $S_t$: current asset price
• $r$: risk-free interest rate
• $\sigma$: volatility
• $\Delta t$: time step size
• $Z$: independent standard normal random variable

This log-Euler scheme, sometimes called the exact discretization for GBM, preserves positivity unconditionally because the exponential function always returns positive values regardless of its input. Since the exponential of any real number is positive, and we multiply this positive quantity by the positive current price, the next price must be positive. This scheme is also exact for geometric Brownian motion rather than just an approximation, making it the preferred method in practice for simulating GBM paths.

We build our simulations by generating standard normal random variables. As we covered in the probability fundamentals, computers produce uniform random numbers $U \sim \text{Uniform}(0, 1)$ through deterministic algorithms called pseudo-random number generators. These uniform random numbers must then be transformed to the desired distribution through mathematical operations.

The most common approaches for converting uniform random variables to standard normal random variables are:

• Box-Muller transform: This elegant method produces two independent standard normal random variables from two independent uniform random variables. Given two independent uniforms $U_1, U_2$, we compute:

  $$\begin{aligned} Z_1 &= \sqrt{-2\ln U_1}\cos(2\pi U_2) \\ Z_2 &= \sqrt{-2\ln U_1}\sin(2\pi U_2) \end{aligned}$$

  where:

  • $Z_1, Z_2$: independent standard normal random variables
  • $U_1, U_2$: independent uniform random variables on $(0, 1)$

  The transformation works by generating a random point in polar coordinates in a cleverly chosen way. The angle $2\pi U_2$ is uniformly distributed around the circle, while the squared radius $-2\ln U_1$ follows an exponential distribution (because $-\ln U$ is exponentially distributed when $U$ is uniform). Converting these polar coordinates to Cartesian coordinates using cosine and sine yields two independent standard normals. This result can be verified by computing the joint density function and showing it equals the product of two standard normal densities.

• Inverse transform: This approach exploits the fundamental theorem connecting cumulative distribution functions to uniform random variables. We apply the inverse of the standard normal CDF:

  $$Z = \Phi^{-1}(U)$$

  where:

  • $Z$: standard normal random variable
  • $\Phi^{-1}$: inverse of the standard normal cumulative distribution function
  • $U$: uniform random variable on $(0, 1)$

  The intuition here is straightforward: if $U$ represents a cumulative probability (a value between 0 and 1), then $\Phi^{-1}(U)$ tells us which value of a standard normal distribution has that cumulative probability. By mapping a uniformly distributed cumulative probability $U$ to its corresponding quantile $Z$, this method exactly reproduces the standard normal distribution. The method is conceptually elegant and directly applicable, but requires numerical approximation of $\Phi^{-1}$ since no closed-form expression exists for this function.

In practice, numerical libraries like NumPy implement highly optimized algorithms (often the Ziggurat method, which is faster than Box-Muller for modern CPU architectures) that we use directly without needing to implement these transformations ourselves.