Abstract
Laplace transform method (LTM) has a lot of applications in the evaluation of European-style options and exotic options without early exercise features. However the Laplace transform methods for pricing American options have unsatisfactory accuracy and suffer from the instability. The aim of this paper is to develop a Laplace transform method for pricing American Strangles options with the underlying asset price following the constant elasticity volatility (CEV) models. By approximating the free boundaries, the Laplace transform is taken on a fixed space region to replace the moving boundaries space. After solving the linear system in Laplace space, Gaver-Stehfest formula (GSF) and hyperbola contour integral method (HCIM) are applied to compute the Laplace inversion. Numerical results show that the LTM-HCIM outperform the LTM-GSF in regard to the accuracy and stability for the option values.
1. Introduction
The Laplace transform methods for option pricing originate from the idea of randomizing the maturity in [1]. The Laplace transform methods are applied in the pricing options without early exercise features: Pelsser [2] for pricing double barrier options, Davydov and Linetsky [3] for pricing and hedging path dependent options under constant elasticity of variance (CEV) models, Sepp [4] for pricing double barrier options under double-exponential jump diffusion models, Cai and Kou [5] for pricing European options under mixed-exponential jump diffusion models, and Cai and Kou [6] for pricing Asian options under hyperexponential jump diffusion models.
This technique is known to work well for options without early exercise features, but, for American-style options, one difficulty has been perceived that the Black-Scholes-Merton PDE only holds where it is optimal to retain the option. Mallier and Alobaidi [7] develop a partial Laplace transform method for pricing American options in which the location of the free boundary for all values of the transform variable has to be determined by solving nonlinear integral equations. The approach is only conceptually appealing but practically ineffective in numerical implementations as for each transform variable one has to solve the nonlinear integral equations. In fact Mallier and Alobaidi [7] do not give the numerical computations.
A simple framework for Laplace transform methods is introduced by Zhu [8] for pricing American options under GBM models. This framework is later developed for evaluation of finite-lived Russian options by Kimura [9], pricing American options under CEV models by Wong and Zhao [10] and by Pun and Wong [11], pricing American options under hyperexponential jump-diffusion model by Leippold and Vasiljević [12], pricing stock loan (essentially an American option with time-dependent strike) by Lu and Putri [13], and pricing American options with regime switching by Ma et al. [14]. But this framework suffers from a drawback that numerical Laplace inversion such as Gaver-Stehfest (GS) and Gaver-Wynn-Rho (GWR) algorithm are not stable. Pun and Wong [11] explain the complex calculation of special functions as one of the possible reasons causing instability.
Zhou et al. [15] develop a new Laplace transform method for solving the free-boundary fractional diffusion equations arising in the American option pricing. By approximating the free boundary, the Laplace transform is taken on a fixed space region to replace the moving boundaries space. Then the hyperbola contour integral method is exploited to restore the option values. The coefficient matrix has theoretically proven to be sectorial. Therefore, the highly accurate approximation and computational stability of the fast Laplace transform method are guaranteed.
In this paper, we develop the new Laplace transform method for solving American Strangles option pricing under CEV model. The Black-Scholes-Merton PDE of Strangles option is defined on the moving region , where is put side and is call side. The idea of the new method is modifying PDE into a new form on certain fixed region . Then performing the Laplace transform leads to a ODE which involves the inverse functions and of and , respectively. Using the finite difference method combined with the approximation of functions and , where the optimal parameters of the approximation are obtained by minimizing the prescribed residual error, we obtain the numerical solution of the ODE in the Laplace space. In the final step, both the inversion Laplace Gaver-Stehfest formula (GSF) and the hyperbola contour integral method (HCIM) are used to recover the option value and the free boundary.
Numerical examples show the inverse Laplace GSF is unstable if the number of discrete integral nodes is great than 20, so HCIM is an effective algorithm for computing Laplace inversion. However, the technics of spectral analysis in [15] cannot be applied over here, because the coefficients of PDEs arising from CEV Strangles option are not constant. This paper gives a convergence theorem by analysing the so-called Laplace transform iteration algorithm.
The remaining parts of this paper are arranged as follows: In Section 2, we develop the Laplace transform method for pricing American CEV Strangles option; in Section 3, we discuss the hyperbola contour integral method to compute inverse Laplace and give some stable conditions for HCIM. In Section 4, some examples are taken to conform the theoretical results of HCIM. Conclusions are given in the final section.
2. Evaluation of American Strangles under CEV Model
Assume that the underlying asset price is governed by CEV (see, e.g., [16]):where is the risk-free interest rate, is the dividend yield, and is standard Brownian motion. represents the local volatility function and can be interpreted as the elasticity of . is the scale parameter fixing the initial instantaneous volatility at time , i.e., . If , the SDE (1) becomes the standard log-normal diffusion model or so-called geometric Brownian motion models (GBM). If , the SDE (1) nests the Cox-Ingersoll-Ross (CIR) model.
Let be the price of an American Strangles position written on an underlying asset with price at time and the payoffThis position is formed using a long put with strike and a long call with strike . Note that . Let , , and the early exercise boundary on the put side be denoted by , and the early exercise boundary on the call side be denoted by . It is known that satisfies the Black-Scholes PDE:with initial conditionand free boundary conditionsChiarella and Ziogas [17] apply Fourier transform technique to derive a coupled integral equation system for the Strangles free boundaries, and then a numerical algorithm is provided to solve this system.
To establish the new Laplace transform method proposed by Zhou et al. [15], we first note some properties of early exercise boundary for American Strangle option under the CEV model (1). It is known from [18–20] that the free-boundary functions and of American Strangles are continuously differentiable on the interval . We call the fact that is a decreasing function of variable , while is an increasing function of variable . DenoteThen we know, and satisfy the following ODE (similar to perpetual American option):
The values of and could be computed by secant method. Define uniform mesh and take larger enough number such that . The discrete form of (10) is as follows:for with . Assuming and , we have the pseudocode as Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||
Theorem 1. With the same initial condition (5) and boundary conditions (6)-(7), the PDE (4) is equivalent to the following PDE:for and , where is the indicator function.
Proof. (i) On the region , , we know from the boundary condition (6) that . Direct calculation shows that satisfies (14). (ii) On the region , , both (4) and (14) have the same forms. (iii) On the region , , we know from the boundary condition (7) that . Direct calculation shows that satisfies (14). Combining (i), (ii), and (iii), we get this theorem.
Let and be the inverse functions of and , respectively. Taking Laplace transformandto the both sides of (14), we geton fixed space region . The initial free boundary conditions (6) and (7) become the following forms:Now, the free boundaries and can be approximated. To this end, we give the fit functions with unknown parameters , for :Redefine uniform mesh with . Equations (18) can be discretized as follows:with index representing the number of space mesh partition. Matrix , solution vector , and RHS vector are defined aswhereDenote and to be the first row and the row of , respectively. and can be obtained by solving linear equation:Then, we reach the expressionsWe find the parameters , () such that conditions (19) and (20) hold for all real values of . A practicable way is to consider a positive sequence which are used for numerical Laplace inversion and find the values of unknown and that optimizing the following objective function:withAfter determining and we can compute by solving linear system (23), and the value of American Strangles can be expressed in terms of the inverse Laplace transformsApplying polynomial interpolations, we can obtain the option value for any values . Moreover, the vector version of inverse Laplace can be formulated as
One of the classical numerical Laplace inversion is the Gaver-Stehfest formula (GSF) (see [21, 22]):where with being taken as an even positive integer. Another effective numerical scheme for inverse Laplace transform is the so-called hyperbolic contour integral method (HCIM). We discuss HCIM in next section and give some examples, in Section 4, to illustrate the advantages of HCIM compared with GSF.
We finally summarize the Laplace transform method (NLTM) for pricing American Strangles as Algorithm 2.
3. Hyperbolic Contour Integral Method for Laplace Inversion
Many computational experiences show that the Gaver-Stehfest inverse Laplace formula (39) is unstable when (see [21, 22]; also see Example 1 in this paper). So, it is very important to develop an appropriate numerical algorithm for restoring the option values. In this section, we discuss how to apply the Hyperbolic contour integral method to compute from .
The inversion of Laplace transform is based on numerical integration of the Bromwich complex contour integral:where and is the convergence abscissa. This means that all the singularities of (with respect to ) lie in the open half-plane . The integral (41) is not well-suited for numerical integration. First, the exponential factor is highly oscillatory on the Bromwich line, , . Second, the transform typically decays slowly as . One strategy for circumventing the slow decay is due to Talbot [23], who suggests that the Bromwich line be deformed into a contour that begins and ends in the left half-plane, such that at each end. On such a contour, the exponential factor in (41) forces a rapid decay of the integrand as . This makes the integral well suited for approximation by the trapezoidal or midpoint rules. Owing to Cauchy’s theorem, such a deformation of contour is permissible as long as no singularities are traversed in the process, provided that uniformly in as .
Because the coefficient matrix in (23) is not Toeplitz matrix, the technics of spectral analysis in [15] cannot be applied over here. To apply the idea of Talbot [23] and Zhou [15], we rewrite the linear system (23) as follows:where is the identity matrix,and is defined by (27). Note that is a Toeplitz matrix. The linear system (42) can be solved by the following iteration algorithm:with .
Let -inner product and the corresponding vector norm be defined bywith being defined by (44). Matrix -norm is induced by the vector -norm, i.e.,
The theorem below shows that the iteration algorithm (46) is convergent if falls outside of a sectorial region.
Theorem 2. (i) The spectrum of defined by (43) lies in the sectorial region for any values of , i.e.,(ii) For defined by (43), we havewith being defined by (49) and .
(iii) Let with being defined by (45) and assume such that . Assuming solve (42), we haveTherefore, is convergent to according to -norm.
(iv) Under the assumptions of (iii), we have
Proof. (i) Matrix defined by (44) is a Toeplitz matrix and the generating function of isFrom the above generating function and the discussion of Pang and Sun [24], we know the spectrum lies the sectorial region for any values of . Since is positive definite, the spectrum also lies in the same sectorial region (see [24, Lemma 3.5]).
(ii) The inequality (50) is proven in [24, Theorem 3.3].
(iii) From linear system (42) and iteration algorithm (46), we have Therefore,andwhich means is convergent to according to -norm.
(iv) From linear equation (42) we have Consequently, which ends the proof.
From Theorem 2, we see the iteration sequence and the limit are analytical if lies out the sectorial region defined by (49). Therefore, the convergence of iteration algorithm (46) requires that is far away from the sectorial region such that . This condition could be satisfied in next discussion.
The Laplace inversion of with respect to can be expressed asThis is the vector version of Laplace inversion (41). Weideman et al. [25] (also see Pang and Sun [24]) suggest to select the hyperbolic contour , which can be parameterized bywhere parameters and set the width and the asymptotic angle of the hyperbolic contour, respectively. Then substituting the contour (60) into (59) givesDiscretization of this integral with uniform node spacing yieldswhere for the trapezoidal rule and for the midpoint rule, is the number of quadrature nodes, are the discretization errors, and is the truncation error. Defineand assume the function to be analytic in the strip with and , then the discretization errors are bounded aswhereand is some vector norm. The truncation error can be approximated by the magnitude of the last term retained, i.e.,We note that the estimates in (66) and (68) are extended in a straightforward manner, from the case being scalar functions (see [24–26]) to the case when takes value in a complex vector space.
To make all the singularities of the integrand in (61) to fall into the sectorial region defined by (49) for a given , the expressions of optimal parameters , , and must be chosen appropriately. By asymptotically matching and , and are derived as follows:and is a free parameter satisfying . With this choice the predicted convergence rate is given bywhereAccording to (70), we see that once the semiangle of the sectorial region (49) is given, then one can find the optimal by maximizing the function . Consequently, the optimal parameters and can be obtained from formulas (69).
Denoting and , we have
Now, we give the convergence condition of the iteration algorithm (46). From [27, Theorem 3.3], the distance between and is given by Inserting the above equation and the expression into following inequality we get Denote . If we take such thatthen and the iteration algorithm (46) is convergent according to norm . Inequality (76) requires that is large enough such that the hyperbolic contour is far away from the sectorial region .
In following Theorem, we give an estimation for and then give the minimum estimation for .
Theorem 3. For matrix being defined by (45), we have the following estimation:where and are defined by (8) in Section 2 andIf is small enough, then
Proof. Firstly, we verified the inequalitywhere represents the vector -norm. Indeed,From definition (45), we have withand . By carefully calculating, we havewhereand for , and . By the Gershgorin circle theorem, each eigenvalue of falls into a certain Gerschgorin disk , i.e.,for some index . Incorporating (84) and (85), we getTherefore,Incorporating (80) and (87), we haveLet , and ; we finally obtain the estimate expression (77).
Finally, we give two remarks as follows.
Remark 4. (i) For the fixed hyperbolic contour parameters and , Theorem 3 and iteration convergence condition (76) show that we must take <