Pricing American-style Parisian down-and-out call options

Applied Mathematics and Computation 305 (2017) 330–347

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Pricing American-style Parisian down-and-out call options Nhat-Tan Le a, Duy-Minh Dang b,∗ a

Department of Mathematics, International University, Vietnam National University, Quarter 6, Linh Trung Ward, Thu Duc District, Ho Chi Minh City, Viet Nam School of Mathematics and Physics University of Queensland, St Lucia, Brisbane 4072, Australia

b

a r t i c l e

i n f o

Keywords: Integral equation approach American-style Parisian options Fourier sine transform Coupled integral equations

a b s t r a c t We propose an integral equation approach for pricing American-style Parisian down-andout call options under the Black–Scholes framework. For this type of options, the knockout feature is activated only if the underlying asset price continuously remains below a pre-determined barrier for a sufficiently long period of time. As such, the corresponding pricing problem becomes a three-dimensional (3-D) free boundary problem, instead of a two-dimensional (2-D) one as is the case of “one-touch” barrier options, and this poses a computational challenge. In our approach , we first reduce the 3-D problem to a 2-D one, and then, by applying the Fourier sine transform to the resulting 2-D problem, we can derive a pair of coupled integral equations governing the option price at any given time in terms of (i) the option price at the barrier and (ii) the optimal exercise boundary at that time. This pair of coupled integral equations can be solved using the Newton–Raphson iterative procedure, after which, the option price, the optimal exercise boundary, and the hedging parameters can be obtained in a straightforward manner. A complexity analysis of the method, together with numerical results, show that the proposed approach is robust and significantly more efficient than existing uniform finite difference methods with Crank–Nicolson timestepping, especially in dealing with spot prices near the barrier. Numerical results are also examined in order to provide new insight into several interesting properties of the option price and the optimal exercise boundary. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Barrier options are very commonly traded path-dependent options in over-the-counter financial markets, as they can be viewed as inexpensive alternatives to vanilla options [1]. However, “one-touch” knock-out options may have an undesirable feature of suddenly losing the embedded options, even if the price of the underlying asset just momentarily breaches the (pre-determined) constant barrier, regardless of how brief the breaching is. Such barrier options are prone to market manipulations, as market participants can deliberately manipulate the underlying asset price to force the cancellation of the option to occur. To partially prevent such market manipulations, Parisian options were first introduced by Chesney et al. [2], with an extended “trigger condition” requiring that the underlying asset must remain continuously below (down-and-out) or above (up-and-out) the barrier for a sufficiently long period of time, referred to as the “option window”, before the knock-out

∗

Corresponding author. E-mail addresses: [email protected] (N.-T. Le), [email protected] (D.-M. Dang).

http://dx.doi.org/10.1016/j.amc.2017.02.015 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

331

feature is activated. Such a requirement certainly makes the market fairer in the sense that it protects the holder of Parisian options from deliberate action taken by the writer. This extended trigger condition, i.e. the “option window”, can also be found in some derivative contracts, such as callable convertible bonds and executive warrants [3]. The “option window” can also be useful in studying an optimal decision to invest in a project when delays are involved [4]. Although Parisian options have been rarely traded after the 2008 financial crisis, they can still be a useful tool in corporate finance [5], credit risk and life insurance [6,7]. Compared to “one-touch” knock-out options, the pricing of Parisian knock-out options is much more challenging, primarily due to the fact that the option prices depend not only on the elapsed time and the underlying asset price, but also on the “excursion time”, i.e. the amount of time the underlying asset price stays continuously above or below the pre-determined barrier [8]. In particular, for American-style Parisian knock-out options, the optimal exercise boundaries also need to be solved simultaneously with the option prices. As a result, the pricing of this type of options is now a 3-D free-boundary problem, which is very challenging to deal with. There have been several pricing approaches for American-style Parisian knock-out options. The partial differential equation (PDE) approach with finite differences is proposed in Haber et al. [9]. More specifically, the authors establish and solve by finite differences a pair of coupled PDE systems governing the price of an American-style Parisian knock-out option. While this method is flexible and easy to implement, there are some deficiencies in the proposed pricing systems, which are pointed out in [10]. Alternatively, in Chesney and Gauthier [11], a probabilistic method to price American-style Parisian knock-out options is proposed. More precisely, the pricing problem is first reduced to the problem of finding the Laplace transform of the distribution of the “Parisian stopping time”. As a result, the option price can be computed by numerically inverting the Laplace transform. However, numerical inversion of Laplace transform may be unstable and sensitive to round-off errors [12,13], and this is the main drawback of this approach. These drawbacks of existing numerical methods for pricing American-style Parisian knock-out options form the motivation for our work. In this paper, we present an integral equation method for pricing American-style Parisian down-and-out call options under the Black–Scholes framework. It turns out that the dimension reduction technique developed in [10] for European Parisian options can be adapted to reduce the number of dimensions of the American Parisian pricing problem from 3-D to 2-D. The major contributions of the paper are mentioned below. • By means of the Fourier sine transform, we convert the resulting 2-D problem into a pair of coupled integral equations that governs the option price at any given time in terms of (i) the option price at the barrier, and (ii) the optimal exercise boundary at that time. • We develop a Newton–Raphson iterative procedure to solve this pair of coupled equations efficiently. After this step the option price, the optimal exercise boundary, and the hedging parameters can be obtained in a straightforward fashion. • We present a complexity analysis of the proposed method. This, together with numerical results, show that the proposed approach is robust and significantly more efficient than existing (uniform) finite difference methods with Crank–Nicolson timestepping, especially in dealing with spot prices near the barrier. The paper is organized as follows. In Section 2, we introduce the PDE systems governing the price of an American-style Parisian down-and-out call option. The solution procedure is presented in Section 3, and the numerical implementation is discussed in Section 4. Section 5 presents selected numerical results to illustrate the significantly higher level of efficiency of the proposed method compared to existing (uniform) finite difference methods with Crank–Nicolson timestepping. Furthermore, through these examples, we demonstrate interesting properties of the price of an American-style Parisian downand-out call, as well as of its optimal exercise boundary. Concluding remarks and future work are given in Section 6. 2. Formulation We assume that the underlying asset price, denoted by S, follows a geometric Brownian motion given by:

dS(t ) = (r − δ )dt + σ dZ. S(t )

(2.1)

Here, r and σ denote the risk-free interest rate and the instantaneous volatility, respectively; δ is a positive constant continuous dividend yield; Z is a standard one-dimensional Brownian motion. We are interested in the valuation problem of an American-style Parisian down-and-out call option written on S, with maturity T and strike E. The knock-out barrier is specified by the pre-determined constant S¯. Also, we denote by J the excursion time, which measures the total time the underlying asset has spent continuously below the barrier S¯. Since the excursion time variable J starts accumulating from zero at the same rate of the elapsed time when S < S¯, it is straightforward to see that



J = 0, dJ = 0 d J = dt

if S ≥ S¯, if S < S¯,

(2.2)

Note that for a Parisian down-and-out option, the knock-out feature is triggered only when S continuously stays below S¯ for a pre-determined period of time, denoted by J¯, usually referred to as the “option window”. One can observe that if J¯ takes on an extreme value, namely zero or a value not smaller than the option maturity T, then an American-style Parisian down-and-out call option degenerates to a much simpler option. Specifically, when J¯ = 0, the option becomes the

332

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

“one-touch” down-and-out counterpart. When J¯ > T , the option becomes its embedded vanilla American option, since it can never be knocked out during its life. In the rest of the paper, we are interested in non-degenerate cases, and hence we consider 0 < J¯ ≤ T . We now make the usual assumption that E > S¯ for a down-and-out call option, because the holder often accepts the loss of his/her option only when the option is out-of-money. From this assumption, we deduce that the optimal exercise price is always greater than S¯ as the option should be optimally exercised only when it is in-the-money. The optimal exercise is therefore a function depending on time only. We denote by Sf (t) the optimal exercise boundary of the option at time t. We note that Sf (t) is not known in advance, and must be solved simultaneously with the option price. We now investigate the pricing domain of the problem. When S < S¯, the excursion time variable J starts to accumulate and the option price depends on three variables, namely S, t and J. The pricing domain in this case, denoted by D1 , can be expressed as

D1 = {(S, t, J )|0 ≤ S ≤ S¯, 0 ≤ t ≤ T , 0 ≤ J ≤ J¯}. When S ≥ S¯, the excursion time variable J is reset to zero and remains zero until S hits S¯ again. Consequently, in this case, the option price depends on two variables, namely S and t. In addition, at any time t, if S ≥ Sf (t) then the option should be exercised immediately, in which case, the option value is equal to the payoff. As a result, we only need to price the option when S¯ ≤ S ≤ S f (t ). Consequently, the pricing domain in this case, denoted by D2 , can be expressed as

D2 = {(S, t )|S¯ ≤ S ≤ S f (t ), 0 ≤ t ≤ T }. It turns out that D1 and D2 can be further simplified. More specifically, in domain D1 , the sub-domain corresponding to t < J, denoted by D3

