ADI schemes for valuing European options under the Bates model

This paper is concerned with the adaptation of alternating direction implicit (ADI) time discretization schemes for the numerical solution of partial integro-differential equations (PIDEs) with application to the Bates model in finance. Three different adaptations are formulated and their (von Neumann) stability is analyzed. Ample numerical experiments are provided for the Bates PIDE, illustrating the actual stability and convergence behaviour of the three adaptations.


Introduction
The traditional asset price model in financial option valuation theory is the geometric Brownian motion. It is well-known that this model assumes constant volatility while the market prices of options clearly indicate varying volatility. For the valuation of options with long maturities, stochastic volatility models like the Heston stochastic volatility model [18] are a common means to introduce such variability. For short maturities, however, pure Brownian motion based models such as the Heston model often require excessively large volatilities to explain the market prices of options. A modern way to resolve this is to incorporate jumps in the asset price model, like the classical Merton jump-diffusion model [36] does. The Bates model [6] combines the Merton model and the Heston model. This asset price model includes both jumps and stochastic volatility and it is therefore popular for valuing options with short as well as long maturities. Under the Bates model, a two-dimensional parabolic partial integro-differential equation (PIDE) can be derived for the values of European-style options with the spatial variables representing the underlying asset price and its instantaneous variance [9]. The differential operator is of the convection-diffusion type. The integral operator, which stems from the jumps, couples all option values in the asset price direction. This paper deals with the numerical solution of the Bates PIDE. We follow the method of lines approach, where the PIDE is first discretized in the spatial variables and the resulting semidiscrete system of ordinary differential equations (ODEs) is solved by applying a suitable implicit time discretization scheme. Due to the nonlocal nature of the integral operator, semidiscretization leads to a large, dense system matrix. For an implicit time discretization scheme one thus faces two challenges: the efficient treatment of the two-dimensional convection-diffusion part and the efficient treatment of the dense integral part. In this paper we focus on operator splitting methods, which handle these two parts separately in an effective manner. The following reviews methods proposed in the literature.
For the time discretization of the two-dimensional partial differential equation (PDE) arising from the Heston model, various alternating direction implicit (ADI) schemes were studied in [23]. ADI and other directional splitting methods, see e.g. [21], are highly efficient. In each time step they yield a sequence of unidirectional linear systems, which can be solved very fast using LU factorization. Since the original papers [13,14,38], the literature on ADI schemes for time-dependent convection-diffusion equations is extensive. In financial option valuation, correlations between the underlying stochastic processes (asset price, variance) give rise to mixed spatial-derivative terms in the pertinent PDEs. First-order ADI schemes have been investigated for such PDEs in [34,35] and second-order ADI schemes in [11,17,23,24,25,26,28,29,37], where the mixed derivative part is always treated in an explicit fashion. Under the Merton model a one-dimensional PIDE for European option values holds. Application of a classical implicit time discretization scheme such as Crank-Nicolson or BDF2 leads in every time step to the solution of a linear system with a dense matrix. An iterative method requiring only multiplications by a dense matrix and solutions to tridiagonal systems can often solve these linear systems with a few iterations [1,12,39,42]. Alternatively, a splitting method for the time discretization of the Merton PIDE was proposed in [2]. This method, of the Peaceman-Rachford type, treats the integral part in an implicit manner and employs FFT to reduce computational cost. An implicit-explicit (IMEX) method that treats the sparse diffusion part implicitly and the dense integral part explicitly was investigated in [10]. The pertinent method is only first-order. Second-order IMEX methods such as the IMEX midpoint and IMEX-CNAB (Crank-Nicolson Adams-Bashforth) schemes were subsequently investigated in [16,33,40] and IMEX schemes of the Runge-Kutta type were examined in [7].
The IMEX-CNAB method has been applied to the Bates PIDE in [41]. For the solution of the linear system in each time step, both LU factorization and an algebraic multigrid method were considered. An adaptive finite difference discretization combined with the IMEX-CNAB method and LU factorization was constructed in [44]. In [4,5] a first-order directional splitting method was used together with Richardson extrapolation with the aim to obtain second-order accuracy. Directional splitting methods of the Strang type were studied in [8,43]. In [31] the authors employed Strang splitting between the convection-diffusion and integrals parts and applied a secondorder ADI scheme to solve the two-dimensional convection-diffusion substeps in this approach; see also [30]. A novel, direct adaptation of ADI schemes to two-dimensional PIDEs has recently been proposed in [32]. Here the integral and mixed derivative parts are treated jointly in an explicit fashion.
For more information on operator splitting methods like ADI and IMEX schemes in finance we refer to [22,27]. A good general reference is [20].
In this paper we consider the adaptation of ADI schemes for the time discretization of two-dimensional PIDEs with the Bates model as the leading example. We formulate three types of adaptations: two of these are novel and one agrees with that from [32]. The obtained operator splitting methods are all of second-order and treat both the integral and mixed derivative parts in an explicit manner. At each time step they require only one or two multiplications with the dense matrix corresponding to the integral part. Hence, the resulting methods are computationally very efficient.
An outline of this paper is as follows. Section 2 describes the Bates model, the PIDE for European option values, and its spatial discretization. In Section 3 the three adaptations of ADI schemes to PIDEs are formulated and their stability is analyzed. Numerical experiments for these adaptations applied to the Bates PIDE are presented in Section 4. The conclusions end the paper in Section 5.

