Double Discretization Difference Schemes for Partial Integrodifferential Option Pricing Jump Diffusion Models

Hindawi Publishing Corporation Abstract and Applied Analysis Volume 2012, Article ID 120358, 20 pages doi:10.1155/2012/120358 Research Article Double Discretization Difference Schemes for Partial Integrodifferential Option Pricing Jump Diffusion Models M.-C. Casab´ an, R. Company, L. J ´ odar, and J.-V. Romero Instituto de Matem´ atica Multidisciplinar, Universitat Polit` ecnica de Val` encia, Camino de Vera s/n, 46022 Valencia, Spain Correspondence should be addressed to L. J´ odar, [email protected].es Received 10 September 2012; Revised 7 November 2012; Accepted 7 November 2012 Academic Editor: Carlos Vazquez Copyright q2012 M.-C. Casab´ an et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. A new discretization strategy is introduced for the numerical solution of partial integrodifferential equations appearing in option pricing jump diffusion models. In order to consider the unknown behaviour of the solution in the unbounded part of the spatial domain, a double discretization is proposed. Stability, consistency, and positivity of the resulting explicit scheme are analyzed. Advantages of the method are illustrated with several examples. 1. Introduction Since empirical studies revealed that the normality of the log returns, as assumed by Black and Scholes, could not capture features like heavy tails and asymmetries observed in market-data log-returns densities 1, a number of models try to explain these empirical observations: stochastic volatility 2,3, deterministic local volatility 4,5, jump diffusion 6,7, and infinite activity L´ evy models 8–11. The two last types of models, discussed in 12and 13, chapters 14, 15allow to calibrate the model to market price of options and reproduce a wide variety of implied volatility skews/smiles. These models are characterized by partial integrodifferential equations PIDEsthat involve a second-order differential operator and a nonlocal integral term that requires specific treatment and presents additional difficulties. In order to solve the PIDE problem numerically, Andersen and Andreasen 14use an unconditionally stable ADI finite difference method and accelerate it using fast Fourier transform FFT.In15–17wavelet methods are applied to infinite activity jump-diffusion 2 Abstract and Applied Analysis models. Interesting analytic-numerical treatments for L´ evy models have been introduced in 18–20. The so-called COS method for pricing European options is presented in 18. This is based on the knowledge of the characteristic function of the jump operator and the close relation of the characteristic function with the series coefficients of the Fourier-cosine expansion of the density function. In 19, an expansion of the characteristic function of local volatility models with L´ evy jumps is developed. The authors in 20derive an analytical formula for the price of European options for any model including local volatility and Poisson jump process by using Malliavin calculus techniques. Various authors apart of 14used finite difference schemes for PIDEs in 21–27. Discretization of the integral term leads to full matrices due to its nonlocal character. Dealing with the integral term several challenges arise, for instance, how to approximate the integral term and how to localizate a bounded computational domain, also the selection of the boundary conditions of the numerical domain and the problem of the double discretization of the differential and integral part of the PIDE. Tavella and Randall in 26used an implicit time discretization and propose a stationary fairly rapid convergent iterative method to solve the full matrix problem quoted above but without a careful numerical analysis. A generalization of this iterative method to price American options is proposed in 25. In the outstanding paper 22the authors propose an explicit-implicit finite difference scheme for solving parabolic PIDEs with possibly singular kernels when the random evolution of the underlying asset is driven by a time-inhomogeneous jump-diffusion process. The authors study stability and convergence of the proposed scheme as well as rates of convergence. However, they use backward or forward difference quotients of only first order depending on the sign of the coefficient of the convection term in order to avoid oscillations. An improvable issue of 22is that in order to approximate the truncated integral term, they assume a particular behavior of the solution outside of the bounded numerical domain. An efficient solution of PIDEs for the jump-diffusion Merton model is proposed in 24with a very efficient treatment of the resulting dense linear system by using a circulant preconditioned conjugate gradient method. However, in 24, they only consider the particular case where the jump sizes have zero mean, μJ0. They also assume a particular behavior of the solution outside the bounded numerical domain. Almendral and Oosterlee 28present an implicit discretization of the PIDE jumpdiffusion model on an uniform grid using finite differences, where a splitting technique combined with FFT is used to accelerate the dense matrix-vector product. The authors also assume a particular behaviour of the solution outside of the bounded numerical domain, in a similar way to 24. In 21afinitedifference method for PIDE associated with the CGMY infinite activity L´ evy model is treated. The equations are discretized in space by the collocation method and in time by an explicit backward differentiation formula. The integral part is transformed into a Volterra equation. After integration by parts and taking advantage of the vanishing derivative behaviour of the payofffunction for large asset values, the authors are able to truncate properly the integral for the case of put and butterfly options. In 27the price of European and American options under PIDE Kou’s jump-diffusion model is solved using finite differences on nonuniform grids, and time stepping is performed using the implicit Rannacher scheme. The evaluation of the integral term is efficient from the computational cost point of view, assuming that the behaviour of the solution for large values of the underlying asset follows the asymptotic behaviour. Abstract and Applied Analysis 3 For the sake of clarity in the presentation we recall that in a jump-diffusion model, the modified stochastic differential equation SDEfor the underlying asset is dS Sμdt σdz η−1dq, 1.1 where Sis the underlying stock price, μis the drift rate, σis the volatility, dz is the increment of Gauss-Wiener process, and dq is the Poisson process. The random variable representing the jump amplitude is denoted by η, and the expected relative jump size is denoted by K Eη−1. The jump intensity of the Poisson process is denoted by λ. Based on the SDE 1.1 the resulting PIDE for a contingent claim VS, tis given by 7,14,29: ∂V ∂t 1 2σ2S2∂2V ∂S2r−λKS∂V ∂S −rλV λ∞ 0 VSη, tgηdη 0,0<S<∞,0≤t<T, 1.2 VS, TfS,0<S<∞,1.3 where ris the risk-free interest rate, the probability density of the jump amplitude is given by gη,andfSis the payofffunction. Merton’s jump-diffusion model assumes that jump sizes are log-normally distributed with mean μJand standard deviation σJ,thatis, gη exp−1/2lnη−μJ/σJ2 σJη√2π.1.4 In this paper we consider Merton’s jump-diffusion model for a vanilla call option with payofffunction fSmaxS−E, 0.1.5 The aim of the paper is the construction and numerical analysis of an explicit finite difference numerical scheme of the PIDE 1.2-1.3,withadifferent treatment of the integral part from previously quoted authors. Instead of assuming long-term information about the solution we perform a full discretization of the integral part, involving the unknown function values in the numerical scheme, discriminating the finite truncation domain and the infinite remaining one. As a consequence, this strategy involves a double discretization with respect to the spatial variable. With respect to the time variable an explicit forward approximation is used. An account of the advantages of this explicit approach has been explained and applied in 30. This paper is organized as follows. Section 2 deals with a transformation of variables in order to eliminate both the advection and the reaction terms of the PIDE 1.2. Then the integral part is split in two parts: a finite integral J1and an infinite one J2, and the last is again transformed into a finite integral. The separation point Aof the two split integrals is becoming a parameter that could be chosen according to the criteria used by 16,22,31. A suitable 4 Abstract and Applied Analysis choice of parameter Ais the one used by other authors when they truncate the numerical domain. For instance in 27, one takes A4E;inSection 6 we take A3E.Section 3 deals with the construction of the numerical scheme and the selection of the numerical domain that always involves the difficulty of the consideration of the boundary conditions. For the case of a PIDE this issue is even more relevant because the values throughout all the unbounded integral domain are unknown. The spatial numerical domain is divided in two parts by the parameter A:0,Aand A, ∞.Inthe0,Adomain, the stepsize discretization is hΔx, consisting in Nequidistributed mesh points Xi,1≤i≤N. The domain A, ∞is transformed into the 0,1by transformation zA/X. In the transformed domain 0,1we consider a stepsize discretization δΔzand Mmesh points with Mδ 1. When the interval 0,1is reversed to the domain A, ∞, the reversed mesh points Xi,N ≤i≤NM−1 become nonuniformly distributed. Hence, the numerical scheme for problem 2.9-2.10is forward in time with time-step discretization kΔτ. The approximation of ∂2U/∂X2is centered in 0,Aof the unique parameter h,andinA, A/δthe approximation of ∂2U/∂X2involves a nonuniform stepsize discretization hjdepending on δand the value of j. The numerical approximation of the integrals is evaluated using trapezoidal quadrature rules with stepsizes hand δ, respectively. The boundary conditions at the boundary of our numerical domain are as follows. At X0 we assume that the solution is zero according to the vanilla call option problem. At our largest finite value considered XA/δ we assume a linear behavior of the solution. This hypothesis has been previously used in 32. In Section 4 sufficient conditions for stability and positivity of the numerical solution are given in terms of the three parameters hΔX, k Δτ, and δΔzas well as of the parameter A. Consistency of the scheme is treated in Section 5.Section 6 includes illustrative examples showing the possible advantages of our new discretization approach. Finally conclusions are shown in Section 7. If vv1,v 2,...,v nTis a vector in Rn, we denote its infinite norm v∞ max{|vj|;i≤j≤n}.Vectorvis said to be nonnegative if vj≥0 for all 1 ≤j≤n. Then we denote v≥0. For a matrix Aaij m×nin Rm×n, we denote by A∞max1≤i≤m{n j1|aij|}. Matrix Ais said to be nonnegative if aij ≥0 for all 1 ≤i≤m, 1≤j≤n, and we denote A≥0. 2. Transformation of the Integrodifferential Problem For the sake of convenience we introduce a transformation of variables to remove both the advection and the reaction terms of the PIDE problem 1.2-1.3. Let us consider the transformation Xexpr−λKT−tS, τ T−t, UX, τexprλT−tVS, t, 2.1 and note that problem 1.2-1.3is transformed into the problem ∂U ∂τ 1 2σ2X2∂2U ∂X2λ∞ 0 UXη,τgηdη, 0<X<∞,0<τ≤T, 2.2 UX, 0fX,0<X<∞.2.3 Abstract and Applied Analysis 5 In order to approximate the integral appearing in 2.2and further discretization it is convenient the change of the variable φXη, ∞ 0 UXη,τgηdη 1 X∞ 0 Uφ,τgφ Xdφ. 2.4 Let us denote JJX, τ∞ 0 Uφ,τgφ Xdφ. 2.5 Taking A>0, let us decompose JJ1J2,J 1A 0 Uφ,τgφ Xdφ, J2∞ A Uφ,τgφ Xdφ, J1J1X, τ, A,J 2J2X, τ, A. 2.6 Following 33, page 201let us consider the substitution zA φ,2.7 into J2, obtaining the expression J2A1 0 UA z,τgA Xz1 z2dz. 2.8 Taking into account 2.4–2.8, the problem 2.2-2.3can be written in the form ∂U ∂τ σ2X2 2 ∂2U ∂X2λ XJ1J2,0<X<∞,0<τ≤T, 2.9 UX, 0fX,0<X<∞.2.10 3. Numerical Scheme Construction In this section a difference scheme for problem 2.9-2.10is constructed. With respect to the time variable, given τwith 0 <τ≤T,letkbe the time-step discretization kΔττ/L,and τnnk,0≤n≤L,withLinteger. With respect to the spatial variable X, given an arbitrary positive fixed A>0, we construct a uniform grid in 0,A, with the spatial step discretization hΔXA/N,withXjjh, 0≤j≤N, being Ninteger. Note that the integral J2X, τ, Agiven by 2.6requires the evaluation of the unknown Uφ,τat points φ∈A, ∞. As this integral has been transformed into 2.8over the interval 6 Abstract and Applied Analysis 0,1for the variable z,see2.7, we consider a uniform mesh of 0,1into Mpoints, of the form zjjδ, 1≤j≤M, where Mis integer, Mδ 1, with M≥3. Taking into account 2.7, for the original variable Xin A, ∞, one has XjA zNM−j ,N≤j≤NM−1,3.1 and since zjjδ, XjA 1−j−Nδ,N≤j≤NM−1.3.2 Thus the spatial domain 0,A∪A, A/δ0,A/δis split into NMpoints Xj⎧ ⎪ ⎨ ⎪ ⎩ jh, 0≤j≤N, A 1−j−Nδ,N≤j≤NM−1,3.3 from those only the first N1 are equidistributed. Let us denote un j≈UXj,τn,0≤j≤NM−1,0≤n≤L, ∂U ∂τ Xj,τn≈un1 j−un j k,1≤j≤NM−1,0≤n≤L. 3.4 For the approximation of ∂2U/∂X2we consider two types of finite differences: ∂2U ∂X2Xj,τn≈un j1−2un jun j−1 h2Δ n j,1≤j≤N−1,3.5 for the internal points of 0,A, and denoting hjXj1−Xj>0, ∂2U ∂X2Xj,τn≈2un j1 hjhjhj−1un j−1 hj−1hjhj−1−un j hjhj−1Δ n j,N≤j≤NM−2, 3.6 for the points Xjlying in A, A/2δ. Note that the discrete operator Δn jhas different expressions depending on the ubication of j,see3.5and 3.6. From the previous approximations, for the internal points we have un1 iun ik 2σ2X2 iΔn ikλ XiJn 1,i Jn 2,i,1≤i≤NM−2,3.7 Abstract and Applied Analysis 7 where Jn 1,i and Jn 2,i are approximations of the composite trapezoidal type of integrals appearing in 2.6,2.8: Jn 1,i ≈A 0 Uφ,τngφ Xidφ, 3.8 Jn 2,i ≈∞ A Uφ,τngφ Xidφ A1 0 UA z,τngA Xiz1 z2dz. 3.9 Let us denote gi,j gXj/Xi. The approximation Jn 1,i takes the form Jn 1,i h⎛ ⎝ N−1  j1 un jgi,j 1 2un Ngi,N⎞ ⎠,1≤i≤NM−2,3.10 where the first term for j0 does not appear due to the null value of the limit of function gηgiven by 1.4as ηtends to zero. On the other hand, considering the assumption that Uφ,τnhas asymptotic linear behaviour as φ→∞z→0and using 1.4it follows that the integrand of 3.9verifies UA z,τngA Xiz1 z2−→ 0,as z−→ 0.3.11 Consequently, the last term of Jn 2,i related to jNMvanishes, and one gets Jn 2,i δ A⎛ ⎝ 1 2un Ngi,NX2 N NM−1  jN1 un jgi,jX2 j⎞ ⎠,1≤i≤NM−2.3.12 The numerical scheme 3.7–3.12needs to incorporate the transformed initial condition u0 ifXimaxXi−E, 0,1≤i≤NM−1,3.13 and the boundary conditions for i0: un 00,0≤n≤L, 3.14 and assuming linear behaviour of the solution for large values of the spatial variable at any time, we have Δn NM−10 and null integral term approximation Jn 1,NM−1Jn 2,NM−10. Hence, considering 3.7for iNM−1, one gets un1 NM−1un NM−1u0 NM−1,0≤n≤L−1.3.15 8 Abstract and Applied Analysis For the sake of convenience in the study of stability we introduce a vector formulation of the scheme 3.7–3.15. Let us consider the vector in RNM−1as Unun 1un 2... u n NM−1t,3.16 and let P∈RNM−1×NM−1be the tridiagonal matrix related to the differential part and defined by P ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ β1γ100 ··· 0 α2β2γ20··· 0 0α3β3γ3··· 0 ......... ......... 0··· αNM−2βNM−2γNM−2 0··· 0αNM−1βNM−1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,3.17 where γ1k 2σ2,α NM−10,β NM−11, γiαik 2h2σ2X2 i,2≤i≤N−1, βi1−k h2σ2X2 i,1≤i≤N−1, αNkσ2A2 hhAδ/1−δ,β N1−kσ2A1−δ δh ,γ Nkσ2A1−δ δhAδ/1−δ, αikσ2X2 i hi−1hihi−1,β i1−kσ2X2 i hihi−1 ,γ ikσ2X2 i hihihi−1,N1≤i≤NM−2. 3.18 Let Bbijbe the matrix in ∈RNM−1×NM−1related to the integral part whose entries bij for each fixed iin 1 ≤i≤NM−2 are defined by bij  ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ kλ Xi hgi,j,1≤j≤N−1, kλ 2Xi hδAgi,N,jN, kλ Xi δ Agi,jX2 j,N1≤j≤NM−1, bNM−1,j 0,1≤j≤NM−1. 3.19 Abstract and Applied Analysis 9 From the previous notation the scheme 3.7–3.15can be written in the form Un1PBUnPBn1U0,0≤n≤L−1, U0fX1fX2··· fXNM−1t. 3.20 4. Positive and Stability of the Numerical Solution Dealing with prices of contracts modeled by PIDE, the solution must be nonnegative. In this section we show that numerical solution provided by scheme 3.7–3.15is conditionally positive and stable. We begin with the following result. Lemma 4.1. With previous notation, assume that stepsizes kΔτ,hΔXin 0,Aand 0<δ≤ 1/3, and δΔzin 0,1satisfy C1k/h2≤1/σ2A2, C2k≤min{δ2/σ21−2δ,δh/σ2A1−δ}. Then matrix Pgiven by 3.17is nonnegative. Proof. From 3.18,for1≤i≤N−1 one has γi>0andαi>0for2≤i≤N−1. On the other hand, for 1 ≤i≤N−1, we have that βi≥0,iffk h2≤1 σ2X2 i .4.1 Thus, under condition C1, condition 4.1holds true. With respect to the nonuniform grid, note that for N≤i≤NM−2, from 3.18one gets that αi>0, γi>0, and αNM−10. From 3.18we also have that βN≥0,iffk≤δh σ2A1−δ.4.2 In order to guarantee the nonnegativeness of the remaining entries of matrix P,letus introduce the function Hihihi−1 X2 i δ2 1−i−N1δ1−i−N−1δ>0,4.3 for N1≤i≤NM−2. With this notation, we have that βiappearing in 3.18satisfies βi≥0,iffk≤Hi σ2.4.4 16 Abstract and Applied Analysis 5 10 15 20 25 30 35 40 45 0 10 20 30 40 50 60 70 80 S Option price k=0.01 k=0.015625 Figure 1: Satisfying and breaking stability conditions. Thus Tn iUOh2Oδ2Ok,5.14 which proves the consistency of the scheme with the PIDE. 6. Numerical Results In the following examples the code was run on Matlab. The first example illustrates that stability conditions of Theorem 4.5 cannot be removed. Example 6.1. Consider the vanilla call option problem 1.2–1.5under Merton jump diffusion model with parameters T1, r0.05, E10, σ0.1, μJ0.5, K0.7, and λ0.1. Taking A3E,δ0.1, and h0.3, Figure 1 shows that when the stability conditions 4.26are satisfied results are good k0.01, while if the stability conditions are broken k0.015625results are unreliables. The next example shows the robustness of our numerical scheme under changes of the jump intensity λof the model. Example 6.2. Taking the same parameters of Example 6.1 apart from λand the stepsize discretizations h0.3, δ0.1, and k0.01, Figure 2 shows the variation of the solution with parameter λ, where λ0 corresponds to the Black-Scholes case. In the next example, the error is the difference between the numerical solution VS, 0 computed by 3.7–3.15and 2.1and the exact solution given by Merton’s formula 7. Example 6.3 shows that the error of the numerical solution with fixed δdecreases with the uniform stepsize hΔXabout the strike E10, while the error close to the truncation Abstract and Applied Analysis 17 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 S Option price λ=0 λ=0.1 λ=0.2 Figure 2: Variation of the jump intensity λ. 0 5 10 15 20 25 −4 −3 −2 −1 0 1 2 3 4 S Absolute errors ×10−3 h=0.6 h=0.3 h=0.1 h=0.075 Figure 3: Absolute errors with several values of hand a fixed δ. separation point A30 remains stationary when hdecreases. This fact agrees with facts illustrated in 28, pages 15-16. Example 6.3. Consider the vanilla call option problem 1.2–1.5under Merton jump diffusion model with parameters T1, r0.05, E10, σ0.1, μJ0, K0.00501, and λ0.1. For A30, k0.001, and δ0.0625, the Figure 3 shows the variation of the absolute error of the solution under changes of the stepsize h. 18 Abstract and Applied Analysis 0 10 20 30 40 50 60 −0.015 −0.01 −0.005 0 0.005 0.01 S Absolute errors δ=0.0625 δ=0.03125 δ=0.01 Figure 4: Absolute errors with several values of δand a fixed h. The next Example 6.4 shows that the errors in the right boundary of the numerical domain when one uses finite difference schemes, quoted by 28, can be reduced with our double spatial discretization by decreasing the stepsize δ. Example 6.4. Taking the problem of Example 6.3 with fixed h0.3, Figure 4 shows the error reduction of the numerical solution about the right boundary of the numerical domain when parameter δdecreases, while the error about the strike remains stationary. 7. Conclusions This work introduces a new discretization strategy for solving partial integrodifferential equations which involves the discretization of the unknown in the unbounded part of the integral. This fact increases the accuracy of the numerical solution in the boundary of the numerical domain as it is shown in Example 6.4. Acknowledgment This paper was supported by the Spanish M.E.Y.C. Grant DPI2010-20891-C02-01. References 1J. Y. Campbell, A. W. Lo, and A. C. MacKinlay, The Econometrics of Financial Markets,Princeton University Press, Princeton, NJ, USA, 1997. 2S. Heston, “A closed-form solution for options with stochastic volatility with applications to bond and currency options,” Review of Financial Studies, vol. 6, no. 2, pp. 327–343, 1993. 3J. Hull and A. White, “The pricing of options with stochastic volatilities,” The Journal of Finance, vol. 42, no. 2, pp. 281–300, 1987. 4J. C. Cox and S. A. Ross, “The valuation of options for alternative stochastic processes,” Journal of Financial Economics, vol. 3, no. 1-2, pp. 145–166, 1976. Abstract and Applied Analysis 19 5B. Dupire, “Pricing with a smile,” Risk Magazine, vol. 1, pp. 18–20, 1994. 6S. G. Kou, “A jump diusion model for option pricing,” Management Science, vol. 48, no. 8, pp. 1086– 1101, 2002. 7R. C. Merton, “Option pricing when the underlying stocks are discontinuous,” Journal of Financial Economics, vol. 3, pp. 1125–2144, 1976. 8O. E. Barndorff-Nielsen, “Processes of normal inverse Gaussian type,” Finance and Stochastics, vol. 2, no. 1, pp. 41–68, 1998. 9E. Eberlein, “Application of generalized hyperbolic L´ evy motions to finance,” in L´ evy Processes—Theory and Applications, O. Barndor-Nielsen, T. Mikosch, and S. Resnick, Eds., pp. 319–336, Birkh¨ auser, Boston, Mass, USA, 2001. 10I. Koponen, “Analytic approach to the problem of convergence of truncated L´ evy ights towards the Gaussian stochastic process,” Physical Review E, vol. 52, no. 1, pp. 1197–1199, 1995. 11D. Madan and F. Milne, “Option pricing with variance gamma martingale components,” Mathematical Finance, vol. 1, no. 4, pp. 39–55, 1991. 12R. Cont and P. Tankov, Financial Modelling with Jump Processes, Chapman & Hall/CRC Financial Mathematics Series, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2004. 13A. Pascucci, PDE and Martingale Methods in Option Pricing, vol. 2 of Bocconi & Springer Series,Springer, Milan, Italy, 2011. 14L. Andersen and J. Andreasen, “Jump-diffusion processes: volatility smile fitting and numerical methods for option pricing,” Review of Derivatives Research, vol. 4, no. 3, pp. 231–262, 2000. 15A.-M. Matache, P.-A. Nitsche, and C. Schwab, “Wavelet Galerkin pricing of American options on L´ evy driven assets,” Quantitative Finance, vol. 5, no. 4, pp. 403–424, 2005. 16A.-M. Matache, T. von Petersdorff, and C. Schwab, “Fast deterministic pricing of options on L´ evy driven assets,” Mathematical Modelling and Numerical Analysis, vol. 38, no. 1, pp. 37–71, 2004. 17A.-M. Matache, C. Schwab, and T. P. Wihler, “Fast numerical solution of parabolic integrodifferential equations with applications in finance,” SIAM Journal on Scientific Computing, vol. 27, no. 2, pp. 369– 393, 2005. 18F. Fang and C. W. Oosterlee, “A novel pricing method for European options based on Fourier-cosine series expansions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 826–848, 2008/09. 19S. Pagliarani, A. Pascucci, and C. Riga, “Adjoint expansions in local L´ evy models,” SSRN eLibrary, 2011. 20E. Benhamou, E. Gobet, and M. Miri, “Smart expansion and fast calibration for jump diffusions,” Finance and Stochastics, vol. 13, no. 4, pp. 563–589, 2009. 21A. Almendral and C. W. Oosterlee, “Accurate evaluation of European and American options under the CGMY process,” SIAM Journal on Scientific Computing, vol. 29, no. 1, pp. 93–117, 2007. 22R. Cont and E. Voltchkova, “A finite difference scheme for option pricing in jump diffusion and exponential L´ evy models,” SIAM Journal on Numerical Analysis, vol. 43, no. 4, pp. 1596–1626, 2005. 23Y. d’Halluin, P. A. Forsyth, and G. Labahn, “A penalty method for American options with jump diffusion processes,” Numerische Mathematik, vol. 97, no. 2, pp. 321–352, 2004. 24E. W. Sachs and A. K. Strauss, “Efficient solution of a partial integro-differential equation in finance,” Applied Numerical Mathematics, vol. 58, no. 11, pp. 1687–1703, 2008. 25S. Salmi and J. Toivanen, “An iterative method for pricing American options under jump-diffusion models,” Applied Numerical Mathematics, vol. 61, no. 7, pp. 821–831, 2011. 26D. Tavella and C. Randall, Pricing Financial Instruments, Wiley, New York, NY, USA, 2000. 27J. Toivanen, “Numerical valuation of European and American options under Kou’s jump-diffusion model,” SIAM Journal on Scientific Computing, vol. 30, no. 4, pp. 1949–1970, 2008. 28A. Almendral and C. W. Oosterlee, “Numerical valuation of options with jumps in the underlying,” Applied Numerical Mathematics, vol. 53, no. 1, pp. 1–18, 2005. 29M. Briani, C. La Chioma, and R. Natalini, “Convergence of numerical schemes for viscosity solutions to integro-differential degenerate parabolic problems arising in financial theory,” Numerische Mathematik, vol. 98, no. 4, pp. 607–646, 2004. 30R. Company, L. J´ odar, E. Ponsoda, and C. Ballester, “Numerical analysis and simulation of option pricing problems modeling illiquid markets,” Computers & Mathematics with Applications, vol. 59, no. 8, pp. 2964–2975, 2010. 31R. Kangro and R. Nicolaides, “Far field boundary conditions for Black-Scholes equations,” SIAM Journal on Numerical Analysis, vol. 38, no. 4, pp. 1357–1368, 2000. 20 Abstract and Applied Analysis 32H. Windcli, P. A. Forsyth, and K. R. Vetzal, “Analysis of the stability of the stability of the linear boundary condition for the Black-Scholes equation,” Journal of Computational Finance,vol.8,no.1,pp. 65–92, 2004. 33P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Computer Science and Applied Mathematics, Academic Press, New York, NY, USA, 2nd edition, 1984. 34P. Linz, Analytical and Numerical Methods for Volterra Equations, vol. 7, Society for Industrial and Applied Mathematics SIAM, Philadelphia, Pa, USA, 1985. 35G. D. Smith, Numerical Solution of Partial Differential Equations, Clarendon Press, Oxford, UK, 3rd edition, 1985.