Finite Difference Methods

Quants who implement valuation models use the finite difference method when they need an option price on a grid, and risk managers use it when they check the calculated prices. It also comes up in the valuation of options that can be exercised before expiry.

Replacing a partial differential equation with differences on a grid

Under the Black-Scholes model, the option price satisfies a partial differential equation. The finite difference method places grid points over the range of prices of interest, replaces the derivatives in the equation with differences between neighboring points, and solves for the price at each point numerically¹. The derivation of the equation is covered in Introduction to Stochastic Calculus (Ito's Lemma).

Let $V$ be the price of a European option, $S$ the price of the underlying asset, and $\tau = T - t$ the time to expiry. The interest rate $r$ and the volatility $\sigma$ are constant, and the underlying pays no dividends. The equation is as follows².

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

Because time is measured as the time to expiry $\tau$, expiry is $\tau = 0$ and today is $\tau = T$. The payoff at expiry is the starting value (the initial condition), and the calculation advances by increasing $\tau$ until it reaches today's price. Changing the variable to $x = \ln(S/K)$ gives the following equation, whose coefficients do not depend on $S$³. $K$ is the strike price.

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

Place a grid with spacing $h$ in the $x$ direction and $\Delta\tau$ in the $\tau$ direction. Let $V_j^n$ be the approximate price at the point $x_j = x_0 + jh$ and time $\tau_n = n\Delta\tau$. The three derivatives are replaced by the following differences⁴.

The order of the whole calculation is as follows. The boundary values in step 2 are explained in the section "Boundary conditions fix the values outside the grid", and the comparison with the exercise value in step 4 in the section "An American option never falls below its exercise value".

1. Put the payoff at expiry on every point of the grid.
2. Put the boundary values at the lower and upper ends (and reset them to the values for that time at each step forward).
3. Advance one time step and find the value at each point.
4. In the calculation on this page, only for an American option, replace each point's value with the larger of that value and the exercise value (this is not a procedure given in the sources).
5. Repeat steps 3 and 4 from expiry to today, then read the value at the point that corresponds to the current price of the underlying as the price.

Three solution methods differ in how time is stepped

The symbols $\theta$, $\mu$, $a$, $b$, and $c$ below are used only on this page. They mean something different from $\theta(t)$ and the mean reversion speed $a$ in the Hull-White model, and from the drift $\mu$ of a stochastic process.

Collect the differences in the price direction at each point $j$ into one operation, written as follows.

$$(\mathcal{L}V)_j = a\,V_{j-1} + b\,V_j + c\,V_{j+1}$$

With $D = \sigma^2/2$ (the coefficient of the second derivative) and $\mu = r - \sigma^2/2$ (the coefficient of the first derivative), the coefficients are $a = D/h^2 - \mu/(2h)$, $b = -2D/h^2 - r$, and $c = D/h^2 + \mu/(2h)$. The parameter $\theta$ chooses whether $\mathcal{L}V$ is evaluated at the earlier time $n$, at the next time $n+1$, or as an average of the two⁵.

$$V_j^{n+1} - \theta\,\Delta\tau\,(\mathcal{L}V^{n+1})_j = V_j^n + (1 - \theta)\,\Delta\tau\,(\mathcal{L}V^n)_j$$

$\theta = 0$ is the explicit method, $\theta = 1$ is the implicit method, and $\theta = 1/2$ is the Crank-Nicolson method⁵. The three differ as follows.

The explicit method avoids solving a system of equations. However, for a diffusion equation it becomes unstable when the time step exceeds $(\Delta x)^2/(2D)$⁶. In the log-price form of the Black-Scholes equation, the coefficient of the second derivative is $\sigma^2/2$, so a guide that looks at the diffusion term only is $\Delta\tau \le h^2/\sigma^2$⁷. Halving the grid spacing $h$ multiplies the number of time steps by 4 and the number of points by 2, so the computation time grows about 8 times⁶.

The implicit and Crank-Nicolson methods solve a tridiagonal system of equations for the values at the next time. A tridiagonal system can be solved in time and memory proportional to the number of points⁸. The implicit method has no time step limit and is unconditionally stable⁹. For smooth initial values, the error of the Crank-Nicolson method is $O(h^2, \Delta\tau^2)$, which shrinks faster than the $O(h^2, \Delta\tau)$ of the other values of $\theta$. For $\theta \ge 1/2$, the $L^2$ norm of the solution does not grow whatever the time step (unconditional stability in the $L^2$ norm), so the time step can be taken in proportion to $h$ and still give both accuracy and efficiency¹⁰.

The explicit formula expresses the value at the next time as a weighted sum of three values at the earlier time.

$$V_j^{n+1} = a\Delta\tau\,V_{j-1}^n + (1 + b\Delta\tau)\,V_j^n + c\Delta\tau\,V_{j+1}^n$$

For a diffusion equation without discounting or advection, the three coefficients lie between 0 and 1 and sum to 1 under the stability condition. They can therefore be read as probabilities, and the explicit method takes the same form as a calculation walking a trinomial tree backward from expiry¹¹. In the formula above, the discounting term makes the coefficients sum to 1 minus $r\Delta\tau$, which is the probabilities multiplied by one period of discounting.

Deep diveHow to check the stability of the explicit and implicit methods

Split the solution of the heat equation $\partial V/\partial\tau = D\,\partial^2 V/\partial x^2$ into components of the wave form $e^{ipx}$. Suppose each component is multiplied by $G$ in one time step. If $|G| \le 1$, the component does not grow, and if $|G| > 1$, it keeps growing. With $\lambda = D\Delta\tau/h^2$, $G$ is as follows for the explicit method ($G_{\mathrm{E}}$) and the implicit method ($G_{\mathrm{I}}$)¹².

$$G_{\mathrm{E}} = 1 - 4\lambda\sin^2\left(\frac{ph}{2}\right)$$ $$G_{\mathrm{I}} = \frac{1}{1 + 4\lambda\sin^2(ph/2)}$$

Because $\sin^2(ph/2)$ can come arbitrarily close to 1, the condition for $|G_{\mathrm{E}}| \le 1$ to hold for every component is $\lambda \le 1/2$. This corresponds to $\Delta\tau \le h^2/(2D)$. For the implicit method, $G_{\mathrm{I}}$ is greater than 0 and at most 1 for any positive $\lambda$, so there is no limit on the time step.

Boundary conditions fix the values outside the grid

Only a finite range can be computed, so boundaries are placed at the lower and upper ends of $x$ and the values there are fixed. For a European option, the boundary values follow from put-call parity. A call approaches $S - Ke^{-r\tau}$ where the underlying price is high and approaches 0 where it is low. A put approaches $Ke^{-r\tau} - S$ at the low end and 0 at the high end¹³. A condition that fixes the value at the boundary in this way is called a Dirichlet boundary condition¹³.

The value at expiry ($\tau = 0$) is the payoff: $\max(S - K, 0)$ for a call option and $\max(K - S, 0)$ for a put option¹⁴. In the common example (next section), the range of $x$ is $-1$ to $1$. In terms of $S$, that is 36.79 to 271.83, which is 5 times the standard deviation of the log price to expiry, $\sigma\sqrt{T} = 0.2$. In the table below, the Crank-Nicolson price with $h = 0.02$ and 200 time steps stays at 8.4235 even when this range is widened to $\pm 1.5$. Narrowing it to $\pm 0.5$ gives 8.4234, a change of 1 in the fourth decimal place¹⁵.

Solving the common example call with the three methods

Consider a fictional example of a call option with underlying price 100, strike price 100, one year to expiry, interest rate 1%, and volatility 20%. The price from the Black-Scholes formula is 8.4333¹⁶.

First, follow one step of the explicit method by hand. Let the grid spacing be $h = 0.1$ and the time step $\Delta\tau = 0.025$ (the year divided into 40 steps). Then $D = 0.2^2/2 = 0.02$ and $\mu = 0.01 - 0.02 = -0.01$. The coefficients of the formula take the following values.

All three weights are 0 or greater, and they sum to 0.99975 ($= 1 - 0.01 \times 0.025$). At the point $x_j = 0.1$ (underlying price 110.5171), advance one time step from the payoff at expiry.

The weighted sum is $0.05125 \times 0 + 0.89975 \times 10.5171 + 0.04875 \times 22.1403 = 0 + 9.4628 + 1.0793 = 10.5421$. At this point, the price 0.025 years before expiry equals the underlying price of 110.5171 minus the present value of the strike, $100 \times e^{-0.01 \times 0.025} = 99.9750$, which is 10.5421¹⁷.

Next, spread the grid over the whole range and solve with the three methods. The table shows the price at the point $S = 100$ ($x = 0$), with the difference from the Black-Scholes price. The grid spacing is $h = 0.02$ (100 intervals). The guide for the explicit method, looking at the diffusion term only, is $\Delta\tau \le h^2/\sigma^2 = 0.01$, that is, 100 time steps or more.

With a time step of 0.025, the time step is 2.5 times the guide limit, and the explicit calculation diverges. Cutting the time step to one fifth hardly changes the Crank-Nicolson difference, so this difference comes mainly from the grid spacing $h$¹⁸. For the implicit method, the difference shrinks as the time step gets smaller.

Next, refine the grid with the Crank-Nicolson method. The time step is $\Delta\tau = h/2$, so halving $h$ also halves the time step.

Each time $h$ is halved, the difference shrinks to about a quarter, which matches the error order $O(h^2, \Delta\tau^2)$ in the table above¹⁹. Other valuations will not necessarily show differences of the same size.

An American option never falls below its exercise value

An American option can also be exercised before expiry. Its value at any time is the larger of the value of holding it and the value of exercising now. The equation for an American put adds to the Black-Scholes equation the condition that the value does not fall below the exercise value $K - S$²⁰. Unlike the European case, the boundary conditions cannot be fixed from the start, and the problem becomes a free boundary problem²¹.

The sources give no procedure for imposing this on a grid. In the calculation on this page, after solving each time step, the value at each point is replaced by the larger of that value and the exercise value (intrinsic value)²². Another approach is the penalty method. It adds a nonlinear term that can be ignored where the option is not exercised and that makes the solution almost equal to the exercise value elsewhere, and solves with Newton's method (a way to solve a nonlinear equation by repeatedly solving linear ones)²³.

The next table shows the put in the common example, solved with the Crank-Nicolson method and the replacement at each time step. The European price is the value solved on the same grid and time steps without the replacement. The time steps are $\Delta\tau = 0.005, 0.0025, 0.0025$²².

The European put from the Black-Scholes formula is 7.4383. The American put is higher than the European put by the value of being able to exercise before expiry. The American put valued on a binomial tree (2,000 steps) is 7.5127, and the finite difference price approaches this value as the grid is refined²².

Choosing between binomial models, the Monte Carlo method, and finite differences

Valuations with few risk factors and early exercise, such as American options, are handled readily by finite differences, while valuations with many underlying assets are described as more efficient with the Monte Carlo method. Both the binomial model and the finite difference method are grid calculations that find the price by working backward from expiry. The binomial model uses a binomial tree. The finite difference method, as seen in the section on stepping through time, can take a form identical to walking a trinomial tree backward when the explicit method is used. The Monte Carlo method approximates the price by the average over many paths built from random numbers. The details are covered in Monte Carlo Simulation. The differences from the finite difference method are as follows.

The difference in dimension shows in the number of points. Consider a fictional example of putting 100 points on each axis. One, two, and three axes give 100, 10,000, and 1,000,000 points, and $d$ axes give $100^d$ points. The curse of dimensionality also affects deterministic methods that use grids³⁰.

Assumptions and limits

The finite difference method is a way to solve an equation numerically. The price is determined by the model that the equation assumes. The Black-Scholes equation assumes that the interest rate and the volatility are constant over the whole period³¹. If this assumption fails, the solved price does not necessarily match the market price.

The numerical error depends on the choice of grid and time step. A payoff that is not smooth is another source of error. The payoff of a call has a kink at the strike, so delta and gamma computed by central differences of the grid values have larger errors than the price itself. The largest gamma error appears at the strike point²⁸. The Crank-Nicolson method is unconditionally stable only in the $L^2$ norm, and when the payoff is not smooth, delta and gamma converge poorly in terms of the maximum error²⁸.

One remedy is Rannacher time-marching, which replaces the first few time steps with half steps of the implicit method. Replacing the first two Crank-Nicolson steps with four half steps of the implicit method has been shown to make the maximum error (the largest error over all grid points, the maximum norm) of the price, delta, and gamma shrink at second order for puts, calls, and digitals³². The calculations in the text use the Crank-Nicolson method without this replacement.

When there are two or more risk factors and the equation contains mixed derivatives (cross derivatives) between different factors, keeping second-order accuracy requires the Craig-Sneyd treatment in how the cross derivative terms are differenced²⁵.

Notes

1. Page 543, "Finite-Difference Methods", of Yuh-Dauh Lyuu's lecture notes (National Taiwan University, May 17, 2006). The notes approximate the function at grid points and solve the equation numerically by substituting differences for derivatives. ↩

2. Section 2 (page 2) of M. B. Giles, "Crank-Nicolson time-marching" (Oxford-Man Institute of Quantitative Finance). The source writes the equation with the time to expiry $\tau \equiv T - t$ and the form after changing to $x \equiv \log S$. It has no dividend term, and the text includes none either. ↩

3. The source uses $x = \log S$. This page uses $x = \ln(S/K)$, which shifts the variable by $\ln K$ and leaves the form of the derivatives unchanged. ↩

4. The forward difference in time and the central differences in space are in Lyuu's notes, pages 547–548 ("Explicit Methods", equations (59)–(61)), and in Sections 1 and 2 of Giles (pages 1–2). Section 2 of Giles uses the central difference $V_{j+1} - V_{j-1}$ for the first derivative. ↩

5. Section 1 of Giles (pages 1–2). The source defines the $\theta$ scheme, with $\theta = 0$ as the explicit method, $\theta = 1$ as the implicit method, and $\theta = 1/2$ as the Crank-Nicolson method. The coefficients of the operator $\mathcal{L}$ on this page are the discretization operator $D$ of Section 2 (page 2), divided by the time step. ↩ ↩²

6. Page 551, "Stability", of Lyuu's notes. The explicit method for the diffusion equation $D\,\partial^2\theta/\partial x^2 - \partial\theta/\partial t = 0$ is unstable unless $\Delta t \le (\Delta x)^2/(2D)$, and halving $\Delta x$ multiplies the computation time by 8. ↩ ↩²

7. The condition obtained by putting the coefficient of the second derivative in the log-price equation, $\sigma^2/2$, into the source's $D$. For $h = 0.02$ in the common example, the condition is $\Delta\tau \le 0.01$. The calculation with 40 steps ($\Delta\tau = 0.025$) diverged, and those with 100 steps or more did not. ↩

8. Pages 553–556, "Implicit Methods", of Lyuu's notes. A tridiagonal system of equations can be solved in $O(N)$ time and $O(N)$ memory. ↩

9. Page 554, "Implicit Methods (continued)", of the same notes. The source calls the implicit formula "unconditionally stable". ↩

10. Section 1 of Giles (page 2). The source shows unconditional stability in the $L^2$ norm for $\theta \ge 1/2$ and the error orders ($O(h^2, k^2)$ for the Crank-Nicolson method, $O(h^2, k)$ for other $\theta$). It says that $k$ can be chosen in proportion to $h$. $k$ is $\Delta\tau$ on this page. ↩

11. Page 552, "Explicit Method and Trinomial Tree", of Lyuu's notes. Under the stability condition the three coefficients lie between 0 and 1 and sum to 1, so the explicit formula is the same as a backward calculation on a trinomial tree. The sum with a discounting term is calculated from the formula on this page. ↩

12. Rutgers University lecture notes "Numerical Solution of PDEs" (page 65, "Example 1: explicit scheme for the heat equation", and page 68, "Example: implicit scheme for the heat equation"). The source's $\sigma$ is the diffusion coefficient, which corresponds to $D$ on this page, and $k$ is the time step. ↩

13. Section 1.5 (page 8) of Holst's lecture notes "Computational Finance" (UC San Diego, 2015). The source sets the call's lower end and the put's upper end to 0 and calls the limits at both ends Dirichlet boundary conditions. ↩ ↩²

14. Section 2 of Giles (page 2). The source sets the initial value of the European call at expiry as $V(S, 0) = \max(S - K, 0)$. The text sets the put's $\max(K - S, 0)$ on the same reasoning. ↩

15. A calculation for the conditions in the text ($h = 0.02$, 200 time steps, Crank-Nicolson method). With the range of $x$ at $\pm 0.5$ the price is 8.4234, and with $\pm 1$ and $\pm 1.5$ it is 8.4235 (rounded at the fifth decimal place). The range of $S$ is 36.79 to 271.83. ↩

16. The value 8.433319 from the Black-Scholes formula ($d_1 = (\ln(S/K) + (r + \sigma^2/2)T)/(\sigma\sqrt{T})$, $d_2 = d_1 - \sigma\sqrt{T}$), rounded at the fifth decimal place. ↩

17. The weights and prices were calculated from the formulas in the text. $100 \times e^{-0.00025} = 99.9750$ and $110.5171 - 99.9750 = 10.5421$. ↩

18. A calculation for the conditions in the text. The price is the value at the grid point $x = 0$ ($S = 100$). The difference is the price rounded at the fifth decimal place minus the Black-Scholes price 8.4333. The explicit price with 40 steps was $-4.86 \times 10^{22}$. ↩

19. A calculation for the conditions in the text. The difference is the price rounded at the fifth decimal place minus 8.4333. The ratios of the absolute differences, from the previous row, are 4.2, 4.0, and 4.1. ↩

20. Page 559, "Numerically Solving the Black-Scholes PDE", of Lyuu's notes. The source gives the American put with the value at expiry $P(S, T) = \max(X - S, 0)$ and, before expiry, the condition $P(S, t) = \max(P(S, t), X - S)$. $X$ is the strike price $K$ on this page, and the source has no procedure for the replacement on a grid. ↩

21. Section 2.4, "American Options as a Free Boundary Problem" (page 17), and the opening (page 9) of Holst's lecture notes. American options become a free boundary problem. The section on how to handle it is unfinished. ↩

22. A calculation for the conditions in the text (put, Crank-Nicolson method). For the American put, after each time step the value is replaced by $\max(\text{value at that point}, K - S)$, with the lower end set to $K - S$ and the upper end to 0. The binomial tree is a CRR tree (2,000 steps) with a price of 7.5127, and the Black-Scholes put is 7.4383. ↩ ↩² ↩³

23. Section 4 of Giles (page 4). A nonlinear penalty term approximates American options. It can be ignored where the option is not exercised and makes the solution almost equal to the exercise value elsewhere. The solution method is Newton iteration. ↩

24. Pages 562–564 ("Monte Carlo Simulation", "The Big Idea") of Lyuu's notes. The notes explain the Monte Carlo method as estimating an expected value by a sample average. ↩

25. Section 4 of Giles (page 4). In higher dimensions solving the system of equations directly costs much more, and approximate factorization by ADI and iterative solvers such as BiCGStab with ILU preconditioning are used. The source states that the cross derivative terms need the Craig-Sneyd treatment to keep second-order accuracy. ↩ ↩²

26. Pages 566–567 ("Accuracy and Number of Replications") of Lyuu's notes. The statistical error of a sample average decreases at the rate $1/\sqrt{N}$ regardless of the dimension $n$, while ordinary numerical integration suffers the curse of dimensionality. The notes say that the Monte Carlo method is more efficient than other procedures for securities that depend on several assets. ↩

27. Page 570, "Pricing American Options", of Lyuu's notes. The standard Monte Carlo method is not suited to American options with early exercise because the exercise point is hard to determine from a single path. A modification by Longstaff and Schwartz (2001) can value them with small bias. ↩

28. Section 2 of Giles (pages 2–3). For a European call, delta and gamma computed by central differences deviate more than the price, and the largest gamma error appears at the strike ($S = 1$ in the source). ↩ ↩² ↩³

29. Pages 571–572 ("Delta and Common Random Numbers") of Lyuu's notes. Computing delta from two separate simulations has a large variance and is not recommended, and using the same random numbers is better. ↩

30. Section 1, Introduction (page 4 of the PDF), of Paweł Przybyłowicz, "Monte Carlo methods and stochastic simulations: from integration to approximation of SDEs" (arXiv:2208.05531). Deterministic algorithms are affected by the curse of dimensionality. ↩

31. Section 1.3 of Holst's lecture notes (page 5). The source lists, among the assumptions of the Black-Scholes-Merton market model, that the interest rate and the volatility are constant on $[0, T]$, calling this "an idealization, but required by the model". ↩

32. Section 3 of Giles (pages 3–4). The source describes Rannacher's method (replacing the first time step with two half steps of the implicit method) and the method shown by Giles and Carter of replacing the first two Crank-Nicolson steps with four half steps. The latter gives second-order accuracy in the maximum norm for the price, delta, and gamma of puts, calls, and digitals. ↩