- Email: [email protected]
A radial basis function based implicit–explicit method for option pricing under jump-diffusion models
Accepted Manuscript A radial basis function based implicit-explicit method for option pricing under jump-diffusion models
Mohan K. Kadalbajoo, Alpesh Kumar, Lok Pati Tripathi
PII: DOI: Reference:
S0168-9274(16)30156-8 http://dx.doi.org/10.1016/j.apnum.2016.08.006 APNUM 3077
To appear in:
Applied Numerical Mathematics
Received date: Revised date: Accepted date:
27 March 2015 17 June 2016 19 August 2016
Please cite this article in press as: M.K. Kadalbajoo et al., A radial basis function based implicit-explicit method for option pricing under jump-diffusion models, Appl. Numer. Math. (2016), http://dx.doi.org/10.1016/j.apnum.2016.08.006
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A radial basis function based implicit-explicit method for option pricing under jump-diffusion models Mohan K. Kadalbajoo, Alpesh Kumar∗, Lok Pati Tripathi Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur 208016, India. Abstract In this article, we present a radial basis function based implicit explicit numerical method to solve the partial integro-differential equation which describe the nature of the option price under jump diffusion model. The governing equation is time semi discrtized by using the implicit-explicit backward difference method of order two (IMEXBDF2) followed by radial basis function based finite difference (RBF-FD) method. The numerical scheme derived for European option is extended for American option by using operator splitting method. Numerical results for put and call option under Merton and Kou models are given to illustrate the efficiency and accuracy of the present method. The stability of time semi discretized scheme is also proved.
Keywords: Radial basis function; Finite difference, Option pricing; Jump-diffusion models; Partial integro-differential equation.
1
Introduction
There is evidence to suggest that the Black Scholes model for stock price behavior does not always model real stock price behavior. Jump can appear at a random time and these jumps can not be captured by the log normal distribution characteristic of the stock price in the Black Scholes model. To overcome the above shortcoming, several models have been proposed in the literature. Among these, the jump diffusion models introduced by Merton[23] and Kou[18] are of the most widely used models. Merton proposed a log-normally distributed process for the jump-amplitudes, whereas Kou suggested log-double-exponentially distributed process. The valuation of option under jump diffusion process requires the solution of a partial integro differential equation containing a non local integral term. There are several numerical methods available in the literature to approximate the above equation. In [1], Almendral and Osterlee presented an implicit second order accurate time discretization with finite difference and finite element spatial discretization on uniform grid. Andersen et al.[2] proposed ∗
corresponding author, Email: [email protected](Alpesh Kumar)
1
an unconditionally stable alternating direction implicit method for its solution. Song Wang et al.[33] developed a fitted finite volume method for jump diffusion process. Their method is based on fitted finite volume method spatial discretization and Crank Nicolson scheme for temporal discretization. More recently, Patidar et al.[24] developed an efficient method for pricing Merton jump diffusion option, combining the spectral domain decomposition method and the Laplace transform method. The scheme proposed by d Halluin et al.[9] required to use an iterative procedure to solve discrete equations. The main difficulty with implicit scheme is due to containing non local integral term in governing equation, which leads to a dense discretization matrix where as fully explicit scheme imposed stability restriction on it. An approach based on implicit-explicit schemes in which integral term is treated explicitly has been proposed by YongHoon Kwon et al.[19] and Briani et al.[3]. More recently Tangman et al.[28] introduced a new scheme called exponential time integration (ETI) scheme to solve the PIDE. In ETI scheme, the time direction of PIDE is directly tackled by a ’one step’ formula, which means temporal discretization is not required. Tangman et al.[28] used the central difference approach with ETI to provide very efficient and second order accurate result. Recently, a new method based on radial basis function (RBF) for approximation of spatial derivative in option pricing equation is under going active research. Application of RBF in one dimension European and American options is given by Hon et al.[13, 14]. Fasshauer et al.[10] solved American option pricing model using penalty method. Golbabai et al.[12], developed an algorithm based on global collocation for jump diffusion process. Bhuruth et al.[25] proposed a radial basis function based differential quadrature rule for spatial descretization with exponential time integration to solve jump diffusion model. In more recent work, Chan et al.[4, 6] used new RBF called cubic spline as basis function to solve PIDE, and show that their scheme is second order accurate. It was recognized that standard approach to solving the radial basis function collocation problem has been ill conditioned due to use of collocation in global sense. Recently many strategies have been developed in the literature to avoid these problems, such as local RBF approach by Lee et al.[21], radial point interpolation method proposed by Liu et al.[22], Shu et al.[27] proposed a local radial basis function-based differential quadrature method and used it to solve two-dimensional incompressible Navier-Stokes equations. Tolstykh[30], Tolstykh and Shirobokov[31], Wright et al.[11, 32] proposed radial basis function finite difference method, the idea is to use radial basis functions with a local collocation as in finite difference mode thereby reducing number of nodes and hence producing a sparse matrix. This technique is further extended by Sanyasiraju et al.[7, 26] for convection diffusion type equations. However these methods have not been extended to solve partial integro differential equation yet. In the present work, we have extended the localization concept proposed by Wright and Fornberg, to solve jump diffusion models. The governing equations are discretized by a three level implicit and explicit time scheme followed by RBF based finite difference method. The paper is organized as follow. In section 2, mathematical models for pricing option with jump diffusion process are given in terms of partial integro-differential equations and provide a brief review of both the Merton and Kou jump diffusion models. Section 3 deals with 2
the construction of three time level implicit explicit scheme to discretize the jump diffusion model. The time semi discrete equation is coupled with radial basis function based finite difference method for spatial discretization. Section 4, provide extension of proposed method for pricing American option by utilizing concept of operator splitting method. In section 5, we give some numerical results for Merton and Kou model and a comparison of the accuracy of our solution with finite difference and finite element method for both American and European options. Finally the paper ends with some conclusive remark in section 6.
2
The mathematical model
In this section, we give brief discussion about the mathematical model for option with jump diffusion process. Consider an asset with the asset price S, then the movement of stock price is modeled by the following stochastic differential equation dS = (ν − κλ)dτ + σdZ + (η − 1)dq S
(2.1)
where, ν is drift rate, τ as the time to maturity, σ represent the constant volatility, dZ is an increment of standard Gauss-Wiener process. The term λ is the intensity of the independent Poisson process dq with 0 with probability 1 − λdτ , dq = 1 with probability λdτ . The expected relative jump size E(η − 1) is denoted by κ, where E[·] is the expectation operator and η − 1 is a impulse function producing jump from S to Sη. Let V (S, τ ) represent the value of a contingent claim that depend on the underlying asset price S with current time τ . Then V (S, τ ) satisfy following backward partial integro differential equation ∞ 1 2 2 ∂ 2V ∂V ∂V + σ S − (r + λ)V + λ + (r − λκ)S V (Sη)g(η)dη = 0, (2.2) ∂τ 2 ∂S 2 ∂S 0 for (S, τ ) ∈ (0, ∞) × (0, T ], where, r is risk free interest rate and g(η) is probability density ∞ function of the jump with amplitude η with properties that ∀η, g(η) ≥ 0 and 0 g(η)dη = 1. The value of V at the expiry date is given by, V (S, T ) = G(S), S ∈ (0, ∞),
(2.3)
where G(S) is the pay-off function for the option contract. Under Merton’s model g(η) is given by the log-normal density 2
(ln η−μJ ) − 1 2σ 2 J g(η) := √ . e 2πσJ η
3
(2.4)
σ2
In this case κ := E(η − 1) = exp(μJ + 2J ) − 1, where μJ and σJ2 are the mean and the variance of jump in return. Under a modified version of Kou’s jump-diffusion model g(η) is the following log-double-exponential density g(η) :=
1 pη1 e−η1 ln(η) H(ln(η)) + qη2 eη2 ln(η) H(− ln(η)) , η
(2.5)
where η1 > 1, η2 > 0, p > 0, q = 1 − p, and H(·) is the Heaviside function. We can show 1 2 that κ : ηpη + ηqη − 1. 1 −1 1 +1 S By using change of variable x = ln( K ), y = ln(η), where K is strike price and letting t = T − τ , computation of the option value requires solving the PIDE ∞ ∂u 1 2 ∂ 2u ∂u σ2 = σ − λκ) − (r + λ)u + λ + (r − u(y, t)f (y − x)dy, (2.6) ∂t 2 ∂x2 2 ∂x −∞ where u(x, t) := V (Kex , T − t) and f (y) = g(ey )ey , under the above transformation the function f (y) can be written as; ⎧ 2 J) ⎨ 1 − (y−μ 2σ 2 √ J e Merton’s model, 2πσJ (2.7) f (y) := ⎩ −η1 y η2 y H(y) + qη2 e H(−y)) Kou’s model. (pη1 e The initial condition and asymptotic behavior of the price of a European put option is described by u(x, 0) = max(K − Kex , 0) Ke−rt − Kex x → −∞, u(x, t) = 0 x → ∞.
(2.8) (2.9)
Other types of initial and boundary conditions can be suitably defines for different types options. In contrast, the American option can be exercised at any time up to the maturity date and can be formulated as the linear complementarity problem(LCP) of the form ⎧ ∂u − Lu ≥ 0, ⎨ ∂t u(x, t) − G(x) ≥ 0, (2.10) ⎩ ∂u ( ∂t − Lu)(u(x, t) − G(x)) = 0, for (x, t) ∈ (0, ∞) × [0, T ) together with the boundary conditions K − Kex x → −∞, u(x, t) = 0 x → ∞, where L is the partial integro differential operator on the right side of (2.6). 4
(2.11)
3
Discretizations
In this section we shall construct an implicit explicit time semi discretization for following PIDE on truncated domain Ω := (xmin , xmax ), ∂u(x, t) = Du(x, t) + λIu(x, t), (x, t) ∈ Ω × [0, T ) ∂t ¯ u(x, 0) = G(x), x ∈ Ω
(3.1) (3.2)
where D is differential operator and I is integral operator of L such that σ2 ∂u 1 2 ∂ 2u + (r − σ − λκ) − (r + λ)u, Du(x, t) = 2 2 ∂x 2 ∂x
(3.3)
u(y, t)f (y − x)dy.
(3.4)
Iu(x, t) =
∞
−∞
In order to approximate the integral operator Iu in (3.4) numerically, we first divide the integral into two parts on Ω and on Ωc := R \ Ω. Now the the integral operator can be split as Iu(x, t) = u(y, t)f (y − x)dy + R(t, x). (3.5) Ω
By using asymptotic behavior of the option the integral over R \ Ω is given by (Ke−rt − Key )f (y − x)dy. R(t, x) :=
(3.6)
Ωc
For the European put option the expression R(t, x) is given as; ⎧
σ2 xmin −x−μJ −σJ2 ⎨ −rt x+μJ + 2J J N Ke N xmin −x−μ − Ke σJ σJ R(t, x) := ⎩K(1 − p)e−rt+η2 (xmin −x) − K(1 − p) η2 e−η2 x+(η2 +1)xmin η2 +1
Merton’s ,
(3.7)
Kou’s.
In the case of American put option above expression can be reformulated as; ⎧
σ2 xmin −x−μJ −σJ2 ⎨ x+μJ + 2J J − Ke Merton’s, N KN xmin −x−μ σJ σJ R(t, x) := ⎩K(1 − p)eη2 (xmin −x) − K(1 − p) η2 e−η2 x+(η2 +1)xmin Kou’s.
(3.8)
η2 +1
where N (·) is the cumulative normal distribution. We use composite trapezoidal rule to approximate the integral term on truncated domain Ω.
3.1
Time discretizations
Let {0 = t0 < t1 < . . . < tN = T ; tn − tn−1 = δt, 1 ≤ n ≤ N } be a partition of the interval [0, T ]. Let us use the notation un := u(x, tn ) then equation (3.1) will be discretized by following implicit explicit scheme,
1 3 n+1 1 n−1 n u = Dun+1 + λI(Eun ), n ≥ 1, − 2u + u (3.9) δt 2 2 5
together with boundary condition un+1 (xmin ) = Ke−rtn+1 − Kexmin un+1 (xmax ) = 0,
(3.10)
where Eun = 2un − un−1 . The above discretization method is called IMEX-BDF2 method with three time levels. In order to use the proposed method we need two initial values on the zeroth and first time level. The value u0 at the zeroth time level is given by initial condition of the model problem, and the value u1 can be obtained by applying the implicit-explicit backward difference method of order one, un+1 − un = Dun+1 + λIun , (3.11) δt together with boundary condition un+1 (xmin ) = Ke−rtn+1 − Kexmin un+1 (xmax ) = 0.
(3.12)
Now we prove the stability of the time discrete scheme (3.9). It can be easily shown that for all u(·, t) ∈ L2 (Ω), t ∈ (0, T ) such that u(x, t) : (x, t) ∈ Ω × [0, T ], u¯(x, t) = 0 : (x, t) ∈ Ωc × [0, T ], the integral operator satisfies the condition
I u¯(·, t) ≤ CI u(·, t) , (3.13) for some constant CI independent of t, where v := ( Ω | v(x) |2 dx)1/2 . Let us suppose that un be the solution of the equation (3.9) and u˜n be the solution of perturbed equation
1 3 n+1 1 n−1 n = D˜ un+1 + λI(E u˜n ) + δ n+1 , n ≥ 1, u˜ − 2˜ u + u˜ (3.14) δt 2 2 together with boundary condition u˜n+1 (xmin ) = Ke−rtn+1 − Kexmin u˜n+1 (xmax ) = 0.
(3.15)
The δ j is the perturbation in the data at any time level tj and what we will examine the effect of these perturbations on the error as we proceed in time. We will show that the error at any time level is bounded by the initial errors at initial time levels(t0 & t1 ; for BDF2) and is a bounded multiple of the introduced perturbations; which establishes stability in time discretization. Now define the error term en := u˜n − un . Since, we are treating boundary conditions exactly in the equations (3.9-3.10) and perturbed equations (3.14-3.15). Hence the error term en will be zero at the boundary points. 6
Theorem 1. (L2 -stability)For sufficiently small δt such that δt <
em 2 ≤ C( e0 2 + e1 2 + max δ j 2 ), 2≤j≤m
1 , (2+4α+2λCI )
∀2≤m≤
T , δt
we have (3.16)
2
(r− σ −λκ)2 −2(r+λ)σ 2
2 where α =| | and C is generic constant depending on the parameter CI , 2σ 2 r, σ, λ, and T . Proof. Relations (3.9-3.10) and (3.14-3.15) implies that the error term en satisfies the following relations (3en+1 − 4en + en−1 ) = Den+1 + λI(Een ) + δ n+1 , (3.17) 2δt en+1 (xmin ) = 0 (3.18) en+1 (xmax ) = 0.
Taking the inner product of the equation (3.17) with en+1 , we obtain;
n+1 (3e − 4en + en−1 ) n+1 ,e = Den+1 , en+1 + λ I(Een ), en+1 + δ n+1 , en+1 , 2δt σ2 σ2 2 n+1 − λκ) en+1
+ (r − − (r + λ) en+1 2 = − en+1 x x ,e 2 2 +λ I(Een ), en+1 + δ n+1 , en+1 (3.19) by simplification of the above relation we have n+1
σ2 (r − − λκ) (3e − 4en + en−1 ) n+1 σ2 2 = − ,e
en+1 − en+1 2 x 2δt 2 σ2 +
(r −
σ2 2
− λκ)2 − 2(r + λ)σ 2 n+1 2
e + λ I(Een ), en+1 2 2σ
+ δ n+1 , en+1 σ2 2
− λκ)2 − 2(r + λ)σ 2 n+1 2 n n+1 ≤
e
+ λ I(Ee ), e 2 n+1 n+12σ + δ ,e (3.20) ≤ α en+1 2 + λ I(Een ), en+1 + δ n+1 , en+1 , (r −
2
(r− σ −λκ)2 −2(r+λ)σ 2
2 |. where α =| 2σ 2 Now by using the relation 2 (3a − 4b + c, a) = a 2 − b 2 + 2a−b 2 − 2b−c 2 + a−2b+c 2 , we have
1 n+1 2
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ α en+1 2 + λ I(Een ), en+1 4δt + δ n+1 , en+1 7
1 n+1 2
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ α en+1 2 + λCI Een
en+1
4δt + δ n+1
en+1
1 n+1 2 1 λCI
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ (α + + ) en+1 2 4δt 2 2 λCI 1
Een 2 + δ n+1 2 + 2 2 1 n+1 2 1 λCI
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ (α + + ) en+1 2 4δt 2 2 +λCI (4 en 2 + en−1 2 ) 1 + δ n+1 2 2
en+1 2 − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ δt[(2 + 4α + 2λCI ) en+1 2 + 2 δ n+1 2 +16λCI en 2 + 4λCI en−1 2 ] After summing up for n between 1 to m − 1, we get m 2
1 2
1
0 2
e − e − 2e − e
0 2
1 2
≤ δt[4λCI e + 20λCI e + 2
m
δ j 2
j=2
+(2 + 4α + 22λCI )
m−1
ej 2 + (2 + 4α + 2λCI ) em 2 ]
j=2
Now consider the δt sufficiently small such that 1 − δt(2 + 4α + 2λCI ) > 0 that is δt < 1 , then above relation imply that (2+4α+2λCI ) m 2
e
0 2
1 2
≤ C( e + e + δt
m
j 2
δ + δt
j=2
m−1
ej 2 )
j=2
≤ C( e0 2 + e1 2 + mδt max δ j 2 + δt
m−1
2≤j≤m
ej 2 ).
j=2
Since mδt ≤ T . Therefore, we have m 2
0 2
1 2
j 2
e ≤ C( e + e + max δ + δt 2≤j≤m
m−1
ej 2 ),
(3.21)
j=2
where C is a generic constant which is independent of mesh length. Now by applying the discrete Gronwall’s inequality, (for more detail see Lemma 3.1 in [17]) we get the result. We solve the resulting time semi discrete scheme by using radial basis function based finite difference method discussed in the next sub section. 8
3.2
Space discretizations
A function Φ : Rd → R is called radial provided there exists a univariate function φ : [0, ∞) → R such that Φ(x) = φ(r), where r = x and · is Euclidean norm on Rd . These functions can be broadly classified into two classes, infinitely smooth and piecewise smooth radial basis functions. The former include a shape parameter , and upon varying this parameter the radial function can vary sharp peak to very flat one. Classical choices of RBF are given in Table 1 with their order, where for any x ∈ R, the symbol x denotes as usual the smallest integer greater than or equal to x. The Gaussian and inverse multiquadric are positive definite functions where as multiquadric are conditionally positive definite functions of order m > 0. dk To derive local RBF-FD approximation of any linear differential operator L := dx k of order Table 1: Examples of radial basis functions and their order
R.B.F. Multiquadric (MQ) Inverse multiquadric (IMQ) Gaussian (GA)
φ(r), r > 0 (1 + ( r)2 )υ , υ > 0, υ ∈ /N 2 υ (1 + ( r) ) , υ < 0, υ ∈ /N −( r)2 e
Order of RBF m = υ m=0 m=0
k at a specific node point xi , in the discretized domain Ω := {x1 , x2 , . . . , xn } containing n number of nodes, consider any subset Ωi containing ni (<< n) nodes in the neighborhood of xi . In RBF-FD approach we are required to compute weights such that; L u(xi ) =
ni
wj u(xj ),
(3.22)
j=1
where ni and wj are known as the stencil size and differential weights respectively. For each node xi ∈ Ω we compute the weights wj on each local support Ωi . In traditional method, generally these nodes are equidistant and the weights are computed using classical polynomial interpolation. At the same time in radial basis function interpolation on randomly distributed nodes is used. Let s(x) be radial basis function interpolant that interpolates function u(x) at the interpolation points contained in Ωi . Then s(x) can be represented by s(x) =
ni
λj φ( x − xj ) +
j=1
l
γj pj (x)
(3.23)
j=1
where · is the Euclidian norm and {pj (x)}lj=1 denote basis of dm−1 , which is space of d-variate polynomials of total degree ≤ m − 1, where m is order of φ. The coefficients λj
9
and γj are evaluated by imposing the following conditions ni
s(xi ) = u(xi ), 1 ≤ i ≤ ni
(3.24)
1 ≤ k ≤ l.
(3.25)
λj pk (xj ) = 0,
j=1
Imposing conditions (3.24-3.25) on s(x) gives a linear system u |Ωi Φ P λ = P O γ O
(3.26)
where Φ := (φij )1≤i,j≤ni := (φ xi − xj )1≤i,j≤ni ∈ Rni ×ni , P := (pj (xi ))1≤i≤ni ,1≤j≤l ∈ Rni ×l . Suppose φ is conditionally positive definite function of order m on Rd and the points Ωi := {xi ∈ Rd ; i = 1, 2, . . . , ni } form (m − 1) unisolvent set of centers. Then the system (3.26) is uniquely solvable. The differential weights wj are calculated by enforcing the linear combination (3.22) should i be exact for RBFs {φ( x − xi )}ni=1 , centered at each of the stencil evaluation points. Hence by using RBFs and polynomial up to degree l, the weights are obtained by solving the following linear system i [L φ( xi − xj )nj=1 ] Φ P w = (3.27) P O v [L pj (xi )lj=1 ] note that the right hand side of equation (3.27) is nothing but [L φ( xi − x1 ) L φ( xi − x2 ) . . . L φ( xi − xni ) | L p1 (xi ) L p2 (xi ) . . . L pl (xi )] . It is obvious from the linear system (3.27) that its size is only ni + l, which is much smaller than the size n + l of global RBF collocation system. Thus the proposed method provides more stable system for a wide value of shape parameter . Algorithm to evaluate an European option n Then the approximated value of the solution unm = u(xm , tn ) denoted by Um can be obtained by following time stepping problem; for n = 0 : N − 1 if n ≤ 1 n+1 n Um −Um n+1 n = D Δ Um + λIΔ Um δt else 1 3 n+1 n n−1 n+1 n U − 2Um + 12 Um + λIΔ (EUm ) = DΔ U m δt 2 m end end where DΔ U is local RBF based space discretization of differential operator Du, obtained by using procedure described above and IΔ is numerical approximation of integral operator defined by equation (3.5). Since the option price problem have non smooth payoff function,there is possibility to have oscillation on the final time level. To overcome this problem, the solution on two time levels is calculated by implicit Euler method. 10
4
American options
The RBF based IMEX method discussed in previous section can be extend to solve LCP (2.10) for American option problem coupled with operator splitting method. The operator splitting method was introduced by Ikonen and Toivanen[15, 16] to evaluate the American put option. The basic idea behind the operator splitting method is the formulation with the auxiliary variable ψ such that ψ = ut − Lu. Now the reformulated LCP (2.10) as ⎧ ∂u − Lu = ψ, ⎪ ∂t ⎪ ⎨ (u(x, t) − G(x)) · ψ = 0, (4.1) u(x, t) − G(x) ≥ 0, ⎪ ⎪ ⎩ ψ ≥ 0, in the region Ω × [0, T ). Let us split the governing equation ut − Lu = ψ on the (n + 1)th time level into two discrete equations as
1 3 ˜ n+1 1 n−1 n n+1 n ˜ = Ψnm , (4.2) Um − 2Um + Um − DΔ Um + λIΔ (EUm ) δt 2 2
1 n−1 1 3 n+1 n n+1 n − DΔ U˜m Um − 2Um + Um + λIΔ (EUm ) = Ψn+1 (4.3) m . δt 2 2 Now the discrete problem for LCM (4.1) is to look for the pair (U n+1 , Ψn+1 ), which satisfy the discrete equations (4.2) and(4.3) and the constrains ⎧ n+1 Um ≥ G(xm ), ⎨ n+1 Ψm ≥ 0, (4.4) ⎩ n+1 n+1 Ψm (Um − G(xm )) = 0. n+1 The first step is to compute the intermediate approximation U˜m by solving the equations n (4.2) with known auxiliary term Ψm . Now the second step of the operator splitting method is to derive a relationship in (4.3) n+1 between Um and Ψn+1 m . To do this rewrite the equation (4.3) using equation (4.2) together n+1 with the constrains in (4.4) as a problem to find the pair (Um , Ψn+1 m ), such that n+1 ˜ n+1 −Um 3 Um n = Ψn+1 m − Ψm , 2 δt (4.5) n+1 n+1 Ψm (Um − G(xm )) = 0,
with the constrains n+1 Um ≥ G(xm ) and Ψn+1 ≥ 0. m
Now by solving the problems (4.5)-(4.6), we get n+1 ˜m U ) (G(xm ), Ψnm + 32 G(xm )− n+1 δt , Ψn+1 ) = (Um m 2δt n n+1 ˜ (Um − 3 Ψm , 0) 11
n+1 if U˜m − 2δt Ψnm ≤ G(xm ), 3 otherwise.
(4.6)
(4.7)
Thus, one can do the second step by solving discrete equation (4.2) with the updating formula (4.7). The above implicit method with three time levels requires the value on previous two 0 time levels. The pair (Um , Ψ0m ) on zeroth time level can be obtained by using initial condition 0 and assign value Ψm = 0 or one can use the algorithm discussed in [20]. To find the pair 1 1 (Um , Ψ1m ) at first level, the first step is to compute the intermediate value U˜m as follows: 1 0 − Um U˜m 1 0 + λIΔ Um = Ψ0m . − DΔ U˜m δt 1 The second step is to find the pair (Um , Ψ1m ) such that ˜1 1 if U˜m − δtΨ0m ≤ G(xm ), (G(xm ), Ψ0m + G(xmδt)−Um ) 1 1 (Um , Ψm ) = 1 (U˜m − δtΨ0m , 0) otherwise.
(4.8)
(4.9)
Algorithm to evaluate an American option for n = 0 : N − 1 if n ≤ 1 n+1 n ˜m U −Um n+1 n = DΔ U˜m + λIΔ Um + Ψnm δt n+1 n+1 Um = max(U˜m − δtΨnm , G(xm )) n+1 ˜ n+1 −Um = Ψnm + Um δt Ψn+1 m else
1 3 ˜ n+1 − 2U n + 1 U n−1 = DΔ U˜ n+1 + λIΔ (EU n ) + Ψn U 2 n+1 Um Ψn+1 m δt
end end
5
m
m
2
m
m
m
m
n+1 = max(U˜m − 2δt Ψnm , G(xm )) 3 n+1 ˜ n+1 −Um = Ψnm + 32 Um δt
Numerical simulation and discussion
In this section, we present several numerical experiments to illustrate the efficiency and accuracy of proposed method. We used the FFT algorithm to compute the product of a dense matrix and a column vector in the discrete integral operator in (3.5). Thus the implicit-explicit method with three time levels requires totally O(M N log2 (M )) operations, where N is the number of time steps and M is the number of spatial steps. Numerical results are presented for both call option and put option under Merton and Kou model using multiquadric (with υ = 12 and hence m = 1) radial basis functions. By keeping the shape parameter fixed, the computational error produced by the numerical schemes was measured against the value of analytical solution or reference solution at specific asset price. In all numerical experiments equispaced grid with ni = 3 is used. In all numerical experiments initial two step are calculated with IMEX-BDF1 method. We used MATLAB R2010 on a PC with 2.93GHz Intel(R) Core(TM)i3 CPU and 4.00GB RAM. Example 1. European option under Merton model with parameters σ = 0.15, r = 0.05, σJ = 0.45, μJ = −0.90, λ = 0.10, T = 0.25, K = 100. 12
Table 2: Numerical results for European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 9.283266 9.284893 9.285282 9.285384 9.285409 9.285416
= 90 Error 2.1520e-03 5.2475e-04 1.3655e-04 3.4508e-05 8.7377e-06 2.5525e-06
S = 100 Value Error 3.120808 2.8218e-02 3.142100 6.9262e-03 3.147302 1.7240e-03 3.148595 4.3054e-04 3.148918 1.0776e-04 3.148998 2.7723e-05
S = 110 Value Error 1.396156 5.0297e-03 1.399881 1.3053e-03 1.400852 3.3346e-04 1.401101 8.5126e-05 1.401165 2.1274e-05 1.401180 5.5654e-06
These parameters are also used by YongHoon Kwon, and Younhee Lee [19] for pricing European put option with truncated domain xmin = −1.5 and xmax = 1.5, and by Y. d’Halluin, P. A. Forsythy, and K. R. Vetzal [9] for European call option. The analytical solution is obtained by using infinite series solution of Merton [23]. The reference values for European put option are 9.285418 at S = 90, 3.149026 at S = 100, and 1.401186 at S = 110 and for European call option are 0.527638 at S = 90, 4.391246 at S = 100, and 12.643406 at S = 110. The numerical results for put option and call option are listed in Tables 2 and 3 respectively Table 3: Numerical results for European call options under the Merton model at different asset price with parameters as given in Example 1. S = 90 M N Value Error 129 25 0.525486 2.1520e-03 257 50 0.527113 5.2475e-04 513 100 0.527501 1.3655e-04 1025 200 0.527604 3.4508e-05 2049 400 0.527629 8.7377e-06 4097 800 0.527635 2.5525e-06
S = 100 Value Error 4.363028 2.8218e-02 4.384320 6.9262e-03 4.389522 1.7240e-03 4.390815 4.3054e-04 4.391138 1.0776e-04 4.391218 2.7723e-05
S = 110 Value Error 12.638376 5.0297e-03 12.642101 1.3053e-03 12.643072 3.3346e-04 12.643321 8.5126e-05 12.643385 2.1274e-05 12.643400 5.5654e-06
with shape parameter = 3.8. The option value and absolute error are given at different asset prices. From these tables, we can observe that the discretization scheme is convergent with quadratic rate, since the ratio in respective error is approximately four. The Greeks help to provide important measurements of an option position’s risks. Greeks are the quantities representing the sensitivity of the price of derivatives. Delta measures 13
Table 4: Numerical result for Delta of European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S = 90 Value Error -0.847817 1.1019e-03 -0.847023 3.0807e-04 -0.846794 7.8244e-05 -0.846735 1.9622e-05 -0.846720 4.9210e-06 -0.846717 1.3125e-06
S = 100 Value Error -0.358125 2.4621e-03 -0.356279 6.1556e-04 -0.355817 1.5392e-04 -0.355702 3.8482e-05 -0.355673 9.6151e-06 -0.355665 2.3771e-06
S = 110 Value Error -0.059047 9.4544e-04 -0.058334 2.3319e-04 -0.058158 5.6766e-05 -0.058115 1.3689e-05 -0.058105 3.4368e-06 -0.058102 8.1211e-07
the rate of change of option value with respect to changes in the underlying asset’s price, where as Gamma represents the rate of change in the delta with respect to changes in the underlying price. For European put option the reference values of Delta are −0.84671538 at S = 90, −0.35566306 at S = 100, and −0.05810123 at S = 110, and the values of Gamma are 0.03486014 at S = 90, 0.04882567 at S = 100, and 0.01212941 at S = 110. Since radial basis functions are sufficiently smooth it helps us to compute some Greeks like Delta and Gamma.It has been shown by Sanyasiraju et al.[26] and Wright et al.[32] that RBF-FD approximate first and second derivative of function with second order accuracy. Using the same parameter values which is employed for results in Table 2, we computed two Greeks viz Delta and Gamma at different asset prices and results obtained are reported in Table 4 and 5. From these tables we can observe that convergence rate of the proposed scheme is quadratic for computing Greeks as well. In Fig. 1 we have plotted the option values, the Delta and the Gamma. The plot shows that the option values and Greeks are very stable and no spurious oscillation occur at or around the strike price. These figures show that Greeks are efficiently evaluated using proposed method. Example 2. European call option under Merton model with parameters σ = 0.20, r = 0, σJ = 0.50, μJ = 0, λ = 0.10, T = 1.0, K = 1. The above parameters are also considered by Bhuruth et al.[25] with truncated domain xmin = −4.0 and xmax = 4.0. In their article, Bhuruth et al. provided a new radial basis function based differential quadrature approach (RBF-DQ). They first used RBF-DQ approach for the spatial discretization followed by exponential time integration scheme for time discretization. Using Merton’s closed form solution, the reference value 0.09413550 at strike price is calculated and the option value and absolute error at strike price are reported in Table 6 with shape parameter = 2.2. We observe that the results obtained with the present method have nice agreement with those given in[5, 25]. Example 3. European option under Kou model with parameters σ = 0.15, r = 0.05, η1 = 3.0465, η2 = 3.0775, λ = 0.10, p = 0.3445, T = 0.25, K = 100.. 14
70
60
Option Value
50
40
30
20
10
0
40
60
80
100
120
140
160
Asset Price
(a)
0
0.06
−0.1 0.05
−0.2
Gamma Value
Delta Value
−0.3 −0.4 −0.5 −0.6
0.04
0.03
0.02
−0.7 −0.8
0.01
−0.9 −1
40
60
80
100
120
140
0
160
Asset Price
40
60
80
100
120
140
160
Asset Price
(b)
(c)
Figure 1: European put option value, Delta and Gamma under the Merton model with parameters as provided in the example 1. These parameters are adopted from [19, 9]. To perform the numerical simulations, we choose the truncated domain as xmin = −1.5 and xmax = 1.5. The analytical solution of the model can be found in Kou[18]. The value of European call option at different asset price are evaluated with six digit after the decimal points by d’Halluin et al.[9]. The reference values for call option is 0.672677 at S = 90, 3.973479 at S = 100, and 11.794583 at S = 110. The reference price for put option is computed using put call parity. The numerical experiment 15
Table 5: Numerical result for Gamma of European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 0.034448 0.034758 0.034835 0.034854 0.034859 0.034860
= 90 Error 4.1186e-04 1.0225e-04 2.5208e-05 6.2791e-06 1.5713e-06 4.1016e-07
S = 100 Value Error 0.049816 9.9002e-04 0.049067 2.4124e-04 0.048886 5.9975e-05 0.048841 1.4976e-05 0.048829 3.7404e-06 0.048827 9.2079e-07
S = 110 Value Error 0.012071 5.8559e-05 0.012117 1.2733e-05 0.012126 3.3105e-06 0.012129 9.0800e-07 0.012129 2.2477e-07 0.012129 6.9091e-08
Table 6: Comparison of different methods for European call options under the Merton model at strike price with parameters as given in Example 2.
M N 65 5 129 10 257 20 513 40
CM[5] Value Error 0.09102 4.01e-03 0.09320 9.32e-04 0.09413 2.72e-04 0.09408 5.39e-05
RBF-DQ[25] Value Error 0.09182 2.32e-03 0.93565 5.70e-04 0.09399 1.38e-04 0.09410 3.02e-05
Present Value 0.092339 0.093865 0.094076 0.094121
Method Error 1.7963e-03 2.7096e-04 5.9464e-05 1.4369e-05
Table 7: Numerical results for European call options under the Kou model at different asset price with parameters as given in Example 3. S = 90 M N Value Error 129 25 0.670284 2.3932e-03 257 50 0.672105 5.7208e-04 513 100 0.672529 1.4783e-04 1025 200 0.672639 3.7620e-05 2049 400 0.672673 3.9368e-06 4097 800 0.672676 1.0867e-06
S = 100 Value Error 3.938910 3.4569e-02 3.964965 8.5136e-03 3.971358 2.1212e-03 3.972948 5.3113e-04 3.973363 1.1569e-04 3.973450 2.9465e-05
S = 110 Value Error 11.787842 6.7414e-03 11.792817 1.7660e-03 11.794129 4.5366e-04 11.794466 1.1692e-04 11.794561 2.2399e-05 11.794577 5.9222e-06
is carried out with shape parameter = 1.7. From Tables 7 and 8, we can observe that the point wise error at different asset price decrease, with convergence rate two. Finally in Fig. 16
70
60
Option Value
50
40
30
20
10
0
40
60
80
100
120
140
160
Asset Price
(a)
0
0.06
−0.1 0.05
−0.2
Gamma Value
Delta Value
−0.3 −0.4 −0.5 −0.6
0.04
0.03
0.02
−0.7 −0.8
0.01
−0.9 −1
40
60
80
100
120
140
0
160
Asset Price
40
60
80
100
120
140
160
Asset Price
(b)
(c)
Figure 2: European put option value, Delta and Gamma under the Kou model with parameters as provided in the example 3. 2, we have plotted the option value and Greeks for put option. These figures shows that Greeks are efficiently evaluated under Kou model as well using proposed method. Example 4. American put option under the Merton model with parameters σ = 0.15, r = 0.05, σJ = 0.45, μJ = −0.90, λ = 0.10, T = 0.25, K = 100. These parameters are consider by d’Halluin, Forsythy, and Labahn [8] and by YongHoon Kwon, and Younhee Lee[20] for Pricing American put option with truncated domain xmin = 17
Table 8: Numerical results for European put options under the Kou model at different asset price with parameters as given in Example 3.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 9.428064 9.429885 9.430309 9.430419 9.430453 9.430456
= 90 Error 2.3932e-03 5.7203e-04 1.4778e-04 3.7571e-05 3.8874e-06 1.0373e-06
S = 100 Value Error 2.696690 3.4569e-02 2.722745 8.5135e-03 2.729138 2.1212e-03 2.730728 5.3108e-04 2.731143 1.1564e-04 2.731230 2.9415e-05
S = 110 Value Error 0.545622 6.7413e-03 0.550597 1.7660e-03 0.551909 4.5361e-04 0.552246 1.1687e-04 0.552341 2.2349e-05 0.552357 5.8728e-06
−1.5 and xmax = 1.5. In [8] authors have developed an implicit discretization coupled with penalty method for pricing such American options. The authors in [20] developed the three time level IMEX method coupled with operator splitting method to solve such problems. The option value and error at different asset price are reported in Table 9 with shape parameter
= 3.8. To obtain error at different asset price in Table 9, we have computed the successive changes in the option value at each asset price. From the above table, we observe that Table 9: Numerical results for American put options under the Merton model at different asset price with parameters as given in Example 4.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S = 90 Value Error 9.998080 10.000259 2.1789e-03 10.003020 2.7610e-03 10.003934 9.1430e-04 10.003835 9.9612e-05 10.003866 3.1405e-05
S = 100 Value Error 3.207266 3.232387 2.5120e-02 3.238920 6.5333e-03 3.240634 1.7132e-03 3.241085 4.5131e-04 3.241207 1.2209e-05
S = 110 Value Error 1.414138 1.418174 4.0361e-03 1.419327 1.1528e-03 1.419660 3.3244e-04 1.419759 9.9174e-05 1.419790 3.0951e-05
point wise errors generally converge with approximately second order accuracy, while near the optimal exercise boundary convergence is less accurate. The above phenomenon is also discussed in [8, 20]. In Fig. 3, we provide the plot of both option values and the early exercise boundary for American put option under the Merton model. Example 5. American put option under the Kou model with parameters σ = 0.15, r = 0.05, η1 = 3.0465, η2 = 3.0775, λ = 0.10, p = 0.3445, T = 0.25, K = 100.. These parameter are also adopted by YongHoon Kwon, and Younhee Lee [20] with same truncated domain and shape parameter. The reference values at different asset prices is 18
70
100
American Payoff 60
98
Holding region
96
Stock Price
Option Value
50
40
30
94
92
20 90
10
0
88
40
60
80
100
120
140
160
Exercise region 0
0.05
0.1
0.15
0.2
0.25
Time (years)
Asset Price
(a)
(b)
Figure 3: American put option value, and optimal early exercise boundary under the Merton model with parameters as provided in the example 4. given by Toivanen [29] are 10.005071 at S = 90, 2.807879 at S = 100, and 0.561876 at S = 110. In Table 10, we get the same scenario as in the case of Merton model. We are Table 10: Numerical results for American put options under the Kou model at different asset price with parameters as given in Example 5. S = 90 M N Value Error 129 25 10.002131 257 50 10.000735 1.3951e-03 513 100 10.005015 4.2794e-03 1025 200 10.005055 3.9945e-05 2049 400 10.005088 3.3276e-05 4097 800 10.005098 9.4650e-06
S = 100 Value Error 2.777966 2.799956 2.1990e-02 2.805766 5.8098e-03 2.807307 1.5414e-03 2.807723 4.1541e-04 2.807834 1.1140e-04
S = 110 Value Error 0.558356 0.560726 2.3706e-03 0.561502 7.7560e-04 0.561750 2.4808e-04 0.561834 8.3698e-05 0.561861 2.7201e-05
getting approximately second order convergence rate, while we have less accuracy near by the early exercise boundary. Finally in Fig. 4, we have plotted the option value and early exercise boundary for American put option under Kou model. Since we are using backward difference method of order two for time semi discretization, and problem have non smooth initial condition, we observe that the increasing the value of ni 19
70
100
American Payoff 60
98
Holding region
96
Stock Price
Option Value
50
40
30
94
92
20 90
10
0
88
40
60
80
100
120
140
160
Exercise region 0
0.05
0.1
0.15
0.2
0.25
Time (years)
Asset Price
(a)
(b)
Figure 4: American put option value, and optimal early exercise boundary under the Kou model with parameters as provided in the example 5. does not seems comparably efficient. Finally cpu times(in seconds) of the method with Merton model for European put option and American put option (example 1 and 4) and with Kou model for for European put option and American put option (example 3 and 5) are given in Table 11. Table 11: CPU times of the proposed methods for different experiment setup.
M
N
129 257 513 1025 2049 4097
25 50 100 200 400 800
Merton model European American option option 0.2260 0.2766 0.0359 0.0617 0.0582 0.0878 0.1602 0.2369 0.5644 0.7946 2.0307 2.7810
20
Kou model European American option option 0.2192 0.2765 0.0512 0.0635 0.0648 0.0913 0.1769 0.2492 0.6078 0.8194 2.1449 2.9971
6
Conclusion
In this article, a RBF based method for solving partial integro-differential equations, which arises when underlying asset follows the jump diffusion process has been presented. After time semi discretization of governing equation, the resulting linear differential equation is solved using RBF based finite difference method, the integral term involved in PIDE is approximated with composite trapezoidal rule. The numerical simulation with European put/call and American put option under Merton and Kou model is carried out. The results are compared against the FD, FEM, and RBF-DQ method available in literature and shown that the proposed method is more accurate and efficient. The proposed method approximate not only the option value but also some of its “Greeks” namely “Delta” and “Gamma” very efficiently.
Acknowledgements The authors would like to thanks the anonymous referees for their valuable comments and suggestion which helped to improve the presentation of the work. This research work of authors Alpesh Kumar, Lok Pati Tripathi is supported by Council of Scientific and Industrial Research, India with grant no. 10-2(5)/2007(ii)-E.U.II and 09/092(0715)/2009-EMR-I respectively.
References [1] A. Almendral, C. W. Oosterlee, Numerical valuation of options with jumps in the underlying, Applied Numerical Mathematics 53 (1) (2005) 1–18. [2] L. Andersen, J. Andreasen, Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing, Review of Derivatives Research 4 (3) (2000) 231–262. [3] M. Briani, R. Natalini, G. Russo, Implicit–explicit numerical schemes for jump–diffusion processes, Calcolo 44 (1) (2007) 33–57. [4] R. Brummelhuis, R. T. Chan, A Radial Basis Function scheme for option pricing in exponential l´evy models, Applied Mathematical Finance 21 (3) (2014) 238–269. [5] P. Carr, A. Mayo, On the numerical evaluation of option prices in jump diffusion processes, European Journal of Finance 13 (4) (2007) 353–372. [6] R. T. L. Chan, S. Hubbert, Options pricing under the one-dimensional jump-diffusion model using the radial basis function interpolation scheme, Review of Derivatives Research 17 (2) (2014) 161–189. [7] G. Chandhini, Y. Sanyasiraju, Local RBF-FD solutions for steady convection–diffusion problems, International Journal for Numerical Methods in Engineering 72 (3) (2007) 352–378.
21
[8] Y. dHalluin, P. A. Forsyth, G. Labahn, A penalty method for american options with jump diffusion processes, Numerische Mathematik 97 (2) (2004) 321–352. [9] Y. dHalluin, P. A. Forsyth, K. R. Vetzal, Robust numerical methods for contingent claims under jump diffusion processes, IMA Journal of Numerical Analysis 25 (1) (2005) 87–112. [10] G. E. Fasshauer, A. Q. M. Khaliq, D. A. Voss, Using meshfree approximation for multiasset American options, Journal of the Chinese Institute of Engineers 27 (4) (2004) 563–571. [11] B. Fornberg, G. Wright, E. Larsson, Some observations regarding interpolants in the limit of flat radial basis functions, Computers & mathematics with applications 47 (1) (2004) 37–55. [12] A. Golbabai, D. Ahmadian, M. Milev, Radial basis functions with application to finance: American put option under jump diffusion, Mathematical and Computer Modelling 55 (3) (2012) 1354–1362. [13] Y. C. Hon, A quasi-radial basis functions method for American options pricing, Computers & Mathematics with Applications 43 (3) (2002) 513–524. [14] Y. C. Hon, X. Z. Mao, A radial basis function method for solving options pricing models, Journal of Financial Engineering 8 (1999) 31–50. [15] S. Ikonen, J. Toivanen, Operator splitting methods for American option pricing, Applied Mathematics Letters 17 (7) (2004) 809–814. [16] S. Ikonen, J. Toivanen, Efficient numerical methods for pricing American options under stochastic volatility, Numerical Methods for Partial Differential Equations 24 (1) (2008) 104–126. [17] M. K. Kadalbajoo, L. P. Tripathi, A. Kumar, Second order accurate IMEX methods for option pricing under merton and kou jump-diffusion models, Journal of Scientific Computing (2015) 1–46. [18] S. G. Kou, A jump-diffusion model for option pricing, Management science 48 (8) (2002) 1086–1101. [19] Y. Kwon, Y. Lee, A second-order finite difference method for option pricing under jump-diffusion models, SIAM Journal on Numerical Analysis 49 (6) (2011) 2598–2617. [20] Y. Kwon, Y. Lee, A second-order tridiagonal method for American options under jumpdiffusion models, SIAM Journal on Scientific Computing 33 (4) (2011) 1860–1872. [21] C. K. Lee, X. Liu, S. C. Fan, Local multiquadric approximation for solving boundary value problems, Computational Mechanics 30 (5-6) (2003) 396–409. 22
[22] X. Liu, K. Tai, Point interpolation collocation method for the solution of partial differential equations, Engineering analysis with boundary elements 30 (7) (2006) 598–609. [23] R. C. Merton, Option pricing when underlying stock returns are discontinuous, Journal of financial economics 3 (1) (1976) 125–144. [24] E. Ngounda, K. C. Patidar, E. Pindza, Contour integral method for European options with jumps, Communications in Nonlinear Science and Numerical Simulation. [25] A. Saib, D. Y. Tangman, M. Bhuruth, A new radial basis functions method for pricing American options under merton’s jump-diffusion model, International Journal of Computer Mathematics 89 (9) (2012) 1164–1185. [26] Y. Sanyasiraju, G. Chandhini, Local radial basis function based gridfree scheme for unsteady incompressible viscous flows, Journal of Computational Physics 227 (20) (2008) 8922–8948. [27] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible navier–stokes equations, Computer Methods in Applied Mechanics and Engineering 192 (7) (2003) 941–954. [28] D. Tangman, A. Gopaul, M. Bhuruth, Exponential time integration and chebychev discretisation schemes for fast pricing of options, Applied Numerical Mathematics 58 (9) (2008) 1309–1319. [29] J. Toivanen, Numerical valuation of European and American options under kou’s jumpdiffusion model, SIAM Journal on Scientific Computing 30 (4) (2008) 1949–1970. [30] A. I. Tolstykh, On using RBF-based differencing formulas for unstructured and mixed structured-unstructured grid calculations, in: Proceedings of the 16th IMACS World Congress, Lausanne, 2000. [31] A. I. Tolstykh, D. A. Shirobokov, On using radial basis functions in a “finite difference mode” with applications to elasticity problems, Computational Mechanics 33 (1) (2003) 68–79. [32] G. B. Wright, B. Fornberg, Scattered node compact finite difference-type formulas generated from radial basis functions, Journal of Computational Physics 212 (1) (2006) 99–123. [33] K. Zhang, S. Wang, Pricing options under jump diffusion processes with fitted finite volume method, Applied Mathematics and Computation 201 (1) (2008) 398–413.
23
Mohan K. Kadalbajoo, Alpesh Kumar, Lok Pati Tripathi
PII: DOI: Reference:
S0168-9274(16)30156-8 http://dx.doi.org/10.1016/j.apnum.2016.08.006 APNUM 3077
To appear in:
Applied Numerical Mathematics
Received date: Revised date: Accepted date:
27 March 2015 17 June 2016 19 August 2016
Please cite this article in press as: M.K. Kadalbajoo et al., A radial basis function based implicit-explicit method for option pricing under jump-diffusion models, Appl. Numer. Math. (2016), http://dx.doi.org/10.1016/j.apnum.2016.08.006
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A radial basis function based implicit-explicit method for option pricing under jump-diffusion models Mohan K. Kadalbajoo, Alpesh Kumar∗, Lok Pati Tripathi Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur 208016, India. Abstract In this article, we present a radial basis function based implicit explicit numerical method to solve the partial integro-differential equation which describe the nature of the option price under jump diffusion model. The governing equation is time semi discrtized by using the implicit-explicit backward difference method of order two (IMEXBDF2) followed by radial basis function based finite difference (RBF-FD) method. The numerical scheme derived for European option is extended for American option by using operator splitting method. Numerical results for put and call option under Merton and Kou models are given to illustrate the efficiency and accuracy of the present method. The stability of time semi discretized scheme is also proved.
Keywords: Radial basis function; Finite difference, Option pricing; Jump-diffusion models; Partial integro-differential equation.
1
Introduction
There is evidence to suggest that the Black Scholes model for stock price behavior does not always model real stock price behavior. Jump can appear at a random time and these jumps can not be captured by the log normal distribution characteristic of the stock price in the Black Scholes model. To overcome the above shortcoming, several models have been proposed in the literature. Among these, the jump diffusion models introduced by Merton[23] and Kou[18] are of the most widely used models. Merton proposed a log-normally distributed process for the jump-amplitudes, whereas Kou suggested log-double-exponentially distributed process. The valuation of option under jump diffusion process requires the solution of a partial integro differential equation containing a non local integral term. There are several numerical methods available in the literature to approximate the above equation. In [1], Almendral and Osterlee presented an implicit second order accurate time discretization with finite difference and finite element spatial discretization on uniform grid. Andersen et al.[2] proposed ∗
corresponding author, Email: [email protected](Alpesh Kumar)
1
an unconditionally stable alternating direction implicit method for its solution. Song Wang et al.[33] developed a fitted finite volume method for jump diffusion process. Their method is based on fitted finite volume method spatial discretization and Crank Nicolson scheme for temporal discretization. More recently, Patidar et al.[24] developed an efficient method for pricing Merton jump diffusion option, combining the spectral domain decomposition method and the Laplace transform method. The scheme proposed by d Halluin et al.[9] required to use an iterative procedure to solve discrete equations. The main difficulty with implicit scheme is due to containing non local integral term in governing equation, which leads to a dense discretization matrix where as fully explicit scheme imposed stability restriction on it. An approach based on implicit-explicit schemes in which integral term is treated explicitly has been proposed by YongHoon Kwon et al.[19] and Briani et al.[3]. More recently Tangman et al.[28] introduced a new scheme called exponential time integration (ETI) scheme to solve the PIDE. In ETI scheme, the time direction of PIDE is directly tackled by a ’one step’ formula, which means temporal discretization is not required. Tangman et al.[28] used the central difference approach with ETI to provide very efficient and second order accurate result. Recently, a new method based on radial basis function (RBF) for approximation of spatial derivative in option pricing equation is under going active research. Application of RBF in one dimension European and American options is given by Hon et al.[13, 14]. Fasshauer et al.[10] solved American option pricing model using penalty method. Golbabai et al.[12], developed an algorithm based on global collocation for jump diffusion process. Bhuruth et al.[25] proposed a radial basis function based differential quadrature rule for spatial descretization with exponential time integration to solve jump diffusion model. In more recent work, Chan et al.[4, 6] used new RBF called cubic spline as basis function to solve PIDE, and show that their scheme is second order accurate. It was recognized that standard approach to solving the radial basis function collocation problem has been ill conditioned due to use of collocation in global sense. Recently many strategies have been developed in the literature to avoid these problems, such as local RBF approach by Lee et al.[21], radial point interpolation method proposed by Liu et al.[22], Shu et al.[27] proposed a local radial basis function-based differential quadrature method and used it to solve two-dimensional incompressible Navier-Stokes equations. Tolstykh[30], Tolstykh and Shirobokov[31], Wright et al.[11, 32] proposed radial basis function finite difference method, the idea is to use radial basis functions with a local collocation as in finite difference mode thereby reducing number of nodes and hence producing a sparse matrix. This technique is further extended by Sanyasiraju et al.[7, 26] for convection diffusion type equations. However these methods have not been extended to solve partial integro differential equation yet. In the present work, we have extended the localization concept proposed by Wright and Fornberg, to solve jump diffusion models. The governing equations are discretized by a three level implicit and explicit time scheme followed by RBF based finite difference method. The paper is organized as follow. In section 2, mathematical models for pricing option with jump diffusion process are given in terms of partial integro-differential equations and provide a brief review of both the Merton and Kou jump diffusion models. Section 3 deals with 2
the construction of three time level implicit explicit scheme to discretize the jump diffusion model. The time semi discrete equation is coupled with radial basis function based finite difference method for spatial discretization. Section 4, provide extension of proposed method for pricing American option by utilizing concept of operator splitting method. In section 5, we give some numerical results for Merton and Kou model and a comparison of the accuracy of our solution with finite difference and finite element method for both American and European options. Finally the paper ends with some conclusive remark in section 6.
2
The mathematical model
In this section, we give brief discussion about the mathematical model for option with jump diffusion process. Consider an asset with the asset price S, then the movement of stock price is modeled by the following stochastic differential equation dS = (ν − κλ)dτ + σdZ + (η − 1)dq S
(2.1)
where, ν is drift rate, τ as the time to maturity, σ represent the constant volatility, dZ is an increment of standard Gauss-Wiener process. The term λ is the intensity of the independent Poisson process dq with 0 with probability 1 − λdτ , dq = 1 with probability λdτ . The expected relative jump size E(η − 1) is denoted by κ, where E[·] is the expectation operator and η − 1 is a impulse function producing jump from S to Sη. Let V (S, τ ) represent the value of a contingent claim that depend on the underlying asset price S with current time τ . Then V (S, τ ) satisfy following backward partial integro differential equation ∞ 1 2 2 ∂ 2V ∂V ∂V + σ S − (r + λ)V + λ + (r − λκ)S V (Sη)g(η)dη = 0, (2.2) ∂τ 2 ∂S 2 ∂S 0 for (S, τ ) ∈ (0, ∞) × (0, T ], where, r is risk free interest rate and g(η) is probability density ∞ function of the jump with amplitude η with properties that ∀η, g(η) ≥ 0 and 0 g(η)dη = 1. The value of V at the expiry date is given by, V (S, T ) = G(S), S ∈ (0, ∞),
(2.3)
where G(S) is the pay-off function for the option contract. Under Merton’s model g(η) is given by the log-normal density 2
(ln η−μJ ) − 1 2σ 2 J g(η) := √ . e 2πσJ η
3
(2.4)
σ2
In this case κ := E(η − 1) = exp(μJ + 2J ) − 1, where μJ and σJ2 are the mean and the variance of jump in return. Under a modified version of Kou’s jump-diffusion model g(η) is the following log-double-exponential density g(η) :=
1 pη1 e−η1 ln(η) H(ln(η)) + qη2 eη2 ln(η) H(− ln(η)) , η
(2.5)
where η1 > 1, η2 > 0, p > 0, q = 1 − p, and H(·) is the Heaviside function. We can show 1 2 that κ : ηpη + ηqη − 1. 1 −1 1 +1 S By using change of variable x = ln( K ), y = ln(η), where K is strike price and letting t = T − τ , computation of the option value requires solving the PIDE ∞ ∂u 1 2 ∂ 2u ∂u σ2 = σ − λκ) − (r + λ)u + λ + (r − u(y, t)f (y − x)dy, (2.6) ∂t 2 ∂x2 2 ∂x −∞ where u(x, t) := V (Kex , T − t) and f (y) = g(ey )ey , under the above transformation the function f (y) can be written as; ⎧ 2 J) ⎨ 1 − (y−μ 2σ 2 √ J e Merton’s model, 2πσJ (2.7) f (y) := ⎩ −η1 y η2 y H(y) + qη2 e H(−y)) Kou’s model. (pη1 e The initial condition and asymptotic behavior of the price of a European put option is described by u(x, 0) = max(K − Kex , 0) Ke−rt − Kex x → −∞, u(x, t) = 0 x → ∞.
(2.8) (2.9)
Other types of initial and boundary conditions can be suitably defines for different types options. In contrast, the American option can be exercised at any time up to the maturity date and can be formulated as the linear complementarity problem(LCP) of the form ⎧ ∂u − Lu ≥ 0, ⎨ ∂t u(x, t) − G(x) ≥ 0, (2.10) ⎩ ∂u ( ∂t − Lu)(u(x, t) − G(x)) = 0, for (x, t) ∈ (0, ∞) × [0, T ) together with the boundary conditions K − Kex x → −∞, u(x, t) = 0 x → ∞, where L is the partial integro differential operator on the right side of (2.6). 4
(2.11)
3
Discretizations
In this section we shall construct an implicit explicit time semi discretization for following PIDE on truncated domain Ω := (xmin , xmax ), ∂u(x, t) = Du(x, t) + λIu(x, t), (x, t) ∈ Ω × [0, T ) ∂t ¯ u(x, 0) = G(x), x ∈ Ω
(3.1) (3.2)
where D is differential operator and I is integral operator of L such that σ2 ∂u 1 2 ∂ 2u + (r − σ − λκ) − (r + λ)u, Du(x, t) = 2 2 ∂x 2 ∂x
(3.3)
u(y, t)f (y − x)dy.
(3.4)
Iu(x, t) =
∞
−∞
In order to approximate the integral operator Iu in (3.4) numerically, we first divide the integral into two parts on Ω and on Ωc := R \ Ω. Now the the integral operator can be split as Iu(x, t) = u(y, t)f (y − x)dy + R(t, x). (3.5) Ω
By using asymptotic behavior of the option the integral over R \ Ω is given by (Ke−rt − Key )f (y − x)dy. R(t, x) :=
(3.6)
Ωc
For the European put option the expression R(t, x) is given as; ⎧
σ2 xmin −x−μJ −σJ2 ⎨ −rt x+μJ + 2J J N Ke N xmin −x−μ − Ke σJ σJ R(t, x) := ⎩K(1 − p)e−rt+η2 (xmin −x) − K(1 − p) η2 e−η2 x+(η2 +1)xmin η2 +1
Merton’s ,
(3.7)
Kou’s.
In the case of American put option above expression can be reformulated as; ⎧
σ2 xmin −x−μJ −σJ2 ⎨ x+μJ + 2J J − Ke Merton’s, N KN xmin −x−μ σJ σJ R(t, x) := ⎩K(1 − p)eη2 (xmin −x) − K(1 − p) η2 e−η2 x+(η2 +1)xmin Kou’s.
(3.8)
η2 +1
where N (·) is the cumulative normal distribution. We use composite trapezoidal rule to approximate the integral term on truncated domain Ω.
3.1
Time discretizations
Let {0 = t0 < t1 < . . . < tN = T ; tn − tn−1 = δt, 1 ≤ n ≤ N } be a partition of the interval [0, T ]. Let us use the notation un := u(x, tn ) then equation (3.1) will be discretized by following implicit explicit scheme,
1 3 n+1 1 n−1 n u = Dun+1 + λI(Eun ), n ≥ 1, − 2u + u (3.9) δt 2 2 5
together with boundary condition un+1 (xmin ) = Ke−rtn+1 − Kexmin un+1 (xmax ) = 0,
(3.10)
where Eun = 2un − un−1 . The above discretization method is called IMEX-BDF2 method with three time levels. In order to use the proposed method we need two initial values on the zeroth and first time level. The value u0 at the zeroth time level is given by initial condition of the model problem, and the value u1 can be obtained by applying the implicit-explicit backward difference method of order one, un+1 − un = Dun+1 + λIun , (3.11) δt together with boundary condition un+1 (xmin ) = Ke−rtn+1 − Kexmin un+1 (xmax ) = 0.
(3.12)
Now we prove the stability of the time discrete scheme (3.9). It can be easily shown that for all u(·, t) ∈ L2 (Ω), t ∈ (0, T ) such that u(x, t) : (x, t) ∈ Ω × [0, T ], u¯(x, t) = 0 : (x, t) ∈ Ωc × [0, T ], the integral operator satisfies the condition
I u¯(·, t) ≤ CI u(·, t) , (3.13) for some constant CI independent of t, where v := ( Ω | v(x) |2 dx)1/2 . Let us suppose that un be the solution of the equation (3.9) and u˜n be the solution of perturbed equation
1 3 n+1 1 n−1 n = D˜ un+1 + λI(E u˜n ) + δ n+1 , n ≥ 1, u˜ − 2˜ u + u˜ (3.14) δt 2 2 together with boundary condition u˜n+1 (xmin ) = Ke−rtn+1 − Kexmin u˜n+1 (xmax ) = 0.
(3.15)
The δ j is the perturbation in the data at any time level tj and what we will examine the effect of these perturbations on the error as we proceed in time. We will show that the error at any time level is bounded by the initial errors at initial time levels(t0 & t1 ; for BDF2) and is a bounded multiple of the introduced perturbations; which establishes stability in time discretization. Now define the error term en := u˜n − un . Since, we are treating boundary conditions exactly in the equations (3.9-3.10) and perturbed equations (3.14-3.15). Hence the error term en will be zero at the boundary points. 6
Theorem 1. (L2 -stability)For sufficiently small δt such that δt <
em 2 ≤ C( e0 2 + e1 2 + max δ j 2 ), 2≤j≤m
1 , (2+4α+2λCI )
∀2≤m≤
T , δt
we have (3.16)
2
(r− σ −λκ)2 −2(r+λ)σ 2
2 where α =| | and C is generic constant depending on the parameter CI , 2σ 2 r, σ, λ, and T . Proof. Relations (3.9-3.10) and (3.14-3.15) implies that the error term en satisfies the following relations (3en+1 − 4en + en−1 ) = Den+1 + λI(Een ) + δ n+1 , (3.17) 2δt en+1 (xmin ) = 0 (3.18) en+1 (xmax ) = 0.
Taking the inner product of the equation (3.17) with en+1 , we obtain;
n+1 (3e − 4en + en−1 ) n+1 ,e = Den+1 , en+1 + λ I(Een ), en+1 + δ n+1 , en+1 , 2δt σ2 σ2 2 n+1 − λκ) en+1
+ (r − − (r + λ) en+1 2 = − en+1 x x ,e 2 2 +λ I(Een ), en+1 + δ n+1 , en+1 (3.19) by simplification of the above relation we have n+1
σ2 (r − − λκ) (3e − 4en + en−1 ) n+1 σ2 2 = − ,e
en+1 − en+1 2 x 2δt 2 σ2 +
(r −
σ2 2
− λκ)2 − 2(r + λ)σ 2 n+1 2
e + λ I(Een ), en+1 2 2σ
+ δ n+1 , en+1 σ2 2
− λκ)2 − 2(r + λ)σ 2 n+1 2 n n+1 ≤
e
+ λ I(Ee ), e 2 n+1 n+12σ + δ ,e (3.20) ≤ α en+1 2 + λ I(Een ), en+1 + δ n+1 , en+1 , (r −
2
(r− σ −λκ)2 −2(r+λ)σ 2
2 |. where α =| 2σ 2 Now by using the relation 2 (3a − 4b + c, a) = a 2 − b 2 + 2a−b 2 − 2b−c 2 + a−2b+c 2 , we have
1 n+1 2
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ α en+1 2 + λ I(Een ), en+1 4δt + δ n+1 , en+1 7
1 n+1 2
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ α en+1 2 + λCI Een
en+1
4δt + δ n+1
en+1
1 n+1 2 1 λCI
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ (α + + ) en+1 2 4δt 2 2 λCI 1
Een 2 + δ n+1 2 + 2 2 1 n+1 2 1 λCI
e − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ (α + + ) en+1 2 4δt 2 2 +λCI (4 en 2 + en−1 2 ) 1 + δ n+1 2 2
en+1 2 − en 2 + 2en+1 − en 2 − 2en − en−1 2 ≤ δt[(2 + 4α + 2λCI ) en+1 2 + 2 δ n+1 2 +16λCI en 2 + 4λCI en−1 2 ] After summing up for n between 1 to m − 1, we get m 2
1 2
1
0 2
e − e − 2e − e
0 2
1 2
≤ δt[4λCI e + 20λCI e + 2
m
δ j 2
j=2
+(2 + 4α + 22λCI )
m−1
ej 2 + (2 + 4α + 2λCI ) em 2 ]
j=2
Now consider the δt sufficiently small such that 1 − δt(2 + 4α + 2λCI ) > 0 that is δt < 1 , then above relation imply that (2+4α+2λCI ) m 2
e
0 2
1 2
≤ C( e + e + δt
m
j 2
δ + δt
j=2
m−1
ej 2 )
j=2
≤ C( e0 2 + e1 2 + mδt max δ j 2 + δt
m−1
2≤j≤m
ej 2 ).
j=2
Since mδt ≤ T . Therefore, we have m 2
0 2
1 2
j 2
e ≤ C( e + e + max δ + δt 2≤j≤m
m−1
ej 2 ),
(3.21)
j=2
where C is a generic constant which is independent of mesh length. Now by applying the discrete Gronwall’s inequality, (for more detail see Lemma 3.1 in [17]) we get the result. We solve the resulting time semi discrete scheme by using radial basis function based finite difference method discussed in the next sub section. 8
3.2
Space discretizations
A function Φ : Rd → R is called radial provided there exists a univariate function φ : [0, ∞) → R such that Φ(x) = φ(r), where r = x and · is Euclidean norm on Rd . These functions can be broadly classified into two classes, infinitely smooth and piecewise smooth radial basis functions. The former include a shape parameter , and upon varying this parameter the radial function can vary sharp peak to very flat one. Classical choices of RBF are given in Table 1 with their order, where for any x ∈ R, the symbol x denotes as usual the smallest integer greater than or equal to x. The Gaussian and inverse multiquadric are positive definite functions where as multiquadric are conditionally positive definite functions of order m > 0. dk To derive local RBF-FD approximation of any linear differential operator L := dx k of order Table 1: Examples of radial basis functions and their order
R.B.F. Multiquadric (MQ) Inverse multiquadric (IMQ) Gaussian (GA)
φ(r), r > 0 (1 + ( r)2 )υ , υ > 0, υ ∈ /N 2 υ (1 + ( r) ) , υ < 0, υ ∈ /N −( r)2 e
Order of RBF m = υ m=0 m=0
k at a specific node point xi , in the discretized domain Ω := {x1 , x2 , . . . , xn } containing n number of nodes, consider any subset Ωi containing ni (<< n) nodes in the neighborhood of xi . In RBF-FD approach we are required to compute weights such that; L u(xi ) =
ni
wj u(xj ),
(3.22)
j=1
where ni and wj are known as the stencil size and differential weights respectively. For each node xi ∈ Ω we compute the weights wj on each local support Ωi . In traditional method, generally these nodes are equidistant and the weights are computed using classical polynomial interpolation. At the same time in radial basis function interpolation on randomly distributed nodes is used. Let s(x) be radial basis function interpolant that interpolates function u(x) at the interpolation points contained in Ωi . Then s(x) can be represented by s(x) =
ni
λj φ( x − xj ) +
j=1
l
γj pj (x)
(3.23)
j=1
where · is the Euclidian norm and {pj (x)}lj=1 denote basis of dm−1 , which is space of d-variate polynomials of total degree ≤ m − 1, where m is order of φ. The coefficients λj
9
and γj are evaluated by imposing the following conditions ni
s(xi ) = u(xi ), 1 ≤ i ≤ ni
(3.24)
1 ≤ k ≤ l.
(3.25)
λj pk (xj ) = 0,
j=1
Imposing conditions (3.24-3.25) on s(x) gives a linear system u |Ωi Φ P λ = P O γ O
(3.26)
where Φ := (φij )1≤i,j≤ni := (φ xi − xj )1≤i,j≤ni ∈ Rni ×ni , P := (pj (xi ))1≤i≤ni ,1≤j≤l ∈ Rni ×l . Suppose φ is conditionally positive definite function of order m on Rd and the points Ωi := {xi ∈ Rd ; i = 1, 2, . . . , ni } form (m − 1) unisolvent set of centers. Then the system (3.26) is uniquely solvable. The differential weights wj are calculated by enforcing the linear combination (3.22) should i be exact for RBFs {φ( x − xi )}ni=1 , centered at each of the stencil evaluation points. Hence by using RBFs and polynomial up to degree l, the weights are obtained by solving the following linear system i [L φ( xi − xj )nj=1 ] Φ P w = (3.27) P O v [L pj (xi )lj=1 ] note that the right hand side of equation (3.27) is nothing but [L φ( xi − x1 ) L φ( xi − x2 ) . . . L φ( xi − xni ) | L p1 (xi ) L p2 (xi ) . . . L pl (xi )] . It is obvious from the linear system (3.27) that its size is only ni + l, which is much smaller than the size n + l of global RBF collocation system. Thus the proposed method provides more stable system for a wide value of shape parameter . Algorithm to evaluate an European option n Then the approximated value of the solution unm = u(xm , tn ) denoted by Um can be obtained by following time stepping problem; for n = 0 : N − 1 if n ≤ 1 n+1 n Um −Um n+1 n = D Δ Um + λIΔ Um δt else 1 3 n+1 n n−1 n+1 n U − 2Um + 12 Um + λIΔ (EUm ) = DΔ U m δt 2 m end end where DΔ U is local RBF based space discretization of differential operator Du, obtained by using procedure described above and IΔ is numerical approximation of integral operator defined by equation (3.5). Since the option price problem have non smooth payoff function,there is possibility to have oscillation on the final time level. To overcome this problem, the solution on two time levels is calculated by implicit Euler method. 10
4
American options
The RBF based IMEX method discussed in previous section can be extend to solve LCP (2.10) for American option problem coupled with operator splitting method. The operator splitting method was introduced by Ikonen and Toivanen[15, 16] to evaluate the American put option. The basic idea behind the operator splitting method is the formulation with the auxiliary variable ψ such that ψ = ut − Lu. Now the reformulated LCP (2.10) as ⎧ ∂u − Lu = ψ, ⎪ ∂t ⎪ ⎨ (u(x, t) − G(x)) · ψ = 0, (4.1) u(x, t) − G(x) ≥ 0, ⎪ ⎪ ⎩ ψ ≥ 0, in the region Ω × [0, T ). Let us split the governing equation ut − Lu = ψ on the (n + 1)th time level into two discrete equations as
1 3 ˜ n+1 1 n−1 n n+1 n ˜ = Ψnm , (4.2) Um − 2Um + Um − DΔ Um + λIΔ (EUm ) δt 2 2
1 n−1 1 3 n+1 n n+1 n − DΔ U˜m Um − 2Um + Um + λIΔ (EUm ) = Ψn+1 (4.3) m . δt 2 2 Now the discrete problem for LCM (4.1) is to look for the pair (U n+1 , Ψn+1 ), which satisfy the discrete equations (4.2) and(4.3) and the constrains ⎧ n+1 Um ≥ G(xm ), ⎨ n+1 Ψm ≥ 0, (4.4) ⎩ n+1 n+1 Ψm (Um − G(xm )) = 0. n+1 The first step is to compute the intermediate approximation U˜m by solving the equations n (4.2) with known auxiliary term Ψm . Now the second step of the operator splitting method is to derive a relationship in (4.3) n+1 between Um and Ψn+1 m . To do this rewrite the equation (4.3) using equation (4.2) together n+1 with the constrains in (4.4) as a problem to find the pair (Um , Ψn+1 m ), such that n+1 ˜ n+1 −Um 3 Um n = Ψn+1 m − Ψm , 2 δt (4.5) n+1 n+1 Ψm (Um − G(xm )) = 0,
with the constrains n+1 Um ≥ G(xm ) and Ψn+1 ≥ 0. m
Now by solving the problems (4.5)-(4.6), we get n+1 ˜m U ) (G(xm ), Ψnm + 32 G(xm )− n+1 δt , Ψn+1 ) = (Um m 2δt n n+1 ˜ (Um − 3 Ψm , 0) 11
n+1 if U˜m − 2δt Ψnm ≤ G(xm ), 3 otherwise.
(4.6)
(4.7)
Thus, one can do the second step by solving discrete equation (4.2) with the updating formula (4.7). The above implicit method with three time levels requires the value on previous two 0 time levels. The pair (Um , Ψ0m ) on zeroth time level can be obtained by using initial condition 0 and assign value Ψm = 0 or one can use the algorithm discussed in [20]. To find the pair 1 1 (Um , Ψ1m ) at first level, the first step is to compute the intermediate value U˜m as follows: 1 0 − Um U˜m 1 0 + λIΔ Um = Ψ0m . − DΔ U˜m δt 1 The second step is to find the pair (Um , Ψ1m ) such that ˜1 1 if U˜m − δtΨ0m ≤ G(xm ), (G(xm ), Ψ0m + G(xmδt)−Um ) 1 1 (Um , Ψm ) = 1 (U˜m − δtΨ0m , 0) otherwise.
(4.8)
(4.9)
Algorithm to evaluate an American option for n = 0 : N − 1 if n ≤ 1 n+1 n ˜m U −Um n+1 n = DΔ U˜m + λIΔ Um + Ψnm δt n+1 n+1 Um = max(U˜m − δtΨnm , G(xm )) n+1 ˜ n+1 −Um = Ψnm + Um δt Ψn+1 m else
1 3 ˜ n+1 − 2U n + 1 U n−1 = DΔ U˜ n+1 + λIΔ (EU n ) + Ψn U 2 n+1 Um Ψn+1 m δt
end end
5
m
m
2
m
m
m
m
n+1 = max(U˜m − 2δt Ψnm , G(xm )) 3 n+1 ˜ n+1 −Um = Ψnm + 32 Um δt
Numerical simulation and discussion
In this section, we present several numerical experiments to illustrate the efficiency and accuracy of proposed method. We used the FFT algorithm to compute the product of a dense matrix and a column vector in the discrete integral operator in (3.5). Thus the implicit-explicit method with three time levels requires totally O(M N log2 (M )) operations, where N is the number of time steps and M is the number of spatial steps. Numerical results are presented for both call option and put option under Merton and Kou model using multiquadric (with υ = 12 and hence m = 1) radial basis functions. By keeping the shape parameter fixed, the computational error produced by the numerical schemes was measured against the value of analytical solution or reference solution at specific asset price. In all numerical experiments equispaced grid with ni = 3 is used. In all numerical experiments initial two step are calculated with IMEX-BDF1 method. We used MATLAB R2010 on a PC with 2.93GHz Intel(R) Core(TM)i3 CPU and 4.00GB RAM. Example 1. European option under Merton model with parameters σ = 0.15, r = 0.05, σJ = 0.45, μJ = −0.90, λ = 0.10, T = 0.25, K = 100. 12
Table 2: Numerical results for European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 9.283266 9.284893 9.285282 9.285384 9.285409 9.285416
= 90 Error 2.1520e-03 5.2475e-04 1.3655e-04 3.4508e-05 8.7377e-06 2.5525e-06
S = 100 Value Error 3.120808 2.8218e-02 3.142100 6.9262e-03 3.147302 1.7240e-03 3.148595 4.3054e-04 3.148918 1.0776e-04 3.148998 2.7723e-05
S = 110 Value Error 1.396156 5.0297e-03 1.399881 1.3053e-03 1.400852 3.3346e-04 1.401101 8.5126e-05 1.401165 2.1274e-05 1.401180 5.5654e-06
These parameters are also used by YongHoon Kwon, and Younhee Lee [19] for pricing European put option with truncated domain xmin = −1.5 and xmax = 1.5, and by Y. d’Halluin, P. A. Forsythy, and K. R. Vetzal [9] for European call option. The analytical solution is obtained by using infinite series solution of Merton [23]. The reference values for European put option are 9.285418 at S = 90, 3.149026 at S = 100, and 1.401186 at S = 110 and for European call option are 0.527638 at S = 90, 4.391246 at S = 100, and 12.643406 at S = 110. The numerical results for put option and call option are listed in Tables 2 and 3 respectively Table 3: Numerical results for European call options under the Merton model at different asset price with parameters as given in Example 1. S = 90 M N Value Error 129 25 0.525486 2.1520e-03 257 50 0.527113 5.2475e-04 513 100 0.527501 1.3655e-04 1025 200 0.527604 3.4508e-05 2049 400 0.527629 8.7377e-06 4097 800 0.527635 2.5525e-06
S = 100 Value Error 4.363028 2.8218e-02 4.384320 6.9262e-03 4.389522 1.7240e-03 4.390815 4.3054e-04 4.391138 1.0776e-04 4.391218 2.7723e-05
S = 110 Value Error 12.638376 5.0297e-03 12.642101 1.3053e-03 12.643072 3.3346e-04 12.643321 8.5126e-05 12.643385 2.1274e-05 12.643400 5.5654e-06
with shape parameter = 3.8. The option value and absolute error are given at different asset prices. From these tables, we can observe that the discretization scheme is convergent with quadratic rate, since the ratio in respective error is approximately four. The Greeks help to provide important measurements of an option position’s risks. Greeks are the quantities representing the sensitivity of the price of derivatives. Delta measures 13
Table 4: Numerical result for Delta of European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S = 90 Value Error -0.847817 1.1019e-03 -0.847023 3.0807e-04 -0.846794 7.8244e-05 -0.846735 1.9622e-05 -0.846720 4.9210e-06 -0.846717 1.3125e-06
S = 100 Value Error -0.358125 2.4621e-03 -0.356279 6.1556e-04 -0.355817 1.5392e-04 -0.355702 3.8482e-05 -0.355673 9.6151e-06 -0.355665 2.3771e-06
S = 110 Value Error -0.059047 9.4544e-04 -0.058334 2.3319e-04 -0.058158 5.6766e-05 -0.058115 1.3689e-05 -0.058105 3.4368e-06 -0.058102 8.1211e-07
the rate of change of option value with respect to changes in the underlying asset’s price, where as Gamma represents the rate of change in the delta with respect to changes in the underlying price. For European put option the reference values of Delta are −0.84671538 at S = 90, −0.35566306 at S = 100, and −0.05810123 at S = 110, and the values of Gamma are 0.03486014 at S = 90, 0.04882567 at S = 100, and 0.01212941 at S = 110. Since radial basis functions are sufficiently smooth it helps us to compute some Greeks like Delta and Gamma.It has been shown by Sanyasiraju et al.[26] and Wright et al.[32] that RBF-FD approximate first and second derivative of function with second order accuracy. Using the same parameter values which is employed for results in Table 2, we computed two Greeks viz Delta and Gamma at different asset prices and results obtained are reported in Table 4 and 5. From these tables we can observe that convergence rate of the proposed scheme is quadratic for computing Greeks as well. In Fig. 1 we have plotted the option values, the Delta and the Gamma. The plot shows that the option values and Greeks are very stable and no spurious oscillation occur at or around the strike price. These figures show that Greeks are efficiently evaluated using proposed method. Example 2. European call option under Merton model with parameters σ = 0.20, r = 0, σJ = 0.50, μJ = 0, λ = 0.10, T = 1.0, K = 1. The above parameters are also considered by Bhuruth et al.[25] with truncated domain xmin = −4.0 and xmax = 4.0. In their article, Bhuruth et al. provided a new radial basis function based differential quadrature approach (RBF-DQ). They first used RBF-DQ approach for the spatial discretization followed by exponential time integration scheme for time discretization. Using Merton’s closed form solution, the reference value 0.09413550 at strike price is calculated and the option value and absolute error at strike price are reported in Table 6 with shape parameter = 2.2. We observe that the results obtained with the present method have nice agreement with those given in[5, 25]. Example 3. European option under Kou model with parameters σ = 0.15, r = 0.05, η1 = 3.0465, η2 = 3.0775, λ = 0.10, p = 0.3445, T = 0.25, K = 100.. 14
70
60
Option Value
50
40
30
20
10
0
40
60
80
100
120
140
160
Asset Price
(a)
0
0.06
−0.1 0.05
−0.2
Gamma Value
Delta Value
−0.3 −0.4 −0.5 −0.6
0.04
0.03
0.02
−0.7 −0.8
0.01
−0.9 −1
40
60
80
100
120
140
0
160
Asset Price
40
60
80
100
120
140
160
Asset Price
(b)
(c)
Figure 1: European put option value, Delta and Gamma under the Merton model with parameters as provided in the example 1. These parameters are adopted from [19, 9]. To perform the numerical simulations, we choose the truncated domain as xmin = −1.5 and xmax = 1.5. The analytical solution of the model can be found in Kou[18]. The value of European call option at different asset price are evaluated with six digit after the decimal points by d’Halluin et al.[9]. The reference values for call option is 0.672677 at S = 90, 3.973479 at S = 100, and 11.794583 at S = 110. The reference price for put option is computed using put call parity. The numerical experiment 15
Table 5: Numerical result for Gamma of European put options under the Merton model at different asset price with parameters as given in Example 1.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 0.034448 0.034758 0.034835 0.034854 0.034859 0.034860
= 90 Error 4.1186e-04 1.0225e-04 2.5208e-05 6.2791e-06 1.5713e-06 4.1016e-07
S = 100 Value Error 0.049816 9.9002e-04 0.049067 2.4124e-04 0.048886 5.9975e-05 0.048841 1.4976e-05 0.048829 3.7404e-06 0.048827 9.2079e-07
S = 110 Value Error 0.012071 5.8559e-05 0.012117 1.2733e-05 0.012126 3.3105e-06 0.012129 9.0800e-07 0.012129 2.2477e-07 0.012129 6.9091e-08
Table 6: Comparison of different methods for European call options under the Merton model at strike price with parameters as given in Example 2.
M N 65 5 129 10 257 20 513 40
CM[5] Value Error 0.09102 4.01e-03 0.09320 9.32e-04 0.09413 2.72e-04 0.09408 5.39e-05
RBF-DQ[25] Value Error 0.09182 2.32e-03 0.93565 5.70e-04 0.09399 1.38e-04 0.09410 3.02e-05
Present Value 0.092339 0.093865 0.094076 0.094121
Method Error 1.7963e-03 2.7096e-04 5.9464e-05 1.4369e-05
Table 7: Numerical results for European call options under the Kou model at different asset price with parameters as given in Example 3. S = 90 M N Value Error 129 25 0.670284 2.3932e-03 257 50 0.672105 5.7208e-04 513 100 0.672529 1.4783e-04 1025 200 0.672639 3.7620e-05 2049 400 0.672673 3.9368e-06 4097 800 0.672676 1.0867e-06
S = 100 Value Error 3.938910 3.4569e-02 3.964965 8.5136e-03 3.971358 2.1212e-03 3.972948 5.3113e-04 3.973363 1.1569e-04 3.973450 2.9465e-05
S = 110 Value Error 11.787842 6.7414e-03 11.792817 1.7660e-03 11.794129 4.5366e-04 11.794466 1.1692e-04 11.794561 2.2399e-05 11.794577 5.9222e-06
is carried out with shape parameter = 1.7. From Tables 7 and 8, we can observe that the point wise error at different asset price decrease, with convergence rate two. Finally in Fig. 16
70
60
Option Value
50
40
30
20
10
0
40
60
80
100
120
140
160
Asset Price
(a)
0
0.06
−0.1 0.05
−0.2
Gamma Value
Delta Value
−0.3 −0.4 −0.5 −0.6
0.04
0.03
0.02
−0.7 −0.8
0.01
−0.9 −1
40
60
80
100
120
140
0
160
Asset Price
40
60
80
100
120
140
160
Asset Price
(b)
(c)
Figure 2: European put option value, Delta and Gamma under the Kou model with parameters as provided in the example 3. 2, we have plotted the option value and Greeks for put option. These figures shows that Greeks are efficiently evaluated under Kou model as well using proposed method. Example 4. American put option under the Merton model with parameters σ = 0.15, r = 0.05, σJ = 0.45, μJ = −0.90, λ = 0.10, T = 0.25, K = 100. These parameters are consider by d’Halluin, Forsythy, and Labahn [8] and by YongHoon Kwon, and Younhee Lee[20] for Pricing American put option with truncated domain xmin = 17
Table 8: Numerical results for European put options under the Kou model at different asset price with parameters as given in Example 3.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S Value 9.428064 9.429885 9.430309 9.430419 9.430453 9.430456
= 90 Error 2.3932e-03 5.7203e-04 1.4778e-04 3.7571e-05 3.8874e-06 1.0373e-06
S = 100 Value Error 2.696690 3.4569e-02 2.722745 8.5135e-03 2.729138 2.1212e-03 2.730728 5.3108e-04 2.731143 1.1564e-04 2.731230 2.9415e-05
S = 110 Value Error 0.545622 6.7413e-03 0.550597 1.7660e-03 0.551909 4.5361e-04 0.552246 1.1687e-04 0.552341 2.2349e-05 0.552357 5.8728e-06
−1.5 and xmax = 1.5. In [8] authors have developed an implicit discretization coupled with penalty method for pricing such American options. The authors in [20] developed the three time level IMEX method coupled with operator splitting method to solve such problems. The option value and error at different asset price are reported in Table 9 with shape parameter
= 3.8. To obtain error at different asset price in Table 9, we have computed the successive changes in the option value at each asset price. From the above table, we observe that Table 9: Numerical results for American put options under the Merton model at different asset price with parameters as given in Example 4.
M N 129 25 257 50 513 100 1025 200 2049 400 4097 800
S = 90 Value Error 9.998080 10.000259 2.1789e-03 10.003020 2.7610e-03 10.003934 9.1430e-04 10.003835 9.9612e-05 10.003866 3.1405e-05
S = 100 Value Error 3.207266 3.232387 2.5120e-02 3.238920 6.5333e-03 3.240634 1.7132e-03 3.241085 4.5131e-04 3.241207 1.2209e-05
S = 110 Value Error 1.414138 1.418174 4.0361e-03 1.419327 1.1528e-03 1.419660 3.3244e-04 1.419759 9.9174e-05 1.419790 3.0951e-05
point wise errors generally converge with approximately second order accuracy, while near the optimal exercise boundary convergence is less accurate. The above phenomenon is also discussed in [8, 20]. In Fig. 3, we provide the plot of both option values and the early exercise boundary for American put option under the Merton model. Example 5. American put option under the Kou model with parameters σ = 0.15, r = 0.05, η1 = 3.0465, η2 = 3.0775, λ = 0.10, p = 0.3445, T = 0.25, K = 100.. These parameter are also adopted by YongHoon Kwon, and Younhee Lee [20] with same truncated domain and shape parameter. The reference values at different asset prices is 18
70
100
American Payoff 60
98
Holding region
96
Stock Price
Option Value
50
40
30
94
92
20 90
10
0
88
40
60
80
100
120
140
160
Exercise region 0
0.05
0.1
0.15
0.2
0.25
Time (years)
Asset Price
(a)
(b)
Figure 3: American put option value, and optimal early exercise boundary under the Merton model with parameters as provided in the example 4. given by Toivanen [29] are 10.005071 at S = 90, 2.807879 at S = 100, and 0.561876 at S = 110. In Table 10, we get the same scenario as in the case of Merton model. We are Table 10: Numerical results for American put options under the Kou model at different asset price with parameters as given in Example 5. S = 90 M N Value Error 129 25 10.002131 257 50 10.000735 1.3951e-03 513 100 10.005015 4.2794e-03 1025 200 10.005055 3.9945e-05 2049 400 10.005088 3.3276e-05 4097 800 10.005098 9.4650e-06
S = 100 Value Error 2.777966 2.799956 2.1990e-02 2.805766 5.8098e-03 2.807307 1.5414e-03 2.807723 4.1541e-04 2.807834 1.1140e-04
S = 110 Value Error 0.558356 0.560726 2.3706e-03 0.561502 7.7560e-04 0.561750 2.4808e-04 0.561834 8.3698e-05 0.561861 2.7201e-05
getting approximately second order convergence rate, while we have less accuracy near by the early exercise boundary. Finally in Fig. 4, we have plotted the option value and early exercise boundary for American put option under Kou model. Since we are using backward difference method of order two for time semi discretization, and problem have non smooth initial condition, we observe that the increasing the value of ni 19
70
100
American Payoff 60
98
Holding region
96
Stock Price
Option Value
50
40
30
94
92
20 90
10
0
88
40
60
80
100
120
140
160
Exercise region 0
0.05
0.1
0.15
0.2
0.25
Time (years)
Asset Price
(a)
(b)
Figure 4: American put option value, and optimal early exercise boundary under the Kou model with parameters as provided in the example 5. does not seems comparably efficient. Finally cpu times(in seconds) of the method with Merton model for European put option and American put option (example 1 and 4) and with Kou model for for European put option and American put option (example 3 and 5) are given in Table 11. Table 11: CPU times of the proposed methods for different experiment setup.
M
N
129 257 513 1025 2049 4097
25 50 100 200 400 800
Merton model European American option option 0.2260 0.2766 0.0359 0.0617 0.0582 0.0878 0.1602 0.2369 0.5644 0.7946 2.0307 2.7810
20
Kou model European American option option 0.2192 0.2765 0.0512 0.0635 0.0648 0.0913 0.1769 0.2492 0.6078 0.8194 2.1449 2.9971
6
Conclusion
In this article, a RBF based method for solving partial integro-differential equations, which arises when underlying asset follows the jump diffusion process has been presented. After time semi discretization of governing equation, the resulting linear differential equation is solved using RBF based finite difference method, the integral term involved in PIDE is approximated with composite trapezoidal rule. The numerical simulation with European put/call and American put option under Merton and Kou model is carried out. The results are compared against the FD, FEM, and RBF-DQ method available in literature and shown that the proposed method is more accurate and efficient. The proposed method approximate not only the option value but also some of its “Greeks” namely “Delta” and “Gamma” very efficiently.
Acknowledgements The authors would like to thanks the anonymous referees for their valuable comments and suggestion which helped to improve the presentation of the work. This research work of authors Alpesh Kumar, Lok Pati Tripathi is supported by Council of Scientific and Industrial Research, India with grant no. 10-2(5)/2007(ii)-E.U.II and 09/092(0715)/2009-EMR-I respectively.
References [1] A. Almendral, C. W. Oosterlee, Numerical valuation of options with jumps in the underlying, Applied Numerical Mathematics 53 (1) (2005) 1–18. [2] L. Andersen, J. Andreasen, Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing, Review of Derivatives Research 4 (3) (2000) 231–262. [3] M. Briani, R. Natalini, G. Russo, Implicit–explicit numerical schemes for jump–diffusion processes, Calcolo 44 (1) (2007) 33–57. [4] R. Brummelhuis, R. T. Chan, A Radial Basis Function scheme for option pricing in exponential l´evy models, Applied Mathematical Finance 21 (3) (2014) 238–269. [5] P. Carr, A. Mayo, On the numerical evaluation of option prices in jump diffusion processes, European Journal of Finance 13 (4) (2007) 353–372. [6] R. T. L. Chan, S. Hubbert, Options pricing under the one-dimensional jump-diffusion model using the radial basis function interpolation scheme, Review of Derivatives Research 17 (2) (2014) 161–189. [7] G. Chandhini, Y. Sanyasiraju, Local RBF-FD solutions for steady convection–diffusion problems, International Journal for Numerical Methods in Engineering 72 (3) (2007) 352–378.
21
[8] Y. dHalluin, P. A. Forsyth, G. Labahn, A penalty method for american options with jump diffusion processes, Numerische Mathematik 97 (2) (2004) 321–352. [9] Y. dHalluin, P. A. Forsyth, K. R. Vetzal, Robust numerical methods for contingent claims under jump diffusion processes, IMA Journal of Numerical Analysis 25 (1) (2005) 87–112. [10] G. E. Fasshauer, A. Q. M. Khaliq, D. A. Voss, Using meshfree approximation for multiasset American options, Journal of the Chinese Institute of Engineers 27 (4) (2004) 563–571. [11] B. Fornberg, G. Wright, E. Larsson, Some observations regarding interpolants in the limit of flat radial basis functions, Computers & mathematics with applications 47 (1) (2004) 37–55. [12] A. Golbabai, D. Ahmadian, M. Milev, Radial basis functions with application to finance: American put option under jump diffusion, Mathematical and Computer Modelling 55 (3) (2012) 1354–1362. [13] Y. C. Hon, A quasi-radial basis functions method for American options pricing, Computers & Mathematics with Applications 43 (3) (2002) 513–524. [14] Y. C. Hon, X. Z. Mao, A radial basis function method for solving options pricing models, Journal of Financial Engineering 8 (1999) 31–50. [15] S. Ikonen, J. Toivanen, Operator splitting methods for American option pricing, Applied Mathematics Letters 17 (7) (2004) 809–814. [16] S. Ikonen, J. Toivanen, Efficient numerical methods for pricing American options under stochastic volatility, Numerical Methods for Partial Differential Equations 24 (1) (2008) 104–126. [17] M. K. Kadalbajoo, L. P. Tripathi, A. Kumar, Second order accurate IMEX methods for option pricing under merton and kou jump-diffusion models, Journal of Scientific Computing (2015) 1–46. [18] S. G. Kou, A jump-diffusion model for option pricing, Management science 48 (8) (2002) 1086–1101. [19] Y. Kwon, Y. Lee, A second-order finite difference method for option pricing under jump-diffusion models, SIAM Journal on Numerical Analysis 49 (6) (2011) 2598–2617. [20] Y. Kwon, Y. Lee, A second-order tridiagonal method for American options under jumpdiffusion models, SIAM Journal on Scientific Computing 33 (4) (2011) 1860–1872. [21] C. K. Lee, X. Liu, S. C. Fan, Local multiquadric approximation for solving boundary value problems, Computational Mechanics 30 (5-6) (2003) 396–409. 22
[22] X. Liu, K. Tai, Point interpolation collocation method for the solution of partial differential equations, Engineering analysis with boundary elements 30 (7) (2006) 598–609. [23] R. C. Merton, Option pricing when underlying stock returns are discontinuous, Journal of financial economics 3 (1) (1976) 125–144. [24] E. Ngounda, K. C. Patidar, E. Pindza, Contour integral method for European options with jumps, Communications in Nonlinear Science and Numerical Simulation. [25] A. Saib, D. Y. Tangman, M. Bhuruth, A new radial basis functions method for pricing American options under merton’s jump-diffusion model, International Journal of Computer Mathematics 89 (9) (2012) 1164–1185. [26] Y. Sanyasiraju, G. Chandhini, Local radial basis function based gridfree scheme for unsteady incompressible viscous flows, Journal of Computational Physics 227 (20) (2008) 8922–8948. [27] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible navier–stokes equations, Computer Methods in Applied Mechanics and Engineering 192 (7) (2003) 941–954. [28] D. Tangman, A. Gopaul, M. Bhuruth, Exponential time integration and chebychev discretisation schemes for fast pricing of options, Applied Numerical Mathematics 58 (9) (2008) 1309–1319. [29] J. Toivanen, Numerical valuation of European and American options under kou’s jumpdiffusion model, SIAM Journal on Scientific Computing 30 (4) (2008) 1949–1970. [30] A. I. Tolstykh, On using RBF-based differencing formulas for unstructured and mixed structured-unstructured grid calculations, in: Proceedings of the 16th IMACS World Congress, Lausanne, 2000. [31] A. I. Tolstykh, D. A. Shirobokov, On using radial basis functions in a “finite difference mode” with applications to elasticity problems, Computational Mechanics 33 (1) (2003) 68–79. [32] G. B. Wright, B. Fornberg, Scattered node compact finite difference-type formulas generated from radial basis functions, Journal of Computational Physics 212 (1) (2006) 99–123. [33] K. Zhang, S. Wang, Pricing options under jump diffusion processes with fitted finite volume method, Applied Mathematics and Computation 201 (1) (2008) 398–413.
23