Fractal approximation of chaos game representations using recurrent iterated function systems

1.   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 integro-differential 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, non-constant 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.

Instead of fixed point iteration, we propose an efficient multigrid method for solving the linear system arising from the discretization of the PIDE using the approach in ^[19]. Multigrid has shown to be a powerful and one of the most efficient numerical techniques for solving elliptic partial differential equations (PDEs) ^[20,21,22]. Its convergence rate is often independent of the mesh size ^[21]. Sophisticated smoothing ^[20,23,24], coarsening, and interpolation ^{[25,26,27,28]} techniques have been developed. The technique has also been extended to convection dominated problems ^{[29,30,31,32]}, and hyperbolic equations ^[33,34,35]. 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).

2.   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 $\sigma$ the volatility. The parameter $\tau$ is the time from expiry; i.e., $\tau = 0$ denotes the expiry time $T$ and $\tau = 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 $\nu$ satisfies

Under the CGMY model ^[2], $\nu$ is defined as follows

where

are the indicator variables. The parameter $C > 0$ is the measure of the overall level of activity, $G \geq 0$ and $M \geq 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 $\nu$. If $Y < 0$, the measure $\nu$ integrates to a finite value yielding a process of finite activity. If $Y \in [0, 1]$, the process displays infinite activity but finite variation since $\int_{|y| < 1} y^2 \nu(dy) < \infty$. If $Y \in (1, 2)$, the process is said to have infinite activity and infinite variation.

