A time multidomain spectral method for valuing affine stochastic volatility and jump diffusion models

A time multidomain spectral method for valuing affine stochastic volatility and jump diffusion models

Journal Pre-proof

A time multidomain spectral method for valuing affine stochastic volatility and jump diffusion models Claude Rodrigue Bambe Moutsinga, Edson Pindza, Eben Mare´ PII: DOI: Reference:

S1007-5704(19)30478-2 https://doi.org/10.1016/j.cnsns.2019.105159 CNSNS 105159

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

11 May 2019 8 November 2019 27 December 2019

´ A Please cite this article as: Claude Rodrigue Bambe Moutsinga, Edson Pindza, Eben Mare, time multidomain spectral method for valuing affine stochastic volatility and jump diffusion models, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.105159

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 • A time spectral domain decomposition method is proposed for valuing affine stochastic volatility and jump diffusion models. • The proposed method is elegant and accurate. • The method is robust as it handles time domain decomposition, plus exponential convergence is offered.

1

A time multidomain spectral method for valuing affine stochastic volatility and jump diffusion models Claude Rodrigue Bambe Moutsinga∗ , Edson Pindza† , Eben Mar´e‡ ∗ † ‡ †

Department of Mathematics and Applied Mathematics, University of Pretoria Pretoria 002, South Africa Achieversklub School of Cryptocurrency and Entrepreneurship, 1 Sturdee Avenue, Rosebank 2196, South Africa

Abstract The general form of existing multivariate models in finance are in based on stochastic volatility and jump diffusion models. These models typically lead to the need of solving systems of stiff Riccati differential equations. In this paper, we propose a time spectral domain decomposition method for solving systems of stiff Riccati differential equations. The technique is applied to solving stiff diffusion model problems found in oil pricing, interest rate and electricity models. Numerical methods show that the present approach is efficient, highly accurate and a good alternative to the existing numerical methods. Keywords: Stochastic volatility, Affine jump diffusion, Correlated state process, Poisson process, Time-spectral method, Domain decomposition. 2010 MSC: 337K10, 44A15, 45K05, 65M12, 65M70.

1. Introduction A quest for efficient solutions for option pricing problems in the financial market has been an ongoing topic of research among academics and practitioners. In the early 70s, Fisher Black and Myron Scholes published an influential paper where they presented the BlackScholes model, which soon became a major source for most financial traders for pricing and hedging options. The model represents a parabolic partial differential equation (PDE) with constant coefficients to be solved, see [3, 7]. In this context, closed-form solutions could still be found in the case of one-dimensional problems with few parameters. However, empirical studies show that volatility of financial assets follows a stochastic process and also that underlying stock prices experience jumps. In this case, option pricing problems result in solving a multi-dimensional partial integro-differential equation (PIDE), which becomes more complicated to solve as the dimension increases. Various numerical solutions have been proposed for solving these problems. In one dimension, classical numerical methods such as the Finite Difference Method (FDM) or Finite Element Method (FEM) [12, 11] were used to solve these PIDEs, nu∗

Correspondence to: Claude Rodrigue Bambe Moutsinga, Email: [email protected]

2

merically. Pindza et al. [23, 24, 25] proposed a robust spectral method to obtain accurate solutions to these PIDEs. Ngounda et al. [21] provided a combined spectral domain decomposition method and an inverse Laplace transform with a Bromwich contour integral to obtain spectral convergence for the pricing of European options. In two dimensions, in’t Hout et al [28] proposed new ADI schemes to solve the pricing of the European option under the Bates model numerically; Zhu et al [34] developed a Legendre quadrilateral spectral element approximation for the Black-Scholes equation to price European options with two underlying assets. Their methodology displayed an exponential convergence. Clift et al [27] derived a modified finite differential method by combining a fixed point iteration scheme with an FFT to avoid dense linear system solutions and improve the efficiency of the existing methods. In three dimensions, computational intensive works were conducted by Taruvinga et al [29], they used an accurate and efficient method of line to compute a three-asset American option with jumps. The above methods are without a doubt very effective in the sense that they provide accurate solutions of the Black-Scholes PIDE directly in a reasonable time. However, as the dimension gets higher, these numerical methods begin to suffer from the curse of dimensionality. A tremendous amount of computation can be reduced if one could avoid the PIDE by taking advantage of the dynamics of the underlying asset and considering volatility models without departing so much from reality, unlike the original Black-Scholes model that considered volatility to be constant. So, in order to maintain a balance between the effectiveness of solution and the solvability of option pricing problems, Heston [13] presents a model with stochastic volatility that includes correlation and the capability to provide a semi-closed form solution for option prices in terms of the characteristic function of the Log-stock price. The solution is more realistic, consistent with the smile and the skewness observed from the empirical data. The model is later on generalised by Duffie et al. [10] to cater for an entire class of affine jump diffusion processes. The main advantage with such processes lies in applying some transform analysis to arrive at a system of stiff Riccati equations that can still be solved analytically although costly, but also numerically. Many approaches have been used to obtain numerical solutions for the stiff Riccati equations, these approaches include finite difference methods and homotopy perturbation methods [20]. Unfortunately these methods tend to require many points to achieve good results and also turn out to lose their accuracy for a larger time scale. This is due to the fact that some components in the solution decay more rapidly than others, thus forcing any numerical method, with a finite region of stability, to use excessively more step sizes in order to get smoothness of solution. This inevitably leads to a decrease of the efficiency and accumulates more machine roundoff error [14]. Huang [15], Dahlquist and Bjorck [8] recommend the use of implicit methods when it comes to stiff problems. For these reasons we introduce a time spectral domain decomposition method for solving a Riccati system of equations, numerically. The literature around spectral methods is rich. These methods have been used in several areas such as computational fluid dynamics Canuto et al. [6], Hussaini and Zang [16], modeling different types of wave motion Komatitsch and Tromp [17], weather forecasting Bourke et al. [4], and finance Pindza et al. [23, 24]. Spectral methods have the particularity 3

