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

Authors: Zhiqiang Zhou; Hongying Wu
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## 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.


---

## Body

## 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[Sp(τ),Sc(τ)], where Sp(τ) is put side and Sc(τ) is call side. The idea of the new method is modifying PDE into a new form on certain fixed region [S_p,S¯c]. Then performing the Laplace transform leads to a ODE which involves the inverse functions τp(S) and τc(S) of Sp(τ) and Sc(τ), respectively. Using the finite difference method combined with the approximation of functions τp(S) and τc(S), 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 Section2, 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]):(1)dSt=r-qStdt+δStβ+1dWt,where r>0 is the risk-free interest rate, q≥0 is the dividend yield, and Wt is standard Brownian motion. σ(S)=δSβ represents the local volatility function and β can be interpreted as the elasticity of σ(S). δ is the scale parameter fixing the initial instantaneous volatility at time t=0, i.e., δ=σ(S0)/S0β≔σ0/S0β. 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(2)fS,t=supτ∈T0,T⁡Ee-rτ-tΠSτ∣St=Sbe the price of an American Strangles position written on an underlying asset with price S at time t and the payoff(3)ΠSτ=max⁡K1-Sτ,0+max⁡Sτ-K2,0.This position is formed using a long put with strike K1 and a long call with strike K2. Note that K1<K2. Let τ=T-t, v(S,τ)=f(S,T-t), and the early exercise boundary on the put side be denoted by Sp(τ), and the early exercise boundary on the call side be denoted by Sc(τ). It is known that v(S,τ) satisfies the Black-Scholes PDE:(4)∂v∂τ=12δ2S2β+2∂2v∂S2+r-qS∂v∂S-rv,Spτ<S<Scτwith initial condition(5)vS,0=max⁡K1-S,0+max⁡S-K2,0,Sp0<S<Sc0and free boundary conditions(6)vS,τ=K1-S,∂∂SvS,τ=-1,S≥Spτ,(7)vS,τ=S-K2,∂∂SvS,τ=1,S≤Scτ.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–20] that the free-boundary functions Sp(τ) and Sc(τ) of American Strangles are continuously differentiable on the interval [0,+∞). We call the fact that Sp(τ) is a decreasing function of variable τ, while Sc(τ) is an increasing function of variable τ. Denote(8)S_p=Sp∞,S¯p=Sp0,S_c=Sc0,S¯c=Sc∞,Then we know(9)S¯p=min⁡K1,rqK1,S_c=max⁡K2,rqK2,S_p>0, S¯c<∞ and v∞(S)=v(S,∞) satisfy the following ODE (similar to perpetual American option):(10)12δ2S2β+2d2v∞dS2+r-qSdv∞dS-rv∞=0,S_p<S<S¯c,(11)v∞S_p=K1-S_p,ddSv∞S_p=-1,(12)v∞S¯c=S¯c-K2,ddSv∞S¯c=1.The values ofS_p and S¯c could be computed by secant method. Define uniform mesh Sj=jΔS,j=0,1,…,N and take larger enough number Smax=SN such that Smax>S¯c. The discrete form of (10) is as follows:(13)δ2Sj2β+22ΔS2-r-qSj2ΔSvj-1-δ2Sj2β+2ΔS2+rvj+δ2Sj2β+22ΔS2+r-qSj2ΔSvj+1=0.for j=1,2,…,N-1 with vj=v∞(Sj). Assuming S¯p=Sj1 and S_c=Sj2, we have the pseudocode as Algorithm 1.Algorithm 1:Algorithm for solving S_p and S¯c.
1:Step1. Let L1=1, H1=j1, L2=j2, H2=N-1. 2:Step2. Do the following loop: 3:while  L1<H1-1 or L2<H2-1  do 4:Let k1 = floor[0.5∗(L1+H1)] and k2 = floor[0.5∗(L2+H2)]. 5:Solve linear system (13) with boundary vk1=K1-Sk1 and vk2=Sk2-K2. 6:% Finite difference approximation (6). 7:if  vk1+1-vk1ΔS<-1  then 8:L1=k1; 9:else 10:H1=k1; 11:end if 12:% Finite difference approximation (7). 13:if  vk2-vk2-1ΔS<1  then 14:L2=k2; 15:else 16:H2=k2. 17:end if 18:end while 19:Step3. Output S_p=Sk1 and S¯c=Sk2.Theorem 1.
With the same initial condition (5) and boundary conditions (6)-(7), the PDE (4) is equivalent to the following PDE:(14)∂v∂τ=12δ2S2β+2∂2v∂S2+r-qS∂v∂S-rv+rK1-qS1S≤Spτ-rK2-qS1S≥Scτ.for τ∈(0,+∞) and S∈(S_p,S¯c), where 1{·} is the indicator function.Proof.
(i) On the regionτ∈(0,∞),  S∈(0,S_p(τ)), we know from the boundary condition (6) that v(S,τ)=K1-S. Direct calculation shows that v(S,τ)=K1-S satisfies (14). (ii) On the region τ∈(0,∞), S∈(S_p(τ),S¯c(τ)), both (4) and (14) have the same forms. (iii) On the region τ∈(0,∞), S∈(S¯c(τ),∞), we know from the boundary condition (7) that v(S,τ)=S-K2. Direct calculation shows that v(S,τ)=S-K2 satisfies (14). Combining (i), (ii), and (iii), we get this theorem.Letτp(S) and τc(S) be the inverse functions of Sp(τ) and Sc(τ), respectively. Taking Laplace transform(15)Lvτz≔∫0∞e-zτvτdτ,(16)L1S≤Spτ=∫0+∞e-zτ1S≤Spτdτ=∫0+∞e-zτ1τ≤τpSdτ=∫0τpSe-zτdτ=-1ze-zτ0τpS=1z1-e-zτpS,and(17)L1S≥Scτ=1z1-e-zτcS,to the both sides of (14), we get(18)zv^=12δ2S2β+2∂2v^∂S2+r-qS∂v^∂S-rv^+max⁡K1-S,0+max⁡S-K2,0+rK1-qS1z1-e-zτpS-rK2-qS1z1-e-zτcSon fixed space region S∈(S_p,S¯c). The initial free boundary conditions (6) and (7) become the following forms:(19)v^S_p,z=1zK1-S_p,∂∂Sv^S_p,z=-1z,(20)v^S¯c,z=1zS¯c-K2,∂∂Sv^S¯c,z=1z.Now, the free boundaries τp(S) and τc(S) can be approximated. To this end, we give the fit functions with unknown parameters pi, αi for i=1,2:(21)τpS≈τpS;p1,α1=0,S≥S¯p,p1-log⁡S-S_pS¯p-S_pα1,S_p<S<S¯p,+∞,S=S_p,(22)τcS≈τcS;p2,α2=+∞,S≥S¯c,p2-log⁡S¯c-SS¯c-S_cα2,S_c<S<S¯c,0,S=S_c.Redefine uniform mesh Sj=S_p+jΔS,j=0,1,…,n+1 with ΔS=(S¯c-S_p)/(n+1). Equations (18) can be discretized as follows:(23)zI-Av^nz=gz,with index n representing the number of space mesh partition. Matrix A, solution vector v^n(z), and RHS vector g(z) are defined as(24)Az=b1c10⋯00a2b2c2⋯000a3b3⋱00⋮⋱⋱⋱⋱⋮000⋱bn-1cn-1000⋯anbnn×n,(25)aj=δ2Sj2β+22ΔS2-r-qSj2ΔS,bj=-δ2Sj2β+2ΔS2-r,cj=δ2Sj2β+22ΔS2+r-qSj2ΔS,(26)v^nz=v^nS1,z,v^nS2,z,…,v^nSn,z⊺,(27)gz;p1,p2,α1,α2=g1,…,gn⊺,where(28)g1=a1K1-S_p,(29)gj=max⁡K1-Sj,0+max⁡Sj-K2,0+1zrK1-qSj1-e-zτpSj;p1,α1-1zrK2-qSj1-e-zτcSj;p2,α2,j=2,…,n-1(30)gn=cnS¯c-K2.Denote a1(z) and an(z) to be the first row and the nth row of (zI-A)-1(z), respectively. a1(z) and an(z) can be obtained by solving linear equation:(31)a1zzI-A=1,0,0,…,0,(32)anzzI-A=0,0,…,0,1.Then, we reach the expressions(33)v^S1,z=a1zgz;p1,p2,α1,α2,v^Sn,z=anzgz;p1,p2,α1,α2.We find the parameters pi, αi (i=1,2) such that conditions (19) and (20) hold for all real values of z>0. A practicable way is to consider a positive sequence {z1,z2,…,zm} which are used for numerical Laplace inversion and find the values of unknown pi and αi that optimizing the following objective function:(34)minp1,p2,α1,α2≥0⁡Fp1,p2,α1,α2=∑i=1mE1i;p1,p2,α1,α2+∑i=1mE2i;p1,p2,α1,α2with(35)E1i;p1,p2,α1,α2=v^S1,zi-v^S0,ziΔS-∂∂Sv^S0,zi2=a1zigzi;p1,p2,α1,α2-1/ziK1-S_pΔS+1zi2,(36)E2i;p1,p2,α1,α2=v^Sn+1,zi-v^Sn,ziΔS-∂∂Sv^Sn+1,zi2=1/ziS¯c-K2-anzigzi;p1,p2,α1,α2ΔS-1zi2.After determining τp(S) and τc(S) we can compute v^n(z) by solving linear system (23), and the value of American Strangles can be expressed in terms of the inverse Laplace transforms(37)vnSj,τ=L-1v^nSj,z,j=1,2,…,n.Applying polynomial interpolations, we can obtain the option value v(S,τ) for any values S∈(0,∞). Moreover, the vector version of inverse Laplace can be formulated as(38)vnτ=L-1v^nz.One of the classical numerical Laplace inversion is the Gaver-Stehfest formula (GSF) (see [21, 22]):(39)vnτ≈vn,Lτ=ln⁡2τ∑k=1LCkLv^nkln⁡2τ,where (40)CkL=-1k+L/2∑j=k+1/2min⁡k,L/2jL/22j!L/2-j!j!j-1!k-j!2j-k!with L 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 Algorithm2.Algorithm 2:NLTM algorithm for pricing American Strangles.
Step  1 . Use Algorithm 1 to compute S_p and S¯c. Step  2 . Compute τp(S) and τc(S) by optimizing (34). Step  3 . Solve the linear system (23) to get v^n(z). Step  4 . Apply Laplace inversion GSF or HCIM to restore American Strangles option valuevn(τ).


