Introduction to the Numerical Solution of Partial Differential - - PDF document

Introduction to the Numerical Solution of Partial Differential Equations in Finance

Claus Munk∗ October 29, 2007

This note gives a short presentation of some numerical techniques for the solution of the type of second-order partial differential equation that appears in many financial asset pricing models. For more information about the techniques discussed here and some alternatives, the interested reader is referred to Ames (1977), Smith (1985), Johnson (1990), Wilmott, Dewynne, and Howison (1993), Thomas (1995), Wilmott (1998), and Tavella and Randall (2000).

1 The problem

Suppose we want to find a function f(x, t), f : S × [0, T] → R, which solves the second-order partial differential equation (PDE) ∂f ∂t (x, t) + µ(x, t)∂f ∂x(x, t) + 1 2σ(x, t)2 ∂2f ∂x2 (x, t) − r(x, t)f(x, t) = 0, (x, t) ∈ S × [0, T), (1) with some terminal condition f(x, T) = F(x), x ∈ S, (2) where F : S → R is a known function. Here S ⊆ R. This is a well-known problem in finance, see e.g. specific problems Black and Scholes (1973), Vasicek (1977), and Cox, Ingersoll, and Ross (1981, 1985), and a more general introduction to the topic in Hull (2006) and Munk (2007). In these financial models f is the unknown price of an asset that depends on an underlying state variable x, which takes values in S with typically S = R or S = R+ ≡ [0, ∞). The function r is the short-term continuously compounded risk-free interest rate, which is also assumed to depend at most on time t and the state variable x. The function F is the payoff function of the asset so that F(x) is the state-dependent payoff of the asset at the maturity date T. Finally, µ is the drift rate of the state variable under the so-called risk-neutral probability measure, while σ denotes the absolute volatility of x.1

∗Department of Business and Economics, University of Southern Denmark, DK-5230 Odense. E-mail: cmu@sam.sdu.dk.

Internet homepage: http://www.sam.sdu.dk/ ∼ cmu

1In other words, the risk-neutral dynamics of the state variable is assumed to be dxt = µ(xt, t) dt + σ(xt, t) dzt, where

z = (zt) is a standard Brownian motion under the risk-neutral probability measure.

1

In some models (i.e. for some assumptions about the functions µ, σ, and r), the relevant PDE can be solved in closed-form for some specific assets (i.e. for some specific payoff functions F). In many

• ther cases, it is not possible to solve the PDE in closed form. Below, we discuss various techniques for

numerically computing approximations to the solution of such a PDE.

2 Discretization of the problem

Our approximative solution techniques are based on a conversion of the problem (1)–(2) to a sequence

• f difference equations that can be solved iteratively starting with the known values at the maturity date

given by (2). For that purpose assume that the x variable can only take on the values xmin ≡ x0, x1, x2, . . . , xJ−1, xJ ≡ xmax, where xj+1 − xj = ∆x for all j, i.e. xj = x0 + j∆x. In particular, ∆x = (xmax − xmin)/J. It is intuitively clear that a good approximation to the unknown solution requires that the probability that x takes on values greater than xmax or smaller than xmin has to be negligible. In some cases the choice of one or both of the boundaries xmin, xmax is natural, but in other problems a subjective choice has to be made. In many situations you care primarily about the value of the unknown function for one specific values

• f the state variable, e.g. f(x∗, 0) where x∗ is the current value of the state variable. In that case the

boundaries should be imposed sufficiently far from x∗ so that the approximate solution f(x∗, 0) is fairly insensitive to small changes in those boundaries. In general, some experimentation with different values

• f xmin and xmax is often useful. Furthermore, assume that the time variable can only take the values

0, ∆t, . . . , N∆t = T (N + 1 time points), i.e. ∆t = T/N. Hence, the state space S × [0, T] is approximated by the lattice {x0, x1, . . . , xJ} × {0, ∆t, . . . , N∆t}. The value of the function f in the lattice node (j, n) corresponding to the x-value xj and the t-value n∆t is denoted by fj,n. Similarly, µj,n denotes µ(xj, n∆t), σ2

j,n denotes σ(xj, n∆t)2, and rjn represents

r(xj, n∆t). The basic idea of the approximative procedure is to replace the partial derivatives in the PDE by

• differences. First consider the partial derivative ∂f

∂x . Two obvious candidates for an approximation of ∂f ∂x (xj, n∆t) are

D+

x fj,n = fj+1,n − fj,n

∆x , (3) D−

x fj,n = fj,n − fj−1,n

∆x , (4) where D+

x is called the forward-looking difference operator and D− x is called the backward-looking differ-

