Abstract
This paper is concerned with the numerical solution of partial integrodifferential equation for option pricing models under a tempered stable process known as CGMY model. A double discretization finite difference scheme is used for the treatment of the unbounded nonlocal integral term. We also introduce in the scheme the Patankar-trick to guarantee unconditional nonnegative numerical solutions. Integration formula of open type is used in order to improve the accuracy of the approximation of the integral part. Stability and consistency are also studied. Illustrative examples are included.
1. Introduction
In Black-Scholes model, it is assumed that the probability distribution of the stock price is lognormal and the instantaneous log return is a geometric Brownian motion. However the market for the options shows that the geometric Brownian model for the underlying asset leads to underprice or overprice for these options [1]. Several models have been proposed to overcome these shortcomings such as models where the volatility follows a stochastic process like Heston model [2] and models incorporating jumps in the underlying asset following Lévy processes [3] and [4, chap. 14, 15]. There are two main features for the Lévy processes: first, models with finite activity, that is, jump diffusion models [5, 6], and, second, models with infinite intensity measure [7–9]. The numerical solution of these models can be achieved through three main techniques: by partial integrodifferential equation (PIDE) methods, as a solution of numerical integration, and by Monte Carlo simulation. Based on the close relation between the characteristic function of the probability density function and the Fourier-cosine expansion, the so-called COS method has been used to obtain the option pricing for European options [10].
In this paper we focus on PIDE approach. The PIDE valuating the option price presents a differential part with reaction, convection, and diffusion terms, while the nonlocal integral part is extended over an infinite or semi-infinite interval. Several finite difference (FD) schemes have been proposed to solve numerically these PIDE problems [11–20]. In order to implement FD methods, there are many challenges to face such as how to treat the unbounded domain for the spatial variable and the possible singularity of the kernel of the integral term. Furthermore, a way of numerical integration has to be chosen and the discretizations of both the integral part and the differential part have to be matched [16, 19].
The model proposed by Carr et al., the so-called CGMY model, is one of the most practical and adaptable Lévy models [9]. In the CGMY model, the diffusions and jumps can be of finite or infinite activity. Its Lévy density is given by where measures the overall level of activity, and measure the skewness, and controls 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 infinite activity but finite variance; that is, . Finally, for , both the activity and variation are infinite. Beside that for one gets the well known Variance Gamma process proposed by Madan and Milne [8] as a special case.
In the innovative paper [16], Cont and Voltchkova provided implicit discretization for the differential part and explicit step for the integral part after truncating its domain. The singularity of the integral kernel and the nonsmoothness of initial conditions are treated using the viscosity solutions and applying this technique to Merton and Variance Gamma models. It is worth mentioning that, in the Variance Gamma model, the integral part is split into singular and nonsingular parts.
An implicit FD method for CGMY model has been proposed by Wang et al. in [19] and a semi-Lagrangian discretization has been used for the drift term. On the other hand the integral part is split into singular and nonsingular parts; in the singular part the singularity has been removed using Taylor expansion; after that the trapezoidal rule is implemented. In order to adapt the unknown option price of the integral part to the computational mesh points, it is approximated using the upwind quadratic interpolation.
The option pricing for jump diffusion models with finite jump intensity has been treated using ADI finite difference method, accelerated by the fast Fourier transformation [13]. In [20], a three-time-level finite difference method is proposed showing a second order convergence rate in the numerical experiments for infinite activity models. However, the authors in [20] focus the interest on computational issues more than the numerical analysis.
An explicit scheme has been used in [15], applying the trapezoidal rule to approximate the integral term after removing the singularity of the kernel and including the unbounded domain using a double discretization technique [14]. Furthermore, the authors provided conditions to guarantee the positivity and stability of the numerical solution.
Dealing with option pricing models, ensuring positive solutions is a necessary requirement. Our objective is to construct a stable and conditionally consistent numerical scheme that guarantees positive solutions for the PIDE governing 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 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. Numerical methods that guarantee the positive of solutions for parabolic equations have been studied in [21, 22] following the idea initiated by Patankar, the so-called Patankar-trick [23].
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 [19]. Then a variable transformation is developed in order to remove both the convection and reaction terms of the differential part. Following partially the ideas about double discretization developed in [14, 15] an explicit scheme is constructed in Section 3. Unconditionally positivity and stability of the numerical solutions are shown 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. The paper ends with a conclusion section.
2. Removing the Reaction and Convection Terms of the PIDE Problem
This section begins with removing the singularity of the kernel of the integral part of PIDE (2). In order to achieve this aim, we split the real line into two regions depending on a parameter ; that is, and ; see [16]. By using Taylor expansion of for in about , one gets a convergency of order ; see [19].
On the other hand, the convection and reaction terms in (2) can be removed, obtaining a simpler PIDE with further numerical advantages, with the following transformation; see [15]: such that Under transformation (5) problem (2)-(3) is approximated to the following form: where In order to match both discretizations of the differential and integral part, let us consider the substitution in (10). Then, the resulting expression can be written in the more compact form where the new kernel takes the form In order to evaluate integral (11) without truncation, following the double discretization technique developed in [14, 15], let us introduce a parameter splitting into . The point can be chosen according to the criteria used by [16, 24, 25] to truncate or split the numerical domain. For instance, in [17] one takes and in [14] one takes . The unbounded integral part related to is transformed to a finite one by means of the change obtaining an integral of the form where and . In particular if then . Since there is a wide class of integrals that can be evaluated using exponential integrals, we recall the definition of the exponential integrals. Let and be continuous (real or complex) variables; the exponential integral of order denoted by is given by [26] From computational point of view, the integrals in (6) and (9) can be evaluated efficiently using exponential integrals [26] and [27, chapter 7], as follows: where denotes the gamma function and is the exponential integral (14). Notice that (17) holds for , while for particular cases and one can find the expression of in [15].
3. Computing the Numerical Solution
In this section a positivity-preserving explicit difference scheme for problem (7)–(13) is constructed. 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 , , is mapped into a nonuniform mesh partition of , , . Hence, we have Let us denote , , , where, in the approximation of the second partial derivative, we use Patankar-trick [23] obtaining two-time-level approximations: and . Then the difference scheme for (7) has the following norm: Note that the first expression of (21) corresponds to spatial zone with uniform discretization, while the second expression of (21) is related to the nonuniform discretization. On the other hand, for the approximation of the integral part of (7), instead of using the trapezoidal rule like in [14, 16, 20], we use a composite four-point integration formula of open type because of the higher order approximation of this rule [28, pp. 92-93]. This higher accuracy comes out because the singularity points of the kernel are not nodes of the integration mesh due to the truncation (see (12)) and the open type nature of the quadrature formula. Thus the approximation of (13) corresponding to the nodes and is given by where . Consequently the corresponding difference scheme for PIDE given by (7) takes the following form: where From (24)-(25), one gets where
In order to complete the difference scheme, we include the initial and boundary conditions. From (8), 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 levels . Thus from (26) for , one gets For the sake of convenience in the study of stability, we now introduce the vector formulation of the scheme (26)–(30) for the sake of the stability. Let us denote the vector in as , and let be the tridiagonal matrix related to the differential part and defined by Let be the matrix in related to the integral part whose entries for each fixed in are defined by , where Hence scheme (26)–(30) can be written in the form
4. Unconditional Positivity and Stability
The numerical solution of scheme (26) is unconditionally nonnegative because all coefficients of (26) and the initial and boundary conditions (28)–(30) are nonnegative.
For the sake of clarity, before studying stability, let us recall the definition of the infinite norm for vectors and matrices. For a given vector such that , the infinite norm of is denoted by and is defined as . For a matrix , the infinite norm of is denoted by and is defined to be .
In order to study the stability of the scheme given by (26)–(30), we first calculate the norm of the matrices and . Since the norm of the matrix is obtained by for , and for , it follows that . On the other hand, Let be the row containing the maximum of (35), as ; then we have Note that the change of variables in (6) gives which coincides with (10) when . Hence from (10), (35), and (37), we conclude that is an approximation for . Thus, for small enough and , one gets [28] Hence, from (35) and (37) independently of the size of matrix .
Since there are many definitions of stability, let us point out that the chosen concept of stability is the one used in [15, pag. 7]. Note that scheme (26)–(30) is equivalent to the vector form scheme (33). After taking norms in (33), Hence, from (40) and the fact that , , Thus the following result has been established.
Theorem 1. The numerical scheme given by (33) is strongly uniform stable.
5. Conditional Consistency of the Scheme
We say that a numerical scheme is consistent with a PIDE if the exact solution of the PIDE approximates well the difference scheme as the stepsizes discretization tend to zero [29, 30].
Let us write (22) in the form let be the value of the exact price of (7) at , let such that , and let us write the PIDE (7) as where is given by (11) and Now we show that the local truncated error at , given by satisfies 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 where Let us denote the maximum of the associated errors , , and by , , and , respectively, where To study the consistency of the integral part, it is convenient to rewrite it in the following form: where In accordance with [29] let us denote the local consistency error of : where By [28, p. 92] we have where in (56) is the fourth derivative with respect to the variable .
Similarly, the local consistency error for the unbounded region is given by where Thus the local truncation error is given by Consequently, the order of the local truncated error is given by and the following result has been established.
Theorem 2. Under condition , , the numerical scheme (42) is consistent with (7).
6. Numerical Examples
Based on the double discretization and Patankar-trick, a difference scheme has been established to obtain a numerical solution for the option price. This solution is guaranteed to be nonnegative and stable. In this section, we illustrate with several examples the behavior of the option price obtained by this scheme using MATLAB. All examples are done using CPU with Microprocessor 3.4 GHz Intel Core i7.
The following example illustrates that the consistency condition cannot be ignored.
Example 1. Here in this example the parameters have been selected as follows: , , , , , , , , , , , , and , for several values of such that , , and . Figure 1 shows that the consistency condition holds for , while, for the other two values, it is broken and the values of the option price become unreliable.
The aim of the following examples is to exhibit the effects of different parameters such as , , and on the variation of the absolute error in two cases: first, when (Variance Gamma case) and, second, for CGMY process when . Also the CPU time is given in seconds (sec.).
In the next example we calculate the associated error with this numerical scheme for the Variance Gamma model as a particular case () of CGMY model for which the exact solution is known [31].
Example 2. For parameters have been selected as follows: , , , , , , , , , , and ; Table 1 shows the variation of the absolute error with with fixed and for two values of and . From Table 1, it is observed that the associated error exhibits the second order convergence rate providing that is small enough in all the cases.
Table 2 reveals the change of the associated error for various values of time stepsize , while for and . Notice that the associated error due to the change of satisfies the expected first order convergence rate .
The aim of Table 3 is to show the sensitivity of the associated error of the option price due to the variation of , for and , while .
Example 3. Here we compare in Tables 4 and 5 in our results with the reference values given in [10, Tables 9, 10] related to accuracy and computational time. We consider the CGMY model for the following parameters: , , , , , , , , , and . Table 4 shows the variation of the associated error for several values of when and , while .
The variation of the associated error for several values of is presented in Table 5 for and , while .
The next example reveals that the double discretization strategy reduces the error near the parameter by changing the stepsize .
Example 4. Consider a call option under the Variance Gamma process with parameters , , , , , , , , , , and . Figure 2 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.
Example 5. Figures 3(a) and 3(b) describe the behavior of the Greek parameters Delta and Gamma for European call option. They exhibit the Greek parameters as functions in the underlying asset and time . The parameters have been chosen as follows: , , , , , , , , , , , , , and .
(a) as a function of and
(b) as a function of and
7. Conclusion
In this paper we develop a double discretization technique for the numerical solution of option pricing PIDE model under CGMY process, previously used in [15]. This technique allows us to consider the information of the solution of the integral part outside of the bounded numerical domain. In this paper we propose a finite difference scheme that introduces two main advantages with respect to [15] such as the unconditional positivity and stability of the proposed numerical scheme as well as a higher accuracy due to the use of more accurate numerical integration.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
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).