On the pricing of multi-asset options under jump-diffusion processes using meshfree moving least-squares approximation

On the Pricing of Multi-Asset Options under Jump-Diffusion Processes using Meshfree Moving Least-Squares Approximation

Journal Pre-proof

On the Pricing of Multi-Asset Options under Jump-Diffusion Processes using Meshfree Moving Least-Squares Approximation Mohammad Shirzadi, Mehdi Dehghan, Ali Foroush Bastani PII: DOI: Reference:

S1007-5704(19)30479-4 https://doi.org/10.1016/j.cnsns.2019.105160 CNSNS 105160

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

24 April 2019 24 October 2019 27 December 2019

Please cite this article as: Mohammad Shirzadi, Mehdi Dehghan, Ali Foroush Bastani, On the Pricing of Multi-Asset Options under Jump-Diffusion Processes using Meshfree Moving LeastSquares Approximation, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.105160

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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. © 2019 Published by Elsevier B.V.

Highlights • Applying the moving least squares approximation for pricing multi-asset European under jump-diffusion models.

• Developing the proposed method for pricing multi-asset option written on L´evy-driven assets with an operator splitting approach.

• Computing hedge parameters of European and American options and early exercise boundary of American options with little extra cost via the proposed method.

• Comparison the proposed method with finite difference method for validating our result.

1

On the Pricing of Multi-Asset Options under Jump-Diffusion Processes using Meshfree Moving Least-Squares Approximation Mohammad Shirzadia , Mehdi Dehghana,∗, Ali Foroush Bastanib,c a

Department of Applied Mathematics, Faculty of Mathematics and Computer Sciences, Amirkabir University of Technology, No. 424 Hafez Ave, Tehran, Iran b Department of Mathematics, Institute for Advanced Studies in Basic Sciences, P.O. Box 45195-1159, Zanjan, Iran c Research Center for Basic Sciences and Modern Technologies (RBST), Institute for Advanced Studies in Basic Sciences (IASBS), P.O. Box 45137-66731, Zanjan, Iran

Abstract The moving least-squares (MLS) approximation is a powerful numerical scheme widely used in the meshfree literature to construct local multivariate polynomial basis functions for expanding the solution of a given differential or integral equation. For partial integro-differential equations arising from the valuation of multi-asset options written on correlated L´evy-driven assets, we propose here an MLS-based collocation scheme in conjunction with implicit-explicit (IMEX) temporal discretization to numerically solve the problem. We apply the method to price both European and American options as well as computing the option hedge parameters. In the case of American options, we use an operator splitting approach to solve the linear complementarity formulation of the problem. Our numerical experiments show the efficiency of the proposed scheme in comparison with some competing approaches, specially finite difference methods. Keywords: Multi-asset option pricing; Moving least-squares method; Jump-diffusion models; Partial integro-differential equations; Linear complementary problem; Meshfree methods. 2010 MSC: 45K05 , 60G51, 60H30, 65M06, 65M70, 35R35.

1. Introduction Recent years have witnessed a tremendous growth and usage of exchange-traded and overthe-counter derivative securities designed to meet various financial needs of market participants like hedging, risk management, speculation and leverage. Among the innovative and tailor-made ∗

Corresponding author. Tel: +982164542503 E-mail addresses: [email protected](M. Dehghan); [email protected](M. Dehghan), [email protected] (M. Shirzadi), [email protected] (A. Foroush Bastani). Preprint submitted to Elsevier

December 30, 2019

products in these markets, the class of multi-asset contingent claims written on several underlying sources of uncertainty has become fairly common investment tools with a variety of pay-off structures and exercise rules [8, 73]. These include but are not limited to options on the maximum and minimum of several assets (max-options and min-options respectively), options on the difference between two prices and an exercise value (spread options), options to exchange one asset for another (out-performance option) and options on an average of prices (basket options). Based on the fact that closed-form analytic valuation formulas are available only in very limited cases under restrictive assumptions on asset price dynamics (e.g. multi-dimensional geometric Brownian motions) or exercise style of the option (e.g. European), the need to develop efficient computational pricing schemes, especially when the underlying assets follow more complex dynamics such as exponential L´evy processes (or more general semi-martingales) has emerged recently (see [13] and the many references therein). The early multi-asset options literature based on the classical Black-Scholes model was pioneered by Margrabe [52] who found an analytic pricing formula for an European-style option to exchange one risky asset for another. This trend was continued by Johnson [35] and Stulz [67] who derived valuation expressions for European options on the maximum and minimum of two assets (see also Johnson [36] for similar results in the multi-asset case). The pricing problem under two or more correlated Brownian diffusion processes is more complicated and must be treated e.g. by numerical partial differential equation (PDE) approaches [59, 74, 75, 53], Monte-Carlo simulations [26] or lattice-based (multinomial) tree methods [7, 20]. These and other numerical approximations are clearly limited to a few underlying assets as the memory requirements and the computational burden increase very rapidly as more underlying assets are considered. Also see [18] for alternative method which is based on a semi–analytical or quasi–numerical approach. In the jump-diffusion framework of Merton [57] and Kou [44], the literature is mostly devoted to single-asset call and put options with a very sparse amount of research work undertaken for multiasset case [12, 39, 51]. Specifically, when the underlying assets follow correlated finite activity jumpdiffusion processes, Clift and Forsyth [12] proposed a finite-difference (FD) approach combining fixed-point iteration and the fast Fourier transform (FFT) for pricing options written on two underlying assets. Kadalbajoo et al. [39] developed a finite element method (FEM) alongside an implicit-explicit (IMEX) time discretization for multi-asset American options under jump–diffusion dynamics. In a recent contribution, Li et al. [51] and Jinghui [34] designed a FEM-based scheme for pricing two-asset European options (with polynomial payoff functions) under L´evy dynamics. For other approximate pricing formulas in the case of multi-asset correlated jump-diffusion processes with idiosyncratic and systematic jump terms, we refer the reader to Xu and Zheng [72] and the references therein. 3