## 3. Hyperbolic Contour Integral Method for Laplace Inversion

Many computational experiences show that the Gaver-Stehfest inverse Laplace formula (39) is unstable when L≥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 vn(τ) from v^n(τ).The inversion of Laplace transform is based on numerical integration of the Bromwich complex contour integral:(41)vS,τ=12πi∫η-i∞η+i∞ezτv^S,zdz,η>η0,where i=-1 and η0 is the convergence abscissa. This means that all the singularities of v^(S,z) (with respect to z) lie in the open half-plane Re(z)≤η0. The integral (41) is not well-suited for numerical integration. First, the exponential factor is highly oscillatory on the Bromwich line, z=η+iy, -∞<y<+∞. Second, the transform v^(S,z) typically decays slowly as y→∞. 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(z)→-∞ at each end. On such a contour, the exponential factor in (41) forces a rapid decay of the integrand as Re(z)→-∞. 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 v^(S,z)→0 uniformly in Re(z)≤η0 as z→∞.Because the coefficient matrixA 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:(42)zI-A~v^nz=Cv^nz+gz,where I is the identity matrix,(43)A~=DB,(44)B=1ΔS2-1120⋯0012-112⋯00012-1⋱00⋮⋱⋱⋱⋱⋮00⋱12-11200⋯012-1,D=δS12β+2δS22β+2⋱δSn2β+2,(45)C=-rr-qS12ΔS0⋯00-r-qS22ΔS-rr-qS22ΔS⋯000-r-qS22ΔS-r⋱00⋮⋱⋱⋱⋱⋮00⋱-r-qSn-12ΔS-rr-qSn-12ΔS00⋯0-r-qSn2ΔS-rand g(z) is defined by (27). Note that B is a Toeplitz matrix. The linear system (42) can be solved by the following iteration algorithm:(46)v^nzl+1=zI-A~-1Cv^nzl+gz,l=0,1,…with [v^n(z)](0)=0.LetD-1-inner product and the corresponding vector norm be defined by(47)v,vD-1=v∗D-1v,vD-1=v,vD-1,∀v∈C,with D being defined by (44). Matrix D-1-norm is induced by the vector D-1-norm, i.e.,(48)AD-1=maxv≠0⁡AvD-1vD-1,∀A∈Cn×n.The theorem below shows that the iteration algorithm (46) is convergent if z falls outside of a sectorial region.Theorem 2.
(i) The spectrumΛ(A~) of A~ defined by (43) lies in the sectorial region Σϵ for any values of ϵ∈0,π/2, i.e.,(49)ΛA~⊆Σϵ=z∈C:arg⁡-z≤ϵ.(ii) For A~ defined by (43), we have(50)zI-A~-1D-1≤11+distz,Σϵ,z∉Σϵ,with Σϵ being defined by (49) and dist(z,Σϵ)=inf{z-ξ:ξ∈Σϵ}.
(iii) LetCD-1≤γ0 with C being defined by (45) and assume dist(z,Σϵ)≥γ such that 0<γ0<1+γ. Assuming v^n(z) solve (42), we have(51)v^nzl+1-v^nzD-1=γ01+γl+1v^nzD-1.Therefore, [v^n(z)](l+1) is convergent to v^n(z) according to ·D-1-norm.
(iv) Under the assumptions of (iii), we have(52)v^nzD-1≤11+γ-γ0gzD-1.Proof.
(i) MatrixB defined by (44) is a Toeplitz matrix and the generating function f(η) of B is(53)fη=1ΔS212eiη-e-0·iη+12e-iη=1ΔS2cos⁡η-1,η∈-π,π.From the above generating function f(η) 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 Λ(A~)=Λ(DB) 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(54)v^nzl+1-v^nzD-1≤zI-A~-1CD-1v^nzl-v^nzD-1≤zI-A~-1D-1CD-1v^nzl-v^nzD-1≤11+distz,ΣϵCD-1v^nzl-v^nzD-1≤γ01+γv^nzl-v^nzD-1,forl=0,1,2,… Therefore,(55)v^nzl+1-v^nzD-1≤γ01+γv^nzl-v^nzD-1≤γ01+γl+1v^nz0-v^nzD-1=γ01+γl+1v^nzD-1,and(56)liml→∞⁡v^nzl+1-v^nzD-1=0,which means [v^n(z)](l+1) is convergent to v^n(z) according to ·D-1-norm.
(iv) From linear equation (42) we have (57)v^nzD-1≤zI-A~-1D-1CD-1v^nzD-1+zI-A~-1D-1gD-1≤γ01+γv^nzD-1+11+γgD-1.Consequently, (58)v^nzD-1≤11+γ-γ0gzD-1,which ends the proof.From Theorem2, we see the iteration sequence {[v^nz]l,l=0,1,2,…} and the limit v^n(z) are analytical if z lies out the sectorial region Σϵ defined by (49). Therefore, the convergence of iteration algorithm (46) requires that z is far away from the sectorial region Σϵ such that (zI-A~)-1D-1CD-1<1. This condition could be satisfied in next discussion.The Laplace inversion ofv^n(z) with respect to z can be expressed as(59)vnτ=12πi∫Γezτv^nzdz.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(60)Γ:zζ=μ1+sin⁡iζ-θ,-∞<ζ<∞,where parameters μ>0 and θ set the width and the asymptotic angle of the hyperbolic contour, respectively. Then substituting the contour (60) into (59) gives(61)vnτ=12πi∫Γezτv^nzdz=12πi∫-∞+∞ezζτz′ζv^nzζdζ.Discretization of this integral with uniform node spacing h yields(62)vnτ=h2πi∑k=-∞+∞ezζkτz′ζkv^nzζk+DE+h+DE-h(63)=h2πi∑k=-LLezζkτz′ζkv^nzζk+DE+h+DE-h+TEhL(64)≔vn,Lτ+DE+h+DE-h+TEhL,where ζk=kh for the trapezoidal rule and ζk=(k+1/2)h for the midpoint rule, L is the number of quadrature nodes, DE±(h) are the discretization errors, and TE(hL) is the truncation error. Define(65)gτζ=12πiezζτz′ζv^nzζand assume the function gτ(ζ+id) to be analytic in the strip d∈(d-,d+) with d+>0 and d-<0, then the discretization errors are bounded as(66)DE+h=Md+e2πd+/h-1,DE-h=Md-e-2πd-/h-1,where(67)Md=∫-∞+∞gτζ+iddζ,d=d+ord-and · is some vector norm. The truncation error TE(hL) can be approximated by the magnitude of the last term retained, i.e.,(68)TEhL=OgτhL,L→∞.We note that the estimates in (66) and (68) are extended in a straightforward manner, from the case gτ(ζ) being scalar functions (see [24–26]) to the case when gτ(ζ) 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 h, μ, and θ must be chosen appropriately. By asymptotically matching DE±(h) and TE(hL), h and μ are derived as follows:(69)h=ρθL,μ=4πθ-π2+2πϵρθLτ,ρθ=cosh-12θ4θ-π+2ϵsin⁡θand θ is a free parameter satisfying θ<π/2-ϵ. With this choice the predicted convergence rate is given by(70)vn,Lτ-vnτ≤Oe-RθL,where(71)Rθ=π2-2πθ-2πϵρθ.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 R(θ). Consequently, the optimal parameters h and μ can be obtained from formulas (69).Denotingzk=z(ζk) and zk′=z′(ζk), we have(72)vn,Lτ=h2πi∑k=-LLezkτz′ζkv^nzk≈vn,l,Lτ≔h2πi∑k=-LLezkτz′ζkv^nzkl.Now, we give the convergence condition of the iteration algorithm (46). From [27, Theorem 3.3], the distance between Σϵ and Γ is given by (73)distΣϵ,Γ≔minP1∈Σφ,P2∈Γ⁡P1-P2=μ1-sin⁡θ.Inserting the above equation and the expression μ=4πθ-π2+2πϵ/ρ(θ)L/τ into following inequality (74)zI-A~-1D-1≤11+distΣϵ,Γwe get (75)zI-A~-1D-1≤11+4πθ-π2+2πϵ1-sin⁡θ/ρθL/τ.Denote γ=μ(1-sin⁡θ)=dist(Σϵ,Γ). If we take L such that(76)Lτ>ρθγ0-14πθ-π2+2πϵ1-sin⁡θ,then 1+γ>γ0 and the iteration algorithm (46) is convergent according to norm ·D-1. Inequality (76) requires that L is large enough such that the hyperbolic contour Γ is far away from the sectorial region Σϵ.In following Theorem, we give an estimation forCD-1 and then give the minimum estimation for L.Theorem 3.
For matrixC being defined by (45), we have the following estimation:(77)CD-1≤γ0≔MS_p,S¯c2ΔS8r2-4rr-qΔS2+3r-q2S_p2,where S¯c and S_p are defined by (8) in Section 2 and(78)MS_p,S¯c=S¯cS_pmaxj⁡Sjβ+1minj⁡Sjβ+1=S¯cS_pβ+2,β≥-1,S¯pS_cβ,β<-1.If ΔS is small enough, then(79)CD-1≤γ0≔3r-qS_pMS_p,S¯c2ΔS=O1ΔS,ΔS→0.Proof.
Firstly, we verified the inequality(80)CD-1≤maxj⁡Sjβ+1minj⁡Sjβ+1C,where · represents the vector l2-norm. Indeed,(81)CD-1=maxv≠0⁡CvD-1vD-1=maxv≠0⁡Cv∗D-1Cvv∗D-1v≤maxv≠0⁡1/minj⁡δSj2β+2Cv∗Cv1/maxj⁡δSj2β+2v∗v=maxj⁡Sj2β+2minj⁡Sj2β+2maxv≠0⁡Cvv=maxj⁡Sjβ+1minj⁡Sjβ+1C.From definition (45), we have C=1/2ΔSD1C1 with(82)C1=-2rΔSS1r-q0⋯00-r-q-2rΔSS2r-q⋯000-r-q-2rΔSS3⋱00⋮⋱⋱⋱⋱⋱00⋱-r-q-2rΔSSn-1r-q00⋯0-r-q-2rΔSSnand D1=DiagS1,S2,…,Sn. By carefully calculating, we have(83)C1⊺C1=cjkn×nwhere(84)cjj=4r2ΔS2Sj2+r-q2,j=1,n,cjj=8r2ΔS2Sj2+r-q2,j=2,3,…,n-1,cj,j-1=cj,j+1=2rr-qΔS2SjSj+1,j=1,2,…,n,cj,j-2=cj,j+2=-r-q2,j=1,2,…,n,and cjk=0 for j=1,2,…,n, k<j-2 and k>j+2. By the Gershgorin circle theorem, each eigenvalue λ(C1⊺C1) of C1⊺C1 falls into a certain Gerschgorin disk Gj, i.e.,(85)λC1⊺C1)∈Gj≔z:z-cjj≤∑k=j-2,k≠jj+2cjkfor some index j. Incorporating (84) and (85), we get(86)λC1⊺C1≤maxj⁡8r2ΔS2Sj2+4rr-qΔS2SjSj+1+3r-q2≤1minj⁡Sj28r2-4rr-qΔS2+3r-q2minj⁡Sj2.Therefore,(87)C1≤λC1⊺C1≤1minj⁡Sj8r2-4rr-qΔS2+3r-q2minj⁡Sj2.Incorporating (80) and (87), we have(88)CD-1≤maxj⁡Sjβ+1minj⁡Sjβ+1C=12ΔSmaxj⁡Sjβ+1minj⁡Sjβ+1D1C1≤12ΔSmaxj⁡Sjβ+1minj⁡Sjβ+1D1C1≤12ΔSmaxj⁡Sjβ+1minj⁡Sjβ+1maxj⁡Sjminj⁡Sj8r2-4rr-qΔS2+3r-q2minj⁡Sj2Let maxj⁡Sj=S¯c, minj⁡Sj=S_p and M(S_p,S¯c)=S¯c/S_pmaxj⁡Sjβ+1/minj⁡Sjβ+1; 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 L=O1/ΔS. (ii) For inverse Laplace formula vn,L(τ) in (72), we solve v^n(z) by directly computing linear equations v^n(z)=(I-A)-1g(z) (i.e., expression (23) in Section 2) if L satisfies inequality (76) with γ0 being defined by (77).Remark 5.
From the upper estimate (52) for v^n(z)D-1, we see the convergence of inverse Laplace dependents on the norm g(z)D-1. Moreover, observing the expression of g(z) defined by (27) and (29), the exponential terms e-zτp(Sj;p1,α1) and e-zτc(Sj;p2,α2) play important roles. Following the strategy of Zhou et al. [15], we split gj(z) into two parts gj(z)=g1j(z)+g2j(z) with (89)g1j=max⁡K1-Sj,0+max⁡Sj-K2,0+1zrK1-qSj-1zrK2-qSj-1zrK1-qSje-zτcSj1τ>τcSj+1zrK2-qSje-zτcSj1τ>τpSjg2j=-1zrK1-qSje-zτcSj1τ<τcSj+1zrK2-qSje-zτcSj1τ<τpSjfor j=1,2,…,n. Therefore, we have splitting expression g(z)=g1(z)+g2(z). ezτg1(z) rapidly decays as Re(z)→-∞ if z falls on the hyperbola contour Γ defined by (60). To make the Laplace transform method convergent, we define an other hyperbola contour Γ~:(90)Γ~:zζ=μ1+sin⁡iζ+θ,-∞<ζ<∞,which begins and ends in the right half-plane, such that Re(z)→+∞ at each end. We select the same parameters ϵ, μ, and θ as those of Γ. Obviously, ezτg2(z) rapidly decays as Re(z)→+∞ if z falls on the hyperbola contour Γ~ defined by (90). The inverse Laplace of v^n(z) with respect to z can be expressed as(91)vnτ=12πi∫Γezτv^nzdz=12πi∫ΓezτzI-A-1g1z+g2zdz=12πi∫ΓezτzI-A-1g1zdz+12πi∫Γ~ezτzI-A-1g2zdz:≈v1n,Lτ+v2n,Lτ,where v1n,L(τ) and v2n,L(τ) have the predicted convergence rate R(θ) defined by formula (70).