ence operator. A simple graph indicates that the central difference operator Dx given by Dxfj,n = fj+1,n − fj−1,n 2∆x , (5) 2

should often give a more precise approximation of ∂f

∂x and hence we will use that. The second-order

derivative ∂2f

∂x2 is approximated by the difference operator D2 x given by

D2

xfj,n = fj+1,n − 2fj,n + fj−1,n

(∆x)2 . (6) Finally, we have to approximate the derivative ∂f

∂t appearing in the PDE (1). Here, the two obvious

choices are a backward-looking difference and a forward-looking difference. In the following two sections we will try both.

3 The explicit finite difference approach

First, let us try a backward-looking difference approximation of the derivative

∂f ∂t , i.e. we replace ∂f ∂t (xj, n∆t) in (1) by

D−

t fj,n = fj,n − fj,n−1

∆t . (7) Substituting the approximations (5), (6), and (7) into the PDE (1) corresponding to the node (j, n) for 0 < j < J, 0 ≤ n < N, we get fj,n − fj,n−1 ∆t + µj,nDxfj,n + 1 2σ2

j,nD2 xfj,n − rj,nfj,n = 0,

(8) i.e. fj,n − fj,n−1 ∆t + µj,n fj+1,n − fj−1,n 2∆x + 1 2σ2

j,n

fj+1,n − 2fj,n + fj−1,n (∆x)2 − rj,nfj,n = 0, which can be rewritten as2 fj,n−1 = αj,nfj−1,n + βj,nfj,n + γj,nfj+1,n, (9) where αj,n = 1 2∆t

• σ2

j,n

(∆x)2 − µj,n ∆x

• ,

βj,n = 1 − ∆t

• rj,n +

σ2

j,n

(∆x)2

• ,

γj,n = 1 2∆t

• σ2

j,n

(∆x)2 + µj,n ∆x

• .

We can now compute an approximation of f(x, 0) by the following simple backward iterative procedure. First, put fj,N = F(xj) for all j in accordance with the terminal condition (2). Then use (9) in order to successively go backwards time step by time step until time 0 is reached. Note that the “new” value fj,n−1 is given explicitly in terms of the already computed values fj−1,n, fj,n, and fj+1,n by (9). Hence, this method for numerical solution of the PDE is called the explicit finite difference method.3

2If you compare this with Equation (17.34) in Hull (2006), you will notice a small deviation. In the last equation on

• p. 423 he should either add one to the time index on the right-hand side (replace his fi,j by fi+1,j) or subtract one from

all the time indices on the left-hand side. His method is therefore not really what other people refer to as the explicit finite difference method, although implication of this deviation for the approximate solution is probably quite small.

3The method is also known as the Euler method.

3

Equation (9) is no good for j = 0 or j = J since the right-hand side will then involve function values

• utside the lattice. In Section 6 we discuss how to determine fJ,n and f0,n. If we only care about a single

value f(x∗, 0) = f(x0 + j∗∆x, 0), e.g. corresponding to the current value x∗ of the state variable, it is sufficient to compute fj,n in the triangle defined by the nodes (j∗, 0), (j∗ + N, N), and (j∗ − N, N). In that case the explicit finite difference method is equivalent to the use of a trinomial tree, cf. Hull and White (1990) or Hull (2006, Ch. 28). It seems obvious to expect and require that a finer lattice (more nodes) should produce a better approximation of the unknown function f(x, t). In particular, we would like the approximate solution to converge to the exact solution when ∆x and ∆t approach zero. However, it can be shown that in

• rder to ensure the convergence of the explicit finite difference method, ∆t has to be very small relative

to ∆x. Loosely speaking, you have to take a lot of time steps to end up with a good approximation to the unknown solution. (If you have implemented the method with too large time steps, it will often be obvious since the results will be clearly unrealistic.) Although every time step only involves simple computations, cf. (9), the total procedure can thus be quite time-consuming. This conflicts with the desire to get the approximate solution fast. Note that the requirement that ∆t is small relative to ∆x is generally not satisfied by the trinomial procedure mentioned above. However, Hull and White (1993) show how to adjust the trinomial procedure to ensure convergence.

4 The implicit finite difference approach

Now let us try to approximate the partial derivative with respect to time by a forward-looking differ- ence, i.e. we replace ∂f

∂t (xj, n∆t) in (1) by

D+

t fj,n = fj,n+1 − fj,n

∆t . (10) Substituting (5), (6), and (10) into the PDE (1) we get the following equation corresponding to the node (j, n) for 0 < j < J, 0 ≤ n < N: fj,n+1 − fj,n ∆t + µj,nDxfj,n + 1 2σ2

