Finite Difference Methods for Option Pricing › The Implicit Finite Difference Scheme

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

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

Spatial derivatives at time $t_n$:

$$\begin{aligned} \frac{\partial V}{\partial S} &\approx \frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} \\[1em] \frac{\partial^2 V}{\partial S^2} &\approx \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{(\Delta S)^2} \end{aligned}$$

where:

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

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

$$\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 - 2V_j^n + V_{j-1}^n}{(\Delta S)^2} \\ &+ rS_j\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} - rV_j^n = 0 \end{aligned}$$

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

$$V_j^{n+1} = V_j^n - \Delta t \left[ \frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1}^n - 2V_j^n + V_{j-1}^n}{(\Delta S)^2} + rS_j\frac{V_{j+1}^n - V_{j-1}^n}{2\Delta S} - rV_j^n \right]$$

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

$$\begin{aligned} V_j^{n+1} &= \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 \\ &\quad + \left[ 1 + \frac{\Delta t}{(\Delta S)^2}\sigma^2 S_j^2 + \Delta t r \right] V_j^n \\ &\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 \end{aligned}$$

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

$$-\alpha_j V_{j-1}^n + (1 + \beta_j) V_j^n - \gamma_j V_{j+1}^n = V_j^{n+1}$$

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:

$$\begin{aligned} \alpha_j &= \frac{\Delta t}{2}\left(\sigma^2 j^2 - rj\right) \\ \beta_j &= \Delta t\left(\sigma^2 j^2 + r\right) \\ \gamma_j &= \frac{\Delta t}{2}\left(\sigma^2 j^2 + rj\right) \end{aligned}$$

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:

$$\mathbf{A}\mathbf{V}^n = \mathbf{V}^{n+1} + \mathbf{b}$$

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:

$$\mathbf{A} = \begin{pmatrix} 1 + \beta_1 & -\gamma_1 & 0 & \cdots & 0 \\ -\alpha_2 & 1 + \beta_2 & -\gamma_2 & \cdots & 0 \\ 0 & -\alpha_3 & 1 + \beta_3 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & -\alpha_{M-1} & 1 + \beta_{M-1} \end{pmatrix}$$

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, and gives 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.