Finite Difference Methods for Option Pricing: PDE Solutions

Derivative pricing uses two major numerical approaches: binomial trees, which discretize the underlying asset's price movements into up and down jumps, and Monte Carlo simulation, which generates random paths of the underlying asset to estimate option values. Both methods are powerful, but they approach the pricing problem from a probabilistic perspective, simulating possible futures and working backward or averaging outcomes.

Finite difference methods take an entirely different approach. Rather than simulating the stochastic process that drives asset prices, they directly attack the partial differential equation (PDE) that governs option values. As shown in the chapter on the Black-Scholes-Merton PDE, the option price $V(S,t)$ satisfies a deterministic equation that can be solved numerically by discretizing both the asset price and time dimensions onto a grid.

This chapter develops three finite difference schemes, each representing a different trade-off between computational simplicity and numerical stability: the explicit scheme, the implicit scheme, and the Crank-Nicolson scheme. Understanding these methods provides essential tools for pricing options where analytical solutions don't exist, including American options with early exercise features and exotic derivatives with complex payoff structures. We'll examine how to set up the computational grid, handle boundary conditions, ensure numerical stability, and implement each scheme in practice.

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:

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:

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:

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:

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, transforms 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.

The choice of grid parameters involves trade-offs that you must carefully consider. Finer grids (larger $M$ and $N$) provide more accurate solutions but require more computation and memory. The grid must be fine enough to capture the option's behavior near the strike price, where values change most rapidly, but extending the grid too far into deep in-the-money or out-of-the-money regions wastes computational resources. The ratio between $\Delta t$ and $(\Delta S)^2$ is particularly important for stability, as we'll see when analyzing each scheme.

The explicit scheme is the most straightforward approach to solving the Black-Scholes PDE numerically. Its simplicity makes it an excellent starting point for understanding finite difference methods, even though its stability limitations often favor other schemes in practice. The explicit scheme approximates the PDE derivatives using values at the current time step to calculate values at the previous time step (remembering that we work backward from expiration).

The central idea of finite difference methods is to replace derivatives with algebraic approximations involving function values at nearby grid points. These approximations come from Taylor series expansions. We approximate the partial derivatives using finite differences. At grid point $(S_j, t_{n+1})$, the approximations are:

Time derivative (backward difference):

where:

• $V_j^n$: option value at grid point $j$ and current time step $n$
• $V_j^{n+1}$: option value at grid point $j$ and next time step $n+1$
• $\Delta t$: time step size

This approximation follows from the definition of a derivative as a limit. When we cannot take an infinitesimal step, we use a small but finite step $\Delta t$ instead. The error in this approximation is proportional to $\Delta t$, meaning we achieve first-order accuracy in time.

First spatial derivative (central difference):

where:

• $V_{j\pm 1}^{n+1}$: option values at adjacent spatial points at time $t_{n+1}$
• $\Delta S$: price step size

The central difference uses function values on both sides of the point of interest, which leads to a symmetric approximation with second-order accuracy. Compared to one-sided differences, the central difference cancels the leading error terms, giving a more accurate estimate of the derivative.

Second spatial derivative (central difference):

where:

• $V_{j\pm 1}^{n+1}$: option values at adjacent spatial points
• $\Delta S$: price step size

This formula for the second derivative can be derived by applying the first derivative approximation twice, or equivalently by using Taylor expansions and solving for the second derivative. The resulting expression has an intuitive interpretation: it measures the curvature of the option value function by comparing the value at a point to the average of its neighbors. Notice that all spatial derivatives use values at time $t_{n+1}$. This makes the scheme "explicit" because we can directly calculate $V_j^n$ from known values at $t_{n+1}$ without needing to solve a system of equations.

With our finite difference approximations in hand, we can now derive the explicit update formula by substituting into the Black-Scholes PDE. Substituting these approximations into the Black-Scholes PDE and rearranging for $V_j^n$:

where:

• $V_j^n$: unknown option value at current time step $n$
• $V_{\cdot}^{n+1}$: known option values at next time step $n+1$
• $S_j$: asset price at grid point $j$
• $\Delta t, \Delta S$: time and price step sizes
• $r, \sigma$: risk-free rate and volatility

Our objective is to isolate $V_j^n$ on one side of the equation, expressing it in terms of known quantities. Multiplying by $\Delta t$ and rearranging to isolate $V_j^n$:

where:

• $V_j^n$: option value at the current time step
• $\Delta t$: time step size
• Terms in parenthesis: discrete approximation of the Black-Scholes operator

The expression in parentheses represents the spatial part of the Black-Scholes operator, evaluated at the known time step. To make the formula more useful for computation, we collect terms by which grid point they involve. Grouping terms by the grid points $V_{j-1}^{n+1}$, $V_j^{n+1}$, and $V_{j+1}^{n+1}$:

where:

• $V_j^n$: option value at the current time step
• $\Delta t$: time step size
• $\Delta S$: price step size
• $\sigma$: volatility
• $r$: risk-free rate
• $S_j$: asset price at grid point $j$
• $V_{\cdot}^{n+1}$: option values at the next time step

The structure of this equation reveals that the option value at any interior grid point is a weighted combination of three values at the next time step: the value at the same price level and the values at adjacent price levels. Using $S_j = j\Delta S$, these bracketed terms simplify to a more compact form:

where:

• $V_j^n$: option value at the current time step $n$
• $V_{j-1}^{n+1}, V_j^{n+1}, V_{j+1}^{n+1}$: known option values at the next time step $n+1$
• $a_j, b_j, c_j$: weighting coefficients for the grid points

The coefficients are:

where:

• $\Delta t$: time step size
• $\sigma$: volatility
• $r$: risk-free rate
• $j$: grid index ($j = S_j/\Delta S$)

The option value at each interior point depends only on three values at the next time step, which we already know from either the terminal condition or previous calculations. Intuitively, this weighted sum resembles an expected value calculation: $a_j$ and $c_j$ represent the weights (pseudo-probabilities) of the asset price moving down or up, while $b_j$ is the weight of staying effectively in the same region, all adjusted for the risk-free discounting. This probabilistic interpretation connects finite difference methods to the binomial tree approach, where we also computed option values as discounted expected values under risk-neutral probabilities. The explicit scheme essentially implements the same logic on a finer grid.

The explicit scheme proceeds as follows:

1. Initialize the grid with the terminal payoff at expiration
2. Apply boundary conditions at $S = 0$ and $S = S_{\max}$
3. Sweep backward through time, calculating each $V_j^n$ from values at $t_{n+1}$

The calculated option price represents the fair value derived from the explicit finite difference grid. Since the grid is relatively fine (5000 time steps), this value closely approximates the theoretical Black-Scholes price, serving as a baseline for verifying the numerical method's accuracy.

The explicit scheme has a critical limitation. It becomes unstable if the time step is too large. This instability arises because the coefficients $a_j, b_j, c_j$ must behave like probabilities to ensure the solution remains bounded. When we interpret these coefficients as weights in an expected value calculation, they should all be non-negative for the scheme to behave sensibly. If $\Delta t$ is too large, $b_j$ becomes negative, causing errors to oscillate and grow exponentially with each time step. Small numerical errors introduced by rounding or discretization get amplified rather than damped, eventually overwhelming the true solution. The stability condition requires:

where:

• $\Delta t$: maximum allowable time step
• $\sigma$: volatility
• $r$: risk-free rate
• $j_{\max}$: maximum grid index ($M$)

In practice, this means:

where:

• $\Delta t$: maximum stable time step
• $\Delta S$: price grid spacing
• $\sigma$: volatility
• $S_{\max}$: maximum price on the grid

This constraint can be severe because the bound depends on the square of the price grid spacing. Halving $\Delta S$ to improve spatial resolution requires dividing $\Delta t$ by four to maintain stability. For typical parameters, the explicit scheme may require thousands of time steps to remain stable, making it computationally expensive despite its conceptual simplicity. This stability restriction motivates the development of implicit methods that remove or relax this constraint.

