An IMEX-Scheme for Pricing Options under Stochastic Volatility Models with Jumps

. Partial integro-diﬀerential equation (PIDE) formulations are often preferable for pricing options under models with stochastic volatility and jumps, especially for American-style option contracts. We consider the pricing of options under such models, namely the Bates model and the so-called stochastic volatility with contemporaneous jumps (SVCJ) model. The nonlocality of the jump terms in these models leads to matrices with full matrix blocks. Standard discretization methods are not viable directly since they would require the inversion of such a matrix. Instead, we adopt a two-step implicit-explicit (IMEX) time discretization scheme, the IMEX-CNAB scheme, where the jump term is treated explicitly using the second-order Adams–Bashforth (AB) method, while the rest is treated implicitly using the Crank–Nicolson (CN) method. The resulting linear systems can then be solved directly by employing LU decomposition. Alternatively, the systems can be iterated under a scalable algebraic multigrid (AMG) method. For pricing American options, LU decomposition is employed with an operator splitting method for the early exercise constraint or, alternatively, a projected AMG method can be used to solve the resulting linear complementarity problems. We price European and American options in numerical experiments, which demonstrate the good eﬃciency of the proposed methods.

B819 [44,49] which automatically generate a sequence of coarser problems. The projected algebraic multigrid (PAMG) method proposed in [52] is a generalization of these algebraic methods for LCPs. This is one of the methods that we employ for the solution of LCPs.
The other method that we will use is the operator splitting method proposed in [25] and employed for the Heston model in [27,28]. This method approximates the LCP as a system of linear equations and a set of simple decoupled linear complementarity problems. A popular alternative approximation is the penalty method which was considered for the Heston model by Zvan, Forsyth, and Vetzal in [54] and since then by many authors for various option pricing models. For the block tridiagonal systems of linear equations many different methods can employed. In this paper, we will show that these systems can be solved very efficiently using a modern sparse direct solver when the time-step is constant. In this case, an LU decomposition of the coefficient-matrix needs to be formed only once. For two-dimensional problems like the underlying one, George showed in his 1973 paper [22] that the decomposition can be formed using O(m 3 ) operations and the solution each time-step requires O(m 2 log m) operations, where m is of the same order as the number of grid-points in both directions. As the number of time-steps is usually also of the same order, i.e., ∼ m, the computational complexity of the time-stepping is greater than the one required by forming the LU decomposition.
Explicit treatment of the jumps leads to matrix-vector multiplications with a full matrix. These multiplications can be performed efficiently by employing FFT based implementations. For a logarithmically uniform grid, the jump matrix under the SVJ model is a Toeplitz matrix, which can be embedded into a circulant matrix, as in [2], for example. The multiplication can then be computed by an FFT and an inverse FFT, which require only O(m log m) operations for each grid line. This approach is more involved for the SVCJ model, where FFTs need to be performed in two directions. Again under the SVCJ model the jump matrix can be embedded into a block circulant matrix with circulant blocks (BCCB); see [18], for example. The FFTs in both directions require O(m 2 log m) operations. Here we describe and employ this approach with the SVCJ model.