of handling technical boundary value problems in the sense that the solution is approximated by a series of polynomials, the particularity of the polynomials reside in the way coefficients are calculated. In the case of the Chebyshev spectral method, the solution is given as a Chebyshev series and the differential equation is shifted to the interval [−1, 1]. Chebyshev spectral methods could be implemented as Galerkin and Tau methods, where the work takes place in the frequency space or as pseudospectral methods where the work takes place in the physical space [5]. Since most softwares are matrix oriented, an effective way to apply the above methods computationally would be in an operatorial way. That is the solution is presented as U = A−1 F , where U is the solution vector at some given points, A is a matrix, and F is some vector. Bhrawy and Alofi [2] introduces the operational matrix to the shifted Chebyshev method to generate an already fast algorithm for fractional integration in the sense that only a small number of shifted Chebyshev polynomials is needed to obtain a satisfactory result. More on this approach can be found in Weideman and Reddy [33], Trefethen [30] with some variants but still based on pseudospectral methods. Although accurate, one main disadvantage with pseudospectral methods is that the differentiation matrix is full, thus making the computation of A−1 a heavy task. However, significant computational savings can be realised via representation in the frequency domain. Trif [31] proposed a solution based on the Tau method, allowing to avoid full matrices. We shall follow in these footsteps. In this paper, we develop a time-spectral domain decomposition method based on the Tau-matrix operational approach using a differentiation matrix method on a time interval divided into disjoint domains. The obtained matrices have the advantage of being sparse upper triangular instead of full as produced by most operational matrix methods. The methodology is applied to solve systems of stiff Riccati equations that are encountered in pricing under the stochastic volatility model even in the presence of jumps. Numerical results show that our method is fast, accurate, and robust as compared to the existing Chebfun method [26], which is one the most well-known and leading packages of Matlab using spectral methods. This work is organised as follows. Section 2 gives a quick layout of the mathematical model and results under affine settings. Section 3 introduces the numerical method, and Section 4 covers the application on case studies. In the last Section 5, we draw a conclusion. 2. The mathematical model setup
Consider the financial market model M = Ω, F, P, (Ft )t≥0 , (St )t≥0 where Ω is the set of all possible outcomes of an experiment known as the sample space, F is the set of all events, i.e., permissible combinations of outcomes, P is a map F −→ [0, 1] that assigns a probability to each event, (Ft )t≥0 is a natural filtration and St a risky underlying asset price process. The triplet (Ω, F, P) is defined as a probability space. Let (Wt )t≥0 be a P-Wiener process, (Zt )t≥0 a pure jump process (a poisson process with parameter λt ), σ(t, Xt ) > 0 the volatility of the underlying asset, and µ(t, Xt ) the drift parameter. We consider a state (underlying) process Xt satisfying the following stochastic differential equation dXt = µ(t, Xt )dt + σ(t, Xt )dWt + dZt .

4

(2.1)

In order to avoid arbitrage, the price Ψ(t, Xt ) at time t of a contingent claim that pays off ΦT at maturity time T ≥ t, under an equivalent martingale measure Q, is given by  RT  Ψ(t, Xt ) = E Q e− 0 r(s,Xs )ds ΦT |Ft , (2.2) where the expectation is taken with respect to the risk-neutral measure Q. The application of Ito differentiation on equation (2.2) (see [20]) yields the partial integro-differential equation (PIDE) of the form Z 2 ∂Ψt ∂Ψt 1 > ∂ Ψt + µt + σt σt − λt Ψt + λt Ψ(xs + z)dν0 (z) = 0. (2.3) ∂t ∂Xt 2 ∂Xt2 R\0 In affine framework (see also [7]), the drift vector µt , the instantaneous covariance matrix σt σt> , the jump intensity λt and the risk free rate rt have an affine dependence on Xt : µt = K 0 + K 1 · X t ,

σt σt> = H0 + H1 · Xt , λt = l0 + l1 · Xt ,

rt = ρ0 + ρ1 · Xt ,

(K0 , K1 ) ∈ Rd × Rd×d

(H0 , H1 ) ∈ Rd×d × Rd×d×d

(l0 , l1 ) ∈ Rd × Rd×d

(ρ0 , ρ1 ) ∈ Rd × Rd×d .

Duffie et al. [10] show that under some conditions the solution of the PIDE (2.3) can be computed via the complex-valued Riccati ordinary differential equations (ODEs) defined as follows: Theorem 2.1. Under technical conditions, if the pay-off function ΦT , at maturity time T , is chosen such that ΦT = euXT , for some u ∈ Rd then Ψ is of the form Ψ(t, Xt ) = eα(t)+β(t)Xt with α and β verifying the following Riccati equation  ∂α  ∂t (t) = ρ0 − K0 · β(t) − 12 β(t)> H0 β(t) − l0 [Λ(β(t)) − 1] 

∂β (t) ∂t

= ρ1 −

K1> β(t)

−

1 β(t)> H1 β(t) 2

(2.4)

,

(2.5)

− l1 [Λ(β(t)) − 1]

with terminal conditions α(T ) = 0 and β(T ) = u. It is without a doubt that equation (2.5) is analytically difficult to solve and in general presents stiffness. We require numerical methods to come to the rescue. In the next section, we introduce the time multidomain spectral method based on an operational matrix.

5

3. Multidomain Spectral method 3.1. Chebyshev Polynomials The Chebyshev polynomial Tn (x) of the first kind is a polynomial in x ∈ [−1, 1] of degree n > 0 defined by the relation: Tn (x) = cos nθ, for x = cos θ ie. Tn (x) = cos(n arccos(x)) 2 b+a Note that by considering the shift mapping s : x → s(x) = b−a x − b−a , we can work on the interval [−1, 1] then apply the inverse shift mapping to allow the catering of any interval [a, b]. The Chebyshev polynomial has the following properties:

T0 (x) = 1, T1 (x) = x, Tn (x) = 2xTn−1 (x) − Tn−2 (x), In turn, this can be expressed in a matrix form as:   1  −2x  1    1  −2x 1     ... ... ...   1 −2x 1

T0 (x) T1 (x) T2 (x) .. . Tn (x)

The zeros of Tn are the points

xk = − cos

(k − 12 )π , n

(3.6) (3.7) (3.8)

n = 2, 3, ... 



      =    

1 −x 0 .. .

k = 1, 2, . . . , n.

0

      

(3.9)

(3.10)

The set {xk }k is termed as collocation points, also called Chebyshev points of the first kind. For any point x ∈ [−1, 1], the set {T0 (x), T1 (x), . . . } is an orthogonal basis according to the weighted inner product defined by: Z 1 f (x)g(x) √ < f, g >= dx (3.11) 1 − x2 −1 where f, g are two continuous functions defined on [−1, 1]. This means that for any polynomial of degree n > 0, there exists a unique set of coefficients {c1 , c2 , ..., cn } such that pn (x) =

