BARYCENTRIC SPECTRAL DOMAIN DECOMPOSITION METHODS FOR VALUING A CLASS OF INFINITE ACTIVITY L´EVY MODELS

. A new barycentric spectral domain decomposition methods algorithm for solving partial integro-diﬀerential models is described. The method is applied to European and butterﬂy call option pricing problems under a class of inﬁnite activity L´evy models . It is based on the barycentric spectral domain decomposition methods which allows the implementation of the boundary conditions in an eﬃcient way. After the approximation of the spatial derivatives, we obtained the semi-discrete equations. The computation of these equations is performed by using the barycentric spectral domain decomposition method. This is achieved with the implementation of an exponential time integration scheme. Several numerical tests for the pricing of European and butterﬂy op- tions are given to illustrate the eﬃciency and accuracy of this new algorithm. We also show that Greek options, such as Delta and Gamma sensitivity, are computed with no spurious oscillation.


1.
Introduction. Under jump models, option pricing problems can be modelled by stochastic processes. Such problems were initially introduced in financial institutions in the late 1960s. The famous stochastic model for the equilibrium condition between the expected return, the option, the stock and the risk-less interest rate is the celebrated Black-Scholes equation. The equation was derived by Black and Scholes in 1973 [5]. However, it is well known that the constant volatility Black-Scholes model is not consistent with market prices. Therefore, more general models for stochastic dynamics of the risky assets have been developed. We can mention the deterministic local volatility functions [8,15], the stochastic volatility [22,24], jump-diffusion [28,33], Lévy [30], and infinite Lévy models [16,41]. Note that infinite activity Lévy models incorporate jumps whose intensity is not a finite measure.
Problems that arise in the form of jump models, as in [28], and Lévy models, as in [30], can be modelled by means of Partial-Integro Differential Equations (PI-DEs). Due to inherent complexity in the modelling equations, one can rarely find closed-form analytical solutions to these models; therefore one has to resort to numerical methods. Briani et al. [6] used the numerical schemes to approximate the viscosity solutions for integro-differential, possibly degenerate, parabolic problems, although the approach required very restrictive conditions for stability. Cont and Volchkova [10] used implicit schemes to treat the differential part. The use of the Crank-Nicolson time stepping for the Partial Differential Equation (PDE) portion and explicit evaluation of the convolution integral term were explicitly tested by Tavella and Randall [47]. However, such an asymmetric treatment of the PDE and the integral part introduces biases in the viscosity solution. Hence, the second order convergence is not always achieved. d'Halluin et al. [14] employed the Crank-Nicolson scheme with the Rannacher time smoothing to the PIDE. They used a fixed point iterative procedure as the system solver and obtained second order convergence. Tangman et al. [46] proposed an improved fourth order spectral discontinuity inclusion method. To get around the non-smooth initial condition, they clustered Chebyshev grid points at the discontinuous points and at boundaries. Pindza et al. [40] suggested that to overcome the problem of discontinuity and differentiability in the payoff at strike prices, the use of a grid refinement is one of the best tools. The approach retains a satisfactory accuracy of the spectral method to numerically solve option pricing problems. Ngounda et al. [35] employed the numerical inverse Laplace transform using the Bromwich contour integral approach to provide fast and accurate results for pricing European options with jumps.
In recent years, pure jump Lévy processes of infinite activity that govern stock market returns have been empirically and theoretically studied [12,19,20,38]. These models are flexible enough to perfectly fit the market financial data and capture the excess kurtosis and skewness arising from the risk-neutral distribution returns. The work done in [7] generalizes the VG model to the Carr, German, Madan and Yor (CGMY) model. The model has a jump component, which follows dynamics that can represent either the finite or the infinite activity of either finite or infinite variation. As seen in [7], a diffusion component is not needed if the infinite activity jump process has finite variation. In the literature, finite difference schemes are mostly used for solving infinite activity Lévy PIDE problems. Wang et al. [52] proposed implicit finite difference (FD) methods for the numerical solution of the CGMY model. Almendral et al. [1] used FD methods to discretise the equation in space by the collocation method and using explicit difference backward schemes focused on the case of infinite activity and finite variation. Very recently, Fakharany et al. [16] developed an efficient finite difference scheme for partial integro-differential models related to European and American option pricing problems under a wide class of infinity Lévy models. This paper proposes a study of a barycentric spectral domain decomposition method algorithm for solving PIDE models related to European and butterfly option pricing problems under a class of infinite activity Lévy models. The method is based on barycentric spectral domain decomposition methods [36], which allow the implementation of the boundary conditions in an efficient way. The system of semi-discretised of ordinary differential equations, obtained after approximation of the spatial derivatives using barycentric spectral domain decomposition methods are solved, using an exponential time integration (ETI) scheme. Furthermore, the paper provides several numerical tests for solving both the diffusion and jump parts. This shows the superiority of the method over the popular Crank-Nicolson method. Various numerical results for the pricing of European and butterfly options are also given to illustrate the efficiency and accuracy of this new algorithm. We show that the Greek options such as Delta and Gamma sensitivity measures, are efficiently computed to high accuracy.
The paper unfolds as follows. In Section 2, we provide a review of Lévy processes in finance and various numerical approaches to price options using Lévy processes. The section also presents some mathematical details of the fundamental approach we use to price European options on infinite activity Lévy processes. In Section 3, we provide a quick review of domain decomposition methods (another technique such as barycentric spectral methods is reviewed in Appendix 6). In Section 4, we show how the continuous PIDE of Section 2 must be discretised using the domain decomposition methods of Section 3 resulting in a system of ordinary differential equations. Next, we discuss the use of exponential time differencing schemes to solve the resulting system of ordinary differential equations. Section 5 provides numerical results to some representative problems which illustrate the advantages of the new methods as compared to more standard finite difference methods. The conclusions are provided in Section 6.

