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) 2 σ ( x, t ) 2 ∂ 2 f ∂f ∂t ( x, t ) + µ ( x, t ) ∂f ∂x ( x, t ) + 1 ∂x 2 ( 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 1 In other words, the risk-neutral dynamics of the state variable is assumed to be dx t = µ ( x t , t ) dt + σ ( x t , t ) dz t , where z = ( z t ) 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 other 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 of 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 x min ≡ x 0 , x 1 , x 2 , . . . , x J − 1 , x J ≡ x max , where x j +1 − x j = ∆ x for all j , i.e. x j = x 0 + j ∆ x . In particular, ∆ x = ( x max − x min ) /J . It is intuitively clear that a good approximation to the unknown solution requires that the probability that x takes on values greater than x max or smaller than x min has to be negligible. In some cases the choice of one or both of the boundaries x min , x max 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 of 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 of x min and x max 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 { x 0 , x 1 , . . . , x J } × { 0 , ∆ t, . . . , N ∆ t } . The value of the function f in the lattice node ( j, n ) corresponding to the x -value x j and the t -value n ∆ t is denoted by f j,n . Similarly, µ j,n denotes µ ( x j , n ∆ t ), σ 2 j,n denotes σ ( x j , n ∆ t ) 2 , and r jn represents r ( x j , 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 ( x j , n ∆ t ) are x f j,n = f j +1 ,n − f j,n D + , (3) ∆ x x f j,n = f j,n − f j − 1 ,n D − , (4) ∆ x 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 D x given by D x f j,n = f j +1 ,n − f j − 1 ,n , (5) 2∆ x 2
should often give a more precise approximation of ∂f ∂x and hence we will use that. The second-order derivative ∂ 2 f ∂x 2 is approximated by the difference operator D 2 x given by x f j,n = f j +1 ,n − 2 f j,n + f j − 1 ,n D 2 . (6) (∆ x ) 2 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 ∂f First, let us try a backward-looking difference approximation of the derivative ∂t , i.e. we replace ∂f ∂t ( x j , n ∆ t ) in (1) by t f j,n = f j,n − f j,n − 1 D − . (7) ∆ t 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 f j,n − f j,n − 1 + µ j,n D x f j,n + 1 2 σ 2 j,n D 2 x f j,n − r j,n f j,n = 0 , (8) ∆ t i.e. f j,n − f j,n − 1 f j +1 ,n − f j − 1 ,n + 1 f j +1 ,n − 2 f j,n + f j − 1 ,n 2 σ 2 + µ j,n − r j,n f j,n = 0 , j,n (∆ x ) 2 ∆ t 2∆ x which can be rewritten as 2 f j,n − 1 = α j,n f j − 1 ,n + β j,n f j,n + γ j,n f j +1 ,n , (9) where � � σ 2 α j,n = 1 (∆ x ) 2 − µ j,n j,n 2∆ t , ∆ x � � σ 2 j,n β j,n = 1 − ∆ t r j,n + , (∆ x ) 2 � � σ 2 γ j,n = 1 (∆ x ) 2 + µ j,n j,n 2∆ t . ∆ x We can now compute an approximation of f ( x, 0) by the following simple backward iterative procedure. First, put f j,N = F ( x j ) 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 f j,n − 1 is given explicitly in terms of the already computed values f j − 1 ,n , f j,n , and f j +1 ,n by (9). Hence, this method for numerical solution of the PDE is called the explicit finite difference method . 3 2 If 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 f i,j by f i +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. 3 The method is also known as the Euler method. 3
Recommend
More recommend