n X

ck Tk (x).

(3.12)

k=0

A Chebyshev approximation of order n > 0 of a function u continuous on an interval [−1, 1] is defined as un (x) =

n X

ck Tk (x)

(3.13)

k=0

= c · T (x) 6

(3.14)

where c = (c0 , c1 , . . . , cn ) is the coefficient vector associated with the approximation un . It is usually termed as the spectral representation of fn . The set of Chebyshev coefficient vectors {c} of continuous functions on [−1, 1] is termed as the frequency space. The vector c is also termed as the spectral representation of u. For simplicity of notation, we shall write u(x) in place of un (x) to denote the Chebyshev approximation of order n of u at x. Another discrete representation of the function u is v = (u(x0 ), u(x1 ), ..., u(xn )), where x = (x0 , x1 , ..., xn ). We shall call v the physical representation of u. On the collocation points (3.10), one writes v(x) = T (x).c ! n n X X (v(x0 ), ..., v(xn )) = ck Tk (x0 ), ..., ck Tk (xn ) k=0

(3.15) (3.16)

k=0

where T is the matrix defined as follows  T0 (x0 ) T1 (x0 )  T0 (x1 ) T1 (x1 )  T = ..  . T0 (xn ) T1 (xn )

... ... .. .

Tn (x0 ) Tn (x1 ) .. .

...

Tn (xn )

Since

    

v = T c ⇒ c = T −1 v From the nature of Tk0 s, the matrix T is sparse and FFT enables it to get T −1 . 3.2. Differentiation in frequency space Recall again Equation (3.13) and differentiate it u0 (x) =

n X

ck Tk0 (x)

(3.17)

k=0

The differentiation of relation (3.8) and relation (3.7) gives T0 = T10 T20 T1 = 2 0 Tn+1 (x) = nTn−1 (x) − 2(1 − x2 )Tn0 (x) 0 0 Tn+1 Tn−1 ie. Tn = − , n = 2, 3, ... 2(n + 1) 2(n − 1)

(3.18) (3.19) (3.20) (3.21)

Inserting this back into (3.17) shows the existence of a matrix D = (dkl )0
c0k Tk (x)

=

n X n X

dkl cl Tk (x)

(3.22)

k=0 l=0

ie. c0 = Dc 7

(3.23)

where c0 stands for the spectral representation of the derivative function u0 . Moreover, D is a sparse upper triangular matrix, with the following properties  for k ≤ l,  dkl = 0, dkl = 0, if l − k is even, (3.24)  dkl = 2k, if l − k is odd.

Applying the above result recursively we get the spectral representation c(p) of the derivative at order p of u as: (3.25) c(p) = Dp c. Consequently, if we consider a general differential equation Au = f of order m with constant coefficients for which the differential operator can be written as A = L + N where L and N are respectively the linear part and the nonlinear part, then the equation can be written as Lu(t) + N u(t) = f (t) m X u(k) (t) = −N u(t) + f (t)

(3.26) (3.27)

k=0

The above equation can be represented in the frequency space as: m X k=0

implying

Dk c = −n + ˜f Ac = f c = A−1 (f)

(3.28)

P k ˜ for some A = m k=0 D ; where n, f are the spectral representation of N u and f respec˜ tively, and f = −n + f. we use the following algorithm 1: Algorithm 1 Pseudo code 1: u0 ← initial solution 2: INITIALIZE A 3: Evaluate N , and f at u0 4: u := A−1 ∗ (N + f ) 5: while ku − u0 k >  do 6: u0 ← u 7: Evaluate N , and f at u0 8: u = A−1 ∗ (N + f ) 9: RETURN u The above method has proven to struggle when T gets larger. For this reason, we suggest a domain decomposition of the interval [0, T ].

8

3.3. Multiple Domain Spectral Method In this section we consider. Ih to be a mesh on the interval [0, T ] and N to be the number of subintervals. Ih := {tn : 0 = t0 < t1 < · · · < tN = T } We denote Λn = [tn−1 , tn ], hn = tn − tn−1 and un (t) as the solution of (3.26) on the n-th element of the partition Ih , namely, un (t) = u(t),

∀t ∈ Λn ,

1 ≤ n ≤ N.

Let Mn > 0 be an integer and consider PMn to be the space of polynomials of order at most Mn built on Λn . We apply the spectral method as described in algorithm 1 to obtain a numerical solution UMn ∈ PMn on Λn . The Time Multidomain Spectral Method on the interval [0, T ] consists of a successive application of the spectral method on each Λn to obtain a global numerical solution UM (t) of (3.26) defined such that: UM (t) = UMn (t),

t ∈ Λn ,

1 ≤ n ≤ N.

where M is taken to be the smallest of the Mn ’s: that is, M = inf 0
i = 0, ..., m − 1.

(3.29)

The overall matrix A of the entire problem is then a diagonal of block matrices A(i) .     (1)  f A(1) 0 c(1)  0 A(2) 0   c(2)   f(2)       (3.30)    ..  =  ..  . . . .   .   .  . . 0 A(m) c(m) f(m)

By inversion of the matrix A(i) on each domain Λi , we obtain c(i) and therefore uMi , which is UM on Λi . In this case a global error can arise and jeopardise the convergence. But thanks be to the following theorem 3.1 that still guarantees an exponential convergence even after discretisation. Theorem 3.1. Assume that u belongs to the broken Sobolev space: u ∈ H 1 (0, T ) and u|Λn ∈ H rn (Λn ), 1 ≤ n ≤ N with integers 2 ≤ rn ≤ Mn + 1, and there exists a constant L ≥ 0 such that for any z1 and z2 , |f (z1 , t) − f (z2 , t)| ≤ L|z1 − z2 |. (3.31) Then for we have

√ 2 2πhmax L ≤ β < 1, ku − UM k2H 1 (0,T ) ≤ cβ T exp(cβ T )

N X i=1

i −2 Mi2−2ri |u|2H ri (Λ) , h2r i

(3.32) (3.33)

where cβ is a positive constant depending only on β.

Where Λn = [tn−1 , tn ], hn = tn − tn−1 and the constant Mn is the order of the Chebyshev polynomial of approximation un defined on Λn . The proof can be found in Ma [19], Wang and Mu [32]. 9