2.
The jump-diffusion model. We assume an arbitrage-free market model with a single risky asset price process S t on [0, T ], where t is any time from the initial value time to the maturity time T of the asset. The model follows an exponential Lévy model of the form on the filtered probability space (Ω, F, {F t } t∈[0,T ] , P) where the Lévy process X t has dynamics given by The parameter δ represents the skewness of the process. The jump process is represented by J t = Nt k=1 Y k , where N t is a Poisson process with intensity λ > 0 and {Y k } k≥1 are independent observations from a jump size variable Y .
Let consider V (S, t) to be the option value with the underlying asset S t . Under the assumption of no arbitrage, the existence of an equivalent martingale measure Q to the probability measure P must hold (Q ∼ P). Therefore the asset price S t can have the form (1), where X t is now given by Equation (2). The value for a European option with strike price K is a discounted expected payoff of V (S, t), such that, where Ψ(S T ) is the payoff function, r is the interest rate and E Q is the expected value of the martingale measure Q . The value of a contingent claim V (S, t) on the underlying asset S then solves the PIDE given by where the operator L is defined as

EDSON PINDZA, FRANCIS YOUBI, EBEN MARÉ AND MATT DAVISON
The function f (y) is the Lévy density function given in Table 1. The boundary and the initial conditions make the difference between American and European style options as well as put, call and other types of options. For European vanilla call option, the initial and the boundary conditions are given by where K is the strike price. A butterfly spread is a neutral strategy that has a combination of a bull spread and a bear spread. It is a limited profit, limited risk options strategy. There are three strike prices (discontinuities) involved in a butterfly spread and it can be constructed using calls or puts. The initial and boundary conditions of butterfly spread options are expressed as where S is the stock price K 1 , K 2 and K 3 are three distinct strike prices such that dζ + max(0, λ)e −α|y| Table 1. Density functions for Lévy Processes model ), singularities are observed in the kernel of integration. The parameters of CGMY model are set to be C − = C + = C > 0, G > 0, M > 0 and Y ∈ [0, 2) [7]. The parameter C indicates the overall level of activity. The parameters G and M are the measures depicting the skewness of the Lévy density such that G = M yields a symmetric distribution. When choosing G = M this leads to skewed distributions. The parameter Y describes the fine structure of the stochastic process. At Y = −1, the CGMY leads to a special case of Kou's double exponential model [16]. Furthermore, for Y = 0, we obtain the variance Gamma process. In the case of Y ∈ (0, 1), infinite activity models with finite variation are obtained, while in the case of Y ∈ [1, 2], infinite activity models with infinite variation are depicted. Other singular kernel of integration Lévy processes exist. Meanwhile, the hyperbolic and generalized hyperbolic (GH) distribution are used to obtain a better estimation for the stock returns [45]. Here, the functions J ν (·) and Y ν (·) are the Bessel functions of the first and second kind, respectively. The Meixner process was introduced in 1998. The process is used when the environment is changing stochastically over the time showing a reliable valuation for some indices such as the Nikkei 225 [44,31].
3. Domain decomposition. Challenges arise when we want to approximate a function with a jump discontinuity by using high order finite difference methods or spectral methods. In general, the jump locations of a function and its derivatives are not explicitly known and computationally expensive methods have to be used to detect the location of these jumps [21]. Meanwhile, in the case of infinite activity Lévy models option pricing problems, these discontinuities are located at strike prices and singularities of the density functions. Therefore, spectral domain decomposition methods [32] can be used to recover the accuracy at discontinuity points. Following [39,32], the domain D = [a, b] can be broken into M sub-domains where each sub-domain has its own set of basis functions and expansion coefficients The notation u (µ) represents the approximation in the µth domain, and the different sub-domains D µ can touch or overlap each other. For example, to solve a second order non-linear elliptic PDE or system of equations in some domain D ⊂ R d with boundary conditions where N and d respectively denote the elliptic operator and mappings, the matching conditions must be satisfied. Hence, each function u (µ) , defined only on the single sub-domain D µ , must fit together to form a smooth solution of (10) over the full domain D. For infinite resolution, the following conditions at the limit must hold [39]: 1. When two sub-domains, D µ and D ν , touch each other (patching), on their intersection surface the function and its derivative must be smooth, hence 2. When two sub-domains, D µ and D ν , overlap each other (overlapping), the functions u (µ) and u (ν) must be identical in D µ ∩ D ν . Since the solution of a PDE is unique, we must prove that, at a boundary of the overlapping domain, In this paper, we restrict ourself to the patching methods. The calculation of the integral part can be estimated over multi-domains as follows In the next section we discretise the PIDE (4) by means of spectral domain decomposition methods.
4. Discretisation of the PIDE. Let us begin this section by transforming the PIDE (4) into a simpler one. Since the kernel of the integral in (4) presents a singularity at y = 0, a useful technique is to split the real line, for an arbitrary small parameter ε > 0, into two regions Ω 1 = [−ε, ε] and Ω 2 = R\Ω 1 (the complementary set of Ω 1 in the real line). The integral on Ω 1 is replaced by a suitable coefficient in the diffusion term of the differential part of (4) obtained by Taylor expansion of V (Se y , τ ) about S, see [11,9]. This coefficient depending on ε is a convergent integral and takes the form Letting τ = T − t, the resulting approximating PIDE from (4) is given by The approximation ofσ 2 in (14) is evaluated using the Clenshaw-Curtis quadrature and it is given byσ where φ k = cos kπ N , and λ k , k = 0, 1, 2, ...., N , are the Chebyshev-Gauss-Lobatto (CGL) nodes and the Clenshaw-Curtis weights [18,50], respectively. The improper integrals χ(ε) and γ(ε) in (16) are approximated using the shifted Laguerre-Gauss quadrature [17]. Under consideration of the change of variables η = −y −ε for y < 0 and η = y − ε for y > 0, the expressions χ(ε) and γ(ε) have the following forms and where Here η k are the roots of the Laguerre polynomial L N (η) of degree N defined by and the weightsλ k , k = 1, 2, . . . N , are determined as in [17] bȳ 4.1. Discretisation of the PIDE on a single domain. We transform the PIDE (15) into a constant coefficient PIDE using the transformation S = Ke x and u(x, τ ) = V (S, t). One obtains where We define the numerical domain by D = [y M , y m ]. The discretised version of (23) is given byu . D 2 and D 1 are matrices with entries defined by (54) and D 0 is the identity matrix. J and ζ are defined in (28) and (29), respectively.
For the sake of convenience in the numerical treatment we rewrite the integral part of (22) as follows and k = y − x.
We use the Clenshaw-Curtis quadrature rule to compute the first integral over the interval [y m , y M ] to obtain where u = [u 0 , u 1 , . . . , u N ] T and J is a (N + 1) × (N + 1) matrix with entries Let g L (x, τ ) and g R (x, τ ) be the left and right boundary conditions of the PIDE (22). Therefore, the second integral over R \ [y m , y M ] is approximated using the shifted Laguerre-Gauss quadrature [17] under consideration of the change of variables η = −y + y m for y < y m and η = y + y M for y > y M . This leads to where K = N k=0λ k e η kf (η k + y M − x).

