Abstract
It is well known that linear complementarity problem (LCP) involving partial integro differential equation (PIDE) arises from pricing American options under Lévy Models. In the case of infinite activity process, the integral part of the PIDE has a singularity, which is generally approximated by a small Brownian component plus a compound Poisson process, in the neighborhood of origin. The PIDE can be reformulated as a fractional partial differential equation (FPDE) under fractional diffusion models, including FMLS (finite moment log stable), CGMY (Carr-Madan-Geman-Yor), and KoBol (Koponen-Boyarchenko-Levendorskii). In this paper, we first present a stable iterative algorithm, which is based on the fractional difference approach and penalty method, to avoid the singularity problem and obtain numerical approximations of first-order accuracy. Then, on the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Numerical tests are performed to demonstrate the effectiveness of the algorithm and the extrapolation method. We believe that this can be used as necessary tools by the engineers in research.
1. Introduction
The holder of American options has the right to exercise at any date prior to maturity and not only at the expiry date, while the holder of European options can exercise his right only at the expiry date. Thus, American options give the holder more freedom; hence, the price of American put is always greater or equal to the price of its European counterpart. There are two types of American options: put option and call option. The put option is a contract that the owner has the right but not the obligation to sell, while the holder of a call option has the right but not the obligation to buy.
Black and Scholes [1] and Merton [2] proposed the classic option pricing model called Black-Scholes model in 1973. In order to fit observed empirical market data better, many modifications to the Black-Scholes model have been proposed to build new models like jump-diffusion models proposed by Kou in [3] and Merton in [4] and infinite jump activity models such as KoBol [5, 6] and CGMY [7]. More general jump-diffusion models with stochastic volatility are considered in [8].
Compared with the jump-diffusion models, infinite activity models like KoBol and CGMY [7, 9] are more flexible to characterize the price process and give a more realistic description of the price process at various time scales. Both time-series and option data indicate that market indices lack a diffusion component and confirm the value of the infinite activity Lévy model [7].
Due to the optimal stopping problem, pricing American options lead to a a linear complementarity problem (LCP) or a variational inequality, which involves a partial integrodifferential equation (PIDE) in the case of models with jumps. Merton et al. [10] proposed a finite difference method for this type of formulation in 1977 under Black-Scholes model. More solutions to the LCP, such as the operator splitting method (OS) [11] under the Black-Scholes model, the PSOR method under the Black-Scholes model, the penalty method [12] under the jump-diffusion models, and the Lagrange multiplier method [13] under the Black-Scholes model, have been studied as well. Generally second-order accurate finite differences are used to discretize the PDE of the LCP for option pricing. In [14–16], fourth-order accurate finite difference methods have been considered for American options.
As we know, pricing American options are a free boundary problem [6]; hence, the free boundary and the price of an option must be solved simultaneously. If the models satisfy the smooth pasting principle, we can use two types of methods to find the moving location of free boundary: front tracking [17–19] and front fixing [20–22]. However smooth pasting principle does not always hold under models with jumps unless the models satisfy certain conditions.
In the case of pricing American options under the infinite activity jump models, how to deal with the integral part of the PIDE is important because it has a singularity in the neighborhood of origin. Almendral [23] approximates the pure-jump process with a compound Poisson process plus a small Brownian component and accelerates the computation by FFT. An integration by parts technique for rewriting the integrodifferential operator in terms of Volterra operators with a weakly singular kernel is proposed. These methods achieve second-order accuracy.
Fractional derivatives provide a more delicate instrument to capture the characteristics of processes or materials in comparison with the integer-order derivatives, which proves to be very useful in many fields [24]. Under certain conditions, the Grünwald-Letnikov definition of fractional derivative is equivalent to the Riemann-Liouville definition of fractional derivative; hence, the Grünwald-Letnikov definition, which is convenient for numerical valuation, is often used for discretization of the fractional derivatives. Due to the instability caused by standard Grünwald-Letnikov estimates, a shifted Grünwald-Letnikov definition, which is unconditionally stable, is proposed [25].
Cartea and del-Castillo-Negrete [26] presented the fractional partial differential equation (FPDE) under some infinite activity Lévy models (FMLS, KoBol, and CGMY) as a generalization to PIDE and gave a finite difference method using the numerical valuation of fractional derivatives. The numerical solutions of three fractional partial differential equations have been compared by Marom and Momoniat [27]. Their methods are of first-order accuracy and are only for evaluation of European options. Second-order accurate approximation schemes for fractional derivatives have been proposed in [28, 29]. The approach in [28] is based on the classical Crank-Nicholson method combined with spatial extrapolation to obtain temporally and spatially second-order accurate numerical estimates for fractional diffusion equation. By using proper shifting parameter, the generalized shifted Grünwald-Letnikov definition has been proven to be second-order accurate [29].
In this paper, our study focuses on the discretization of fractional derivatives, spatial extrapolation, and design of stable iterative algorithms for pricing American options in the scheme of FPDE. We first present a stable iterative algorithm, which is based on the fractional difference approach and penalty method, to avoid the singularity problem and obtain numerical approximations of first-order accuracy. On the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Discretization of Fractional derivatives results in dense matrices which make the computation expensive, so the fast Fourier transform (FFT) method can be employed to decrease the computation cost to gain an almost linear complexity from using direct solution.
The outline of this paper is as follows. In Section 2, linear complementarity problem and FPDE arising from pricing American options are introduced. In Section 3, we present the finite difference time and space discretization method. Section 4 describes our method used to solve the systems of linear equations and linear complementarity problems for American options. Numerical experiments are given in Section 5, and finally Section 6 contains conclusions. The main contributions of this paper are the discretization method in Section 3, the algorithm proposed in Section 4, and the numerical experiments in Section 5.
2. Mathematical Model for American Options
Because the value of an American call on an asset paying no dividends is equal to the value of the European call with same strike and maturity [2], only American put options are considered here.
To introduce the mathematical model, let us first recall the Riemann-Liouville definition [26] which is the most widely known definition of the fractional derivatives: where is an arbitrary order, is a positive integer, and is the boundary.
Under certain conditions, Grünwald-Letnikov definition is equivalent to the Riemann-Liouville definition; hence, we can use Grünwald-Letnikov definition to evaluate the Riemann-Liouville derivatives of the LCP for pricing American options. Grünwald-Letnikov definition [24] is more convenient for numerical valuation than the Riemann-Liouville definition.
Using the known definition of the binomial coefficients we have the Grünwald-Letnikov definitions of fractional derivatives where is an arbitrary order and is the boundary.
It has been proved that the standard Grünwald-Letnikov definition has stability problem if being used to approximate fractional derivate in the space-fractional advection-dispersion equation. Consider the following: So the shifted Grünwald-Letnikov definition was proposed to approximate Liouville fractional derivative : where is a nonnegative integer. is of first-order accuracy.
Generally, is set to 1. We have shifted Grünwald-Letnikov definitions for the left and right fractional derivatives: where is an arbitrary order and is the boundary.
In the following, a LCP involving FPDE for pricing American put options under fractional diffusion model is given. We denote the value of American put option by , the exercise price by , the expiry time by , the risk free interest rate by , and the underlying asset price by . The solution of the following LCP with additional constraints gives the price of the option.
Using the fractional partial differential equation (FPDE) given in [26] for pricing European options we have the following LCP for pricing American options: where and is the fractional order.
For KoBol model, we have where
For CGMY model, we have where
The final value is given by , where is the payoff function for the option contract. For a put option, it is . We should note that the function is the free boundary which divides into two parts. If , then we have the equation , which means that the price of the option is the exercise price. If , then equation means that the price is the solution of the above FPDE.
3. Space and Time Discretization
The above mathematical model involving the fractional derivatives assumes an infinite domain, . In order to apply finite difference methods, we will need to discretize and truncate the domain of . For and finite time interval , the terminal condition is given by . Further, the boundary conditions are and . By previous work in the field [11, 30], we take at four times the exercise price of the option to ensure that the errors generated by the truncation are small enough so as to be negligible. In fact, the choice of is simpler than that of because it is easy to choose proper so as to keep far from the moving free boundary, hence, the errors caused by the left truncation are small enough too.
For the space discretization, we use a uniform grid on interval with grid points. The grid step size is denoted by . Let . Then we denote space grid point by . Analogous to the space discretization, a uniform time step is used to get gird points on interval . Let . We denote time grid point by . In order to perform a finite difference approximation of the partial derivatives, we introduce a simplifying piece of notation along with the discrete approximation
3.1. Singularity of the Fractional Derivatives
Using the above truncation, we actually assume that and . As we know, is singular at the lower, , boundary and is singular at the upper, , boundary. To understand the nature of the singularity, we take the left fractional derivative as an example.
Expanding in Taylor series around , we have For , at least, the first term on the right-hand side is singular. Truncation of domain is a must for discretization; so singularity is an important problem that we will clarify in the following discussion on left and right fractional derivatives seperately.
The LCP (8) has a function which is the free boundary at time . Generally, the free boundary is far from the lower, , boundary. The price of American option is equal to , when , so we only need to deal with those fractional derivatives which satisfy the condition ; hence, no singularity problems exist for left fractional derivatives.
As to the right fractional derivatives, the singularity exists at the upper, , boundary if we try to discretize the derivative. Taking into account the boundary condition and the assumption , the equality holds in such a case.
3.2. First-Order Approximation for the Fractional Derivatives
Using the shifted Grünwald-Letnikov definition of fractional derivatives we have for which give the first-order accurate approximations.
3.3. First-Order Approximation for the Integer-Order Derivatives
There is a first-order space derivative in the LCP (8) to be approximated and we have many choices to use to obtain approximation of order 1 or the higher.
For first-order approximation, we have the forward difference or the backward difference
In this paper, we choose the backward difference method.
3.4. Crank-Nicolson Discretization
Assuming is the matrix arising from the finite difference scheme, a general time stepping scheme for pricing American options is where and is the free boundary at time step . For above we recover the explicit forward Euler scheme. gives the implicit backward Euler method. We have the Crank-Nicolson method when .
3.5. Matrix Formulation
Using (18), (20), and (21), we can now have the discretization of the full FPDE (7): where which can be efficiently computed by FFT.
Define matrices , , and such that Equation (22) then can be rewritten in matrix form as follows:
If is set to , this fractional Crank-Nicholson discretization approach will produce a method with a local truncation error, that is, .
3.6. Spatially Second-Order Approximation Obtained by Extrapolation
In numerical analysis, extrapolation is used to improve the rate of convergence of a sequence. Based on the spatially first-order accurate numerical approximation in 3.5, extrapolation can be employed to improve the low order of spatial convergence. First, we state below the result obtained in [28].
Proposition 1. Let , such that all derivatives of up to order belong to . For any integer , define the shifted Grünwald operator by
Then, one has for some constants independent of , , and the fact that
uniformly in .
Writing (27) as
where and are independent of , and (28) shows that the truncation error in the shifted Grünwald finite difference approximation is . Now, one can use the Richardson extrapolation method in the following way. First, the Crank-Nicholson method is applied for grid size and again for grid size . This gives two solutions of first-order accuracy for the spatial dimension. We denote by the solution of grid size and by the solution of gird size . Then, one has
where and are real numbers. Taking , it is obvious that is spatially 2nd-order accurate.
Based on the classical Crank-Nicholson method combined with spatial extrapolation, we obtain a solution with local truncation error .
The extrapolation method can be used at each time step before the maturity of options; because we can obtain two solutions of first-order estimates for diffenent spatial grid size at each time step, then second-order approximation can be computed at each time step.
4. The Penalty Method
In order to apply penalty method, we need to replace the following equation in LCP: with where is the payoff upon the exercise and in the limit as the penalty parameter , the solution satisfies . Let be the matrix given by the above penalty item. Consider the following: Then the matrix form for the penalty method is obtained using ((25), (31), (32)). Consider the following:
4.1. Iterative Algorithm
At each time step, we need to solve (33) to obtain the valuation of American options at grid points of . Option values at grid points of are given by payoff function . Iterative methods are introduced in Algorithm 1.
| ||||||||||||||||||
5. Numerical Results
In this section, we consider an American put option. A series of tests were carried out to study the effectiveness of our method under KoBol model.
In Table 1, column “change” shows the convergence of our approximations along with the temporal and spatial grid getting finer. Convergence rates of second-order estimates obtained by extrapolation are significantly improved, compared with those of first-order approximations. Hence, the test results confirm the effectiveness of our method.
Figure 1 is an example of first- and second-order approximations under same model conditions.
Figure 2 gives the evolution of pricing process for second-order approximation at all time steps.
6. Conclusions
In this paper, starting from the LCP involving FPDE, we first present the algorithms of first-order approximations for American options under fractional diffusion models. Then, on the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Numerical results in Section 6 prove that the algorithm is feasible and the extrapolation method is effective.
Using the Richardson extrapolation method, we even can obtain numerical approximations of higher order. Hence, the algorithm and extrapolation method form a useful tool box for the engineer. The use of fractional derivatives discretization conveniently avoids the singularity problem in the integral part of the PIDE. With the development of the fractional derivatives, we could see more possibilities of its application in finance.
Finally, as for American options, how to extend our method to other models and get combined with other existing methods might be interesting topics. In order to improve computation performance, new technology like parallel computing or new iterative algorithm having quicker convergence rate are important in further study.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was funded by Chinese NSF (Grant no. 91230109).