Bates model and its semidiscretization
We consider in this paper the Bates model [6], which merges the Heston stochastic volatility model [18] and the Merton jump-diffusion model [36] for the asset price S τ ≥ 0. Here the instantaneous variance V τ ≥ 0 follows a CIR process with meanreversion level η > 0, mean-reversion rate κ > 0 and volatility-of-variance σ > 0. The underlying Brownian motions for S τ and V τ have correlation ρ ∈ [−1, 1]. The jump process is a compound Poisson process with intensity λ > 0. The probability of a jump from S τ to S τ Y with random variable Y > 0 is defined by the log-normal probability density function where γ and δ are given real constants with δ > 0 that are equal to the mean and standard deviation, respectively, of the normal random variable log(Y ). Let r ≥ 0 be the risk-free interest rate and let ε denote the mean relative jump size, given by ε = e γ+δ 2 /2 − 1. Consider an European-style option with maturity time T . Under the Bates model, if S T −t = s and V T −t = v, then the option value u(s, v, t) satisfies the PIDE whenever s > 0, v > 0, 0 < t ≤ T . The initial condition for (2.1) is where φ denotes the pay-off function, which specifies the value of the option at expiry. In this paper we shall consider vanilla put options. The pay-off function is then φ(s) = max{K − s, 0} with given strike price K.
For feasibility of the numerical solution, the unbounded spatial domain is truncated to [0, S max ] × [0, V max ] with sufficiently large S max and V max chosen such that the error caused by this truncation is negligible. We pose the Dirichlet and Neumann boundary conditions and at the boundary v = 0 it is assumed that the PIDE (2.1) holds, compare [15]. For the spatial discretization on the computational domain [0, S max ] × [0, V max ] a smooth, nonuniform Cartesian grid is chosen with meshes in the s-and v-directions of the type considered in [17]. This grid has a large fraction of points in the neighbourhood of the important location (s, v) = (K, 0). Let the grid points in the s-direction be denoted by 0 = s 0 < s 1 < · · · < s m 1 = S max and the grid points in the v-direction by For the partial derivatives of u with respect to s we apply the standard central finite differences and For the partial derivatives of u with respect to v the analogous finite differences are considered, with the exception of the boundary v = 0, where we apply the forward formula For the mixed derivative ∂ 2 u/∂s∂v the central nine-point formula is used that is obtained by successive application of the standard central finite differences for the first derivative in the s-and v-directions.
For the discretization of the integral term in (2.1), we employ linear interpolation for u between consecutive grid points in the s-direction, and u(s, v, t) is approximated to be zero whenever s ≥ S max . After some calculations, this leads to the quadrature rule where for d = 0, 1 with erf(·) denoting the error function. This same approximation of the integral term was used in [39], for example. Here we adopt the convention log(0) = −∞ and erf(∞) = 1.
Collecting the semidiscrete approximations to the grid point values Here A is a given m×m-matrix and G(t) is a given m-vector with m = (m 1 −1)(m 2 +1). The vector G(t) depends on the Dirichlet boundary data in (2.3). Let I denote the m × m-identity matrix. The matrix A can be written as where D represents the convection-diffusion part of (2.1) and J represents the integral.
. It can be shown that all entries of J are nonnegative and all its row sums are bounded from above by one. In particular this implies that the spectrum of J belongs to the complex unit disk.
The ODE system (2.9) is complemented with an initial vector U(0) that is given by the values of the pay-off function φ at the spatial grid points where, in view of the nonsmoothness of φ at the strike K, the values at the grid points (s i , v j ) with s i the point lying closest to K are replaced by the cell average value 3 Operator splitting methods