4.2.
Discretisation of the PIDE on multi sub-domains. On each sub-domain D ν , the PIDE can be written as where and Next, we discretise the PIDE (22) in the numerical domain D = [x min , x max ] by means of the spectral domain decomposition method described in Section 3. To this end, we divide the domain D into M sub-domains such that The discretised version of (30) is given bẏ where and I (ν) u = ζ(τ ), which incorporates the boundary conditions. Let L ν = A ν + B ν + I, Equation (35) becomeṡ The solution on the whole domain D is given bẏ where (40) Note that when two sub-domains, D ν and D ν+1 , touch each other, we apply the continuation conditions of the form These last equations are reduced to ODEs when the boundary conditions of a call or put option are imposed on each sub-domains. The obtained ODEs are solved by using Exponential time differencing (ETD) schemes.

4.3.
Exponential time differencing schemes. We consider to solve the system of ODEs (39) where L is either a block dense diagonal matrix or a dense matrix depending on the number of domain taken into consideration, using exponential time differencing methods. Exponential time differencing (ETD) schemes are known as an alternative to implicit methods for solving stiff systems of ODEs [42,25,43]. These methods rely on a fast and stable computation of ϕ-functions i.e., functions of the form (e z − 1)/z. The computation of these functions depends significantly on the structure and the range of eigenvalues of the linear operator and the dimensionality of the semi-discretised PDE. Unfortunately, for spectral methods the linear parts have eigenvalues approaching zero, which lead to complications in the computation of the coefficients. Saad [42], and Hochbruck and Lubich [23] introduced Krylov methods to compute ϕ-functions. Kassam and Trefethen [25] used Cauchy integral representation on a circle for a stable computation of ϕ-functions.
Other evaluations of exponential and related ϕ-matrix functions follow the idea of Schmelzer and Trefethen [43]. This method is based on computing optimal rational approximations to the matrix functions on the negative real axis using the Carathéodory-Fejér procedure [51], The system of ODE (42) can be integrated explicitly on the interval [0 T ] to give The following lemma provides the background for the time-stepping procedure for the evaluation of (44).
has the solution The proof of the above lemma can be found in Nielsen and Wright [34]. The computation of the matrix functions ϕ is obtained by means of the Krylov projection algorithm [34].

