Two-Factor Finite Differences: When One Dimension Is Not Enough

Chapter 79 is about what happens when your problem has two random factors. Convertible bonds with stochastic interest rates. Exotic options with stochastic volatility. Barrier options where you model both the stock and the rates. These are real problems that real desks face, and they need two-factor numerical methods.

The jump from one factor to two is not trivial. Your grid goes from two dimensions (S and t) to three dimensions (S, r, and t). Memory requirements multiply. The stability constraint for the explicit method gets worse. And the implicit methods require new tricks.

But it is doable. And Chapter 79 gives you three approaches: the straightforward explicit method (easy but slow), ADI (clever and fast), and Hopscotch (clever in a different way).

The Two-Factor PDE

The general two-factor equation looks like:

$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma_S^2 S^2 \frac{\partial^2 V}{\partial S^2} + \rho \sigma_S \sigma_r S \frac{\partial^2 V}{\partial S \partial r} + \frac{1}{2}\sigma_r^2 \frac{\partial^2 V}{\partial r^2} + \text{drift terms} - rV = 0$$

Wilmott writes this concretely with S (asset) and r (interest rate) because abstract notation is harder to follow. But everything applies to any pair of variables.

The key new term is the cross derivative, the mixed partial with respect to both S and r. This term comes from the correlation between the two random factors. When correlation is zero, the cross derivative vanishes and life is simpler. When it is not zero, you have to deal with it, and it makes some methods harder.

The Three-Dimensional Grid

The grid now has three axes: S indexed by i (from 0 to I), r indexed by j (from 0 to J), and time indexed by k. The option value is V_ij^k. At expiry, you fill in the payoff for all combinations of S and r. Then you march backward in time.

Boundary conditions must be specified on all four edges of the S-r rectangle at each time step: S = 0, S = S_max, r = 0, and r = r_max. These depend on the contract. For a convertible bond, you might have the bond value approaching zero for large r and large S, the conversion value for large S, and specific behavior at S = 0.

The Explicit Method in Two Dimensions

The explicit method extends naturally. At each point (i, j), you approximate all the derivatives using values at time step k and calculate the value at time step k+1. The stencil now uses nine points (the center, four direct neighbors, and four diagonal neighbors).

The cross derivative is approximated as:

$$\frac{\partial^2 V}{\partial S \partial r} \approx \frac{V_{i+1,j+1} - V_{i-1,j+1} - V_{i+1,j-1} + V_{i-1,j-1}}{4 \delta S \delta r}$$

This uses the four corner points and has the nice property that it is accurate to O(delta_S^2, delta_r^2).

Wilmott shows a VB code fragment for a convertible bond. Points worth noting:

• Central differences for the S derivative (delta).
• Upwind differencing for the r derivative. Interest rate volatility is usually small, making the problem nearly hyperbolic in the r direction. Central differences can be unstable there.
• The time stepping loop visits all interior (i, j) points.
• After time stepping, the conversion constraint is applied: V = max(V, conversion_rate * S).

The code is straightforward. You could implement it in an afternoon if you already understand the one-factor version.

Stability: Worse Than Before

The bad news: the stability constraint is even more severe than in one dimension. For a pure diffusion problem with no correlation, the constraint becomes:

$$\delta t \left( \frac{a}{\delta S^2} + \frac{d}{\delta r^2} \right) \le \frac{1}{2}$$

where a and d are the diffusion coefficients. This is the sum of two terms, each of which would individually constrain the time step. Together, they make delta_t even smaller.

With correlation, things get messier. The cross derivative term can cause additional instability issues.

Computation Time

Wilmott derives the scaling for the explicit method in d dimensions. If you want accuracy epsilon:

• Asset step: delta_S proportional to sqrt(epsilon) (since error is O(delta_S^2))
• Time step: delta_t proportional to epsilon (since error is O(delta_t) and stability links delta_t to delta_S^2)
• Number of grid points in each direction: proportional to 1/sqrt(epsilon)
• Total computation time for one option:

$$\text{Time} \sim \epsilon^{-(1 + d/2)}$$

For d=1: Time scales as epsilon^(-3/2). For d=2: Time scales as epsilon^(-2). The jump to two factors makes things significantly more expensive for the same accuracy.

For a portfolio of M options: Time scales as M * epsilon^(-(1 + d/2)). But you get all the basic Greeks for free (delta, gamma with respect to both variables). Vega requires a separate run.

Early exercise barely adds to the computation. It is just a comparison and replacement at each grid point, same as in one dimension.

ADI: Alternating Direction Implicit