j,nD2 xfj,n − rj,nfj,n = 0,

(11) i.e. fj,n+1 − fj,n ∆t + µj,n fj+1,n − fj−1,n 2∆x + 1 2σ2

j,n

fj+1,n − 2fj,n + fj−1,n (∆x)2 − rj,nfj,n = 0, which can be rewritten as aj,nfj−1,n + bj,nfj,n + cj,nfj+1,n = fj,n+1, (12) 4

where aj,n = −1 2∆t

• σ2

j,n

(∆x)2 − µj,n ∆x

• ,

bj,n = 1 + ∆t

• rj,n +

σ2

j,n

(∆x)2

• ,

cj,n = −1 2∆t

• σ2

j,n

(∆x)2 + µj,n ∆x

• .

Again we want to apply a backward iterative procedure beginning with the known function values at the terminal date, i.e. fj,N = F(j∆x). Suppose we know fj,n+1 for all j and thus want to compute fj,n for all j. In contrast to the explicit approach we cannot directly compute fj,n from values already known. But since (12) has to hold for all j = 1, . . . , J − 1, we have a system of linked equations involving the unknown function values fj,n. More precisely, we have J−1 equations in the J+1 unknowns f0,n, . . . , fJ,n. Therefore we need to add two more equations (linearly independent of the other equations) or to fix the values of two of the unknowns. In particular, if we add two equations of the form aJ,nfJ−1,n + bJ,nfJ,n = dJ,n+1 (13) and b0,nf0,n + c0,nf1,n = d0,n+1, (14) the full system of equations will have a particularly nice structure, which simplifies the solution. In order to see this, write the equations (9), (13), and (14) in matrix form as                              b0,n c0,n . . . a1,n b1,n c1,n . . . a2,n b2,n c2,n . . . . . . ... ... ... . . . . . . ... ... ... . . . . . . aJ−1,n bJ−1,n cJ−1,n . . . aJ,n bJ,n                                                           f0,n f1,n f2,n . . . . . . fJ−1,n fJ,n                              =                              d0,n+1 d1,n+1 d2,n+1 . . . . . . dJ−1,n+1 dJ,n+1                              . (15) Here dj,n+1 = fj,n+1 for j = 1, . . . , J − 1. Given fj,n+1 for all j, we can compute fj,n for all j by solving the equation system (15). The matrix is tridiagonal, which simplifies the solution of the system, as will be discussed in Section 7. In Section 6 we discuss how to come up with equations of the form (13) and (14) at the boundaries of the lattice. The method described above is called the implicit finite difference method.4 Since it involves

4The method is also known as the backward Euler method.

5

solving a system of equations in every time step, it is more complicated to implement than the explicit finite difference method. On the other hand, it can be shown that it is not necessary to use very small time steps relative to ∆x to ensure convergence of the approximate solution generated by the implicit method to the true solution. Hence, the implicit method is considered superior to the explicit method.

5 The Crank-Nicolson approach

Equation (8) is the key equation in the explicit finite difference method. If we replace n by n + 1, it reads fj,n+1 − fj,n ∆t + µj,n+1Dxfj,n+1 + 1 2σ2

j,n+1D2 xfj,n+1 − rj,n+1fj,n+1 = 0.

(16) In the implicit finite difference method the key equation is given by (11), i.e. fj,n+1 − fj,n ∆t + µj,nDxfj,n + 1 2σ2

j,nD2 xfj,n − rj,nfj,n = 0.

(17) The so-called Crank-Nicolson method is based on taking “the average” of the two equations (16) and (17), which gives fj,n+1 − fj,n ∆t + 1 2 {µj,n+1Dxfj,n+1 + µj,nDxfj,n} + 1 2 1 2σ2

j,n+1D2 xfj,n+1 + 1

2σ2

j,nD2 xfj,n

• − 1

2 {rj,n+1fj,n+1 + rj,nfj,n} = 0. (18) After substitution of the difference operators and some manipulations, we arrive at the equation Aj,nfj−1,n + Bj,nfj,n + Cj,nfj+1,n = −Aj,n+1fj−1,n+1 + B∗

j,n+1fj,n+1 − Cj,n+1fj+1,n+1,

(19) where Aj,n = 1 4∆t

• σ2

j,n

(∆x)2 − µj,n ∆x

• ,

Bj,n = −1 − 1 2∆t

• σ2

j,n

(∆x)2 + rj,n

• ,

Cj,n = 1 4∆t

• σ2

j,n

(∆x)2 + µj,n ∆x

• ,

B∗

j,n = −1 + 1

2∆t

• σ2