The stability requirement explains why we used 5,000 time steps in our implementation. Using fewer steps would cause the solution to oscillate wildly and diverge from the true option value.

The implicit scheme overcomes the stability limitation by evaluating the spatial derivatives at the current time step $t_n$ rather than the future time step $t_{n+1}$. This small change significantly improves numerical stability. Rather than computing values at time $t_n$ directly from known values at $t_{n+1}$, the implicit scheme requires solving for all values at $t_n$ simultaneously, since the update formula for each point involves unknown values at neighboring points at the same time level.

The implicit scheme uses the same time derivative approximation as the explicit scheme but evaluates the spatial derivatives at a different time level. The implicit scheme uses:

Time derivative (backward difference):

where:

• $V_j^n$: option value at grid point $j$ and current time step $n$
• $V_j^{n+1}$: option value at grid point $j$ and next time step $n+1$
• $\Delta t$: time step size

Spatial derivatives at time $t_n$:

where:

• $V_{j\cdot}^n$: option values at current time $t_n$ (unknowns)
• $\Delta S$: price step size

The crucial distinction from the explicit scheme is that these spatial derivatives use values at time $t_n$, which are the unknowns we seek to determine. This creates a coupling between neighboring grid points at the same time level: to find $V_j^n$, we need $V_{j-1}^n$ and $V_{j+1}^n$, but these values are also unknown.

Because all values at time $t_n$ are coupled together, we cannot compute them one at a time. Instead, we must solve a system of linear equations simultaneously. Substituting the approximations into the Black-Scholes PDE:

where:

• $V_{\cdot}^n$: unknown option values at current time $t_n$
• $V_j^{n+1}$: known option value at next time $t_{n+1}$

We multiply by $\Delta t$ and isolate $V_j^{n+1}$:

Grouping the unknown terms at $t_n$ by spatial index:

Defining coefficients $\alpha_j, \beta_j, \gamma_j$ to simplify notation, we obtain the linear system:

where:

• $V_{\cdot}^n$: unknown values at current time $t_n$
• $V_j^{n+1}$: known value at next time $t_{n+1}$
• $\alpha_j, \beta_j, \gamma_j$: coefficients for the implicit system

The coefficients are:

where:

• $\Delta t$: time step
• $\sigma$: volatility
• $r$: risk-free rate
• $j$: spatial index corresponding to $S_j$

Unlike the explicit scheme, we cannot directly calculate $V_j^n$ from $V_j^{n+1}$. Instead, we must solve a system of linear equations at each time step. For the $M-1$ interior points, this forms a tridiagonal matrix equation:

where:

• $\mathbf{A}$: tridiagonal coefficient matrix
• $\mathbf{V}^n$: vector of option values at time $t_n$ (unknowns)
• $\mathbf{V}^{n+1}$: vector of known option values at time $t_{n+1}$
• $\mathbf{b}$: vector incorporating boundary conditions

The matrix $\mathbf{A}$ has the form:

where:

• $\alpha_j, \beta_j, \gamma_j$: coefficients defined previously
• $M$: number of spatial grid intervals

The tridiagonal structure arises because each interior equation involves only three unknowns: the value at the current grid point and its two immediate neighbors. This sparse structure is extremely important computationally because it allows efficient solution algorithms. By coupling all $V_j^n$ values in a single system, the implicit method allows information to propagate across the entire grid instantaneously, mirroring the infinite propagation speed of the diffusion equation. In contrast, the explicit scheme propagates information only one grid point per time step, which is why it requires the stability constraint. This global dependency ensures errors are damped rather than amplified, providing unconditional stability. Tridiagonal systems can be solved very efficiently using the Thomas algorithm (a specialized form of Gaussian elimination), requiring only $O(M)$ operations per time step, making the implicit scheme practical even though it requires solving a linear system.

The implicit method yields a price very close to the explicit method but achieves this using only 500 time steps, one-tenth of the effort required for the explicit scheme. This efficiency demonstrates the practical benefit of unconditional stability.