The explicit method is simple but slow. A full two-factor Crank-Nicolson would be fast but requires solving a large system where every (i,j) point is linked to every other. The matrix is not tridiagonal. It has a banded structure that is expensive to solve.

ADI splits the time step into two halves and solves each half with implicitness in only one direction.

First half-step (k to k+1/2): Explicit in S, implicit in r. For each fixed i, you solve a tridiagonal system in j. This is no harder than a one-factor implicit solve, you just have I separate systems to solve.

Second half-step (k+1/2 to k+1): Implicit in S, explicit in r. For each fixed j, you solve a tridiagonal system in i.

Each half-step alternates which direction is treated implicitly. Hence the name: Alternating Direction Implicit.

The result:

• Stable for all time steps (no stability constraint)
• Error: O(delta_t^2, delta_S^2, delta_r^2) (same as Crank-Nicolson)
• Each half-step is a set of one-dimensional tridiagonal solves (fast and simple)

The main drawback: the cross derivative term does not decompose cleanly into S-only and r-only terms. There are workarounds but they add complexity. When correlation is zero, ADI is perfect. When correlation is significant, you need to be more careful with the implementation.

Hopscotch: A Different Approach

The Hopscotch method is named after the children’s game because of how it jumps around the grid. The idea is creative.

Color the grid like a checkerboard. Points where i+j is odd get circles. Points where i+j is even get squares.

First sweep (odd time steps): Calculate all circle points explicitly using the standard formula with values from time step k. Then, in a second sweep, calculate all square points. But here is the trick: the square points use the just-calculated k+1 values at the circle points (which are their neighbors). So the formula looks implicit, but because the circle values are already known, no system of equations needs to be solved.

Second sweep (even time steps): Reverse the roles. Squares go first (explicit), then circles use the updated values.

The result:

• Stable for all time steps (despite looking partly explicit)
• Error: O(delta_t, delta_S^2, delta_r^2) (first-order in time, so less accurate than ADI)
• No simultaneous equations to solve (everything is computed pointwise)

The catch: when there is a cross derivative term, the neat explicitness breaks down and you need to solve small systems. For uncorrelated problems, Hopscotch is great. For correlated problems, ADI is usually better.

When to Use What

Here is the practical decision tree for two-factor problems:

Prototyping or simple problems: Use the explicit method. Easy to code, easy to debug. Accept the slow computation.

Production with no correlation: ADI or Hopscotch, both work well. ADI is more accurate (second-order in time). Hopscotch is simpler to code.

Production with correlation: ADI with appropriate handling of the cross derivative. More work to implement but worth it for the stability and accuracy.

Three or more factors: Finite differences start to struggle. Consider Monte Carlo instead. The explicit method can technically handle any number of dimensions, but the stability constraint and memory requirements become prohibitive.

Wilmott notes that simple explicit finite differences are similar to binomial trees. So the methods in this chapter can never do worse than the binomial method. ADI and Hopscotch can do much better.

Practical Example: Convertible Bonds

Convertible bonds are the poster child for two-factor methods. You have a stock that is random (equity risk) and an interest rate that is random (rate risk). The bond can be converted into stock at the holder’s option, adding an early-exercise feature.

The two-factor explicit method handles all of this:

• The PDE has terms for equity diffusion, rate diffusion, correlation, and both drifts
• The interest rate model can be yield-curve fitted (time-dependent drift)
• The conversion constraint is applied at each time step: V = max(V, conversion_rate * S)
• Boundary conditions handle what happens at S = 0, large S, r = 0, and large r

The code fragment in the chapter is maybe 30 lines of VB, including the conversion constraint. It is not elegant, but it works and it is readable.

The Takeaway

Two-factor finite differences are a natural extension of the one-factor methods. The explicit method extends with minimal effort but the stability constraint gets worse. ADI cleverly splits each time step into two half-steps, each of which is just a bunch of one-dimensional solves. Hopscotch avoids simultaneous equations entirely by alternating explicit and implicit treatment across the grid.

The main enemy is the cross derivative term from correlation. When correlation is zero, everything is clean. When it is not, ADI needs modifications and Hopscotch loses its simplicity.

For three or more factors, Wilmott recommends switching to Monte Carlo. Finite differences scale badly in high dimensions. But for two factors, which covers a huge range of practical problems (convertible bonds, stochastic volatility options, two-asset products), the methods in this chapter are the right tool.

Next up: Monte Carlo simulation, the method that works when dimensions get too high for grids.

Previous post: Advanced Finite Differences: Implicit, Crank-Nicolson, and More

Next post: Monte Carlo Simulation: The Random Path to Option Prices