Adaptation of ADI schemes to PIDEs
For the efficient temporal discretization of the semidiscrete system (2.9) we study in this paper operator splitting methods. To this purpose, the matrix A is decomposed into four matrices, Here A 1 , respectively A 2 , corresponds to the finite difference discretization of all spatial derivatives in the s-direction, respectively v-direction, minus half of the (r+λ)I matrix.
Next, A (D) 0 represents the finite difference discretization of the mixed derivative term in (2.1) and A (J) 0 = λJ represents the finite difference discretization of the integral term.
We note that incorporating (r + λ)I in the above manner turns out to be beneficial for stability. The vector G(t), containing the boundary contributions, is decomposed analogously to A. It is convenient to define for 0 ≤ t ≤ T and V ∈ R m , The basis for the operator splitting approach under investigation in this paper is the class of ADI schemes, which is well-established in the literature for option valuation PDEs, without integral term. As a prominent representative of this class we select the Modified Craig-Sneyd (MCS) scheme, which was introduced in [29].
Let θ > 0 be a given parameter. Let step size ∆t = T /N with given integer N ≥ 1 and temporal grid points t n = n∆t for n = 0, 1, . . . , N. If the integral term is absent, and the MCS scheme generates successively, in a one-step manner, approximations U n to U(t n ) for n = 1, 2, . . . , N by The MCS scheme has classical order of consistency (that is, for fixed nonstiff ODEs) equal to two for any parameter value θ. The particular choice θ = 1 2 yields the Craig-Sneyd (CS) scheme, introduced in [11].
In the MCS scheme the F 0 part is treated in an explicit fashion and the F 1 and F 2 parts are treated in an implicit fashion. In each step of the scheme four linear systems need to be solved, involving the two matrices I − θ∆tA j for j = 1, 2. Since these two matrices are essentially tridiagonal (up to permutation), the linear systems can be solved very efficiently by using LU factorizations, which can be computed upfront. Notice further that Y 0 and Y 0 can be combined in (3.2), so that just two evaluations of F 0 are required per time step.
We study three adaptations of the MCS scheme to PIDEs. The first adaptation is straightforward. It consists of (3.2) with F 0 = F (J) 0 +F (D) 0 as defined above. The integral term is then conveniently treated in an explicit way, along with the mixed derivative term. This type of adaptation has recently been proposed in [32], where the (related) Hundsdorfer-Verwer ADI scheme, instead of the MCS scheme, was considered.
The second adaptation of the MCS scheme to PIDEs is obtained by treating the integral term in a one-step, explicit, second-order fashion directly at the beginning of each time step: We mention that the scheme defined by the first three lines of (3.3) followed by U n = Y 2 has recently been studied, in a different context, in [3] for problems without integral term.
In the third adaptation of the MCS scheme to PIDEs, the integral term is treated in a two-step, explicit, second-order fashion directly at the start of each time step: It can be shown that each of the adaptations (3.2), (3.3), (3.4) of the MCS scheme retains classical order of consistency equal to two for any value θ. The mixed derivative and integral terms are always handled in an explicit manner. In (3.2) and (3.3) the integral term is dealt with in an explicit trapezoidal rule fashion, whereas in (3.4) it is treated in a two-step Adams-Bashforth fashion. An advantage of the adaptation (3.4) is that only one multiplication with the dense matrix A (J) 0 is required per time step, instead of two for (3.2) and (3.3). Treating the integral term in a twostep Adams-Bashforth fashion has recently been advocated in [40,41] where IMEX multistep schemes, without directional splitting, were considered for option valuation PIDEs.
In the following section we analyze the stability of the three adaptations of the MCS scheme.
Clearly, if w 0 = 0, then R, S, T 1 reduce to the same rational function and T 0 vanishes.
We investigate the stability of the three processes (3.6), (3.7), (3.8). By definition of A j it holds for j = 1, 2 that where µ j represents an eigenvalue of the matrix corresponding to the finite difference discretization of all spatial derivatives in the j-th spatial direction. Let z j = µ j ∆t for j = 1, 2. In the stability analysis of ADI schemes for two-dimensional convectiondiffusion equations with mixed derivative term, the following condition on z 0 , z 1 , z 2 plays a key role: where ℜ designates the real part. It is satisfied [28] in the model (von Neumann) framework, with semidiscretization by second-order central finite differences on uniform, Cartesian grids and periodic boundary condition. Consider the following, related condition on w 0 , z 0 , z 1 , z 2 : For the integral term there holds λ 0 = λν, where ν ∈ C represents an eigenvalue of the matrix J . A natural assumption on ν is |ν| ≤ 1, compare Section 2. For w 0 , z 0 , z 1 , z 2 , z 1 , z 2 defined above, we have the following useful lemma.
In the literature on the stability analysis of the MCS scheme pertinent to PDEs with mixed derivative terms, a variety of favourable results has been derived on the stability bound |R(0, z 0 , z 1 , z 2 )| ≤ 1 under the condition Obviously, (3.15) is weaker than (3.14). Hence, by virtue of Lemma 3.1, one directly arrives at positive stability results for the adaptation (3.2) pertinent to two-dimensional PIDEs with mixed derivative term. The following theorem gives two main results, which are relevant to diffusion-dominated and convection-dominated problems, respectively.
Letting successively ξ → ∞ and x → ∞, the right-hand side tends to the value 1−1/θ 2 . The modulus of this must be bounded from above by 1, which yields θ ≥ 1 2 √ 2. (b) Taking w 0 = −x ≤ 0, an analogous derivation as in (a) leads to the requirement that 1 + 1/θ 2 must be bounded from above by 1, which obviously does not hold.
Clearly, in Theorem 3.3 for the adaptation (3.3) the popular values θ = 1 3 and θ = 1 2 are excluded. We remark that the condition w 0 ∈ R, w 0 ≥ 0 is relevant to the case where the eigenvalues ν of J are real and lie in the interval [0, 1]. This is often seen to be numerically fulfilled in our applications, compare Section 4. For the adaptation (3.4), stability is determined by power-boundedness of the companion matrix This is established by studying the roots of the characteristic polynomial For any given point (w 0 , z 0 , z 1 , z 2 ) the recurrence relation (3.8) is stable if and only if the root condition holds: both roots ζ of P have a modulus of at most one, and those with modulus equal to one are simple. If T 0 = T 0 (w 0 , z 0 , z 1 , z 2 ) and T 1 = T 1 (w 0 , z 0 , z 1 , z 2 ) are real, then by the Schur criterion the root condition is satisfied if and only if and there is no multiple root with modulus equal to one. (a) If θ ≥ 1 3 , then (3.8) is stable whenever w 0 , z 0 , z 1 , z 2 ∈ R satisfy (3.14) and w 0 ≥ 0.
Proof. From (3.11), (3.12) it is directly clear that Assume |R 0 | ≤ 1. Then, by Lemma 3.5, and the remainder of the proof is the same as that of Theorem 3.6.