## 4. 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)-(7) are taken as T=1(year), r=0.05(1/year), q=0.1($), K1=1($), K2=1.5($), S0=1.5, δ=σ0/S0β and Smax=4K2. The values of parameters in HCIM are taken as θ=0.8 and ϵ=0.4. We take the number of space partition N=3000 to compute S_p and S¯c and take n=3000 to optimise τp(S) and τc(S). The positive values zi=1+3i/m with m=6 in expressions (35) and (36). The size of space mesh and the time mesh is set as N=2000, Δt=T/2000 for FDM.Example 1.
Table1 lists the American Strangles option prices with σ0=0.4, different values of β and different values of L. From the table, we see LTM-GSF is unstable when L>20, while the values of LTM-HCIM is very close to ones of FDM for L=10,20,26. So, LTM-HCIM outperforms the LTM-GSF formula in regard to the accuracy and stability for the option values.Table 1
Values of American Strangles at timet=0 (i.e., τ=T) with K1=1, K2=1.5, r=0.05, q=0.1, S0=1.5, σ0=0.4, different values of β, and different values of L. The CPU time is about 20s, 22s, and 680s for LTM-GSF, LTM-HCIM, and FDM, respectively.
β L=10 L = 20 L = 26 FDM LTM-GSF LTM-HCIM LTM-GSF LTM-HCIM LTM-GSF LTM-HCIM -0.50 0.260862 0.260864 0.297905 0.260864 18.447801 0.260864 0.258805 -0.25 0.253815 0.253800 0.284163 0.253800 22.070192 0.253800 0.251154 0 0.247986 0.247964 0.261926 0.247964 24.667532 0.247964 0.244310 0.25 0.245120 0.245090 0.239312 0.245090 -0.098661 0.245090 0.238169 0.50 0.236749 0.236703 0.217691 0.236724 -9.087067 0.236724 0.232615Example 2.
This example illustrates the convergence rate of LTM-HCIM for the numbern of space mesh and the number L for HCIM. Table 2 lists the space convergence rate with L=20 and fixed N=3000. The column labeled “Value” shows the option values at stock price S0=1.5; the column labeled “Err” is defined by(92)Errn=ΔSnvn,Lτ-Iv2n,Lτ,where I[v2n,L(τ)] is the linear interpolation of v2n,L(τ) from 2n-mesh to n-mesh, · is the l2 vector norm and ΔS(n)=S¯c-S_p/n+1, and the convergence rate is defined by(93)Raten=log2⁡ErrnErr2n.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 solutionvn,L(τ) with n=1599,L=30 as the reference option prices, i.e., vn(τ)≈v1599,30(τ), and the LTM-HCIM time-errors are defined by(94)ErrL=ΔS1600v1599,Lτ-v1599,30τ,for different values of L. Figure 1 plots the LTM-HCIM errors in log-scale v.s. L. From this Figure we see the numerical convergence rates are consistent with the predicted rates (95)vn,Lτ-vnτ≤Oe-RθL,for moderate values of L.Table 2
Values of American Strangles at timet=0 (i.e., τ=T) with K1=1, K2=1.5, r=0.05, q=0.1, S0=1.5, L=20, different values of σ0, different values of β, and different space partition n.
n σ 0 = 0.2, β=0.5 σ 0 = 0.4, β=-0.5 Value Err Rate Value Err Rate 24 0.150660 3.3277e-003 0.6061 0.262272 2.7302e-003 1.5515 49 0.149394 2.1862e-003 0.5367 0.260271 9.3141e-004 0.9603 99 0.149262 1.5070e-003 0.5154 0.260335 4.7871e-004 1.0611 199 0.149211 1.0543e-003 0.5076 0.260319 2.2943e-004 0.9601 399 0.149187 7.4158e-004 0.5037 0.260274 1.1793e-004 0.9812 799 0.149191 5.2302e-004 0.5018 0.260276 5.9737e-005 0.9920 1599 0.149189 3.6937e-004 — 0.260275 3.0034e-005 — 3199 0.149189 — — 0.260275 — —Figure 1
LTM-HCIM errors in log-scale v.s.L for American CEV Strangles option at τ=1. Left: log scale errors with σ0=0.2, β=0.5; Right: log scale errors with σ0=0.4, β=-0.5.


## 5. 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.



---

*Source: 5908646-2018-09-04.xml*