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

1. IntroductionThe 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 Vasiljevic [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 [[S.sub.p]([tau]),[S.sub.c]([tau])], where [S.sub.p]([tau]) is put side and [S.sub.c]([tau]) is call side. The idea of the new method is modifying PDE into a new form on certain fixed region [[[S.bar].sub.p], [[bar.S].sub.c]]. Then performing the Laplace transform leads to a ODE which involves the inverse functions [[tau].sub.p](S) and [t.sub.c](S) of [S.sub.p]([tau]) and [S.sub.c]([tau]), respectively. Using the finite difference method combined with the approximation of functions [[tau].sub.p](S) and [[tau].sub.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 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]):

d[S.sub.t] = (r-q) [S.sub.t]dt + [delta][S.sup.[beta]+1.sub.t]d[W.sub.t], (1)

where r > 0 is the risk-free interest rate, q [greater than or equal to] 0 is the dividend yield, and [W.sub.t] is standard Brownian motion. [sigma](S) = [delta][S.sup.[beta]] represents the local volatility function and [beta] can be interpreted as the elasticity of [sigma](S). [delta] is the scale parameter fixing the initial instantaneous volatility at time t = 0, i.e., [delta] = [sigma]([S.sub.0])/[S.sup.[beta].sub.0] [??] [[sigma].sub.0]/[S.sup.[beta].sub.0]. If [beta] = 0, the SDE (1) becomes the standard log-normal diffusion model or so-called geometric Brownian motion models (GBM). If [beta] = -1/2, the SDE (1) nests the Cox-Ingersoll-Ross (CIR) model.

Let

[mathematical expression not reproducible] (2)

be the price of an American Strangles position written on an underlying asset with price at time and the payoff

[PI] ([S.sub.[tau]]) = max ([K.sub.1] - [S.sub.[tau]], 0) + max ([S.sub.[tau]] - [K.sub.2],0). (3)

This position is formed using a long put with strike [K.sub.1] and a long call with strike [K.sub.2]. Note that [K.sub.1] < [K.sub.2]. Let [tau] = T - t, v(S, [tau]) = f (S, T - t), and the early exercise boundary on the put side be denoted by [S.sub.p]([tau]), and the early exercise boundary on the call side be denoted by [S.sub.c]([tau]). It is known that v(S, [tau]) satisfies the Black-Scholes PDE:

[mathematical expression not reproducible] (4)

with initial condition

v (S, 0) = max ([K.sub.1] - S,0) + max (S - [K.sub.2],0), [S.sub.p] (0) < S < [S.sub.c] (0) (5)

and free boundary conditions

v (S,[tau]) = [K.sub.1] - S, [partial derivative]/[partial derivative]S v (S, [tau]) = -1; S [greater than or equal to] [S.sub.p] ([tau]), (6)

v (S,[tau]) = S - [K.sub.2], [[partial derivative]/[partial derivative]]v(S, [tau]) = 1, S [less than or equal to] [S.sub.c] ([tau]). (7)

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 [S.sub.p]([tau]) and [S.sub.c]([tau]) of American Strangles are continuously differentiable on the interval [0,+[infinity]). We call the fact that [S.sub.p]([tau]) is a decreasing function of variable [tau], while [S.sub.c]([tau]) is an increasing function of variable [tau]. Denote

[[S.bar].sub.p] = [S.sub.p]([infinity]), [[bar.S].sub.p] = [S.sub.p](0), [[S.bar].sub.c] = [S.sub.c] (0), [[bar.S].sub.p] = [S.sub.c] ([infinity]), (8)

Then we know

[[bar. S].sub.p] = min([K.sub.1], [r/q] [K.sub.1]), [[S.bar].sub.c] = max ([K.sub.2], [r/q] [K.sub.2]), (9)

[[S.bar].sub.p] > O, [[bar.S].sub.c] < [infinity] and [v.sup.[infinity]] (S) = v(S, [infinity]) satisfy the following ODE (similar to perpetual American option):

[mathematical expression not reproducible], (10)

[mathematical expression not reproducible], (11)

[mathematical expression not reproducible], (12)