4. Applications and numerical results In this section, we apply our method to different problems found in finance and test the convergence and the efficiency of the proposed method against the existing Chebfun method [26]. In all numerical experiments, when analytical solutions are not available, we choose the solution from ODE15s with relative and absolute tolerances to be 10−14 as the benchmark solution as one of the appropriate Matlab packages for stiff problems [22]. The error E is the maximal error given by: ||E|| = ||SolBenchmark − SolN umerical ||∞ .

(4.34)

4.1. Oil price model The model for oil pricing is derived as follows. Let St be the spot price of crude oil and Vt the stochastic variance. We allow the log-stock price to be driven by p d ln St = µt dt + Vt dWts + dZts (4.35)

The stochastic variance is considered to be a square root process. In absence of jumps, the mean reverting toward the parameter v is the long run mean of V . The process reverts at a speed controlled by the parameter κ and σv is the volatility of the volatility. It measures the responsiveness to diffuse volatility shocks. The volatility is driven by the following dynamics: p ˜v (4.36) dVt = κ(v − Vt )dt + σv Vt dW t with corr[dWts , dWtv ] = ρ and that the jump process Zts of St arrives exponentially with rate λt and with size normally distributed N (¯ µs , σs ). Therefore the model is expressed as: 





 p  1 0 p d = dt + Vt dWt + dZt (4.37) ρσv 1 − ρ2 σv   Yt Considering the state process Xt = , equation (4.37) can then be written as: Vt Yt Vt

r − λt ms − 21 Vt κ(v − Vt )



dXt = µt dt + σt dWt + dZt where µt =



r − λt ms κv



+



0 − 21 0 −κ



Xt

σt σt>

and

=



1 ρσv ρσv σv2



Vt

Referring to affine settings in (2.5), we see that: ρ0 = r, ρ1 = (0, 0), l0 = λt , l1 = (0, 0), K0 =



r − λms κv



, K1 =



0 − 12 0 −κ



, H0 =



10

0 0 0 0



and H1 =

"

# . 0 0 .. 1 ρσv . . 0 0 .. ρσv σv2

The oil price is given by φt = eα(t)+β(t) x = eα(t)+β1 (t)yt +β2 (t)Vt β (t) = St 1 eα(t)+β2 (t)Vt where α and β = (β1 , β2 ) are real valued functions that satisfy the Riccati equation  ˙ = r − (r − λms )β1 (t) + κvβ2 (t) − λ(Λ(β(t)) − 1)  α(t) β˙1 (t) = 0 , (4.38)  ˙ 1 1 2 1 2 2 β2 (t) = 2 β1 (t) + κβ2 (t) − 2 β1 (t) − β1 (t)β2 (t)ρσv − 2 β2 (t)σv together with terminal conditions

α(T ) = 0 and β(T ) = (u, 0). R and Λ is defined as Λ(c) = R2 exp (c z)dν(z). In other words, Λ(c) is the moment generating function of the random vector Z taken at c. The risk neutral condition imposes 2 σs