D3 = {(S, t, J )|0 ≤ S ≤ S¯, 0 ≤ t ≤ J, 0 ≤ J, ≤ J¯}, see Fig. C.1, can be cut off, because in this region, the elapsed time t is less than the excursion time J, which contradicts (2.2). Furthermore, in the sub-domain corresponding to t > T − J¯ + J, denoted by D4 ,

D4 = {(S, t, J )|0 ≤ S ≤ S f (t ), t > T − J¯ + J, 0 ≤ J ≤ J¯}, there is not enough time left for J to reach J¯, and hence, the option can never be knocked out, i.e. it degenerates to its embedded vanilla American option. As such, sub-domain D4 can also be removed from D1 and D2 . The pricing domains of those non-degenerate cases can now be elegantly reduced to (Fig. C.1).





I = (S, t, J )|0 ≤ S ≤ S¯, J ≤ t ≤ T ∗ + J, 0 ≤ J ≤ J¯ ,





II = (S, t )|S¯ ≤ S ≤ S f (t ), 0 ≤ t ≤ T ∗ , where T ∗ = T − J¯.

(2.3)

Denote by V1 (S, t, J) and V2 (S, t) the price of an American-style Parisian down-and-out call option in pricing domains I and II, respectively. We now establish a pair of coupled PDE systems governing V1 (S, t, J) and V2 (S, t). In domain I, it can be shown that V1 (S, t, J) is governed by a modified Black–Scholes equation [9]:

∂ V1 ∂ V1 + + LV1 = 0, ∂t ∂J

(2.4)

where

LV =

σ 2 S2 ∂ 2V ∂V + ( r − δ )S − rV, 2 ∂ S2 ∂S

(2.5)

with the following boundary conditions:

V1 (0, t, J ) = 0,

(2.6)

lim V1 (S, t, J ) = 0,

(2.7)

V1 (S¯, t, J ) = V1 (S¯, t, 0 ).

(2.8)

J→J¯

Condition (2.6) is standard for a call option when the underlying price is zero. Condition (2.7) follows from the fact that when J reaches J¯, the option is knocked out and becomes worthless. Condition (2.8) is the so-called “reset condition”which indicates that J is reset to zero every time S hits S¯ from below. To ensure the no arbitrage opportunity, we require that the option Delta be continuous across the barrier S¯ (from domain I to domain II), i.e.

∂ V1 ¯ ∂V (S, t, 0 ) = 2 (S¯, t ). ∂S ∂S Condition (2.9) is hereinafter referred to as the Delta condition.

(2.9)

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

333

In domain II, V2 (S, t) satisfies the classical Black–Scholes equation

∂ V2 + LV2 = 0, ∂t

(2.10)

where the operator L is defined in (2.5). We recall that, as noted earlier, in t > T ∗ + J, the American-style Parisian call option is identical to its embedded vanilla American call, as there is not enough time left for J to reach J¯. As such, the terminal condition for V2 is given by

V2 (S, T ∗ ) = CA (S, T ∗ ),

(2.11)

T∗ )

where CA (S, denotes the price of the embedded vanilla American call option at the underlying asset price S and time t = T ∗ , where T ∗ = T − J¯ as defined in (2.3). Besides the terminal condition, a set of boundary conditions along the S direction is also needed to solve for V2 . The continuity of the option price across the barrier S¯ requires

V2 (S¯, t ) = V1 (S¯, t, 0 ).

(2.12)

Also, the two necessary conditions for determining the optimal exercise boundary of an American-style Parisian down-andout call option are:

V2 (S f (t ), t ) = S f (t ) − E,

∂ V2 (S (t ), t ) = 1. ∂S f

(2.13)

Eqs. (2.4)–(2.13) constitute a pair of coupled PDE systems governing the value of an American-style Parisian down-andout call option at any underlying price S, any excursion time J, and any time t before the expiration T. The coupled PDE systems can be summarized in (2.14)–(2.15):

⎧ ∂ V1 ∂ V1 ⎪ ⎪ + + LV1 = 0, ⎪ ⎪ ∂t ∂J ⎪ ⎨lim V (S, t, J ) = 0, 1 A1 J→J¯ and ⎪V (0, t, J ) ⎪ = 0, 1 ⎪ ⎪ ⎪ ⎩ V1 (S¯, t, J ) = V1 (S¯, t, 0 ), Delta condition :

⎧ ∂V 2 ⎪ ⎪ ∂ t + LV2 ⎪ ⎪ ⎪ ∗ ⎪ ⎪ ⎨V2 (S, T )

= 0, = CA (S, T ∗ ),

A2 V2 (S f (t ), t )

= S f (t ) − E, ⎪ ⎪ ⎪ ⎪ ⎪ ∂ V2 (S f (t ), t ) = 1, ⎪ ⎪ ⎩ ∂S V2 (S¯, t ) = V1 (S¯, t, 0 ).

∂ V1 ¯ ∂V (S, t, 0 ) = 2 (S¯, t ), ∂S ∂S

(2.14)

(2.15)

where A1 and A2 are defined on the regions I and II, respectively. 3. Solution procedure 3.1. Transformation to a 2-D PDE system By taking advantage of the shape of the pricing domain, the above 3-D systems (2.14), (2.15) can be further simplified to a pair of 2-D systems. We adapt the transformation technique in [10,14]. More specifically, to transform (2.14), (2.15) to a 2-D system, we only need to deal with the system governing V1 , because the system governing V2 is already in 2-D. The idea is √ ∂V ∂ V1 ∂ V1 to replace the sum of the two partial derivatives appearing in A1 , namely + , by the directional derivative 2 1 , ∂t ∂J √ √ ∂v which represents the instantaneous rate of change of the function V1 at the point (t, J) in the direction of ( 2/2, 2/2 ). √ v Furthermore, the rescaling u = √ can be applied to absorb the constant 2. After these steps, the PDE systems (2.14), 2 (2.15) become

⎧ ∂ V1 ⎪ ⎪ = 0, ⎪ ∂ u + LV1 ⎪ ⎪ ⎨lim V (S, u; t ) = 0, 1 B1 u→J¯ and ⎪V (0, u; t ) ⎪ = 0, 1 ⎪ ⎪ ⎪ ⎩ ¯ V1 (S, u; t ) = W (t + u ), Delta condition :

⎧ ∂V 2 ⎪ + LV2 ⎪ ⎪ ∂ t ⎪ ⎪ ∗ ⎪ ⎪ ⎨V2 (S, T )

= 0, = CA (S, T ∗ ),

B2 V2 (S f (t ), t )

= S f (t ) − E, ⎪ ⎪ ⎪ ∂ V 2 ⎪ ⎪ (S f (t ), t ) = 1, ⎪ ⎪ ⎩ ∂S ¯ V2 (S, t ) = W (t ),

∂ V1 ¯ ∂V (S, 0; t ) = 2 (S¯, t ). ∂S ∂S

(3.1)

334

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

The systems B1 and B2 are defined on:







B1 = (S, t, u ) 0 ≤ S ≤ S¯, 0 ≤ t ≤ T ∗ , 0 ≤ u ≤ J¯ ,







B2 = (S, t ) S¯ ≤ S ≤ S f (t ), 0 ≤ t ≤ T ∗ , respectively. Here,

W (t ) = V1 (S¯, 0; t ),

(3.2)

T∗ ),

