Submitted to Operations Research manuscript (Please, provide the manuscript number!)

Simulation of tempered stable L´evy bridges and its applications Kyoung-Kuk Kim Department of Industrial and Systems Engineering, Korea Advanced Institute of Science and Technology, Daejeon 305-701, South Korea, [email protected]

Sojung Kim Department of Mathematical Sciences, Korea Advanced Institute of Science and Technology, [email protected]

We consider tempered stable L´evy subordinators and develop a bridge sampling method. An approximate conditional PDF given the terminal values is derived with stable index less than one, using the double saddlepoint approximation. We then propose an acceptance-rejection algorithm based on the existing gamma bridge and the inverse Gaussian bridge as proposal densities. Its performance is comparable to existing sequential sampling methods such as Devroye (2009) and Hofert (2011) when generating a fixed number of observations. As applications, we consider option pricing problems in L´evy models. First, we demonstrate the effectiveness of bridge sampling when combined with adaptive sampling under finite-variance CGMY processes. Second, further efficiency gain is achieved in terms of variance reduction via stratified sampling. Key words : bridge sampling, L´evy process, saddlepoint approximation, tempered stable subordinator History : This paper was first submitted on May 19, 2014 and has been with the authors for 12 months for 4 revision.

1. Introduction When it comes to simulating sample paths of a stochastic process, one fundamental idea applied to a Brownian motion is bridge sampling. This refers to the procedure such that we first generate a skeleton of a Brownian motion and fill in the details as needs arise by utilizing the readily available conditional distributions. It is contrasted with the typical sequential sampling, moving forward in time. And the underlying process, Brownian bridge, has found numerous applications such as statistical inference, large deviations, and option pricing just to name a few. From the simulation point of view, bridge sampling provides us with the control over the coarseness of a simulated 1

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

2

process, the freedom to combine with sequential sampling, and other possibilities to make use of variance reduction techniques and low-discrepancy methods (Glasserman 2003). Naturally, the extension of such bridge sampling to more advanced stochastic processes has received a growing amount of attention. In this paper, we consider a L´evy process conditioned on the terminal values, that is, (Xt )0≤t≤T given X0 and XT . This generalizes the notion of Brownian bridge and thus is called L´evy bridge. Along with various useful applications of Brownian bridge, there have been extensive studies on the probabilistic properties of such conditional distributions of L´evy processes. Concentrating on rather practical aspects, in this work, we are particularly interested in simulating sample paths of L´evy bridge processes, that is, bridge sampling for L´evy processes. Despite its potential usefulness, simulation methods for L´evy bridges have not been developed except for some special processes. A few well-known examples include gamma bridge, inverse Gaussian (IG) bridge, and squared Bessel bridge. Those processes allow explicit or semi-explicit forms of probability density functions (PDF). Bridge sampling methods for gamma and variance gamma (VG) processes are developed in Ribeiro and Webber (2004), Avramidis, L’Ecuyer, and Tremblay (2003) and Avramidis and L’Ecuyer (2006). For IG and normal inverse Gaussian (NIG) processes, see Ribeiro and Webber (2003). The paper applies the MSH method in Michael et al. (1976) and the authors utilize quite a unique transformation of an IG bridge, resulting in large efficiency gains. For a squared Bessel bridge, its Laplace transform is presented in Pitman and Yor (1982) and a series expansion is available as well. On the other hand, simulation algorithms for diffusion bridges have received much attention over the last fifteen years. Early attempts at bridge sampling for diffusions are based on the MetropolisHastings algorithm and a discrete-time approximation of diffusion processes in Roberts and Stramer (2001) and Durham and Gallant (2002); see also Lin et al. (2010) for an extension. Beskos et al. (2006) introduce a remarkable exact sampling algorithm of general diffusion processes. Their proposed method accepts a sampled skeleton of a standard Brownian motion with some probability

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

3

proportional to the Radon-Nikodym derivative between the laws of the Brownian motion and the diffusion bridge. Lastly, there is a simple approximation scheme for diffusion bridges of Bladt and Sørensen (2014), using two diffusions, one moving forward in time and the other backward, and constructing an approximate bridge. Except for those mentioned above, L´evy bridge sampling methods are hardly known. One of the difficulties comes from the fact that the distribution of a L´evy process is mostly specified via its characteristic function (CHF) or the associated L´evy measure. Therefore, the PDF is only available via numerical inversion or infinite series expansion, which renders some of the prior algorithms impractical. In this respect, the (approximate) PDF or related probabilistic features are reported in the literature. Some representative works include R¨ uschendorf and Woerner (2002) and FigueroaL´opez and Houdr´e (2009). In Figueroa-L´opez and Tankov (2014), the authors obtain small-time asymptotic behaviors of the exit probabilities of a L´evy process out of certain intervals, using L´evy bridges. Unfortunately, such a small-time density expansion is not easily applicable from the simulation perspective because an approximating quantity fails to be a legitimate PDF in general. We, however, note that they describe an interesting bridge sampling algorithm for mid-points if the PDF of a L´evy process is known and unimodal, with an application to stopping times. For more general instances, there are some approximation schemes such as the forward-backward method for Markov bridges in Asmussen and Hobolth (2012) and the beta approximation method for L´evy processes in Glasserman and Kim (2008). We consider stable and tempered stable processes, which have infinite jump activity and are widely used for heavy-tailed models. We propose a simulation algorithm for the L´evy bridges of stable and tempered stable subordinators along with efficient simulation schemes for the stochastic processes based on or derived from them. In fact, their bridge processes have the same law thanks to exponential tilting when the stable index is less than one. For this purpose, we obtain an approximate conditional PDF given the end-points via saddlepoint techniques. The resulting closed form conditional PDF enables us to consider an acceptance-rejection algorithm in which two known proposal densities are utilized, namely the gamma bridge and the IG bridge.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

4

The saddlepoint method obtains asymptotic expansions by approximating a contour integral in the complex plane near a saddlepoint, via the steepest descent method. It is first introduced in statistics by Daniels (1954) to provide an approximate PDF of the mean of i.i.d. random variables, and later the tail probabilities of the sample mean are derived in Lugannani and Rice (1980). In financial applications, many authors have applied saddlepoint methods to option pricing under various models such as affine jump-diffusions or credit models. The reader is referred to Rogers and Zane (1999), Carr and Madan (2009), Glasserman and Kim (2009),Yang et al. (2006), Dembo et al. (2004), and Gordy (2002) etc. The performance and advantages of the proposed bridge sampling scheme are demonstrated via several numerical studies. We particularly consider option pricing problems under L´evy models. The main messages from our numerical tests can be summarized as follows: • It is comparable in cost and accuracy to three existing sequential simulation methods if we

generate a fixed number of observations. The first method is an acceptance-rejection method based on the representation of a stable random variable in Chambers et al. (1976), the second is the method of Devroye (2009) with uniformly bounded complexity. The last under consideration is the fast rejection method proposed by Hofert (2011). • It allows adaptive sampling for pricing path-dependent options under finite variation tem-

pered stable processes (e.g., CGMY processes). Such an adaptive technique extended from Becker (2010) results in remarkable savings in simulation costs. • It makes it possible to use stratified sampling techniques for pricing path-dependent options

beyond VG and IG processes, which can lead to significant variance reduction. • It is combined with sequential sampling utilizing the strengths of both, and such a hybrid

method is applied for least square Monte Carlo (LSMC) methods for American options valuation under asset dynamics with subordinated Brownian motions. This avoids the massive data storage requirements for LSMC pricing without sacrificing precision. (See the e-companion to this paper.)

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

5

In addition to the pricing of path-dependent options, it is worth noting that there are prospective applications of our approach to statistical inference and the information-based asset pricing framework. Inference based on incomplete data, the so-called missing data problem, is a fundamental issue in statistics. By treating the data observed at a discrete set of points as a missing data problem and generating intermediate values conditioned on endpoints, one can exploit inference tools designed for continuous paths. On the other hand, the essential idea of the information-based asset pricing is to model the flow of information in markets as stochastic processes. For a market factor ST revealed at time T , the corresponding information process (ξtT )0≤t≤T conditioned on ξT T = ST that market participants have about ST is generated. Brody et al. (2007, 2008) adopt Brownian bridge and gamma bridge as information processes; Hoyle et al. (2011) use more general L´evy bridges for such information processes. The rest of this paper is organized as follows. Section 2 provides more detailed background on stable and tempered stable processes as well as on saddlepoint methods. In Section 3, we explain some useful properties of a tempered stable L´evy bridge and derive its approximate PDF. Then we develop a bridge sampling scheme later in the section. In Section 4, numerical tests are conducted and reported under finite variation tempered stable processes. Stratified sampling for variance reduction is considered as well. Section 5 concludes.

2. Preliminaries 2.1. Stable and tempered stable processes A L´evy process (Xt )t≥0 is a stochastically continuous process with stationary and independent increments starting at 0. For every t, Xt has an infinitely divisible distribution. Conversely, if µ is an infinitely divisible distribution then there exists a L´evy process (Xt )t≥0 such that the distribution of X1 equals µ. We first define a stable distribution on the positive real line and the L´evy process corresponding to the stable distribution. The general definition of a stable process is given in EC.1.3 of the e-companion to this paper; here we focus on the case 0 < α < 1.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

6

Definition 1. Let µ be an infinitely divisible probability measure on (R+ , B (R+ )). It is called α-stable, denoted by S(α, c), with an index of stability 0 < α < 1 and c > 0 if it has the CHF Z exp


eiux − 1 ν(dx) , u ∈ R,

(1)

R

with the L´evy measure of the form ν(dx) =

c x1+α

1(0,∞) (x)dx.

(2)

The corresponding nonnegative L´evy process Xt is called a stable (L´evy) subordinator, whose distribution is S(α, ct) for each t > 0. The PDF of S(α, c) is only known as an infinite series of the following form: fS(α,c) (x) =

X (−1)n−1

1 π (−cΓ(−α))

1/α n∈N

n!

sin(nπα)Γ(nα + 1)

!−nα−1

x (−cΓ(−α))

1/α

.

(3)

For some α, such as α = 1/3 or 1/2, the closed form expression is known in terms of some special functions. Nevertheless, an exact simulation method for a stable distribution S(α, c) exists through the representation (Chambers et al. 1976, Kawai and Masuda 2011): 1/α sin(αU + ϑ)  cos((1 − α)U − ϑ)  1−α α S(α, c) = − cΓ(−α) , 1/α (cos U ) E 

(4)

where ϑ = πα/2, U has a Unif(−π/2, π/2) distribution, and E is exponentially distributed with mean 1 independently of U . Now, let us consider an exponentially tempered version of a stable subodinator, that is, ν(dx) = R eθx νS (dx) for some θ that satisfies |x|≥1 eθx νS (dx) < ∞ where νS (dx) is the L´evy measure (2). Definition 2. An infinitely divisible distribution η on (R+ , B (R+ )) is called one-sided tempered stable distribution, denoted by TS(α, λ, c), with parameters c, λ ∈ (0, ∞) and α ∈ [0, 1) if its CHF is given by (1) with the L´evy measure ν(dx) =

c −λx e 1(0,∞) (x)dx. x1+α

The L´evy process (Xt )t≥0 associated with η is called a tempered stable (L´evy) subordinator.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

7

Clearly, the distribution of Xt is given by TS(α, λ, ct) for each t > 0. The parameter c represents the intensity of jumps, and α determines the relative importance of small jumps for the trajectories of the process. Here, α = 0 case is included for more generality as in Cont and Tankov (2004). The new exponential tempering parameter λ controls the decay rate of large jumps. The properties of a one-sided tempered distribution relative to a stable distribution have been investigated (for example, see Rosi´ nski 2007). A one-sided tempered stable distribution also admits a series expansion of the PDF on the positive real line, by the following relation to the PDF fS(α,c) (x): h i fTS(α,λ,c) (x) = exp − λx − cΓ(−α)λα fS(α,c) (x).

(5)

Based on this relation, TS(α, λ, c) can be simulated exactly by the acceptance-rejection sampling, which is described in Algorithm 1. Algorithm 1 Exact sampling algorithm for a tempered stable subordinator [SSR] 1: repeat 2:

Generate X ∼ S(α, ct) from (4) and an independent U ∼ Unif(0, 1)

3:

until U ≤ e−λX .

4:

return X

Algorithm 1, which we call simple stable rejection (SSR), has the complexity of O (exp (tλ/α)). Its performance deteriorates quickly for large t, λ, or small α. Devroye (2009) develops an exact double rejection method that is uniformly bounded in complexity over all parameter ranges to address such possible low acceptance rates. Recently, Hofert (2011) suggests a fast rejection algorithm that obtains a logarithmic improvement in complexity O (tλ/α) in comparison to SSR. A tempered stable subordinator can be naturally extended to the whole real line. Definition 3. For fixed parameters c+ , c− , λ+ , λ− ∈ (0, ∞) and α+ , α− ∈ [0, 1), an infinitely divisible distribution η on (R, B (R)) is called a tempered stable distribution, denoted by TS (α+ , λ+ , c+ ; α− , λ− , c− ) , if η = η + ∗ η − , where η + = TS (α+ , λ+ , c+ ) and η − = υ e with υ =

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

8

TS (α− , λ− , c− ) and denoting the dual of υ by υ e(B) = υ(−B) for B ∈ B (R). Its L´evy measure ν is given by  ν(dx) =

c+ x1+α+

e

−λ+ x

 c− −λ− |x| 1(0,∞) (x) + 1+α− e 1(−∞,0) (x) dx. |x|