ms = Λ(1, 0) − 1, that is ms = eµs + 2 (see Duffie et al. [10]). µs , σs ), we get Considering the jump distribution to be of Gaussian type, that is Zts = N (¯ Λ(β(s)) = Λ(u, β2 (s)) = eu¯µs +

u2 2 σ 2 s

s + 1 σ2 2 s

and ms = Λ(1, 0) − 1 = eµ

Operating the change in time τ = T − t the system of equations reduces to  α(τ ˙ ) = −r + (r − λms )β1 (τ ) − κvβ2 (τ ) + λ(g(τ ) − 1) , β˙2 (τ ) = − 12 β1 (τ ) − κβ2 (τ ) + 21 β12 (τ ) + ρσv β1 (τ )β2 (τ ) + 21 σv2 β22 (τ )

(4.39)

with α(0) = 0

β2 (0) = 0

β1 (τ ) = u

The exact solution of the above equation is: −a (1 − e−γτ ) 2γ − (γ + b) (1 − e−γτ )    2 γ+b γ+b −γτ α(τ ) = −rτ + r uτ − κv τ + 2 ln 1 − (1 − e ) σv2 σv 2γ − λ(1 + ms u)τ + λg(τ )

β2 (τ ) =

s + 1 σ2 2 s

where γ = b2 − aσv2 , and g(τ ) = τ eµ

(4.40)

(4.41)

.

For numerical illustration, we choose the parameters to have the following values, as illustrated in [18] : µs = −1.8521, σs = 1.9578, λ = 0.0294, ρ = −0.0083, r = 0, ms = 0.0642.

κ = 0.0118,

11

v = 4.6978,

σv = 0.1655294,

We compare our numerical solution with the exact solution and evaluate the error as the L∞ distance between the exact and the numerical solution. The plots in Figure 1a and 1b show that the exact and numerical solutions for variable α and β are in good agreement. We also run a comparison with another appropriate well-known numerical method, here we mention Chebfun, which is also based on spectral Chebyshev polynomials. As we vary the number of collocation points for both methods, we record in Figure 2a the error on variable α. It shows that as the number of collocation points gets larger (n > 8) TMDSM achieves better accuracy quicker than Chebfun. The same remark also holds for variable β. In figure 2b a plot of the efficiency of the TMDSM and of Chebfun is provided. It is interesting to see that we achieve accuracy of order 10−15 within 0.005s, whereas Chebfun would take close to 0.02s to achieve the same order of accuracy on the same problem with the same processor Core I5-8th Gen. This is mostly due to the structure of the matrices inherited from the differentiation in the frequency space rather than in the physical space, which was mentioned to be sparse upper triangular whereas Chebfun differentiation matrices are full, see Figure 3. More elaborately, Figure 3a exhibits two upper triangles together with lines, it represents the plot of matrix operator A. The first triangle cares for the first variable x together with its boundary condition, and the second triangle is for the second variable y also with its boundary condition. In reality these triangles contain trips of zeros, due to the relation (3.24). For a case of 100 collocation points for instance, the operator A is a (200 × 200)-matrix which normally has 40, 000 elements but in our case the total number of non-zero elements is 5398 and the rest are zeros. The sparsity coefficient is thus about 86.5%, reducing the amount of calculation considerably, A−1 . As for the structure of the Chebfun matrix in 3b we have a totally dense matrix. 0

0

αexact αTDMSM

−1

−1

−3

−1.5

−4

−2

−5

−2.5

10

20

τ

30

40

βTDMSM

−0.5

−2

−6 0

βexact

−3 0

50

(a) α

10

20

τ

30

(b) β

Figure 1: exact and numerical solutions plot for T = 50 with 2 domains

12

40

50

−2

−6

10

10

TMDSM Chebfun

TMDSM Chebfun

−4

10

−8

10 −6

10

−10

10 Log10E

Log10E

−8

10

−10

10

−12

10

−12

10

−14

10

−14

10

−16

10

−16

0

10 20 30 Number of collocation pts

10

40

(a) Convergence

0

0.005

0.01 CPU time (s)

0.015

0.02

(b) Efficiency

Figure 2: Convergence and efficiency of TMDSM vs Chebfun on α.

0 20 40 60 80 100 120 140 160 180 200 0

50

100 nz = 5398

150

200

(a) 1-Domain TMDSM matrix

(b) Chebfun matrix

Figure 3: Plots of the structure of the underlying matrix A for 1-Domain TMDSM vs Chebfun

13

4.2. Interest rate model Short rate models are helpful in pricing of assets, such as bonds. It has been shown that most short rate time series experience means reverting phenomenon that varies with time. It also appears that interest rates rt converge rather quickly towards a time varying mean level µt , which itself reverts more slowly toward an unconditional mean θ. One reason for this mean reversion is due to the inflation factor that manifests itself in the required nominal interest rate, as mentioned by Andersen et al. [1]. Another key feature is that theory predicts that an unexpected information arrival should induce a discontinuity in the return process. An announced shift in monetary policy, for instance, can induce the appearance of jumps in the interest rate process. Plus having jumps as part of the modeling can improve model description, calibration as jumps can also accommodate outliers in the short rate distribution. For this example we will restrict to jumps, which are normally distributed and uncorrelated. Thus, we consider the interest rate r to be driven by the following process: p drt = (κ1 µ ¯t − λ¯ µJ − κ1 r) dt + Vt dWtr + dZtr

(4.42)

where Wtr , Ztr are respectively a Brownian process and a pure jump process (with √ constant intensity λ) associated with the stock price process. The stock price volatility Vt is stochastic with the following dynamics, p dVt = κ2 (v − Vt )dt + σv Vt dWtv

where Wtr , Wtv are independent Brownian motions. Note in this model that the volatility is mean reverting, with mean v and κ2 as mean reversion coefficient, σv is the volatility of the volatility. Also we consider the mean value µ ¯t of the drift to be stochastic. This is due to some unusual changes, by governments, in their monetary policies. Let’s write the dynamics of µ ¯t to be driven by the following equation: √ d¯ µt = κ3 (¯ v−µ ¯t )dt + σm µ ¯t dWtm . (4.43) As it stands, a good candidate to consider as a state process is the short rate process.      rt 1 0 κ1 µ ¯t − λ¯ µJ − κ1 rt p  dt + Vt  0 σv κ2 v − κ2 Vt d  Vt  =  µ ¯t 0 0 κ3 v¯ − κ3 µ ¯t

Xt = (rt , Vt , µ ¯t ) where rt is  0 0  dWt + dZt σm

where Wt = (Wtr , Wtv , 0) is a three-dimensional Brownian motion and Zt = (Ztr , 0, 0) is a three-dimensional Poisson process in t with mean arrival rate λt . As such the entire process Xt jumps only when rt jumps. Hence the intensity arrival λt = λ and the multivariate distribution of the jump size of Zt is identical to the distribution of Ztr . Referring to affine settings from section 2 and if we consider a financial claim that pays off ΨT = eu¯·XT at time T for some constant u¯ = (u1 , u2 , u3 ) ∈ R3 then theorem 2.1 ensures that its price φt at time 0 < t < T is of the form φt = eα(t)+β(t) x 14

(4.44)

where α and β = (β1 , β2 , β3 ) are deterministic functions of t that satisfy the Riccati equation (4.38), which yields the following:  α(t) ˙ = 0 − (−λr µJ )β1 (t) − κ2 v(β2 (t)) − κ3 ηβ3 (t) − λr (Λ(β(t)) − 1)    β˙ (t) = 1 + κ β 1 1 1 , (4.45) ˙ β2 (t) = κ2 β2 (t) − 21 β12 (t) − 12 σv2 β22 (t)    ˙ 2 2 β3 (t) = −κ1 β1 (t) + κ3 β3 (t) − 12 σm β3 (t) together with terminal conditions

α(T ) = 0 and β(T ) = (0, 0, 1). Operate the rescaling time τ = T − t. Also, it is not difficult to see that the solution of β1 is: β1 (τ ) =

1 −κ1 τ 1 [e − 1] i.e. β1 (t) = [eκ1 (t−T ) − 1]. κ1 κ1

(4.46)

The simultaneous equations (4.45) then result in the following initial value problem  ˙ ) = −λr µJ β1 (τ ) + κ2 vβ2 (τ ) + κ3 ηβ3 (τ ) + λr (Λ(β(τ )) − 1)  α(τ β˙2 (τ ) = −κ2 β2 (τ ) + 21 β12 (τ ) + 12 σv2 β22 (τ ) , (4.47)  ˙ 1 2 2 β3 (τ ) = κ1 β1 (τ ) − κ3 β3 (τ ) + 2 σm β3 (τ ) with initial conditions

α(0) = 0, Zts

and β(0) = (0, 0, 1).

Once again for the case of a jump distribution of Zts being of Gaussian type, that is = N (µs , σs ) then, Λ(β(s)) = Λ(β1 , β2 (s)) = eµs β1 (s)+

2 β 2 (s) σs 1 2

For numerical implementation let’s consider the same set of parameters (see 4.48) used by Andersen et al. [1]. As such the problem presents a moderate stiffness factor. A similar problem with a severe stiffness factor is encountered in bond pricing [14]. κ1 = 1.7887, κ2 = 1.7895, κ3 = 0.2792, v = 0.000052, σv = 0.0110, η = 0.0525,

σm = 0.0459, σJ = 0.0016,

µJ = 0, λr = 3.2688.

(4.48) (4.49)

Equation (4.47) is a system of nonlinear ordinary differential equations of order 3. It can be written as: Au + N u = f (4.50) Using our Time Multidomain Spectral Method described in Section 3, we transport the equation in the frequency space and it becomes Ac = f

15

(4.51)

where c and f are spectral representations of the unknown solution vector u = (α, β2 , β3 ) and the nonlinear part f − N u, respectively. In addition, the matrix A is of the form:   D κ2 vIn κ3 ηIn  0 A =  0 D − κ2 In 0 0 D − κ3 In where In is the identity matrix of order n and D is the differentiation matrix as defined in (3.23). The nonlinear part we will write as:   −λµJ β1 (τ ) − lr [Λ(β(τ )) − 1] 1 . β (τ )2 − 21 σv2 β22 (τ ) N = 2 1 1 2 2 κ1 β1 (τ ) 2 σm β3 (τ )

A plot of the benchmark solution generated from ODE15s (with absolute tolerance 10−14 ) and TMDSM is provided in Figure 4 for each variable α, β2 and β3 . 0.02

0.09 αODE15s

0.01

0.08

αTMDSM

0

1 β3 ODE15s

β

2 ODE15s

β2 TMDSM

β3 TMDSM

0.5

0.07 0

−0.01

0.06

−0.02

0.05

−0.5

−0.03

0.04

−1

−0.04

0.03

−0.05

0.02

−0.06

0.01

−1.5

−0.07 0

1

2

τ

3

(a) α

4

5

0 0

−2

1

2

τ

3

4

5

−2.5 0

(b) β2

1

2

τ

3

4

5

(c) β3 (τ )

Figure 4: plot of the 3 variables for T = 5

As we vary the number of collocation points from n = 4, 8, 16, 32, 64, 128, 256 and 512, we record in table 1 the error on each variable α, β2 and β3 . The results show again that it only requires a few points to reach the solution with high precision. The table shows that we achieve an error of about 10−8 within 16 points. The same holds also for the other two variables. This is in line with Figure 2a observed earlier. In other words, the spectral convergence still holds even in a higher dimension case. n=4 n=8 n = 16 n = 32 n = 64 n = 128 α 2.223E-5 1.526E-7 7.591E-9 4.070E-10 2.383E-11 1.398E-12 β2 6.058E-4 8.195E-6 3.367E-7 1.768E-8 1.024E-9 6.183E-11 β3 2.890E-3 1.972E-5 9.704E-7 5.204E-8 3.042E-9 1.781E-10

n = 256 8.40E-14 3.839E-13 1.078E-12

n = 512 5.551E-15 2.703E-14 6.90E-14

Table 1: Convergence of the error of α, β2 and β3

Given that Chebfun returns the solution in 1.31 seconds, we also record in Figure 5 the accuracy achieved as the running time increases. The graph shows that for the very same 16

Efficiency on α

−4

10

−6

10

−8

Log10Error

10

−10

10

−12

10

−14

10

−16

10

0.2

0.4

0.6

0.8

1

1.2

CPU time in (s)

1.4

1.6

1.8

Figure 5: Efficiency of 2-Domain TMDSM for α

problem, the TMDSM achieves an accuracy of order 10−13 within 0.06 seconds. Again this is mostly due to the structure of the matrices produced from our differentiation for the system, compared to those from Chebfun which are full, see figure 6. Let’s increase T to 100 and introduce different partitions in 1, 2, and 4 subintervals for a total number of Chebyshev collocation points varying from 256, 512, 1024, 2048, and 4096. Let’s record the error as well as the CPU running time in the table 2. Clearly the T where N is the number of discretisation of our interval [0, T ] is uniform, that is hn = h = N subintervals. We also consider Mn = M to be constant since we generate the same number of Chebyshev points in each subinterval. Moreover it is not difficult to see that our function f here abides to the Lipschitz conditions. Therefore we should expect an exponential decay of the error. N = 1 CPU Error N = 2 CPU Error N = 4 CPU Error

n = 256 0.450 6.136E-9 0.397 2.158E-8 0.394 7.572E-8

n = 512 1.653 3.802E-10 0.765 1.327E-9 0.501 4.538E-9

n = 1024 8.875 2.368E-11 2.883 7.960E-11 1.160 2.779E-10

n = 2048 66.589 1.485E-12 16.870 5.117E-12 5.18 1.672E-11

n = 4096 607 5.337E-13 156 4.993E-13 45 1.068E-12

Table 2: Convergence and Efficiency of TMDSM with 1, 2 and 4 domains at T = 100.

The table 2 shows that as the number of collocation points gets larger (here n > 1000) on the interval [0, T ], the TMDSM tends to suffer. Indeed the matrix A gets very large, making inversion a complicated task. But if the structure of A gets more porous (lower sparsity factor) then the spectral method would still be capable of handling even larger problems without losing much in accuracy. For instance, when the total number of points is 4092 it 17

0

0

100

100

200

200

300

300

400

400

500

500

600 0

100

200

300 nz = 31396

400

500

600 0

600

(a) 1-Domain TMDSM matrix

100

200

300 nz = 16392

400

500

600

(b) 2-Domain TMDSM matrix

0

100

200

300

400

500

600 0

100

200

300 nz = 8884

400

500

600

(c) 4-Domain TMDSM matrix

(d) Chebfun matrix

Figure 6: Plots of the underlying matrix A for 1, 2, 4-Domain TMDSM vs Chebfun

takes 607 seconds (nearly 10min) for TMDSM to deliver the solution whereas it will take the same TMDSM only 45 seconds when the interval [0, 100] is split into four subintervals [0, 25] ∪ [25, 50] ∪ [50, 75] ∪ [75, 100]. This is a tremendous gain in computational time. Indeed the structure of the linear operator A plays an important role. We plot the structure of the matrix A in figure 6 using 200 collocation points. The matrix is 600 × 600 and for the case of one domain the total number of non-zero elements is 31396 (see Figure 6a) implying a sparsity of 91.2% whereas when splitting the domain into four sub-domains the number of non-zero terms reduces to 8884 (see Figure 6c) resulting in a sparsity of 97.5%. It is also remarkable to notice that the overall error did not suffer, it remained in the order of 10−12 . The method is stable. The accuracy depends not on the number of subintervals but on the total number of Chebyshev collocation points.

18

4.3. electricity pricing under affine process Energy commodity prices are important to model as their cost of storage is usually very high. Modeling their dynamics with affine jump diffusion looks reasonable in order to capture salient features of energy commodity prices, [9]. For instance, electricity cannot be stored or inventoried economically once generated. The electricity demand and supply in a bulk electricity power network has to be balanced continuously so as to prevent the network from collapsing. The most noticeable behaviour of energy commodities is the mean-reverting aspect. That is, when the price of the commodity is high, its supply tends to increase, thus putting a downward pressure on the price. When the spot price is low, the supply of commodity tends to decrease thus providing an upward lift to the price. Another feature is the presence of jumps and spikes. These happen when a massive storage of commodities is not economically available and the demand exhibits low elasticity. Also, a forced outage of a major power plant or sudden uprise of demand would either cause the supply curve to shift to the left or lift up the demand curve, therefore causing a price jump. For electricity, spikes can be observed when an unexpected breakdown of some power plant occurs, the spot price can experience an abrupt upward increase, but as the lost generator gets restored, the price falls back quickly to their normal range [9]. A suitable equation that can capture the abovementioned features is:      √ Vt 0 0 Xt κ1 θ1 − κ1 Xt p √ √ d  Vt  =  κ2 θ2 − κ2 Vt  dt +  ρ1 σ2 √Vt 1 − ρ21 σ2 Vt 0  dWt + dZt Yt κ3 θ3 − κ3 Yt ρ2 σ3 Vt 0 σ3 

(4.52)

where Xt is the log stock price, Vt is the volatility, and Yt is the product supply of raw material, Wt = (Wtx , Wtv , Wty ) is a three-dimensional Brownian motion and Zt = (Ztx , Ztv , 0) is a three-dimensional Poisson process in t with mean arrival rate λ = λx + λv + λc . The process can then be written as: dXt = µt dt + σt dWt + dZt

(4.53)

and the electricity price is of the form: φt = eα(t)+β(t) x

(4.54)

where α and β = (β1 , β2 , β3 ) are real valued functions that satisfy the Riccati equation α˙ = r −

3 X i=1

1 κi θi βi − σ32 β32 − λx [Λ(β)] 2

1 β˙ 1 = κ1 β1 − σ32 β32 2 1 β˙ 2 = κ2 β2 − [β12 + σ22 β22 + ρ22 σ32 β32 ] − λ2 [Λ(β) − 1] 2 −ρ1 σ2 β1 β2 − ρ2 σ3 β1 β3 − ρ1 ρ2 σ2 σ3 β2 β3 β˙ 3 = κ3 β3 19

(4.55) (4.56)

(4.57) (4.58)

together with terminal conditions α(T ) = 0 and β(T ) = (0, 0, 1). This problem can only be solved by using the numerical method. Operating again the rescaling τ = T − t the system becomes α˙ = −r +

3 X i=1

1 κi θi βi + σ32 β32 + λx [Λ(β) − 1] 2

1 β˙ 1 = −κ1 β1 + σ32 β32 2 1 β˙ 2 = −κ2 β2 + [β12 + σ22 β22 + ρ22 σ32 β32 ] + λ2 [Λ(β) − 1] 2 +ρ1 σ2 β1 β2 + ρ2 σ3 β1 β3 + ρ1 ρ2 σ2 σ3 β2 β3 ˙ β3 = −κ3 β3

(4.59) (4.60)

(4.61) (4.62)

with initial condition: α(0) = 0 and β(0) = (0, 0, 1). As in the previous problem we will compute the numerical solution using the set of parameters considered in [9]: κ1 = 2.17 κ2 = 3.5 κ3 = 1.8 θ1 = 3.2 θ2 = 0.85 θ3 = 0.87 σ2 = 0.8 σ3 = 0.54 ρ1 = 0.25 ρ2 = 0.2 λ1 = 6.43 λ2 = 5 µ11 = 0.23 µ12 = 0.22 µ21 = −0.14. The associated differentiation matrix is:  D −k1 θ1 −k2 θ2 −k3 θ3  0 D + k1 0 0 A=  0 0 D + k2 0 0 0 0 D + k3 the corresponding nonlinear part of the equation is:  −r + lr [Λ(β) − 1] 1 2 2  σ β 2 3 3 N =  λ2 [Λ(β) − 1] + N3 0

   

   

where N3 = 21 [β12 + σ22 β22 + ρ22 σ32 β32 ] + ρ1 σ2 β1 β2 + ρ2 σ3 β1 β3 + ρ1 ρ2 σ2 σ3 β2 β3 . Once again we run the spectral method and compare the solution against the ODE15s solution. The resulting error is captured in Figure 7a It can be seen that with a number of Chebyshev collocation points n = 16 the accuracy is relatively good. A log-log plot of the efficiency of the two methods (TMDSM and Chebfun) is provided in figure 7b. As in the previous example, the TMDSM achieves better accuracy in quicker time Chebfun. Once again the sparsity of the matrix A, (see figure 8) inherited from the 20

−4

−2

10

10 α β1

chebfun

β2

−6

10

TMDSM

−4

10

β3 −8

−6

10 Log10E

Log10E

10

−10

10

−12

−10

10

10

−14

10

−8

10

−12

1

10

2

3

10 10 number of collocation points

10

4

10

(a) Convergence of TMDSM vs Chebfun

−2

10

−1

10

0

10 Log10CPU time

1

10

2

10

(b) Efficiency of TMDSM vs Chebfun

Figure 7: Comparison of the convergence and efficiency of TMDSM against Chebfun

sparsity of the differentiation matrix D obtained in each subinterval, is the main cause of such difference in efficiency. Here we considered 200 collocation points generating therefore 800 × 800 matrix. Using a single domain, the number of non-zero terms is 41994 resulting in a sparsity factor of 93.4%. With a two domain, Figure 8b shows 21988 non-zero terms in A, that is a sparsity of 96.5%. For the case of four domains, the sparsity increases to 98.1%. 5. Conclusion We have designed a spectral method that accommodates differential equations arising from financial models of the affine type. The affine structure of these financial models is used to avoid solving the PIDE but rather to solve a system of Riccati equations. The proposed numerical method to solve theses Riccati equations presents an operational matrix based on Chebyshev polynomials and the solution is obtained in the frequency domain. In doing so, the original problem is transformed into an iterative system of algebraic equations making it easier to solve. In addition the time domain is discretised into small domains allowing for spectral convergence to still hold. Three numerical examples are implemented and solutions are compared to numerical solutions from Chebfun, a leading package available for solving differential equations using the spectral method. The result from the error calculation shows that our method still maintains its spectral convergence even as we increase the time space interval. In addition, the method shows robustness and competitiveness compared to other numerical methods on the field, here ODE15s and Chebfun. The method can be applied in the entire class of affine models with jumps. Consequently it can be applied to the pricing of many other financial derivatives such as swaptions, coupon-bearing bonds, captions, currency options, etc., to mention a few. For further possible applications it would be interesting to see how TMDSM would perform when affine structure is not taken into account, that is, directly on the PIDE.

21

0

0

100

100

200

200

300

300

400

400

500

500

600

600

700

700

800 0

200

400 nz = 41994

600

800 0

800

(a) 1-Domain TMDSM matrix

200

400 nz = 21988

600

800

(b) 2-Domain TMDSM matrix

0 100 200 300 400 500 600 700 800 0

200

400 nz = 11976

600

800

(c) 4-Domain TMDSM matrix

(d) Chebfun matrix

Figure 8: Plots of the underlying matrix A for 1, 2, 4-Domain TMDSM vs Chebfun

Declaration of Competing Interest This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue. The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript. 6. References References [1] Torben G Andersen, Luca Benzoni, and Jesper Lund. Stochastic volatility, mean drift, and jumps in the short-term interest rate. Northwestern University, Chicago, 2004. 22

[2] AH Bhrawy and AS Alofi. The operational matrix of fractional integration for shifted chebyshev polynomials. Applied Mathematics Letters, 26(1):25–31, 2013. [3] T Bj¨ork. Arbitrage Theory in Continuous Time. Oxford University Press, second edition, 2004. [4] WILLIAM Bourke, Bryant McAvaney, KAMAL Puri, and R Thurling. Global modeling of atmospheric flow by spectral methods. Methods in computational physics, 17:267–324, 1977. [5] John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001. [6] Claudio Canuto, M Yousuff Hussaini, Alfio Quarteroni, A Thomas Jr, et al. Spectral methods in fluid dynamics. Springer Science & Business Media, 2012. [7] Rama Cont and Peter Tankov. Financial Modeling With Jump Processes. Chapman and Hall/CRC, 2004. [8] Germund Dahlquist and Ake Bjorck. Numerical Methods in Scientific Computing: Volume 1, volume 103. Siam, 2008. [9] Shijie Deng. Stochastic models of energy commodity prices and their applications: Meanreversion with jumps and spikes. University of California Energy Institute Berkeley, 2000. [10] D Duffie, J Pan, and K Singleton. Transform analysis and asset pricing for affine jumpdiffusions. Econometrica, (68):1343–1376, 2000. [11] Matthias Ehrhardt and Ronald E Mickens. A fast, stable and accurate numerical method for the black–scholes equation of american options. International Journal of Theoretical and Applied Finance, 11(05):471–501, 2008. [12] Houde Han and Xiaonan Wu. A fast numerical method for the black–scholes equation of american options. SIAM Journal on Numerical Analysis, 41(6):2081–2095, 2003. [13] S.L. Heston. A closed-form solution for options with stochastic volatility with applications to bond and currency options. The Review of Financial Studies, 6(2):327–343, 1993. [14] Shirley J Huang and Jun Yu. On stiffness in affine asset pricing models. Available at SSRN 676203, 2004. [15] Shirley Jun Ying Huang. Implementation of general linear methods for stiff ordinary differential equations. The University of Auckland, New Zealand, 2005. [16] M Yousuff Hussaini and Thomas A Zang. Spectral methods in fluid dynamics. Annual review of fluid mechanics, 19(1):339–367, 1987.

23

[17] Dimitri Komatitsch and Jeroen Tromp. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical journal international, 139 (3):806–822, 1999. [18] Karl Larsson and Marcus Nossman. Jumps and stochastic volatility in oil prices: Time series evidence. Energy Economics, 33(3):504–514, 2011. [19] Heping Ma. Chebyshev–legendre spectral viscosity method for nonlinear conservation laws. SIAM Journal on Numerical Analysis, 35(3):869–892, 1998. [20] Claude Rodrigue Bambe Moutsinga, Edson Pindza, and Eben Mare. Homotopy perturbation transform method for pricing under pure diffusion models with affine coefficients. Journal of King Saud University-Science, 30(1):1–13, 2018. [21] Edgard Ngounda, Kailash C Patidar, and Edson Pindza. Contour integral method for european options with jumps. Communications in Nonlinear Science and Numerical Simulation, 18(3):478–492, 2013. [22] D Omale, PB Ojih, and MO Ogwo. Mathematical analysis of stiff and non-stiff initial value problems of ordinary differential equation using matlab. International journal of scientific & engineering research, 5(9):49–59, 2014. [23] E Pindza, KC Patidar, and E Ngounda. Robust spectral method for numerical valuation of european options under merton’s jump-diffusion model. Numerical Methods for Partial Differential Equations, 30(4):1169–1188, 2014. [24] Edson Pindza, Kailash C Patidar, and Edgard Ngounda. Implicit-explicit predictorcorrector methods combined with improved spectral methods for pricing european style vanilla and exotic options. Electronic Transactions on Numerical Analysis, 40:268–293, 2013. [25] Edson Pindza, Francis Youbi, Eben Mar´e, and Matt Davison. Barycentric spectral domain decomposition methods for valuing a class of infinite activity l´evy models. Discrete & Continuous Dynamical Systems-S, 12(3):625–643, 2019. [26] R. B. Platte and N. Trefethen. Chebfun: A new kind of numerical computing. Progress in Industrial Mathematics at ECMI 2008. Mathematics in Industry, (15):69–87, 2010. [27] Simon S.Clift and Peter A.Forsyth. Numerical solution of two asset jump diffusion models for option valuation. Applied Numerical Mathematics, (58):743–782, 2008. [28] Karel J.in ’t Hout and Jari Toivanen. Adi schemes for valuing european options under the bates model. Applied Numerical Mathematics, (130):1343–1376, 2018. [29] B. Taruvinga, B. Kang, and C.S. Nikitopoulos. Quantitative finance research centre, university of technology, sydney. Applied Numerical Mathematics, (394):1–43, 2018. [30] Lloyd N Trefethen. Spectral methods in MATLAB, volume 10. Siam, 2000.

24

[31] Damian Trif. Matrix based operational approach to differential and integral problems. 2011. [32] Zhong-qing Wang and Jun Mu. A multiple interval chebyshev-gauss-lobatto collocation method for ordinary differential equations. Numerical Mathematics: Theory, Methods and Applications, 9(4):619–639, 2016. [33] J.A. Weideman and S.C. Reddy. A matlab differentiation matrix suite. ACM Transactions on Mathematical Software (TOMS), 26(4):465–519, 2000. [34] Wuming Zhu and David A. Kopriva. A spectral element approximation to price european options. ii. the black-scholes model with two underlying assets. Journal of Scientific Computing, (39):323339, 2009.

25