5.
Numerical results. In this section, we numerically solve the PIDE discretised in Section 4. Two options are used to compare the accuracy of the CGMY, Meixner, and GH models on the financial PIDE. We refer as Example 1, the case of European call options and as Example 2, the case of butterfly call options. Example 1. We consider a European call option with K = 50, T = 0.5, r = 0.05, q = 0, σ = 0.2, ε = 0.1, x min = −3 and x max = 1. The parameters K, T, r, q, σ respectively represent the strike price, the time of maturation, the interest rate, the continuous dividend, and the volatility of the underlying asset price. The parameters for Lévy models used in this example and arbitrary chosen for computation purpose are given in Table 2. Next, we discretise the PIDE (22)  cal domain D = [x min , x max ] by the means of the spectral domain decomposition method described in Section 3 such that where The domain D is divided into four sub-domains. Figure 1 (a) represents the matrix structure of L in Equation (39), using finite difference and a naive spectral methods for spatial discretisation of the PIDE (22). Note that the matrix in Figure 1 (a) applies for both the finite difference and the naive spectral method discretisation. In the case of FD methods, the matrix L is full due to the discretisation of the non-local part (32). Figure 1 (b) and (c) represent the structure of the matrix L in Equation 39, obtained by means of spectral domain decomposition methods where the continuity conditions (41) were applied. These block diagonal matrices are preferred to full matrices in a sense that they reduce the number of unknowns and the computational time in solving the linear system (42). It is important to note that Figure 1 Table 2. their Greeks (∆ = ∂V ∂S and Γ = ∂ 2 V ∂S 2 ) under CGMY, Meixner, and GH Lévy models. The Greeks measure the sensitivity of the option value with respect to the variations in the asset price and the parameters associated with the model [47]. In practice, accurate approximations to Greeks are needed for hedging purposes. In order to   Table 2.
show the superiority of the SDDM over the FDM, we perform some numerical experiments on the case where Equation (2) is a geometric Brownian motion leading to a standard Black-Scholes equation. To show the efficiency of the present method we report the maximal norm error ME = max |U Analytical − U N umerical | between analytical solutions (available for the Black-Scholes model) and numerical solutions.
In Figure 3, we investigate the trade-off between computational time (CPU time) and the accuracy as the number of grid points (N T ) increases for European call options. Clearly the SDDM is faster and more accurate than the FDM and achieves spectral convergence as expected.
For the general PIDE (5), we report the accuracy of our numerical scheme by means of absolute error AE = |U Benchmark − U N umerical | where U Benchmark and U N umerical represent the benchmark solution computed with N = 150 (the number of grid points in each sub-domain) and the numerical solution, respectively. Table 3 shows the benchmark solution values at S = {40, 50, 60} for different Lévy models. We vary the number of grid points N and compute the absolute errors (AE) for each Lévy model. Table 4 represents all the computed AE values for each Lévy model with different values of N and S. In all the cases, the number of grid points is chosen to be N = 100.
We observe a very rapid decrease of the AE as the number of grid points N increases. Note that the approximations of order 10 −4 , 10 −5 , 10 −7 and 10 −10 in Table 4 are in general difficult to attain with standard finite difference, finite element, and finite volume methods.  Table 4. Absolute errors (AE) of the benchmark and the European call option apply to the CGMY, Meixner and GH processes models with different values of N and S for the parameters in Table  2.
Example 2. In this subsection, we investigate the performance of our proposed method for valuing European butterfly options under Lévy models at the strike prices K 1 = 40, K 2 = 50 and K 3 = 60 using the parameters presented in Table  2 and in Example 1. In this particular case, we need to divide the domain at five different points, namely three different strike prices (K 1 , K 2 and K 3 ), and at two singularities (−ε and ε) present in the kernel of the integral (26)). Figure 4 represents numerical solutions of European butterfly call options and their Greeks (∆ = ∂V ∂S and Γ = ∂ 2 V ∂S 2 ) under the CGMY, Meixner, and GH Lévy models.   Table 2.   Table 6. Absolute errors (AE) of the benchmark and the European call option apply to the CGMY, Meixner, and GH processes models with different values of N and S for the parameters in Table  2.
In Figure 5 we investigate the trade-off between computational time and the accuracy as the number of grid points increases for European vanilla butterfly call options. Clearly the SDDM is faster and more accurate than the FDM and achieves spectral convergence as expected. Table 5 shows the benchmark prices of butterfly call options. Table 6 depicts the AE between benchmark prices and numerical solutions of each model with different values of N and S. We observe a very rapid convergence in the case of European butterfly spread option, which has five regions of singularity. Our approach allows a high resolution of grids around the strike prices K 1 , K 2 and K 3 , and at two singularities −ε and ε present in the kernel of the integral (26). Once again, we obtain approximations of order 10 −4 , 10 −5 , 10 −7 , and 10 −10 in Table 6. 6. Conclusion. We presented a spectral domain decomposition method coupled with the exponential time integrator (44) for pricing European call and European butterfly call options for a class of infinite activity Lévy models, including the CGMY, Meixner, and GH models. Our method produced fast and very accurate results for Example 1 and 2 as compared to FDM. Furthermore, the computed Greeks were free of spurious oscillations. Currently, we are investigating our approach to solve multi-asset Lévy models.
In this section, we review the concept of interpolation, differentiation matrix, and quadrature rule in a barycentric spectral method framework in one domain and multi-domains.
A.1. Spectral barycentric interpolation. The review done by [3] on the Lagrange interpolation and the barycentric formula shows the importance of the the type of spatial discretisation adapted for the performance of spectral methods. At first, a polynomial u N (x) is considered to be found among the vector space of all polynomials of degree N such that u N (x j ) = u j with j = 0, ....., N . The result can be written in Lagrange form as ( [29]) with Lagrange polynomial φ j corresponding to the node x j having the property The evaluation of (48) requires an order O(N 2 ) additions and multiplications. In addition, the presence of instability in the numerical computation is certain. Therefore, modifications of (48) are required to overcomes those disadvantages. Hence, Berrut and Trefethen [4] proposed a more stable barycentric formula of (48) that allows the computation of u N (x) in O(N ) operations. The new formula is written as where w 0 , w 1 , . . . , w N are called barycentric weights. It is well known that an adequate choice of nodes contributes in a decisive manner to improve the accuracy of the approximation and reduce the computational effort. In this paper, a valuable selection is the well-known Chebyshev-Gauss-Lobatto (CGL) nodes x k = cos( kπ N ), k = 0, 1, 2, ...., N , [18,50], which are not uniformly distributed but clustered near the endpoints. The corresponding set of barycentric weights is w 0 = c/2, w k = (−1) k c, k = 1, . . . , N − 1, and w N = (−1) N c/2 for some non-zero constant c [2]. Note that for every set of points {x k }, there is a unique set of barycentric weights {w k }. More details are given in [4] to obtain (50).
A.2. Differentiation matrix. An important task in pseudo-spectral methods is to compute the derivatives u (p) (x) in terms of the values of u(x) at the collocation points x k . A common, practical and efficient way is the use of differentiation matrix which allows one to perform the numerical differentiation in a straightforward way in terms of matrix vector products. For u(x) sufficiently smooth and p a positive integer number, one can approximate the p th derivative of u(x) by u (p) where the entries of the p-order differentiation matrix D (p) are given by

