Abstract
This paper deals with the numerical analysis of PIDE option pricing models with CGMY process using double discretization schemes. This approach assumes weaker hypotheses of the solution on the numerical boundary domain than other relevant papers. Positivity, stability, and consistency are studied. An explicit scheme is proposed after a suitable change of variables. Advantages of the proposed schemes are illustrated with appropriate examples.
1. Introduction
The hypothesis that asset prices behave according to the geometric Brownian motion when one derives the option prices is inconsistent with market prices [1]. This drawback has been overcome using Lévy process models [2–9] allowing the calibration of the model to the option market price and the reproduction of a wide variety of implied volatility skews/smiles; see [10] and [11, chapters 14, 15]. Among the Lévy process models, it is remarkable to distinguish these with finite activity, that is, jump diffusion models [2, 3], and those where the intensity of the jumps is not a finite measure [4–9]. These models are characterized by the fact that option price is given by the solution of a partial integrodifferential equation (PIDE) involving a second-order differential operator part and a nonlocal integral term that presents additional difficulties. In [12] wavelet methods are applied to infinite Lévy models. Monte Carlo approaches are developed by [13, 14]. Interesting analytic-numerical treatments have been introduced in [15–17]. The so-called COS method for pricing European options is presented in [15]. This is based on the knowledge of the characteristic function and its relation with the coefficients of the Fourier-cosine expansion of the density function. In [16], an expansion of the characteristic function of local volatility models with Lévy jumps is developed. The authors in [17] derive an analytical formula for the price of European options for any model including local volatility and Poisson jump process by using Malliavin calculus techniques.
Many authors used FD schemes for solving these PIDE problems [18–28]. Dealing with FD methods for such PIDEs, the following challenges should be addressed, for instance, how to approximate the integral term and how to localize a bounded computational domain in order to consider relevant information like large jumps. In addition, the possible singularities of the integral kernel should be carefully treated [18, 19].
The nonlocal character of the integral part involves a dense discretization matrix. In the outstanding paper [18], Cont and Voltchkova presented an explicit-implicit method (explicit into the integral part and implicit into the differential one) to obtain the numerical approximation of viscosity solutions for European and barrier options. An improvable issue of [18] is that, in order to approximate the truncated integral term, they assume a particular behavior of the solution outside of the bounded numerical domain. This last drawback is experienced by most of the authors; see [22, 26, 27].
Implicit FD methods for the numerical solution of the CGMY model have been used by Wang et al. [19] who proposed an implicit time stepping method avoiding dense linear systems but involving the iteration methods drawbacks of the implicit methods such as ungranted positivity. They also assume that, for large enough values of , the solution behaves like Black-Scholes.
In [21], the authors use an unconditional ADI FD method and accelerate it using fast Fourier transform (FFT) for jump diffusion models with finite jump intensity. Tavella and Randall in [24] use an implicit time discretization and propose a stationary rapid convergent iterative method to solve the full matrix problem quoted above, but with poor numerical analysis. A generalization of their iterative method to price American options is proposed in [25].
One of the most relevant and versatile Lévy models is the one proposed by Carr et al., the so-called CGMY model [8], that belongs to the family of KoBoL models [9]. It is considered a prototype of the general class of models with jumps and enjoys widespread applicability. The CGMY model allows diffusions and jumps of both finite and infinite activity. The CGMY Lévy density is given by where , , , and . The parameter allows controling the fine structure of asset return distribution. For , the Lévy process is of finite activity; that is, the measure is finite; . For , it is of infinity activity but finite variance; that is, . Finally, for , both the activity and variation are infinite. Note that, for , one gets the well known Variance Gamma process proposed by Madan and Seneta [29] as a particular case. So CGMY model is an extension of the Variance Gamma model [7].
The authors in [22] use FD methods discretizing the equation in space by the collocation method and using explicit difference backward schemes focused on the case of infinite activity and finite variation.
Very recently [28] proposed an efficient three-time-level finite difference method for the infinite activity Lévy model. Second-order convergence rate is shown in numerical experiments although the numerical analysis of the method is not developed.
The aim of this paper is the construction and numerical analysis of an explicit finite difference scheme of the PIDE for the CGMY model with measure given by (1):
Here is the option price depending on the underlying asset , the time is the volatility parameter, and and are the risk-free interest and the continuous dividend paid by the asset, respectively. The payoff function for a vanilla call option is given by where is the strike price.
Like [20, 23] for jump diffusion models we transform the original PIDE problem in order to remove the convection term to avoid possible numerical oscillations. With respect to the singularity of the integral kernel quoted above, the jump component in the neighborhood of log jump size zero is approximated by using a Taylor expansion, like [18, 19].
The selection of the boundary conditions of the numerical domain, the discretization of the infinite domain of the integration, and matching the discretization of both the differential and the integral part are important challenges. Some authors, like those of [18], assume a particular behavior of the solution outside of the bounded numerical domain. In order to weaken these hypotheses we do not truncate the infinite integral and we use a nonuniform partition of the complete unbounded domain, allowing a proper matching of the discretizations of the differential and integral parts by assuming asymptotic linear behavior of the solution. This strategy involves a double discretization with two spatial stepsize parameters that will allow a better flexibility to improve the approximation in different zones of the domain.
This paper is organized as follows. In Section 2, the integral part of (2) is approximated in a neighborhood of to obtain a new PIDE integral part extended outside a neighborhood of . Then a variable transformation is developed in order to remove both the convection and reaction terms of the differential part. Following the idea developed in [20] for the finite activity case, a double discretization explicit scheme for solving PIDE problem (2) is constructed in Section 3 for the infinite activity case affording the new challenges. Positivity and stability of the numerical solutions given by the scheme are studied in Section 4. Consistency of the scheme is treated in Section 5. In Section 6, some illustrative numerical examples show the advantages of the new discretization approach showing how the double discretization allows flexible improvement of the accuracy in different zones of the domain.
If is a vector in , we denote its infinite norm . Vector is said to be nonnegative if for all ; then we denote . For a matrix in , we denote . Matrix is said to be nonnegative if for all , , and we denote .
The exponential integrals have a major role in evaluating important class of integrals. Let and be continuous (real or complex) variables; the exponential integral of order denoted by is given by [30]
2. Transformation of the PIDE Problem
We begin this section by removing the singularity of the kernel of the integral term of PIDE (2). Let and let us split the real line into two regions and . For the term in , taking Taylor expansion for about one gets Taking into account (1) the integral part of (2) can be written as where the integrals can be evaluated with high accuracy using the exponential integrals ([30, 31], chapter 7). Let us denote where denotes the gamma function and is the exponential integral. Considering the expression (1) and the exponential integral (5) for , one gets Notice that (10) holds for . For the particular case where , one gets while for , we have For the remaining integrals in (8), we have Hence, the problem (2) takes the following form: where .
In order to remove the convection and reaction terms from (15), let us introduce the following transformation of variables: Hence the problem (15) is approximated by the following form: where Finally in order to combine both discretizations of the differential and integral part, we use to change the integrand as follows: where . For evaluating the integrals in all the positive real line, let us introduce a parameter that separates into . The point can be chosen according to the criteria used by [18, 32, 33] to truncate the numerical domain. For instance, in [26] one takes and in [20] one takes . To evaluate the integrals related to , they are transformed to finite integrals by using the substitution consequently, obtaining integrals of the form where , . In particular if then . Hence, the problem (17) takes the form
3. Numerical Scheme Construction
In this section a difference scheme for the problem (21)-(22) is designed. For the time variable, given , let be the time-step discretization and , , with integer. With respect to the spatial variable and for an arbitrary fixed , we divide the interval into equal intervals with a spatial step , with , . Note that the unbounded domain is transformed into by the above quoted change . Thus a uniform distributed mesh partition of the interval of the form , , and is mapped into a nonuniform mesh partition of , , . Hence, we have Let us denote , , , whereand . With respect to the approximation of (19), note that for each we need to evaluate two integrals corresponding to and , denoted by , . Let be the biggest with such that and let be the first such that . Then the expression (19) for the point has the following form: Then we apply the trapezoidal rule for the integrals over and because of (20) and using the first mean value theorem for integrals [34, page 1063], the two remaining integrals are approximated by Let us denote Depending on the location of for each with , we approximate given by (26) in the following form.
Case 1 (). Note that in this case , and thus is approximated by . Also one has in the domain of the integral and is approximated by , taking into account (20) for . Thus
Case 2 (). As and , the approximation of becomes
Case 3 (). Here and the approximation of is given by Assuming that tends to zero at least linearly as tends to zero one has by (1) and (28). On the other hand, assuming linear behavior of the solution for large values of , the integrand of (20) , as . Thus, both the terms involving and do not appear in the expressions of (29)–(31). Taking into account (25)–(31) the resulting difference scheme for the PIDE problem (21) takes the form In order to obtain a complete difference scheme, we include the initial and boundary conditions. From (22), we have On the other hand, for a vanilla call option the boundary condition for is and by assuming the linear behavior of the solution for large values of the spatial variable, we have and thus and the null integral term approximation , for all time level . Thus from (32) for , one gets For the sake of convenience to study the stability, we now introduce the vector formulation of scheme (32)–(35). Let us denote the vector in as and let be a tridiagonal matrix in related to the differential part, defined by where Let be the matrix in related to the integral part whose entries for each fixed in are defined by where Hence scheme (32)–(35) can be written in the form
4. Positivity and Stability of the Numerical Solution
The price of contracts modelled by PIDE must be a nonnegative value. Our objective here is to demonstrate that the solution of scheme (32)–(35) is conditionally nonnegative and stable.
First we study the positivity of the matrix . The following lemma has been proved in [20].
Lemma 1. With previous notation, assume that stepsizes , in and , in , satisfy (C1); (C2).
Then matrix given by (37) is nonnegative.
Note that as the matrix defined by (40)-(39) is always nonnegative, from Lemma 1 and (43) starting from nonnegative initial vector , the following result is established.
Theorem 2. With the hypotheses and notation of Lemma 1, the solution of scheme (32)–(35) is nonnegative if the initial values , .
The next result will be used below to guarantee stability.
Lemma 3. Let matrices and be defined by (37)–(39), and let ; then the following results hold. (1)Under conditions (C1) and (C2) of Lemma 1, . (2), where is defined by (14).
By [20, lemma 2], part 1 is proved. Since the norm of is given by if denotes the row containing the maximum of (44), one gets the elements of the summation in (45) are given by (40)-(39). To upper bound (45), we apply the change of variables in (8), resulting in which coincides with (19) when . Hence from (19), (45), and (46), we conclude that is an approximation for . Thus, for small enough and , one gets [35] Hence independently of the value of the size of matrix .
There are many definitions of stability in the literature; we recall our concept of stability in the next definition.
Definition 4. Let be a numerical solution of the PIDE (21), (22) computed from scheme (32)–(35) with stepsizes in , in , and in . Let be the corresponding vector form; that is, ; of (43). One says that is strongly uniformly stable if where is independent of , , , and .
If the property (49) is satisfied for appropriate relationships between the stepsizes , and , then one says that the strong uniform stability is conditional.
Theorem 5. With the previous notation, the numerical solution of scheme (32)–(35) is strongly uniformly stable if one satisfies the condition together with
Proof. Note that scheme (32)–(35) is equivalent to the vector form scheme (43). Under condition (50), by Lemma 3 one gets, after taking norms in (43), Hence, from (51) and that , , Thus the conditional strong uniform stability is established.
5. Consistency
A numerical scheme is consistent with a PIDE if an exact theoretical solution of the PIDE approximates well the difference scheme as the stepsizes discretization tend to zero [36, 37].
After recalling the definition of consistency, let us write (32) in the form Let us denote as the value of the theoretical solution of (21) and let such that , and let us write the PIDE (21) as where The local truncated error at is defined by In order to prove the consistency, we must show that Assuming that is twice continuously partially differentiable with respect to and four times partially differentiable with respect to , and using Taylor’s expansion about , it follows that In accordance with [36, page 101] let us denote the local consistency error of (see (29)) by By (26) and (29), the local consistency error for is given by where From the first mean value theorem for integrals [34, page 1063], one gets and since it follows that where Analogously, Let , , and be defined as where the second derivatives appearing in (67) are taken with respect to the variable for and , and with respect to the variable for . From the expression of the error of the trapezoidal rule, [35, page 54], (59)–(67), one gets Summarizing, one gets Thus, showing the consistency of the scheme with PIDE.
6. Numerical Results
In the following examples the Matlab program has been used. The first example reveals the effect of Yor parameter on the option price.
Example 6. Consider the vanilla call option problem (2)–(4) under CGMY process with parameters , , , , , , , , , , , and . Figure 1 exhibits the variation of the option price versus the underlying asset at various values of Yor parameter.
The next example illustrates the importance of positivity conditions given by Lemma 1.
Example 7. Here in this example the parameters have been selected as follows: , , , , , , , , , , , , and . Positivity conditions hold for , while for , the positivity conditions are broken and the values of the option price become unreliable as shown in Figure 2.
Next example shows that the variation of the absolute and relative error of the solution in light of the stability and positivity conditions hold at the strike for two cases: first, for several values of the stepsize discretization and second, for different values of the parameter .
Example 8. Consider the European call option for CGMY process with the following values: , , , , , , , , and , for several values of Yor parameter , and .
We consider the evaluation of the price option at the strike and . Table 1 reveals the deviation between our numerical solutions and the reference values used in [15, tables 8–10] for different stepsizes and fixed . Notice that the numerical solution exhibits the expected second-order convergence rate.
Table 2 shows the deviation for several values of , while .
In the next example, we consider the Variance Gamma model as a particular case () of CGMY model for which the exact solution is known [38].
Example 9. Consider a call option under the Variance Gamma process with parameters , , , , , , , , , , and . Figure 3 displays the associated error of the numerical solution for several values of the stepsize .
The next example shows that the double discretization strategy reduces the error near the parameter by changing the stepsize .
Example 10. Consider Example 9 with fixed ; Figure 4 shows the variation of the error of the numerical solution for various values of . Notice that the error decreases near the right boundary of the numerical domain by decreasing the stepsize , while the error near the strike remains stationary.
Acknowledgments
This work has been partially supported by the European Union in the FP7-PEOPLE-2012-ITN program under Grant Agreement no. 304617 (FP7 Marie Curie Action, Project Multi-ITN STRIKE-Novel Methods in Computational Finance) and by the Spanish M.E.Y.C. Grant DPI2010-20891-C02-01.