j,n

(∆x)2 + rj,n

• .

Again we can perform successive backwards iterations, where the right-hand side of (19) is known and the left-hand side involve the unknowns fj−1,n, fj,n, and fj+1,n. Again we need to add two equations of the form AJ,nfJ−1,n + BJ,nfJ,n = dJ,n+1 and B0,nf0,n + C0,nf1,n = d0,n+1 6

in order to “complete the system.” The resulting system of equations has a tridiagonal structure as for the implicit method. Just as for the implicit method, convergence of the Crank-Nicolson method does not require par- ticularly small time steps. The Crank-Nicolson method will typically converge faster than the implicit method and will thus typically give a more precise approximation of the unknown solution for the same lattice.

6 How to handle the boundaries

In some cases it seems reasonable to fix the value of the asset at the boundaries of the lattice, i.e. in the nodes (J, n) and (0, n). As an example, consider the pricing of a European put option on a stock in the Black-Scholes-Merton model where the state variable x is the price of the underlying stock, and the interest rate r is assumed constant. Since the stock price is lognormally distributed in that model, the obvious choice is to set the lower bound at xmin = 0 and to let f0,n = Ke−r[T −n∆t], where K is the exercise price of the put option. If the stock price should reach zero, it will stay zero, so that the payoff

• f the put option will be K. Conversely, we can set fJ,n = 0 since the put option will be worthless for

very high values of the stock. In many situations, however, it seems difficult to fix values of the unknown function f(x, t) on the boundaries, but it may be reasonable to impose conditions on the derivatives of the function at the

• boundaries. For some assets it may be a reasonable approximation to impose the condition that the

second-order derivative ∂2f

∂x2 is zero at the upper bound. If, furthermore, the first-order derivative ∂f ∂x at

the upper bound is approximated using a one-sided, backward-looking difference as in (4), the implicit finite difference method results in fJ,n+1 − fJ,n ∆t + µJ,n fJ,n − fJ−1,n ∆x − rJ,nfJ,n = 0, which can be rewritten as ∆t ∆xµJ,nfJ−1,n +

• 1 + rJ,n − ∆t

∆xµJ,n

• fJ,n = fJ,n+1.

(20) Note that this equation is of the form (13).

7 How to solve a tridiagonal system of equations

With the implicit finite difference method and the Crank-Nicolson method a number of tridiagonal systems of equations like (15) have to be solved. While a mathematician would attack the problem by computing the inverse of the matrix, this is not the most efficient way to solve the system of equations

• n the computer. A more efficient solution can be implemented as follows. First, perform a Gauss

7

elimination:5 for j = 0, 1, . . . , J − 1: aj+1,n := aj+1,n/bj,n, bj+1,n := bj+1,n − aj+1,ncj,n. Second, perform a forward substitution: for j = 0, 1, . . . , J − 1: dj+1,n+1 := dj+1,n+1 − aj+1,ndj,n+1. Third, perform a backward substitution: fJ,n := dJ,n+1/bJ,n, for j = J − 1, J − 2, . . . , 1, 0: fj,n := [dj,n+1 − cj,nfj+1,n] /bj,n. Note that if the matrix is the same at all time steps, i.e. the coefficients aj,n, bj,n, and cj,n are independent

• f n, it is only necessary to do the Gauss elimination once. Also note when several systems of equations

involving the same matrix—but different right-hand sides—a forward and a backward substitution have to be done for each right-hand side, while the Gauss elimination is the same in all cases. This is relevant if several assets are to be priced in the same model.

8 American-style derivatives

For a derivative asset of the American type the pricing function f(x, t) has to satisfy the PDE (1) with terminal condition (2) and an extra boundary condition of the form f(˜ x(t), t) = F(˜ x(t), t), where F(x, t) gives the payoff if the asset is exercised at time t if a situation where the state variable has the value x. The extra condition is to be understood as follows. If the value of the state variable reaches a certain critical level ˜ x(t), it is optimal to exercise the derivative asset pre-maturely and cash in the payoff F(˜ x(t), t). However, the critical values ˜ x(t) are not known in advance, but has to be computed as a part of the solution. This can be done in the backward iterative procedure by checking whether early exercise is optimal—just as in a binomial tree for option pricing. Suppose we know the value of the American-style derivative at time (n + 1)∆t, i.e. we know fj,n+1 for all j. First, compute a value fj,n for all j as explained in the previous sections. Then replace that value with max(fj,n, F(xj, n∆t)). Next, we can go backwards another time step.

5The symbol := should be interpreted as the assignment operator as it would be understood by a computer. As an