The implicit scheme is unconditionally stable, meaning it produces bounded solutions for any choice of $\Delta t$ and $\Delta S$. This property allows you to use much larger time steps, dramatically reducing computation. Although each time step requires solving a linear system, the overall efficiency often surpasses the explicit scheme because far fewer time steps are needed.

The stability comes from the implicit treatment of the spatial derivatives. Mathematically, when we analyze the amplification factor of the numerical scheme, the implicit treatment ensures that errors decay rather than grow over time. Any errors introduced in one time step are damped rather than amplified in subsequent steps. The price we pay for this stability is the need to solve a linear system at each time step, but for tridiagonal systems this cost is minimal.

The Crank-Nicolson scheme represents a middle ground between explicit and implicit methods. It averages the spatial derivatives at times $t_n$ and $t_{n+1}$, achieving second-order accuracy in time while maintaining unconditional stability. This combination of accuracy and stability makes Crank-Nicolson the method of choice for most practical applications.

Instead of evaluating spatial derivatives entirely at $t_n$ (implicit) or $t_{n+1}$ (explicit), Crank-Nicolson uses the average of both evaluations. This averaging can be motivated by the trapezoidal rule for numerical integration: when integrating the time derivative, we approximate the integral of the spatial operator by averaging its values at the beginning and end of the time interval. The scheme uses:

where:

• $V_{\cdot}^n$: values at current time step
• $V_{\cdot}^{n+1}$: values at next time step
• $\Delta S$: price step size

and similarly for the first derivative. This averaging can be viewed as taking a half-step with the explicit scheme and a half-step with the implicit scheme. The explicit part uses known information from the future time step, while the implicit part ensures stability by coupling values at the current time step.

Substituting the approximations into the Black-Scholes PDE:

where:

• $\mathcal{L}_j V$: spatial operator at node $j$, equal to $\frac{1}{2}\sigma^2 S_j^2 \frac{\partial^2 V}{\partial S^2} + rS_j\frac{\partial V}{\partial S} - rV$
• $V^n, V^{n+1}$: option values at current and next time steps

The spatial operator $\mathcal{L}_j$ encapsulates all the terms involving spatial derivatives and the discounting term. Applying this operator at both time levels and averaging creates the characteristic Crank-Nicolson structure. Rearranging terms by time step:

where:

• LHS: implicit terms at current time $t_n$
• RHS: explicit terms at next time $t_{n+1}$
• $\alpha_j, \beta_j, \gamma_j$: Crank-Nicolson coefficients

The left side involves unknown values at $t_n$, forming a tridiagonal system just as in the implicit scheme. The right side uses known values at $t_{n+1}$, which we can compute explicitly before solving the system. The coefficients are:

where:

• $\Delta t$: time step
• $\sigma$: volatility
• $r$: risk-free rate
• $j$: spatial index

Notice that these coefficients are exactly half of the corresponding implicit scheme coefficients, reflecting the averaging of the spatial operator between two time levels. The left side involves unknown values at $t_n$, while the right side uses known values at $t_{n+1}$. This again leads to a tridiagonal system at each time step.

The Crank-Nicolson price is consistent with the other methods. By averaging the explicit and implicit approaches, it maintains stability while improving accuracy, effectively reducing the bias introduced by time discretization.

The Crank-Nicolson scheme achieves $O((\Delta t)^2)$ accuracy in time, compared to $O(\Delta t)$ for both the explicit and implicit schemes. This means that halving the time step reduces the time discretization error by a factor of four rather than two. The improved accuracy comes from the averaging, which cancels the leading-order error terms in the Taylor expansion. For most applications, Crank-Nicolson provides the best balance of accuracy, stability, and computational efficiency.

Proper boundary conditions are essential for finite difference methods. They determine the option value at the edges of the computational grid, where the standard update formulas cannot be applied directly because they would require values outside the grid.

At expiration ($t = T$), the option value equals its payoff. This is the fundamental condition that initializes our backward-in-time computation:

where:

• $V(S,T)$: option value at expiration
• $S$: stock price at expiration
• $K$: strike price