Algorithm 1: Algorithm for solving [[bar.S].sub.p] and [[bar.S].sub.c]. 1: Step 1. Let [L.sub.1] = 1, [H.sub.1] = [j.sub.1], [L.sub.2] = [j.sub.2], [H.sub.2] = N - 1. 2: Step 2. Do the following loop: 3: while [L.sub.1] < [H.sub.1] - 1 or [L.sub.2] < [H.sub.2] - 1 do 4: Let [k.sub.1] = floor[0.5 * ([L.sub.1] + [H.sub.1])] and [k.sub.2] = floor[0.5 * ([L.sub.2] + [H.sub.2])]. 5: Solve linear system (13) with boundary [mathematical expression not reproducible] 6: % Finite difference approximation (6). 7: if [mathematical expression not reproducible] then 8: [L.sub.1] = [k.sub.1]; 9: else 10: [H.sub.1] = [k.sub.1]; 11: end if 12: % Finite difference approximation (7). 13: if [mathematical expression not reproducible] then 14: [L.sub.2] = [k.sub.2]; 15: else 16: [H.sub.2] = [k.sub.2]. 17: end if 18: end while 19: Step 3. Output [mathematical expression not reproducible].

The values of [[S.bar].sub.p] and [[bar.S].sub.p] could be computed by secant method. Define uniform mesh [S.sub.j] = j[DELTA]S, j = 0, 1, ..., N and take larger enough number [S.sub.max] = [S.sub.N] such that [S.sub.max] > [[bar.S].sub.c]. The discrete form of (10) is as follows:

[mathematical expression not reproducible] (13)

for j = 1,2, ..., N-1 with [v.sub.j] = [v.sup.[infinity]]([S.sub.j]). Assuming [mathematical expression not reproducible] and [mathematical expression not reproducible], 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:

[mathematical expression not reproducible] (14)

for [tau] [member of] (0, + [infinity]) and S [member of] ([[S.bar].sub.p], [[bar.S].sub.c]), where [1.sub.{x}] is the indicator function.

Proof, (i) On the region [tau] [member of] (0, [infinity]), S [member of] (0, [[S.bar].sub.p]([tau])), we know from the boundary condition (6) that v(S, [tau]) = [K.sub.1] - S. Direct calculation shows that v(S, [tau]) = [K.sub.1] -S satisfies (14). (ii) On the region [tau] [member of] (0, [infinity]), S [member of] ([[bar.S].sub.p]([tau]), [[bar.S].sub.c]([tau])), both (4) and (14) have the same forms, (iii) On the region [tau] [member of] (0, [infinity]), S [member of] ([[bar.S].sub.c]([tau]), [infinity]), we know from the boundary condition (7) that v(S, [tau]) = S - [K.sub.2]. Direct calculation shows that v(S, [tau]) = S-[K.sub.2] satisfies (14). Combining (i), (ii), and (iii), we get this theorem.

Let [[tau].sub.p](S) and [[tau].sub.c](S) be the inverse functions of [S.sub.p]([tau]) and [S.sub.c]([tau]), respectively. Taking Laplace transform

[mathematical expression not reproducible], (15)

[mathematical expression not reproducible] (16)

and

[mathematical expression not reproducible] (17)

to the both sides of (14), we get

[mathematical expression not reproducible] (18)

on fixed space region S [member of] ([[S.bar].sub.p], [[bar.S].sub.c]). The initial free boundary conditions (6) and (7) become the following forms:

[mathematical expression not reproducible], (19)

[mathematical expression not reproducible], (20)

Now, the free boundaries [[tau].sub.p](S) and [[tau].sub.c](S) can be approximated. To this end, we give the fit functions with unknown parameters [p.sub.i], [[alpha].sub.i] for i = 1, 2:

[mathematical expression not reproducible], (21)

[mathematical expression not reproducible], (22)

Redefine uniform mesh [S.sub.j] = [[S.bar].sub.p] + j[DELTA]S, j = 0, 1, ..., n + 1 with [DELTA]S = ([[bar.S].sub.c] - [[S.bar].sub.p])/(n + 1). Equations (18) can be discretized as follows:

(zI - A) [[??].sup.n] (z) = g (z), (23)

with index n representing the number of space mesh partition. Matrix A, solution vector [[??].sup.n](z), and RHS vector g(z) are defined as

[mathematical expression not reproducible], (24)

[mathematical expression not reproducible], (25)

[mathematical expression not reproducible], (26)

g (z; [p.sub.1], [p.sub.2], [[alpha].sub.1], [[alpha].sub.2]) = [[[g.sub.1], ..., [g.sub.n]].sup.T], (27)