2.1.   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 $[-\infty, \infty]$ 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-\Delta y/2, y_j+\Delta y/2]$, $j = 0, 1, \ldots, N-1$, where $\Delta y = (y_{\max}-y_{\min})/N$. Note that the Lévy measure, $\nu(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_{\Omega_0}$, $I_{\Omega_1}$ and $I_{\Omega_2}$, corresponding to the three regions

Notice that $\nu(y)$ is smooth on $\Omega_2$ but exhibits singular behaviours on $\Omega_0$ and $\Omega_1$. Using the discretization techniques introduced in ^[19] to handle the singularity, the integral can be approximated by

where

and $\hat{\gamma}(y_j)$ is the value of the following integral $\gamma(y_j)$ computed by a quadrature rule:

As mentioned above, $\nu(y)$ is smooth on $\Omega_2$. The quadrature rule can be just any standard numerical integration method such as the composite trapezoidal rule. On $\Omega_1$, a special quadrature rule is needed to capture the singular behaviour of $\nu(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_j^2$ 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$, $\ldots$, $S_m$] be a set of grid points for $S$ and $\tau^n$ be the $n$-th time step. Also let $V_{\ell}^n \approx V(S_{\ell}, \tau^n)$ be the approximate value at $(S_{\ell}, \tau^n)$. The semi-Lagrangian method traces node $S_{\ell}$ at $\tau = \tau^{n+1}$ along the characteristic trajectory back to $\tau = \tau^n$ at a point $S_{{\ell}^*}$, which in general is not a grid node $S_j$. Let $V_{{\ell}^*}^n$ be the option value at $S_{{\ell}^*}$, which is approximated by interpolating the nearby values of $V_{\ell}^n$. In order to achieve second order accuracy, a quadratic interpolation $\Phi^{n+1}$ is used and the approximate value is denoted by $\Phi^{n+1}V_{\ell}^n$. Note that the superscript $n+1$ for $\Phi^{n+1}$ emphasizes that the interpolation operator depends on time $\tau^{n+1}$ since the location to be interpolated, $S_{{\ell}^*}$, is derived from $S_{\ell}$ at $\tau^{n+1}$.

Applying a fully-implicit timestepping with step size $\Delta \tau$ along the characteristic direction gives:

where

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 $(\mathcal{Q} V^{n+1})_{\ell}$ can be expensive in general. However, it can be accurately approximated by an FFT method ^[12], and the approximate value is denoted by $(\mathcal{B}V^{n+1})_{\ell}$. Applying central finite differencing to discretize $(\mathcal{L} V^{n+1})_{\ell}$, the fully discrete system can be written as:

where $\alpha_{\ell}$ and $\beta_{\ell}$ are given by

for $\ell = 2, \ldots, m-1$. When $\ell = 1$, we set $\alpha_{\ell} = \beta_{\ell} = 0$ at $S_1 = 0$; and when $\ell = 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.

2.2.   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 $A = A_1-A_2$, where

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, $\rho$, is bounded by ^[19]:

Although the fixed point iteration is convergent, the convergence rate can be very slow. Note that $\lambda$ is a function of the density function $\nu(y)$. It can be estimated by

for $Y > 0$. As a result,

as $\Delta y$ approaches 0.

3.   Multigrid method

As shown in (2.11), the convergence of the fixed point iteration can be very slow when $Y > 0$ and $\Delta 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, $\hat{V}_h^k \approx V^{n+1}$, on the fine grid. The first step is pre-smoothing. A small number of the fixed point iterations (typically $\nu_1$ = 1 or 2) is applied to compute an approximate solution of (2.7),

where $\mathcal{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 \tilde{V}_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 = \mathcal{R} r_h$ is the restriction of $r_h$ given by a restriction operator $\mathcal{R}$. Here we use full weighting for $\mathcal{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 $\mathcal{P}$ is a prolongation operator which interpolates the coarse grid error $e_H$ onto the fine grid. Linear interpolation is used for $\mathcal{P}$; i.e.,

The last step is to apply smoothing again,

where $\nu_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.

3.1.   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 $\overline{V}(x, \tau) = V(e^x, \tau)$. Then (2.6) can be written as the following PIDE with constant coefficients,

where $\overline{V} \otimes \hat{\gamma}$ denotes a correlation product; i.e. $(\overline{V} \otimes \hat{\gamma})_i = \sum_{j = 0}^{N-1} \bar{V}(x_i, y_j) \hat{\gamma}(y_j)$. The PIDE is discretized similarly as described in Section 2.1 where the unknowns are now $\overline{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 $\tau^{n+1}$ is traced back to $x_{j^*}$ at time $\tau^n$. The value of $\overline{V}^n_{j^*} \equiv \overline{V}(x_{j^*}, \tau^n)$ is computed by an upwind quadratic interpolation. After rearranging terms, (3.3) becomes

where $\alpha = (\sigma^2 + \bar{\sigma})/(2 (\Delta x)^2)$. Then, the corresponding fixed point iteration (2.9) is given in component form by

where $\overline{V}^{n+1, k}$ is an approximate solution to $\overline{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 $\overline{V}_x$ in (3.2). Without loss of generality, assuming the coefficient is negative, we have

where $\psi$'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

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). Let $E^{k} = \overline{V}^{n+1, k}-\overline{V}^{n+1}$ be the error of $\overline{V}^{n+1, k}$. Let $C^k_{\mu}$ and $G_{\mu}$ be the symbols of $E^{k}$ and $\hat{\gamma}$, respectively, where $\mu$ is the Fourier mode; i.e. $C^k_{\mu}$ and $G_{\mu}$ are the Fourier transform of $E^k$ and $\hat{\gamma}$, respectively ^[20]. Also, let $h = 1/N$ be the mesh size and $i = \sqrt{-1}$. Then (3.7) becomes

Note the negative subscript of $G_{-\mu}$ 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 $\eta_j = \sum_{m = j-p-2}^{j-p} \psi_{j, m} e^{-(j-m)\mu \pi h i}$.

Let $\hat{S}_{\mu}$ be the symbol of the smoother. Then one can easily show that $\hat{S}_{\mu} = C_{\mu}^{k+1}/C_{\mu}^k$; i.e.,

It has been proved that $|\eta_j| \leq 1$ [40, Section 5.2]. Furthermore, since $G_{\mu}$ is the symbol of $\hat{\gamma}$, by (2.5), we have $|G_{-\mu}| \leq \lambda$. Hence, we can bound $\hat{S}_{\mu}$ as

Let $\zeta$ be the smoothing factor of the fixed point iteration, which is defined as

Considering the high frequency Fourier modes; i.e., $N/2 \leq \mu \leq N$, it follows from (3.11) that

To bring more insight into the formula, we estimate the upper bound for $\zeta$ by noting that

Considering the limit where $\Delta x \approx 0$ and keeping $\Delta \tau$ fixed, we have

In the Lévy process under the CGMY model, we always have $Y < 2$. Hence for small $\Delta x$, the smoothing factor for the fixed point iteration is small, which is fundamental to achieve an effective multigrid method.

3.2.   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, $\hat{\mathcal{M}}$, is a block diagonal matrix with each block a $2 \times 2$ matrix, $\hat{\mathcal{M}}_{\mu}$, which is given by

where

Here $c_{\mu} = \cos(\mu \pi h)/2$, $s_{\mu} = \sin(\mu \pi h)/2$, and $h = \Delta x$.

By (3.10), we have

where $\tilde{\sigma} = (\sigma^2+\bar{\sigma})/2$. Also, we have assumed the interest rate $r$ to be zero, just to simplify the algebraic expression.

A similar calculation yields

where $\eta_{j'}$ is defined similarly as $\eta_j$ on the coarse grid. Details can be found in Appendix A.

In order to prove convergence, we assume the following inequalities hold:

Note that $\tilde{\sigma} = O(1)$ and hence the right-hand side of assumption (3.17) is also $O(1)$. In practice, $\Delta \tau$ is typically chosen as $O(h)$. Thus assumption (3.17) can be easily satisfied. By (2.10), $\lambda h^2 = O(h^{2-Y})$. Since $Y < 2$, for sufficiently small $h$, $\lambda h^2$ will satisfy the bound in assumption (3.18).

Lemma 3.1. Let $\hat{\mathcal{M}}_{\mu}$ be defined as (3.12). Assume the inequalities (3.17) and (3.18) hold. Then

for any $\mu$.

Proof. The proof of the lemma is very technical. Here we will outline the main steps of the proof. We refer the interested reader to Appendix A for details.

The idea is to derive a bound for $|\hat{\mathcal{M}}_{\mu}(1, 1)|$ + $|\hat{\mathcal{M}}_{\mu}(1, 2)|$ and $|\hat{\mathcal{M}}_{\mu}(2, 1)|$ + $|\hat{\mathcal{M}}_{\mu}(2, 2)|$, where $\hat{\mathcal{M}}_{\mu}(i, j)$ is the $(i, j)$ entry of the $2 \times 2$ matrix $\hat{\mathcal{M}}_{\mu}$. Combining (3.13), (3.14), (3.15), and (3.15), we have

By direct computation and that facts that $|G_{\mu}| = |G_{-\mu}| \leq \lambda$ and $|\eta_j| \leq 1$, one can show that

By assumptions (3.17) and (3.18), one can further show that

where $\xi = 4\tilde{\sigma}/\lambda h^2$. The case for $|\hat{\mathcal{M}}_{\mu}(2, 1)| + |\hat{\mathcal{M}}_{\mu}(2, 2)| < 1$ is similar and hence is omitted here.

Theorem 3.2. Assume the inequalities (3.17) and (3.18) hold. Let $\rho(\mathcal{M})$ be the spectral radius of the iteration matrix $\mathcal{M}$. Then

Hence, the two-grid method converges and the upper bound is independent of the grid size and the value $Y$.

Proof. By lemma 3.1, $\|\hat{\mathcal{M}}_{\mu}\|_{\infty} < 1$ for any $\mu$. Furthermore, the proof of the lemma shows that

The last inequality comes from assumption (3.18). The same bound for $|\hat{\mathcal{M}}_{\mu}(2, 1)| + |\hat{\mathcal{M}}_{\mu}(2, 2)|$ can be derived similarly. Thus

As a result, $\rho(\mathcal{M}) = \rho(\hat{\mathcal{M}}) \leq \|\hat{\mathcal{M}}\|_{\infty} < 0.75$.

4.   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 \approx 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 $\Delta \tau = 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 $\Delta \tau$ 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 $\eta_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 of $\hat{\mathcal{M}}$. Instead, we compute an upper bound for the smoothing factor and the infinity norm of $\hat{\mathcal{M}}$ 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.

5.   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 \approx 0$) but very slowly when the kernel becomes more singular ($Y \approx 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.

Acknowledgments

This author was supported by the Natural Sciences and Engineering Research Council of Canada.

Conflict of interest

The author declares no conflict of interest.

Appendix

A.   Two-grid convergence analysis

Here we provide the technical details for the analysis in Section 3.2.

Lemma A.1. The Fourier symbols for the fine and coarse grid operators are given by

Proof. By (3.4), the fine grid operator can be written as

with the assumption that $r = 0$. Applying the local Fourier analysis similar to (3.8), we have

noting that $\alpha = \tilde{\sigma})/h^2$. The formula for $\hat{L}_{N-\mu}^h$ can be easily derived from $\hat{L}_{\mu}^h$ which is omitted here.

Similarly, the coarse grid operator can be written as

where $j'^*$ is defined similarly as $j^*$ on the coarse grid. Applying the local Fourier analysis, we obtain

Lemma A.2. Let $\hat{\mathcal{M}}_{\mu}$ be defined as (3.12). Assume $r = 0$ and the inequalities (3.17) and (3.18) hold. Then

for any $\mu$.

Proof. By (3.12), (3.13), (3.14), (3.15), and (3.15), we have

Thus, the (1, 1) and (1, 2) entries are given by

We first consider $\hat{\mathcal{M}}_{\mu}(1, 1)$. Using the formulas (3.15), (3.15), and (3.14) for $\hat{L}_{\mu}^h(1, 1)$, $\hat{L}_{\mu}^H$, and $\hat{\mathcal{S}}_{\mu}$, respectively, we obtain

where $\tilde{\sigma} = (\sigma^2 + \bar{\sigma})/2$. Since $G_{\mu}$ is the Fourier symbol of $\hat{\gamma}$, by (2.5), we have $|G_{\mu}| \leq \lambda$. Also, $|\eta_j| \leq 1$ as proved in [40], Section 5.2]. Hence

Write the complex number $G_{\mu} \eta_{j'}$ as $a + b \, i$ where $i = \sqrt{-1}$. Then we have

Note that $a^2 + b^2 \leq \lambda^2$ since $|G_{\mu} \eta_{j'}| \leq \lambda$. Thus $\lambda - a \geq 0$. Since all the terms on the right-hand side of (A.4) are non-negative, it can be lower bounded by

Combining (A.3) and (A.5), we have

Consider the first term on the right-hand side of (A.6). Since $0.5 \leq c_{\mu}^2 \leq 1$ and $4 \tilde{\sigma}s_{\mu}^2 \geq 0$, we have

By assumption (3.17), equation (A.7) becomes

Let $\xi = 4 \tilde{\sigma}/ \lambda h^2$. Then, we have

For the second term on the right-hand side of (A.6), we have

By assumption (3.18), we have

Also by assumption (3.18), we have

Multiplying (A.10) and (A.11) yields

which implies

Combining (A.9) and (A.12), we obtain

We will bound the third term on the right-hand side of (A.6) later. Now, consider $\hat{\mathcal{M}}_{\mu}(1, 2)$, the (1, 2) entry of $\hat{\mathcal{M}}_{\mu}$. By (A.2), (3.15), (3.15), and (3.14), we have

By a similar calculation as above, we can bound $|\hat{\mathcal{M}}_{\mu}(1, 2)|$ as:

The first term on the right-hand side can be bounded as:

where the last inequality follows from (A.8). The second term on the right-hand side of (A.14) can be bounded similarly as:

where the last inequality follows from (A.12). Combining (A.6)-(A.16), we obtain

By simple calculus, the maximum of the second term occurs when $\sin^2(\mu \pi h) = 1$. Hence it can be bounded by $1/(1+\frac{1}{2}\xi)$. As a result, we have

The case for the second row of $\hat{\mathcal{M}}_{\mu}$ is similar and we omit the corresponding calculations here.