We call the L´evy process associated with η a tempered stable process. Remark 1. In Cont and Tankov (2004), generalized tempered stable processes for parameters c+ , c− , λ+ , λ− ∈ (0, ∞) and α+ , α− < 2 are defined. For α+ , α− ∈ [0, 1), which we consider in this paper, X is a finite-variation process making infinitely many jumps in a finite horizon. For α+ , α− ∈ [1, 2), we have an infinite-variation process. As a special case, when c+ = c− , α+ = α− , it is called a CGMY process in Carr et al. (2002). For more information on tempered stable processes, we refer the reader to K¨ uchler and Tappe (2011), Cont and Tankov (2004), Sato (1999) and the references therein. Remark 2. The term ‘tempered stable process’ sometimes refers to a more general class of processes that contains the version in Remark 1 as a subclass. For example, see Rosi´ nski (2007) and Grabchak (2012). In this broader context, the tempered stable processes that we consider correspond to ‘exponentially tilted stable processes’.

2.2. Saddlepoint approximations for conditional probabilities Daniels (1954) introduces the saddlepoint method to statistics in order to approximate the PDF of the mean of n i.i.d. random variables. In what follows, we present the result in the case n = 1 but for multivariate random vectors. Suppose X = (X1 , · · ·, Xm )> is a random vector with a non-degenerate distribution in Rm . With all continuous components Xi , let f (x) be the PDF at each x ∈ Rm . With χ the support of the random vector X, we define its joint moment generating function (MGF) as 
 M (u) = E exp u> X =

Z χ


exp u> x f (x)dx

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

9

for all values of u = (u1 , · · ·, um ) ∈ Rm such that the integral converges. We assume that the maximal convergence set S for M contains a neighborhood of 0 ∈ Rm . Clearly, S is a convex set. Let Iχ ⊂ Rm denote the interior of the convex hull of χ and K(u) be the cumulant generating function (CGF) of f defined by K(u) = ln M (u). Then, the saddlepoint approximation to unknown f (x) for the continuous vector X is given in Iχ ⊂ Rm as fˆ(x) =

  1 > ˆ exp K(ˆ u ) − u x , (2π)m/2 |K 00 (ˆ u)|1/2

x ∈ Iχ

(6)

ˆ is the unique solution in S to the m-dimensional saddlepoint equation K 0 (ˆ where u u) = x. A natural approach to obtain a saddlepoint approximation for a conditional PDF is to use two saddlepoint approximations, one for the joint PDF and another for the marginal as in (6). The idea is introduced in Daniels (1958) and the explicit formula is presented in Barndorff-Nielsen and Cox (1979). To be specific, let (X, Y) be a random vector having a non-degenerate distribution in Rm with dim(X) = mx , dim(Y) = my , and mx + my = m. With all components continuous, suppose there is a joint PDF f (x, y) with support (x, y) ∈ χ ⊂ Rm . The double saddlepoint approximation for f (y|x), the conditional PDF of Y at y given X = x, is expressed as fˆ(y|x) = (2π)

my /2

× exp