In order to discretize high-dimentional parabolic partial integro-differential equations (PIDEs) arising from multi-asset options, a natural step is to go beyond the conventional mesh-based FD and FE methods and examine the efficiency of mesh-free numerical approaches which were introduced into the numerical analysis field in 70’s and 80’s (see e.g. Fasshauer [23] for a comprehensive review in the field). A hallmark of mesh-free discretization methods is their ability to deal with high-dimensional problems by exploiting a scattered point cloud (without any specified connection structure between the nodes) distributed in the problem domain and its boundary to approximate the derivative and integral terms appearing in the considered equation. Among the popular meshfree methods developed and widely used in the area of scientific computing, the class of radial basis function (RBF) collocation methods [17] and moving leastsquares (MLS) approximation schemes [58, 63] have gained prominence as being fast and accurate PDE solvers available in the literature. These and other mesh-free discretization techniques have been applied recently with great success in single-asset and multi-asset European and American option pricing problems on a variety of asset return distributions (see e.g. [4, 9, 24, 25, 28, 37, 38, 41, 42, 62, 65] for RBF collocation and [43, 70] for MLS approximation. The origins of MLS approximation could be traced back to some ideas in spatial statistics (e.g. optimal least-squares prediction or kriging [54]). It also appeared in the non-parametric statistics literature as a robust tool for local weighted regression and smoothing of scatter-plots [11] (see also the original work of Shepard [66] in the multivariate interpolation context). The general form of the method in its current usage is due to Lancaster and Salkauskas [49] while a through analysis of the interpolation, smoothing and derivative approximation properties of the method are due to Levin [50]. In this paper, we propose an MLS-based approximation scheme followed by direct collocation at interior and boundary nodes to price a collection of multi-asset European and American options under jump-diffusion processes. In this respect, we first show that the proposed method approximates the arbitrage-free price of European options efficiently and accurately with an easily implementable algorithm. As a byproduct, we show that the proposed scheme is capable to approximate the option Greeks (e.g. Delta and Gamma) alongside the price with little extra cost. In the second part of the paper, we present an efficient pricing engine for American multi-asset contingent claims based on operator-splitting approach which solves the linear complementarity sub-problems at each time step by an efficient iterative procedure (see [32, 68] for more details). This will enable us to find the location of the early exercise boundary at each time level easily and efficiently. Similar to the European case, we are able to compute the option hedge parameters conveniently. Comparison with a FEM based approach presented in Li et al. [51] and also two FD based 4

discretization schemes (one based on the work of Kwon and Lee [45, 46] and another one presented in the Appendix) shows that the proposed method gives better performance in terms of error and ease of implementation both in single-asset and multi-asset European and American cases. The reminder of this paper is organized as follows. In Section 2, we derive the PIDE for multi-asset European option pricing problem under jump–diffusion processes and present the corresponding localization error estimates and defining boundary conditions for the resulting problem. Section 3 contains the mathematical foundations of the MLS approximation and Section 4 presents the details of our discretization scheme composed of an IMEX time marching combined with MLS collocation in space for valuation of European option contracts. Extension of the framework presented in Section 4 to price American-style options based on the incorporation of an operator splitting idea is considered in 5. The details of computing the hedge parameters are investigated in Section 6 while Section 7 is devoted to extensive numerical experiments on a variety of single-asset and multi-asset options to demonstrate the effectiveness of the proposed methodology. We conclude the paper by presenting some final remarks on the domain of applicability of the presented MLS based scheme to other problems in computational finance. 2. Multi-Dimensional Exponential L´ evy Dynamics In the remainder, we consider a basket of d ≥ 1 risky assets, S = (S1 , ..., Sd ), whose log-returns

are modeled by a L´evy process with state-space Rd (see e.g. Reich et al. [61] for more details).

Let (Ω, F, (Ft )t≥0 , Q) be a filtered probability space in which Ω is the set of all possible events,

(Ft )t≥0 is the filtration produced by the sigma-algebra of price vector S(t) = (S1 (t), ..., Sd (t)) and Q is a multi-dimensional risk-neutral probability measure. Let us also assume that the (risk-neutral) dynamics of the price vector, S under Q is given by Si (t) = Si (0) exp(rt + Xi (t)),

i = 1, 2, · · · , d,

(2.1)

in which S(0) = (S1 (0), S2 (0), · · · , Sd (0)) denotes the vector of initial asset prices, r is the risk-free

interest rate and X(t) = (X1 (t), ..., Xd (t)) is a d-dimensional time-homogeneous jump-diffusion

(L´evy) process starting from zero with the characteristic triplet (Σ, Γ, ν). Here Σ is a symmetric positive semi-definite covariance matrix, Γ denotes the expected return vector of the underlying assets and ν is a (L´evy) measure on Rd , satisfying Z ν({0}) = 0 and min(1, |z|2 )ν(dz) < ∞.

(2.2)

the process X(t) could be expressed as # " Z h i  h T i  1 exp(iθT z) − 1 − iθT z1{|z|≤1} ν(dz) . E eiθ X(t) = exp t − θT Σθ + iθT Γ + 2 Rd

(2.3)

Rd

It is well known (based on the L´evy-Khinchine theorem [13]) that the characteristic exponent of

5

In order to obtain a detailed stochastic representation for the risk-neutral dynamics of the i-th asset, let us assume that the (L´evy) measure, ν, satisfies the extra condition Z ezi ν(dz) < ∞.

(2.4)

|z|>1

Now, each component, eXi , will be a martingale with respect to the filtration (F)t≥0 , if and only if Σi,i + Γi + 2

Z

Rd



 ezi − 1 − zi 1|z|≤1 ν(dz) = 0.

(2.5)

So, the risk-neutral dynamics of each underlying asset, Si (t) is given by Si (t) = Si (0) +

Z

0

t

rSi (u−)du +

d Z X j=1

t

Aij Si (u−)dWj (u) +

Z tZ 0

0

R

(ezi − 1)Si (u−)JXi (du × dzi ), (2.6)

in which A is the Choleski factor of Σ (i.e. Σ = AAT ), W (t) = (W1 (t), · · · , Wd (t)) is a ddimensional Brownian motion process with correlation matrix Σ and JXi (·) is the Poisson random measure describing the jumps of Xi . Remark 2.1. In the sequel, we assume that the jump components of the asset price processes are statistically independent, in the sense that ν(dz) = ν1 (dz1 ) × · · · × νd (dzd ).

(2.7)

This assumption implies that the asset prices never jump together and its imposition will ease the presentation of the material by reducing multi-dimensional integrals to a sum of univariate integrals. 2.1. Deriving the PIDE Describing the Price of Multi-Asset Options Based on the classical theory of risk-neutral valuation, the price of an European multi-asset option at time t with maturity T , strike price K and pay-off function VT (S) is given by   V (t, S) = EQ e−r(T −t) VT (S) S(t) = S .

(2.8)

    If we now assume that V (t, S) ∈ C1,2 (0, T ) × Rd+ ∩ C0 [0, T ] × (Rd+ ∪ {0}) , then V (t, S) will

solve the backward Kolmogorov equation of the form

6

 d d X  1 X ∂2V ∂V ∂V   +r − rV + Σi,j Si Sj Si    ∂t 2 ∂Si ∂Sj ∂Si  i,j=1 i=1  d Z h X ∂V i zi zi  + νi (dzi ) = 0, V (t, S , · · · , S e , · · · , S ) − V (t, S , · · · , S ) − (e − 1)S  1 i 1 i d d   ∂Si  i=1 R    V (T, S1 , · · · , Sd ) = VT (S1 , · · · , Sd ),

(2.9)

defined on the spatio-temporal domain [0, T ] × Rd+ (see e.g. [61]). In the above PIDE, νi ’s,

i = 1, 2, · · · , d are absolutely continuous L´evy measure components defined on R with densities λi fi (z) where λi is the jump intensity and fi (z) depends on the model being used. The two most common choices are the normal distribution  (z − µi )2  1 J fi (z) = i √ exp − , 2 σJ 2π 2σJi

(2.10)

with given µJ and σJ in the Merton’s model [57] and the double-exponential distribution fi (z) = pλi+ exp(−λi+ z)1z≥0 + (1 − pi )λi− exp(−λi− z)1z<0 ,

(2.11)

with λi− > 0, λi+ ≥ 1 and 0 ≤ pi ≤ 1 in the Kou’s model [44].

The pay-off function, VT (S1 , · · · , Sd ) in equation (2.9) takes different forms based on the option

type, among them we mention the following: • Basket call/put option

VT (S1 , · · · , Sd ) =

 d X +    w S − K , (call on basket)  i i  i=1

d  +  X    K − w S , i i 

(2.12)

(put on basket)

i=1

in which

x+

= max{x, 0} and wi denotes the weight of i-th asset.

• Max call/put option

  +  max(S1 , .., Sd ) − K , (call on max)   VT (S1 , · · · , Sd ) =  K − max(S , .., S ) + , (put on max) 1 d

(2.13)

• Min call/put option

  +  min(S1 , .., Sd ) − K , (call on min)   VT (S1 , · · · , Sd ) =  K − min(S1 , .., Sd ) , (put on min) 7

(2.14)

• Spread option on two assets ( VT (S1 , S2 ) = (S1 − S2 − K)+ , (call on spread) VT (S1 , S2 ) = VT (S1 , S2 ) = (K − S1 − S2 )+ , (put on spread)

(2.15)

In order to simplify the formulation of the PIDE (2.9), we apply change of variables (see e.g. [45]) of the form   τ = T − t, ,      x = ln( Si ), i Si (0)  u(τ, x) = exp(−rτ )V (t, S(0) exp(x)),      h(x) = V (S(0) exp(x)), T

(2.16)

in which x = (x1 , · · · , xd )T . Equation (2.9) could now be re-written as uτ = L[u] := D[u] + I[u],

(τ, x) ∈ [0, T ] × Rd ,

(2.17)

where D[u] := I[u] =

d X

d

i,j=1

d Z X i=1

X ∂u ∂2u + αi − Λu, Σi,j ∂xi ∂xj ∂xi

R

i=1

(2.18)

u(τ, x1 , · · · , xi + y, · · · , xd )νi (dy),

R with ξi := R (ey − 1)fi (y)dy, α = (α1 , ..., αd )T = (r − 21 σ12 − λ1 ξ1 , ..., r − 21 σd2 − λd ξd )T and P Λ = di=1 λi . The pricing of a multi-asset European option could now be expressed as the solution of the following multi-dimensional Cauchy problem uτ = D[u] + I[u] on (τ, x) ∈ (0, T ] × Rd ,

u(0, x) = h(x), x ∈ Rd .

(2.19)

2.2. Localization and Boundary Conditions In order to numerically solve PIDE (2.19), we first need to localize the problem domain, Rd into a bounded computational region, ΩM := {x ∈ Rd : kxk ≤ M } for some positive M and some

norm k · k on Rd . For making the problem well-defined on the new bounded domain, we need to

specify suitable boundary conditions on the boundary of ΩM by defining a (continuous) function g(τ, x), x ∈ ∂ΩM for τ ∈ (0, T ]. Among the possible choices for the function g, the two most

important ones are g(τ, x) ≡ 0 (extending the solution by zero outside ΩM ) and g(τ, x) = h(x)

(extending the solution by pay-off value outside ΩM ). As shown in Cont and Voltchkova [14], the

second choice will result in asymptotically close approximations to the exact solution at infinity. They also showed that these two choices will cause the total truncation error, eM = u − uM (uM 8

being the solution of the localized problem) to exponentially decrease with increasing the size of ΩM . There exists also a third proposal to choose g(τ, x) by considering it as a linear combination of exponential functions with coefficients determined by substituting g(τ, x) into the PIDE (2.19). In the special case of an European basket put option, g(τ, x) is given by g(τ, x) := max(K − S0

d X

wi exp(xi + rτ ), 0),

i=1

(τ, x) ∈ ∂ΩM .

(2.20)

As shown in [51, 34], this choice will result in a total truncation error exponentially decreasing with a rate depending on M as given in Theorem 2.2 below. We first recall that L2 (Ω) denotes the space of square integrable functions on Ω and the Sobolev spaces W k,p (Ω) are defined by W k,p (Ω) = {u ∈ Lp (Ω) : Dα u ∈ Lp (Ω) ∀|α| ≤ k}, equipped with the norm   1 p  P p α , 1 ≤ p < ∞, |α|≤k kD ukLp (Ω) kukW k,p (Ω) =  max|α|≤k kDα ukL∞ (Ω) , p = ∞.

(2.21)

(2.22)

In the above expressions, α = (α1 , · · · , αd )T is a multi-index with non-negative integer components αi , i = 1, · · · , d and Dα is defining by Dα u(x) = in which |α| = [2]).

∂ |α| u(x) ∂xα11 · · · ∂xαdd

Pd

i=1 αi .

(2.23)

We also use the notation H k (Ω) = W k,2 (Ω) in the remainder (see e.g.

Theorem 2.2. (Li et al. [51]) There exist constants C = C(T ) > 0 and α = α(M ) > 0, such that

keM (τ, .)k2L2 (ΩM/2 )

+

Z

0

τ

keM (s, .)k2H 1 (ΩM/2 ) ds ≤ C exp(−αM ).

(2.24)

3. Mesh-Free MLS Approximation Moving least-squares approximation is a universal approach to reconstruct smooth multidimensional surfaces from a finite set of scattered data nodes in which (in contrast to ordinary least-squares) both the weight functions and coefficients are considered to be functions of the space variable. In this section, we present a general overview of the method and some of its features which will be used in the sequel. 9

Consider a predetermined set of scattered data nodes, X = {x1 , · · · , xn } ⊆ Ω ⊂ Rd (in a

neighborhood of x) with assigned values, u = {u(x1 ), · · · , u(xn )} coming from a scalar function u : Ω → R. For a quasi-interplant of the form Pu (x) =

n X

x ∈ Rd ,

u(xi )ψi (x),

i=1

(3.1)

the goal of MLS approximation is to find the expansion coefficients ψi , i = 1, · · · , n in (3.1) such

that Pu becomes the best approximation out of Πdm (the space of d-variate polynomials of total degree at most m) for u at given x. Using the Backus-Gilbert formulation of the MLS scheme [6], our aim is to find the generating functions ψi (·), i = 1, 2, · · · , n pointwise at evaluation point x by minimizing the quadratic functional n

1X 2 ψi (x)wi (x), 2

(3.2)

i=1

subject to the linear constraint n X

p(xi )ψi (x) = p(x),

i=1

∀p ∈ Πdm ,

(3.3)

which is known as polynomial reproduction condition (see e.g. [22] or [71] for more details). In the above expression, wi (·)’s, i = 1, · · · , n, are positive weight functions centered at the point xi which vanish outside a distance δ from its center point.

Using the Lagrange multipliers λ(x) = [λ1 (x), · · · , λQ (x)]T with Q =

ψi (x) could be computed as

m+d d




, each function

Q

ψi (x) =

1 X λj (x)pj (xi ), wi (x) j=1

i = 1, · · · , n,

(3.4)

in which p(x) = [p1 (x), · · · , pQ (x)]T forms a basis for Πdm . The vector λ(x) in (3.4) solves the

linear system

G(x)λ(x) = p(x),

(3.5)

in which G (known as the Gram matrix) is defined componentwise by Gjk (x) = hpj (x), pk (x)iφ(x) :=

n X

pj (xi )pk (xi )φi (x) = PT DP,

i=1

j, k = 1, 2, · · · , Q,

(3.6)

where D = diag[φ1 (x), · · · , φn (x)] ∈ Rn×n with φi (x) = 1/wi (x) and P = [Pij ] = pi (xj ) ∈ Rn×m , i = 1, · · · , Q, j = 1, · · · , n. Using the above notation, the generator function Ψ(x) = [ψ1 (x), · · · , ψn (x)]T could be obtained as

Ψ(x) = DPλ(x) = DPG−1 p(x) = DP(PT DP)−1 p(x). 10

(3.7)

Now, if we define u = [u(x1 ), · · · , u(xn )]T and use the expression (3.1), the MLS approximation could be written as n X Pu (x) = u(xi )ψi (x) = ΨT (x)u,

(3.8)

i=1

and its β-order derivative could be obtained as Dβ Pu (x) =

n X

u(xi )Dβ ψi (x).

(3.9)

i=1

Based on the fact that the MLS generating function ψi (·) for each i depends on all polynomial basis functions, pj (·), j = 1, ..., m and also the weight function, φi (·), we need the following lemma which gives us the necessary conditions for the order of smoothness of the MLS approximation: Lemma 3.1. (Wendland [71]) If the weight functions φi (·) for each i = 1, 2, · · · , Q possess k ∈ N continuous derivatives, then the MLS approximation (3.8) will belong to C k (Ω).

4. MLS Approximation for Multi-Asset European Options The moving least-squares approximation could be employed in the numerical solution of PDEs and PIDEs in different guises. When it is used in the context of a Galerkin scheme for discretizing the weak form of the governing equations, the MLS shape functions will play the role of trial and test functions and this will result in the element-free Galerkin family of schemes (we refer the interested reader to [5] and the references therein). On the other hand, if this approximation scheme is applied to the strong form of the equation in a collocation context, the resulting scheme will be called a finite point method (see e.g. [21]). The error analysis of the finite point method has been studied in e.g. Cheng and Cheng [10] for general elliptic PDEs. In the remainder, we use an implicit-explicit (IMEX) temporal discretization based on finitedifferences in the time direction and a collocation approach based on MLS shape functions in the space variable for pricing European options given by the solution of PIDE   (τ, x) ∈ (0, T ] × ΩM ,  uτ = L[u],  u|τ =0 = h(x), x ∈ ΩM ,    u| ∂ΩM = g(x), τ ∈ (0, T ].

(4.1)

4.1. Temporal Discretization

Consider a partition of [0, T ] into sub-intervals of the form [τn−1 , τn ] for n = 1, 2, · · · , Nt with

0 = τ0 < τ1 < · · · < τNt = T using a fixed time step-size, ∆τ =

define our approximations to the exact solution at τ = τn and τ  1 1 un (x) := u(τn , x), un+ 2 (x) := u(τn , x) + u(τn+1 , x) . 2 11

T Nt . For each n = 1, 2, · · · , Nt , = τn +τ2n+1 respectively by

we

(4.2)

Applying the well-known IMEX temporal discretization (see e.g. [14]) to solve (4.1), we arrive at the expression 1 un+1 (x) − un (x) = D[un+ 2 (x)] + I[un (x)], ∆τ

(4.3)

where re-arranging the terms gives us H+ un+1 (x) = H− un (x),

un = [un (x1 ), · · · , un (xn )].

(4.4)

The operators H+ and H− are defined respectively by ∆τ ∆τ H+ := 1 − D =1− 2 2

d X

i,j=1

! d X ∂2 ∂ + −Λ , Σi,j αi ∂xi ∂xj ∂xi

∆τ ∆τ H− := 1 + D + ∆τ I = 1 + 2 2

d X

Σi,j

i,j=1

i=1 ∂2

∂xi ∂xj

+

d X i=1

! ∂ − Λ + ∆τ I. αi ∂xi

(4.5)

In order to solve the elliptic partial differential equation (4.4), we use a collocation-based finite point method described in the reminder. 4.2. Spatial Approximation For approximating the function un+1 (x) in (4.4) at time level n + 1, we propose to use an MLS expansion of the form n

u (x) =

N X

un (xi )ψi (x),

(4.6)

i=1

in which N = Ni + Nb and X = X1 ∪ X2 is the union of internal and boundary nodes with

X1 = {x1 , ..., xNi } ⊂ Ω and X2 = {xNi +1 , ..., xN } ⊂ ∂Ω. Using a collocation approach (see e.g. [21]) based on the nodes X1 and X2 , we arrive at a system of linear equations (H+ Ψ)un+1 = (H− Ψ)un ,

(4.7)

where [H+ Ψ]ij = H+ ψi (xj ) and [H− Ψ]ij = H− ψi (xj ), i, j = 1, · · · , n. Starting from the initial

condition, u0 = h(x), we can solve (4.7) in a forward recursive manner and the boundary conditions

could be imposed through the iteration at each time step. In order to construct the entries of matrices H− and H+ in (4.7), we use the expansion (2.23)

to evaluate the derivative terms in the equation. For computing I[u] in equation (4.5), we first truncate the entire real line to [−A, A] and then use a trapezoidal integration scheme with step size ∆x := 2A/N to arrive at I[u] ≈

d Z X i=1

A

−A

u(τ, x1 , · · · , xi + y, · · · , xd )νi (dy) ≈ 12

Kr d X X

i=1 j=Kl

νij ui+j ,

(4.8)

in which Kl and Kr are constant parameters for which 1 1 [−A, A] ⊂ [(Kl − )∆x, (Kr + )∆x], 2 2

(4.9)

and ξi ≈

Kr X

j=Kl

(e

yj

− 1)νij ,

νij =

Z

(j+1/2)∆x

νi (dy).

(4.10)

(j−1/2)∆x

In order to have an estimate of the incurred error during the space-discretizatin phase of the proposed numerical scheme, we state and prove Theorem 4.1 below. In its statement, we will need the concept of fill distance, hX,Ω , defined by hX,Ω = sup min kx − xj k2 .

(4.11)

x∈Ω 1≤j≤N

Theorem 4.1. Let un (x) and uh,n (x) be the exact and approximate solutions of the equation (4.7), respectively. Then, we have kun (x) − uh,n (x)k ≤ C Cond(HΦ−1 )hQ−1 X,Ω , h,n ku (x)k

(4.12)

in which C is a constant, H = H−1 + H− , [Φ]ij = ψi (xj ), i, j = 1, · · · , N and the condition number,

Cond(·) is computed using the corresponding matrix norm.

Proof. The result follows from the fact that wi (x) ∈ CQ (ΩM )∩W Q,∞ (ΩM ) and un (x) ∈ CQ+1 (ΩM )∩

H Q+1 (ΩM ) and employing Theorem 3 in [10] with s = 2.



In the next section, we apply the proposed method to price multi-asset American options under multi-dimensional L´evy processes. 5. MLS Approximation for Multi-Asset American Options The MLS-based IMEX scheme discussed in the previous section could be extended to the American option case in a straightforward manner. It is well-known (see e.g. [33]) that the price of an American option satisfies a variational inequality of the form     uτ − L[u] ≥ 0, u − h(x) ≥ 0,    (u − L[u])(u − h(x)) = 0, τ

(5.1)

for (τ, x) ∈ (0, T ] × Ω in which the operator, L and the payoff function, h(·) are as defined in the European case (2.17) and (2.16), respectively. We must note in the passing that this general 13

formulation appears in many other contexts in sciences and engineering (see e.g. [15, 16, 27, 29, 30, 47, 48, 55, 56, 60]). In order to solve the above variational inequality formulation of the problem, different approaches have been proposed in the literature (see e.g. [1, 64] for detailed expositions). If we discretize the partial differential operator by finite differences or finite elements, we arrive at a linear complementarity problem (LCP) in each time step which could now be solved by an iterative scheme like the projected successive overrelaxation (PSOR) method [64]. As an efficient alternative, we present the details of an operator splitting method due to Ikonen and Toivanen [31] in the sequel (see e.g. [3, 32, 69] for more details). 5.1. Operator Splitting Method The basic idea in the operator splitting approach is to define an auxiliary function ϕ = D[u] − I[u] and then rewriting the LCP formulation (5.1) as   ϕ = ∂u − D[u] − I[u], ∂τ  ϕ ≥ 0, u − h ≥ 0, ϕ · (u − h) = 0,

∂u ∂τ

−

(5.2)

in (0, T ] × Ω with the initial condition h(·). We now could split the first equation in (5.2) on the (n + 1)-th time level into the following two intermediate steps o u en+1 − un n h u en+1 + un i − D + I[un ] = Φn , ∆τ 2

un+1

− ∆τ

un

n hu o en+1 + un i − D + I[un ] = Φn+1 . 2

(5.3)

So the solution of (5.2) reduces to finding a pair (un+1 , Φn+1 ) satisfying (5.3) under the conditions

un+1 ≥ h(X),

Φn+1 ≥ 0,

Φn+1 (un+1 − h(X)) = 0.

(5.4)

In order to find the relationship between un+1 and the supplemental term Φn+1 , we rewrite the two expressions in (5.3) under the constraints (5.4) as  n+1 −u en+1  u = Φn+1 − Φn , ∆τ  n+1 n+1 Φ (u − h(X)) = 0,

(5.5)

with the restrictions un+1 ≥ h(X),

Φn+1 ≥ 0.

(5.6) 14

Now the solution of (5.5) in conjunction with the constraints expressed in (5.6) could be obtained as  

h(X) − u en+1 ), if u en+1 − ∆τ Φn ≤ h(X), (un+1 , Φn+1 ) = ∆τ  n+1 (e u − ∆τ Φn , 0), Otherwise. (h(X), Φn +

(5.7)

The final step of the operator splitting approach is to solve the first equation in (5.3) with the update formula (5.7). This implicit method needs the value of the option on the previous time level. The pair, (u0 , Ψ0 ), on the initial level could be obtained using the initial condition and setting Ψ0 = 0 (see e.g. [46]). 5.2. Treatment of the Early Exercise Boundary In order to determine a fair price for an American option contract, the early exercise provision must be taken into account by specifying an optimal exercise boundary, S(τ ), from which to ¯ ) or determine whether or not to exercise the option at time t (by checking whether S(τ ) ≤ S(τ ¯ ), respectively). In terms of the log-price variable, the last expression is equivalent to S(τ ) > S(τ ¯ ), which solves the functional equation finding, S(τ   ¯ )) = K − S0 exp S(τ ¯ ) . u(τ, S(τ

¯ ) at each time level n, by solving the (generally nonlinear) equation It is possible to approximate S(τ

Gn (x) = un (x) − K + x = 0,

(5.8)

with an iterative Newton-like scheme of the form x(m+1) = x(m) −

Gn (x(m) ) , (Gn )0 (x(m) )

(5.9)

in which Gn (x) and its derivative, (Gn )0 (x), are computed as n

G (x) =

n X i=1

n 0

(G ) (x) =

ψi (x)un (xi ) − K + x,

n X

ψi0 (x)un (xi ) + 1,

i=1

where un (xj ), j = 1, ..., n denote the option price value in the log-price variable.

15

6. Computing Option Greeks as Hedge Parameters In the options trading world, the term “Greeks” (including delta, gamma, vega, theta and many others as special cases) refers to quantities providing information on the sensitivity of option’s price to factors such as the underlying asset’s price, volatility, interest rate, time, etc. They are used to hedge the option position’s risks in the case of adverse movements in the market prices and assess its potential rewards against risks. The Greek delta, ∆, measures the sensitivity of an option’s theoretical value to a change in the price of the underlying asset. It is defined by ∆ :=

∂V , ∂S

(6.1)

and is normally represented as a number between minus one and one. It indicates how much the value of an option should change when the price of the underlying asset rises by one dollar. On the other hand, the Geek gamma, Γ, describes the rate of change of the considered option’s delta with respect to the underlying asset’s price and is defined by ∂2V . ∂S 2

Γ :=

(6.2)

In the case that the MLS expansion (3.9) of the price is at hand, we do not require any extra computational cost to obtain the delta and gamma values in the proposed method. With the change of variables (2.16), delta is given by ∆ :=

∂V 1 ∂u = , ∂S S ∂x

(6.3)

and gamma is computed as Γ :=

1  ∂ 2 u ∂u  ∂2V = − . ∂S 2 S 2 ∂x2 ∂x

(6.4)

d n X 1 X ∂ψj uj , Si ∂xi

(6.5)

Now using the MLS approximation (3.9), we can get the approximation ∆=

i=1

j=1

for the delta Greek and Γ=

d n X ∂ψj  1 X  ∂ 2 ψj uj , − 2 2 ∂xi S ∂x i i i=1 j=1

(6.6)

for the gamma Greek. It should be mentioned that similar strategies could be exploited for measuring other Greek quantities. 16

7. Numerical Experiments In this section, we present the results of our numerical simulations which illustrate the accuracy and efficiency of the proposed MLS-based collocation scheme for pricing multi-asset European and American option contracts. In this respect, we consider four test cases which span an acceptable range of parameter values, payoff structures and L´evy measures. For each test problem, we compare our results with those obtained from other approximation schemes, sometimes from the existing literature and sometimes with the newly developed approximations. In all test problems, we consider the space of polynomials of degree at most 2 (i.e. Π12 := span{1, x, x2 } in the univariate case and Π22 := span{1, x, y, x2 , xy, y 2 } in the two-variate case).

We will also use a regular point distribution on the spatial domain with mesh spacing h as our

collocation points with Nx and Ny denoting the number of points in x and y directions, respectively. We employ the C4 Wendland’s compactly supported weight function φ(r) = (1 − r)6+ (35r2 +

18r +3) where the support is selected to be δ = 2h in all of our experiments. The truncated domain will be set to [−1.5, 1.5] in the single-asset case and [−1.5, 1.5] × [−1.5, 1.5] for two-asset options. To measure the accuracy, we consider the following maximum error E = max |uk − u ˜k |,

(7.1)

1≤k≤N

where u ˜ and uk denote the approximate and exact solution at kth node, respectively. Test Problem 1 (Single-Asset American Put Option) As a first instance, we consider an American option under both Merton’s and Kou’s dynamics. The parameters are chosen to be σ = 0.15, r = 0.05, σJ = 0.45, µJ = −0.9, λ = 0.10, T = 0.25, K = 100,

(7.2)

for the Merton’s case and σ = 0.15, r = 0.05, λ− = 3.0775, λ+ = 3.0465, p = 0.3445, λ = 0.10, T = 0.25, K = 100, (7.3) for the Kou’s dynamics which are the same as those used in Kwon and Lee [46].

17

Table 1: The numerical results obtained for American put option prices under Merton’s jump-diffusion model.

S = 90

S = 100

S = 110

Nx

Nt

MLS

FDM [46]

MLS

FDM [46]

MLS

FDM [46]

128

25

10.011116

9.997971

3.199700

3.198242

1.413177

1.410850

256

50

10.000442

10.000196

3.230860

3.231112

1.418567

1.417464

512

100

10.003078

10.003226

3.238671

3.238672

1.419675

1.419167

1024

200

10.003910

10.004099

3.240638

3.240649

1.419804

1.419641

2048

400

10.003850

10.003840

3.241112

3.241126

1.419773

1.419763

2.8316e-06

1.8000e-5

4.2627e-05

1.2500e-04

2.0569e-05

3.9999e-05

Error

In Tables 1 and 2, we have reported the maximum error in option prices under Merton’s and Kou’s model for S = 90, 100 and S = 110, where the reference values are 10.003822 at S = 90, 3.241251 at S = 100, and 1.419803 at S = 110 in the Merton model and 10.005071 at S = 90, 2.807879 at S = 100, and 0.561876 at S = 110 in the Kou model [46]. Table 2: The numerical results obtained for American put option prices under Kou’s jump-diffusion model.

S = 90

S = 100

S = 110

Nx

Nt

MLS

FDM [46]

MLS

FDM [46]

MLS

FDM [46]

128

25

10.01090

10.001336

2.767558

2.767419

0.555976

0.553315

256

50

10.000388

10.000180

2.797691

2.797994

0.560942

0.553315

512

100

10.004940

10.005150

2.805261

2.797994

0.561853

0.561266

1024

200

10.005046

10.005150

2.807196

2.807297

0.561891

0.561710

2048

400

10.005077

10.005129

2.807698

2.807744

0.561830

0.561831

6.0009e-07

5.8000e-05

6.4332e-05

1.3500e-04

8.0323e-05

4.5000e-05

Error

Figure 1 shows the payoff function and option surface at the maturity for the range S ∈ [60, 140]

in the two models. We also have plotted the Delta and Gamma Greeks at the maturity in the same range in Figure 2. Figure 3 shows the early exercise boundary of the considered option under the two models. Test Problem 2 (Two-Asset Polynomial Option) Here, we consider a polynomial option with the payoff function VT (S1 , S2 ) := (S1 + S2 )2 having an exact analytic solution (see e.g. [51]). We use the parameters reported in Table 3 (also used in [51]) and apply our collocation scheme for different correlation structures between the assets. 18

Table 4 contains the computed option prices for various combinations of asset volatilities and correlations. We have used the values Nx = Ny = 31 as the resolution of spatial grid and Nt = 100 for time discretization. Comparison with the analytic solution (in both cases of positive and negative correlations) shows that in almost all cases considered (except the two rows corresponding to (σ1 , σ2 ) = (0.1, 0.2) and (0.2, 0.2)), the presented MLS method yields more accurate values than the FEM-based method of [51]. We also have plotted the option price surface in the polynomial option case in Figure 4. Table 3: Parameters used in test problem (2) for polynomial option.

Parameters

Values

Mean jump sizes (µ1 , µ2 ) Mean jump volatilities

(-0.9,-0.9)

(σJ1 , σJ2 )

(0.45,0.45)

Jump intensities (λ1 , λ2 )

(0.1,0.1)

Underlying prices (S1 , S2 )

(40,40)

Strike price (K)

80

Interest rate (r)

0.05

Maturity time (T )

0.1

Table 4: The numerical results obtained for polynomial option prices.

ρ

0.3

-0.3

σ1

σ1

Analytic Solution [51]

FEM [51]

Relative Error

MLS

Relative Error

0.1

0.1

6.4899e+03

6.5695e+03

0.0122

6.5209e+03

0.0047

0.1

0.2

6.4958e+03

6.4785e+03

0.0026

6.5268e+03

0.0047

0.1

0.3

6.5050e+03

6.4434e+03

0.0094

6.5361e+03

0.0047

0.2

0.2

6.5026e+03

6.4977e+03

0.0007

6.5221e+03

0.0029

0.2

0.3

6.5128e+03

6.4611e+03

0.0079

6.5323e+03

0.0029

0.3

0.3

6.5239e+03

6.4646e+03

0.0090

6.5350e+03

0.0017

0.1

0.1

6.4899e+03

6.5695e+03

0.0122

6.5190e+03

0.0044

0.1

0.2

6.4958e+03

6.4785e+03

0.0026

6.5229e+03

0.0041

0.1

0.3

6.5050e+03

6.4434e+03

0.0094

6.5303e+03

0.0038

0.2

0.2

6.5026e+03

6.4977e+03

0.0007

6.5144e+03

0.0018

0.2

0.3

6.5128e+03

6.4611e+03

0.0079

6.5207e+03

0.0012

0.3

0.3

6.5239e+03

6.4646e+03

0.0090

6.5177e+03

0.0009

Test Problem 3 (Two-Asset European Basket Put Option) As the third test case, we consider the pricing of a two-asset European basket put option with the parameters listed in Table 5. 19

Table 5: Parameters of an European basket put option.

Parameters

Values

Diffusion volatility (σ1 , σ2 )

(0.1, 0.1)

Mean jump size (µ1 , µ2 )

(-0.9,-0.9)

Mean jump volatility

(σJ1 , σJ2 )

(0.45,0.45)

Jump intensity (λ1 , λ2 )

(0.1,0.1)

Correlation (ρ)

0.3

Weights (w1 , w2 )

(0.5,0.5)

Strike price (K)

80

Interest rate (r)

0.05

Maturity time (T )

0.25

In Table 6, we present the results of our computations for a variety of spatial grid spacings (Nx ×Ny )

and also a finite difference scheme proposed in this paper (see the Appendix at the end of the paper). The MLS collocation scheme converges as the grid spacing increases which is in accordance with the theoretical results in Section 4. In Table 7 we also compute order of convergence, defined as log2

kuh k , kuh/2 k

of the proposed method at (S1 , S2 ) = (60, 60). From Theorem 4.1 we know that the

order of the convergence is depend on hX,Ω as well as Cond(HΦ−1 ). Since Cond(HΦ−1 ) depends on hX,Ω , as hX,Ω → 0 the Cond(HΦ−1 ) increases and effect the term hQ−1 X,Ω and we can not say anything about the order of convergence.

Table 6: Result for European two-assets basket put option.

MLS (11 × 11)

MLS (21 × 21)

MLS (31 × 31)

MLS (41 × 41) 18.1763

18.4884

(60,70)

12.7992

13.2185

13.2920

13.3521

13.6118

(60,80)

8.4817

8.5058

8.5420

8.5671

8.7336

(60,90)

5.2510

4.2961

4.2108

4.1397

4.1475

(70,60)

12.8843

13.3302

13.4004

13.4578

13.6118

(80,60)

8.6406

8.6876

8.7162

8.7339

8.7336

(90,60)

5.4512

4.4908

4.3976

4.3173

4.1475

(70,70)

8.1769

8.5909

8.5967

8.6551

8.7328

(80,80)

0.3028

0.8212

1.1044

1.2540

1.2135

(S1 , S2 ) (60,60)

17.7834

18.0366

18.1302

20

FDM

Table 7: The convergence rate of error estimate obtained for Example 3.

h

MLS

Error

Order

0.15

17.7834

0.0380

—

0.075

18.0366

0.0240

0.20

0.375

18.1302

0.0190

0.10

0.0187

18.1763

0.1160

0.07

In order to have a feeling about the behaviour of the error term in our proposed scheme, we have plotted the price surface in conjunction with the incurred error at each point in Figure 5. We have used a FDM-based scheme using Nx × Ny = 101 × 101 points as our benchmark solution. We also

have plotted the Delta and Gamma Greeks for this option in Figure 6.

Test Problem 4 (Two-Asset American Basket Put Option) As the last example in this section, we consider an American two-asset basket put option with parameters listed in Table 5. We report the results obtained from applying our proposed scheme as well as the ones coming from a finite difference scheme in Table 8. We have used three different asset correlations (ρ) and in each case, different volatility combinations (σ1 , σ2 ) are examined. All of the results correspond to (S1 , S2 ) = (80, 80). Again we observe a good agreement between the obtained results and also the benchmark solution. Table 8: The numerical results obtained for American basket put option prices using MLS and FDM schemes.

ρ = 0.3

ρ = 0.0

σ1

σ2

MLS

FDM

MLS

FDM

0.1

0.1

1.3359

1.3321

1.2249

0.1

0.2

2.0180

1.9650

0.1

0.3

2.7271

0.2

0.2

0.2 0.3

ρ = −0.3 MLS

FDM

1.2070

1.1076

1.0782

1.8400

1.7991

1.6466

1.6084

2.6515

2.5440

2.4586

2.3525

2.2465

2.6229

2.5122

2.3596

2.2345

2.0607

1.9078

0.3

3.2809

3.1381

2.9453

2.8000

2.5672

2.4187

0.3

3.8989

3.7125

3.4680

3.2793

2.9912

2.7875

Figure 7 concerns the surface plot of the corresponding American basket put option obtained correspondingly from the MLS and FDM approximations. The Delta and Gamma Greeks with ρ = −0.3 and ρ = 0.3 are also plotted in Figure 8 and Figure 9, respectively. Last but not least, Figure 10 demonstrates the early exercise boundary which splits the (S1 , S2 ) price space into three

different regions corresponding to exercise with respect to S1 , exercise with respect to S2 and the holding region. As it could be seen from this figure, there exists an intersection point of the early 21

exercise boundaries. It has been proved in [40] that for any (S1 , S2 ) below this intersection point, it is optimal to exercise the option. 8. Conclusion In this study, we proposed and analyzed a mesh-free method based on moving least squares collocation to discretize multi-dimensional partial integro-differential equations arising from multiasset option pricing problems under L´evy processes. We employed the method to price European as well as American options and to compute the hedging parameters such as Delta and Gamma. We compared the performance of the MLS collocation scheme with a finite-difference discretization where the results show advantages of the new method in terms of computational costs and ease of implementation. As possible extensions, this method could be employed when the dynamics of the underlying asset follow more complex processes such as infinite activity L´evy or regime-switching jump-diffusion models. Acknowledgements The authors would like to thank the editor and two anonymous referees for providing valuable and helpful suggestions and comments on the paper. The third author also gratefully acknowledges support by the Institute for Advanced Studies in Basic Sciences (IASBS) Research Council under grant No. G2018IASBS22614. Conflict of Interest The authors declare that they have no conflict of interest. References [1] Y. Achdou, O. Pironneau, Computational Methods for Option Pricing, vol. 30, Siam, 2005. [2] K. Atkinson, W. Han, Theoretical Numerical Analysis, vol. 39, Springer, 2005. [3] L. V. Ballestra, L. Cecere, A fast numerical method to price American options under the Bates model, Computers & Mathematics with Applications 72 (5) (2016) 1305–1319. [4] L. V. Ballestra, G. Pacelli, Pricing European and American options with two stochastic factors: A highly efficient radial basis function approach, Journal of Economic Dynamics and Control 37 (6) (2013) 1142–1167. [5] T. Belytschko, Y. Y. Lu, L. Gu, Element-free Galerkin methods, International Journal for Numerical Methods in Engineering 37 (2) (1994) 229–256. [6] L. Bos, K. Salkauskas, Moving least-squares are Backus-Gilbert optimal, Journal of Approximation Theory 59 (3) (1989) 267–275. [7] P. P. Boyle, J. Evnine, S. Gibbs, Numerical evaluation of multivariate contingent claims, The Review of Financial Studies 2 (2) (1989) 241–250.

22

[8] O. Brockhaus, Equity Derivatives and Hybrids, Palgrave Macmillan, 2016. [9] 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. [10] R. Cheng, Y. Cheng, Error estimates for the finite point method, Applied Numerical Mathematics 58 (6) (2008) 884–898. [11] W. S. Cleveland, Robust locally weighted regression and smoothing scatterplots, Journal of the American Statistical Association 74 (368) (1979) 829–836. [12] S. S. Clift, P. A. Forsyth, Numerical solution of two asset jump diffusion models for option valuation, Applied Numerical Mathematics 58 (6) (2008) 743–782. [13] R. Cont, P. Tankov, Financial Modelling with Jump Processes, CRC Press, 2003. [14] R. Cont, E. Voltchkova, A finite difference scheme for option pricing in jump diffusion and exponential L´evy models, SIAM Journal on Numerical Analysis 43 (4) (2005) 1596–1626. [15] S. Dafermos, Sensitivity analysis in variational inequalities, Mathematics of Operations Research 13 (3) (1988) 421–434. [16] D. Damircheli, M. Bhatia, Solution approaches and sensitivity analysis of variational inequalities, in: AIAA Scitech 2019 Forum, 2019, p. 0977. [17] M. Dehghan, A. Shokri, A numerical method for solution of the two–dimensional sine–Gordon equation using the radial basis functions, Mathematics and Computers in Simulation 79 (3) (2008) 700–715. [18] M. Dehghan, S. Pourghanbar, Solution of the Black–Scholes equation for pricing of barrier option, Z. Naturforsch., A, 66a(5), (2011) 289–296. [19] D. J. Duffy, Finite Difference Methods in Financial Engineering: a Partial Differential Equation approach, John Wiley & Sons, 2013. [20] N. Ekvall, A lattice approach for pricing of multivariate contingent claims, European Journal of Operational Research 91 (2) (1996) 214–228. [21] G. E. Fasshauer, Approximate moving least-squares approximation for time–dependent PDEs, in: WCCM V, Fifth World Congress on Computational Mechanics (http://wccm. tuwien. ac. at), HA Mang, FG Rammerstorfer, and J. Eberhardsteiner (eds.), Vienna University of Technology, 2002. [22] G. E. Fasshauer, Approximate moving least-squares approximation with compactly supported radial weights, in: Meshfree Methods for Partial Differential Equations, Springer, 2003, pp. 105–116. [23] G. E. Fasshauer, Meshfree Approximation Methods with MATLAB, vol. 6, World Scientific, 2007. [24] G. E. Fasshauer, A. Q. M. Khaliq, D. A. Voss, Using meshfree approximation for multi-asset American options, Journal of the Chinese Institute of Engineers 27 (4) (2004) 563–571. [25] A. Foroush Bastani, Z. Ahmadi, D. Damircheli, A radial basis collocation method for pricing American options under regime-switching jump-diffusion models, Applied Numerical Mathematics 65 (2013) 79–90. [26] P. Glasserman, Monte Carlo methods in financial engineering, vol. 53, Springer Science & Business Media, 2013. [27] C. D. Ha, Application of degree theory in stability of the complementarity problem, Mathematics of Operations Research 12 (2) (1987) 368–376. [28] M. Haghi, R. Mollapourasl, M. Vanmaele, An RBF–FD method for pricing American options under jump– diffusion models, Computers & Mathematics with Applications 76 (10) (2018) 2434–2459. [29] P. T. Harker, J.-S. Pang, Existence of optimal solutions to mathematical programs with equilibrium constraints, Operations Research Letters 7 (2) (1988) 61–64. [30] P. T. Harker, J.-S. Pang, Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications, Mathematical Programming 48 (1-3) (1990) 161–220.

23

[31] S. Ikonen, J. Toivanen, Operator splitting methods for American option pricing, Applied Mathematics Letters 17 (7) (2004) 809–814. [32] 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. [33] P. Jaillet, D. Lamberton, B. Lapeyre, Variational inequalities and the pricing of American options, Acta Applicandae Mathematica 21 (3) (1990) 263–289. [34] Z. Jinghui, Multi-asset Option Pricing with L´evy Process, Ph.D. thesis (2009). [35] H. Johnson, The pricing of complex options, Unpublished manuscript. [36] H. Johnson, Options on the maximum or the minimum of several assets, Journal of Financial and Quantitative Analysis 22 (3) (1987) 277–283. [37] M. K. Kadalbajoo, A. Kumar, L. P. Tripathi, Application of the local radial basis function-based finite difference method for pricing American options, International Journal of Computer Mathematics 92 (8) (2015) 1608–1624. [38] M. K. Kadalbajoo, A. Kumar, L. P. Tripathi, A radial basis function based implicit–explicit method for option pricing under jump-diffusion models, Applied Numerical Mathematics 110 (2016) 159–173. [39] M. K. Kadalbajoo, L. P. Tripathi, A. Kumar, An error analysis of a finite element method with IMEX-time semidiscretizations for some partial integro-differential inequalities arising in the pricing of American options, SIAM Journal on Numerical Analysis 55 (2) (2017) 869–891. [40] J. Kay, M. Davison, H. Rasmussen, The early exercise region for Bermudan options on two underlyings, Mathematical and Computer Modelling 50 (9-10) (2009) 1448–1460. [41] S.-M.-M. Kazemi, M. Dehghan, A. Bastani, On a new family of radial basis functions: Mathematical analysis and applications to option pricing, Journal of Computational and Applied Mathematics 328 (2018) 75–100. [42] L. Khodayari, M. Ranjbar, A numerical method to approximate multi-asset option pricing under exponential L´evy model, Computational Economics 50 (2) (2017) 189–205. [43] Y. Kim, H.-O. Bae, H. K. Koo, Option pricing and Greeks via a moving least square meshfree method, Quantitative Finance 14 (10) (2014) 1753–1764. [44] S. G. Kou, A jump-diffusion model for option pricing, Management Science 48 (8) (2002) 1086–1101. [45] 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. [46] Y. Kwon, Y. Lee, A second-order tridiagonal method for American options under jump-diffusion models, SIAM Journal on Scientific Computing 33 (4) (2011) 1860–1872. [47] J. Kyparisis, Sensitivity analysis framework for variational inequalities, Mathematical Programming 38 (2) (1987) 203–213. [48] J. Kyparisis, Solution differentiability for variational inequalities, Mathematical Programming 48 (1-3) (1990) 285–301. [49] P. Lancaster, K. Salkauskas, Surfaces generated by moving least squares methods, Mathematics of Computation 37 (155) (1981) 141–158. [50] D. Levin, The approximation power of moving least-squares, Mathematics of Computation of the American Mathematical Society 67 (224) (1998) 1517–1531. [51] X. Li, P. Lin, X.-C. Tai, J. Zhou, Pricing two-asset options under exponential l´evy model using a finite element method, arXiv preprint arXiv:1511.04950. [52] W. Margrabe, The value of an option to exchange one asset for another, The Journal of Finance 33 (1) (1978) 177–186. [53] J. Mart´ın-Vaquero, A. Khaliq, B. Kleefeld, Stabilized explicit Runge–Kutta methods for multi-asset American

24

options, Computers & Mathematics with Applications 67 (6) (2014) 1293–1308. [54] G. Matheron, The intrinsic random functions and their applications, Advances in Applied Probability 5 (3) (1973) 439–468. [55] L. McLinden, An analogue of Moreaus proximation theorem, with application to the nonlinear complementarity problem, Pacific Journal of Mathematics 88 (1) (1980) 101–161. [56] N. Megiddo, M. Kojima, On the existence and uniqueness of solutions in nonlinear complementarity theory, Mathematical Programming 12 (1) (1977) 110–130. [57] R. C. Merton, Option pricing when underlying stock returns are discontinuous, Journal of Financial Economics 3 (1-2) (1976) 125–144. [58] D. Mirzaei, R. Schaback, M. Dehghan, On generalized moving least squares and diffuse derivatives, IMA Journal of Numerical Analysis 32 (2012) 983–100. [59] B. F. Nielsen, O. Skavhaug, A. Tveito, Penalty methods for the numerical solution of American multi-asset option problems, Journal of Computational and Applied Mathematics 222 (1) (2008) 3–16. [60] J. Pang, Solution differentiability and continuation of Newton’s method for variational inequality problems over polyhedral sets, Journal of Optimization Theory and Applications 66 (1) (1990) 121–135. [61] N. Reich, C. Schwab, C. Winter, On Kolmogorov equations for anisotropic multivariate L´evy processes, Finance and Stochastics 14 (4) (2010) 527–567. [62] A. Safdari-Vaighani, A. Heryudono, E. Larsson, A radial basis function partition of unity collocation method for convection–diffusion equations arising in financial applications, Journal of Scientific Computing 64 (2) (2015) 341–367. [63] R. Salehi, M. Dehghan, A moving least square reproducing polynomial meshless method, Applied Numerical Mathematics, 69 (2013) 34–58. [64] R. Seydel, R. Seydel, Tools for Computational Finance, vol. 3, Springer, 2006. [65] V. Shcherbakov, E. Larsson, Radial basis function partition of unity methods for pricing vanilla basket options, Computers & Mathematics with Applications 71 (1) (2016) 185–200. [66] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings of the 1968 23rd ACM National Conference, ACM, 1968, pp. 517–524. [67] R. Stulz, Options on the minimum or the maximum of two risky assets: analysis and applications, Journal of Financial Economics 10 (2) (1982) 161–185. [68] J. Toivanen, Numerical valuation of European and American options under Kou’s jump-diffusion model, SIAM Journal on Scientific Computing 30 (4) (2008) 1949–1970. [69] J. Toivanen, et al., Application of operator splitting methods in finance, in: Splitting Methods in Communication, Imaging, Science, and Engineering, Springer, 2016, pp. 541–575. [70] S. Vahdati, D. Mirzaei, The finite points approximation to the PDE problems in multi-asset options, CMES: Computer Modeling in Engineering & Sciences 109 (3) (2015) 247–262. [71] H. Wendland, Scattered Data Approximation, vol. 17, Cambridge University Press, 2004. [72] G. Xu, H. Zheng, Approximate basket options valuation for a jump-diffusion model, Insurance: Mathematics and Economics 45 (2) (2009) 188–194. [73] P. G. Zhang, Exotic Options: A Guide to Second Generation Options, World Scientific, 1997. [74] R. Zhang, Q. Zhang, H. Song, An efficient finite element method for pricing American multi-asset put options, Communications in Nonlinear Science and Numerical Simulation 29 (1) (2015) 25–36. [75] R. Zvan, P. A. Forsyth, K. Vetzal, A finite volume approach for contingent claims valuation, IMA Journal of Numerical Analysis 21 (3) (2001) 703–731.

25

Appendix In this appendix, we present a θ-weighted finite-difference approximation scheme for the numerical solution of PIDE (4.1). For simplicity of presentation, we consider the two-asset case (see Duffy [19] for more general schemes in higher dimensional cases). In this respect, we first consider a mesh grid on the truncated domain (0, T ] × [−M, M ] × [−M, M ] and for positive integers Nt , Nx1 and Nx2 , we define ∆t = T /Nt , ∆x1 = 2M/Nx1 and ∆x2 = 2M/Nx2 . Also, we define x1h = −M + h∆x1

and x2k = −M + k∆x2 for h = 0, 1, · · · , Nx1 and k = 0, 1, · · · , Nx2 . For each tn = n∆t with an

integer n = 0, 1, · · · , Nt , we define unh,k := u(tn , x1h , x2k ) and consider the difference approximations

 ∂u n 1  θ unh−1,k − unh+1,k + (1 − θ) un+1 − un+1 , ' h+1,k h+1,k ∂x1 h,k 2∆x1



 ∂u n 1  n+1 θ unh,k+1 − unh,k−1 + (1 − θ) un+1 − u , ' h,k+1 h,k−1 ∂x2 h,k 2∆x2



1  ∂ 2 u n n+1 n+1 n+1 n n n , + u − 2u + (1 − θ) u + u − 2u θ u ' h−1,k h,k h+1,k h−1,k h,k h+1,k ∂x21 h,k ∆x21



∂ 2 u n 1  n+1 n+1 n+1 n n n , + u − 2u + (1 − θ) u + u − 2u θ u ' h,k−1 h,k h,k+1 h,k−1 h,k h,k+1 ∂x22 h,k ∆x22 
∂ 2 u n 1 θ unh+1,k+1 − unh+1,k−1 − unh−1,k+1 + unh−1,k−1 ' ∂x1 ∂x2 h,k 4∆x1 ∆x2

n+1 n+1 n+1 +(1 − θ) un+1 h+1,k+1 − uh+1,k−1 − uh−1,k+1 + uh−1,k−1




(8.1) ,

for a given constant θ ∈ [0, 1]. If we choose θ = 0, θ = 1 and θ = 1/2 the obtained schemes are called respectively as implicit Euler, Explicit Euler and Crank-Nicolson schemes.

The operator D in equation (2.18) could now be estimated at each point (tn , x1h , x2k ) by e n . Dunh,k ≈ Du h,k

(8.2)

In order to numerically approximate the integral term I in (2.18), we use the same trapezoidal

rule as in (4.8) and (4.10). Finally, using the expression n un+1 ∂u n h,k − uh,k = , ∂t h,k ∆t

(8.3)

the approximate solution unh,k could be obtained at each grid point (tn , x1h , x2k ).

26

Merton

40 35

35

30

30

25

25

20

20

15

15

10

10

5

5

0 60

70

80

90

100

110

120

130

Kou

40

Payoff Option value

0 60

140

Payoff Option value

70

80

90

S

100

110

120

130

140

S

Figure 1: American put option value for single-asset American option under Kou and Merton’s models.

Merton

0 -0.1

-0.1

-0.2

-0.2

-0.3

-0.3

-0.4

-0.4

-0.5

-0.5

-0.6

-0.6

-0.7

-0.7

-0.8

-0.8

-0.9

-0.9

-1 60

70

80

90

100

Kou

0

110

120

130

-1 60

140

70

80

90

S

110

120

130

140

110

120

130

140

S

Merton

0.07

100

Kou

0.06

0.06

0.05

0.05 0.04 0.04 0.03 0.03 0.02 0.02 0.01

0.01

0 60

70

80

90

100

110

120

130

0 60

140

S

70

80

90

100

S

Figure 2: Delta and Gamma values for single-asset American option under Kou and Merton’s models.

27

Merton

100

98

Early exercise boundary

Early exercise boundary

98

96

94

92

90

88

Kou

100

96

94

92

90

0

0.05

0.1

0.15

0.2

88

0.25

T

0

0.05

0.1

0.15

0.2

0.25

T

Figure 3: Optimal early exercise boundary for single-asset American put option under consideration models.

Figure 4: The surface of polynomial option with the parameter of Table 3.

28

Figure 5: European two-assets basket option values in conjunction with error estimate.

Figure 6: European two-assets Delta and Gamma of basket option.

29

Figure 7: American two-asset basket option surfaces with MLS and FDM methods.

Figure 8: American two-asset Delta of basket option with correlation parameters ρ = −0.3, 0.3.

30

Figure 9: American two-asset Gamma of basket option with correlation parameters ρ = −0.3, 0.3.

120

110

Hold option 100

S2

Exercise region K-S1 90

80

Exercise region K-S2

70 Intersection point

60 60

70

80

90

100

110

120

S1

Figure 10: Early exercise boundary for two-asset American put option.

31