and CA (S, as introduced in (2.11), is the price of the embedded vanilla American call at the underlying price S and time t = T ∗ , where T ∗ = T − J¯. 3.2. Non-dimensionalization In order to make the newly-derived 2-D PDE systems simpler, we shall now non-dimensionalize the PDE systems by introducing dimensionless variables, as suggested in [1,15]

S x = ln , S¯

τ=

σ2 2

( T ∗ − t ),

l = (J¯ − u )

σ2 2

;

constants

σ2

L = J¯

2

,

γ=

2r

σ2

,

q=

k 2δ , k = γ − q − 1, α = − , β = −α 2 − γ ; 2 σ2

and unknown functions C1 (x, l; τ ), C2 (x, τ ), U(τ ) and xf (l; τ ) defined by:

C1 (x, l; τ ) = S¯−1 e−α x−β l V1 (S, u; t ), C2 (x, τ ) = S¯−1 e−α x−βτ V2 (S, t ),

U (τ ) = S¯−1W (t, ) x f (τ ) =

ln S f (t ) . S¯

Under this change of variables, the systems become dimensionless, and are given by:

B3

⎧ ∂ C1 ⎪ ⎪ (x, l; τ ) ⎪ ⎪ ⎪ ⎨C∂(l x, 0; τ ) 1

∂ C1 = (x, l; τ ), ∂ x2 2

= 0,

⎪ lim C1 (x, l; τ ) = 0, ⎪ ⎪ x→−∞ ⎪ ⎪ ⎩ C1 (0, l; τ ) = e −β l U ( τ + l − L ) , Delta condition : eβ L

and

⎧ ∂C ⎪ ⎪ 2 (x, τ ) ⎪ ⎪ ∂τ ⎪ ⎪ ⎪ C2 (x, 0 ) ⎪ ⎨ B4 C2 (0, τ ) ⎪ ⎪ ⎪ C2 (x f (τ ), τ ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ C2 (x (τ ), τ ) ∂x f

∂ 2C2 (x, τ ), ∂ x2 = g( x ) ,

=

= g1 ( τ ),

(3.3)

= g2 ( x f ( τ ), τ ), = g3 ( x f ( τ ), τ ),

∂ C1 ∂C (0, L; τ ) = eβτ 2 (0, τ ). ∂x ∂x

(3.4)

Here, the systems B3 and B4 are defined on



σ2 ∗ B3 = (x, τ , l ) − ∞ < x ≤ 0, 0 ≤ τ ≤ T , 0≤l≤L , 2 

σ2 ∗ B4 = (x, τ ) 0 ≤ x ≤ x f (τ ), 0 ≤ τ ≤ T , , 2

and U (τ ) = eβ LC1 (0, L; τ ). Datum g(x), g1 (τ ), g2 (xf (τ ), τ ) and g3 (xf (τ ), τ ) are given by:

g(x ) = S¯−1 e−α xCA (S¯ex , T ∗ ), g1 (τ ) = e−βτ U (τ ),

(3.5)

E g2 (x f (τ ), τ ) = e(1−α )x f (τ )−βτ − e−α x f (τ )−βτ , S¯ E g3 (x f (τ ), τ ) = (1 − α )e(1−α )x f (τ )−βτ + α e−α x f (τ )−βτ . S¯ Note that in (3.3), (3.4), C1 and C2 are still coupled via U(τ ) and xf (τ ), hence they can not be solved directly. The underlying idea of our solution procedure is as follows. • First, we find integral representations for C1 and C2 in terms of the unknown functions U(τ ) and xf (τ ). This can be achieved by solving C1 and C2 from (3.3), (3.4) separately, as if they were not coupled. We then revert change of variables used previously in Section 3.2 to find the corresponding integral representations for V1 (S, u; t) and V2 (S, t) in terms of the unknown function W(t), defined in (3.2) and the unknown optimal exercise boundary Sf (t).

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

335

• The next step is to find a pair of coupled integral equations that governs Sf (t) and W(t), and solve for Sf (t) and W(t) at discrete time points. • Finally, use these discrete values of Sf (t) and W(t) to approximate V1 (S, u; t) and V2 (S, t). In the remainder of the paper, we discuss these steps in detail. 3.3. Integral representations of V1 (S, u; t) and V2 (S, t) The system B3 is clearly a classical heat problem in a semi-finite domain and has been well-studied in a number of text books [16–19]. Using therein results to solve the system B3 to get C1 , together with reverting change of variables to the original variables S and t, we obtain the following integral representation for V1 (S, u; t) in terms of the unknown function W(t), defined in (3.2):

V1 (S, u; t ) =

(ln S¯ − ln S )  S α √ S¯ σ 2π

J¯ u

S¯ )2 W (t + v ) β σ22 (v−u )− (2lnσS−ln 2 (v−u ) e d v. ( v − u )3



(3.6)

The system B4 is a heat problem in a finite time-dependent domain, and hence, it is much more complicated to solve. This system, however, resembles the governing system of the price of an American down-and-out call, which is discussed in Le et al. [1] via the Fourier sine transform. We will therefore adopt this approach to solve for C2 in the system B4 . To this end, we define by H(x) the Heaviside function

H (x ) =

⎧ ⎪ ⎨1

if x > 0,

1/2 if x = 0,

⎪ ⎩

(3.7)

if x < 0.

0

As shown in Appendix A, an integral representation of the solution C2 of the system B4 can be shown to be





x f (0 )

(x−u )2 ( x+u )2 g( u ) e− 4τ − e− 4τ du √ 2 πτ 0

τ x −(S¯e − E )1x=x f (τ ) (x ) + M(x, τ ) + Q(x, τ , ξ , x f (ξ ))dξ ,

S¯eα x+βτ H (x f (τ ) − x )C2 (x, τ ) = S¯eα x+βτ

(3.8)

0

where

1 1x=x f ( τ ) ( x ) =

if x = x f (τ ),

2 0

if x = x f (τ ),

and

 M(x, τ ) = S¯ex−qτ N

 −e

2α x

x − x f ( 0 ) + 2 ( 1 − α )τ √ 2τ



S¯e−qτ −x N

 Q(x, l, ξ , x f (ξ )) = qS¯ex−q(τ −ξ ) N

− Ee





− e2α x qS¯e−q(τ −ξ )−x N

 −γ ( τ −ξ )

x − x f (0 ) − 2ατ N √ 2τ



− Ee

N

−x − x f (0 ) − 2ατ N √ 2τ



2 (τ − ξ )

 − E γ e −γ ( τ −ξ ) N

−x − x f (ξ ) + 2(1 − α )(τ − ξ )







x − x f ( ξ ) − 2α ( τ − ξ )



(3.9)



2 (τ − ξ )

2 (τ − ξ )

−x − x f (ξ ) − 2α (τ − ξ )







−γ τ

x − x f (ξ ) + 2(1 − α )(τ − ξ )



+



−γ τ

−x − x f (0 ) + 2(1 − α )τ √ 2τ

and

− Eγ e





2 (τ − ξ )

x2 xS¯U (ξ ) eβ (τ −ξ )+α x− 4(τ −ξ ) .  2 π ( τ − ξ )3

(3.10)

Transforming (3.8), (3.9), (3.10) back to the original variables S and t, we obtain the following formula for V2 (S, t) in terms of Sf (t) and the unknown function W(t), defined in (3.2),





H (S f (t ) − S )V2 (S, t ) = −(S − E )1S=S f (t ) (S ) + R(S, t ) + M S, T ∗ − t, S f (T ∗ ) +

t

T∗

Q (S, t, v, S f (v ))dv,

(3.11)

336

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

where

  x α eβ σ22 (T ∗ −y)C (v, T ∗ )  (ln x−ln v)2 (ln x+ln v−2 ln S¯ )2 A e 2σ 2 (T ∗ −y) − e 2σ 2 (T ∗ −y) d v,  v S¯ vσ 2π ( T ∗ − y )   x 2α  S¯2 M (x, y, z ) = M1 (x, y, z ) − M1 , y, z , x S¯   α  x 2α  S¯2 x Q (x, y, z, w ) = Q1 (x, y, z, w ) − Q1 , y, z, w + K (x, y, z ). R(x, y ) =

S f (0 )

S¯

(3.12)

S¯

x

and T ∗ = T − J¯ as defined in (2.3). Here, M1 , Q1 and K are defined by

M1 (x, y, z ) = xe−δ y N (d1 (x, y, z )) − Ee−ry N (d2 (x, y, z ))

Q1 (x, y, z, w ) = xδ e−δ (z−y ) N (d1 (x, z − y, w )) − Ere−r (z−y ) N (d2 (x, z − y, w )), 2 S¯ )2 (ln x − ln S¯ )W (z ) − (2lnσx2−ln +β σ2 (z−y ) (z−y ) K (x, y, z ) = e ,  σ 2π ( z − y )3

(3.13)

with

d1 (x, y, z ) = and

ln x − ln z + (r − δ + σ 2 /2 )y , √ σ y

 IS=S f (t ) (S ) =

1 2 0

if

S = S f (t ),

if

S = S f (t ).

d2 (x, y, z ) = d1 (x, y, z ) − σ

√ y.

(3.14)

(3.15)

3.4. Coupled integral equations for Sf (t) and W(t) From (3.6)–(3.11), we have obtained integral representations of V1 (S, u; t) and V2 (S, t) in terms of Sf (t) and W(t). The next step is to derive a pair of coupled integral equations governing Sf (t) and W(t). The first one can be easily derived by substituting S = S f (t ) into (3.11). After some simplifications, we obtain the first integral equation governing the values of Sf (t) and W(t):





S f (t ) − E = M S f (t ), T ∗ − t, S f (T ∗ ) + R(S f (t ), t ) +

T∗

t

Q (S f (t ), t , v, S f (v ))dv,

(3.16)

T ∗ = T − J¯ as defined in (2.3), and M(·), R(·), and Q(·) are defined in (3.12). Note that W(·) is hidden in K(·) defined in (3.13). The second integral equation governing the values of Sf (t) and W(t) can be derived from the Delta condition (3.1). To ∂ V1 ∂V this end, we first need to calculate the Delta () of the option when S is below or above S¯, i.e. (S, u; t ) and 2 (S, t ), ∂S ∂S respectively. ∂V

Proposition 3.1. The Delta of the option when S < S¯, i.e. ∂ S1 (S, u; t ), is given by

√

 S )2 ∂ 2W (t + J¯) S α β σ22 (J¯−u )− (lnσ 2S¯−ln (J¯−u ) V1 (S, u; t ) = e  ¯ ∂S Sσ π (J¯ − u ) S ⎧ ⎪ ⎨

√

⎩ 2α Sα −1 +∞ ⎪ + √ e π S¯α 0

− 12

 S α √2σ √J¯−u

− Proof. See Appendix B

S¯

√

πS

0



e

η

−

S¯−ln S + ln√ σ J¯−u

(ln S¯−ln S )2 2σ 2 v2

2

2 ¯ + β (ln S−ln S ) 2 S¯−ln S 2 η+ ln√ σ J¯−u

β

+ 2 σ 2 v2

!

⎫ ⎪ ⎬ ⎪ ⎭

⎛

⎞ ( ln S¯ − ln S )2 ⎠dη W ⎝t + u + S¯−ln S 2 σ 2 (η + ln√ ) ¯ σ

βW (t + u + v2 ) +

2

σ

J−u

"

W (t + u + v2 ) dv. 2

(3.17)



From the formula (3.17), we can obtain

∂ V1 ¯ α (S, 0; t ) = W (t ) + ∂S S¯ ∂V



∂ V1 ¯ (S, 0; t ) straightforwardly as follows: ∂S

√ √ √J¯ ! " 2W (t + J¯) β σ 2 J¯ σ 2 2 σ2 2 e 2 − √ eβ 2 u β W (t + u2 ) + 2 W (t + u2 ) du.  σ S¯ π 0 S¯σ π J¯

(3.18)

The calculation of ∂ S2 (S, t ), on the other hand, is very similar to that of the Delta of the American down-and-out call options, which was presented in the appendix A4 in [1]. The main difference lies on the appearance of the term R(S, t)

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

337

in the formula of V2 (S, t). Fortunately, its derivative can be easily calculated simply using the chain rule. Therefore, the following proposition can be obtained. ∂V

Proposition 3.2. The Delta of the option when S ≥ S¯, i.e. ∂ S2 (S, t ), is given by

∂ ˜ (S, T ∗ − t, S f (T ∗ )) + V (S, t ) = R˜(S, t ) + K˜1 (S, t ) + M ∂S 2 where,

T∗

t

∀S < S f (t ),

L(S, t, v, S f (v ))dv,

(3.19)

  ln x − ln u du α− 2 ∗ σ (T − y ) S¯ uα +1 σ 2π (T ∗ − y )  

S f (0 ) σ2 ∗ ln S¯ )2 ln x + ln u − 2 ln S¯ CA (u, T ∗ )xα −1 eβ 2 (T −y ) − (ln x2+σln2 (uT−2 ∗ −y ) − e du, α−  σ 2 (T ∗ − y ) S¯ uα +1 σ 2π (T ∗ − y ) ⎛ ⎛ ⎞2 ⎞ β (ln x−ln S¯)2 S¯ 2 √

− 12 (η+ ln√x−ln ) x−ln S¯ )2 σ T ∗ −y 2(η+ ln√ 2α xα −1 +∞ ⎜ ln x − ln S¯ ⎟ ∗ −y σ T ⎠ ⎠e K˜1 (x, y ) = √ W ⎝y + ⎝ dη, S¯ π S¯α 0 σ (η + σln√x−ln ) T ∗ −y R˜(x, y ) =

+ +

S f (0 )

CA (u, T ∗ )xα −1 eβ



 x α √2σ √T ∗ −y S¯ xσ

√

πx 0 √ 2W (T ∗ ) 

λxλ−1 S¯λ

− (ln2x−ln∗ u )

+ β2 σ 2 v2

! " 2 β W ( y + v2 ) + 2 W ( y + v2 ) d v σ ¯ 2

S) β σ2 (T ∗ −y )− (σln2x(−ln T ∗ −y ) 2

e





Q1

2

2σ (T −y )

S¯2 , y, z x

 +

S¯2 , y, z, w x



 x 2α−2 S¯ +

 ˜1 M

 x λ−2 S¯

 Q˜1



S¯2 , y, z , x



S¯2 , y, z, w , x



 ˜ 1 (x, y, z ) = e−δ y N (d1 (x, y, z )) + M

 Q˜1 (x, y, z, w ) = e−δ (z−y )

e

(ln x−ln S¯ )2 2σ 2 v2

2α x2α −1 M1 S¯2α

L(x, y, z, w ) = Q˜1 (x, y, z, w ) − Here

e−

 x α

π (T ∗ − y ) S¯

˜ (x, y, z ) = M ˜ 1 (x, y, z ) − M

σ 2 (T ∗ −y ) 2

˜ (d1 (x, y, z ))  E N 1− , √ σ y z

δ N (d1 (x, z − y, w )) +



˜ (d1 (x, z − y, w ) )  Er N δ− , √ σ z−y w

with

1 − x2 ˜ (x ) = √ N e 2 , and d1 (x, y, z ), d2 (x, y, z ) are defined in 3.14. 2π Based on Proposition 3.2, we now can calculate

∂ V2 ¯ (S, t ) by ∂S

∂ ˜ (S¯, T ∗ − t, S f (T ∗ )) + V (S¯, t ) = R˜(S¯, t ) + K˜1 (S¯, t ) + M ∂S 2

t

T∗

L(S¯, t, v, S f (v ))dv,

(3.20)

˜ , L are defined as in Proposition 3.2, and T ∗ = T − J¯ as defined in (2.3). As a result, from (3.18) and (3.20), where R˜, K˜1 , M using the Delta condition (3.1), we can derive the second integral equation governing the values of Sf (t) and W(t)

α S¯

√ √ √J¯ ! " 2W (t + J¯) β σ 2 J¯ σ 2 2 σ2 2 e 2 − √ eβ 2 u β W (t + u2 ) + 2 W (t + u2 ) du  σ S¯ π 0 S¯σ π J¯

T∗ ˜ (S¯, T ∗ − t, S f (T ∗ )) + = R˜(S¯, t ) + K˜1 (S¯, t ) + M L(S¯, t, v, S f (v ))dv.

W (t ) +

(3.21)

t

It should be noted that once Sf (t) and W(t) are approximated by numerically solving the coupled integral Eqs. (3.16) and (3.21) at discrete time points, they can be used to approximate the value of an American-style Parisian down-and-out call when S is below or above S¯ via formulas (3.6) and (3.11), respectively. Similarly, the hedging parameters Delta, Gamma, Theta, Vega and Rho, typically referred to as the Greeks, can the also be approximated by differentiating the formulas (3.6) and (3.11) with respect to the relevant parameters. For instance, the value of Delta associated with S below or above S¯ can be easily approximated by using formulas (3.17) and (3.19), respectively. Therefore, it is very important to develop an efficient numerical scheme for solving the coupled integral Eqs. (3.16)–(3.21) at discrete time points. This is the focus of the next section.

338

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

4. Numerical approximations of V1 (S, ·; 0) and V2 (S, 0) 4.1. Numerical procedure We now describe a numerical procedure to approximate W(t) and Sf (t) at discrete time points. These values are then used to compute V1 (S, ·; 0) (for J > 0) and V2 (S, 0) (for J = 0) via formulas (3.6) and (3.11), respectively. p Without loss of generality, we can assume that T ∗ = J¯, where p and q are positive integer numbers. With this assumpq tion, we can construct the following uniform partition for the interval [0, T∗ ]

{T ∗ = t1 > t2 > · · · > tm > tm+1 = 0},

(4.1)

J¯ where t j = T ∗ − ( j − 1 )h with h = , j = 1, . . . , (m + 1 ), n = qz, and m = pz for some positive integer number z. n We aim to approximate W(tj ) and Sf (tj ) at each j ∈ {1, 2, . . . , m + 1}, starting from j = 1. We note that, at j = 1, W (T ∗ ) = V1 (S¯, 0; T ∗ ), as defined in (3.2), which is the same as CA (S¯, T ∗ ) the price of the embedded vanilla American option at time T∗ . This quantity can be approximated using a numerical method. To find Sf (T∗ ), we make use of the following corollary. Corollary 4.1. Denote by SVf (t ), 0 ≤ t ≤ T, the time-t optimal exercise boundary of the embedded vanilla American call option.1 Then

S f (T ∗ ) = SVf (T ∗ ). Proof. See Appendix C. The quantity

SVf (T ∗ )

 can be computed using a numerical method, and by Corollary 4.1, we also have Sf (T∗ ). In our exper-

iments, we use the integral equation methods of [20,21] to compute both CA (S, T∗ ) and SVf (T ∗ ). We now focus on a numerical procedure to approximate W(tj ) and Sf (tj ) at each discrete time point in (4.1). At each time t j+1 , j = 1, . . . , m, given previously computed W(tj ), we compute approximations to S f (t j+1 ) and W (t j+1 ) using the Newton

method. To this end, we denote by S(fk ) (t j ) and W(k) (tj ), k = 0, . . . , approximations to Sf (tj ) and W(tj ) at the kth Newton iteration, respectively. Once all approximations W(tj ) and Sf (tj ), j = 1, . . . , (m + 1 ), are determined, they can be used to approximate V1 (S, ·; 0) and V2 (S, 0) in a straightforward fashion, using the formulas (3.6) and (3.11), respectively. Similarly, we can also approximate the Delta of the option via formulas (3.17) and (3.19). In Algorithm 4.1, we present the computational

Algorithm 4.1 Algorithm to approximate V1 (S, ·; 0) and V2 (S, 0) using m time steps for the time period [0, T∗ ], given W (T ∗ ) = CA (S, T ∗ ). set W (t1 ) = W (T ∗ ); for j = 1, 2, . . . , m + 1 do 3: set W (0 ) (t j+1 ) = W (t j ); for k = 0, 1, 2 . . . , do 4: 5: substitute W (k ) (t j+1 ) into (3.16) and solve for S(fk ) (t j+1 ) using the Newton–Raphson scheme; 1:

2:

6:

(k ) (k ) S f (t j+1 )−S(fk−1) (t j+1 ) (k−1 ) (t j+1 ) < tol and W (t j+1()k−W if k > 0 and < tol then (k ) ) W (t j+1 ) S (t j+1 ) f

7: 8: 9: 10: 11: 12: 13:

break from the iteration; end if substitute S(fk ) (t j+1 ) into (3.21) and solve for W (k+1 ) (t j+1 ) using the Newton-Raphson scheme; end for set S f (t j+1 ) = S(fk+1 ) (t j+1 )and W (t j+1 ) = W (k+1 ) (t j+1 ); end for approximate V1 (S, ·; 0 ) and V2 (S, 0 )using formulas (3.6) and (3.11), respectively;

algorithm to compute V1 (S, ·; 0) and V2 (S, 0), given W (t1 ) ≡ W (T ∗ ) = CA (S¯, T ∗ ). We conclude this subsection by highlighting that, as clearly detailed in Algorithm 4.1, our approach does not involve discretizing either the S or the J variable as required by the finite difference approach, which is clearly an important computational advantage of our approach over the finite difference approach. More specifically, approximations to W(tj ) and Sf (tj ) at discrete time points for only J = 0 are first obtained, and then used in approximating V1 (S, ·; 0) and V2 (S, 0) via formulas (3.6) and (3.11), respectively. Furthermore, due to the Parisian feature, a grid-base approach typically does not handle well the knock-out feature for spot prices near the barrier, i.e. when S ↑ S¯, and J ↑ J¯. However, the proposed method does not have this issue. As a result, it is expected that the proposed integral equation method will be even more efficient than finite 1

We note that the underlying is a dividend-paying stock.

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

339

difference methods in this case. We will illustrate these advantages of the proposed integral equation method by numerical examples in Section 5. 4.2. Computational complexity We now discuss the computational complexity of the proposed integral equation method. We note that Algorithm 4.1 only covers the time period [0, T∗ ]. As input to this algorithm, as mentioned in the previous subsection, we need to compute both CA (S, T∗ ) and SVf (T ∗ ). We choose to compute these quantities via the integral equation methods

of [20,21] developed for American options. Specifically, these methods involve first approximating values of SVf (t ) at discrete time points in [T∗ , T] via the Newton method, then evaluating CA (S, T∗ ) using these values. At each timestep, and for each Newton iteration, the major cost involves evaluating the function value and its derivative. In this process, a total of two integrals need to be approximated, one for the function value and one for the derivative. The evaluation of CA (S, T∗ ) involves only one integral. The reader is referred to [20,21] for details of the computational steps. We now focus on the complexity of Algorithm 4.1, which covers the time period [0, T∗ ]. Similar to the period [T∗ , T], we need to evaluate a number of integrals. Unless otherwise stated, all integrals, including those arising in computing CA (S, T∗ ) as discussed above, are approximated by the composite Gauss–Legendre or Gauss–Laguerre rules [22]. Specifically, for integrals with smooth integrands on finite domains, the composite Gauss–Legendre quadrature rule is used. For singular % T∗ integrals, such as t K (S(fk ) (t j+1 ), t j+1 , v )dv which has a singularity at v = t j+1 , we first apply a suitable variable transform j+1

to convert them into a semi-infinite domain, and then apply the composite Gauss–Laguerre quadrature rule on the resulting integrals. We denote by g the number of points used for these rules. A breakdown of the major costs required by Algorithm 4.1 is as follows: 1. A construction of abscissae and weights for the g-point Gauss–Legendre and Laguerre quadrature rules entails a cost of approximately 2g2 (flops)2 . 2. Steps 2-12: these are for computing numerical approximations to S f (t j+1 ) and W (t j+1 ) at each of the m discrete time points. The computational cost of Steps 2–12 can be approximated by

cost-per-iteration × (total number of Newton iterations over [0, T ∗ ]) (flops). We now focus on each iteration. (a) At each iteration, to compute S(fk ) (t j+1 ) using the Newton–Raphson scheme (Step 5), the major computational cost involves computing the function value and its derivative. To compute the function value, examination of quantities appearing in Step 5 shows that a total of 4 integrals, namely one in R(S(fk ) (t j+1 ), t j+1 ) % T∗ and three in t Q (S(fk ) (t j+1 ), t j+1 , v, S f (v ))dv, need to be evaluated using the g-point Gauss quadrature rules. For j+1

each integral evaluation, there are approximately g additions, and g multiplications between the Gauss weights and the values of the integrand evaluated at the Gauss points. Furthermore, examination of the associated integrands show that the average cost for evaluating an integrand at a Gauss point is about 20 (flops). Thus, the cost for evaluating the function value is about

4(integrals) × 22g(flops/integral) = 88g (flops). For evaluating the derivative, 6 integrals need approximating, totalling 132g (flops). As a result, the cost per iteration for computing S(fk ) (t j+1 ) is approximately

88g + 132g = 220g (flops). (b) At each iteration, the computational cost for computing W (k+1 ) (t j+1 ) (Step 9) primarily involves evaluating a total of 6 integrals for the function value (this is achieved after some simplifications). For the function derivative, no integral needs approximating. Thus, the cost per iteration for computing W (k+1 ) (t j+1 ) is about

6(integrals) × 22g(flops/integral) = 132g (flops). As a result, the cost-per-iteration is approximately 220g + 132g = 352g (flops). 3. Step 13: approximately requires 88 mg (flops), taking into account that there are m time steps, and the cost per time step is approximately 88 g (flops) (4 integrals need to be evaluated, one for V1 (S, ·; 0) and three for V2 (S, 0), at the cost of about 22 g (flops) each). Thus, the computational cost of Algorithm 4.1 can be approximated by

2g2 + 352g × (total number of Newton iterations over [0, T ∗ ]) + 88mg(flops). 2

A flop is one addition, or one multiplication, or one division of two floating-point numbers.

340

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

ˆ be the number of time steps used for this Now we approximate the computational cost for the time period [T∗ , T]. Let m time period. First, we analyze the computational complexity in approximating SVf (t ) at discrete time points in [T∗ , T] via the Newton method. Similarly to the above analysis, the cost for this step can be approximated by

cost-per-iteration × (total number of Newton iterations over [T ∗ , T ]) (flops), where the cost-per-iteration is about

2(integrals) × 22g(flops/integral) = 44g (flops). (Recall that two integrals per iteration with the cost of approximately 22 g (flops) each.) The cost for evaluating CA (S, T∗ ) ˆ g (flops), taking into account that there is only one integral that needs to be evaluated at the cost of is approximately 22 m ˆ g (flops). about 22 m In summary, the total computational cost of the proposed method can be modelled by

total cost ≈ 2g2 + 352g × (total number of Newton iterations over [0, T ∗ ]) + 88mg ∗

ˆ g (flops). +44g × (total number of Newton iterations over [T , T ]) + 22m

(4.2) (4.2)

We conclude by highlighting that, as illustrated later in Section 5, the average number of iterations per time step required for achieving the stopping criterion for the Newton method with a tight tolerance of tol = 10−5 is relatively small, only 2–3 iterations per time step, and is independent of the time step size used, which is a desirable property. 5. Numerical examples and discussions In this section, we first provide selected examples to validate the proposed method and compare our numerical method with existing uniform finite difference methods in terms of efficiency, i.e. accuracy versus complexity. We then illustrate some interesting properties of the price and the optimal exercise boundary of an American-style Parisian down-and-out call option. The parameters for the tests are S¯ = {75, 95}, E = $100, T = 1(year ), J¯ = 1/15(years ), σ = 40%, r = 6%, δ = 10%. We denote by τ = T − t the time to maturity. We refer to the proposed integral equation method as the “IE method”. Unless otherwise stated, in all experiments for the IE method, we use the same uniform step of size 1/200 in the t-direction, and for the Newton method, we use tol = 10−5 . 5.1. Validation We now provide selected examples to validate the IE method. The benchmark solutions are obtained by solving the systems (2.14) directly using a uniform finite difference method with Crank–Nicolson timestepping and Rannacher smoothing [23], referred to as the “FD-CN” method, with very fine stepsizes. Here, the Rannacher smoothing involves running a fully implicit scheme for just a few time steps (3 in our experiments) near expiration, and then switching to Crank–Nicolson for all remaining time steps. It is well-established that the Rannacher smoothing can eliminate oscillations created by discontinuous payoffs (for example, see [24,25]). In Table C.1, as an illustrating example, we present computed option prices obtained by the IE and the FD-CN methods for different initial values of S(0) and J for the case S¯ = 95. In this test, for the benchmark solutions, we use the FD-CN method with a very fine computational grid having t = J = 1/3600, and  log(S ) = 1/1800. We also report the relative errors between the prices obtained by the FD-CN and those obtained by the IE method, defined as

relative error =

|VIE (S(0 ), 0, J ) − VFD-CN (S(0 ), 0, J )| , VFD-CN (S(0 ), 0, J )

where V· (S(0), 0, J) denotes the option price obtained by the corresponding method. From Table C.1, it is clear that the option prices obtained by the IE method agree well with the benchmark values, with relative errors less than 1.1%. This excellent agreement to the benchmark solutions shows the validity of the proposed IE method. 5.2. Efficiency comparison with the FD-CN method In this subsection, we compare the IE method with the FD-CN method in terms of efficiency, i.e. accuracy per unit of cost, which is a more accurate measure than accuracy or cost alone. Since we compare the efficiency of two methods applied to option pricing, it is important to determine the computational cost of each method. The computational cost of the IE method can be approximated as detailed in Section 4.2 via formula (4.2). With respect to the FD-CN method, discretizations of the S and J variables via central finite differences are used, and the Crank–Nicolson timestepping method is employed. At each timestep, the FD-CN method involves solving a total of “# J points” tridiagonal linear systems of size “# S points” × “# S points”, where “# J points” and “# S points” respectively denote the number of discretization points in the J- and S-directions. The cost for solving such a tridiagonal linear system is approximately “# S points” (flops). The

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

341

American constraint is imposed explicitly at each J discrete value. Hence, the total computation cost of the FD-CN method can conservatively be modelled by

total cost = # timesteps × # J steps × # S points (flops).

(5.1)

In Table C.2, we present selected results for an efficiency comparison study between the IE and the FD-CN methods. For the IE method, we separately report the numbers of time steps used by the method over [0, T∗ ] and [T∗ , T], where T ∗ = T − J¯, ˆ in the computational cost formula (4.2). The total numbers of Newton iterations required which respectively are m and m for convergence, denoted by “# Newton iter.”, for [0, T∗ ] and [T∗ , T] are also reported separately. The total computational costs of the IE method for different grid sizes are then approximated using formula (4.2). For illustration purposes, we also report the computational time (“run time”) in seconds. The programs were implemented using Matlab 2015a, and run on a Surface Pro 3 with Intel Dual Core i7-4650U @ 1.70 GHz 2.30 GHz processor and 4 GB RAM. We make the following observations. • The prices (V) obtained by the IE and FD-CN methods appear to converge to the same values, with the point-wise relative errors at the finest grid levels being less than 0.05% for both cases, namely (S = 100, J = 0 ) and (S = 94.9, J = 1/15 − 0.001 ). These convergence results could also serve as a validation test for accuracy of the IE method. • More importantly, the IE method appears to be significantly more efficient than the FD-CN method. For the case of S = 100 and J = 0, as an illustration, to achieve one decimal digit of accuracy, i.e. to get V = 10.1, the FD-CN method needs between 1.9 × 108 flops (which gives V = 10.5902) and 1.5 × 109 flops (which gives V = 10.1741), while the IE methods needs less than 5.1 × 106 flops (which gives V = 10.1307). That is, to achieve one decimal digit of accuracy, the cost required by the FD-CN method appears to be between 38 times and 300 times more than that required by the IE method. In other words, for this case, the IE method is between 38 times and 300 times more efficient than the FD-CN method. For the case of S = 94.9 and J = 1/15 − 0.001, i.e. S ↑ S¯ and J ↑ J¯, which is a challenging case due to high possibility of knock-out, it is observed that the IE method is even more efficient than the FD-CN method. To achieve just one decimal digit of accuracy, the IE method requires less than 9.1 × 106 flops (which gives V = 6.4773), whereas the FD-CN method requires between 1.5 × 109 flops (which gives V = 6.8151) and 1.2 × 1010 flops (which gives V = 6.4798). This makes the IE method between about 168 times and 1318 times more efficient than the FD-CN method in this case. Finally we note that in the efficiency comparison, we use computational complexity (flop counts) instead of computational times. We argue that using computational complexity results in a more accurate efficiency comparison because computational times are implementation-dependent. Nonetheless, it is obvious from Table C.2 that under either measure, the IE method is significantly more efficient than the FD-CN method. As well-known in the computational finance literature, the lattice/grid-based methods cannot handle the knock-out feature well, especially for asset prices near the barrier. As such, the observed efficiency issues that the FD-CN method has in the case S ↑ S¯ and J ↑ J¯ is expected. These results indicate that the proposed IE method are much more robust and efficient than existing finite difference methods in dealing with American Parisian options.

5.3. Optimal exercise boundary In Fig. C.2(a) and (c), we plot, as a function of time to expiry τ , the optimal exercise boundaries of (i) an American-style Parisian down-and-out call option, (ii) its embedded vanilla American call option, and (iii) its American down-and-out call counterpart. Fig. C.2(b) plots zooms to provide a clearer comparison at τ ≈ J¯. It is clear from these figures that, at τ ≤ J¯, the optimal exercise boundary of the Parisian option coincides with that of its embedded vanilla, as expected from Corollary 4.1. Furthermore, when τ > J¯ ≈ 0.067, the optimal exercise boundary of the Parisian option always lies between those of its embedded vanilla option (upper) and its down-and-out counterpart (lower). The above observation is indeed a result of the fact that the knock-out risk inherent in the contract of an American-style Parisian option is higher than that of its embedded vanilla American option (zero knock-out risk), but lower than that of its “one-touch” barrier American option counterpart. More specifically, compared to its vanilla American option counterpart, a Parisian knock-out option has an extra risk that it may be knocked out, i.e. become worthless, before expiry. However, if compared with its barrier option counterpart, the risk of being “knocked out” of the Parisian option is still lower, as its knock-out feature is harder to be activated. By comparing Fig. C.2(a) and (c), we also observe that the optimal exercise boundary of the American-style Parisian option increases when the barrier S¯ decreases. More specifically, when S¯ decreases from $95 to $75, the optimal exercise price of the Parisian option increases towards that of its embedded vanilla American option, which, of course, does not depend on S¯. This can be explained as follows. A smaller value of S¯ would not only reduce the likelihood of the knockout event (down-and-out feature), but also would reduce the loss for the option holder if this happens. In other words, the knock-out feature has less effect on the Parisian option as S¯ decreases, and thus, the difference between the optimal exercise boundaries of the American-style Parisian option and the embedded vanilla American option becomes smaller in this case.

342

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

5.4. Option price In Fig. C.3(a), we plot, as a function of time to maturity τ , the option prices at S = S¯ = 95 of (i) an American-style Parisian down-and-out call option, (ii) its embedded vanilla American call at the same S. Fig. C.3(b) plots zooms to provide a clearer comparison at τ ≈ J¯. Other parameters are the same with the previous test. It is clear from Fig. C.3(a) that, similar to its embedded vanilla American option, the price of the Parisian option increases when time to maturity increases. This is because with a larger τ , the option holder has more chances to optimally exercise the option. Fig. C.3(a),(b) also demonstrate that the price of the American-style Parisian option is always less than or at most equal to that of its embedded vanilla option. This is due to the extra risk of being “knocked out”. In particular, when τ < J¯ ≈ 0.067, this knock-out risk disappears since there is not enough time left for J to reach J¯. As a result, the two prices are equal, as evident in Fig. C.3(b). In Fig. C.4(a) and (b), we plot, as functions of S, the prices of (i) an American-style Parisian down-and-out call option at different fixed values of J, (ii) its embedded vanilla American option, and (iii) its American down-and-out call counterpart, when S¯ = $95 and S¯ = $75, respectively. Similar to observations made for optimal exercise boundaries in Fig. C.2, the price of the American-style Parisian option always lies between those of its embedded vanilla American option and its barrier option counterpart. In addition, it is clear from both Fig. C.4(a) and (b) that the price of the American-style Parisian option is an increasing function of the underlying asset price. This is indeed expected, as an increase in S always brings more benefits to the option holder: it reduces the risk of losing the option and increases the chance of optimally exercising the option. By contrast, when S decreases below the barrier, the American-style Parisian option becomes out-of-money. Even worse, J then starts to accumulate at the same rate of the passing time and the risk of losing the option increases, which devaluates the option. This is clearly demonstrated in both Fig C.4(a) and (b): the price of the Parisian option associated with a higher J is lower than that associated with a lower J. It can be further observed from these figures that only the option price curve associated with J = 0 passes the barrier smoothly, but not those associated with non-zero J values. This is because the Delta condition (2.15) ensures the continuity of the Delta of the option only at J = 0, but not at other non-zero J values. 6. Conclusion and future work In this paper, we present an integral equation approach for pricing American-style Parisian down-and-out call options. The key step of the method is to reduce the pricing problem to a pair of coupled integral equations. Once these integral equations are effectively solved numerically using the Newton–Raphson iterative procedure, the option prices, the optimal exercise prices, and the hedging parameters can be obtained easily. Through selected numerical examples, we validate the proposed numerical methods, and discuss interesting properties of the price and the optimal exercise boundary of Americanstyle Parisian down-and-out call options. Moreover, an efficiency analysis show that the proposed method is significantly more efficient than existing finite difference methods. Extensions of this work may include incorporating time-dependent rebates for the knock-out feature. In this case, the method of [1] can be adapted and combined with the proposed method in this paper. In addition, convertible bonds and callability features in convertible bonds [26], or recoveries in structural models of defaultable bonds [27], are expected to be handled by extensions of the techniques developed in this paper. We plan to investigate these applications in the near future. Another interesting research direction is to extend our approach to a more realistic option pricing model, such as a jumpdiffusion model with local and/or stochastic volatility. Recently published works on an integral equation approach for pricing American options under a jump diffusion model, such as [28,29], can be used as a good starting point for this extension. We leave this as a future work. Appendix A. Proof of (3.8) In this appendix, we sketch a proof for the result (3.8). The main idea is to apply the Fourier sine transform to solve the system B4 . To this end, the x-domain of B4 which is a finite domain, first needs to be extended to a semi-infinite domain first by expressing the PDE as:

H (x f (τ ) − x )

∂ C2 ∂ 2C (x, τ ) = H (x f (τ ) − x ) 22 (x, τ ) ∂τ ∂x

(A.1)

where H(x) is the Heaviside function, defined as in (3.7). As a result of applying the Fourier sine transform with respect to x, the PDE system B4 can be reduced to the following linear first-order ODE in the Fourier sine space:

∂ C&2 (ω, τ ) + ω2C&2 (ω, τ ) = G(ω, τ ) ∂τ where

G(ω, τ ) = sin



     ωx f (τ ) g3 (x f (τ ), τ ) − ω cos ωx f (τ ) g2 (x f (τ ), τ ) + ωg1 (τ ) + x f (τ )g2 (x f (τ ), τ ) sin ωx f (τ ) ,

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

with initial condition C&2 (ω, 0 ) =

% x f (τ ) 0

2 C&2 (ω, τ ) = C&2 (ω, 0 )e−ω τ +

τ

343

g(x ) sin(ωx )dx. The solution of this initial-value ODE can be easily solved as: 2 e −ω ( τ −ξ ) G ( ω , ξ ) d ξ .

0

(A.2)





As C&2 (ω, τ ) denotes the Fourier sine transform of H x f (τ ) − x C2 (x, τ ), from (A.2), we obtain:





2 H x f (τ ) − x C2 (x, τ ) = Fs−1 {C&2 (ω, τ )} = Fs−1 {C&2 (ω, 0 )e−ω τ } + Fs−1



τ 0



2 e −ω ( τ −ξ ) G ( ω , ξ ) d ξ .

(A.3)

Applying the Convolution theorem for the Fourier Sine transform, we obtain:



Fs−1

=

= The calculation of Fs−1 result is:

Fs−1







l

0

%

+∞ 0



'

x2

e− 4τ 2 {C&2 (ω, 0 )e−ω τ } = Fs−1 Fs H (x f (0 ) − x )g(x ) Fc √

πτ



(x−u )2 ( x+u )2 1 H ( x f ( 0 ) − u ) g( u ) √ e− 4τ − e− 4τ du 2 πτ

x f (0 ) 0





(x−u )2 ( x+u )2 g( u ) e− 4τ − e− 4τ du √ 2 πτ

τ −ω2 (τ −ξ ) e G ( ω , ξ )d ξ

(A.4)

is the same as that presented in the appendix A2 in [1]. The obtained

0

2 e −ω ( τ −ξ ) G ( ω , ξ ) d ξ

= −(S¯ex − E )1x=x f (τ ) (x ) + M(x, τ ) +

τ 0

Q(x, τ , ξ , x f (ξ ))dξ ,

where functions M and Q are defined as in (3.9) and (3.10), respectively. The result (3.8) can be now obtained straightforwardly. Appendix B. Proof of Proposition 3.1 ∂V

Note that V1 (S, u; t) has removable singularities at (S¯, u ). In this case, in order to calculate ∂ S1 (S, u; t ) when S closes to S¯, we first need to remove these singularities by using the following variable transformation:

ξ=

ln S − ln S¯ . √ σ v−u

As a result,

V1 (S, u; t ) =

+∞

 S α √2 S¯

S¯−ln S ln√ σ J¯−u

√ e

π

2

− ξ2 + β2



ln S¯−ln S

2

ξ



 (ln S¯ − ln S )2 W t +u+ dξ . σ 2ξ 2

By using the Leibniz integral rule, we can obtain:

√   S α √2 +∞ ln S¯ − ln S − ξ 2 + β  ln S¯−ln S 2 S¯ )2 ∂ 2W (t + J¯) S α β σ22 (J¯−u )− (lnσ 2S−ln ξ (J¯−u ) − V1 (S, u; t ) = e e 2 2 √  2 ln√ S¯−ln S ¯ ¯ ∂S S ξ π S S ¯ Sσ π ( J − u ) σ J¯−u      2 2 ¯ ¯ (ln S − ln S ) 2 t + u + (ln S − ln S ) · βW t + u + + W dξ σ 2ξ 2 σ2 σ 2ξ 2   2 

+∞ √ 2 S 2α Sα −1 − ξ2 + β2 ln S¯−ln (ln S¯ − ln S )2 ξ + e W t + u + dξ . √ α ln√ S¯−ln S σ 2ξ 2 π S¯ ¯ σ

(B.1)

J−u

By using variable transformations η = ξ − formula (3.17).

S¯−ln S ln √ σ J¯−u

and v =

ln S¯−ln S

σξ

for the first and second integral in (B.1), we obtain the

Appendix C. Proof of Corollary 4.1 In order to prove Corollary 4.1, we need to take the limit as t tends to T∗ of both sides of the integral Eq. (3.16). It is clear that

lim∗

t→T

t

T∗

Q (S f (t ), t , v, S f (v ))dv = 0.

344

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

Fig. C1. Pricing domain of the American-style Parisian down-and-out call options. Table C1 Option prices obtained by the IE and the FD-CN methods. S(0)

J

FD-CN

IE

Relative error

110 105 100 95 90 85

0 0 0 0 1/30 1/30

16.9778 13.4799 10.1741 6.9182 3.0171 0.8008

17.0058 13.5157 10.1788 6.9406 3.0491 0.8101

0.16% 0.27% 0.05% 0.32% 1.00% 1.10%

In addition, we have:

 

lim N d1 S f (t ), T ∗ − t , S f (T ∗ )

t→T ∗

 

lim N d1

t→T ∗



 

= lim∗ N d2 S f (t ), T ∗ − t , S f (T ∗ ) t→T

S¯2 , T ∗ − t, S f (T ∗ ) S f (t )



 

= lim∗ N d2 t→T



=

1 , 2

S¯2 , T ∗ − t, S f (T ∗ ) S f (t )

 = 0.

Therefore, it can be easily proved that





lim∗ M S f (t ), T ∗ − t , S f (T ∗ ) =

t→T

 1 S f (T ∗ ) − E . 2

The term R(Sf (t), t) can be expressed as: R(S f (t ), t ) = I1 (S f (t ), t ) − I2 (S f (t ), t ), where

I1 (S f (t ), t ) = I2 (S f (t ), t ) =

S f (T ∗ ) S¯

S f (T ∗ ) S¯

 

S f (t )

v S f (t )

v

α α

eβ

σ 2 (T ∗ −t ) 2



CA (v, T ∗ )

vσ 2π ( T ∗ − t ) eβ

σ 2 (T ∗ −t ) 2



CA (v, T ∗ )

vσ 2π ( T ∗ − t )

e

e

−

−

(ln S f (t )−ln v )2 2σ 2 (T ∗ −t )

d v,

(ln S f (t )+ln v−2 ln S¯ )2 2σ 2 (T ∗ −t )

d v.

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

170

345

123.5

Vanilla

Vanilla 160 150

Parisian

Parisian

Sf

Sf

140 123

Barrier

130

Barrier 120 110 100

0

0.2

0.4

0.6

τ

0.8

1

122.5 0.066

(a) S¯ = $95

0.067

0.068

τ

(b) S¯ = $95 - zoom in

170

Vanilla 160

Parisian

150

Barrier

Sf

140 130 120 110 100

0

0.2

0.4

0.6

τ

0.8

1

(c) S¯ = $75 Fig. C2. Optimal exercise boundaries of an American-style Parisian down-and-out call, its embedded vanilla American call option, and its American downand-out call counterpart.

1.95

12

Vanilla 1.94

8

option price (V)

option price (V)

10

Vanilla

6

Parisian 4

1.92

1.91

2 0

Parisian

1.93

0

0.2

0.4

τ

0.6

(a) option price

0.8

1

1.9 0.066

0.067

τ

0.068

(b) option price-zoom in

Fig. C3. Option prices at S = S¯ = $95, as function of time to maturity τ , of an American-style Parisian down-and-out call and its embedded vanilla American option.

346

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347 Table C2 Efficiency comparison study between the IE and FD-CN methods. Prices at time-0 with the barrier S¯ = 95. The computational cost for the IE and FD-CN methods are approximated by formulas (4.2) and (5.1), respectively. FD-CN S(0)

J

# Time steps

# S points

# J points

Price (V)

Total cost (flops)

Time (seconds)

100

0

450 900 1800 3600

225 450 900 1800

30 60 120 240

7.5228 10.9559 10.5902 10.1741

3.0 × 106 2.4 × 107 1.9 × 108 1.5 × 109

2.5 14.8 91.3 650.8

94.9 (↑ S¯)

1/15 -0.001 (↑ J¯)

900 1800 3600 7200

450 900 1800 3600

60 120 240 480

7.7542 7.3482 6.8151 6.4798

2.4 × 107 1.9 × 108 1.5 × 109 1.2 × 1010

14.9 91.3 650.8 4668.1

J

# Time steps

# Gauss

# Newton iter.

Price (V)

Total cost

Time

(flops)

(seconds)

IE S

∗

100

0

94.9 (↑ S¯)

1/15 - 0.001 (↑ J¯)

∗

∗

∗

[0, T ] (m)

[T , T] (mˆ )

points (g)

[0, T ]

[T , T]

70 140 280 70

5 10 20 5

57 57 57 57

235 415 730 235

17 31 58 17

10.1307 10.1765 10.1788 6.6425

140 280

10 20

57 57

415 730

31 58

6.4773 6.4788

16

× × × ×

106 106 107 106

8.3 24.8 71.3 8.3

9.1 × 106 1.6 × 107

24.8 71.3

5.1 9.1 1.6 5.1

14

Vanilla Vanilla

12

12

Barrier

option price (V)

option price (V)

14

10 8

J=0

6

J = 1/30 4

8 6

J=0 4

J = 1/20 2

2 0 40

60

80

100

Barrier

10

0 40

asset price

J = 1/30

60

J = 1/20 80

100

asset price

(a) S¯ = $95

(b) S¯ = $75

Fig. C4. Comparison between prices of an American-style Parisian down-and-out call, at different J, with that of its embedded option.

Using a variable transformation u =

I1 (S f (t ), t ) =

ln S f (t )−ln S¯

σ

√

2(T ∗ −t )

ln S f (t )−ln S f (T ∗ )

σ

√

ln S f (t ) − ln v

, I1 (Sf (t), t) can be transformed to  σ 2 (T ∗ − t )

√ √ e−u 2 ∗ ∗ ∗ √ CA (S f (t )e−σ u 2(T −t ) , T ∗ )eα uσ 2(T −t )+βσ /2(T −t ) du.

2(T ∗ −t )

2

π

Therefore,

lim∗ I1 (S f (t ), t ) = CA (S f (T ∗ ), T ∗ )

t→T

0

+∞

e−u 1 √ du = CA (S f (T ∗ ), T ∗ ). 2 π 2

Similarly, we can prove that lim I2 (S f (t ), t ) = 0. Now, by taking the limit as t tends to T∗ of both sides of the integral t→T ∗

Eq. (3.16), we obtain:

S f (T ∗ ) − E = CA (S f (T ∗ ), T ∗ ). This implies that S f (T ∗ ) = SVf (T ∗ ), where SVf (t ) denotes the optimal exercise price of the embedded vanilla option at time t. This completes our proof.

N.-T. Le, D.-M. Dang / Applied Mathematics and Computation 305 (2017) 330–347

347

References [1] N.-T. Le, S.-P. Zhu, X. Lu, An integral equation approach for the valuation of american-style down-and-out calls with rebates, Comput. Math. Appl. 71 (2) (2016) 544–564. http://dx.doi.org/10.1016/j.camwa.2015.12.013. [2] M. Chesney, M. Jeanblanc-Picque, M. Yor, Brownian excursions and Parisian barrier options, Adv. Appl. Probab. 29 (1) (1997) pp.165–184. [3] M. Dai, Y.K. Kwok, Knock-in American options, J. Futur. Mark. 24 (2) (2004) 179–192. [4] L. Gauthier, Excursions height-and length-related stopping times, and application to finance, Adv. Appl. Probab. 34 (4) (2002) 846–868. [5] C. Bernard, O.L. Courtois, F. Quittard-Pinon, A new procedure for pricing Parisian options, J. Deriv. 12 (4) (2005) 45–53. [6] F. Moraux, Valuing corporate liabilities when the default threshold is not an absorbing barrier, SSRN Working Paper Series (2002) https://ssrn.com/ abstract=314404. [7] A. Chen, M. Suchanecki, Default risk, bankruptcy procedures and the market value of life insurance liabilities, Insur.: Math. Econ. 40 (2) (2007) 231–255. [8] A. Dassios, S. Wu, Perturbed Brownian motion and its application to Parisian option pricing, Finance Stoch. 14 (3) (2010) 473–494. [9] R.J. Haber, P.J. Schönbucher, P. Wilmott, Pricing Parisian options, J. Deriv. 6 (3) (1999) 71–79. [10] S.-P. Zhu, W.-T. Chen, Pricing Parisian and Parasian options analytically, J. Econ. Dyn. Control 37 (4) (2013) 875–896. [11] M. Chesney, L. Gauthier, American Parisian options, Finance Stoch. 10 (4) (2006) 475–506. [12] Y.K. Kwok, D. Barthez, An algorithm for the numerical inversion of Laplace transforms, Inverse Problems (5) (1989) 1089–1095. [13] A.H.-D. Cheng, P. Sidauruk, Y. Abousleiman, Approximate inversion of the Laplace transform, Mathematica J. 4 (2) (1994) 76–81. [14] G.F. Carrier, C.E. Pearson, Partial Differential Equations: Theory and Technique, Academic Press, (London) LTD, 1976. [15] N.-T. Le, D.M. Dang, T.-V. Khanh, A decomposition approach via Fourier sine transform for valuing American knock-out options with time-dependent rebates, J. Comput. Appl. Math. 317 (2017) 652–671. [16] J. Kevorkian, Partial differential equations, Texts in Applied Mathematics, 35, second ed., Springer-Verlag, New York, 20 0 0. [17] D.J. Duffy, Finite difference methods in financial engineering, Wiley Finance Series, John Wiley & Sons, Ltd., Chichester, 2006. [18] C. Constanda, Solution techniques for elementary partial differential equations, CRC Press, Boca Raton, FL, 2010. [19] H. Hattori, Partial differential equations, World Scientific Publishing, 2013, ISBN 9789814407564. [20] C. Chiarella, A. Kucera, A. Ziogas, A Survey of the Integral Representation of American Option Prices, Technical Report, Quantitative Finance Research Centre, University of Technology, Sydney, 2004. [21] I. Kim, The analytic valuation of American options, Rev. Financial Stud. 3 (4) (1990) 547–572. [22] P.K. Kythe, M.R. Schaferkotter, Handbook of Computational Methods for Integration, Chapman and Hall/CRC, 2014, ISBN 1584884282. [23] R. Rannacher, Finite element solution of diffusion problems with irregular data, Numerische Mathematik 43 (1984) 309–327. [24] M.B. Giles, R. Carter, Convergence analysis of Crank–Nicolson and Rannacher time-marching, J. Comput. Finance 9 (2005) 563–576. [25] P. Forsyth, K. Vetzal, Quadratic convergence of a penalty method for valuing American options, SIAM J. Sci. Comput. 23 (2002) 2096–2123. [26] Y.K. Kwok, K.W. Lau, Pricing algorithms for options with exotic path-dependence, J. Deriv. 9 (1) (2001) 28–38. [27] C. Bernard, O. Le Courtois, F. Quittard-Pinon, A new procedure for pricing Parisian options, J. Deriv. 12 (4) (2005) 45–53. [28] G.H. Cheang, C. Chiarella, A. Ziogas, The representation of american options prices under stochastic volatility and jump-diffusion dynamics, Quantitative Finance 13 (2) (2013) 241–253. [29] C. Chiarella, A. Ziogas, American call options under jump-diffusion processes–a Fourier transform approach, Appl. Math. Finance 16 (1) (2009) 37–79.