Iterative Methods for Pricing American Options under the Bates Model

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Iterative Methods for Pricing American Options under the Bates Model Salmi, Santtu; Toivanen, Jari; von Sydow, Lina Salmi, S., Toivanen, J., & von Sydow, L. (2013). Iterative Methods for Pricing American Options under the Bates Model. Procedia Computer Science, 18, 11361144. https://doi.org/10.1016/j.procs.2013.05.279 2013 Procedia Computer Science 18 ( 2013 ) 1136 – 1144 1877-0509 © 2013 The Authors. Published by Elsevier B.V. Selection and peer review under responsibility of the organizers of the 2013 International Conference on Computational Science doi: 10.1016/j.procs.2013.05.279 International Conference on Computational Science, ICCS 2013 Iterative Methods for Pricing American Options under the Bates Model Santtu Salmia, Jari Toivanena,b,∗, Lina von Sydowc aDepartment of Mathematical Information Technology, P.O. Box 35 (Agora), FI-40014 University of Jyv¨askyl¨a, Finland bInstitute for Computational and Mathematical Engineering, Stanford University, Stanford, CA 94305, USA cDepartment of Information Technology, Uppsala University, Box 337, 751 05 Uppsala, Sweden Abstract We consider the numerical pricing of American options under the Bates model which adds log-normally distributed jumps for the asset value to the Heston stochastic volatility model. A linear complementarity problem (LCP) is formulated where partial derivatives are discretized using finite differences and the integral resulting from the jumps is evaluated using simple quadrature. A rapidly converging fixed point iteration is described for the LCP, where each iterate requires the solution of an LCP. These are easily solved using a projected algebraic multigrid (PAMG) method. The numerical experiments demonstrate the efficiency of the proposed approach. Furthermore, they show that the PAMG method leads to better scalability than the projected SOR (PSOR) method when the discretization is refined. Keywords: American option; Bates model; Finite difference method; Iterative method; Linear complementarity problem 1. Introduction In this paper we consider the numerical pricing of American options. Such options can be exercised prior to the date of maturity which leads to a free-boundary problem. This is in contrast to European options that can only be exercised on the date of maturity leading to an easier problem to solve. Since trading of options has grown to a tremendous scale during the last decades the need for accurate and effective numerical option pricing methods is obvious. The most common options give the holder either the right to sell (put option) or buy (call option) the underlying asset for the strike price. A mathematical model to describe the behavior of the underlying asset is needed to compute the option price. Many such models of varying complexity exist. Typically, more complicated models reproduce more realistic paths of the underlying asset and are hence better to give accurate option prices but they also make the numerical pricing process more challenging. The most commonly used model is the BlackScholes model [1], which assumes the value of the underlying asset to follow a geometric Brownian motion. In the Merton model [2] log-normally distributed jumps are added to the Black-Scholes model, while in the Kou model [3] the jumps are log-doubly-exponentially distributed. By making the volatility a stochastic quantity the Heston model is derived [4], while the Bates model [5] combines the Merton model with the Heston model by adding log-normally distributed jumps to the latter one. Finally, the correlated jump model [6] also lets the volatility jump in the Bates model. ∗Corresponding author. E-mail address: toiv[email protected]. Available online at www.sciencedirect.com Open access under CC BY-NC-ND license. 1137 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 One way to price options is to employ a Monte-Carlo type solver that simulates the behavior of the underlying asset using the model employed and then compute discounted mean values. Such methods are known to have nonfavorable convergence properties and the treatment of the early exercise feature is nontrivial which is why we use another approach here. We formulate a linear complementarity problem (LCP) for a partial (integro-) differential equation (P(I)DE) operator for the price, discretize the P(I)DE, and then solve the resulting LCPs. Several methods have been proposed for solving the resulting LCPs. The Brennan and Schwartz algorithm [7] is a direct method for pricing American options under the Black-Scholes model; see also [8]. Numerical methods for pricing under the Heston model have been developed in [9], [10], [11], [12], [13], [14], for example. The treatment of the jumps in the Merton and Kou models have been studied in [15], [16], [17], [18], [19], [20], for example. Pricing under the Bates model has been considered in [21], [22] and under the correlated jump model in [23]. In this paper, we price American call options under the Bates model. The spatial partial derivatives in the resulting partial integro-differential operator are discretized using a seven-point finite difference stencil and the integral term is discretized using a simple quadrature rule. The Rannacher scheme [24] is employed for the time stepping. We solve the resulting LCPs by employing a fixed point iteration described and analyzed in [25] where each iteration requires the solution of an LCP. These are solved using a projected multigrid (PAMG) method which was recently introduced in [26]. The numerical experiments demonstrate that the proposed method is orders of magnitude faster than the projected successive overrelaxation (PSOR) method. The outline of the paper is the following. The Bates model and an LCP formulation for an American call option is described in Sect. 2. In Sect. 3 the discretization of the LCPs is introduced and the iterative method to solve them is proposed in Sect. 4. Numerical experiments are presented in Sect. 5 and conclusions are given in Sect. 6. 2. Option Pricing Model Here, we consider the Bates model [5] that combines the Merton jump model [2] and the Heston stochastic volatility model [4]. It describes the behavior of the asset value sand its variance yby the coupled stochastic differential equations ds =(μ−λξ)sdt +√ysdw1+(J−1)sdn, dy =κ(θ−y)dt +σ√ydw2. Here μis the growth rate of the asset value, κis the rate of reversion to the mean level of y,θis the mean level of y, and σis the volatility of the variance y. The two Wiener processes w1are w2have the correlation ρ. The Poisson arrival process nhas the rate λand the jump size Jis taken from a distribution f(J)=1 √2πδJexp −[ln J−(γ−δ2/2)]2 2δ2, where γand δdefine the mean and variance of the jump. The mean jump ξis given by ξ=exp(γ)−1. For simplicity, from now on we assume that the market prices of the volatility and jump risks are zero. Applying the Feynman-Kac formula to the Bates model we arrive at the following PIDE 0=∂u ∂τ −1 2ys2∂2u ∂s2−ρσys ∂2u ∂s∂y−1 2σ2y∂2u ∂y2−(r−q−λξ)s∂u ∂s−κ(θ−y)∂u ∂y+(r+λ)u−λ∞ 0u(Js,y,τ)f(J)dJ =∂u ∂τ −a11 ∂2u ∂s2−a12 ∂2u ∂s∂y−a22 ∂2u ∂y2−a1∂u ∂s−a2∂u ∂y+(r+λ)u−λ∞ 0u(Js,y,τ)f(J)dJ =:Lu, (1) where uis the price of a European option, τ=T−tis the time to expiry and qis the dividend yield. The initial condition for (1) is defined by u=g(s,y), where gis the payofffunction which gives the value of option at the expiry. In the following, we consider only call options. A similar approach can be also applied for put options. The payofffunction for a call option with the strike price Kis g(s,y)=max(s−K,0). 1138 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 For the computations, the unbounded domain is truncated to (s,y,τ)∈(0,S)×(0,Y)×(0,T] (2) with sufficiently large Sand Y. The price uof an American option under the Bates model satisfies an LCP ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Lu ≥0, u≥g, (Lu)( u−g)=0. (3) We impose the boundary conditions u(0,y,τ)=g(0,y), u(S,y,τ)=g(S,y),y∈(0,Y), uy(s,Y,τ)=0,s∈(0,S). Beyond the boundary s=S, the price uis approximated to be the same as the payoffg, that is, u(s,y,τ)=g(s,y) for s≥S. On the boundary y=0, the LCP (3) holds and no additional boundary condition needs to be posed. 3. Discretization We will compute approximate prices uon a space-time grid defined by the grid points (xi,yj,τ k), 0 ≤i≤m, 0≤j≤n,0≤k≤l. In space we use a uniform grid with grid steps Δs=S/min the s-direction and Δy=Y/nin the y-direction. We start by introducing a semidiscrete approximation for u(si,yj,τ),0≤i≤m,0≤j≤n. For the non cross-derivatives in (1) we use standard second-order, centered finite difference approximations. In this paper, we assume that the correlation ρis negative. Due to the cross-derivative, we use a seven-point finite difference stencil. A similar stencil has been described in [21], [22]. For a positive correlation ρ, a suitable seven-point stencil is given in [10], [11]. The cross-derivative is approximated by ∂2u ∂s∂y(si,yj,τ)≈1 2ΔsΔy2u(si,yj,τ)−u(si−1,yj+1,τ)−u(si+1,yj−1,τ) −(Δs)2∂2u ∂s2(si,yj,τ)−(Δy)2∂2u ∂y2(xi,yj,τ). (4) Due to the additional derivative terms in (4), we define modified coefficients for ∂2u ∂s2and ∂2u ∂y2as ˜a11 =a11 +1 2 Δs Δya12,and ˜a22 =a22 +1 2 Δy Δsa12. To avoid positive weights in the computational stencil when the convection dominates the diffusion, we add some artificial diffusion according to ˆa11 =min ˜a11,−1 2b1Δs,1 2b1Δsand ˆa22 =min ˜a22,−1 2b2Δy,1 2b2Δy. This is equivalent to using a combination of one-sided and central differences for the convection part. The resulting matrix is an M-matrix with nonpositive off-diagonals and positive diagonal. It is strictly diagonally dominant when r+λ>0. 1139 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 The integral term in (1) at each grid point siis denoted by Ii. We start by making the change of variable J=ez, to obtain Ii=∞ 0 u(Jsi,y,τ)f(J)dJ =∞ −∞ u(ezxi,y,τ)p(z)dz, where pis the probability density function of the normal distribution with mean γ−δ2/2 and variance δ2given by p(z)=1 √2πδ exp −[z−(γ−δ2/2)]2 2δ2. Then we decompose Iiinto one integral over the computational domain defined in (2) and one integral over the remainder of the interval. The first part is then divided on the spatial grid so that we get Ii= n−1  j=0 Ii,j+∞ ln sn−ln si g(ezsi,y)p(z)dz,(5) where Ii,j=ln sj−ln si ln sj+1−ln si u(ezsi,y,τ)p(z)dz.(6) The price function u(s,y,τ) needs to be approximated between each grid point pair (si,si+1) in order to define approximate values of Ii,j. For this, we use a piecewise linear interpolation u(s,y,τ)≈si+1−s si+1−si u(si,y,τ)+s−si si+1−si u(si+1,y,τ) (7) for s∈[si,si+1]. Using (7) in (6) we get Ii,j≈eγ 2erf si,j+1−δ2/2 δ√2−erf si,j−δ2/2 δ√2αjxi+1 2erf si,j+1+δ2/2 δ√2−erf si,j+δ2/2 δ√2βjxi, where erf(·) is the error function, si,j=ln sj−ln si−γ, αj=u(sj+1,y,τ)−u(sj,y,τ) sj+1−sj ,and βj=u(sj,y,τ)sj+1−u(sj+1,y,τ)sj sj+1−sj . The spatial discretization leads to a semi-discrete LCP ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ uτ+Au +a≥0, u≥g, (uτ+Au +a)T(u−g)=0, where Ais an (m+1)(n+1) ×(m+1)(n+1) matrix, ais a vector resulting from the second term in (5), uand gare vectors containing the grid point values of the price uand the payoffg, respectively. In the above LCP, the inequalities hold componentwise. For the temporal discretization we use the Rannacher scheme [24]; see also [27]. The first four time steps are performed with the implicit Euler method with the time step Δτ/2, and then the rest of the time steps are performed with the Crank-Nicolson method with the time step Δτ, where Δτ=T/(l−2). Thus, the time grid is defined by τk=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ k 2(l−2) T,k=0,1,2,3, k−2 l−2T,k=4,5,...,l. The purpose of a few Euler steps in the beginning of the time-stepping process is to damp oscillatory components of the solution. The discretization in time leads to the solution of the following sequence of LCPs: LCP(B(k+1),u(k+1),b(k+1),g),(8) 1140 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 where u(k)denotes the vector uat the kth time step. Here LCP(B,u,b,g) denotes the linear complementarity problem ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ (Bu −b)≥0, u≥g, (Bu −b)T(u−g)=0. For the first four time steps k=0,1,2,3, the implicit Euler method leads to B(k+1) =I+1 2ΔτAand b(k+1) =u(k)−1 2Δτa.(9) For the rest of the time steps k=4,5,...,l−1, the Crank-Nicolson method leads to B(k+1) =I+1 2ΔτAand b(k+1) =I−1 2ΔτAu(k)−Δτa.(10) 4. The solution of LCPs The projected SOR method (PSOR) for LCPs was introduced by Cryer in [30]. The method performs successive over relaxed corrections for the components of the solution vector combined with a projection when a component violates the early exercise constraint. For pricing American options this methods has been discussed in the books [31], [29], for example. The method is fairly simple to implement, but typically the number of iterations grows substantially when the discretization is refined. Thus, it is not usually efficient when fairly accurate option prices are sought. In this paper we will employ PSOR to (8) for comparison, see Sect. 5. Here we will focus on an iterative scheme introduced in [22] and [25]. Let Bdenote the matrix B(k+1) in (9) or (10) associated with the LCP (8). It has a regular splitting [28] B=T−J, where Jis a block diagonal matrix with full diagonal blocks resulting from the integral term and Tis the rest which is a block tridiagonal matrix. Based on this splitting, the first two authors of this paper proposed a fixed point iteration for LCPs in [22], [25]. It is a generalization of an iteration for linear systems described in [29] and applied in [15], [18], [20]. The fixed point iteration for LCP(B,u,b,g) reads LCP(T,uj+1,b+Juj,g),j=0,1,... (11) Each iteration requires the solution of an LCP with the block tridiagonal Tand the multiplication of a vector by J. Below we will describe and compare PSOR and PAMG to solve these LCPs. Based on a convergence result in [25] and the properties of the discretization, we can easily see that the reduction of the l∞-norm of the error in each iteration of (11) is proportional to Δτλ. In general Δτλ is much less than one yielding that the iteration converges very rapidly. In practice, typically only a few iterations are needed to reach sufficient accuracy for practical purposes. One way to solve the LCPs in (11) is to use PSOR. Since it is not an efficient method for refined discretizations we will only use it here for comparison and instead make use of a projected algebraic multigrid method (PAMG) introduced in [26]. With a well designed multigrid method, the number of iterations does not grow with refined discretizations. For extensive literature on this see, the book [32], for example. For solving LCPs Brandt and Cryer introduced a projected full approximation scheme (PFAS) multigrid method in [33]. Another multigrid method for similar problems was described in [34]. The PFAS method was used to price American options under stochastic volatility by Clarke and Parrott in [9], and Oosterlee in [13]. Some alternative approaches employing multigrid methods for option pricing have been considered in [35], [36], [37]. Reisinger and Wittum described a projected multigrid (PMG) method for LCPs which resembles more closely a classical multigrid method for linear problems in [38]. This method has been used to price American options in [38], [11]. The above mentioned methods are so-called geometrical multigrid methods which means that the spatial operators are discretized on sequence of grids. Furthermore, transfer operators between grids need to be implemented. 1141 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 The geometrical multigrid method can be implemented with some effort especially when the computational domain is a rectangle like in this case, but it is not a black-box method to which one can just give the matrix and vectors defining the LCP. An algebraic multigrid (AMG) method [39], [40] builds the coarse problems and the transfer operators automatically using the properties of the matrix. Recently, Toivanen and Oosterlee generalized an AMG method for LCPs and called the resulting method as the projected algebraic multigrid (PAMG) method [26]. Its treatment of LCPs in the coarser levels resemble the one in the PMG method [38]. The PAMG method is easy to use and efficient [26]. Below we present the algorithms for one iteration of PSOR and PAMG respectively. Algorithm One iteration of PSOR(B,u,b,g)Algorithm One iteration of PAMG(B,u,b,g) for i=1,...,dimBif coarsest level then ri=bi−dimB j=1Bijujsolve LCP(B,u,b,g) ui=ui+ωri/Bii else ui=max(ui,gi) PS(B,u,b,g) end for uc=0 rc=˜ R(b−Bu) gc=ˆ R(g−u) PAMG(Bc,uc,rc,gc) x=x+˜ Pxc PS(B,u,b,g) end if Here ˜ Rand ˆ Rdenote the restriction operators for the solution of the LCP and its constraint respectively. The prolongation for the LCP is denoted by ˜ P. Finally PS is a smoother for the LCP. For more details on these operators, see [26]. Finally, we summarize our algorithm to numerically price American options under the Bates model. Note that PSOR or PAMG form inner iterations to the outer LCP-iteration. In the next section we will see that for PAMG, both the outer and inner iteration-count is kept very low for each time-step. Algorithm Discretize (3) with (1) giving (8) with (9) and (10) for k=1,...,l(Time-stepping) for j=1,2,...until convergence (LCP-iteration) Solve (11) using PSOR or PAMG end for end for 5. Numerical Experiment In our numerical example, we price American call options. The parameters for the Bates model are the same as in [22] and they are defined below. 1142 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 Parameter Notation Value Risk free interest rate r0.03 Dividend yield q0.05 Strike price K100 Correlation between the price and variance processes ρ-0.5 Mean level of the variance θ0.04 Rate of reversion to the mean level κ2.0 Volatility of the variance σ0.25 Jump rate λ0.2 Mean jump γ-0.5 Variance of jump δ0.4 The computational domain is (x,y,τ)∈[0,400] ×[0,1] ×[0,0.5]. For the PSOR method, the relaxation parameter ω=1.5 is used. In Table 1 we report the numerical results. The table has the following columns: •Grid (m,n,l) defines the number of grid-points in x,y, and τto be m,n, and l, respectively. •LCP iter. gives the average number of LCP iterations on each time step. •PSOR/PAMG iter. gives the average number of inner PSOR or PAMG iterations for solving one LCP. •Error gives the root mean square relative error given by error =⎡⎢⎢⎢⎢⎢⎢⎣ 1 5 5  i=1u(xi,θ,T)−U(xi,θ,T) U(xi,θ,T)2⎤⎥⎥⎥⎥⎥⎥⎦ 1/2 , where x=(80,90,100,110,120)T. The reference prices Ugiven in [22] at (xi,θ,T), i=1,2,...,5 are 0.328526, 2.109397, 6.711622, 13.749337, 22.143307. They were computed using a componentwise splitting method on the grid (4096,2048,514). •Ratio is the ratio of the consecutive root mean square relative errors. •CPU gives the CPU time in seconds on a 2.0 GHz Intel Core i7 PC using one thread. For the PAMG method, the CPU time includes the AMG initialization time. For the iterations, we use the stopping criterion   ¯ rj  2≤10−6b2, where ¯ rjis the reduced residual for the LCP iterations and the pure PSOR iterations respectively. It is defined by ¯ rj i=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Buj i−biif uj i>gi 0ifuj i=gi. For the inner PSOR/PAMG iterations it is defined similarly with Tand the associated vectors instead of B,uj, and b. The multiplication by the matrix Jis the most expensive operation in the iteration. In order to perform it efficiently with the LCP iterations, we collected all n+1 multiplications corresponding to all x-grid lines together and then performed the resulting matrix-matrix multiplication using the optimized GotoBLAS library [41]. In Table 1, roughly second-order accuracy is observed with the proposed discretization as the ratio is about four on average. On finer grids, only two LCP iterations are required to satisfy the stopping criterion. With the coarsest grid (64,32,10), the LCP iterations with the PSOR and PAMG methods require the same amount of time while the pure PSOR is twice slower. On finer grids, the speed differences become large and the number of PSOR iterations roughly doubles with each refinement. On the finest grid (1024,512,130), the LCP iteration with the PAMG method is about 12 times faster than the LCP iteration with the PSOR method, and it is 150 times faster than the pure PSOR iteration. 1143 Santtu Salmi et al. / Procedia Computer Science 18 ( 2013 ) 1136 – 1144 Table 1. The numerical results with five different space-time discretizations. method Grid (m,n,l) LCP iter. PSOR/PAMG iter. error ratio CPU PSOR (64,32,10) 47.6 0.11765 0.09 (128,64,18) 50.2 0.04068 2.89 1.24 (256,128,34) 109.6 0.00740 5.49 42.00 (512,256,66) 216.8 0.00227 3.26 1608.39 (1024,512,130) 396.6 0.00038 5.93 251758.87 LCP iter. (64,32,10) 2.4 35.5 0.11759 0.04 with PSOR (128,64,18) 2.0 44.4 0.04064 2.89 0.28 (256,128,34) 2.0 101.0 0.00738 5.50 5.19 (512,256,66) 2.0 212.4 0.00226 3.27 107.91 (1024,512,130) 2.0 417.0 0.00038 5.96 1673.90 LCP iter. (64,32,10) 2.5 2.9 0.11760 0.04 with PAMG (128,64,18) 2.1 3.5 0.04064 2.89 0.26 (256,128,34) 2.0 3.5 0.00738 5.50 2.10 (512,256,66) 2.0 3.3 0.00226 3.26 17.85 (1024,512,130) 2.0 2.6 0.00038 5.93 133.11 6. Conclusions In this paper we considered a linear complementarity problem (LCP) with a partial integro-differential operator for pricing American options under the Bates model. For the partial derivatives and integral we employed finite differences and simple quadrature respectively. In the numerical experiments, the discretizations are roughly second-order accurate in both space and time. We proposed a rapidly converging iteration for solving LCPs at each time step. In each such iteration, an LCP with a sparse matrix needs to be solved. We demonstrated that these problems can be efficiently and easily solved with a projected algebraic multigrid method. With finer discretizations this approach leads to an order or several orders of magnitude faster method than using the projected SOR method. References [1] F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Political Economy 81 (1973) 637–654. [2] R. C. Merton, Option pricing when underlying stock returns are discontinuous, J. Financial Econ. 3 (1976) 125–144. [3] S. G. Kou, A jump-diffusion model for option pricing, Management Sci. 48 (8) (2002) 1086–1101. [4] S. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financial Stud. 6 (1993) 327–343. [5] D. S. Bates, Jumps and stochastic volatility: Exchange rate processes implicit Deutsche mark options, Review Financial Stud. 9 (1) (1996) 69–107. [6] D. Duffie, J. Pan, K. Singleton, Transform analysis and asset pricing for affine jump-diffusions, Econometrica 68 (6) (2000) 1343–1376. [7] M. J. Brennan, E. S. Schwartz, The valuation of American put options, J. Finance 32 (1977) 449–462. [8] S. Ikonen, J. Toivanen, Pricing American options using LU decomposition, Appl. Math. Sci. 1 (49-52) (2007) 2529–2551. [9] N. Clarke, K. Parrott, Multigrid for American option pricing with stochastic volatility, Appl. Math. Finance 6 (1999) 177–195. [10] S. Ikonen, J. Toivanen, Componentwise splitting methods for pricing American options under stochastic volatility, Int. J. Theor. Appl. Finance 10 (2) (2007) 331–361. [11] S. Ikonen, J. Toivanen, Efficient numerical methods for pricing American options under stochastic volatility, Numer. Methods Partial Differential Equations 24 (1) (2008) 104–126. [12] K. Ito, J. Toivanen, Lagrange multiplier approach with optimized finite difference stencils for pricing American options under stochastic volatility, SIAM J. Sci. Comput. 31 (4) (2009) 2646–2664. [13] C. W. Oosterlee, On multigrid for linear complementarity problems with application to American-style options, Electron. Trans. Numer. Anal. 15 (2003) 165–185. [14] R. Zvan, P. A. Forsyth, K. R. Vetzal, Penalty methods for American options with stochastic volatility, J. Comput. Appl. Math. 91 (2) (1998) 199–218. [15] A. Almendral, C. W. Oosterlee, Numerical valuation of options with jumps in the underlying, Appl. Numer. Math. 53 (1) (2005) 1–18. [16] L. Andersen, J. Andreasen, Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing, Rev. Deriv. Res. 4 (3) (2000) 231–262. [17] R. Cont, E. Voltchkova, A finite difference scheme for option pricing in jump diffusion and exponential L´ evy models, SIAM Numer. Anal. 43 (4) (2005) 1596–1626. [18] Y. d’Halluin, P. A. Forsyth, K. R. Vetzal, Robust numerical methods for contingent claims under jump diffusion processes, IMA J. Numer. Anal. 25 (1) (2005) 87–112.