Finite Difference Methods for Option Pricing › The Explicit Finite Difference 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):

$$\frac{\partial V}{\partial t} \approx \frac{V_j^{n+1} - V_j^n}{\Delta t}$$

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):

$$\frac{\partial V}{\partial S} \approx \frac{V_{j+1}^{n+1} - V_{j-1}^{n+1}}{2\Delta S}$$

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):

$$\frac{\partial^2 V}{\partial S^2} \approx \frac{V_{j+1}^{n+1} - 2V_j^{n+1} + V_{j-1}^{n+1}}{(\Delta S)^2}$$

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$:

$$\begin{aligned} \frac{V_j^{n+1} - V_j^n}{\Delta t} &+ \frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1}^{n+1} - 2V_j^{n+1} + V_{j-1}^{n+1}}{(\Delta S)^2} \\ &+ rS_j\frac{V_{j+1}^{n+1} - V_{j-1}^{n+1}}{2\Delta S} - rV_j^{n+1} = 0 \end{aligned}$$

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$:

$$\begin{aligned} V_j^n = V_j^{n+1} + \Delta t \bigg( &\frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1}^{n+1} - 2V_j^{n+1} + V_{j-1}^{n+1}}{(\Delta S)^2} \\ &+ rS_j\frac{V_{j+1}^{n+1} - V_{j-1}^{n+1}}{2\Delta S} - rV_j^{n+1} \bigg) \end{aligned}$$

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}$:

$$\begin{aligned} V_j^n &= \left[ \frac{\Delta t}{2(\Delta S)^2}\sigma^2 S_j^2 - \frac{\Delta t}{2\Delta S}rS_j \right] V_{j-1}^{n+1} \\ &\quad + \left[ 1 - \frac{\Delta t}{(\Delta S)^2}\sigma^2 S_j^2 - \Delta t r \right] V_j^{n+1} \\ &\quad + \left[ \frac{\Delta t}{2(\Delta S)^2}\sigma^2 S_j^2 + \frac{\Delta t}{2\Delta S}rS_j \right] V_{j+1}^{n+1} \end{aligned}$$

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:

$$V_j^n = a_j V_{j-1}^{n+1} + b_j V_j^{n+1} + c_j V_{j+1}^{n+1}$$

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:

$$\begin{aligned} a_j &= \frac{\Delta t}{2}\left(\sigma^2 j^2 - rj\right) \\ b_j &= 1 - \Delta t\left(\sigma^2 j^2 + r\right) \\ c_j &= \frac{\Delta t}{2}\left(\sigma^2 j^2 + rj\right) \end{aligned}$$

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:

$$\Delta t \leq \frac{1}{\sigma^2 j_{\max}^2 + r}$$

where:

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

In practice, this means:

$$\Delta t \leq \frac{(\Delta S)^2}{\sigma^2 S_{\max}^2}$$

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.