example, the declaration aj+1,n := aj+1,n/bj,n means that the variable aj+1,n is assigned a new value equal to its current value divided by the value of the variable bj,n. Hence, this new value of aj+1,n should be used in all later calculations until a new value is assigned to the variable.

8

The procedure just outlined is first to solve the system of equations and then adjust the solution values fj,n for all j. For the implicit method and the Crank-Nicolson method the early exercise feature can be included in a more efficient way by checking for early exercise during the solution of the equation

• system. This can be obtained by adjusting the backward substitution slightly:

fJ,n := max (dJ,n+1/bJ,n, F(xJ, n∆t)) , for j = J − 1, J − 2, . . . , 1, 0: fj,n := max ([dj,n+1 − cj,nfj+1,n] /bj,n, F(xj, n∆t)) . Now, when the value of fj,n is computed, the term fj+1,n on the right-hand side has already been modified to take into account the possibility of early exercise in the node (j + 1, n). The computed value of fj,n is thus based on a better approximation of fj+1,n.

9 Assets with intermediate payments

It can be shown that the price f(x, t) of an asset giving a continuous dividend yield q(x, t) has to satisfy the PDE ∂f ∂t (x, t) + µ(x, t)∂f ∂x(x, t) + 1 2σ(x, t)2 ∂2f ∂x2 (x, t) − (r(x, t) − q(x, t))f(x, t) = 0, (x, t) ∈ S × [0, T) with appropriate boundary conditions. The extra term will not cause any problems in the finite difference methods discussed above. Next consider an asset giving discrete payments at predetermined points in time, say F(x(τi), τi) at time τi, i = 1, 2, . . . , n, where 0 < τ1 < · · · < τn = T. For example, this is relevant for a coupon

• bond. The price f(x, t) of such an asset has to satisfy the PDE in between payment dates, i.e. in all

intervals (τi, τi+1). At a given payment date the price of the asset will drop by an amount exactly equal to the size of the payment. We can price such an asset using one of the backward iterative procedures described in earlier sections with the following small adjustment for each intermediate payment F(x, τi). Suppose that τi ∈ [n∆t, (n+1)∆t) and that the values fj,n+1 for all j have already been computed. Then the values at time n∆t are determined by first computing fj,n as described earlier and then replacing fj,n by fj,n + F(xj, τi)e−rj,n[τi−n∆t].

References

Ames, W. F. (1977). Numerical Methods for Partial Differential Equations (Second ed.). Academic Press. Black, F. and M. Scholes (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy 81(3), 637–654. Cox, J. C., J. E. Ingersoll, Jr., and S. A. Ross (1981). A Re-examination of Traditional Hypotheses about the Term Structure of Interest Rates. The Journal of Finance 36(4), 769–799. 9

Cox, J. C., J. E. Ingersoll, Jr., and S. A. Ross (1985). A Theory of the Term Structure of Interest

• Rates. Econometrica 53(2), 385–407.

Hull, J. and A. White (1990). Valuing Derivative Securities Using the Explicit Finite Difference Method. Journal of Financial and Quantitative Analysis 25(1), 87–100. Hull, J. and A. White (1993). One-Factor Interest-Rate Models and the Valuation of Interest-Rate Derivative Securities. Journal of Financial and Quantitative Analysis 28(2), 235–254. Hull, J. C. (2006). Options, Futures, and Other Derivatives (6th ed.). Prentice-Hall, Inc. Johnson, C. (1990). Numerical Solution of Partial Differential Equations by the Finite Element Method. Cambridge University Press. Longstaff, F. A. and E. S. Schwartz (1992). Interest Rate Volatility and the Term Structure: A Two- Factor General Equilibrium Model. The Journal of Finance 47(4), 1259–1282. Munk, C. (2007). Introduction to Fixed-Income Analysis. Lecture notes, University of Southern Den- mark. Smith, G. D. (1985). Numerical Solution of Partial Differential Equations (Third ed.). Oxford Applied Mathematics and Computing Science Series. Oxford University Press. Tavella, D. and C. Randall (2000). Pricing Financial Instruments: The Finite Difference Method. Wiley. Thomas, J. W. (1995). Numerical Partial Differential Equations, Volume 22 of Texts in Applied Math-

• ematics. Springer-Verlag.

Vasicek, O. (1977). An Equilibrium Characterization of the Term Structure. Journal of Financial Economics 5, 177–188. Wilmott, P. (1998). Derivatives: The Theory and Practice of Financial Engineering. Oxford University Press. Wilmott, P., J. Dewynne, and S. Howison (1993). Option Pricing: Mathematical Models and Compu-

• tation. Oxford Financial Press.

10