Numerical examples
In this section we examine by numerical experiments the stability and convergence behaviour of the three adaptations (3.2), (3.3), (3.4) of the MCS scheme in the application to the semidiscretized Bates PIDE described in Section 2. We consider the four parameter sets for the Bates model and European put option listed in Table 1. Case I has been considered in e.g. [41,44] and Case II in [4,5]. Case III stems from the recent monograph [30]. Case IV is new and has been constructed as a specifically challenging example, where in particular λT is large and the Feller condition 2κη > σ 2 is violated. For the MCS scheme the common parameter values θ = 1 3 and θ = 1 2 are selected. Recall that the latter choice yields the Craig-Sneyd (CS) scheme. In addition, we study the reduced version of each of (3.2), (3.3), (3.4) that ends directly with the computation of Y 2 and sets U n = Y 2 . These can be viewed as adaptations of the well-known Douglas (Do) scheme (cf. e.g. [19,20,28]) and will be applied with θ = 1 2 . We investigate the global temporal discretization error at time t = T = N∆t, where index k is such that U k (T ) and U N, k correspond to the spatial grid point (s i , v j ). Clearly, the temporal error under consideration is defined via the maximum norm. The set {(s, v) : 1 2 K < s < 3 2 K, 0 < v < 1} forms a natural region of interest for practical applications. For the truncation of the spatial domain, S max = 8K and V max = 5 are taken, and the semidiscrete solution vector U(T ) in e is approximated to high accuracy by applying a suitable time stepping method with a very small step size. Figure 1 shows the eigenvalues of the matrix J in the four cases of Table 1 when m 1 = 2m 2 = 200. They all lie in or are close to the real interval [0, 1] in Cases I, II, III, whereas they are dispersed in the complex unit disk in Case IV.   Table 1.
each figure concerns the i-th adaptation (i = 1, 2, 3). The results corresponding to the MCS scheme with θ = 1 3 and θ = 1 2 are indicated by green and blue circles, respectively, and those corresponding to the Do scheme with θ = 1 2 by red diamonds. In Figures 2, 3 we observe that in all but one instance the temporal errors always remain below a moderate value and are (essentially) monotonically decreasing as N increases. In Case IV the second adaptation yields large errors for each ADI scheme if N 50. We attribute this to a lack of unconditional stability, which may reveal itself when λT is large. Except for this instance, for each given adaptation, the MCS scheme with θ = 1 3 always yields temporal errors that are smaller than, or approximately equal to, those for the CS and Do schemes, and it always shows a smooth, second-order convergence behaviour. Additional numerical experiments on finer spatial grids, such as with m 1 = 2m 2 = 400, show that the temporal errors for all adaptations of the MCS scheme with θ = 1 3 are visually identical, indicating that second-order convergence holds in a favourable stiff sense.
In Cases I, II the temporal errors obtained with the CS and Do schemes are often relatively large for N 100. This is due to the nonsmoothness of the initial function and can be resolved by using backward Euler damping (Rannacher time stepping). This is computationally expensive, however, and does not alter our above conclusion.
Finally, we note that Figures 2, 3 suggest that the temporal error constant increases as λT increases. Additional numerical experiments, where λ and T are varied and all else is kept fixed, confirm this.

Conclusion
In this paper we have studied three adaptations of the MCS scheme to two-dimensional PIDEs that are second-order consistent in the classical ODE sense. They treat both the integral term and mixed derivative term explicitly and require per time step the solution of four linear systems with essentially tridiagonal matrices. The first and second adaptations, (3.2) and (3.3), handle the integral term in a one-step fashion and employ two multiplications per time step with the dense matrix corresponding to this term. The third adaptation, (3.4), deals with the integral term in a two-step fashion and uses just one multiplication per time step with the pertinent dense matrix. In view of the stability results obtained in Section 3 and numerical experiments for the Bates PIDE performed in Section 4, we recommend the third adaptation of the MCS scheme with parameter value θ = 1 3 , directly followed by the first adaptation of this scheme.