- Email: [email protected]
Implicit–explicit Runge–Kutta methods for financial derivatives pricing models
European Journal of Operational Research 171 (2006) 991–1004 www.elsevier.com/locate/ejor
Implicit–explicit Runge–Kutta methods for financial derivatives pricing models Javier de Frutos
*
Departamento de Matema´tica Aplicada, Universidad de Valladolid, C/Prado de la Magdalena S/N, 47005 Valladolid, Spain Available online 10 March 2005
Abstract Implicit–explicit Runge–Kutta methods are investigated for application to financial derivatives pricing models in the partial differential equations approach. The methods are showed to be an alternative to other existing procedures for the numerical valuation of American type contracts. We follow the method of lines in order to have a numerical method that can be used with a variety of state variable discretizations including finite elements, finite differences and finite volume methods. Some numerical experiments are presented. 2005 Elsevier B.V. All rights reserved. Keywords: Finance; American options; Implicit–explicit Runge–Kutta methods; Finite element method
1. Introduction A wide variety of derivative securities traded in exchange markets such as call or put options on dividend paying stocks, foreign currency options, callable bonds, among others, are American type contracts. An American contract gives the holder the right to exercise and exit the contract at any time before the expiration date specified in the contract. The problem of valuating such contracts is considerably more difficult than the corresponding European contract where the holder can only exercise at expiration. In fact there is no closed form analytical formulas for valuing American contracts even in the simpler cases. Numerical methods are then necessary for pricing such contracts. Derivative securities can usually be defined as optimal stopping problems [18] of the form V ðx; tÞ ¼ sup Ex ½gðs; X s Þ;
x 2 X;
s2Tðt;T Þ
*
Tel.: +34 983 42 31 81; fax: +34 983 42 30 13. E-mail address: [email protected]
0377-2217/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2005.01.013
ð1Þ
992
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
associated with a system of stochastic differential equations dX t ¼ bðX t ; tÞ dt þ rðX t ; tÞ dW t ;
ð2Þ
and a reward function g, where X is a possibly unbounded subset of Rn and Tðt; T Þ is the set of stopping times. Here, Ex denotes the mathematical expectation with respect the probability law of the process Xt, starting at X0 = x. In many cases the reward function is chosen to be the discounted value of a given payoff w. Z s gðs; X s Þ ¼ exp rðs; X s Þ ds wðs; X s Þ: t
The value function is representable deterministically as the solution of a system of time dependent variational inequalities [2,9,12] in terms of the characteristic differential operator associated with (2). In this formulation the problem is similar to a parabolic, possibly degenerate, obstacle problem in which the state variables play the role of spatial degrees of freedom. Most practical numerical methods for the solution of the equivalent linear complementarity problem use well known finite differences approaches, [10,23,24]. In [26,8] a penalty method coupled with a finite volume discretization of spatial like variables and finite differences in time are proved to be a useful and very efficient alternative for this kind of problems. The main idea is to change the set of variational inequalities into a set of nonlinear partial differential equations through a penalization of the restriction [2,17,20]. As it is claimed in [8], the main advantage of the penalty method is that it is easily generalizable to multidimensional problems and a variety of contracts. Furthermore, it can be used coupled with any type of spatial (state) discretization, finite differences, finite volume or finite element methods. In this paper we are concerned with implicit–explicit Runge–Kutta methods for the numerical integration of the system of ordinary differential equations which arises after spatial discretization of the penalty approximation to the linear complementarity problem. That is, we use the method of lines to decouple the spatial like and temporal discretization allowing the use of efficient stiffly accurate time integrators others than simple finite differences. We remark that, in the so-called continuation region, one has to numerically solve a system of convection–diffusion–reaction partial differential equations. This problem qualifies as stiff, see [11]. So, stiffly accurate time integrators are necessary to have efficient numerical methods. Note that efficiency is a crucial issue in order to design numerical methods of practical use in general, possibly multidimensional, financial problems. Surprisingly, it seems that there has not been any attempt to use efficient, stiffly accurate time integrators other than the first-order backward Euler method in the financial literature. The only exception is [19], where the authors propose to use a numerical differentiation formula with quasiconstant step size and variable order. Implicit–explicit Runge–Kutta methods are a particular instance of additive Runge–Kutta methods [7]. It has been proved in [5], see also [1,14], that being easily implementable, implicit–explicit Runge–Kutta methods are competitive with standard stiffly accurate time integrators for convection–diffusion–reaction equations. To our knowledge this is the first time that they have been successfully applied to the numerical solution of time dependent variational inequalities. Recently, implicit–explicit finite difference methods has been successfully used for European and barrier options in jump diffusion and exponential Le´vy models in [6]. In the following, for clarity in the exposition, we restrict ourselves to a simple one dimensional model problem, the American put option, although the numerical methods is intended to be general enough to be applied to more complex financial problems. In Section 2 we state the model problem and its penalty approximation. Implicit–explicit Runge–Kutta time integrators coupled with a finite element discretization are presented in Section 3. The next section is devoted to show the qualities of the method proposed by means of some numerical experiments. We finish the paper with some conclusions and final remarks.
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
993
2. The model problem Consider an asset with price S which is supposed to follow a log-normal evolution process dS ¼ lS ds þ rS dW ;
ð3Þ
where l is the expected rate of return, r is the volatility, dW is the increment of a Wiener process and s denotes time. We are concerned with the problem of determining the fair value V(S, s) of an American put option in which the holder can exercise at any time, receiving a payoff specified in the contract as w(S, s) = max(E S, 0), where E is the strike. Under the usual assumptions on the market, the problem can be formally reduced to finding V(S, s) such that V t LBS V P 0;
0 6 S 6 1; 0 < t 6 T ;
V ðS; tÞ P wðS; tÞ;
0 6 S 6 1; 0 < t 6 T ;
ðV t LBS V ÞðV wÞ ¼ 0; V ðS; 0Þ ¼ wðSÞ;
0 6 S 6 1; 0 < t 6 T ;
0 6 S 6 1:
ð4Þ ð5Þ ð6Þ ð7Þ
Here, LBS denotes the well-known Black–Scholes operator, 1 LBS V ¼ r2 S 2 V SS þ rSV S rV : 2
ð8Þ
r is the risk free interest rate and t = T s is the time to expiration of the contract. The problem has a unique solution which is smooth in the continuation region and is continuously differentiable globally [12,2, Chapter 3]. Note that no boundary conditions are imposed as corresponds to a singular problem. From a practical point of view the solution of (4)–(7) can be approximated by means of a sequence of problems posed over bounded increasing domains. More precisely, let us consider a sequence of domains ½S km ; S kM ð0; 1Þ with S km ! 0 and S kM ! 1. Let Vk be the solution of the localized problem S km 6 S 6 S kM ; 0 < t 6 T ;
V kt LBS V k P 0;
S km 6 S 6 S kM ; 0 < t 6 T ;
V k ðS; tÞ P wðS; tÞ;
ðV kt LBS V k ÞðV k wÞ ¼ 0; V k ðS km ; tÞ ¼ E S km ; V k ðS kM ; tÞ ¼ 0;
0 < t 6 T;
0 < t 6 T;
V k ðS; 0Þ ¼ wðSÞ; Let X ¼ ½S l ; S u
S km 6 S 6 S kM ; 0 < t 6 T ;
S km 6 S 6 S kM :
ðS km ; S jM Þ k
ð10Þ ð11Þ ð12Þ ð13Þ ð14Þ
a fixed approximation domain. Then from [16,12,13],
lim max kV ð; tÞ V ð; tÞkL1 ðXÞ ¼ 0:
k!1 t2½0;T
ð9Þ
ð15Þ
We note that the boundary conditions (12) and (13) are essentially arbitrary but nevertheless financially reasonable. Furthermore, in [16] it is shown that every choice of the boundary conditions in an increasing exhausting sequence of bounded domains leading to well posed systems can be used to approximate V(S, t) in a given fixed bounded approximation domain. However, the choice of boundary conditions surely affects the quality of the approximation. The rationale for the choice (12), (13) is as follows: First it is clear that the
994
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
value of the option must be zero when the asset value tends to infinity. The boundary condition in the neighborhood of S = 0 is slightly more subtle. A wrong value will generate a boundary layer near S = 0 that is, obviously, not financially acceptable. From the complementarity condition (6) one has that, at S = 0, either Vt + rV = 0 or V(0, t) = w(0, t) = E. Using (7), in the first case it is found that V(0, t) = exp(rt)E < E for all time t > 0. Then, the constraint (5) forces V(0, t) = E and in consequence (12) for S km small enough. Other boundary conditions have been used in the literature. In [21] linearity type boundary conditions of the form VSS = 0, together with some other useful alternatives, are advocated for options with piecewise linear payoff when finite differences are used as discretization tool. However, this type of boundary conditions are more difficult to implement in the context of a finite element method. 2.1. The penalty problem There are several approaches that can be used to solve the localized problem (9)–(14). Here, we shall use the penalization technique that has been proved to be an efficient way to treat some financial problems, see [8]. To this end, let us start by choosing a smooth function u with uðS km Þ ¼ E S km , uðS kM Þ ¼ 0, and otherwise arbitrary. For > 0, we seek Vk,, the unique solution of the (regularized) penalty problem k; U k; þ b ðU k; Þ ¼ f ; t LBS U
U k; ðS km ; tÞ ¼ U k; ðS kM ; tÞ ¼ 0;
S km 6 S 6 S kM ; 0 < t 6 T ;
ð17Þ
0 < t 6 T;
U k; ðS; 0Þ ¼ w ðSÞ ¼ wðSÞ uðSÞ;
ð16Þ
S km 6 S 6 S kM ;
ð18Þ
where f = ut + LBSu, w*(S) = w(S) u(S) is the modified payoff and b is the penalty term, b ðuðSÞÞ ¼
1 minðuðSÞ w ðSÞ; 0Þ:
ð19Þ
Problem (16)–(18) has a unique solution with derivatives up to second-order in Lp ðS km ; S kM Þ for t > 0, see for example, [9, Theorem 10.4.2]. Furthermore, for small enough, there exists a constant C > 0, such that [3, Theorem 4.1] max kV k ð; tÞ ðU k; ð; tÞ þ uð; tÞÞkL1 ðXÞ 6 C:
ð20Þ
t2½0;T
Besides, the American restriction (10) is also satisfied to OðÞ terms, that is, V k; ðS; tÞ ¼ U k; ðS; tÞ þ uðS; tÞ P wðSÞ C;
S 2 ½S km ; S kM
ð21Þ
for some constant C > 0 independent of .
3. The discrete problem We follow the method of lines and discretize first the asset value and then time. One of the advantages of this methodology is that decoupling the discretization problem into a spatial-like part and a temporal one allows an easier treatment of the specific characteristic of each problem. Here, we discretize the asset value by means of the linear finite element method. In order to simplify somewhat the exposition we take S km ¼ 0 in the rest of the paper. For this particular problem it can be observed that in practice all results of the previous section remain valid.
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
995
Let us consider a partition 0 ¼ S 0 < S 1 < < S J þ1 ¼ S kM of the computational domain, and let us denote h = max06i6J(Si+1 Si). We consider the space of test functions n o X h ¼ / 2 Cð½0; S kM Þj/j½S i ;S iþ1 2 P1 ; /ð0Þ ¼ /ðS kM Þ ¼ 0 ; where P1 denotes the space of polynomials of degree less than or equal to one in the asset variable S. The semidiscrete problem consists in finding a function uh : [0, T] ! Xh such that Z k Z Sk Z Sk Z Sk Z Sk M M M M r2 S M 2 uh;t / dS þ S uh;S /S dS ðr r2 Þ Suh;S / dS þ r uh / dS þ b ðuh Þ/ dS 2 0 0 0 0 0 Z S kM f / dS 8/ 2 X h : ð22Þ ¼ 0
We remark that although other variational formulations are possible (see for example [9,12,16,22]), from a computational point of view this one has the merit of being numerically efficient and directly related to the financial problem. Let /j 2 Xh, j = 1, . . . , J, be defined by /j(Sk) = djk. It is clear that Xh = span{/j, j = 1, . . . , J}. Therefore, writing uh ðS; tÞ ¼
J X
U j ðtÞUj ðSÞ;
j¼1
the variational equations (22) are equivalent to the following system of ordinary differential equations: M
dU h r2 þ KU h ðr r2 ÞCU h þ rMU h þ Mb ðU h Þ ¼ MF h ; dt 2 T
T
ð23Þ T
where U h ¼ ½U 1 ; . . . ; U J , F h ¼ ½f ðS 1 Þ; . . . ; f ðS J Þ and b ðU h Þ ¼ ½b ðU 1 Þ; . . . ; b ðU J Þ . The mass, M, stiffness, K, and convection, C, matrices are defined by Z S M J M¼ /i /j dS ; 0
K¼
Z
i;j¼1 SM
S 2 /i;S /j;S dS
J
0
C¼
Z
J
SM
; i;j¼1
S/i;S /j dS 0
: i;j¼1
For the semidiscrete approximation (22) it is possible to show that for each t > 0, the following error bound holds [4, Chapter 8]: kU k; ð; tÞ uh ð; tÞkL1 ðXÞ 6 Ch2 j log hj; and in consequence, taking (20) into account, kV k ð; tÞ ðuh ð; tÞ þ /ð; tÞÞkL1 ðXÞ 6 Cð þ h2 j log hjÞ:
ð24Þ
3.1. Implicit–explicit time discretizations The system of ordinary differential equations (23) can be solved by means of an standard ODE time integrator. The main difficulty when dealing which this kind of system is that the use of explicit time integrators
996
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
is inefficient because the system becomes stiffer as the asset value mesh is refined. As a consequence the time step must be heavily reduced in order to fulfill the severe stability restriction inherent to explicit methods. On the other hand, if a stiffly accurate, implicit time integrator is chosen, one has to solve implicit nonlinear equations with nonsmooth Jacobian. Furthermore, the convective terms introduce additional numerical difficulties. This could be one of the reasons why very few implicit codes other than mid-point rule can be found in the literature. To our knowledge perhaps the only exception is [19] where a backward differentiation type code is proposed. However, multistep based methods are reported in [25], see also [8], to perform badly with American type financial problems as shout options. Although it is not clear in the financial literature if this is also the case for simpler contracts, the use of one step time integrators present some practical advantages. Among others, there exist easily implementable methods with variable time steps and automatic error control that are efficient for convection–diffusion type equations that typically arise in financial problems. The main idea in an implicit–explicit time discretization is to use different methods for the different terms rolled into a single composite method that could treat more efficiently each one of the terms appearing in (23). To this end let us write (23) in the form M
dU h ¼ LðU h Þ þ NðU h Þ; dt
ð25Þ
where r2 KU h Mb ðU h Þ; 2 NðU h Þ ¼ ðr r2 ÞCU h rMU h þ MF h : LðU h Þ ¼
The partition (25) has been chosen such that L(Uh) collects the stiffer terms appearing in the equations (diffusive and penalty terms in our case) and N(Uh) gathers the less stiff terms (convective and source terms). The simplest implicit–explicit Runge–Kutta method is a combination of the well known explicit and implicit Euler rules. Let us suppose that at time tn < T a numerical approximation Unh U h ðtn Þ has been found. Then the approximation U hnþ1 to Uh(tn+1) = Uh(tn + Dt) is the solution of the system of algebraic equations n MU hnþ1 ¼ MU nh þ DtLðU nþ1 h Þ þ DtNðU h Þ:
ð26Þ
Note that as the most stiff terms, diffusion and penalty ones, are treated implicitly by means of a L-stable method (Backward Euler method in this simple case), only a mild stability restriction Dt 6 C independent of asset grid size has to be imposed, see [5]. This is the main advantage of implicit–explicit methods over pure explicit approaches. The price to be paid is that at each time step, a nonlinear system of equations has to be solved. However, as the convection term in (25) is treated explicitly, only a linear solver for symmetric systems is needed when (26) is solved by means of a Newton like iteration procedure. This represents usually an advantage over fully implicit methods such as multistep BDF formulas. For the numerical solution of (26), following [8], we propose to use the following fixed point Newton-like iteration: Given an initial iterate Y0, for example, Y 0 ¼ U nh , for k = 0,1, . . . , M þ DtK þ DtMBk Y kþ1 ¼ MU nh þ DtNðU nh Þ þ DtMBk w ; ð27Þ T
where w ¼ ½w ðS 1 Þ; . . . ; w ðS J Þ and Bk is the diagonal matrix defined by ( 1= if ðY ki Þj < w ðS j Þ; k Bjj ¼ 0 if ðY ki Þj P w ðS j Þ: The iteration is stopped when the difference between two consecutive iterates is small enough. We refer to Section 4 for the practical details concerning the implementation of the nonlinear iteration (27).
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
997
We note that (27) is, formally, the Newton method applied to the nonlinear system (26) in which matrix B plays the role of the Jacobian of the (nonsmooth) nonlinear penalty term. In [8] it is proved that if a lumped (diagonal) mass matrix [15] is used in (27), the iterates converge monotonically to the unique solution of (26). In our experiments we have found that with consistent mass matrix, (27) still converges in a small number of iterations. 3.2. Implicit–explicit Runge–Kutta methods The implicit–explicit Euler method is a first-order method, that is, the error decreases proportionally to time step Dt and in consequence is somewhat inefficient. However, high order one step methods sharing the same ideas that the implicit–explicit Euler methods are available in the literature. In an implicit–explicit Runge–Kutta method, the equations that describe the step tn # tn+1 when applied to (25) take the form MY 1 ¼ MU nh ; MY i ¼
MU nh
ð28Þ
þ Dtaii LðY i Þ þ Dt
i1 X
aij LðY j Þ þ
j¼2
MU hnþ1 ¼ MU nh þ Dt
sþ1 X
bi LðY i Þ þ
i¼2
i1 X
! ^ai;j NðY j Þ ;
2 6 i 6 s þ 1;
ð29Þ
j¼1 sþ1 X
! ^ bi NðY i Þ ;
ð30Þ
i¼1
b ¼ ð^aij Þ where Yi denotes the internal stages. The coefficients of the method, that is, the matrices A = (aij), A T and the vectors bT = (bi) and ^ b ¼ ð^ bi Þ are computed as to satisfy the so-called order conditions (to a certain attainable order depending on the number of stages s + 1) and some stability restrictions. Usually, the coefficients are chosen such as the following simplifying assumptions are satisfied: ci ¼
i X j¼2
aij ¼
i1 X
^ aij ;
bi ¼ ^ bi :
j¼2
We refer to [5], see also [11], for a review of the order conditions and a study of the stability properties of implicit–explicit Runge–Kutta methods. Several efficient implicit–explicit Runge–Kutta methods exist in the literature. We refer for example to [1,14], where an exhaustive comparison among different possibilities are reported. Each stage of the method is an implicit–explicit Euler like step, see Eq. (29), that is treated as has been previously described by means of an iteration similar to (27). Note that only an efficient procedure to solve symmetric positive definite linear systems is needed. Furthermore, the method is easily generalizable to other models were a special stabilization procedure of the convective terms could be needed. We do not pursue further this issue in this paper as this is not the case for the simple model of the American put option. We remark that the method (28)–(30) can be seen as pair of Runge–Kutta methods specified by the T b ^ coefficients A, bT, c and A, b , c, respectively. The first method is implicit and has to be applied to the stiff terms, in this case the diffusive and penalty terms; the second one is explicit and is applied only to the convective and source terms. Moreover, if the implicit method is L-stable [11] it is possible to construct an explicit companion method in such a way that the full implicit–explicit procedure needs only a mild stability restriction of the type Dt 6 C with C independent of the mesh size, see [5]. This is a necessary condition in order to have an efficient numerical procedures for partial differential equations. Besides, it is also a very convenient property when one is dealing with nonsmooth initial data as (18), avoiding the need of special
998
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
starting procedures required when one uses a non strongly stable procedure as for example the implicit midpoint rule (Crank–Nicolson method) [8]. Finally we remark that since the asset and time discretizations are decoupled, the same time integrator can be used in conjunction with other spatial like discretization such as finite differences, finite volume methods or even high order finite elements. 4. Numerical experiments In this section we consider the third-order method (LIRK3) with s + 1 = 4 stages constructed in [5]. The implicit part of the method is based in the well known L-stable, three stage, third-order SDIRK [11]. The coefficients of the explicit part of the method are computed in such a way that the method have quasioptimal stability properties when applied to convection–diffusion equations, see [5, Section 4], like the Black–Scholes equation. The method is implemented in variable-step mode, so that the time step is automatically chosen by the code depending on an estimation of the local errors. To this end, we have considered the second-order formula: e nþ1 ¼ MU n þ Dt b2 ðLðY 2 Þ þ NðY 2 ÞÞ þ b3 ðLðY 3 Þ þ NðY 3 ÞÞ þ b4 ðLð U e nþ1 Þ þ NðY 4 ÞÞ ; MU h h h where Y2, Y3 and Y4 are the internal stages of the third-order scheme. The second-order approximation e nþ1 can be interpreted as the solution generated with a second-order implicit–explicit Runge–Kutta U h e nþ1 defined by the coefficients a5;i ¼ ^a5;i ¼ bi , 1 6 i 6 3, a5;4 ¼ method with an additional stage Y 5 ¼ U h ^ a5;5 ¼ 0, a5;5 ¼ ^ a5;4 ¼ b4 . Each time a new approximation U nþ1 is computed, the local error is estimated by h e nþ1 k : ERR ¼ kU hnþ1 U 1 h The step is accepted if ERR 6 TOL, where TOL is the precision required by the user. If ERR > TOL the step is rejected and the computation is restarted from the previously accepted approximation using a smaller time step. In both cases the new time step size is selected according to the rule 1=3 !! TOL Dtnew ¼ min hlim; Dtold min FACMAX; FAC ; ð31Þ ERR where ‘‘hlim’’ is the maximum admissible step size (in our code hlim = T/5) and FAC and FACMAX are safety factors that limit the variation of step size in order to make the code more robust. In our code we choose FAC ¼ 0:9;
FACMAX ¼ 1:5
if the step has been accepted (EST 6 TOL), and FAC ¼ 0:7;
FACMAX ¼ 1
after a rejection (ERR > TOL). Each stage, a nonlinear problem of the form (26) has to be solved. To this end, we have used the Newton like iteration (27), stopping the iteration when max
jY kþ1 Y kj j j maxð1; jY kj jÞ
6 103 TOL;
where Y kj , j = 1, . . . , denote the components of the iterant Yk in (27).
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
999
In order to control the cost of the iterations in the solution of the nonlinear systems, we have also decided to adopt a conservative policy and limit the maximum number of iterations in each stage. If, in the computation of one stage, the nonlinear iteration has not converged after four iterations, the full time step is rejected and the computation restarted with time step halved. We have tested the method for several values of the parameters. Here, we present the results obtained for the particular values r = 0.2, r = 0.1 and E = 100. The conclusions for other values were similar to the ones reported here. In all the numerical experiments we have taken S kM ¼ 200. For this computational domain the error caused by the arbitrary boundary conditions was found to be small enough compared with the errors made in the discretization process. Fig. 1 shows the time step history for a single run using a grid with 160 nodes and setting TOL = 103 and = 104. In this experiment time to expiry was T = 1. As we can see the code selects small time steps at the beginning of the computation and smoothly increments the time step as the computation proceeds. We remark that the initial condition (18) has only piecewise smooth derivatives although is smoother for t > 0 due to the parabolic character of the Black–Scholes equation. The small steps are the response of the code to the lack of regularity at t = 0. Once the solution has been smoothed, the code increments the time step until a default in the convergence of the nonlinear iterations (27) is detected. Then, the code automatically reduces the time step, see Fig. 1 near T = 0.4, and continues the time integration increasing again the time step according to the rule (31). Note that no special procedure is needed to deal with nonsmooth initial conditions as is the case with other time integrators of common use as, for example, mid-point rule, see [8]. The time step history was found to depend only on the time integrator tolerance and is essentially independent of the asset grid size, see below. In the next test we have checked the behaviour of the method with respect to the penalty parameter. We have monitored the maximum relative error in enforcing the American constraint by the penalty method by means of the quantity AMERR ¼ max
S2½0;S kM
maxðw ðSÞ uh ðSÞ; 0Þ : maxð1; w ðSÞÞ
ð32Þ
Table 1 shows the results obtained for several values of the discretization parameters, number of nodes and tolerance in the step control procedure of the time integrator, varying the size of the penalization parameter . We observe that the size of the relative American error is proportional to , in good agreement
0.12
0.1
0.08
0.06
0.04
0.02
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 1. Time step selection history. r = 0.2, r = 0.1, K = 100, T = 1, grid with 160 nodes, TOL = 103, = 104.
1000
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
Table 1 American put option Nodes
2
10 103 104 105
TOL
American error
Ratio
Steps
· · · ·
3
10 104 104 105
www 8.00 7.62 11.07
17 17 17 27
7.57 8.70 1.04 9.64
· · · ·
103 104 104 106
www 8.70 8.36 10.79
25 25 30 32
7.56 8.63 9.42 1.25
· · · ·
103 104 105 105
www 8.76 9.16 7.54
42 46 50 52
160 160 160 160
2
10 102 102 102
7.56 9.45 1.24 1.12
102 103 104 105
320 320 320 320
103 103 104 105
102 103 104 105
640 640 640 640
104 104 104 104
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, 640 nodes, = 106.
with bound (21). Furthermore, the error is nearly independent of the discretization parameters. It is interesting to see that the number of time steps and, as a consequence, the cost of time integration is nearly independent of the penalization parameter. This is a consequence of the strong stability properties of the time integrator. The slight increment in the number of time steps observed with decreasing values of is due to the conservative policy adopted for time step control. Each time a slow convergence of the iterates in (27) is detected, the time step is halved. The number of time step increases because of this somewhat drastic reduction in the time step. In order to check the accuracy of the implicit–explicit Runge–Kutta method in this application we performed several runs for different values of the parameter TOL. The number of nodes was set to be 640 in all the reported experiments so that the numerical solutions are nearly free of spatial discretization error. The results of the experiment are shown in Table 2. The most precise result corresponds to TOL = 106 and has
Table 2 American put option TOL
Steps
2
10 103 104 105 106
V (1 0 0)
18 37 50 91 153
Difference 2.13 · 2.40 · 3.01 · 2.00 · www
3.04907 3.06794 3.07003 3.07032 3.07034
Ratio
2
www 8.88 7.97 15.05 www
10 103 104 105
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, 640 nodes, = 106.
Table 3 American put option Nodes
TOL
Steps
V (1 0 0)
Difference
40 80 160 320
102 102 103 103
17 17 30 30
3.10001 3.07736 3.07262 3.07097
2.965 7.019 2.274 5.267
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, TOL = 106, = 106.
· · · ·
102 103 103 104
Ratio www 4.2 3.1 4.3
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
1001
50
45
40
35
30
25
20
15
10
5
0 50
60
70
80
90
100
110
120
130
140
150
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-0.9
-1 50
60
70
80
90
100
110
120
130
140
150
90
100
110
120
130
140
150
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
-0.01 50
60
70
80
Fig. 2. Top: Option value (V), middle: delta (VS), bottom: gamma (VSS). Grid with 160 nodes, TOL = 103, = 104, T = 0.25.
been taken as a reference solution in order to measure the efficiency of the numerical method. For each value of the tolerance we show in Table 2 the number of steps, the computed price of the option at
1002
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
S = 100 and the difference between computed prices and the reference solution. This difference is an indicator of the error due to the temporal discretization. Last column shows the ratio of consecutive differences. We observe that the error is of the same order as the prescribed tolerance. We remark that in a step size control procedure, the objective is to have an error in the numerical solution as close as possible to the user prescribed tolerance. In this way, the user is able to optimize, for a given precision, the number of time steps and in consequence the CPU time required by the computation. The high values in the last column of Table 2 nicely show a very rapid convergence as the tolerance is reduced. Note also that each time the time integrator error is reduced by an order of magnitude, the number of steps is increased by a factor less than two. So, in spite of the some times alleged lack of regularity of the American put problem, by means of a stiffly accurate, third-order time integrator, the error can be drastically reduced with only a moderate increment in the cost of the computation. It is clear that, in a practical computation, an equilibration of the errors coming from asset and temporal discretizations is needed. We present in Table 3 the results obtained by minimizing the cost of the computation for a given error. The value of the time integrator tolerance showed in each row is the best value in the sense that a smaller value does not reduce the error while it increases the number of steps and in consequence the cost of the computation. In all runs the value of the penalization parameter was taken to be = TOL/10. As we have done before we take as an indicator of the error in the computation the difference with respect to the numerical solution computed with 640 nodes, TOL = 106 and = 107. As we can see, the ratio between two consecutive differences is very near 4 as corresponds to a second-order method. Furthermore, thanks to the high order temporal discretization, the error in the asset discretization can be reduced while maintaining a small number of time steps. To sum up, the numerical experiments show that the error in the proposed method decouples in three components, the error due to the penalty term, the temporal discretization error and the asset discretization error which can be controlled independently. The high order strongly stable method used for the time integration makes the global procedure very efficient providing rapid convergence with respect to the time integration tolerance and quadratic convergence with respect to the number of degrees of freedom used in the asset discretization. The sensitivity of the option value to variations in asset price and parameters of the model is measured by the so called Greeks [21,23]. In practice, besides the option price, accurate approximations to the Greeks are needed for edging purposes [23]. Fig. 2 shows value, delta (oV/oS) and gamma (o2V/oS2) of the option for the same values of the parameters as before and time to expiry T = 0.25. We can observe that the numerically calculated Greeks are free of spurious oscillations which are typically present in the numerical solution of American type contracts, specially in the computed gamma, due to the discontinuity of the second derivative in the optimal exercise price. This unwanted behaviour of the gamma has been reported in [8] where a second-order mid-point rule (Crank–Nicolson method) is used as time integrator. There, the authors claim that the oscillations disappear if the first-order, fully implicit, Backward Euler method is used instead. Note that the implicit–explicit Runge–Kutta method is a strongly stable method that shares the dissipative properties of implicit Euler method while it retains the efficiency derived from the third-order accuracy and variable step implementation.
5. Conclusions We have presented a third-order variable step implicit–explicit Runge–Kutta method which is showed to be an efficient alternative to other existing procedures for the numerical valuation of American put options. The method is coupled with a linear finite element discretization for the asset variable and treats the American restriction by the penalty method. The computed Greeks are free of the spurious oscillations which are
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
1003
typically present when no strongly stable time integrators, such as Crank–Nicolson method, are used, while the efficiency inherent to high order methods is retained. Using the method of lines to state the full discretization it is possible to completely decouple asset discretization from time integration. This allows the time integrator to be used with other asset discretizations such as the finite difference or finite volume method. We remark that, for the same reason, the extension to multifactor problems is also straightforward.
Acknowledgments Research supported by DGI-MCYT under grant BFM2001-2138 (cofinanced by FEDER funds) and by JCYL under project VA044/03. This research was completed while the author was visiting professor at GERAD and HEC Montre´al partially supported by grant SEEU PR2002-0223.
References [1] U.M. Ascher, S.J. Ruuth, R.J. Spiteri, Implicit–explicit Runge–Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics 25 (1997) 151–167. [2] A. Bensoussan, J.L. Lions, Applications des Ine´quations Variationnelles en Controˆle Stochastique, Dunod, Paris, 1978. [3] M. Boman, A posteriori error analysis in the maximum norm for a penalty finite element method for the time-dependent obstacle problem, Chalmers Finite Element Center, Preprint 2000-12, January 2001. [4] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2002. [5] M.P. Calvo, J. de Frutos, J. Novo, Linearly implicit Runge–Kutta methods for advection–reaction–diffusion equations, Applied Numerical Mathematics 37 (2001) 535–549. [6] R. Cont, E. Voltchkova, Finite difference methods for option pricing in jump-diffusion and exponential Le´vy models, Rapport interne CMAP 513 (September) (2003). [7] G.J. Cooper, A. Sayfy, Additive Runge–Kutta methods for stiff ordinary differential equations, Numerische Mathematik 36 (1981) 431–445. [8] P.A. Forsyth, K.R. Vetzal, Quadratic convergence for valuing American options using a penalty method, SIAM Journal on Scientific Computing 23 (2002) 2095–2122. [9] A. Friedman, Stochastic Differential Equations and Applications, 2 vol., Academic Press, New York, 1975. [10] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer, New York, 1984. [11] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems, Springer, New York, 1991. [12] P. Jaillet, D. Lamberton, B. Lapeyre, Variational inequalities and the pricing of American options, Acta Applicandae Mathematicae 21 (1990) 263–289. [13] R. Kangro, R. Nicolaides, Far field conditions for Black–Scholes equations, SIAM Journal on Numerical Analysis 38 (2000) 1357–1368. [14] C.A. Kennedy, M.H. Carpenter, Additive Runge–Kutta schemes for convection–diffusion–reaction equations, Applied Numerical Mathematics 44 (2003) 131–181. [15] X. Jiang, R.H. Nochetto, Effect of numerical integration for elliptic obstacle problems, Numerische Mathematik 67 (1994) 501– 512. [16] M.D. Marcozzi, On the approximation of optimal stopping problems with applications to financial mathematics, SIAM Journal on Scientific Computing 22 (2001) 1865–1884. [17] R.H. Nochetto, Sharp L1-error estimates for semilinear elliptic problems with free boundaries, Numerische Mathematik 54 (1988) 243–255. [18] B. Øksendal, Stochastic Differential Equations, Springer, New York, 2000. [19] I. Sapariuc, M.D. Marcozzi, J.E. Flaherty, A numerical analysis of variational valuation techniques for derivative securities, preprint, July 2003. [20] R. Scholz, Numerical solution of the obstacle problem by the penalty method. Part II: Time-dependent problems, Numerische Mathematik 49 (1986) 255–268. [21] D. Tavella, C. Randall, Pricing Financial Instruments. The Finite Difference Method, Wiley, New York, 2000.
1004
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
[22] C. Va´zquez, An upwind numerical approach for an American and European option pricing model, Applied Mathematics and Computation 97 (1998) 273–286. [23] P. Wilmott, Derivatives, Wiley, New York, 1998. [24] P. Wilmott, J. Dewynne, S. Howison, Option Pricing, Oxford Financial Press, Oxford, 1993. [25] H. Windcliff, P.A. Forsyth, K.R. Vetzal, Shout options: A framework for pricing contracts which can be modified by the investor, Journal of Computational and Applied Mathematics 134 (2001) 213–241. [26] R. Zvan, P.A. Forsyth, K.R. Vetzal, Penalty methods for American options with stochastic volatility, Journal of Computational and Applied Mathematics 91 (1998) 199–218.
Implicit–explicit Runge–Kutta methods for financial derivatives pricing models Javier de Frutos
*
Departamento de Matema´tica Aplicada, Universidad de Valladolid, C/Prado de la Magdalena S/N, 47005 Valladolid, Spain Available online 10 March 2005
Abstract Implicit–explicit Runge–Kutta methods are investigated for application to financial derivatives pricing models in the partial differential equations approach. The methods are showed to be an alternative to other existing procedures for the numerical valuation of American type contracts. We follow the method of lines in order to have a numerical method that can be used with a variety of state variable discretizations including finite elements, finite differences and finite volume methods. Some numerical experiments are presented. 2005 Elsevier B.V. All rights reserved. Keywords: Finance; American options; Implicit–explicit Runge–Kutta methods; Finite element method
1. Introduction A wide variety of derivative securities traded in exchange markets such as call or put options on dividend paying stocks, foreign currency options, callable bonds, among others, are American type contracts. An American contract gives the holder the right to exercise and exit the contract at any time before the expiration date specified in the contract. The problem of valuating such contracts is considerably more difficult than the corresponding European contract where the holder can only exercise at expiration. In fact there is no closed form analytical formulas for valuing American contracts even in the simpler cases. Numerical methods are then necessary for pricing such contracts. Derivative securities can usually be defined as optimal stopping problems [18] of the form V ðx; tÞ ¼ sup Ex ½gðs; X s Þ;
x 2 X;
s2Tðt;T Þ
*
Tel.: +34 983 42 31 81; fax: +34 983 42 30 13. E-mail address: [email protected]
0377-2217/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2005.01.013
ð1Þ
992
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
associated with a system of stochastic differential equations dX t ¼ bðX t ; tÞ dt þ rðX t ; tÞ dW t ;
ð2Þ
and a reward function g, where X is a possibly unbounded subset of Rn and Tðt; T Þ is the set of stopping times. Here, Ex denotes the mathematical expectation with respect the probability law of the process Xt, starting at X0 = x. In many cases the reward function is chosen to be the discounted value of a given payoff w. Z s gðs; X s Þ ¼ exp rðs; X s Þ ds wðs; X s Þ: t
The value function is representable deterministically as the solution of a system of time dependent variational inequalities [2,9,12] in terms of the characteristic differential operator associated with (2). In this formulation the problem is similar to a parabolic, possibly degenerate, obstacle problem in which the state variables play the role of spatial degrees of freedom. Most practical numerical methods for the solution of the equivalent linear complementarity problem use well known finite differences approaches, [10,23,24]. In [26,8] a penalty method coupled with a finite volume discretization of spatial like variables and finite differences in time are proved to be a useful and very efficient alternative for this kind of problems. The main idea is to change the set of variational inequalities into a set of nonlinear partial differential equations through a penalization of the restriction [2,17,20]. As it is claimed in [8], the main advantage of the penalty method is that it is easily generalizable to multidimensional problems and a variety of contracts. Furthermore, it can be used coupled with any type of spatial (state) discretization, finite differences, finite volume or finite element methods. In this paper we are concerned with implicit–explicit Runge–Kutta methods for the numerical integration of the system of ordinary differential equations which arises after spatial discretization of the penalty approximation to the linear complementarity problem. That is, we use the method of lines to decouple the spatial like and temporal discretization allowing the use of efficient stiffly accurate time integrators others than simple finite differences. We remark that, in the so-called continuation region, one has to numerically solve a system of convection–diffusion–reaction partial differential equations. This problem qualifies as stiff, see [11]. So, stiffly accurate time integrators are necessary to have efficient numerical methods. Note that efficiency is a crucial issue in order to design numerical methods of practical use in general, possibly multidimensional, financial problems. Surprisingly, it seems that there has not been any attempt to use efficient, stiffly accurate time integrators other than the first-order backward Euler method in the financial literature. The only exception is [19], where the authors propose to use a numerical differentiation formula with quasiconstant step size and variable order. Implicit–explicit Runge–Kutta methods are a particular instance of additive Runge–Kutta methods [7]. It has been proved in [5], see also [1,14], that being easily implementable, implicit–explicit Runge–Kutta methods are competitive with standard stiffly accurate time integrators for convection–diffusion–reaction equations. To our knowledge this is the first time that they have been successfully applied to the numerical solution of time dependent variational inequalities. Recently, implicit–explicit finite difference methods has been successfully used for European and barrier options in jump diffusion and exponential Le´vy models in [6]. In the following, for clarity in the exposition, we restrict ourselves to a simple one dimensional model problem, the American put option, although the numerical methods is intended to be general enough to be applied to more complex financial problems. In Section 2 we state the model problem and its penalty approximation. Implicit–explicit Runge–Kutta time integrators coupled with a finite element discretization are presented in Section 3. The next section is devoted to show the qualities of the method proposed by means of some numerical experiments. We finish the paper with some conclusions and final remarks.
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
993
2. The model problem Consider an asset with price S which is supposed to follow a log-normal evolution process dS ¼ lS ds þ rS dW ;
ð3Þ
where l is the expected rate of return, r is the volatility, dW is the increment of a Wiener process and s denotes time. We are concerned with the problem of determining the fair value V(S, s) of an American put option in which the holder can exercise at any time, receiving a payoff specified in the contract as w(S, s) = max(E S, 0), where E is the strike. Under the usual assumptions on the market, the problem can be formally reduced to finding V(S, s) such that V t LBS V P 0;
0 6 S 6 1; 0 < t 6 T ;
V ðS; tÞ P wðS; tÞ;
0 6 S 6 1; 0 < t 6 T ;
ðV t LBS V ÞðV wÞ ¼ 0; V ðS; 0Þ ¼ wðSÞ;
0 6 S 6 1; 0 < t 6 T ;
0 6 S 6 1:
ð4Þ ð5Þ ð6Þ ð7Þ
Here, LBS denotes the well-known Black–Scholes operator, 1 LBS V ¼ r2 S 2 V SS þ rSV S rV : 2
ð8Þ
r is the risk free interest rate and t = T s is the time to expiration of the contract. The problem has a unique solution which is smooth in the continuation region and is continuously differentiable globally [12,2, Chapter 3]. Note that no boundary conditions are imposed as corresponds to a singular problem. From a practical point of view the solution of (4)–(7) can be approximated by means of a sequence of problems posed over bounded increasing domains. More precisely, let us consider a sequence of domains ½S km ; S kM ð0; 1Þ with S km ! 0 and S kM ! 1. Let Vk be the solution of the localized problem S km 6 S 6 S kM ; 0 < t 6 T ;
V kt LBS V k P 0;
S km 6 S 6 S kM ; 0 < t 6 T ;
V k ðS; tÞ P wðS; tÞ;
ðV kt LBS V k ÞðV k wÞ ¼ 0; V k ðS km ; tÞ ¼ E S km ; V k ðS kM ; tÞ ¼ 0;
0 < t 6 T;
0 < t 6 T;
V k ðS; 0Þ ¼ wðSÞ; Let X ¼ ½S l ; S u
S km 6 S 6 S kM ; 0 < t 6 T ;
S km 6 S 6 S kM :
ðS km ; S jM Þ k
ð10Þ ð11Þ ð12Þ ð13Þ ð14Þ
a fixed approximation domain. Then from [16,12,13],
lim max kV ð; tÞ V ð; tÞkL1 ðXÞ ¼ 0:
k!1 t2½0;T
ð9Þ
ð15Þ
We note that the boundary conditions (12) and (13) are essentially arbitrary but nevertheless financially reasonable. Furthermore, in [16] it is shown that every choice of the boundary conditions in an increasing exhausting sequence of bounded domains leading to well posed systems can be used to approximate V(S, t) in a given fixed bounded approximation domain. However, the choice of boundary conditions surely affects the quality of the approximation. The rationale for the choice (12), (13) is as follows: First it is clear that the
994
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
value of the option must be zero when the asset value tends to infinity. The boundary condition in the neighborhood of S = 0 is slightly more subtle. A wrong value will generate a boundary layer near S = 0 that is, obviously, not financially acceptable. From the complementarity condition (6) one has that, at S = 0, either Vt + rV = 0 or V(0, t) = w(0, t) = E. Using (7), in the first case it is found that V(0, t) = exp(rt)E < E for all time t > 0. Then, the constraint (5) forces V(0, t) = E and in consequence (12) for S km small enough. Other boundary conditions have been used in the literature. In [21] linearity type boundary conditions of the form VSS = 0, together with some other useful alternatives, are advocated for options with piecewise linear payoff when finite differences are used as discretization tool. However, this type of boundary conditions are more difficult to implement in the context of a finite element method. 2.1. The penalty problem There are several approaches that can be used to solve the localized problem (9)–(14). Here, we shall use the penalization technique that has been proved to be an efficient way to treat some financial problems, see [8]. To this end, let us start by choosing a smooth function u with uðS km Þ ¼ E S km , uðS kM Þ ¼ 0, and otherwise arbitrary. For > 0, we seek Vk,, the unique solution of the (regularized) penalty problem k; U k; þ b ðU k; Þ ¼ f ; t LBS U
U k; ðS km ; tÞ ¼ U k; ðS kM ; tÞ ¼ 0;
S km 6 S 6 S kM ; 0 < t 6 T ;
ð17Þ
0 < t 6 T;
U k; ðS; 0Þ ¼ w ðSÞ ¼ wðSÞ uðSÞ;
ð16Þ
S km 6 S 6 S kM ;
ð18Þ
where f = ut + LBSu, w*(S) = w(S) u(S) is the modified payoff and b is the penalty term, b ðuðSÞÞ ¼
1 minðuðSÞ w ðSÞ; 0Þ:
ð19Þ
Problem (16)–(18) has a unique solution with derivatives up to second-order in Lp ðS km ; S kM Þ for t > 0, see for example, [9, Theorem 10.4.2]. Furthermore, for small enough, there exists a constant C > 0, such that [3, Theorem 4.1] max kV k ð; tÞ ðU k; ð; tÞ þ uð; tÞÞkL1 ðXÞ 6 C:
ð20Þ
t2½0;T
Besides, the American restriction (10) is also satisfied to OðÞ terms, that is, V k; ðS; tÞ ¼ U k; ðS; tÞ þ uðS; tÞ P wðSÞ C;
S 2 ½S km ; S kM
ð21Þ
for some constant C > 0 independent of .
3. The discrete problem We follow the method of lines and discretize first the asset value and then time. One of the advantages of this methodology is that decoupling the discretization problem into a spatial-like part and a temporal one allows an easier treatment of the specific characteristic of each problem. Here, we discretize the asset value by means of the linear finite element method. In order to simplify somewhat the exposition we take S km ¼ 0 in the rest of the paper. For this particular problem it can be observed that in practice all results of the previous section remain valid.
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
995
Let us consider a partition 0 ¼ S 0 < S 1 < < S J þ1 ¼ S kM of the computational domain, and let us denote h = max06i6J(Si+1 Si). We consider the space of test functions n o X h ¼ / 2 Cð½0; S kM Þj/j½S i ;S iþ1 2 P1 ; /ð0Þ ¼ /ðS kM Þ ¼ 0 ; where P1 denotes the space of polynomials of degree less than or equal to one in the asset variable S. The semidiscrete problem consists in finding a function uh : [0, T] ! Xh such that Z k Z Sk Z Sk Z Sk Z Sk M M M M r2 S M 2 uh;t / dS þ S uh;S /S dS ðr r2 Þ Suh;S / dS þ r uh / dS þ b ðuh Þ/ dS 2 0 0 0 0 0 Z S kM f / dS 8/ 2 X h : ð22Þ ¼ 0
We remark that although other variational formulations are possible (see for example [9,12,16,22]), from a computational point of view this one has the merit of being numerically efficient and directly related to the financial problem. Let /j 2 Xh, j = 1, . . . , J, be defined by /j(Sk) = djk. It is clear that Xh = span{/j, j = 1, . . . , J}. Therefore, writing uh ðS; tÞ ¼
J X
U j ðtÞUj ðSÞ;
j¼1
the variational equations (22) are equivalent to the following system of ordinary differential equations: M
dU h r2 þ KU h ðr r2 ÞCU h þ rMU h þ Mb ðU h Þ ¼ MF h ; dt 2 T
T
ð23Þ T
where U h ¼ ½U 1 ; . . . ; U J , F h ¼ ½f ðS 1 Þ; . . . ; f ðS J Þ and b ðU h Þ ¼ ½b ðU 1 Þ; . . . ; b ðU J Þ . The mass, M, stiffness, K, and convection, C, matrices are defined by Z S M J M¼ /i /j dS ; 0
K¼
Z
i;j¼1 SM
S 2 /i;S /j;S dS
J
0
C¼
Z
J
SM
; i;j¼1
S/i;S /j dS 0
: i;j¼1
For the semidiscrete approximation (22) it is possible to show that for each t > 0, the following error bound holds [4, Chapter 8]: kU k; ð; tÞ uh ð; tÞkL1 ðXÞ 6 Ch2 j log hj; and in consequence, taking (20) into account, kV k ð; tÞ ðuh ð; tÞ þ /ð; tÞÞkL1 ðXÞ 6 Cð þ h2 j log hjÞ:
ð24Þ
3.1. Implicit–explicit time discretizations The system of ordinary differential equations (23) can be solved by means of an standard ODE time integrator. The main difficulty when dealing which this kind of system is that the use of explicit time integrators
996
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
is inefficient because the system becomes stiffer as the asset value mesh is refined. As a consequence the time step must be heavily reduced in order to fulfill the severe stability restriction inherent to explicit methods. On the other hand, if a stiffly accurate, implicit time integrator is chosen, one has to solve implicit nonlinear equations with nonsmooth Jacobian. Furthermore, the convective terms introduce additional numerical difficulties. This could be one of the reasons why very few implicit codes other than mid-point rule can be found in the literature. To our knowledge perhaps the only exception is [19] where a backward differentiation type code is proposed. However, multistep based methods are reported in [25], see also [8], to perform badly with American type financial problems as shout options. Although it is not clear in the financial literature if this is also the case for simpler contracts, the use of one step time integrators present some practical advantages. Among others, there exist easily implementable methods with variable time steps and automatic error control that are efficient for convection–diffusion type equations that typically arise in financial problems. The main idea in an implicit–explicit time discretization is to use different methods for the different terms rolled into a single composite method that could treat more efficiently each one of the terms appearing in (23). To this end let us write (23) in the form M
dU h ¼ LðU h Þ þ NðU h Þ; dt
ð25Þ
where r2 KU h Mb ðU h Þ; 2 NðU h Þ ¼ ðr r2 ÞCU h rMU h þ MF h : LðU h Þ ¼
The partition (25) has been chosen such that L(Uh) collects the stiffer terms appearing in the equations (diffusive and penalty terms in our case) and N(Uh) gathers the less stiff terms (convective and source terms). The simplest implicit–explicit Runge–Kutta method is a combination of the well known explicit and implicit Euler rules. Let us suppose that at time tn < T a numerical approximation Unh U h ðtn Þ has been found. Then the approximation U hnþ1 to Uh(tn+1) = Uh(tn + Dt) is the solution of the system of algebraic equations n MU hnþ1 ¼ MU nh þ DtLðU nþ1 h Þ þ DtNðU h Þ:
ð26Þ
Note that as the most stiff terms, diffusion and penalty ones, are treated implicitly by means of a L-stable method (Backward Euler method in this simple case), only a mild stability restriction Dt 6 C independent of asset grid size has to be imposed, see [5]. This is the main advantage of implicit–explicit methods over pure explicit approaches. The price to be paid is that at each time step, a nonlinear system of equations has to be solved. However, as the convection term in (25) is treated explicitly, only a linear solver for symmetric systems is needed when (26) is solved by means of a Newton like iteration procedure. This represents usually an advantage over fully implicit methods such as multistep BDF formulas. For the numerical solution of (26), following [8], we propose to use the following fixed point Newton-like iteration: Given an initial iterate Y0, for example, Y 0 ¼ U nh , for k = 0,1, . . . , M þ DtK þ DtMBk Y kþ1 ¼ MU nh þ DtNðU nh Þ þ DtMBk w ; ð27Þ T
where w ¼ ½w ðS 1 Þ; . . . ; w ðS J Þ and Bk is the diagonal matrix defined by ( 1= if ðY ki Þj < w ðS j Þ; k Bjj ¼ 0 if ðY ki Þj P w ðS j Þ: The iteration is stopped when the difference between two consecutive iterates is small enough. We refer to Section 4 for the practical details concerning the implementation of the nonlinear iteration (27).
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
997
We note that (27) is, formally, the Newton method applied to the nonlinear system (26) in which matrix B plays the role of the Jacobian of the (nonsmooth) nonlinear penalty term. In [8] it is proved that if a lumped (diagonal) mass matrix [15] is used in (27), the iterates converge monotonically to the unique solution of (26). In our experiments we have found that with consistent mass matrix, (27) still converges in a small number of iterations. 3.2. Implicit–explicit Runge–Kutta methods The implicit–explicit Euler method is a first-order method, that is, the error decreases proportionally to time step Dt and in consequence is somewhat inefficient. However, high order one step methods sharing the same ideas that the implicit–explicit Euler methods are available in the literature. In an implicit–explicit Runge–Kutta method, the equations that describe the step tn # tn+1 when applied to (25) take the form MY 1 ¼ MU nh ; MY i ¼
MU nh
ð28Þ
þ Dtaii LðY i Þ þ Dt
i1 X
aij LðY j Þ þ
j¼2
MU hnþ1 ¼ MU nh þ Dt
sþ1 X
bi LðY i Þ þ
i¼2
i1 X
! ^ai;j NðY j Þ ;
2 6 i 6 s þ 1;
ð29Þ
j¼1 sþ1 X
! ^ bi NðY i Þ ;
ð30Þ
i¼1
b ¼ ð^aij Þ where Yi denotes the internal stages. The coefficients of the method, that is, the matrices A = (aij), A T and the vectors bT = (bi) and ^ b ¼ ð^ bi Þ are computed as to satisfy the so-called order conditions (to a certain attainable order depending on the number of stages s + 1) and some stability restrictions. Usually, the coefficients are chosen such as the following simplifying assumptions are satisfied: ci ¼
i X j¼2
aij ¼
i1 X
^ aij ;
bi ¼ ^ bi :
j¼2
We refer to [5], see also [11], for a review of the order conditions and a study of the stability properties of implicit–explicit Runge–Kutta methods. Several efficient implicit–explicit Runge–Kutta methods exist in the literature. We refer for example to [1,14], where an exhaustive comparison among different possibilities are reported. Each stage of the method is an implicit–explicit Euler like step, see Eq. (29), that is treated as has been previously described by means of an iteration similar to (27). Note that only an efficient procedure to solve symmetric positive definite linear systems is needed. Furthermore, the method is easily generalizable to other models were a special stabilization procedure of the convective terms could be needed. We do not pursue further this issue in this paper as this is not the case for the simple model of the American put option. We remark that the method (28)–(30) can be seen as pair of Runge–Kutta methods specified by the T b ^ coefficients A, bT, c and A, b , c, respectively. The first method is implicit and has to be applied to the stiff terms, in this case the diffusive and penalty terms; the second one is explicit and is applied only to the convective and source terms. Moreover, if the implicit method is L-stable [11] it is possible to construct an explicit companion method in such a way that the full implicit–explicit procedure needs only a mild stability restriction of the type Dt 6 C with C independent of the mesh size, see [5]. This is a necessary condition in order to have an efficient numerical procedures for partial differential equations. Besides, it is also a very convenient property when one is dealing with nonsmooth initial data as (18), avoiding the need of special
998
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
starting procedures required when one uses a non strongly stable procedure as for example the implicit midpoint rule (Crank–Nicolson method) [8]. Finally we remark that since the asset and time discretizations are decoupled, the same time integrator can be used in conjunction with other spatial like discretization such as finite differences, finite volume methods or even high order finite elements. 4. Numerical experiments In this section we consider the third-order method (LIRK3) with s + 1 = 4 stages constructed in [5]. The implicit part of the method is based in the well known L-stable, three stage, third-order SDIRK [11]. The coefficients of the explicit part of the method are computed in such a way that the method have quasioptimal stability properties when applied to convection–diffusion equations, see [5, Section 4], like the Black–Scholes equation. The method is implemented in variable-step mode, so that the time step is automatically chosen by the code depending on an estimation of the local errors. To this end, we have considered the second-order formula: e nþ1 ¼ MU n þ Dt b2 ðLðY 2 Þ þ NðY 2 ÞÞ þ b3 ðLðY 3 Þ þ NðY 3 ÞÞ þ b4 ðLð U e nþ1 Þ þ NðY 4 ÞÞ ; MU h h h where Y2, Y3 and Y4 are the internal stages of the third-order scheme. The second-order approximation e nþ1 can be interpreted as the solution generated with a second-order implicit–explicit Runge–Kutta U h e nþ1 defined by the coefficients a5;i ¼ ^a5;i ¼ bi , 1 6 i 6 3, a5;4 ¼ method with an additional stage Y 5 ¼ U h ^ a5;5 ¼ 0, a5;5 ¼ ^ a5;4 ¼ b4 . Each time a new approximation U nþ1 is computed, the local error is estimated by h e nþ1 k : ERR ¼ kU hnþ1 U 1 h The step is accepted if ERR 6 TOL, where TOL is the precision required by the user. If ERR > TOL the step is rejected and the computation is restarted from the previously accepted approximation using a smaller time step. In both cases the new time step size is selected according to the rule 1=3 !! TOL Dtnew ¼ min hlim; Dtold min FACMAX; FAC ; ð31Þ ERR where ‘‘hlim’’ is the maximum admissible step size (in our code hlim = T/5) and FAC and FACMAX are safety factors that limit the variation of step size in order to make the code more robust. In our code we choose FAC ¼ 0:9;
FACMAX ¼ 1:5
if the step has been accepted (EST 6 TOL), and FAC ¼ 0:7;
FACMAX ¼ 1
after a rejection (ERR > TOL). Each stage, a nonlinear problem of the form (26) has to be solved. To this end, we have used the Newton like iteration (27), stopping the iteration when max
jY kþ1 Y kj j j maxð1; jY kj jÞ
6 103 TOL;
where Y kj , j = 1, . . . , denote the components of the iterant Yk in (27).
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
999
In order to control the cost of the iterations in the solution of the nonlinear systems, we have also decided to adopt a conservative policy and limit the maximum number of iterations in each stage. If, in the computation of one stage, the nonlinear iteration has not converged after four iterations, the full time step is rejected and the computation restarted with time step halved. We have tested the method for several values of the parameters. Here, we present the results obtained for the particular values r = 0.2, r = 0.1 and E = 100. The conclusions for other values were similar to the ones reported here. In all the numerical experiments we have taken S kM ¼ 200. For this computational domain the error caused by the arbitrary boundary conditions was found to be small enough compared with the errors made in the discretization process. Fig. 1 shows the time step history for a single run using a grid with 160 nodes and setting TOL = 103 and = 104. In this experiment time to expiry was T = 1. As we can see the code selects small time steps at the beginning of the computation and smoothly increments the time step as the computation proceeds. We remark that the initial condition (18) has only piecewise smooth derivatives although is smoother for t > 0 due to the parabolic character of the Black–Scholes equation. The small steps are the response of the code to the lack of regularity at t = 0. Once the solution has been smoothed, the code increments the time step until a default in the convergence of the nonlinear iterations (27) is detected. Then, the code automatically reduces the time step, see Fig. 1 near T = 0.4, and continues the time integration increasing again the time step according to the rule (31). Note that no special procedure is needed to deal with nonsmooth initial conditions as is the case with other time integrators of common use as, for example, mid-point rule, see [8]. The time step history was found to depend only on the time integrator tolerance and is essentially independent of the asset grid size, see below. In the next test we have checked the behaviour of the method with respect to the penalty parameter. We have monitored the maximum relative error in enforcing the American constraint by the penalty method by means of the quantity AMERR ¼ max
S2½0;S kM
maxðw ðSÞ uh ðSÞ; 0Þ : maxð1; w ðSÞÞ
ð32Þ
Table 1 shows the results obtained for several values of the discretization parameters, number of nodes and tolerance in the step control procedure of the time integrator, varying the size of the penalization parameter . We observe that the size of the relative American error is proportional to , in good agreement
0.12
0.1
0.08
0.06
0.04
0.02
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 1. Time step selection history. r = 0.2, r = 0.1, K = 100, T = 1, grid with 160 nodes, TOL = 103, = 104.
1000
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
Table 1 American put option Nodes
2
10 103 104 105
TOL
American error
Ratio
Steps
· · · ·
3
10 104 104 105
www 8.00 7.62 11.07
17 17 17 27
7.57 8.70 1.04 9.64
· · · ·
103 104 104 106
www 8.70 8.36 10.79
25 25 30 32
7.56 8.63 9.42 1.25
· · · ·
103 104 105 105
www 8.76 9.16 7.54
42 46 50 52
160 160 160 160
2
10 102 102 102
7.56 9.45 1.24 1.12
102 103 104 105
320 320 320 320
103 103 104 105
102 103 104 105
640 640 640 640
104 104 104 104
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, 640 nodes, = 106.
with bound (21). Furthermore, the error is nearly independent of the discretization parameters. It is interesting to see that the number of time steps and, as a consequence, the cost of time integration is nearly independent of the penalization parameter. This is a consequence of the strong stability properties of the time integrator. The slight increment in the number of time steps observed with decreasing values of is due to the conservative policy adopted for time step control. Each time a slow convergence of the iterates in (27) is detected, the time step is halved. The number of time step increases because of this somewhat drastic reduction in the time step. In order to check the accuracy of the implicit–explicit Runge–Kutta method in this application we performed several runs for different values of the parameter TOL. The number of nodes was set to be 640 in all the reported experiments so that the numerical solutions are nearly free of spatial discretization error. The results of the experiment are shown in Table 2. The most precise result corresponds to TOL = 106 and has
Table 2 American put option TOL
Steps
2
10 103 104 105 106
V (1 0 0)
18 37 50 91 153
Difference 2.13 · 2.40 · 3.01 · 2.00 · www
3.04907 3.06794 3.07003 3.07032 3.07034
Ratio
2
www 8.88 7.97 15.05 www
10 103 104 105
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, 640 nodes, = 106.
Table 3 American put option Nodes
TOL
Steps
V (1 0 0)
Difference
40 80 160 320
102 102 103 103
17 17 30 30
3.10001 3.07736 3.07262 3.07097
2.965 7.019 2.274 5.267
Computational domain [0, 200]. T = 0.25, r = 0.1, r = 0.2, E = 100, TOL = 106, = 106.
· · · ·
102 103 103 104
Ratio www 4.2 3.1 4.3
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
1001
50
45
40
35
30
25
20
15
10
5
0 50
60
70
80
90
100
110
120
130
140
150
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
-0.7
-0.8
-0.9
-1 50
60
70
80
90
100
110
120
130
140
150
90
100
110
120
130
140
150
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
-0.01 50
60
70
80
Fig. 2. Top: Option value (V), middle: delta (VS), bottom: gamma (VSS). Grid with 160 nodes, TOL = 103, = 104, T = 0.25.
been taken as a reference solution in order to measure the efficiency of the numerical method. For each value of the tolerance we show in Table 2 the number of steps, the computed price of the option at
1002
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
S = 100 and the difference between computed prices and the reference solution. This difference is an indicator of the error due to the temporal discretization. Last column shows the ratio of consecutive differences. We observe that the error is of the same order as the prescribed tolerance. We remark that in a step size control procedure, the objective is to have an error in the numerical solution as close as possible to the user prescribed tolerance. In this way, the user is able to optimize, for a given precision, the number of time steps and in consequence the CPU time required by the computation. The high values in the last column of Table 2 nicely show a very rapid convergence as the tolerance is reduced. Note also that each time the time integrator error is reduced by an order of magnitude, the number of steps is increased by a factor less than two. So, in spite of the some times alleged lack of regularity of the American put problem, by means of a stiffly accurate, third-order time integrator, the error can be drastically reduced with only a moderate increment in the cost of the computation. It is clear that, in a practical computation, an equilibration of the errors coming from asset and temporal discretizations is needed. We present in Table 3 the results obtained by minimizing the cost of the computation for a given error. The value of the time integrator tolerance showed in each row is the best value in the sense that a smaller value does not reduce the error while it increases the number of steps and in consequence the cost of the computation. In all runs the value of the penalization parameter was taken to be = TOL/10. As we have done before we take as an indicator of the error in the computation the difference with respect to the numerical solution computed with 640 nodes, TOL = 106 and = 107. As we can see, the ratio between two consecutive differences is very near 4 as corresponds to a second-order method. Furthermore, thanks to the high order temporal discretization, the error in the asset discretization can be reduced while maintaining a small number of time steps. To sum up, the numerical experiments show that the error in the proposed method decouples in three components, the error due to the penalty term, the temporal discretization error and the asset discretization error which can be controlled independently. The high order strongly stable method used for the time integration makes the global procedure very efficient providing rapid convergence with respect to the time integration tolerance and quadratic convergence with respect to the number of degrees of freedom used in the asset discretization. The sensitivity of the option value to variations in asset price and parameters of the model is measured by the so called Greeks [21,23]. In practice, besides the option price, accurate approximations to the Greeks are needed for edging purposes [23]. Fig. 2 shows value, delta (oV/oS) and gamma (o2V/oS2) of the option for the same values of the parameters as before and time to expiry T = 0.25. We can observe that the numerically calculated Greeks are free of spurious oscillations which are typically present in the numerical solution of American type contracts, specially in the computed gamma, due to the discontinuity of the second derivative in the optimal exercise price. This unwanted behaviour of the gamma has been reported in [8] where a second-order mid-point rule (Crank–Nicolson method) is used as time integrator. There, the authors claim that the oscillations disappear if the first-order, fully implicit, Backward Euler method is used instead. Note that the implicit–explicit Runge–Kutta method is a strongly stable method that shares the dissipative properties of implicit Euler method while it retains the efficiency derived from the third-order accuracy and variable step implementation.
5. Conclusions We have presented a third-order variable step implicit–explicit Runge–Kutta method which is showed to be an efficient alternative to other existing procedures for the numerical valuation of American put options. The method is coupled with a linear finite element discretization for the asset variable and treats the American restriction by the penalty method. The computed Greeks are free of the spurious oscillations which are
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
1003
typically present when no strongly stable time integrators, such as Crank–Nicolson method, are used, while the efficiency inherent to high order methods is retained. Using the method of lines to state the full discretization it is possible to completely decouple asset discretization from time integration. This allows the time integrator to be used with other asset discretizations such as the finite difference or finite volume method. We remark that, for the same reason, the extension to multifactor problems is also straightforward.
Acknowledgments Research supported by DGI-MCYT under grant BFM2001-2138 (cofinanced by FEDER funds) and by JCYL under project VA044/03. This research was completed while the author was visiting professor at GERAD and HEC Montre´al partially supported by grant SEEU PR2002-0223.
References [1] U.M. Ascher, S.J. Ruuth, R.J. Spiteri, Implicit–explicit Runge–Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics 25 (1997) 151–167. [2] A. Bensoussan, J.L. Lions, Applications des Ine´quations Variationnelles en Controˆle Stochastique, Dunod, Paris, 1978. [3] M. Boman, A posteriori error analysis in the maximum norm for a penalty finite element method for the time-dependent obstacle problem, Chalmers Finite Element Center, Preprint 2000-12, January 2001. [4] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2002. [5] M.P. Calvo, J. de Frutos, J. Novo, Linearly implicit Runge–Kutta methods for advection–reaction–diffusion equations, Applied Numerical Mathematics 37 (2001) 535–549. [6] R. Cont, E. Voltchkova, Finite difference methods for option pricing in jump-diffusion and exponential Le´vy models, Rapport interne CMAP 513 (September) (2003). [7] G.J. Cooper, A. Sayfy, Additive Runge–Kutta methods for stiff ordinary differential equations, Numerische Mathematik 36 (1981) 431–445. [8] P.A. Forsyth, K.R. Vetzal, Quadratic convergence for valuing American options using a penalty method, SIAM Journal on Scientific Computing 23 (2002) 2095–2122. [9] A. Friedman, Stochastic Differential Equations and Applications, 2 vol., Academic Press, New York, 1975. [10] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer, New York, 1984. [11] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems, Springer, New York, 1991. [12] P. Jaillet, D. Lamberton, B. Lapeyre, Variational inequalities and the pricing of American options, Acta Applicandae Mathematicae 21 (1990) 263–289. [13] R. Kangro, R. Nicolaides, Far field conditions for Black–Scholes equations, SIAM Journal on Numerical Analysis 38 (2000) 1357–1368. [14] C.A. Kennedy, M.H. Carpenter, Additive Runge–Kutta schemes for convection–diffusion–reaction equations, Applied Numerical Mathematics 44 (2003) 131–181. [15] X. Jiang, R.H. Nochetto, Effect of numerical integration for elliptic obstacle problems, Numerische Mathematik 67 (1994) 501– 512. [16] M.D. Marcozzi, On the approximation of optimal stopping problems with applications to financial mathematics, SIAM Journal on Scientific Computing 22 (2001) 1865–1884. [17] R.H. Nochetto, Sharp L1-error estimates for semilinear elliptic problems with free boundaries, Numerische Mathematik 54 (1988) 243–255. [18] B. Øksendal, Stochastic Differential Equations, Springer, New York, 2000. [19] I. Sapariuc, M.D. Marcozzi, J.E. Flaherty, A numerical analysis of variational valuation techniques for derivative securities, preprint, July 2003. [20] R. Scholz, Numerical solution of the obstacle problem by the penalty method. Part II: Time-dependent problems, Numerische Mathematik 49 (1986) 255–268. [21] D. Tavella, C. Randall, Pricing Financial Instruments. The Finite Difference Method, Wiley, New York, 2000.
1004
J. de Frutos / European Journal of Operational Research 171 (2006) 991–1004
[22] C. Va´zquez, An upwind numerical approach for an American and European option pricing model, Applied Mathematics and Computation 97 (1998) 273–286. [23] P. Wilmott, Derivatives, Wiley, New York, 1998. [24] P. Wilmott, J. Dewynne, S. Howison, Option Pricing, Oxford Financial Press, Oxford, 1993. [25] H. Windcliff, P.A. Forsyth, K.R. Vetzal, Shout options: A framework for pricing contracts which can be modified by the investor, Journal of Computational and Applied Mathematics 134 (2001) 213–241. [26] R. Zvan, P.A. Forsyth, K.R. Vetzal, Penalty methods for American options with stochastic volatility, Journal of Computational and Applied Mathematics 91 (1998) 199–218.