ˆ )| |K 00 (ˆ u, v 00 |Kuu (ˆ u0 , 0)|

−1/2

  ˆ) − u ˆ >x − v ˆ > y − K(ˆ ˆ> K(ˆ u, v u0 , 0) − u 0x

(7)

ˆ ) solves the set of m equations K 0 (ˆ ˆ) = for (x, y) ∈ Iχ . Here, the m-dimensional saddlepoint (ˆ u, v u, v ˆ 0 is the mx -dimensional saddlepoint for the denominator that solves K 0 (ˆ (x, y), and u u0 , 0) = x. 00 ) Here, K 0 (Ku0 ) is the gradient with respect to both components u and v (u only), and K 00 (Kuu

is the corresponding Hessian. For the comprehensive treatments of saddlepoint approximations, good references are Butler (2007) and Jensen (1995).

3. Bridge Sampling Methods In this section, we study the bridge process of a stochastic process (Xt (α, λ, c))t≥0 associated with TS(α, λ, c), and we call it a tempered stable L´evy bridge. Note that gamma bridge and IG bridge

10

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

correspond to α = 0 and α = 1/2, respectively. Let us first describe some of interesting properties of a tempered stable L´evy bridge. Later we will propose an efficient bridge sampling algorithm. All proofs developed in this section are deferred to EC.1 of the e-companion. It is well known that a tempered stable subordinator possesses the following scaling property: For every r > 0, rXt (α, λ, c) has the same law as Xrα t (α, λ/r, c). Accordingly, the resulting bridge process preserves the scaling property as well and we record it in the next proposition. Proposition 1. Let Xt (α, λ, c) be a tempered stable subordinator. For any r > 0, we have d



( Xt | XT = z) =

 1 Xrα t (α, λ/r, c) Xrα T (α, λ/r, c) = rz . r

Another interesting feature for such a process is that it is identical in law to the bridge process for a stable subordinator, again named a stable L´evy bridge. Lemma 1. The distribution of a stable L´evy bridge is equal to that of its tempered version. Thus throughout the rest of the paper, our illustration will be based on a tempered stable process for which any bridge sampling method is applicable to stable L´evy bridges as well. Additionally, we can use known distributional features of a stable L´evy bridge for its tempered counterpart. For instance, we are able to compute conditional moments for the latter process, based on the theory of multi-dimensional stable processes. In the following proposition, we prove that the conditional expectation given initial and end points is simply a linear interpolation of those points. Proposition 2. Let Xt (α, λ, c) be a tempered stable subordinator. For 0 < t1 < t2 , given Xt2 = z the conditional expectation is linear in z, that is, E[Xt1 |Xt2 = z] = (t1 /t2 )z.

3.1. Approximate conditional PDF Thanks to the scaling property presented above, it is sufficient to consider tempered stable L´evy subordinators with the moment condition E[Xt ] = t, from which we have a two-parameter family. For notational convenience, we denote Xt (α, κ) for the two-parameter tempered stable L´evy

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

11

subordinator associated with a one-sided tempered stable distribution, say TS(α, κ), whose L´evy measure is 1 ρ(x)dx = Γ(1 − α)



1−α κ

1−α

e−(1−α)x/κ dx, x1+α

and the MGF is given by E[e

uXt

tl(u)

]=e



,

1−α ∀u ∈ −∞, κ



where the Laplace exponent l(·) is  1−α  α  α  1−α 1−α 1 1−α − −u . l(u) = α κ κ κ

(8)

In this parametrization, λ = (1 − α)/κ and c = λ1−α /Γ(1 − α). The first four central moments are computed as K1 = E[Xt ] = t, K2 = Var[Xt ] = κt, K3 = (2 − α)κ2 t/(1 − α), and K4 = (3 − α)(2 − α)κ3 t/(1 − α)2 . Consider a two-dimensional random vector (Xt2 , Xt1 ) with 0 < t1 < t2 , where Xt = Xt (α, κ) has the CGF tl(u) as in (8). Then the joint CGF K(u, v) of (Xt2 , Xt1 ) is computed as follows:   K(u, v) = ln E euXt2 +vXt1   = ln E e(u+v)Xt1 +uXt2 −t1

(by stationary increment property)

= t1 l(u + v) + (t2 − t1 )l(u) (by independent increment property)    t2 λ − t2 −t1 λ1−α (λ − u)α − t1 λ1−α (λ − u − v)α , α ∈ (0, 1); α α α =
  − t1 log 1 − u+v − t2 −t1 log 1 − u
, α = 0, κ κ κ κ

(9)

with λ = (1 − α)/κ. Note that K(u, v) is well-defined in an open neighborhood of (0, 0). In order to derive an approximate conditional PDF of Xt1 given Xt2 = z, denoted by fXt1 |Xt2 (x|z), we apply the double saddlepoint density approximation (7). Theorem 1. Let Xt be a tempered stable L´evy subordinator associated with TS(α, κ), α ∈ [0, 1), and κ > 0. Then the saddlepoint approximation fˆXt1 |Xt2 (x|z) for the conditional PDF fXt1 |Xt2 (x|z) conditioned on Xt2 = z, where 0 < t1 < t2 , is as follows. If α ∈ (0, 1), then 2−α  1   2(1−α)  1 1 t1 (t2 − t1 ) 2(1−α) z ˆ fXt1 |Xt2 (x|z) = ·√ N (z; α, κ; t1 , t2 ) t2 (z − x)x " (2πκ1 )# 1 1 2 (1 − α) t1 1−α (t2 − t1 ) 1−α t2 1−α × exp − + − α , α α ακ x 1−α (z − x) 1−α z 1−α

(10)

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

12

where N (z; α, κ; t1 , t2 ) is a normalizing constant. If α = 0, then for some normalizing constant N 0 t2 +1

fˆXt1 |Xt2 (x|z) =

t2 −t1 t1 t2 t2κ 2 1 1 · ·√ (z − x) κ −1 x κ −1 z 1− κ . t1 1 0 t −t N (z; α, κ; t1 , t2 ) (t − t ) 2 κ 1 + 12 t κ + 2 2πκ 2 1 1

Recall the two special cases of α = 0 (gamma bridge) and α = 1/2 (IG bridge) whose exact conditional densities fXt1 |Xt2 (x|z) are known from Ribeiro and Webber (2003, 2004): fXGB (x|z) = t |Xt 1

2

1 Beta

t1

t1 t2 −t1 , κ κ


x κ −1 (z − x)

t2 −t1 κ −1

t2

z 1− κ ,

(11)

and fXIGB (x|z) t1 |Xt2

t1 (t2 − t1 ) 1 · =√ t2 2πκ



z (z − x)x

3/2

   1 t1 2 (t2 − t1 )2 t2 2 exp − . + − 2κ x z−x z

(12)

Somewhat surprisingly, the exact densities and the approximate densities from Theorem 1 coincide in these special cases. Remark 3. An unscaled version of the approximate conditional PDF fˆXt1 |Xt2 (x|z) associated 1 − (1−α)

with TS(α, λ, c) can be easily computed by simply replacing κ with (1 − α)(Γ(1 − α)c)

. This

indicates that the conditional law depends only on the parameters α and c, as expected from the scaling property. Remark 4. The same procedure can be conducted for the CGF of a stable process, and it gives the equivalent result as in Theorem 1. This is also expected from Lemma 1. Remark 5. The double saddlepoint approximation above agrees with the conditional density fˆXt1 (x)fˆXt2 −t1 (z − x)/fˆXt2 (z), where we use the one-dimensional single saddlepoint approximation (6) for the PDF fXt (x) with the CGF K(u) = tl(u) as in (8): " # 1 1 2 2(1−α) (1−α) 1 t (1 − α)t (1 − α)x (1 − α) t fˆXt (x) ∝ √ · 2−α exp − − · . α ακ κ ακ 2πκ x 2(1−α) x (1−α) Note that the double saddlepoint approximation is determined from the associated joint CGF while the one-dimensional version is derived directly from the univariate CGF. The two approaches are not always seen to agree. However, it is not difficult to see that the two conditional approximate PDFs are the same for L´evy processes due to their stationary and independent increment properties.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

9

3 α=0 (Gamma) α=0.2 α=0.5 (IG) α=0.8

7

α=0 (Gamma) α=0.2 α=0.5 (IG) α=0.8

2.5 coditional density f(x|z)

8

coditional density f(x|z)

13

6 5 4 3

2

1.5

1

2 0.5 1 0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

(i) Figure 1

0.5 x

0.6

0.7

0.8

0.9

1

(ii)

Approximate conditional PDF fˆXt1 |Xt2 (x|z) given z = 1 with t2 = 1/2 and κ = 0.2505: (i) t1 = 1/3, and (ii) t1 = 1/4 (symmetric case).

Remark 6. An extension of Theorem 1 to a multivariate random vector is given in EC.2.1 of the e-companion. Figure 1 shows the graphs of the approximate conditional PDF fˆXt1 |Xt2 (x|z) for z = 1 in α where parameters are given as t2 = 1/2 and κ = 0.2505. The graphs for α = 0, 1/2 are the exact conditional densities of gamma bridge and IG bridge. As seen in Figure 1(i), the density achieves its maximum on [z/2, z) when t2 − t1 ≤ t1 ; we can also observe from the explicit form (10) that it attains the maximum on (0, z/2] when t2 − t1 > t1 . In the right panel (ii), on the other hand, we see the symmetric shape of conditional densities when t1 = t2 /2. The graphs give us a clue in choosing proposal densities for the acceptance-rejection algorithm that we design in Section 3.2. Even though an error analysis of saddlepoint approximation is difficult, it is still possible to test the accuracy of the approximation in Theorem 1 numerically. The exact conditional PDF can be obtained from Lemma 1 and the PDF (3) of a stable process. With the new parametrization (α, κ), the PDF (3) is re-written as 1 X (−1)n−1 sin(nπα)Γ(nα + 1) fXt (u) = π n∈N n!



1−α κ

1−α

t α

!n u−nα−1 .

(13)

Our MATLAB implementation retains enough terms in the infinite sum in (13) so that there is essentially no difference in numerical outcomes. Figure 2 compares two conditional PDFs where

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

14

the parameters are given as α = 0.3, κ = 0.2505, z = 1, t1 = 1/4, and t2 = 1/2. The average of the relative errors with respect to the true conditional PDF 1 − fˆXt1 |Xt2 (x|z)/fXt1 |Xt2 (x|z) over the range of x, i.e. (0, z), is 0.9%.

1.4 Exact PDF SPA PDF

coditional density f(x|z)

1.2 1 0.8 0.6 0.4 0.2 0

Figure 2

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Comparison of the approximate conditional PDF with the exact PDF for parameters α = 0.3, κ = 0.2505, z = 1, t1 = 1/4, and t2 = 1/2.

We repeat the same exercise for different κ and t1 /t2 values as they turn out to be important factors that affect the performance of the approximation based on extensive numerical tests. Figure 3(i) and (ii) plot the average and maximum of relative errors over κ and t1 /t2 , respectively. While κ varies, the other parameters are fixed as done in Figure 2. When we compute the relative errors, we set the step size of x equal to 0.001. While t1 /t2 varies, κ is set to be 1 and the rest of the parameters remain the same. In Figure 3(i), the average relative errors are less than 4% but the maximum errors tend to increase as κ increases. On the other hand, in Figure 3(ii) it is shown that the averages are more or less small (less than 4%) in the most of the region except for the extreme cases where t1 /t2 is close to 0 or 1. This implies that sampling at middle points would produce the highest accuracy. This phenomenon naturally leads us to the idea of using a bisection method when generating the entire trajectories of a process of interest in order to guarantee good performance, and this is exactly what we propose in Algorithm 4.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

15

Another important parameter under consideration is α, however, the behavior of the series representation is quite unstable for the parameter domain (0, 1), making the overall comparison of two approaches impractical. Nevertheless, we observe from many individual trials that the saddlepoint method yields a high accuracy near α = 1/2.

0.2 Maximum of RE Average of RE

Relative error

0.15

0.1

0.05

0

0.5

1

1.5 κ

2

2.5

(i)

Relative error

0.4

Maximum of RE Average of RE

0.3 0.2 0.1 0

0.1

0.2

0.3

0.4

0.5 t1/t2

0.6

0.7

0.8

0.9

(ii) Figure 3

Average and maximum relative errors with α = 0.3, z = 1, and t2 = 1/2: (i) κ varies with t1 = 1/4 (ii) t1 /t2 varies with κ = 1.

3.2. Sampling algorithms To sample from (Xt1 |Xt2 = z), utilizing fˆXt1 |Xt2 (x|z) in Theorem 1, we propose an acceptancerejection algorithm that uses the known PDFs of the gamma and IG bridges as proposal densities. More precisely, when α < 1/2, we let our proposal density be the conditional PDF of a gamma process (α = 0) of (11). If α > 1/2, then we choose to use the conditional PDF of an IG process (α = 1/2) of (12). For these proposal densities, an acceptance-rejection method becomes implementable thanks to the following result.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

16

Proposition 3. The ratio of fˆXt1 |Xt2 (x|z) relative to the gamma bridge density (11) if α ∈ (0, 1/2) or to the IG bridge density (12) if α ∈ (1/2, 1) is uniformly bounded on (0, z). The ratios of fˆXt1 |Xt2 (x|z) and the proposal densities appear as constant times some auxiliary functions, say g1 (x) and g2 (x) as in EC.1.5. Let us denote the respective upper bounds of those functions by Cg1 and Cg2 . We first begin with the bridge sampling algorithm for IG subordinators, introduced in Ribeiro and Webber (2003) in Algorithm 2. Then in Algorithm 3, our sampling algorithm for tempered stable L´evy bridges is presented. Algorithm 2 Inverse Gaussian Bridge sampling algorithm [IGB] 1: Generate Q ∼ χ2 (1) q 3 2 t −t zκ zκ 1) 1) 4Q (t2t−t Q2 and s2 = 2: Compute s1 = 2t 1 + Q 2t + (t2 −t 2 − 2t (t −t ) t2 1 1 2 1 1 zκ 1

1

t2 −t1 (1 + s1 ) ÷ t1

3:

Generate Ber ∼ Bernoulli(p1 ) with p1 =

4:

Set R ← Ber · s1 + (1 − Ber) · s2 and X ←

5:

return X

t2 t2 −t1 ( t1 t1

(t2 −t1 )2 t2 1 s1

+ s1 )

z 1+R

In applying Algorithm 3 which we call [TSLB] hereafter, it is crucial to compute tight upper bounds Cg1 and Cg2 for efficient simulation. Unfortunately, analytical computations for the maxima of g1 (x) and g2 (x) on the interval (0, z) turn out to be very difficult albeit their explicit expressions. Any numerical methods to compute Cg1 and Cg2 can be applied. One straightforward method is to compute the maximum values at discrete points on the half interval [0, t2 /2] or [t2 /2, t2 ], depending on the relative size of t1 versus t2 , and then we can interpolate them. A bisection method can be an alternative to find a smaller interval that contains the maximum. Numerical experiments show that the efficiency of such computational methods for Cg1 and Cg2 varies according to parameter settings. We may further consider tabulating those upper bounds in terms of z and δ = t1 /t2 via the scaling property of a bridge process in Proposition 1. More specifically, suppose that we have a −1/α

table of Cg1 or Cg2 given parameter set (α, λ, c) at t2 = 1. For general t2 , let r = t2 . Then    1 1 d Xrα δt2 (α, λ/r, c) Xrα t2 (α, λ/r, c) = z Xδt2 (α, λ, c) Xt2 (α, λ, c) = z = r r

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

17

Algorithm 3 Tempered Stable L´evy Bridge sampling algorithm [TSLB] 1: if 0 < α < 12 then repeat

2: 3:

Compute Cg1 given z, t1 , t2

4:

Generate B ∼ Beta

5:

Set X ← zB

t1 t2 −t1 , κ κ




and U ∼ Unif(0, 1)

until U ≤ g1 (X)/Cg1

6: 7:

return X

8:

else if

1 2

< α < 1 then

repeat

9: 10:

Compute Cg2 given z, t1 , t2

11:

Generate Y from [IGB] and U ∼ Unif(0, 1)

12:

Set X ← Y until U ≤ g2 (X)/Cg2

13: 14:

return X

15:

end if  =

 1 Xδ (α, λ/r, c) X1 (α, λ/r, c) = rz . r

(14)

Since the parameter λ does not affect the bridge process, one can sample from (Xδ |X1 = rz) by 1/α

using an interpolated bound from the table and then multiply it by t2 . When t1 = t2 /2, which commonly happens in financial applications, it suffices to compute a one-dimensional table for z only. Remark 7. In order to improve the efficiency of [TSLB], we may approximate the bridge process as the conditional expectation by Proposition 2 when z is very small. Since 0 ≤ (Xt1 |Xt2 = z) ≤ z, the error |(Xt1 |Xt2 = z) − zt1 /t2 | does not exceed z. If z is very small and additionally time points t1 and t2 are small, the samples from the gamma bridge and IG bridge proposals may

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

18

not be distinguishable from zero numerically. Therefore, sampling from [TSLB] is relatively timeconsuming, and using the conditional expectation could be a reasonable alternative for speed-up. The expected number of iterations Mi , (i = 1, 2) of [TSLB] can be computed as (i) 0 < α < 21 :

fˆXt1 |Xt2 (x|z)

t1 t2 − t1 = Beta , GB ft1 |t2 (x|z) κ κ

M1 = sup x∈(0,z)

(ii)

1 2





t2

z κ −1

Cg1 , N

< α < 1: M2 = sup x∈(0,z)

√   Cg2 t2 2πκ − 3 t22 2 = z exp − . IGB ft1 |t2 (x|z) t1 (t2 − t1 ) 2κz N

fˆXt1 |Xt2 (x|z)

Here, the normalizing constant N is given by Z N=

z



x(z − x)

−1− p2

 exp −

0

1 p(1 + p)κ



t1+p (t2 − t1 )1+p 1 + xp (z − x)p

 dx.

Figure 4 plots three-dimensional shaded surfaces of M1 and M2 over a rectangular grid (t1 /t2 , κ) such that 0 < t1 /t2 < 1 with step size 0.01 and 0 < κ < 2 with step size 0.05. Other parameters are

12

Expected number of iterations

Expected number of iterations

fixed as z = 1 and t2 = 1; α = 0.3 for M1 and 0.7 for M2 .

10 8 6 4 2 0 2 1.5 1 κ

0.5 0

0

0.2

0.4

0.6 t1 / t2

(i) Figure 4

0.8

1

5 4 3 2 1 2 1.5 1 κ

0.6

0.5 0

0

0.2

0.8

1

0.4 t / t 1 2

(ii)

Expected number of iterations of [TSLB] algorithm as t1 /t2 and κ vary where z = 1, t2 = 1: (i) M1 with α = 0.3 (ii) M2 with α = 0.7.

First of all, Mi rapidly increases as t1 /t2 goes to the extremes 0 or 1. In these cases, an existing sequential simulation algorithm is recommended, considering that the accuracy of the approximate

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

19

conditional PDF is low as well. A guideline for the specific choice of such an algorithm is briefly described in EC.3. However, except for these extremes, the number of iteration of [TSLB] is mostly less than 2 in expectation. In particular, when t1 /t2 = 1/2, the average M1 and M2 over κ ∈ (0, 2) with step size 0.05 are 1.2658 and 1.3636, respectively. The behaviors of Mi ’s with respect to other parameters κ and α are also observed. In the same figure, we see that the expected number of iterations tends to increase as κ increases, however, the growth rate is quite low compared to their sensitivities with respect to t1 /t2 near extremes. As for the stability index α, Figure 5 shows the average values of Mi over κ ∈ (0, 2) with step size 0.05 versus the parameter α. As expected from Figure 1, Mi ’s increase in α but again they stay below 2 in a large region of α.

1.5

4.5

1.45 4 3.5

1.35

Averaged M2 over κ

Averaged M1 over κ

1.4

1.3 1.25 1.2 1.15

3 2.5 2

1.1 1.5 1.05 1 0.05

0.1

0.15

0.2

0.25 α

0.3

0.35

0.4

0.45

1 0.55

0.6

0.65

(i) Figure 5

0.7

0.75 α

0.8

0.85

0.9

0.95

(ii)

Averaged Mi of [TSLB] algorithm over κ ∈ (0, 2) as α varies, where z = 1, t2 = 1 and t1 /t2 = 1/2 are fixed: (i) M1 with 0.05 < α < 0.45 (ii) M2 with 0.55 < α < 0.95.

In principle, saddlepoint methods for conditional distributions are applicable to any L´evy process as long as its Laplace transform is available. Additionally, if a process of interest belongs to some parametric class, one of which has a closed form PDF, then a non-Gaussian base would likely result in excellent performance for the process as mentioned in EC.2.2. However, in most cases, it is difficult to obtain the saddlepoint explicitly or even to approximate it. Moreover, often

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

20

the complexity of the formula (7) makes the task of sampling from the conditional PDF highly nontrivial.

3.3. Discussion on performance of the proposed method For the rest of this section, we make a short discussion on the relative performance of [TSLB] in comparison to existing simulation algorithms for tempered stable distributions. In particular, we consider the three algorithms mentioned below Algorithm 1. Following the guideline in EC.3, one can find the best performing algorithm for a given set of parameters. This combined exact algorithm is denoted by hybrid sequential algorithm, or HSQ in short. Let Xt be a tempered stable L´evy subordinator associated with TS(α, κ). We generate the same number of sample paths (Xt1 , Xt2 , · · · , Xtd ) at given time points {t1 , . . . , td } with constant step size ∆ = ti+1 − ti . When simulating a trajectory via bridge sampling (BS), we first generate the terminal value Xtd using HSQ or the numerical Laplace inversion of a tempered stable distribution. At the given end point, we fill the intermediate samples by bisecting time points. Algorithm 4 describes the BS procedure with a possible tabulation where a time grid with d = 2m for some integer m is given. We compare computational times to generate the same number of trajectories of (X∆ , · · · , Xd∆ ) for a fixed ∆. We refer the reader to EC.3 for detailed numerical outcomes. To sum up, we observe that the HSQ algorithm is more efficient than BS without tabulation for generating a given number of points in general. However, bridge sampling can be an attractive alternative when we tabulate bounding constants, and the computational burden of this tabulation is small. In our experiments, BS is comparable in speed to HSQ without losing its high accuracy. Furthermore, BS will boost its strength when variance reduction techniques or low-discrepancy methods can be adopted. Moreover, the advantages of BS stand out in settings where we need to fill in intermediate values, as is illustrated in the next section.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

21

Algorithm 4 Dyadic generation of a sample path of a tempered stable L´evy subordinator via [TSLB] without or with tabulation 1: At given time points 0 = t0 < t1 < · · · < td with d = 2m , 2:

X0 ← 0 and generate Xtd ∼ TS(α, (1 − α)/κ, td ((1 − α)/κ)(1−α) /Γ(1 − α))

3:

for l ← 1 to m do

4:

n ← 2m−l

5:

for j ← 1 to 2l−1 do

6:

i ← (2j − 1)n

7:

if the pre-computed table of [zvec ; Cvec ] bounds for (X1/2 |X1 ) is given then

8:

Set r ← (ti+n − ti−n )−1/α

9:

Find C that is interpolated with r(Yti+n − Yti−n ) given the table [zvec , Cvec ]

10:

Set Cgi ← C for i = 1, 2 in line 3 and 10 of [TSLB]

11:

Generate Y from [TSLB] with t1 = 1/2, t2 = 1, and z = r(Yti+n − Yti−n )

12:

Set Y ← Y /r

13:

else

14:

Set t1 = ti − ti−n , t2 = ti+n − ti−n and z = Yti+n − Yti−n

15:

Generate Y from [TSLB] with t1 , t2 , and z

16:

end if

17:

Set Xti ← Xti−n + Y

18: 19:

end for end for

4. Option Pricing in L´ evy Models In this section, we consider option pricing problems when the underlying asset price dynamics is governed by tempered stable L´evy processes. Poirot and Tankov (2006) show that under an appropriate equivalent probability measure a tempered stable process becomes a stable process. This provides a fast Monte Carlo algorithm for computing the expectation of European contingent

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

22

claims under a tempered stable process. This method, however, does not allow access to the entire trajectory of the process, and thus we need a different sample-path generation method to deal with path-dependent financial contracts. More specifically, in Section 4.1, we extend the adaptive sampling scheme developed in Becker (2010) to finite variation tempered stable processes for pricing path-dependent options. This leads to remarkable efficiency gains in terms of computational costs. Significant amount of variance reduction is observed by stratified sampling in Section 4.2.

4.1. Adaptive bridge sampling under finite variation tempered stable (CGMY) processes The asset price is assumed to follow a geometric L´evy process: h i St = S0 exp (w + r)t + Xt where r is a risk free interest rate and w is the constant chosen so that the discounted value of a portfolio invested in the asset is a martingale. One popular model of asset price dynamics based on L´evy models is two-sided tempered stable processes associated with TS(α+ , λ+ , c+ ; α− , λ− , c− ) in Definition 3. Such a process, say Xt , can be written as the difference of two independent L´evy processes, that is, Xt = Xt+ − Xt− where Xt+ and Xt− are tempered stable subordinators associated with TS (α+ , λ+ , c+ ) and TS (α− , λ− , c− ), respectively. This includes the widely known finite-variation CGMY processes (Carr et al. 2002). In the six-parameter setting, it is further constrained that α+ = α− and c+ = c− . Avramidis and L’Ecuyer (2006) develop efficient Monte Carlo algorithms for pricing pathdependent options under the VG model, based on the representation of a VG process as the difference of two increasing gamma processes. (Hence the name difference-of-gamma bridge sampling or simply DGBS.) The authors obtain a pair of estimators whose expectations are monotone in increasing resolution, and the true option price lies in between. When unbiased price estimates are expensive to obtain, this pair provides upper and lower bounds without having to simulate all the process trajectories.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

23

In principle, as described in Section 6 of Avramidis and L’Ecuyer (2006), this method is extendable to option-pricing models driven by any L´evy process whose paths have finite variation. However, it is noted as a major difficulty for such an extension that there has been no bridge sampling algorithm for a more general class of L´evy processes such as CGMY processes. In this sense, we may extend the algorithm DGBS to Xt based on our bridge sampling scheme. A direct extension, however, seems undesirable since BS loses its relative merits with respect to HSQ as the monitoring step size decreases to zero. Nevertheless, we tested this idea and confirmed the accuracy of BS through the pricing of path-dependent options described below. We do not report the results in the paper to economize on space. Although such a direct extension of DGBS does not yield any relative advantages, bridge sampling is ideally suited to adaptive versions of DGBS that can lead to large savings in simulation costs. Becker (2010) has developed such an adaptive algorithm for VG processes, and we now extend the algorithm to CGMY processes. The key idea of adaptive sampling is to exclude the intervals from a given partition of the time dimension that neither contribute to the minimum nor to the maximum of a sample path. The author then presents adaptive sampling procedures for generating the infimum and supremum of VG processes and for pricing of lookback and barrier options described in Figures 2 and 3 in Becker (2010). Numerical studies show considerable reduction in computational efforts and memory requirements by reducing the number of bridge sampling steps. In this subsection, we adopt this adaptive DGBS and extend the whole idea to the option pricing under CGMY processes. To be specific, let Xt = Xt+ − Xt− be a CGMY process where Xt+ ∼ TS(Y, G, C) and Xt− ∼ TS(Y, M, C) are independent tempered stable L´evy subordinators, h i n o and let µ = w + r where w = − ln E exp(X1 ) = −CΓ(−Y ) (G − 1)Y − GY + (M + 1)Y − M Y . For t1 < t < t2 , we can see that L , µt1 + Xt+1 − Xt−2 + (t2 − t1 )µ− ≤ µt + Xt ≤ µt1 + Xt+2 − Xt−1 + (t2 − t1 )µ+ , U by the monotonicity of the subordinators with µ− = min(0, µ) and µ+ = max(0, µ). If T is the current set of grid points during the execution of the sampling procedure, we clearly have that

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

24

maxt∈T Xt ≤ maxt∈[0,T ] Xt and mint∈T Xt ≥ mint∈[0,T ] Xt . Then, for a particular interval [t1 , t2 ] of the current partition, if mint∈T Xt ≤ L (≤) U ≤ maxt∈T Xt holds, the interval neither contains the minimum nor the maximum of the path. Therefore, one can exclude the interval from further subdivisions. In what follows, we demonstrate numerical experiments for lookback and barrier options whose discounted payoffs are given as follows with maturity T and the monitoring interval size ∆ = T /d: h  i 1. floating strike lookback call VL = e−rT E ST − max S0 , S∆ , . . . , Sd∆ ; h i  2. up-and-in call barrier VB = e−rT E max 0, ST − K 1{max{S0 ,S∆ ,...,Sd∆ }>B} . Table 1 reports the reference parameter set, denoted by (a), that is obtained by calibrating the model to the market prices of the options on the S&P 500 index for June 2, 2003 in Carr et al. (2005). In addition, a risk free interest rate r = 0.0548, initial stock price S0 = 100, strike K = 100 Table 1

The results of the calibration of CGMY model on SPX (Carr et al. 2005). Parameter set (a)

T

C

G

M

Y

1.0453 0.3251 3.7103 18.4460 0.6029

and barrier level B = 120 are used. In the parameter set (a), we take SSR algorithm as the representative of HSQ, since it works most efficiently. In fact, BS without tabulation cannot beat SSR if a fixed number of points per path are generated. However, the adaptive sampling scheme greatly reduces the number of simulated points by BS. For an illustration, see Figure 6. We plot the computing times of SSR, adaptive BS, and adaptive BS with tabulation, as the observation time d increases. The sample size for price estimation is fixed at 104 . The computing time of SSR increases exponentially, while that of adaptive BS grows slowly. Furthermore, adaptive BS with tabulation shows outstanding performances for larger d values. Here, upper bounds are tabulated with the step size 0.01 for z. The speed-up from adaptive sampling is more prominent for barrier options as shown in the right panel of Figure 6. On the other hand, there is the associated trade-off between bias and variance. In adaptive BS, biases are caused by saddlepoint approximation in BS and the minimum approximation due to

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

120

SSR adaptive BS adaptive BS w/ tabulation

100

Computing time for barrier option

Computing time for lookback option

120

80

60

40

20

0 6

Figure 6

25

7

8

9

10

11

12

SSR adaptive BS adaptive BS w/ tabulation

100

80

60

40

20

0 6

7

8

9

log2 (d)

log2 (d)

(i)

(ii)

10

11

12

Computing times (in seconds) for pricing (i) lookback and (ii) barrier options versus the observation times d with the sample size 104 under the parameter set (a).

the adaptive scheme. The true option prices are computed using 109 Monte Carlo simulation trials by SSR algorithm. This is used to check the biases generated from adaptive BS with and without tabulation at that particular maturity. Figure 7 illustrates the estimated mean squared errors (denoted by MSE, given by the sum of the squared standard error and the squared bias) versus the computing times for the sample size 104 , 105 and 106 under the parameter set (a) when d = 27 . Note that biases caused by adaptive BS does not show a marked difference and that adaptive BS dominates HSQ in terms of MSE. 4.2. Stratified sampling Stratification is often applied to Monte Carlo simulation in option pricing, yielding successful effective speed-ups in convergence and variance reduction. For example, the use of stratified sampling for the VG process and the NIG process with the gamma bridge and the IG bridge is studied in Ribeiro and Webber (2003, 2004). In this section, we discuss the construction of such stratified samples for the terminal values of two-sided tempered stable processes, apply it to the path-dependent options 1–2, and examine its effects on variance reduction. The final value XT of the tempered stable subordinator can also be easily generated by numerical Laplace inversion. A reliable and efficient technique is explained in Glasserman and Liu (2010) based

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

26

0.14

0.14

SSR Adaptive BS Adaptive BS w/ tabulation

0.1

0.1

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0 0

0.5

1

1.5

2

2.5

3

SSR Adaptive BS Adaptive BS w/ tabulation

0.12

MSE

MSE

0.12

0 −1

3.5

−0.5

0

log10 (time)

0.5

1

(i) Figure 7

1.5

2

2.5

3

log10(time)

(ii)

Estimated mean squared errors versus computation times (in seconds) for the (i) Lookback and (ii) Barrier options for the sample size 104 , 105 and 106 under the parameter set (a).

on the analysis of Abate and Whitt (1992). The sampling algorithm generates U = Unif[0, 1] and then finds x > 0 such that F (x) = U with F the distribution function of XT , which is pre-computed and tabulated following Abate and Whitt (1992). We may use spline or linear interpolation for x. In more detail, the distribution F with its Laplace transform ϕ is computed as M

F (x) ≈

hx 2 X sin(hkx) + Re(ϕ)(−ikh), π π k=1 k

h=

2π . x + u

Here, u is set equal to µXT + mσXT where µXT is the mean of XT , σXT is its standard deviation, and m is some integer suitably chosen to achieve sufficient accuracy (e.g. m = 10). With the help of the Laplace inversion method, sample paths with stratified final values can be generated as follows for a given number of simulation trials n and length of time steps d: 1. generate a stratified sample (ˆ ui , vˆi ), i = 1, · · · , n, from the unit hypercube of dimension 2; 2. set XT+ i = F+−1 (ˆ ui ) and XT− i = F−−1 (ˆ vi ) for each i = 1, · · · , n where F± is the distribution of XT± respectively; 3. generate entire trajectories

n o Xt+j i

j=1,··· ,d−1

fied terminal values for each i by Algorithm 4.

and

n o Xt−j i

j=1,··· ,d−1

conditioned on the strati-

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

27

In Step 1, let K1 and K2 be the numbers of intervals of equal length along each dimension and thus each stratum of the hypercube is of the form 

l1 − 1 l1 , K1 K1





l2 − 1 l2 , × K2 K2

 ,

lj ∈ {1, 2, · · · , Kj } for j = 1, 2.

The related sampling probabilities are 1/K1 K2 for each region. A random vector (ˆ ui , vˆi ) uniformly distributed over this stratum for each i is simply  ui , vˆi ) = (ˆ

l1 − 1 + U1 l2 − 1 + U2 , K1 K2



where U1 , U2 are independent uniform random variables in the unit interval. Now we present the numerical results of stratified sampling under the parameter set (a). The sample size is n and the number of strata K1 and K2 are both set to be 10. We fix d = 26 and the estimated standard deviations (STD) of the resulting price estimates turn out not to vary substantially over d. Figure 8 shows a significant amount of additional variance reduction for BS with terminal stratification. It appears that the estimated STD of all the path-dependent option values are greatly reduced. Note that, in Figure 8, we deliberately use the log sample size (not the computation time) on the horizontal axis and STD on the vertical axis. This is to identify the effect of stratified sampling for the additional variance reduction. Remark 8. Stratification can also be combined with sequential sampling by stratifying at an initial non-zero time point. However, we find from numerical tests that there is little gain from such initial stratification. This is anticipated since the most key features of the path in option pricing often depends on the asset values at the maturity. Although not reported, effective variance reductions are observed for different maturities, monitoring frequencies, barrier levels, and model parameters. The effect of stratification for normal tempered stable processes is discussed as well in EC.4. Remark 9. In EC.5 of the e-companion, we present another application of bridge sampling. In short, a hybrid method that benefits from both sequential and bridge sampling is proposed. It is

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

28

0.35

0.35

SSR Adaptive BS Adaptive BS − Stratification

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0

4

4.5

5

5.5

6

SSR Adaptive BS Adaptive BS − Stratification

0.3

STD

STD

0.3

0

4

4.5

log10 (sample size)

(i) Figure 8

5

5.5

6

log10 (sample size)

(ii)

Estimated standard deviation due to stratification versus the number of simulation trials for the (i) Lookback, and (ii) Barrier options under the parameter set (a).

applied to least square Monte Carlo methods for American option valuation under subordinated Brownian motions, in order to avoid large memory requirements. Numerical results are also reported in EC.5.

5. Conclusion We showed that the approximate PDF of a tempered stable subordinator conditioned on the two observed endpoints using saddlepoint methods gives an accurate bridge distribution for a wide range of parameters. We proposed an acceptance-rejection method to simulate tempered stable L´evy bridges with the gamma bridge and the IG bridge as proposal densities. Special properties of the process such as scaling and linear conditional moments were exploited to enable further speed-ups. The suggested bridge sampling algorithm was applied to the pricing of path-dependent options for a finite-variation CGMY process. Through the suggested adaptive sampling scheme, extensive numerical tests demonstrated the accuracy and effectiveness of bridge sampling. However, the computational burden for evaluating upper bounds in bridge sampling still remains although it can be alleviated via tabulation. Diverse applications of bridge sampling are possible. One particular example that we implemented was stratified sampling and we observed a large amount of efficiency gains and variance

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

29

reduction. In addition, one might combine bridge sampling with randomized quasi-Monte Carlo method. Optimal stopping time problems are also potentially related; in particular adaptive sampling via bridge sampling is plausible by choosing the next step size of the simulated process adaptively near the stopping boundary. There are other important research topics related to the proposed method. First, one can consider an extension to infinite-variation tempered stable processes with α > 1, for which only an approximate sequential sampling method exists. Second, it is potentially useful to have an approximate conditional density and the bridge sampling scheme for statistical inference of L´evy processes. One can investigate simulation-based solutions for missing data problems as done for diffusions, or possible applications to the Monte Carlo EM algorithm and Bayesian methods.

Acknowledgments The authors would like to thank Prof. Paul Glasserman for his helpful comments and the workshop participants of Stochastic Processes and their Statistics held in 2013, Okinawa, for their feedback. The authors are also very grateful for many valuable comments from two anonymous reviewers and the Associate Editor which helped them improve the manuscript greatly. K. Kim’s work was supported by the Basic Science Research Program through the National Research Foundation of Korea funded by the Ministry of Education (NRF-2014R1A1A2054868). The author names are in alphabetical order.

References Abate J, Whitt W (1992) The Fourier-series method for inverting transforms of probability distributions. Queueing Systems, 10(1):5–88. Asmussen S, Hobolth A (2012) Markov bridges, bisection and variance reduction. In Monte Carlo and Quasi-Monte Carlo Methods 2010, edited by L. Plaskota, H. Wo´zniakowski, Springer-Verlag, Berlin. Avramidis AN, L’Ecuyer P (2006) Efficient Monte Carlo and Quasi-Monte Carlo option pricing under the variance gamma model. Management Science, 52(12):1930–1944. Avramidis AN, L’Ecuyer P, Tremblay PA (2003) Efficient simulation of gamma and variance-gamma processes. In Proceedings of the 2003 Winter Simulation Conference, edited by S. Chick, P. J. S´anchez, D. Ferrin, and D. J. Morrice. IEEE Press, Piscataway, NJ, 319–326.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

30

Barndorff-Nielsen OE, Cox DR (1979) Edgeworth and saddlepoint approximations with statistical applications. Journal of the Royal Statistical Society, Series B, 41(3):279–312. Barndorff-Nielsen OE, Shephard N (2001) Normal modified stable processes. Theory of Probability and Mathematical Statistics, 65:1–19. Becker M (2010) Unbiased Monte Carlo valuation of lookback, swing and barrier options with continuous monitoring under variance gamma models. Journal of Computational Finance, 31(4):35–61. Beskos A, Papaspiliopoulos O, Roberts GO, Fearnhead P (2006)

Exact and computationally efficient

likelihood-based estimation for discretely observed diffusion processes. Journal of the Royal Statistical Society, Series B, 68(3):333–382. Bladt M, Sørensen M (2014) Simple simulation of diffusion bridges with application to likelihood inference for diffusions. Bernoulli, 20(2):645–675. Brody DC, Hughston LP, Macrina A (2007) Beyond hazard rates: A new framework for credit-risk modelling. In Advances in Mathematical Finance, edited by Elliot, R., Fu, M., Jarrow, R. and Yen, JuYi, Applied and Numerical Harmonic Analysis Series, Festschrift volume in honour of Dilip Madan, Birkh¨ auser/Springer. Brody DC, Hughston LP, Macrina A (2008) Information-based asset pricing. International Journal of Theoretical and Applied Finance, 11(1):107–142. Butler RW (2007) Saddlepoint Approximations with Applications. Cambridge University Press, Cambridge. Carr P, Geman H, Madan D, Yor M (2002) The fine structure of asset returns: An empirical investigation. Journal of Business, 75(2):305–333. Carr P, Geman H, Madan D, Yor M (2005) Pricing options on realized variance. Finance and Stochastics, 9(4):453–475. Carr P, Madan D (2009) Saddlepoint methods for option pricing. Journal of Computational Finance, 13(1): 49–61. Chambers JM, Mallows CL, Stuck BW (1976) A method for simulating stable random variables. Journal of the American Statistical Association, 71(354):340–344.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

31

Cont R, Tankov P (2004) Financial Modeling with Jump Processes. Chapman and Hall/CRC, Boca Raton, FL. Daniels HE (1954) Saddlepoint approximations in statistics. Annals of Mathematical Statistics, 25(4): 631–650. Daniels HE (1958) Discussion of “The regression analysis of binary sequences” by D. R. Cox. Journal of the Royal Statistical Society, Series B, 20(2):236–238. Dembo A, Deuschel JD, Duffie D (2004) Large portfolio losses. Finance and Stochastics, 8(1):3–16. Devroye L (2009) Random variate generation for exponentially and ploynomially tilted stable distributions. ACM Transactions on Modeling and Computer Simulation, 19(4):18:1–20. Durham GB, Gallant AR (2002) Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes. Journal of Business & Economic Statistics , 20(3):297–316. Figueroa-L´ opez JE, Houdr´e C (2009) Small-time expansions for the transition distributions of L´evy processes. Stochastic Processes and their Applications, 119(11):3862–3389. Figueroa-L´ opez JE, Tankov P (2014) Small-time asymptotics of stopped L´evy bridges and simulation schemes with controlled bias. Bernoulli, 20(3):1126–1164. Glasserman P (2003) Monte Carlo Methods in Financial Engineering. Springer-Verlag, New York. Glasserman P, Kim K (2008) Beta approximations for bridge sampling. In Proceedings of the 2008 Winter Simulation Conference, edited by Mason SJ, Hill RR, M¨onch L, Rose O, Jefferson T, Fowler JW. IEEE Press, Piscataway, NJ, 569–577. Glasserman P, Kim K (2009) Saddlepoint approximations for affine jump-diffusion models. Journal of Economic Dynamics & Control, 33(1):15–36. Glasserman P, Liu Z (2010) Sensitivity estimates from characteristic funtions. Operations Research, 58(6): 1611–1623. Gordy MB (2002) Saddlepoint approximation of CreditRisk+ . Journal of Banking & Finance, 26(7):1335– 1353. Grabchak M (2012) On a new class of tempered stable distributions: moments and regular variation. Journal of Applied Probability, 49(4):1015–1035.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

32

Hirsa A, Madan D (2004) Pricing American options under variance gamma. Journal of Computational Finance, 7(2):63–80. Hofert, M (2011) Sampling exponentially tilted stable distributions. ACM Transactions on Modeling and Computer Simulation, 22(1):3:1–11. Hoyle E, Hughston LP, Macrina A (2011) L´evy random bridges and the modelling of financial information. Stochastic Processes and their Applications, 121(4):856–884. Jensen JL (1995) Saddlepoint Approximations. Oxford University Press, Oxford. Kaishev V, Dimitrova D (2009) Dirichlet bridge sampling for the variance gamma process: Pricing pathdependent options. Management Science, 55(3):483–496. Kawai R, Masuda H (2011) On simulation of tempered stable random variates. Journal of Computational and Applied Mathematics, 235(8):2873–2887. K¨ uchler U, Tappe S (2011) Tempered stable distributions and applications to financial mathematics. Working paper. Lin M, Chen R, Mykland P (2010) On generating Monte Carlo samples of continuous diffusion bridges. Journal of the American Statistical Association, 105(490):820–838. Lugannani R, Rice S (1980) Saddlepoint approximation for the distribution of the sum of independent random variables. Advances in Applied Probability, 12(2):475–490. Michael J, Schucany W, Haas R (1976) Generating random variates using transformations with multiple roots. The American Statistician, 30(2):88–90. Pitman J, Yor M (1982) A decomposition of Bessel bridges. Zeitschrift f¨ ur Wahrscheinlichkeitstheorie und Verwandte Gebiete, 59(4):425–457. Poirot J, Tankov P (2006) Monte Carlo option pricing for tempered stable (CGMY) processes. Asia-Pacific Financial Markets, 13(4):327–344. Ribeiro C, Webber N (2003) A Monte Carlo method for the normal inverse Gaussian option valuation model using an inverse Gaussian bridge. Working paper. Ribeiro C, Webber N (2004) Valuing path dependent options in the variance-gamma model by Monte Carlo with a gamma bridge. Journal of Computational Finance, 7(2):81–100.

Kim and Kim: Simulation of tempered stable L´ evy bridges Article submitted to Operations Research; manuscript no. (Please, provide the manuscript number!)

Roberts GO, Stramer O (2001)

33

On inference for partially observed nonlinear diffusion models using

Metropolis-Hastings algorithms. Biometrika, 88(3):603–621. Rogers LCG, Zane O (1999) Saddlepoint approximations to option prices. Annals of Applied Probability, 9(2):493–503. Rosi´ nski J (2007) Tempering stable processes. Stochastic Processes and their Applications , 117(6):677–707. R¨ uschendorf L, Woerner J (2002) Expansion of transition distributions of L´evy processes in small time. Bernoulli, 8(1):81–96. Sato K (1999) L´evy Processes and Infinitely Divisible Distributions. Cambridge University Press, Cambridge. Skovgaard IM (1987) Saddlepoint expansions for conditional distributions. Journal of Applied Probability, 24(4):875–887. Yang J, Hurd TR, Zhang X (2006) Saddlepoint approximation method for pricing CDOs. Journal of Computational Finance, 10(1):1-20.

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

ec1

E-Companion of the paper titled “Simulation of tempered stable L´ evy bridges and its applications” This E-Companion consists of four sections. EC.1 provides all proofs of the results developed in Section 3 including the proof of the main results, Theorem 1 and Proposition 3. EC.2 is then devoted to supplementary results followed from Section 3. In EC.3, the performance of the proposed bridge sampling method is compared to that of the existing sequential methods when generating a fixed number of observations or a skeleton of a tempered stable subordinator. EC.4 provides additional numerical results for path-dependent option pricing as well as terminal stratification under normal tempered stable processes. Finally, EC.5 illustrates a hybrid sampling scheme applied to least-square Monte Carlo method for American option pricing.

EC.1. Proofs of the results in Section 3 EC.1.1. Proof of Proposition 1 By the scaling property of Xt , we have the relationship fTS(α,λ,ct) (x) = rfTS(α,λ/r,crα t) (rx). Then, for s < t < T , h i h i P Xt ∈ dy |XT = z, Xs = x = P Xt−s (α, λ, c) ∈ d(y − x)|XT −s (α, λ, c) = z − x fTS(α,λ,c(t−s)) (y − x)fTS(α,λ,c(T −t)) (z − y) dy fTS(α,λ,c(T −s)) (z − x) rfTS(α,λ/r,crα (t−s)) (r(y − x))rfTS(α,λ/r,crα (T −t)) (r(z − y)) = dy rfTS(α,λ/r,crα (T −s)) (r(z − x))   1 1 = P Xrα (t−s) (α, λ/r, c) ∈ d(y − x) Xrα (T −s) (α, λ/r, c) = z − x r r   1 1 1 = P Xrα t (α, λ/r, c) ∈ dy Xrα T (α, λ/r, c) = z, Xrα s (α, λ/r, c) = x . r r r =

 EC.1.2. Proof of Lemma 1 For a tempered stable L´evy subordinator (Xt )t≥0 , let fXt1 |Xt2 (x|z) be the PDF of Xt1 conditioned on Xt2 = z, where t2 > t1 ≥ 0. It is straightforward to get fXt1 |Xt2 (x|z) =

fTS(α,λ,ct1 ) (x)fTS(α,λ,c(t2 −t1 )) (z − x) fS(α,ct1 ) (x)fS(α,c(t2 −t1 )) (z − x) = fTS(α,λ,ct2 ) (z) fS(α,ct2 ) (z)

ec2

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

due to the relation (5).



EC.1.3. Multi-dimensional stable processes and Proof of Proposition 2 There are several equivalent definitions of stable distribution, see Samorodnitsky and Taqqu (1994) e.g., and here we provide one with the CHF for its use. Definition EC.1. A random variable X is said to have a stable distribution if there are parameters 0 < α < 2, σ ≥ 0, −1 ≤ β ≤ 1, and µ ∈ R such that its CHF has the following form: 
    exp −σ α |u|α 1 − iβ sign(u) tan( πα ) + iµu , if α 6= 1; 2 φX (u) = 
  exp −σ |u| 1 + iβ 2 sign(u) ln |u| + iµu , if α = 1. π

(EC.1)

A multivariate stable distribution, as the distribution of a stable random vector, is described in terms of a finite measure Γ on the unit sphere of Rd and a shift vector µ0 . Definition EC.2. A d-dimensional random vector X = (X1 , · · ·, Xd ) with 0 < α < 2 is said to be α-stable if there exists a finite measure Γ on the unit sphere of Rd and a vector µ0 ∈ Rd such that  i h R
 πα 0 α  exp − ) Γ(ds) + i(u, µ ) , if α 6= 1; | (u, s) | 1 − i sign(u, s) tan( 2 Sd φX (u) = h i   exp − R |(u, s)| 1 + i 2 sign(u, s) ln |(u, s)|
Γ(ds) + i(u, µ0 ) , if α = 1 π S d

where (u, s) denote a inner product. Here, Γ is called a spectral measure and the pair (Γ, µ0 ) is unique, called a spectral representation. For a two-dimensional stable random vector X = (X1 , X2 ) with µ0 = 0, called a strictly stable distribution, we have  Z    πα   φX (u, v) = exp − |us1 + vs2 |α 1 − i tan sign(us1 + vs2 ) Γ(ds) , 2 S2

(EC.2)

where s = (s1 , s2 ). Each component Xi of X has a marginal stable distribution with the CHF (EC.1), where µ = 0, Z

α

|si | Γ(ds)

σ = σi =

1/α ,

S2

and 1 β = βi = α σi

Z

|si |α sign(si )Γ(ds). S2

ec3

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

In Hardin et al. (1991), the authors give conditions for the existence of the conditional moment     E |X2 |p |X1 = x for p > α and obtain an explicit form of E X2 |X1 = x as a function of x, which is given in the next theorem for a non-symmetric stable subordinator. Theorem EC.1 (Hardin et al. (1991)). Let (X1 , X2 ) be an α-stable process with 0 < α < 1 and a spectral measure Γ(ds) on the unit circle S2 in R2 . If it holds that, for some ν > 1 − α, Z S2

Γ(ds) < ∞, |s1 |ν

(EC.3)

then we have h

i

Z

E X2 |X1 = x =

s2 |s1 |α−1 sign(s1 )Γ(ds)

S2

x . σ1α

(EC.4)

For the proof of Proposition 2, let Xt be a stable subordinator associated with S(α, c) in Definition 1. Applying its L´evy measure (2), X1 has the CHF (EC.1) with µ = 0, σ = [−cΓ(−α) cos(πα/2)]1/α , and β = 1. Then we consider a two-dimensional stable subordinator X = (Xt2 , Xt1 ) with 0 < t1 < t2 where Xt = Xt (α, c) for 0 < α < 1. Its joint CHF is h i φX (u, v) = E exp i(uXt2 + vXt1 ) = φXt1 (u + v)φXt2 −t1 (u) n o h = exp − σ α |u + v |α t1 + |u|α (t2 − t1 ) n oi +iσ α tan(πα/2) sign(u + v)|u + v |α t1 + sign(u)|u|α (t2 − t1 ) .

(EC.5)

Matching the real part of (EC.5) with that of the spectral representation (EC.2) for X = (X1 , X2 ) = (Xt2 , Xt1 ), we have Z

n o |us1 + vs2 |α Γ(ds) = σ α |u + v |α t1 + |u|α (t2 − t1 ) .

(EC.6)

S2

First we show that X = (Xt2 , Xt1 ) satisfies the standard assumption (EC.3) in Theorem EC.1. This is verified from the following proposition. Proposition EC.1 (Samorodnitsky and Taqqu (1994)). The spectral measure Γα of an αstable random vector X is concentrated on a finite number of points on the unit sphere Sd if only if (X1 , · · ·, Xd ) can be expressed as a linear transformation of independent α-stable random variables, say X = AY, where Y1 , · · ·, Yd are independent α-stable and A is a d × d matrix.

ec4

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Since we have 









 Xt2  d  Xt1 + Xt2 −t1   1  = = Xt1 Xt1 1





1   Xt1    0 Xt2 −t1

by Proposition EC.1, the spectral measure Γ(ds) on S2 is of the form Γ(ds) = Pn with i=1 ai < ∞. In fact, we can specify the spectral measure as follows.

Pn

i=1 ai δ(si1 ,si2 ) (ds)

Lemma EC.1. The spectral measure Γ of X = (Xt2 , Xt1 ) is concentrated on two points (1, 1) and (1, 0). More explicitly, Γ = σ α t1 δ(1,1) + σ α (t2 − t1 )δ(1,0) . Proof.

Assume that Γ = a1 δ(1,1) + a2 δ(1,0) . Due to the uniqueness of the spectral measure, it is

enough to show that this choice of Γ is adequate. Equating the CHF (EC.5) to (EC.2) using our choice of Γ yields a1 = σ α t1 and a2 = σ α (t2 − t1 ).



By Lemma EC.1, the standard assumption in (EC.3) is satisfied. Thus, it remains to compute R the coefficient in Proposition 2. Let σ α (u, v) = S2 |us1 + vs2 |α Γ(ds). Then Z 1 ∂σ α (u, v) ∂ 1 · = |us1 + vs2 |α Γ(ds) α ∂v α S2 ∂v Z = |us1 + vs2 |α−1 s2 sign(us1 + vs2 )Γ(ds), S2

so that Z 1 ∂σ α (u, v) · = |s1 |α−1 s2 sign(s1 )Γ(ds) = σ1α a = σ α t2 a α ∂v S2 u=1,v=0 by (EC.4). The interchange of the order of integration and differentiation follows from the dominated convergence theorem since α |us1 + (v + h)s2 |α − |us1 + vs2 |α hs2 ≤ 1 1 + − 1 · |us1 + vs2 |α h h us1 + vs2 1 hs2 · |us1 + vs2 |α ≤ h us1 + vs2 = |s2 ||us1 + vs2 |α−1 for us1 + vs2 6= 0, which is integrable by Proposition EC.1. The second inequality uses the fact ||1 + z |α − 1| ≤ |z | for all z ∈ R. Then by (EC.6) " # n o 1 ∂σ α (u, v) 1 ∂ · = · σ α |u + v |α t1 + |u|α (t2 − t1 ) α ∂v α ∂v u=1,v=0

and this leads to a = t1 /t2 . This completes the proof.



= σ α t1 u=1,v=0

ec5

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

EC.1.4. Proof of Theorem 1 First, consider the case α ∈ (0, 1). Using the joint CGF (9), the saddlepoint equations K 0 (ˆ u, vˆ) = (z, x) and Ku0 (ˆ u0 , 0) = z have explicit solutions of the form   1 ! t2 − t1 1−α 1−α 1− u ˆ= , κ z−x 1 !   1−α   1 t1 1−α t2 − t1 1−α vˆ = − , and κ z−x x 1 !   1−α 1−α t2 u ˆ0 = . 1− κ z Therefore, plugging this solution into the formula (7), we have the following saddlepoint approximation fˆXt1 |Xt2 (x|z) of fXt1 |Xt2 (x|z): 1 fˆXt1 |Xt2 (x|z) ∝ √ 2π



× exp

h

|K 00 (ˆ u, vˆ)| 00 (ˆ |Kuu u0 , 0)|

−1/2



i K(ˆ u, vˆ) − u ˆz − vˆx − K(ˆ u0 , 0) − u ˆ0 z .

(EC.7)

Straightforward computations then yield the desired result. When α = 0, we also solve the saddlepoint equations with the joint CGF (9), which lead to u ˆ=κ−

t2 − t1 , κ(z − x)

vˆ =

t2 − t1 t1 − , κ(z − x) κx

and

u ˆ0 = κ −

A similar procedure using (EC.7) gives the expression in the statement.

t2 . κz



EC.1.5. Proof of Proposition 3 Let p = α/(1 − α). First, in the case of 0 < α < 1/2 or equivalently 0 < p < 1, the ratio can be written as fˆXt1 |Xt2 (x|z) fXGB (x|z) t |Xt 1

2

=

B1 (z; α, κ; t1 , t2 ) g1 (x), N (z; α, κ; t1 , t2 )

where    1+p   1 t1 t2 − t1 t1 (t2 − t1 ) 2 p + tκ2 t1+p 2 B1 = √ Beta , z2 exp κ κ t2 p(1 + p)κz p 2πκ and t − κ1 − p 2

g1 (x) = x

(z

t2 −t1 p − x)− κ − 2

 exp −

1 p(1 + p)κ



t1 1+p (t2 − t1 )1+p + xp (z − x)p

 .

ec6

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Finding the maximum of g1 (x) over (0, z) is equivalent to finding the minimum of h1 (y) over (0, 1) where h1 (y) = a log y + b log(1 − y) + cy −p + d(1 − y)−p + e with t2 − t1 p b= + , κ 2

t1 p a= + , κ 2

t1+p 1 c= , p(1 + p)κz p

(t2 − t1 )1+p d= , p(1 + p)κz p

 e=

 t2 + p log z. κ

Note that   p/2 p/2 √ 1−y y cy −p + d(1 − y)−p = c + d ≥ 2 cd, (y(1 − y))−p/2 y 1−y from which we have √ min h1 (y) ≥ (a ∧ b) log y(1 − y) + 2 cd(y(1 − y))−p/2 + e y∈(0,1)    √    √ 2/p   2(a∧b) log p cd + 1 + e, if p cd < 14 ; p a∧b a∧b ≥   (a ∧ b) log( 1 ) + √cd 2p+1 + e, otherwise. 4

Second, in the case of 1/2 < α < 1 or equivalently p > 1, we need to find the maximum of g2 (x) over (0, z) where g2 (x) = x

− p−1 2

(z − x)

− p−1 2

 exp −

1 p(1 + p)κ



t1 1+p (t2 − t1 )1+p + xp (z − x)p



1 + 2κ



t21 (t2 − t1 )2 + x z−x

 .

Similarly as above, we consider the minimum of h2 (y) over (0, 1) where h2 (y) = a log y(1 − y) + by −p + c(1 − y)−p − dy −1 − e(1 − y)−1 + f, with a=

p−1 , 2

b=

t1+p 1 , p(1 + p)κz p

c=

(t2 − t1 )1+p , p(1 + p)κz p

d=

t21 , 2κz

e=

(t2 − t1 )2 , 2κz

and lastly f = (p − 1) log z. Then, we observe that there exists a small 0 > 0 such that h2 (y) = a log y(1 − y) + O(y −p ) + O((1 − y)−p ) on [0 , 1 − 0 ]c because y −p  1/y and (1 − y)−p  1/(1 − y) for y values near 0 and near 1, respectively. This implies that h2 (y) is bounded below on [0 , 1 − 0 ]c . On the other hand, the function is bounded and continuous on [0 , 1 − 0 ]. Hence, the desired result follows.



ec7

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

EC.2. Supplementary results in Section 3 EC.2.1. An extension of Theorem 1 The corollary below is an extension of 1 to the multivariate random vector (Xt1 , · · ·, Xtn−1 ) given Xtn = xn where 0 < t1 < · · · < tn . Its proof follows in a straightforward fashion from the fact that its multi-dimensional CGF K(u1 , · · · , un ) is given by n n X X (tk − tk−1 ) l ui K(u1 , · · · , un ) = k=1

! with

t0 = 0.

i=k

Corollary EC.1. Let Xt be a tempered stable L´evy subordinator associated with TS(α, κ), α ∈ [0, 1), and κ > 0. Then the saddlepoint approximation fˆXt1 ,···,Xtn−1 |Xtn (x1 , · · ·, xn−1 |xn ) for the conditional PDF fXt1 ,···,Xtn−1 |Xtn (x1 , · · ·, xn−1 |xn ) conditioned on Xtn = xn , where 0 < t1 < · · · < tn , is as follows: fˆXt1 ,···,Xtn−1 |Xtn (x1 , · · ·, xn−1 |xn ) 2−α  Qn  2(1−α)  1  (tk − tk−1 ) 2(1−α) xn − n−1 k=1 Qn ∝ (2πκ) 2 tn k=1 (xk − xk−1 ) " ( n )# 1 1 (1 − α)2 X (tk − tk−1 ) 1−α tn 1−α × exp − − α α 1−α ακ (x − x ) xn 1−α k k−1 k=1 for α ∈ (0, 1), and n−1 fˆXt1 ,···,Xtn−1 |Xtn (x1 , · · ·, xn−1 |xn ) ∝ (2πκ)− 2 Q

×

n Y

tn +1 2

tnκ n k=1 (tk

(xk − xk−1 )

− tk−1 )

tk −tk−1 −1 κ

tk −tk−1 +1 κ 2 tn

xn 1− κ

k=1

for α = 0, where x0 = 0 and t0 = 0. EC.2.2. Non-Gaussian-based saddlepoint approximation The saddlepoint expansions which involve the PDF φ and the CDF Φ of a standard normal random variable are said to be Gaussian-based. However, it is known that a greater accuracy might be obtained by saddlepoint expansions with a different distributional base rather than the standard normal. In fact, non-Gaussian-based saddlepoint approximations, developed in Wood et al. (1993),

ec8

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

are known to attain better precision for some heavy-tailed distributions for which the Gaussianbased method does not work properly in certain settings. In what follows, we present an IG-based saddlepoint approximation for tempered stable processes. Although this is a reasonable choice given the known density form (12) at α = 1/2, the IG-based conditional density (EC.12) does not stand out as being more efficient. The average and maximum errors are almost the same as those of the Gaussian-based one, while the latter allows a more explicit form and handy algorithms. Saddlepoint approximations that have involved a non-normal density base is mostly developed in Wood et al. (1993). Consider the target random variable X with a PDF f (x) and a CGF K(u), and suppose that the base function has a PDF λ(η) and a CGF L(s). Then the (λ, L)-based saddlepoint density approximation f˜(x) of f (x) is given by s f˜(x) = λ(ˆ η)

L00 (ˆ s(ˆ η )) , K 00 (ˆ u)

(EC.8)

where sˆ(η) is the saddlepoint root of L0 (ˆ s) = η and u ˆ(x) is the root of K 0 (ˆ u) = x. The value ηˆ is an optimally tilted η value chosen so that ηˆ solves L(ˆ s(η)) − sˆ(η)η = K(ˆ u) − u ˆx

(EC.9)

in η. The mapping x 7→ ηˆ determined by the solution in (EC.9) is a well-defined smooth bijection. The right side of (EC.9) is a function of x that implicitly depends on x through u ˆ that solves K 0 (ˆ u) = x. The maximum of K(ˆ u) − u ˆx is 0 and occurs when u ˆ = 0 and x = E[X]. The mapped value for E[X] is ηˆ = L0 (0) that occurs when 0 = sˆ(ˆ η ). In all the other cases, K(ˆ u) − u ˆx < 0 and each negative value is associated with two values of x, say x− and x+ , for which x− < E[X] < x+ . The mapping associates x− 7→ ηˆ− and x+ 7→ ηˆ+ where ηˆ− < L0 (0) < ηˆ+ and ηˆ− , ηˆ+ are two solutions to (EC.9). We refer the reader to Butler (2007) for more details. For our application, let the base function be a two-parameter family of an IG distribution IG(µ, σ 2 ) with its density function h i µ2/3 −3/2 µ λ(x; µ, σ 2 ) = √ x exp − 2 (x − µ)2 . 2σ x 2πσ

ec9

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

√ This is in fact TS(1/2, λ, c) with λ = µ/σ 2 and c = µ3/2 / 2πσ in our parameter setting. The

IG distribution is the most prominent among the base functions for heavy-tailed distributions. Moreover, it is also a special case of a tempered stable distribution that has an exact closed formula for the PDF, at least when α = 1/2. We follow the parameter setting and formulas in Chapter 16, Butler (2007) and in Wood et al. (1993). The scale invariance below enables us to reduce one parameter. Proposition EC.2 (Wood et al. (1993)). Suppose the base distribution Λa,b (z) is a location and scale transformed CDF of Λ(z), i.e. Λa,b (z) = Λ((z − a)/b). Then, provided b > 0, the resulting formulas do not depend on a or b. As in Butler (2007), let µ = θ and σ 2 = θ3 and then the PDF and CGF become   1 1 −3/2 2 exp − 2 (η − θ) , λ(η; θ) = √ η 2θ η 2π

1 L(s) = − θ

r

1 − 2s. θ2

With K(u) = tl(u) in (8), matching the standardized cumulants for the base function with the target density gives the optimal choice of the parameter θ written as ˆ x) = θ(t,

{K 000 (ˆ u)}

s

2

3 {K 00 (ˆ u)}

3+ω ˆ (t, x)

3

2

!−1

3

{K 00 (ˆ u)} !−1 α √ (2 − α) κx 2(1−α) 3+ω ˆ (t, x) 1 (1 − α)t 2(1−α)

α

=

{K 000 (ˆ u)}

(2 − α)2 κx 1−α 1

(1 − α)2 t 1−α

(EC.10)

with v 1 !u   1−α u 2(1 − α) t t ω ˆ (t, x) = sign 1 − x κ

1

t 1 − α t 1−α x− + α α α x 1−α

! .

ˆ the root of (EC.9), ηˆ, can be expressed as For this optimal choice of θ, 2 ˆ ˆ x) + θ(t, x) ηˆ(t, x) = θ(t, 2

s 2

ω ˆ (t, x) + ω ˆ (t, x) ω ˆ (t, x)2 +

and the saddlepoint sˆ for the CGF L(s) is 1 sˆ (ˆ η (t, x)) = 2

1

1 − 2 ˆ x) ηˆ(t, x)2 θ(t,

! .

4 ˆ x) θ(t,

! ,

(EC.11)

ec10

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Then, the IG-based saddlepoint approximation f˜Xt1 |Xt2 (x|z) for the conditional PDF fXt1 |Xt2 (x|z) for 0 < t1 < t2 can be computed from (EC.8): 1 f˜Xt1 |Xt2 (x|z) ∝ √ 2πκ



t1 (t2 − t1 ) t2

1  2(1−α) 

z (z − x)x

2−α  2(1−α)

  exp −

 2  ˆ 1 , x)  η ˆ (t , x) − θ(t 1 1

ˆ 1 , x)2 ηˆ(t1 , x) θ(t 2   2  ˆ 2 , z)  ˆ 2 − t1 , z − x)  ηˆ(t2 , z) − θ(t ηˆ(t2 − t1 , z − x) − θ(t  − + (EC.12)  ˆ 2 − t1 , z − x)2 ηˆ(t2 − t1 , z − x) ˆ 2 , z)2 ηˆ(t2 , z)  θ(t θ(t  2 

ˆ x) and ηˆ(t, x) are as in (EC.10) and (EC.11), respectively. where θ(t,

EC.3. Numerical performance in comparison to existing sequential methods In this section, we provide numerical studies to examine the efficiency of [TSLB] algorithm by comparing it with HSQ. To specify which algorithm is used for the hybrid scheme, we here classify HSQ into three different algorithms: the simple stable rejection (SSR) described in Algorithm 1, Hofert’s algorithm in Hofert (2011), and Devroye’s algorithm in Devroye (2009). SSR and Hofert’s algorithm perform better as ∆(1 − α)/κ/α gets smaller. (Recall our parametrization λ = (1 − α)/κ.) If ∆(1 − α)/κ/α < 1, SSR is the fastest algorithm and Hofert’s algorithm becomes identical to SSR. On the other hand, if ∆(1 − α)/κ/α > 1, Hofert’s and Devroye’s algorithms are competitors to each other. In Hofert (2011), the author suggests a simple guideline; use Devroye’s algorithm if ∆(1 − α)/κ/α > 2.7 and Hofert’s otherwise. As an illustration, we compare computational times to generate 104 number of trajectories (X∆ , · · · , X8∆ ) with ∆ = 1 and varying α and κ in Figure EC.1. The fastest algorithm in each parameter set (α, κ) is colored in the figure among the four algorithms: SSR, Hofert, Devroye, and BS with tabulation. When α = 0.5, Xt becomes an inverse Gaussian process and the exact bridge sampling (IGB) is implemented. To provide a guideline for deciding when to apply which, the complexity of Hofert is written in each cell corresponding to α and κ. Here, in our numerical experiments, we adopt Devroye’s algorithm for generating the terminal value in Algorithm 4 since it performs best for a large time step. We first note that, due to the extra burden of computing bounding constants in BS without tabulation, this method does not compete well with Hofert’s or Devroye’s. However, BS with

ec11

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

tabulation is superior as shown in most of the cells in the upper left part. More specifically, BS with tabulation performs best where the complexity of Hofert is greater than approximately 1.5 even though Devroye’s algorithm is fastest for parameter sets with very high complexities. In the remaining cells, SSR usually dominates all the others. Although not presented in the figure, Hofert’s algorithm works more efficiently than Devroye’s algorithm when its complexity is less than approximately 2.3 in our experiments. In addition, the computation times of Hofert, Devroye, BS without tabulation, and BS with tabulation are reported in Figure EC.1 for the borderline cells indicated by ∗. 𝛼\𝜅

0.1

0.2

0.3

0.4

0.5

0.6

0.9

1

2

0.05 190

95

63.3 47.5

38

31.7 27.1 23.8 21.1

19

9.5

0.1

45

9

4.5

90

30

22.5

18

15

0.15 56.7 28.3 18.9 14.2 11.3 9.44 0.2 0.25 0.3

40

20

13.3

10

8

30

15

10

7.5

6

0.7

3

4

6.33 4.75 3

2.25

5

6

7

8

9

10

3.8

3.17 2.71 2.38 2.11 1.9*

1.8

1.5* 1.29 1.13

12.9 11.3

10

8.1

7.08

6.3

5

4.44

4

2

1.33*

1

0.8

0.67 0.57

4.29 3.75 3.33

3

1.5*

1

0.75

0.6

0.5

6.67 5.71 5

0.8

1

𝛼

Hofert

Devroye

BS

BS with tabulation

0.05

0.36

0.38

1.62

0.28

0.1

0.41

0.43

1.68

0.3

0.15

0.36

0.45

1.93

0.25 0.24

0.9

5.67 2.83 1.89 1.42* 1.13 0.94 0.81 0.71 0.63 0.57 0.5

Computation times at the border *

0.44

0.4

0.43 0.38 0.33

0.3

23.3 11.7 7.78 5.83 4.67 3.89 3.33 2.92 2.59 2.33* 1.17 0.78 0.58 0.47 0.39 0.33 0.29 0.26 0.23

0.2

0.33

0.47

2.23

0.25

0.35

0.49

2.51

0.2

0.35 18.6 9.29 6.19 4.64 3.71

3.1

2.65 2.32 2.06 1.86* 0.93 0.62 0.46 0.37 0.31 0.27 0.23 0.21 0.19

0.3

0.53

0.54

2.57

0.21

0.4

2.5

2.14 1.88 1.67 1.5* 0.75

0.25 0.21 0.19 0.17 0.15

0.35

0.41

0.53

2.64

0.23

0.45 12.2 6.11 4.07 3.06 2.44 2.04 1.75 1.53 1.36* 1.22 0.61 0.41 0.31 0.24

15

0.2

0.4

0.38

0.52

2.61

0.25

0.5

0.17 0.14 0.13 0.11

0.45

0.35

0.5

2.87

0.29

0.55

0.27

0.48

3.65

0.23

0.6

0.33

0.51

3.7

0.22

0.65

0.33

0.5

3.62

0.21

0.7

0.35

0.5

3.6

0.21

0.75

0.48

0.51

3.71

0.23

0.8

0.59

0.54

4.11

0.22

0.18 0.09 0.06 0.04 0.04 0.03 0.03 0.02 0.02 0.02

0.85

0.41

0.48

4.25

0.23

0.9 1.11* 0.56 0.37 0.28 0.22 0.19 0.16 0.14 0.12 0.11 0.06 0.04 0.03 0.02 0.02 0.02 0.01 0.01 0.01

0.9

0.24

0.41

4.38

0.26

0.95 0.53* 0.26 0.18 0.13 0.11 0.09 0.08 0.07 0.06 0.05 0.03 0.02 0.01 0.01 0.01 0.01 0.01 0.01 0.01

0.95

0.15

0.29

6.98

0.26

10

7.5

5

5

3.33

3.75

2.5

3

2

1.67 1.43 1.25 1.11

1

0.5

0.5

0.33 0.25

0.55 8.18 4.09 2.72 2.05 1.64 1.36 1.17* 1.02 0.91 0.82 0.41 0.27 0.6

0.2

0.3

0.2

0.9

0.77 0.67

0.6

4.29 2.14 1.43* 1.07 0.86 0.71 0.61 0.54 0.48

0.17 0.15 0.14 0.12

0.16 0.14 0.12

6.67 3.34 2.22 1.67 1.33* 1.11 0.95 0.83 0.74 0.67 0.33 0.22 0.17 0.13 0.11

0.65 5.38 2.69 1.79 1.35* 1.08 0.7

0.38

0.1

0.1

0.1

0.09 0.08

0.08 0.07 0.07

0.54 0.27 0.18 0.13 0.11 0.09 0.08 0.07 0.06 0.05 0.3

0.21 0.14 0.11 0.09 0.07 0.06 0.05 0.05 0.04

0.75 3.34 1.67* 1.11 0.83 0.67 0.56 0.48 0.42 0.37 0.33 0.17 0.11 0.08 0.07 0.06 0.05 0.04 0.04 0.03 0.8

2.5* 1.25 0.83 0.63

0.5

0.42 0.36 0.31 0.27 0.25 0.13 0.08 0.06 0.05 0.04 0.04 0.03 0.03 0.03

0.85 1.76* 0.88 0.59 0.44 0.35 0.29 0.25 0.22

SSR

Hofert

Figure EC.1

Devroye

BS

0.2

BS with tabulation

The fastest algorithm for sampling 104 number of trajectories (X∆ , · · · , X8∆ ) among SSR, Hofert, Devroye, and BS with tabulation, where the number in each cell is the complexity of Hofert. The computation times (in seconds) in the borderline cells, indicated as ∗, are presented.

In this experiment, a bound table is generated offline by parameterizing z where the length of z vector is set to be 300, and a bound for each realized z value is obtained by linear interpolation. Since we bisect time points, it suffices to compute a one-dimensional table for z and we exploit the scaling property (14). The computing time for off-line tabulation is not included in our comparison.

ec12

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

This is because the tabulation time becomes relatively negligible as the number of sample paths increases and the computing time is approximately 0.1 seconds only. All the tests are implemented in C language (visual studio premium 2013) on a desktop computer with processor Intel Core i5 CPU 760. The accuracy of BS is numerically checked for the parameter sets under consideration by comparing the empirical distributions and Q-Q plots generated by the four algorithms. We note here that the graphs are visually indistinguishable among the four algorithms for all the parameter sets examined. In our applications, we examine the accuracy of BS via pricing of complex pathdependent options. In conclusion, BS algorithm with tabulation is an attractive candidate for generating a skeleton of a tempered stable subordinator. Both speed and accuracy of BS are comparable to that of the existing exact algorithms, HSQ. The specific decision rule may differ as ∆ and the number of time steps vary. However, at least in the examined conditions, BS with tabulation works more efficiently than the others when the complexity of Hofert is greater than approximately 1.5.

EC.4. Path-dependent option pricing under normal tempered stable processes In this section, we assume that Xt is a normal tempered stable process, which is a time-changed Brownian motion subordinated by a tempered stable subordinator. In other words, we define h i St = S0 exp (w + r)t + Xt ,

Xt = Xt (θ, σ, α, κ) = BYt

where Bt = θt + σWt , Wt is a standard Wiener process, and Yt is a tempered stable subordinator associated with TS(α, κ). Bridge sampling for this process, therefore, consists of two parts: one for Yti ’s using an algorithm such as Algorithm 4 and one for Xti ’s using the usual Brownian bridge sampling. For further information about Xt such as the characteristic exponent or the L´evy measure, we refer the reader to Cont and Tankov (2004) or Barndorff-Nielsen and Shephard (2001). In the model, w is found via the characteristic exponent so that   α  κ(θ + σ 2 /2) 1−α 1− 1− . w = − ln E exp(X1 ) = κα 1−α h

i

ec13

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

We consider three different path-dependent options with maturity T and the monitoring interval size ∆ = T /d: 1. floating strike lookback put h i  VL = e−rT E max S0 , S∆ , . . . , Sd∆ − Sd∆ ; 2. fixed strike Asian call h  i VA = e−rT E max 0, AT − K ,

d

AT =

1 X Sk∆ ; d + 1 k=0

3. up-and-in call barrier h i  VB = e−rT E max 0, Sd∆ − K 1{max{S0 ,S∆ ,...,Sd∆ }>B} . Two parameter sets (b) α = 0.3, κ = 0.01 and (c) α = 0.6, κ = 0.01 are chosen as representatives of the left edge of the parameter boxes in Figure EC.1. Other reference parameters for the normal tempered stable process are r = 0.0548, θ = −0.2859, and σ = 0.1927, which are borrowed from Hirsa and Madan (2004). Also, S0 = K = 100, B = 120, T = 2, and d = 8. As a benchmark, the true option prices at T = 2 are computed by SSR with the sample size 109 . Table EC.1 compares the prices of options 1–3 computed by Devroye’s algorithm, which is the fastest SQ under parameter sets considered, and by BS with tabulation. Here, n is again the number of simulation trials and computing times are reported in seconds. We note that a great amount of variance reduction can be achieved by terminal stratification. Generation of stratified samples at the final time step is similar to the procedure in Section 4.2. Again, variance reduction is most notable for barrier options. In addition, the variability in the options’ payoffs can be reduced by stratifying the terminal values of normal tempered stable processes. We can generate sample paths with stratified final values as follows for a given number of simulation trials n and length of time steps d: 1. generate a stratified sample (ˆ ui , vˆi ), i = 1, · · · , n, from the unit hypercube of dimension 2; 2. set YTi = F −1 (ˆ ui ) for each i = 1, · · · , n where F is the distribution of YT ;

ec14

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Table EC.1

Option prices, standard errors, and mean squared errors for the normal tempered stable process with the reference parameter sets (b) α = 0.3, κ = 0.01 and (c) α = 0.6, κ = 0.01.

Parameter set n

104

105

106

n

104

105

106

(b)

Lookback

Asian

Devroye BS-tabulation Devroye BS-tabulation

Barrier Devroye

BS-tabulation

Price

12.3817

12.3384

8.5391

8.5832

15.1380

15.1010

Std

0.1220

0.1218

0.1132

0.1131

0.2260

0.2241

Bias

0.0981

0.0548

0.0985

0.0544

0.0465

0.0835

RMSE

0.1565

0.1336

0.1501

0.1255

0.2307

0.2391

Time

0.326

0.355

0.346

0.351

0.334

0.368

Price

12.2518

12.2482

8.6352

8.6342

15.1930

15.1756

Std

0.0381

0.0380

0.0357

0.0358

0.0715

0.0714

Bias

0.0318

0.0354

0.0024

0.0034

0.0085

0.0089

RMSE

0.0497

0.0520

0.0359

0.0359

0.0720

0.0720

Time

3.203

3.262

3.132

3.322

3.239

3.528

Price

12.2761

12.2890

8.6437

8.6439

15.1923

15.1756

Std

0.0121

0.0121

0.0113

0.0113

0.0226

0.0226

Bias

0.0075

0.0054

0.0061

0.0063

0.0078

0.0089

RMSE

0.0142

0.0133

0.0129

0.0130

0.0239

0.0243

Time

33.633

35.312

33.754

35.426

33.742

35.341

BS-tabulation

Devroye

BS-tabulation

(c)

Devroye BS-tabulation Devroye

Price

12.3516

12.1396

8.5587

8.6092

15.0157

15.1626

Std

0.1206

0.1207

0.1146

0.1123

0.2270

0.2242

Bias

0.0699

0.1421

0.078

0.0274

0.1669

0.0199

RMSE

0.1394

0.1864

0.1386

0.1156

0.2818

0.2251

Time

0.3473

0.251

0.347

0.253

0.335

0.264

Price

12.2909

12.3018

8.6412

8.6184

15.1455

15.1750

Std

0.0381

0.0381

0.0359

0.0358

0.0715

0.0716

Bias

0.0092

0.0202

0.0046

0.0182

0.0370

0.0075

RMSE

0.0392

0.0431

0.0362

0.0402

0.0805

0.0720

Time

3.286

2.4415

3.365

2.543

3.315

2.471

Price

12.2706

12.2668

8.6425

8.6289

15.1941

15.1881

Std

0.0121

0.0121

0.0113

0.0113

0.0226

0.0226

Bias

0.0110

0.01482

0.0058

0.0078

0.0115

0.0055

RMSE

0.0164

0.0191

0.0128

0.0137

0.0253

0.0232

Time

33.2544

25.7301

33.654

25.801

33.679

25.728

ec15

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

3. set XTi = θYTi + σ

p YTi Z i where Z i = Φ−1 (ˆ vi ) for Φ standard normal cumulative distribution

function; n o 4. generate entire trajectories Ytij

j=1,··· ,d−1

n o at tj by Algorithm 4 and Xtij

j=1,··· ,d−1

at Ytij by

Brownian bridge sampling conditioned on the stratified terminal values for each i Taking the equiprobable strata K1 = K2 = 10, Figure EC.2 shows the significant amount of variance reduction for BS with terminal stratification. It appears that the estimated STDs of all the pathdependent option values are greatly reduced, and in particular, it is most significant for barrier option.

0.25

0.25 (b) Devroye (b) BS − Strat (c) Devroye (c) BS − Strat

0.2

0.25 (b) Devroye (b) BS − Strat (c) Devroye (c) BS − Strat

0.2

0.2

STD

0.15

STD

0.15

STD

0.15

0.1

0.1

0.1

0.05

0.05

0.05

0

4

5

log10 (sample size)

Figure EC.2

6

(b) Devroye (b) BS − Strat (c) Devroye (c) BS − Strat

0

4

5

log10 (sample size)

6

0

4

5

6

log10 (sample size)

Additional reduction of the estimated standard deviation due to stratification versus the number of simulation trials n for the Lookback, Asian, Barrier options from left to right under the parameter sets (b) and (c).

EC.5. American option pricing under normal tempered stable processes In this section, we consider the same asset model based on normal tempered stable processes as in the previous section. Before further discussion, we first present the backward generation algorithm of normal tempered stable processes in Algorithm 5 which we will utilize together with American option pricing. However, it is worth noting that dyadic generation is preferable in terms of efficiency and accuracy if applicable.

ec16

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Algorithm 5 Backward generation of a sample path of a normal tempered stable process via the Brownian bridge and [TSLB] 1: At given time points 0 = t0 < t1 < · · · < td , 2:

Y0 ← 0; X0 ← 0

Generate Ytd ∼ TS(α, (1 − α)/κ, td ((1 − α)/κ)(1−α) /Γ(1 − α)) and Z ∼ N(0, 1) p 4: Set Xtm ← θYtm + σ Ytm Z 3:

5: 6: 7:

for i ← 1 to d do Generate Ytd−i from [TSLB] with t1 = td−i , t2 = td−i+1 , and z = Yt2 .  (Y  −Yt )Yt t Generate Z ∼ N 0, d−i+1Yt d−i d−i σ 2 d−i+1

8: 9:

Set Xtd−i ←

Ytd−i Ytd−i+1

Xtd−i+1 + Z

end for

American option pricing is a challenging task due to its early exercise feature. It is often approximated by assuming that the option can be exercised at finitely many time points. Then at each time step, the intrinsic value, the payoff if exercised at that time step, and the continuation value, the conditional expected payoff if not exercised, are compared to determine an optimal early exercise policy. Such comparisons are done backward in time. The LSMC methodology developed by Longstaff and Schwartz (2001) is one of the most widely adopted methods to estimate the continuation value. The description of LSMC, which utilizes least square regression, is referred to the original paper of Longstaff and Schwartz (2001). The convergence of LSMC is rigorously proved in Stentoft (2004). From the computational point of view, one practical concern is a data management problem because the procedure requires to generate and store the n × d matrix of the simulated paths with sample size n and number of time steps d. Its purpose is to evaluate the continuation value, which is compared with the intrinsic value, in each step of a backward dynamic program. This often limits the accuracy of price estimates given computational capacity. However, the difficulty can be resolved if bridge sampling of the underlying asset process is available. In this case, the comparison

ec17

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

of the intrinsic value and the continuation value at each step can be performed by storing simulated data only at two time points. This is due to the fact that the LSMC needs only the current asset price to estimate the continuation value. Recently, Pellegrino and Sabino (2015) combine diffusion bridge construction with the LSMC, which enables a large size simulation without storing the whole trajectories of an underlying asset. We therefore adopt our bridge sampling to enhance the LSMC based American option pricing under normal tempered stable processes. However, it should be recalled that the accuracy of the proposed bridge sampling deteriorates when the time ratio is close to 0 or 1. This could cause problems in applying Algorithm 5 when the discrete exercise-time points are close together. Thus, we instead propose a hybrid sampling method that combines HSQ and BS in order to reduce the computational burden in implementing the LSMC. The central idea of this approach is to strike the balance between the unbiasedness of HSQ and the efficiency of BS. More precisely, for a given time grid (t1 , · · · , td ) and a computational budget, one generates and stores the n × ds matrix of the underlying paths at the time grid (tdb , t2db , · · · , tds db ) where ds = d/db is the number of time points to sample from HSQ. For each k = ds − 1 to 0, the intermediate values between tkdb and t(k+1)db are filled recursively by Algorithm 5. Therefore, the total memory requirement is reduced by approximately a factor of db . See Figure EC.3 for an illustration. Sequential sampling

t0

t db

∙∙∙ ∙∙∙

t 2 db

t( ds 1) db

t d s db

∙∙∙ Backward bridge sampling

Figure EC.3

The hybrid construction of a sample path at d time points with ds number of HSQ and db number of BS where d = ds × db .

We check the benefit of hybrid sampling in the following numerical example. The valuation of an American put option is considered with payoff max(0, K − ST ), strike K = 100 and maturity T = 1. Early exercise is possible only at 50 equally distributed time points including the maturity. The

ec18

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

basis functions in the LSMC are the first nine Laguerre polynomials. This setting is adopted from Longstaff and Schwartz (2001). Note that the number of basis functions is selected not to have rank deficiency when computing the inverse matrix in the procedure of LSMC. We continue to assume r = 0.0548, θ = −0.2859, and σ = 0.1927. We let the volatility of tempered stable subordinator be the same as that of the VG model, i.e., κ = 0.2505 (Avramidis and L’Ecuyer 2006, Kaishev and Dimitrova 2009). Lastly, α = 0.3 and we denote this parameter setting by (d). Numerical tests for American puts are implemented in MATLAB (version R 2013a). Table EC.2 compares the classical LSMC method (denoted by LSMC) and the LSMC combined with hybrid sampling (denoted by LSMC-hybrid) while varying the initial stock price S0 when the sample size n = 105 and ds = 5. The true prices are estimated by LSMC with 108 simulation trials using SSR. Table EC.2

American option prices, estimated STDs, and MSEs for the normal tempered stable process

computed using LSMC and LSMC-hybrid when n = 105 , d = 50, and ds = 5 under the parameter set (d). S0

80 LSMC

90

100

LSMC-hybrid

LSMC

LSMC-hybrid LSMC LSMC - hybrid

Price 19.9354

19.8542

11.9069

11.8894

7.3206

7.1996

STD

0.0098

0.0082

0.0343

0.0336

0.0319

0.0312

Bias

0.0323

0.0488

0.0442

0.0616

0.0654

0.0557

MSE

0.0011

0.0024

0.0031

0.0049

0.0053

0.0041

Table EC.3 compares the same option prices, varying κ and σ. Here, we report the relative differences (PLSMC − PLSMC-hybrid )/PLSMC × 100 in percentage, instead of computing true prices, where Pmethod denotes the option price computed by the method. For various κ and σ values, we notice that the percentage differences of the put prices from LSMC and LSMC-hybrid are between 0.1% and 0.9%. The optimal ds and db may differ as parameters change. Under the parameter setting (d), we let d = 256 and do the same experiments by increasing ds from 1 to 32. From this exercise, it is

ec19

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

Table EC.3

American option prices, estimated STDs, and percentage differences for the normal tempered

stable process computed using LSMC and LSMC-hybrid when n = 105 , d = 50, and ds = 5 varying κ and σ under the parameter set (d) σ κ

0.5

0.5

1

1.25

LSMC

LSMC-hybrid

LSMC

LSMC-hybrid

LSMC

LSMC-hybrid

LSMC

LSMC-hybrid

Price

16.4483

16.2975

24.5177

24.6316

33.9261

33.8523

44.6802

44.6236

STD

0.0632

0.0631

0.0816

0.0817

0.0952

0.0949

0.1029

0.1026

RD(%)

1

0.75

0.9168

0.4648

0.2174

0.1266

Price

15.7662

15.6868

23.2012

23.3266

33.1135

33.056

48.0758

48.3182

STD

0.0671

0.0672

0.0835

0.083

0.0948

0.0946

0.0988

0.0985

RD(%)

Table EC.4

0.5031

0.5402

0.1737

0.5043

American option prices, estimated STDs, and MSEs for the normal tempered stable process

computed using LSMC and LSMC-hybrid when n = 105 and d = 256, varying ds under the parameter set (d) S0

80 LSMC

ds

90

LSMC-hybrid 32

16

LSMC 8

100

LSMC-hybrid 32

16

LSMC 8

LSMC-hybrid 32

16

8

Price 19.9532 19.9525 19.9433 19.9181 11.8845 11.8840 11.8618 11.7913 7.2521 7.2540 7.2332 7.2256 STD

0.0041

0.0035

0.0035

0.0029

0.0342

0.0332

0.0339

0.0348

0.0317 0.0309 0.0311 0.0320

Bias

0.0275

0.0281

0.0374

0.0625

0.0881

0.0886

0.1108

0.1813

0.0205 0.0185 0.0394 0.0469

MSE

0.0007

0.0008

0.0014

0.0039

0.0089

0.0089

0.0134

0.0341

0.0014 0.0013 0.0025 0.0032

seen that db value for LSMC-hybrid not greater than 20 could be a reasonable choice, yielding MSEs comparable to LSMC. We report this experiment briefly in Table EC.4. Finally, we want to emphasize that such a gain is very important since it assures the convergence of the numerical scheme in actual implementations. For example, if d = 256, then MATLAB cannot perform LSMC if n > 105 while LSMC-hybrid can easily perform the task with ds = 16.

ec20

e-companion to Kim and Kim: Simulation of tempered stable L´ evy bridges

References Avramidis AN, L’Ecuyer P (2006) Efficient Monte Carlo and Quasi-Monte Carlo option pricing under the variance gamma model. Management Science, 52(12):1930–1944. Hardin Jr C, Samorodnitsky G, Taqqu MS (1991) Nonlinear regression of stable random variables. The Annals of Applied Probability, 1(4):481–634. Kaishev V, Dimitrova D (2009) Dirichlet bridge sampling for the variance gamma process: Pricing pathdependent options. Management Science, 55(3):483–496. Longstaff, F. A., Schwartz, E. S. (2001) Valuing American options by simulation: A simple least-squares approach. Review of Financial Studies, 14(1):113–147. Pellegrino T, Sabino P (2015) Enhancing least squares Monte Carlo with diffusion bridges: An application to energy facilities. Quantitative Finance, 15(5):761–772. Samorodnitsky G, Taqqu MS (1994) Stable Non-Gaussian Random Processes. Chapman & Hall/CRC, Boca Raton, FL. Stentoft L (2004) Convergence of the least squares Monte Carlo approach to American option valuation. Management Science, 50(9):1193–1203. Wood ATA, Booth JG, Butler RW (1993) Saddlepoint approximations to the CDF of some statistics with nonnormal limit distributions. Journal of the American Statistical Association, 88(422):680–686.