Finite Difference Methods for Option Pricing › The Black-Scholes PDE on a Grid

The foundation of finite difference methods is the Black-Scholes PDE, which is derived using Itô's lemma and the principle of no-arbitrage pricing. Recall that this partial differential equation emerges from the requirement that a perfectly hedged portfolio of an option and its underlying asset must earn the risk-free rate. The equation captures how the option value $V(S,t)$ evolves as time passes and as the underlying asset price $S$ fluctuates. For an option with value $V(S,t)$ on an underlying asset with price $S$, the PDE is:

$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0$$

where:

• $V(S,t)$: option price as a function of asset price $S$ and time $t$
• $r$: risk-free interest rate
• $\sigma$: volatility of the underlying asset
• $\frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2}$: diffusion term capturing the effect of volatility
• $rS\frac{\partial V}{\partial S}$: drift term reflecting the risk-neutral expected return
• $-rV$: discounting term reflecting the time value of money

Each term in this equation has a clear financial interpretation that helps us understand why the option value must satisfy this relationship. The time derivative $\frac{\partial V}{\partial t}$ captures how the option value changes purely due to the passage of time, with all else held constant. The second derivative term, sometimes called the gamma term, reflects the convexity of the option's payoff and how volatility affects value. When an option has positive gamma, increased volatility raises its expected payoff under risk-neutral pricing. The first derivative term represents the delta exposure, capturing how the option value responds to small changes in the underlying price. Finally, the discounting term ensures that the option value is expressed in present-value terms, consistent with no-arbitrage pricing. This equation, combined with appropriate boundary conditions and a terminal condition at expiration, uniquely determines the option value.

Working directly with the Black-Scholes PDE can be cumbersome because the coefficient $\frac{1}{2}\sigma^2 S^2$ in front of the second derivative varies with $S$. This spatial dependence of the coefficients means that the diffusion rate changes across the price grid, which can complicate numerical analysis and introduce additional discretization challenges. A common transformation simplifies the problem by converting to log-price coordinates. The key insight is that if the stock price follows geometric Brownian motion, then the log of the stock price follows arithmetic Brownian motion with constant volatility. Let $x = \ln(S)$, so $S = e^x$. After applying the chain rule to transform the derivatives:

$$\begin{aligned} \frac{\partial V}{\partial S} &= \frac{\partial V}{\partial x}\frac{\partial x}{\partial S} && \text{(chain rule)} \\ &= \frac{\partial V}{\partial x}\frac{1}{S} && \text{(derivative of } x = \ln S \text{)} \\ &= \frac{1}{S}\frac{\partial V}{\partial x} && \text{(rearrange)} \end{aligned}$$

where:

• $V$: option value
• $S$: asset price
• $x$: log-price coordinate ($x = \ln S$)

The transformation of the first derivative follows directly from the chain rule. For the second derivative, we must apply the product rule and chain rule together, recognizing that $\frac{1}{S}$ itself depends on $x$. The calculation yields:

$$\frac{\partial^2 V}{\partial S^2} = \frac{1}{S^2}\left(\frac{\partial^2 V}{\partial x^2} - \frac{\partial V}{\partial x}\right)$$

where:

• $\frac{\partial^2 V}{\partial S^2}$: second derivative (gamma) in price coordinates
• $\frac{\partial^2 V}{\partial x^2}, \frac{\partial V}{\partial x}$: derivatives in log-price coordinates

Notice how the transformation introduces a correction term: the second derivative in log-coordinates does not simply equal $S^2$ times the second derivative in price coordinates. This correction arises because the log transformation is nonlinear, stretching the coordinate system differently at different price levels. Substituting these transformed derivatives into the Black-Scholes PDE and simplifying yields a much simpler equation:

$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 \frac{\partial^2 V}{\partial x^2} + \left(r - \frac{1}{2}\sigma^2\right)\frac{\partial V}{\partial x} - rV = 0$$

where:

• $V$: option value
• $t$: time
• $\sigma$: volatility
• $r$: risk-free rate
• $x$: log-price coordinate

This transformed equation has constant coefficients, making it more amenable to finite difference approximations because the same discretization scheme applies uniformly across the entire spatial domain. The coefficient $r - \frac{1}{2}\sigma^2$ in front of the first derivative is the risk-neutral drift of the log-price process, which should be familiar from the study of the Black-Scholes-Merton model. To keep things clear and practical, we'll focus on the original $S$-coordinate formulation. This makes setting boundary conditions and interpreting results more intuitive.

To solve the PDE numerically, we must replace the continuous price and time domains with discrete grid points. This process, called discretization, turns the continuous partial differential equation into a system of algebraic equations that a computer can solve. For the asset price $S$, we create a grid from $S_{\min} = 0$ to $S_{\max}$, where $S_{\max}$ is chosen large enough that the option value is negligible for out-of-the-money options. For time, we work backward from expiration $T$ to the current time $t = 0$, which aligns with the natural direction for solving parabolic PDEs with terminal conditions.

The grid points are defined as:

• Price grid: $S_j = j \cdot \Delta S$ for $j = 0, 1, \ldots, M$, where $\Delta S = S_{\max}/M$
• Time grid: $t_n = n \cdot \Delta t$ for $n = 0, 1, \ldots, N$, where $\Delta t = T/N$

Let $V_j^n$ denote the option value at grid point $(S_j, t_n)$. This notation uses subscripts for spatial indices and superscripts for temporal indices, a convention that helps distinguish the two dimensions of the problem. Our goal is to compute the option values at $n = 0$ (current time) given the terminal condition at $n = N$ (expiration). The backward-in-time solution process mirrors the logic of risk-neutral valuation: we know the option's payoff at expiration and work backward to determine its present value.