EDSON PINDZA, FRANCIS YOUBI, EBEN MARÉ AND MATT DAVISON
Welfert [53] proposed the following hybrid formula to compute the entries of the p-th order differentiation matrix, In the above, the entries of the first and the second differentiation matrices D (1) and D (2) are given as ( [53]): where j, k = 0, 1, . . . , N . These differentiation matrices are more stable, with respect to rounding errors, than direct evaluation [53]. In view of the above setting, the rational collocation interpolation (50) enjoys exponential convergence property as stated in the following theorem.
The constant ρ refers to the Bernstein ellipse E ρ in which u is analytically continuable. For large regions in the complex plane in which u is analytical, the interpolation polynomials u N has a faster convergence to u as N grows (and the same holds for the respective derivatives).
A.3. Spectral barycentric quadrature. Let u be a real function and N + 1 the number of sampled points. The approximation or linear quadrature rule I ≈ N j=0 λ j u j of the integral I = b a u(x)dx can be achieved by two different ways: • non-uniform grids, e.g., Gauss or Clenshaw-Curtis quadrature; • uniform or equispaced points, e.g., the Newton-Cotes, trapezoidal rule, Simpson, or Boole.
In this article, we consider the Clenshaw-Curtis quadrature. We now show how the replacement of polynomial by linear rational interpolation leads to quadrature formulas which allow arbitrarily large numbers of equispaced nodes. Clearly, every linear interpolation formula trivially yields a linear quadrature rule. For a barycentric rational interpolant, we have: where is the integral of the j th Lagrange fundamental rational function. There are two different ways to compute (57): • One can use direct rational quadrature. The technique consists of applying existing quadrature rules such as Gauss-Legendre and Clenshaw-Curtis [13,49], which are known to perfectly approximate the integrals in (57). • On the other hand, we can apply the indirect rational quadrature. This will produce the integral I = b a u(x)dx through the solution of an ordinary differential equation, see for example [26]. The Clenshaw-Curtis quadrature formula has the following convergence property.
Theorem A.2 ( [13,49]). Let f be analytic function in [−1, 1] and analytically continuable with |f (z)| < M in the closed ellipse E ρ . The error in I N (f ), the Clenshaw-Curtis quadrature of degree N to I(f ) will decay geometrically with the bound In other words, Proof. See [13]. The main advantage with the Clenshaw-Curtis quadrature rule is that its weights and nodes can be computed efficiently via a fast Fourier transform (FFT) in only (O (N ln N )) operations.