Multigrid method for pricing European options under the CGMY process

: We propose a fast multigrid method for solving the discrete partial integro-di ﬀ erential equations (PIDEs) arising from pricing European options when the underlying asset is driven by an inﬁnite activity L´evy process. We consider the CGMY model whose kernel singularity gets worse when the parameter Y approaches two. Due to the integral term, the discretization matrix is dense. In order to obtain an e ﬃ cient multigrid method, we apply a ﬁxed point iteration as a smoother for multigrid. In each smoothing step, we only need to solve a sparse matrix corresponding to the di ﬀ erential operator and compute a matrix-vector product involving the integral operator by a fast Fourier transform (FFT). We prove that the ﬁxed point iteration smoother is e ﬀ ective reducing the high frequency components. Moreover, we also prove a two-grid convergence of the multigrid method by a local mode analysis. We demonstrate the e ﬀ ectiveness of the multigrid method by solving the option pricing equation under the CGMY model with ﬁnite and inﬁnite variation processes.


Introduction
The standard Black-Scholes model for pricing options assumes the underlying asset price satisfies the geometric Brownian motion model. However, it has been observed in practice that market prices often have large jumps that do not follow the geometric Brownian motion. Jump diffusion models have been proposed to more accurately capture the actual market behavior; see e.g., [1]. Recently, models based on the Lévy process have become popular in the financial literature [1][2][3][4][5][6][7]. Option pricing under the exponential Lévy process with finite activity [5,[8][9][10][11][12] and infinite activity [13][14][15][16][17] has been extensively studied. Under this model, pricing of the European option requires solving a partial integrodifferential equation (PIDE), which consists of differential and integral terms in the equation.
Numerical solution of PIDEs is challenging since the discretization of the integral term will give rise to a dense matrix. In addition, the singularity of the Lévy measure inside the jump integral may lead to a less accurate numerical approximate solution. (More discussion about the singularity of the jump integral will be given in the next section.) In the literature, the jump process can be classified into finite activity, infinite activity with finite variation, and infinite activity with infinite variation, where the singularity increases from the first to the later cases. In [14], the jump integral was split into local and nonlocal parts with the local term computed using implicit time stepping and the nonlocal term computed using explicit time stepping. In [13,18], an integration by parts technique was used to transform the integral part into a weakly singular Volterra equation. A collocation method was then applied to the CGMY (Carr-Geman-Madan-Yor) process [2] in the case of infinite activity and finite variation. In [16], a wavelet Galerkin finite element technique was applied to an infinite activity process.
In this paper, we are concerned with the numerical solution method proposed by Wang et al. [19]. This approach is able to achieve second order accuracy under an infinite activity and finite variation process, and a better than first order accurate for an infinite variation process. Moreover, it is applicable to different types of options such as European and American, non-constant coefficient PIDEs, nonconstant volatility, and various types of initial and boundary conditions. In each time step, a fixed point iteration is performed to compute the option value. In this approach, an FFT is used to efficiently compute the jump integral and more importantly, it avoids inverting a full dense matrix arising from the discretization of the jump integral. However, the convergence of the fixed point iteration can be very slow. In fact, it can be proved that the convergence rate deteriorates as the number of grid points increases (Y larger than 1) or the singularity of the Lévy measure increases. It will be made more precise in the next section.
Recently, it has been applied to solving American option equations [36]. Multigrid methods have been developed for option pricing in the literature, for instance in [37][38][39]. However these methods do not consider the CGMY or any jump models.
The success of a multigrid technique depends on smoothing and its interaction with the coarse grid correction. Our idea is to apply the fixed point iteration as a smoother in the multigrid framework, rather than to use it as a solver, for solving the PIDE arising from option pricing. We will perform a smoothing analysis and show that the fixed point iteration can effectively reduce high frequency errors. Moreover, using a local mode analysis, we provide a two-grid convergence analysis of the method. Direct discretization is used for the coarse grid matrices. Linear interpolation and full weighting are used for inter-grid transfer.
We remark that the proposed multigrid method can be used with the penalty method for pricing American options. The discretized linear system is almost the same as that of the European option with some of the diagonal entries modified as a result of the early exercise constraint. Due to the similarity, we will focus only on European option in this paper.
In Section 2, we describe the PIDE for option pricing under the Lévy process and derive the linear system resulting from a finite difference method. In Section 3, we present a multigrid method for solving the linear system. Moreover, we provide a smoothing analysis for the proposed smoother and a convergence analysis for the two-grid method. Finally, numerical results are given in Section 4 to demonstrate the fast convergence of the proposed multigrid method by solving PIDEs under a Lévy process with finite activity, infinite activity with finite variation, and infinite activity with infinite variation. Compared with the fixed point iteration which may take tens to hundreds or even thousands of iterations to converge, our proposed multigrid method converges within 10 iterations as the Y parameter of the CGMY model increases from 0 (mildly singular) to close to 2 (highly singular).

Model equation
In Black-Scholes modelling with jumps, the value of the European option price, V, satisfies the following partial integro-differential equation (PIDE) [1]: where r is the risk-free interest rate, q the continuous dividend yield, and σ the volatility. The parameter τ is the time from expiry; i.e., τ = 0 denotes the expiry time T and τ = T denotes the start of the option. The evolution of the underlying risky asset price S is driven by a Lévy process whose Lévy measure ν satisfies |y|<1 y 2 ν(dy) < ∞, |y|≥1 ν(dy) < ∞.
Under the CGMY model [2], ν is defined as follows are the indicator variables. The parameter C > 0 is the measure of the overall level of activity, G ≥ 0 and M ≥ 0 control the rate of exponential decay on the left and right of the Lévy density, and Y < 2 describes the behavior of the Lévy density in the neighborhood of zero where the density tends to infinity. The parameter Y is of interest mathematically and numerically since its value determines the regularity of ν. If Y < 0, the measure ν integrates to a finite value yielding a process of finite activity. If Y ∈ [0, 1], the process displays infinite activity but finite variation since |y|<1 y 2 ν(dy) < ∞. If Y ∈ (1, 2), the process is said to have infinite activity and infinite variation.

Discretization
There have been a number of discretization methods for solving PIDEs [13,14,16,18]. Here, we will consider the method proposed in [19], which is easy to implement and applicable to non-constant coefficients. To be self contained, we present the major steps here and refer the readers to [19] for details. Meanwhile, notations are introduced which will be used in the smoothing analysis in Section 3.1.
The jump integral in (2.1) defined on [−∞, ∞] is first approximated by a finite domain [y min , y max ] with appropriately chosen y min and y max so that the error of approximation is negligibly small. The interval is divided into N subintervals [y j − ∆y/2, y j + ∆y/2], j = 0, 1, . . . , N − 1, where ∆y = (y max − y min )/N. Note that the Lévy measure, ν(y) in (2.3) goes to infinity as y goes to 0, and it is only second moment integrable in general; see (2.2). The computation of the approximate integral is splitted into three parts: I Ω 0 , I Ω 1 and I Ω 2 , corresponding to the three regions , , Notice that ν(y) is smooth on Ω 2 but exhibits singular behaviours on Ω 0 and Ω 1 . Using the discretization techniques introduced in [19] to handle the singularity, the integral can be approximated bȳ andγ(y j ) is the value of the following integral γ(y j ) computed by a quadrature rule: j y j +∆y/2 y j −∆y/2 y 2 ν(y)dy; if y j ∈ Ω 1 y j +∆y/2 y j −∆y/2 ν(y)dy; if y j ∈ Ω 2 0; if y j ∈ Ω 0 .
As mentioned above, ν(y) is smooth on Ω 2 . The quadrature rule can be just any standard numerical integration method such as the composite trapezoidal rule. On Ω 1 , a special quadrature rule is needed to capture the singular behaviour of ν(y). The idea is to choose quadrature weights such that the integral is exact for polynomials of up to degree 3. Thus the accuracy of the quadrature rule is fourth order which is needed to compensate for the 1/y 2 j factor when y j is near 0. Combining (2.1) and (2.4), the equation to be solved can be written as We see that the complex integral is simplified to a discrete sum. A semi-Lagrangian method is used to discretize the first order term. Let [S 1 , . . ., S m ] be a set of grid points for S and τ n be the n-th time step. Also let V n ≈ V(S , τ n ) be the approximate value at (S , τ n ). The semi-Lagrangian method traces node S at τ = τ n+1 along the characteristic trajectory back to τ = τ n at a point S * , which in general is not a grid node S j . Let V n * be the option value at S * , which is approximated by interpolating the nearby values of V n . In order to achieve second order accuracy, a quadratic interpolation Φ n+1 is used and the approximate value is denoted by Φ n+1 V n . Note that the superscript n + 1 for Φ n+1 emphasizes that the interpolation operator depends on time τ n+1 since the location to be interpolated, S * , is derived from S at τ n+1 .
Applying a fully-implicit timestepping with step size ∆τ along the characteristic direction gives: We note that the Crank-Nicolson scheme can also be used to obtain second order accuracy in time.
However, since there is no difference in applying the multigrid method for solving the discrete linear system arising from the two timestepping schemes, we will focus on the fully-implicit discretization. The computation of (QV n+1 ) can be expensive in general. However, it can be accurately approximated by an FFT method [12], and the approximate value is denoted by (BV n+1 ) . Applying central finite differencing to discretize (LV n+1 ) , the fully discrete system can be written as: where α and β are given by , for = 2, . . . , m − 1. When = 1, we set α = β = 0 at S 1 = 0; and when = m, we set V n+1 m equal to the relevant Dirichlet boundary condition.
In matrix form, the discrete system can be rewritten as: where I is the identity matrix, L is defined such that and B is a dense matrix whose entries satisfy the following properties [12] The matrix-vector product BV n+1 can be computed efficiently by using an FFT.

Fixed point iteration
Each time step requires solving the linear system (2.7), where the coefficient matrix, denoted by A, is dense since B is dense. In [19], the authors proposed a fixed point iteration which splits A into The solution is computed by solving iteratively where V n+1,k is an approximate solution to V n+1 after k iterations. This method only involves inverting a sparse matrix A 1 and the matrix-vector product with the dense matrix B, which can be computed efficiently.
It can be proved that the fixed point iteration converges and the convergence rate, ρ, is bounded by [19]: Although the fixed point iteration is convergent, the convergence rate can be very slow. Note that λ is a function of the density function ν(y). It can be estimated by for Y > 0. As a result, as ∆y approaches 0.

Multigrid method
As shown in (2.11), the convergence of the fixed point iteration can be very slow when Y > 0 and ∆y is small. This arises in the infinite activity case for the Lévy process. The situation gets worse when Y approaches 2 in the infinite variation case.
We propose an efficient multigrid method for solving the linear system (2.7). The idea is to apply the fixed point iteration as a smoother in the multigrid context. While the convergence rate of the fixed point iteration can be slow, especially when Y gets close to 2, it turns out that the fixed point iteration is an effective smoother which damps the high frequency.
A V-cycle multigrid method is used to solve the PIDE (2.6). We remark that other cycles such as W-cycle, F-cycle, or full multigrid can also be used. But they are generally more expensive. We find that V-cycle is the more effective option in this case.
In a two-grid setting, the algorithm consists of three major steps. We start with an initial guess for the solution,V k h ≈ V n+1 , on the fine grid. The first step is pre-smoothing. A small number of the fixed point iterations (typically ν 1 = 1 or 2) is applied to compute an approximate solution of (2.7), where S denotes the fixed point iteration given by (2.9), and A h and f h are the fine grid matrix and right-hand side from (2.7), respectively. The role of the smoother is to damp away the high frequency components of the error, thus resulting in a smooth error. The smoothing property is analyzed in Section 3.1. The second step is coarse grid correction. Note that the error, e h , on the fine grid satisfies the equation, where r h = f h − A hṼh is the residual on the fine grid. Since the error is smooth after pre-smoothing, we can solve it accurately and efficiently on the coarse grid with grid size H, where A H is the discretization matrix on the coarse grid and r H = Rr h is the restriction of r h given by a restriction operator R. Here we use full weighting for R; i.e., Note that Dirichlet conditions are used for the boundary points (see [19]). Thus the residual values are zero on the boundary points which will be used in the restriction near the boundary.
After e H is obtained, we update the fine grid solution by where P is a prolongation operator which interpolates the coarse grid error e H onto the fine grid. Linear interpolation is used for P; i.e., The last step is to apply smoothing again, where ν 2 is usually 1 or 2.
In a multigrid setting, the coarse grid equation (3.1) is solved by applying the two-grid algorithm recursively.
The success of the multigrid method hinges on the effectiveness of the smoother. In the next section, we perform a smoothing analysis for the fixed point iteration.

Smoothing analysis
To simplify the analysis, we assume that an evenly spaced grid along the log S coordinate is used together with periodic boundary conditions. By a change of variable x = log S , we define V(x, τ) = V(e x , τ). Then (2.6) can be written as the following PIDE with constant coefficients, where V ⊗γ denotes a correlation product; i.e. (V ⊗γ) i = N−1 j=0V (x i , y j )γ(y j ). The PIDE is discretized similarly as described in Section 2.1 where the unknowns are now V n+1 . More precisely, we have where the index j * is defined such that when applying the semi-Lagrangian discretization, the grid point x j at time τ n+1 is traced back to x j * at time τ n . The value of V n j * ≡ V(x j * , τ n ) is computed by an upwind quadratic interpolation. After rearranging terms, (3.3) becomes where α = (σ 2 +σ)/(2(∆x) 2 ). Then, the corresponding fixed point iteration (2.9) is given in component form by is an approximate solution to V n+1 after k fixed point iterations. Note that the choice of the grid points for the upwind quadratic interpolation depends on the sign of the coefficient for V x in (3.2). Without loss of generality, assuming the coefficient is negative, we have where ψ's are the Lagrangian basis functions and x j * lies between x i−p−1 and x i−p for some integer p. Substituting (3.6) into (3.5), we obtain V n+1,k+1 j The fixed point iteration is used as smoother in our multigrid method. We will estimate the smoothing factor of the fixed point iteration by applying the local Fourier analysis [20] to (3.7).
be the error of V n+1,k . Let C k µ and G µ be the symbols of E k andγ, respectively, where µ is the Fourier mode; i.e. C k µ and G µ are the Fourier transform of E k andγ, respectively [20]. Also, let h = 1/N be the mesh size and i = √ −1. Then (3.7) becomes C k+1 Note the negative subscript of G −µ is a result of the correlation operator. The left-hand side of (3.8) can be simplified as Similarly, the right-hand side of (3.8) can be written as where η j = j−p m= j−p−2 ψ j,m e −( j−m)µπhi . LetŜ µ be the symbol of the smoother. Then one can easily show thatŜ µ = C k+1 µ /C k µ ; i.e., It has been proved that |η j | ≤ 1 [40, Section 5.2]. Furthermore, since G µ is the symbol ofγ, by (2.5), we have |G −µ | ≤ λ. Hence, we can boundŜ µ as |Ŝ µ | = |∆τG −µ ||η j | |1 + (r + λ + 4α sin 2 (µπh/2))∆τ| ≤ λ∆τ 1 + (r + λ + 4α sin 2 (µπh/2))∆τ . (3.11) Let ζ be the smoothing factor of the fixed point iteration, which is defined as Considering the high frequency Fourier modes; i.e., N/2 ≤ µ ≤ N, it follows from (3.11) that To bring more insight into the formula, we estimate the upper bound for ζ by noting that Considering the limit where ∆x ≈ 0 and keeping ∆τ fixed, we have In the Lévy process under the CGMY model, we always have Y < 2. Hence for small ∆x, the smoothing factor for the fixed point iteration is small, which is fundamental to achieve an effective multigrid method.

Two-grid convergence analysis
In the previous section, we have shown by a local Fourier analysis that the fixed point iteration smoother is effective damping high frequency errors. In this section, we will perform a two-grid convergence analysis based on the local Fourier analysis. For easy exposition, we assume there is no post-smoothing.
The iteration matrix of the two-grid iteration can be written as where L h and L H are the fine and coarse grid operator, respectively, P is the prolongation, R is the restriction, and S is the smoother. For two-grid analysis [21], it is customary to pair up the low and high frequency modes such that the Fourier transformed matrix,M, is a block diagonal matrix with each block a 2 × 2 matrix,M µ , which is given bŷ Here c µ = cos(µπh)/2, s µ = sin(µπh)/2, and h = ∆x. By (3.10), we haveŜ (3.14) whereσ = (σ 2 +σ)/2. Also, we have assumed the interest rate r to be zero, just to simplify the algebraic expression.
A similar calculation yieldŝ where η j is defined similarly as η j on the coarse grid. Details can be found in Appendix A.

Numerical results
In this section, we present results of numerical experiments to demonstrate the effectiveness of the proposed multigrid method for pricing European options under the CGMY model. The PIDE is discretized by the semi-Lagrangian method in the spatial domain and implicit Euler's method for time. The discrete linear system of equations is solved by the proposed multigrid method. We also include numerical results using a fixed point iteration for comparison.
In the multigrid procedure, a V-cycle is used with one pre-and one post-smoothing where the smoother is a fixed point iteration. The coarse grid operators are obtained from direct discretization with linear interpolation and full-weighting restriction for intergrid transfer. The stopping criterion is when the relative residual l 2 -norm was less than 10 −8 .
As explained in Section 2, the fixed point iteration (cf. Section 2.1) used in the literature converges slowly when the grid size is small and becomes worse for Y ≈ 2. We will show the multigrid convergence results where Y changes from 0 to 1.98.
Example 1: In this example, Y = 0. In this case, the CGMY process is also known as variance gamma. The input parameters for the PIDE and CGMY model are given in Table 1.  Table 2 presents the convergence of the proposed multigrid method (MG) for different grid sizes. In this example, the total number of time steps is 500, and so ∆τ = 10 −3 , which is typical for pricing option values with reasonable accuracy. Multigrid converges very fast in only 2 iterations. The number of fixed point iterations (FP) is relatively constant which shows that the convergence is independent of the grid size under a Variance Gamma process, as explained in (2.11). Example 2: In this example, Y = 0.6442. When 0 < Y < 1, the Lévy process has infinite activity but finite variation. The input parameters are given in Table 3. The convergence results are shown in Table 4. In this case, the number of fixed point iteration increases for smaller grid size, supported by the formula in (2.11) with Y > 0. The multigrid convergence remains constant, independent of the grid size. In this example, we also use a larger ∆τ to test its effect on the convergence. In practice, a large time step size may be necessary when pricing options with long expiry dates; e.g., T = 10. However, a larger time step size will generally result in a more ill-condition discrete linear system. Table 4 shows that the fixed point iteration numbers increase dramatically, by a factor between 5 to 7. The number of multigrid iteration only increases mildly.
Example 3: When 1 < Y < 2, this corresponds to a Lévy process with infinite activity and infinite variation. Table 6 shows the convergence results for pricing a call option for Y = 1.4 and Y = 1.98, with input model parameters given in Table 5. When Y = 1.4, the discrete problem becomes much more difficult to solve, as can be seen from the number of fixed point iteration which increases from hundreds to thousands as the grid size decreases. The multigrid iteration remains very efficient with mesh independent convergence rate. The case where Y = 1.98 is numerically challenging since the kernel function of the Lévy process becomes very singular. The fixed point iteration converges very slowly. The multigrid iteration deteriorates slightly from 4 to 9, but still converges very quickly.
Example 4: In this example, we illustrate the convergence behaviour of multigrid, fixed point iteration and BiCGstab with the fixed point iteration as preconditioner [19]. The results are shown in Figure 1.
Here we consider the case where Y = 1.4. The convergence behaviour is similar for other values of Y and they are not shown here. We can see that multigrid converges quickly with a constant rate and the convergence is essentially unchanged for different grid sizes. Fixed point iteration converges at a much slower rate and the rate of convergence deteriorates as the grid size h decreases. The convergence of preconditioned BiCGstab is somewhere between the two. It is much better than fixed point iteration but not as good as multigrid. Also, the convergence is not as smooth and uniform as the other two, and it is mesh dependent. Table 7 shows the rates of convergence for multigrid and fixed point iteration, which are computed by taking the average rates of the last 3 iterations. The convergence of BiCGstab is not very uniform and we do not include its rate of convergence here.
In Sections 3.1 and 3.2, we presented a local Fourier analysis to analyze the smoother and the two-grid convergence. It would be interesting to show how the theoretical analysis results compare with the actual convergence. However, both the smoothing factor and two-grid convergence factor include a quantity η j (3.9) which comes from the interpolation operator during the semi-Lagrangian discretization. This quantity is difficult to compute and it varies at different grid points, grid sizes and time steps. As such, we were unable to compute the smoothing factor and the spectral radius ofM. Instead, we compute an upper bound for the smoothing factor and the infinity norm ofM as described in (3.11), (3.19) and (3.20). They are shown in columns 3 and 4 in Table 7. Since they are upper bounds, they do not match very well with the actual convergence of multigrid. Nevertheless, we can see that the smoothing factor decreases with h as discussed in Section 3.1 and the two-grid factor is bounded by 0.75 as shown in Theorem 3.2. We also remark that the discrepancy may be due to the fact that the underlying PIDE is a nonlocal operator where the local mode analysis may not be the right tool for analysis.

Conclusion
In this paper, we have proposed a fast multigrid method for solving the discrete equations from a finite difference discretization of the partial integro-differential equations arising from pricing European options under the CGMY process. The fixed point iteration originally proposed in [19] converges well when the density kernel is not singular (Y ≈ 0) but very slowly when the kernel becomes more singular (Y ≈ 2). While the fixed point iteration is not an efficient solver by itself, we have proved that it is effective removing high frequency error components. By using the fixed point iteration as a smoother for multigrid, we have obtained a fast multigrid solver which converges within 10 iterations for different values of Y ranging from 0 to 1.98. In addition, we have provided a two-grid convergence analysis for the multigrid method.
The current paper focuses primarily on European options. For American options under the CGMY process, the corresponding PIDE is nonlinear with a linear complimentarity condition. The proposed multigrid method would not be directly applicable since it is designed for solving linear systems. An interesting idea is to combine the techniques in this paper with a nonlinear multigrid method such as Full Approximation Scheme (FAS), which will be a future work. Another future work is to perform better local mode analysis of the method, such as three-grid analysis, in order to obtain improved quantitative results.