where

[g.sub.1] = [a.sub.1] ([K.sub.1] - [[bar.S].sub.p]), (28)

[mathematical expression not reproducible] (29)

[g.sub.n] = [c.sub.n] ([[bar.S].sub.c] - [K.sub.2]), (30)

Denote [a.sub.1] (z) and [a.sub.n](z) to be the first row and the nth row of [(zI - A).sup.-1] (z), respectively. [a.sub.1] (z) and [a.sub.n](z) can be obtained by solving linear equation:

[a.sub.1] (z)(zI - A) = [1,0,0, ... ,0], (31)

[a.sub.n] (z)(zI - A) = [0,0, ..., 0,1]. (32)

Then, we reach the expressions

[mathematical expression not reproducible] (33)

We find the parameters [p.sub.i], [[alpha].sub.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 [[z.sub.1], [z.sub.2], ..., [z.sub.m]] which are used for numerical Laplace inversion and find the values of unknown [p.sub.i] and [[alpha].sub.i] that optimizing the following objective function:

[mathematical expression not reproducible] (34)

with

[mathematical expression not reproducible] (35)

[mathematical expression not reproducible] (36)

Algorithm 2: NLTM algorithm for pricing American Strangles. Step 1. Use Algorithm 1 to compute [[S.bar].sub.p] and [[bar.s].sub.c]. Step 2. Compute [[tau].sub.p] (S) and [[tau].sub.c] (S) by optimizing (34). Step 3. Solve the linear system (23) to get [[??].sup.n](z). Step 4. Apply Laplace inversion GSF or HCIM to restore American Strangles option value [v.sup.n]([tau]).

After determining [[tau].sub.p](S) and [[tau].sub.c](S) we can compute [[??].sup.n](z) by solving linear system (23), and the value of American Strangles can be expressed in terms of the inverse Laplace transforms

[v.sup.n] ([S.sub.j], [tau]) = [[Laplace].sup.-1] [[[??].sup.n] ([S.sub.j], z)], j = 1,2, ..., n. (37)

Applying polynomial interpolations, we can obtain the option value v(S, [tau]) for any values S [member of] (0, [infinity]). Moreover, the vector version of inverse Laplace can be formulated as

[v.sup.n] ([tau]) = [[Laplace].sup.-1] [[[??].sup.n] (z)]. (38)

One of the classical numerical Laplace inversion is the Gaver-Stehfest formula (GSF) (see [21, 22]):

[mathematical expression not reproducible] (39)

where

[mathematical expression not reproducible], (40)

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 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 L [greater than or equal to] 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 [v.sup.n]([tau]) from [[??].sup.n]([tau]).

The inversion of Laplace transform is based on numerical integration of the Bromwich complex contour integral:

[mathematical expression not reproducible], (41)

where i = [square root of -1] and [[eta].sub.0] is the convergence abscissa. This means that all the singularities of [??](S, z) (with respect to z) lie in the open half-plane Re(z) [less than or equal to] [[eta].sub.0]. The integral (41) is not well-suited for numerical integration. First, the exponential factor is highly oscillatory on the Bromwich line, z = q + 1y, -[infinity] < y < + [infinity]. Second, the transform [??](S,z) typically decays slowly as [absolute value of y] [right arrow] [infinity]. One strategy for circumventing the slow decay is due to Talbot [23], who suggests that the Bromwich line be deformed into a contour [GAMMA] that begins and ends in the left half-plane, such that Re(z) [right arrow] [infinity] at each end. On such a contour, the exponential factor in (41) forces a rapid decay of the integrand as Re(z) [right arrow] -[infinity]. 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 [??](S,z) [right arrow] 0 uniformly in Re(z) [less than or equal to] [[eta].sub.0] as [absolute value of z] [right arrow] [infinity].

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:

[mathematical expression not reproducible], (42)

where I is the identity matrix,

[mathematical expression not reproducible], (43)

[mathematical expression not reproducible] (44)

[mathematical expression not reproducible] (45)

and 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:

[mathematical expression not reproducible] (46)

with [[[[??].sup.n](z)].sup.(0)] = 0.

Let [D.sup.-1]-inner product and the corresponding vector norm be defined by

[mathematical expression not reproducible] (47)

with D being defined by (44). Matrix [D.sup.-1]-norm is induced by the vector [D.sup.-1]-norm, i.e.,

[mathematical expression not reproducible] (48)

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 [??] defined by (43) lies in the sectorial region [[SIGMA].sub.[epsilon]] for any values of [epsilon] [epsilon] (0, [pi]/2), i.e.,

[mathematical expression not reproducible]. (49)

(ii) For [??] defined by (43), we have

[mathematical expression not reproducible] (50)

with [[SIGMA].sub.[epsilon]] being defined by (49) and dist(z, [[SIGMA].sub.[epsilon]] = inf {[absolute value of z-[xi]] : [xi] [member of] [[SIGMA].sub.[epsilon]]}.

(iii) Let [mathematical expression not reproducible] with C being defined by (45) and assume dist(z, [[SIGMA].sub.[epsilon]]) [greater than or equal to] [gamma] such that 0< [[gamma].sub.0] < 1 + [gamma]. Assuming [[??].sup.n](z) solve (42), we have

[mathematical expression not reproducible] (51)

Therefore, [[[[??].sup.n](z)].sup.(l+1)] is convergent to [[??].sup.n](z) according to [mathematical expression not reproducible]-norm.

(iv) Under the assumptions of (iii), we have

[mathematical expression not reproducible] (52)

Proof. (i) Matrix B defined by (44) is a Toeplitz matrix and the generating function f([eta]) of B is

[mathematical expression not reproducible] (53)

From the above generating function f([eta]) and the discussion of Pang and Sun [24], we know the spectrum [LAMBDA](B) lies the sectorial region [[SIGMA].sub.[epsilon]] for any values of [member of] [member of] (0, [pi]/2). Since D is positive definite, the spectrum [LAMBDA]([??]) = [LAMBDA](DB) also lies in the same sectorial region [[SIGMA].sub.[epsilon]] (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

[mathematical expression not reproducible] (54)

Therefore,

[mathematical expression not reproducible] (55)

[mathematical expression not reproducible] (56)

which means [[[[??].sup.n](z)].sup.(l+1)] is convergent to [[??].sup.n](z) according to [mathematical expression not reproducible]-norm.

(iv) From linear equation (42) we have

[mathematical expression not reproducible] (57)

Consequently,

[mathematical expression not reproducible] (58)

which ends the proof.

From Theorem 2, we see the iteration sequence {[[??].sup.n] (z)](l) l = 0,1,2, ...} and the limit [[??].sup.n](z) are analytical if z lies out the sectorial region [[SIGMA].sub.[epsilon]] defined by (49). Therefore, the convergence of iteration algorithm (46) requires that z is far away from the sectorial region [[SIGMA].sub.[epsilon]] such that [mathematical expression not reproducible]. This condition could be satisfied in next discussion.

The Laplace inversion of [[??].sup.n](z) with respect to z can be expressed as

[mathematical expression not reproducible]. (59)

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 r, which can be parameterized by

[mathematical expression not reproducible], (60)

where parameters [mu] > 0 and [theta] set the width and the asymptotic angle ofthe hyperbolic contour, respectively. Then substituting the contour (60) into (59) gives

[mathematical expression not reproducible] (61)

Discretization of this integral with uniform node spacing h yields

[mathematical expression not reproducible] (62)

[mathematical expression not reproducible] (63)

[mathematical expression not reproducible], (64)

where [[zeta].sub.k] = kh for the trapezoidal rule and [[zeta].sub.k] = (k + 1/2)h for the midpoint rule, L is the number of quadrature nodes, [DE.sub.[+ or -]](h) are the discretization errors, and TE(hL) is the truncation error. Define

[mathematical expression not reproducible] (65)

and assume the function [g.sub.r]([[zeta].sub.k] + id) to be analytic in the strip d [member of] ([d.sub.-], [d.sub.+]) with [d.sub.+] > 0 and [d.sub.-] < 0, then the discretization errors are bounded as

[mathematical expression not reproducible] (66)

where

[mathematical expression not reproducible] (67)

and [parallel]*[parallel] is some vector norm. The truncation error TE(hL) can be approximated by the magnitude of the last term retained, i.e.,

[mathematical expression not reproducible]. (68)

We note that the estimates in (66) and (68) are extended in a straightforward manner, from the [g.sub.[tau]]([zeta]) being scalar functions (see [24-26]) to the case when [g.sub.[tau]]([zeta]) takes value in a complex vector space.

To make all the singularities of the integrand in (61) to fall into the sectorial region [[SIGMA].sub.[epsilon]] defined by (49) for a given [epsilon] > 0, the expressions of optimal parameters h, [mu], and [theta] must be chosen appropriately. By asymptotically matching [parallel][DE.sub.[+ or - ]](h)[parallel] and [parallel]TE(hL)[parallel], h and [mu] are derived as follows:

[mathematical expression not reproducible]. (69)

and [theta] is a free parameter satisfying [theta] < [pi]/2 - [epsilon]. With this choice the predicted convergence rate is given by

[mathematical expression not reproducible], (70)

where

R([theta]) = [[pi].sup.2]-2[pi][theta]-2[pi][epsilon]/[rho]([theta]). (71)

According to (70), we see that once the semiangle [epsilon] of the sectorial region (49) is given, then one can find the optimal [theta] by maximizing the function R([theta]). Consequently, the optimal parameters h and [mu] can be obtained from formulas (69).

Denoting [z.sub.k] = z([[zeta].sub.k]) and [z'.sub.k] = z ([[zeta].sub.k]), we have

[mathematical expression not reproducible] (72)

Now, we give the convergence condition of the iteration algorithm (46). From [27, Theorem 3.3], the distance between [[SIGMA].sub.[epsilon]] and [GAMMA] is given by

[mathematical expression not reproducible]. (73)

Inserting the above equation and the expression [mu] = ((4[pi][theta]- [[pi].sup.2] + 2[pi][epsilon])/[rho]([theta]))(L/[tau]) into following inequality

[mathematical expression not reproducible] (74)

we get

[mathematical expression not reproducible] (75)

Denote [gamma] = [mu](1 - sin d) = dist([[SIGMA].sub.[epsilon]], [GAMMA]). If we take L such that

[mathematical expression not reproducible] (76)

then 1 + [gamma] > [[gamma].sub.0] and the iteration algorithm (46) is convergent according to norm [mathematical expression not reproducible]. Inequality (76) requires that L is large enough such that the hyperbolic contour [GAMMA] is far away from the sectorial region [[SIGMA].sub.[epsilon]].

In following Theorem, we give an estimation for [mathematical expression not reproducible] and then give the minimum estimation for L.

Theorem 3. For matrix C being defined by (45), we have the following estimation:

[mathematical expression not reproducible] (77)

where [[bar.S].sub.c] and [[S.bar].sub.p] are defined by (8) in Section 2 and

[mathematical expression not reproducible] (78)

If [DELTA]S is small enough, then

[mathematical expression not reproducible] (79)

Proof. Firstly, we verified the inequality

[mathematical expression not reproducible], (80)

where [parallel]*[parallel] represents the vector [l.sub.2]-norm. Indeed,

[mathematical expression not reproducible] (81)

From definition (45), we have C = (l/2[DELTA]S)[D.sub.1][C.sub.1] with

[mathematical expression not reproducible] (82)

and [D.sub.1] = Diag[[S.sub.1], [S.sub.2], ..., [S.sub.n]. By carefully calculating, we have

[mathematical expression not reproducible] (83)

where

[mathematical expression not reproducible] (84)

and [c.sub.jk] = 0 for j = 1, 2, ..., n, k < j - 2 and k > j + 2. By the Gershgorin circle theorem, each eigenvalue [lambda]([C.sup.T.sub.1][C.sub.1]) of [C.sup.T.sub.1][C.sub.1] falls into a certain Gerschgorin disk [G.sub.j], i.e.,

[mathematical expression not reproducible] (85)

for some index j. Incorporating (84) and (85), we get

[mathematical expression not reproducible] (86)

Therefore,

[mathematical expression not reproducible] (87)

Incorporating (80) and (87), we have

[mathematical expression not reproducible] (88)

Let [max.sub.j][S.sub.j] = [[bar.S].sub.c], [min.sub.j][S.sub.j] = [[S.bar].sub.p] and M([[S.bar].sub.p], [[bar.S].sub.c]) = ([[bar.S].sub.c]/ [[S.bar].sub.p])([max.sub.j][S.sup.[beta]+1.sub.j]/[min.sub.j][S.sup.[beta]+1.su b.j]); we finally obtain the estimate expression (77).

Finally, we give two remarks as follows.

Remark 4. (i) For the fixed hyperbolic contour parameters [theta] and [epsilon], Theorem 3 and iteration convergence condition (76) show that we must take L = O(1/[DELTA]S). (ii) For inverse Laplace formula [v.sup.n,L]([tau]) in (72), we solve [[??].sup.n](z) by directly computing linear equations [[??].sup.n](z) = [(I - A).sup.-1] g(z) (i.e., expression (23) in Section 2) if L satisfies inequality (76) with [[gamma].sub.0] being defined by (77).

Remark 5. From the upper estimate (52) for [mathematical expression not reproducible], we see the convergence of inverse Laplace dependents on the norm [mathematical expression not reproducible]. Moreover, observing the expression of g(z) defined by (27) and (29), the exponential terms [mathematical expression not reproducible] and [mathematical expression not reproducible] play important roles. Following the strategy of Zhou et al. [15], we split [g.sub.j](z) into two parts [g.sub.j](z) = [g.sub.1j](z) + [g.sub.2j]z) with

[mathematical expression not reproducible] (89)

for j = 1,2, ..., n. Therefore, we have splitting expression g(z) = [g.sub.1](z)+ [g.sub.2](z). [e.sup.Z[T.sub.g1](z) rapidly decays as Re(z) [right arrow] -[infinity] if z falls on the hyperbola contour r defined by (60). To make the Laplace transform method convergent, we define an other hyperbola contour [??]:

[??]: z([zeta] = [mu][1 + sin (1[zeta] + [theta])], [infinity] < [zeta] < [infinity], (90)

which begins and ends in the right half-plane, such that Re(z) [right arrow] +[infinity] at each end. We select the same parameters [epsilon], [mu], and [theta] as those of [GAMMA]. Obviously, [mathematical expression not reproducible] rapidly decays as Re(z) [right arrow] +[infinity] if z falls on the hyperbola contour [??] defined by (90). The inverse Laplace of [[??].sup.n](z) with respect to z can be expressed as

[mathematical expression not reproducible] (91)

where [v.sup.n,L.sub.1]([tau]) and [v.sup.n,L.sub.2]([tau]) have the predicted convergence rate R([theta]) 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($), [K.sub.1] = 1($), [K.sub.2] = 1.5($), [S.sub.0] = 1.5, [delta] = [[sigma].sub.0]/[S.sup.[beta].sub.0] and [S.sub.max] = 4[K.sub.2]. The values of parameters in HCIM are taken as [theta] = 0.8 and [epsilon] = 0.4. We take the number of space partition N = 3000 to compute [[S.bar].sub.p] and [[bar.S].sub.c] and take n = 3000 to optimise [[tau].sub.p](S) and [[tau].sub.c](S). The positive values [z.sub.i] = 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, [DELTA]t = T/2000 for FDM.

Example 1. Table 1 lists the American Strangles option prices with [[sigma].sub.0] = 0.4, different values of [beta] 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 LTMGSF formula in regard to the accuracy and stability for the option values.

Example 2. This example illustrates the convergence rate of LTM-HCIM for the number n 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 [S.sub.0] = 1.5; the column labeled "Err" is defined by

[mathematical expression not reproducible], (92)

where I[[v.sup.2n,L]([tau])] is the linear interpolation of [v.sup.2n,L]([tau]) from 2n-mesh to n-mesh, [parallel]*[parallel] is the [l.sub.2] vector norm and [DELTA]S(n) = ([[bar.S].sub.c] - [[S.bar].sub.p])/(n + 1), and the convergence rate is defined by

Rate (n) = [log.sub.2] [Err (n)/Err (2n)]. (93)

From Table 2, we see the convergence rate of LTM-HCIM is 0.5 and 1.0 for [beta] = 0.5 and [beta] = -0.5, respectively. It is very interesting to investigate the relationship between the convergence rates and the values of [beta] and we leave this topic to the future.

We regard the computational solution [v.sup.n,L]([tau]) with n = 1599, L = 30 as the reference option prices, i.e., [v.sup.n]([tau]) [approximately equal to] [v.sup.1599,30]([tau]), and the LTM-HCIM time-errors are defined by

Err (L) = [square root of ([DELTA]S (1600))] [parallel][v.sup.1599, L] ([tau]) - [v.sup.1599,30] ([tau])[parallel], (94)

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

[mathematical expression not reproducible], (95)

for moderate values of L.

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 g[infinity]d 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.

https://doi.org/10.1155/2018/5908646

Data Availability

All Matlab codes supporting the results of this study are available from the corresponding author on request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors' Contributions

Hongying Wu carried out the experiments in Section 4 and Zhiqiang Zhou made the main contributions to the other sections. All authors read and approved the final manuscript.

Acknowledgments

The work was supported by Natural Science Foundation of Hunan Province of China (no. 2017JJ2113).

References

[1] P. Carr, "Randomization and the American put," Review of Financial Studies, vol. 11, no. 3, pp. 597-626,1998.

[2] A. Pelsser, "Pricing double barrier options using Laplace transforms," Finance and Stochastics, vol. 4, no. 1, pp. 95-104,2000.

[3] D. Davydov and V. Linetsky, "Pricing and hedging path-dependent options under the CEV process," Management Science, vol. 47, no. 7, pp. 949-965, 2001.

[4] A. Sepp, "Analytical pricing of double-barrier options under a double-exponential jump diffusion process: applications of Laplace transform," International Journal of Theoretical and Applied Finance, vol. 7, no. 2, pp. 151-175, 2004.

[5] N. Cai and S. G. Kou, "Option pricing under a mixed-exponential jump diffusion model," Management Science, vol. 57, no. 11, pp. 2067-2081, 2011.

[6] N. Cai and S. Kou, "Pricing Asian options under a hyperexponential jump diffusion model," Operations Research, vol. 60, no. 1, pp. 64-77, 2012.

[7] R. Mallier and G. Alobaidi, "Laplace transforms and American options," Applied Mathematical Finance, vol. 7, no. 4, pp. 241256, 2000.

[8] S.-P. Zhu, "A new analytical approximation formula for the optimal exercise boundary of American put options," International Journal of Theoretical and Applied Finance, vol. 9, no. 7, pp. 1141-1177, 2006.

[9] T. Kimura, "Valuing finite-lived Russian options," European Journal of Operational Research, vol. 189, no. 2, pp. 363-374, 2008.

[10] H. Y. Wong and J. Zhao, "Valuing American options under the CEV model by Laplace-Carson transforms," Operations Research Letters, vol. 38, no. 5, pp. 474-481, 2010.

[11] C. S. Pun and H. Y. Wong, "CEV asymptotics of American options," Journal of Mathematical Analysis and Applications, vol. 403, no. 2, pp. 451-463, 2013.

[12] M. Leippold and N. Vasiljevic, "Pricing and Disentanglement of American Puts in the Hyper-Exponential Jump-Diffusion Model," SSRN Electronic Journal, 2015, http://ssrn.com/abstract =2571208.

[13] X. Lu and E. R. Putri, "Finite maturity margin call stock loans," Operations Research Letters, vol. 44, no. 1, pp. 12-18, 2016.

[14] J. Ma, Z. Zhou, and Z. Cui, "Hybrid Laplace transform and finite difference methods for pricing American options under complex models," Computers and Mathematics with Applications, vol. 74, no. 3, pp. 369-384, 2017.

[15] Z. Zhou, J. Ma, and H.-w. Sun, "Fast Laplace transform methods for free-boundary problems of fractional diffusion equations," Journal of Scientific Computing, vol. 74, no. 1, pp. 49-69, 2018.

[16] J. C. Cox, "The constant elasticity of variance option pricing model," The Journal of Portfolio Management, vol. 22, pp. 16-17, 1996.

[17] C. Chiarella and A. Ziogas, "Evaluation of American strangles," Journal of Economic Dynamics and Control, vol. 29, no. 1-2, pp. 31-62, 2005.

[18] X. Chen and J. Chadam, "A mathematical analysis of the optimal exercise boundary for American put options," SIAM Journal on Mathematical Analysis, vol. 38, no. 5, pp. 1613-1641, 2006/07.

[19] Y. K. Kwok, Mathematical Models of Financial Derivatives, Springer Verlag, Singapore, Singapore, 1998.

[20] M. Xu and C. Knessl, "On a free boundary problem for an American put option under the CEV process," Applied Mathematics Letters, vol. 24, no. 7, pp. 1191-1198, 2011.

[21] H. Stehfest, "Numerical inversion of Laplace transforms," Communications of the ACM, vol. 13, no. 1, pp. 47-49,1970.

[22] P. P. Valko and J. Abate, "Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion," Computers & Mathematics with Applications, vol. 48, no. 3-4, pp. 629-636,2004.

[23] A. Talbot, "The accurate numerical inversion of Laplace transforms," Journal of the Institute of Mathematics and Its Applications, vol. 23, no. 1, pp. 97-120,1979.

[24] H.-K. Pang and H.-W. Sun, "Fast numerical contour integral method for fractional diffusion equations," Journal of Scientific Computing, vol. 66, no. 1, pp. 41-66, 2016.

[25] J. A. Weideman and L. N. Trefethen, "Parabolic and hyperbolic contours for computing the Bromwich integral," Mathematics of Computation, vol. 76, no. 259, pp. 1341-1356, 2007.

[26] J. A. Weideman, "Improved contour integral methods for parabolic PDEs," IMA Journal of Numerical Analysis (IMAJNA), vol. 30, no. 1, pp. 334-350, 2010.

[27] J. Ma and Z. Zhou, "Convergence analysis of iterative Laplace transform methods for the coupled PDEs from regime-switching option pricing," Journal of Scientific Computing, vol. 75, no. 3, pp. 1656-1674,2018.

Zhiqiang Zhou (iD) and Hongying Wu

School of Mathematics and Finance, Xiangnan University, Chenzhou 423000, China

Correspondence should be addressed to Zhiqiang Zhou; zqzhou@2014.swufe.edu.cn

Received 24 March 2018; Accepted 13 August 2018; Published 4 September 2018

Academic Editor: Jorge E. Macias-Diaz

Caption: Figure 1: LTM-HCIM errors in log-scale v.s. L for American CEV Strangles option at [tau] = 1. Left: log scale errors with [[sigma].sub.0] = 0.2, [beta] = 0.5; Right: log scale errors with [[sigma].sub.0] = 0.4, [beta] = -0.5.

Table 1: Values of American Strangles at time t = 0 (i.e., [tau] = T) with [K.sub.1] = 1, [K.sub.2] = 1.5, r = 0.05, q = 0.1, [S.sub.0] = 1.5, [[sigma].sub.0] = 0.4, different values of [beta], 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 [beta] LTM-GSF LTM-HCIM LTM-GSF LTM-HCIM -0.50 0.260862 0.260864 0.297905 0.260864 -0.25 0.253815 0.253800 0.284163 0.253800 0 0.247986 0.247964 0.261926 0.247964 0.25 0.245120 0.245090 0.239312 0.245090 0.50 0.236749 0.236703 0.217691 0.236724 L = 26 [beta] LTM-GSF LTM-HCIM FDM -0.50 18.447801 0.260864 0.258805 -0.25 22.070192 0.253800 0.251154 0 24.667532 0.247964 0.244310 0.25 -0.098661 0.245090 0.238169 0.50 -9.087067 0.236724 0.232615 Table 2: Values of American Strangles at time t = 0 (i.e., [tau] = T) with [K.sub.1] = 1, [K.sub.2] = 1.5, r = 0.05, q = 0.1, [S.sub.0] = 1.5, L = 20, different values of [[sigma].sub.0], different values of [beta], and different space partition n. n Value [[sigma].sub.0] = 0.2, Rate Value [beta] = 0.5 Err 24 0.150660 3.3277e-003 0.6061 0.262272 49 0.149394 2.1862e-003 0.5367 0.260271 99 0.149262 1.5070e-003 0.5154 0.260335 199 0.149211 1.0543e-003 0.5076 0.260319 399 0.149187 7.4158e-004 0.5037 0.260274 799 0.149191 5.2302e-004 0.5018 0.260276 1599 0.149189 3.6937e-004 -- 0.260275 3199 0.149189 -- -- 0.260275 n [[sigma].sub.0] = 0.4, [beta] = -0.5 Err Rate 24 2.7302e-003 1.5515 49 9.3141e-004 0.9603 99 4.7871e-004 1.0611 199 2.2943e-004 0.9601 399 1.1793e-004 0.9812 799 5.9737e-005 0.9920 1599 3.0034e-005 -- 3199 -- --

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Zhou, Zhiqiang; Wu, Hongying |

Publication: | Discrete Dynamics in Nature and Society |

Date: | Jan 1, 2018 |

Words: | 6438 |

Previous Article: | Transient Analysis of Speed-Varying Rotor with Uncertainty Based on Interval Approaches. |

Next Article: | Investigating a Coupled Hybrid System of Nonlinear Fractional Differential Equations. |