Laplace Transform Method for Pricing American CEV Strangles Option with Two Free Boundaries

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.


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 mixedexponential 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 2 Discrete Dynamics in Nature and Society 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.

Evaluation of American Strangles under CEV Model
Assume that the underlying asset price is governed by CEV (see, e.g., [16]): where > 0 is the risk-free interest rate, ≥ 0 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 = 0, i.e., = ( 0 )/ 0 fl 0 / 0 . If = 0, the SDE (1) becomes the standard log-normal diffusion model or so-called geometric Brownian motion models (GBM). If = −1/2, 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 payoff This position is formed using a long put with strike 1 and a long call with strike 2 . Note that 1 < 2 . Let = − , V( , ) = ( , − ), 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 V( , ) satisfies the Black-Scholes PDE: with initial condition and free boundary conditions Chiarella 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][19][20] that the freeboundary functions ( ) and ( ) of American Strangles are continuously differentiable on the interval [0, +∞). We call the fact that ( ) is a decreasing function of variable , while ( ) is an increasing function of variable . Denote Then we know Discrete Dynamics in Nature and Society 3 1: Step 1. Let 1 = 1, 1 = 1 , 2 = 2 , 2 = − 1.
Let ( ) and ( ) be the inverse functions of ( ) and ( ), respectively. Taking Laplace transform Discrete Dynamics in Nature and Society and to the both sides of (14), we get on 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 = 1, 2: Redefine uniform mesh = + Δ , = 0, 1, . . . , +1 with Δ = ( − )/( + 1). Equations (18) can be discretized as follows: with index representing the number of space mesh partition. Matrix A, solution vectork ( ), and RHS vector g( ) are defined as Denote a 1 ( ) and a ( ) to be the first row and the ℎ row of ( I − A) −1 ( ), respectively. a 1 ( ) and a ( ) can be obtained by solving linear equation: Then, we reach the expressionŝ We find the parameters , ( = 1, 2) such that conditions (19) and (20) hold for all real values of > 0. A practicable way is to consider a positive sequence { 1 , 2 , . . . , } which are used for numerical Laplace inversion and find the values Step 1. Use Algorithm 1 to compute and .
Step 4. Apply Laplace inversion GSF or HCIM to restore American Strangles option value k ( ).
of unknown and that optimizing the following objective function: After determining ( ) and ( ) we can computek ( ) by solving linear system (23), and the value of American Strangles can be expressed in terms of the inverse Laplace transforms Applying polynomial interpolations, we can obtain the option value V( , ) for any values ∈ (0, ∞). 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.

Hyperbolic Contour Integral Method for Laplace Inversion
Many computational experiences show that the Gaver-Stehfest inverse Laplace formula (39) is unstable when ≥ 20 (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 k ( ) fromk ( ).
The inversion of Laplace transform is based on numerical integration of the Bromwich complex contour integral: where i = √ −1 and 0 is the convergence abscissa. This means that all the singularities ofV( , ) (with respect to ) lie in the open half-plane Re( ) ≤ 0 . The integral (41) is not wellsuited for numerical integration. First, the exponential factor is highly oscillatory on the Bromwich line, = + i , −∞ < < +∞. Second, the transformV( , ) 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 Re( ) → −∞ at each end. On such a contour, the exponential factor in (41) forces a rapid decay of the integrand as Re( ) → −∞. 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 thatV( , ) → 0 uniformly in Re( ) ≤ 0 as | | → ∞.
Because the coefficient matrix A 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 I is the identity matrix, and g( ) is defined by (27). Note that B is a Toeplitz matrix.
The linear system (42) can be solved by the following iteration algorithm: with [k ( )] (0) = 0. Let D −1 -inner product and the corresponding vector norm be defined by with D being defined by (44). Matrix D −1 -norm is induced by the vector D −1 -norm, i.e., Discrete Dynamics in Nature and Society 7 The theorem below shows that the iteration algorithm (46) is convergent if falls outside of a sectorial region.

Proof. (i) Matrix B defined by (44) is a Toeplitz matrix and the generating function ( ) of B is
From the above generating function ( ) and the discussion of Pang and Sun [24], we know the spectrum Λ(B) lies the sectorial region Σ for any values of ∈ (0, /2). Since D is positive definite, the spectrum Λ(Ã) = Λ(DB) also lies in the same sectorial region Σ (see [24,Lemma 3.5
The Laplace inversion ofk ( ) with respect to can be expressed as This 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 by where parameters > 0 and set the width and the asymptotic angle of the hyperbolic contour, respectively. Then substituting the contour (60) into (59) gives 8 Discrete Dynamics in Nature and Society Discretization of this integral with uniform node spacing ℎ yields where = ℎ for the trapezoidal rule and = ( + 1/2)ℎ for the midpoint rule, is the number of quadrature nodes, ± (ℎ) are the discretization errors, and (ℎ ) is the truncation error. Define and assume the function ( + i ) to be analytic in the strip ∈ ( − , + ) with + > 0 and − < 0, then the discretization errors are bounded as where and ‖⋅‖ 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][25][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 > 0, 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 < /2 − . With this choice the predicted convergence rate is given by where According 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 = ((4 − 2 + 2 )/ ( ))( / ) into following inequality we get Denote = (1 − sin ) = dist(Σ , Γ). If we take such that then 1 + > 0 and the iteration algorithm (46) is convergent according to norm ‖ ⋅ ‖ D −1 . 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 ‖C‖ D −1 and then give the minimum estimation for . Theorem 3. For matrix C being defined by (45), we have the following estimation: Discrete Dynamics in Nature and Society 9 where and are defined by (8) in Section 2 and

Numerical Examples
In this subsection, we list some Strangles option values computed by the Laplace transform method and hyperbola contour integral method (labeled by "LTM-HCIM"). These results are compared with those obtained by LTM using inverse Laplace GSF (labeled by "LTM-GSF") and FDM. The parameters in model (4) Table 2 lists the space convergence rate with = 20 and fixed = 3000. The column labeled "Value" shows the option values at stock price 0 = 1.5; the column labeled "Err" is defined by where I[k 2 , ( )] is the linear interpolation of k 2 , ( ) from 2 -mesh to -mesh, ‖ ⋅ ‖ is the ℓ 2 vector norm and Δ ( ) = ( − )/( + 1), and the convergence rate is defined by From Table 2, we see the convergence rate of LTM-HCIM is 0.5 and 1.0 for = 0.5 and = −0.5, respectively. It is very interesting to investigate the relationship between the convergence rates and the values of and we leave this topic to the future. We regard the computational solution k , ( ) with = 1599, = 30 as the reference option prices, i.e., k ( ) ≈ for different values of . Figure 1 for moderate values of .

Conclusions
An efficient Laplace transform method was developed for pricing American Strangles option under CEV diffusion. Numerical examples show LTM is a fast algorithm to solve PDEs with free boundaries and the hyperbola contour integral method has higher numerical accuracy and good stability for numerical Laplace inversion. The method developed here can be extended to the pricing of other American options, such as American options with regime switching and with jump-diffusion process.

Data Availability
All Matlab codes supporting the results of this study are available from the corresponding author on request.
other sections. All authors read and approved the final manuscript.