Governing equations.
First, we consider a model with stochastic volatility and jumps in returns described by Bates [5]. Under this model the behavior of the asset value s and its variance v is described by the coupled stochastic differential equations Here μ B is the drift rate defined by μ B = r − q − λξ B , where r ≥ 0 is the riskfree interest rate and q ≥ 0 is the continuous dividend yield. The jump process J is a compound Poisson process with intensity λ > 0 and J + 1 has a log-normal distribution p(y) with the mean in log y being γ and the variance in log y being δ 2 , i.e., the probability density function is given by . The parameter ξ B is defined by ξ B = e γ+δ 2 /2 − 1. The variance v has mean level θ, κ is the rate of reversion to the mean level of v, and σ is the volatility of the variance v. The two Wiener processes w 1 are w 2 have the correlation ρ. Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php By combining derivative pricing arguments from [11,19] for the Bates model, we can obtain the PIDE (formulated in forward time) where u is the price of a European option and τ = T − t is the time to expiry. The operators L B D and L B I are defined as the differential part (including the term −(r+λ)u) and the integral part of (2.2), respectively. The initial condition for (2.2) is defined by where g is the pay-off function which gives the value of the option at the expiry. Next, we allow for jumps in the volatility and study SVCJ. Then we have the correlation between jumps in returns and variance. The two-dimensional jump process (J s , J v ) is an R × R + -valued compound Poisson process with intensity λ > 0. The distribution of the jump size in variance is assumed to be exponential with mean ν. Conditional on a jump of size z v in the variance process, J s + 1 has a log-normal distribution p(z s , z v ) with the mean in log z s being γ + ρ J z v . This gives a bivariate probability density function defined by p(z s , . As in [15,18], we assume that the price of a European option under the SVCJ model can be obtained as the solution to the PIDE ∂u ∂τ = 1 2 vs 2 ∂ 2 u ∂s 2 + ρσvs (2.5) The initial condition is defined by (2.3). For computational reasons we truncate the unbounded domain to (s, v, τ ) ∈ (0, s max )×(0, v max )×(0, T ]. We impose the boundary conditions u(0, v, τ) = e −rτ g(0) and ∂u ∂v (s, v max , τ) = 0. For put options we impose u(s max , v, τ) = 0. Furthermore, for integrations we extend u for v > v max as u(s, v, τ ) = u(s, v max , τ). For put options, we extend u for s > s max as u(s, v, τ ) = 0. On the boundary v = 0, the PIDEs (2.2) and (2.5) themselves can be posed as a boundary condition; see [16] for a discussion on this boundary and its treatment.

Linear complementarity problem for American options.
In the following, L D and L I denote the differential operator L B D or L S D , and the integral operator L B I or L S I , respectively. With this notation both (2.2) and (2.5) can be expressed as

Discretization.
3.1. Computational grid. We will use a computational grid that is uniform in τ and nonuniform in s and v. The spatial grid is denoted (s i , v j ), i = 1, . . . , m s , j = 1, . . . , m v . We employ grid generating functions s : There are many ways to choose the functions s and v. Here we use the quadratic functions with the coefficients a, b, c, and d defined by the following conditions. These quadratic functions offer one of the simplest ways to construct refined grids. Also other grid generating functions can be used equally well with the numerical methods described in the following. In order to have the end points at s max and v max the conditions s(1) = s max and v(1) = v max have to hold. We choose s to map the point p K to the strike price K, i.e., s(p K ) = K. When this point satisfies p K < K/s max the grid is refined in the interval [0, K]. A similar condition can be defined also for v, but instead we require the grid to be α times finer at 0 than at v max . This leads to the condition v (0) = αv (1). These four conditions define the coefficients a, b, c, and d. At a grid point s i = s((i − 1)/(m s − 1)) the grid step to the left is denoted by Δs − i = s i − s i−1 and the grid step to the right is denoted by Δs + i = s i+1 − s i . The grid steps around v j are denoted in the same way.

IMEX-discretization in time.
In [47] a number of IMEX-discretization methods for option pricing problems under jump-diffusion models are examined. The authors found that the IMEX-CNAB method produced the smallest error among the methods they studied. This method can be seen as a modification of the popular Crank-Nicolson method with a second-order accurate explicit treatment for the integral operator. For these reasons, we will employ and promote this method here.
We start by taking a small even number 2Ñ of Rannacher-style smoothing Euler steps [42,23] with the time-step Δτ/2 given by This is followed by IMEX-CNAB steps defined by (3.3)

Discretization of the differential operator.
Here we construct a sevenpoint finite difference discretization for the differential operator (3.4) Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php We start with the mixed derivative and assume that the correlation ρ is negative. A similar seven-point finite difference discretization can also be constructed for a positive ρ; see [26,27], for example. An alternative approach to construct finite difference stencils is considered in [29].
Using a weighted average of the approximations in (3.6) we obtain Using (3.7) in (3.4), we obtain the approximation Second-order finite difference discretizations of the second-derivatives are defined by Similarly, second-order finite differences for the first-order derivatives are defined by When a finite difference method is employed, spurious oscillations might occur when the discretization matrix of −L D does not lead to a so-called M -matrix. Sufficient conditions for an M -matrix are that it is strictly diagonally dominant with positive diagonal elements and it has nonpositive off-diagonal elements. To ensure that the off-diagonals are to be nonpositive, we add artificial diffusion when necessary. With the second-order differences (3.10) and (3.11), this leads to the modified diffusion coefficientsc When the original coefficients c ss and c vv are positive, this addition of artificial diffusion is equivalent to using one-sided first-order finite differences to discretize a part of the first-order derivatives. Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php In order not to add excessive diffusion to the discretization, it is desirable that the coefficients c ss and c vv are nonnegative. For a general weight w ∈ [0, 1], this leads to the bounds for the grid step sizes Δs − i and Δs + i . These bounds can be quite stringent when the correlation ρ approaches −1.
Note that the diagonal element of the discretization matrix is r + λ minus the sum of the off-diagonal elements. Thus, the matrix is strictly diagonally dominant with positive diagonals when r + λ > 0 and the off-diagonal elements are nonpositive.

Discretization of the integral operator.
3.4.1. Bates' model. The integral term can be written as , make the change of variable ζ = z + x and study the value of the integral at the point x i , The interval (x min , x max ) is chosen to be so large that the impact of the two last integrals of (3.15) is negligible. For put options x max = log s max is a natural choice asū(ζ, v, τ) = 0 for ζ ≥ x max due to the extension of u for s > s max . The first part of the integral is evaluated using the trapezoidal quadrature rule on an equidistant grid in x with spacing Δx and m x grid-points in (x min , x max ) giving Note that the above approximation includes the additional terms Δx 2ū (x min , v, τ)p(x min − x i ) and Δx 2ū (x max , v, τ)p(x max − x i ), which are also assumed to be negligible. Computing allĪ i , i = 1, . . . , m x , can be accomplished by the matrix-vector multiplication defined byĪ The computation ofĪ = T mxū can be performed by first computingĨ = C 2mx−1ũ , whereũ = ū 1ū2 · · ·ū mx−1ūmx 0 · · · 0 T . ThenĪ is given by the m x first elements inĨ.
The circulant matrix C 2mx−1 can be decomposed as We summarize the computation of the integral in (3.14) as follows: • Interpolate values from the s i grid-points to equidistant points x i between x min and x max . • Compute and embed the matrix T mx into a circulant matrix C Mx .
• ComputeĪ using the algorithm described above using FFTs.
• InterpolateĪ i to the original grid-points s i .

SVCJ model. The integral term can be written as
Now, make the changes of variables ζ = x + z x and η = v + z v and study the value of the integral at the point (x i , v j ), Similarly to the treatment of the integral in the Bates model, we choose a rectangle (x min , x max ) × (v j ,v max ), which is large enough so that integrating over it gives a sufficiently good approximation for (3.18). The first part of the integral is evaluated Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php using the two-dimensional generalization of the trapezoidal rule on an equidistant grid in x and in v with spacing Δv andm v grid-points giving (3.19) By definingĪ = Ī 1Ī2 · · ·Ī mx−1Īmx T ,ū = ū 1ū2 · · ·ū mx−1ūmx T , we can computeĪ byĪ = Tm v mxū , where Tm v mx is a block-Toeplitz matrix with Toeplitz blocks (BTTB-matrix) defined by where   We summarize the computation of the integral in (3.18) as follows: • Interpolate values from grid-points in (s, v) to equidistant points x i between x min and x max and equidistant points v j between 0 andv max . • Compute the BTTB-matrix Tm v mx and embed this matrix in a BCCB-matrix C Mv Mx . • ComputeĪ using the algorithm described above using FFTs.
• InterpolateĪ i,j to the original grid in (s, v).

Numerical experiments.
In this section we price European and American put options under the Bates and the SVCJ models. We compare two alternative approaches to solving the discretized systems: the algebraic multigrid (AMG) method, and the LU decomposition method. The operator splitting (OS) method [25,28] is employed to enable the LU decomposition method for American options as well. Similarly, the projected algebraic multigrid (PAMG) method is adopted instead of the AMG for American options. These combinations lead to a total of eight different cases.
The used AMG implementation is based on the Ruge-Stüben algorithm [44] and the PAMG method is described in [52]. This algorithm generates automatically the coarse grid problems using only the left-hand side matrix. In the considered option prices problems, each coarsening reduces the size of the system to be between onethird and one-half. For the used grids the number of levels generated by the multigrid methods varies between 7 and 16. The AMG methods employ a multigrid V-cycle. The smoother for moving downwards and upwards is one (projected) Gauss-Seidel iteration over all points, followed by one more iteration over the so-called fine-points (F-points). For the LU decompositions, we employed the UMFPACK version 5.6.1 [14]. For FFTs, we employed FFTW version 3.3.2 [21]. We performed the experiments on a Mac laptop with 2 GHz Intel Core i7 processor and 8 GB of memory.
The shared model parameters employed in the numerical experiments are listed in Table 1. In addition, under the SVCJ model the correlation between jumps is set as ρ J = −0.5, and the mean of the exponentially distributed jump sizes in variance is set as ν = 0.2. The truncation boundary in the s-direction is s max = 4K = 400. The truncation boundary in the v-direction for the Bates model is v max = 0.5 and for the SVCJ model v max = 3. These truncations are chosen experimentally so that the error caused by them should not be larger than the discretization errors on the considered grids. The coefficients for the grid generating function s in (3.1) are defined by the parameter p K = 7/16. The coefficients for the grid generating function v are defined by α = 2 for the Bates model and α = 4 for the SVCJ model. The truncation boundaries for the integrations are defined by x min = log s 2 − 1 8 (log s max − log s 2 ), x max = log s max , andv max = 5 4 v max . The number of grid-points in the equidistant integration grid is twice as many in each direction compared to the computational grid. Reference prices, listed in Table 2, were computed using the (P)AMG method on fine grids defined by m s = 4097, m v = 2049, N = 513 for the Bates model and by Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php m s = 4097, m v = 4097, N = 513 for the SVCJ model. Figure 1 illustrates that the pricing errors of the (P)AMG and LU(+OS) approaches are almost identical. Option prices, ratios of consecutive errors, average iterations, and CPU times for the European options under the Bates and SVCJ models are listed in Tables 3 and 5 for the AMG approach, and in Tables 4 and 6 for the LU decomposition approach. Similarly, the numerical results for American options under the Bates and SVCJ models are reported in Tables 7 and 9 for the PAMG approach, and in Tables 8 and 10 for the LU+OS approach. The seven refinements in the plots of Figure 1 correspond to the seven grids described in Tables 3-10.
The ratios of consecutive errors are computed using the l 2 norm with respect to the reference prices. The ratios are, on average, of second order. Thus, this clearly suggests that the IMEX scheme is also second-order accurate as expected. With the LU+OS approach the accuracy is slightly better for the American option prices at 90 and 100 on the finest grids. The operator splitting method yields more accurate prices near the early exercise boundary in our experiments. Due to this the ratios for the LU+OS approach are better on the finest grids and it seems to converge faster on the refinement levels 6 and 7 in Figure 1.
As the grid is refined, the LU decomposition approach is slightly faster than the (P)AMG method. In both cases, however, penny-accurate prices can be obtained in a fraction of a second under the Bates model, and in about six seconds under the SVCJ model. On finer grids computing prices under the SVCJ model is up to 15 times more time consuming than under the Bates model. The two main reasons for this is the need to perform 2d FFTs instead of 1d FFTs and larger computational domain/grids for the SVCJ model. Under the SVCJ model, performing the integrations require an order of magnitude more time than solving the discretized PDEs or associated LCPs. Thus, the proposed IMEX scheme is much faster than any iterative procedure requiring multiple integrations per time step. Figure 1 and Tables 2-10 show that the prices of American options are slightly Downloaded 10/06/14 to 171.64.163.125. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php less accurate than for the European options on the same grid. The ratio of the l 2 errors of American and European options is larger on finer grids. For the (P)AMG and LU(+OS) approaches these ratios are less than 2 and 1.3, respectively. On this respect the methods behave very similarly under the Bates and SVCJ models. Thus, the discretizations are also rather accurate for American options. The CPU times for pricing European and American options are essentially the same in the tables. Thus, the PAMG and LU+OS methods offer very efficient way to solve the LCPs resulting from American options.

Conclusions.
In this work we have considered the numerical pricing of options under the Bates model and the SVCJ model. For the time discretization we employed the second-order accurate IMEX-CNAB scheme which treats the differential operator implicitly and the integral term explicitly. This way we avoid having to solve systems of equations with dense matrices. The matrix-vector multiplications arising from the explicit jump terms can be computed efficiently using FFTs.
European options lead to linear systems of equations, which under the IMEX-CNAB scheme can be solved efficiently by employing AMG method or LU decomposition. For American options a common approach is to formulate the pricing problem as an LCP. Here we solve these LCPs by employing either PAMG or operator splitting method in combination with LU decomposition. Both these methods were employed to price European/American options under the Bates/SVCJ model, rendering a total setup of eight combinations.
Numerical experiments show that the (P)AMG and LU(+OS) approaches produce almost identical prices. Since the operator splitting method does not essentially reduce the accuracy in the solution and the LU decomposition can be precomputed prior to the time-stepping, the LU(+OS) methodology turns out to be surprisingly accurate and fast. With the proposed methods pricing American options is essentially equally fast and only a bit less accurate than pricing European options. Thus, the methods are especially efficient for pricing American options. To the best of our knowledge, this is the first paper where American options under the SVCJ model are priced.