An importance sampling-based smoothing approach for quasi-Monte Carlo simulation of discrete barrier options

Accepted Manuscript

An Importance Sampling-Based Smoothing Approach for Quasi-Monte Carlo Simulation of Discrete Barrier Options Fei Xie, Zhijian He, Xiaoqun Wang PII: DOI: Reference:

S0377-2217(18)30883-X https://doi.org/10.1016/j.ejor.2018.10.030 EOR 15422

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

27 April 2017 29 August 2018 12 October 2018

Please cite this article as: Fei Xie, Zhijian He, Xiaoqun Wang, An Importance Sampling-Based Smoothing Approach for Quasi-Monte Carlo Simulation of Discrete Barrier Options, European Journal of Operational Research (2018), doi: https://doi.org/10.1016/j.ejor.2018.10.030

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • An importance sampling method is proposed to smooth payoffs of barrier options.

AC

CE

PT

ED

M

AN US

CR IP T

• A good path generation ordering is proposed to reduce the effective dimension. Numerical studies show the advantage of the proposed method in quasi-Monte Carlo.

1

ACCEPTED MANUSCRIPT

An Importance Sampling-Based Smoothing Approach for Quasi-Monte Carlo Simulation of Discrete Barrier Options Fei Xiea , Zhijian Heb,∗, Xiaoqun Wanga a Department

of Mathematical Sciences, Tsinghua University, Beijing 100084, China of Mathematics, South China University of Technology, Guangzhou 510641, China

CR IP T

b School

Abstract

ED

M

AN US

Handling discontinuities in financial engineering is a challenging task when using quasi-Monte Carlo (QMC) method. This paper develops a so-called sequential importance sampling (SIS) method to remove multiple discontinuity structures sequentially for pricing discrete barrier options. The SIS method is a smoothing approach based on importance sampling, which yields an unbiased estimate with reduced variance. However, removing discontinuities still may not recover the superiority of QMC when the dimensionality of the problem is high. In order to handle the impact of high dimensionality on QMC, one promising strategy is to reduce the effective dimension of the problem. To this end, we develop a good path generation method with the smoothed estimator under the Black-Scholes model and models based on subordinated Brownian motion (e.g., Variance Gamma process). We find that the order of path generation influences the variance of the SIS estimator, and show how to choose optimally the first generation step. As confirmed by numerical experiments, the SIS method combined with a carefully chosen path generation method can significantly reduce the variance with improved rate of convergence. In addition, we show that the effective dimension is greatly reduced by the combined method, explaining the superiority of the proposed procedure from another perspective. The SIS method is also applicable for general models (with the Euler discretization). The smoothing effect of the SIS method facilitates the use of general dimension reduction techniques in reclaiming the efficiency of QMC.

CE

PT

Keywords: simulation, importance sampling, quasi-Monte Carlo, discrete barrier option, smoothing 2010 MSC: 65C05, 65D30, 91G20, 91G60

AC

1. Introduction

Barrier options are popular financial instruments in practice. They can be broadly categorized into two types: “knock-out” options and “knock-in” options (see Derman and Kani, 1996, 1997). Thanks to the higher probability of expiring worthless than the corresponding options without ∗ Corresponding

author. Tel.: +86 13450361878. Email addresses: [email protected] (Fei Xie), [email protected] (Zhijian He), [email protected] (Xiaoqun Wang) Preprint submitted to European Journal of Operational Research

October 17, 2018

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

barriers, barrier options are usually cheaper. This paper concentrates on pricing discretely monitored barrier options in the quasi-Monte Carlo (QMC) framework. We propose an importance sampling (IS) based method to enhance the efficiency of QMC. In respect of pricing barrier options, most models assume continuous monitoring of the barrier, i.e., if the barrier is reached at any instant during the life of the option, a knock-out or knock-in is deemed to occur. The closed-form solutions for continuously monitored barrier option pricing are available (see, e.g., Merton, 1973; Heynen and Kat, 1994a,b). However, many contracts embedded with barrier options in practice are featured by discrete monitoring, i.e., the barrier crossing is monitored only at some prescribed times, such as daily closing (see Birge and Linetsky, 2007). Many asset classes, such as equity, commodity and interest rates, on which the barrier options are usually written, are discretely monitored at a regular frequency (see Lian et al., 2017). These concerns contribute to the preference of discretely monitored barrier options to continuous monitored ones. Unfortunately, the analytical results for continuously monitored barrier options cannot be extended to the discrete ones. Essentially, simple analytical solutions for the prices of discretely monitored barrier options do not yet exist, even under the Black-Scholes (BS) model. In the literature, there have been a range of researches contributing to the techniques for pricing discrete barrier options. Broadie et al. (1997) introduced a continuity correction for approximate pricing of discrete barrier options under the BS model. Dia and Lamberton (2011) and Fuh et al. (2013) extended the continuity correction to jump-diffusion models with the assumption that the correction factor is independent of the jump part, while the assumption was found questionable later (see Lian et al., 2017). There also exist some semi-analytical methods for pricing discrete barrier options, mainly including the fast Hilbert transform-based approach (see Feng and Linetsky, 2008), the Wiener-Hopf method (see Fusai et al., 2006; Green et al., 2010; Fusai et al., 2016) and the Fourier cosine series method (see Fang and Oosterlee, 2009; Lian et al., 2017). As for the numerical methods for pricing discrete barrier options, lattice methods (see Duan et al., 2003; Dai and Lyuu, 2010) and the finite difference methods (see Zhu and De Hoog, 2010; Golbabai et al., 2014) are among the most popular ones. However, these methods may not be useful if stochastic processes used to model the underlying variables are too complex or high dimensional. Simulation-based methods thus become essential tools for pricing barrier options. They have a wider scope and are more portable (i.e., less model dependent). Monte Carlo (MC) simulation with n samples has a root mean squared error O(n−1/2 ) independent of the dimension of the problem. Glasserman and Staum (2001) proposed an MC sampling scheme to advance the path conditional on one-step survival at each monitoring date for pricing barrier options. This method produces a substantial improvement compared to plain MC sampling, especially when the probability of knock-out is high. QMC method, on the other hand, has the advantage to accelerate the MC rate to O(n−1 (log n)d ) for d-dimensional integrands satisfying appropriate conditions (Niederreiter, 1992). QMC method has attracted a great attention in financial engineering (L’Ecuyer, 2009; Glasserman, 2004). Achtsis et al. (2013a,b) developed a conditional QMC sampling scheme compatible with the linear transformation (LT) algorithm (Imai and Tan, 2006) for pricing knock-out barrier options under the BS model and the Heston model. He and Wang (2017) then extended the method proposed by Achtsis et al. (2013b) to general financial options. Avramidis and L’Ecuyer (2006) combined the gamma bridge sampling with QMC to improve the efficiency of pricing path-dependent options under the Variance Gamma (VG) model, in which large efficiency improvement for barrier options was presented. Pricing barrier options is a challenging task in QMC due to two kinds of difficulties. First, the payoffs of barrier options are discontinuous. Direct use of QMC for discontinuous integrands 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

may bring small or no efficiency gain compared to MC (Wang and Tan, 2013), since QMC is designed for smooth functions. Second, the efficiency of QMC may depend highly on the dimension of the problem. He and Wang (2015) found a convergence rate for QMC integration with discontinuous functions, which deteriorates quickly as the dimension d goes up. To handle the discontinuity, we develop a so-called sequential importance sampling (SIS) method to price discrete barrier options. The SIS method can completely remove the discontinuities in the payoff functions, and thus is expected to improve the QMC accuracy. To handle the high dimensionality, we aim at finding a path generation method (PGM) that can reduce the effective dimension of the SIS estimator. Reducing the effective dimension is a promising strategy in improving the QMC efficiency (Imai and Tan, 2006). Taking the knock-out barrier options as examples, we intensively study how to design the optimal PGM to maximize the efficiency under the BS model and the VG model. By combining the SIS method with the optimal PGM, the integrand not only becomes smoother, but also has smaller effective dimension. As a result, the convergence rate of the QMC method increases. There are some related work on handling discontinuous payoffs with high dimensionality in the QMC literature. One strategy is to realign the discontinuities such that they are QMC-friendly (see, e.g., Wang and Tan, 2013; He and Wang, 2014; Imai and Tan, 2014). A common disadvantage of this approach is that the discontinuities are still involved in the resulting integrands despite of the realignments. Another branch of strategies is to use a combination of smoothing methods and dimension reduction methods. To this end, He and Wang (2017) proposed a so-called variable push-out (VPO) method under the assumption that the indicator function involved in the payoff function is variable separable. Weng et al. (2017) then developed a general dimension reduction method to combine with the VPO method. However, the validation of the variables separation condition required for the VPO method depends on both the model and the payoff. For pricing barrier options, although the VPO method delivers a continuous estimator, the estimator still involves kink structures (i.e., max functions). This is why the efficiency gain of VPO is minor for barrier options (He and Wang, 2017). The SIS estimator that we propose can have unlimited smoothness, which is much smoother than the VPO estimator. The SIS method can be easily applied to a variety of asset price models, only if the asset prices can be generated by the inverse transform method step by step. Commonly used models, such as the BS model and the class of asset price models based on subordinated Brownian motion, are all within this range. For general models with the Euler discretization, the method is also applicable. The application range of the proposed method is not restricted to barrier options. It may be extended to other types of options with a similar structure in their payoffs. The remainder of this paper is organized as follows. Section 2 introduces the basic idea of QMC integration. The conditioning on one-step survival method proposed by Glasserman and Staum (2001) is revisited in Section 3. In Section 4, we present a general smoothing method based on IS. In Section 5, we analyze how the variance reduction can be maximized by choosing an optimal first generation step under the BS model and the VG model as an illustrative example of the subordinated Brownian motion. In Section 6, numerical experiments on pricing knockout barrier options are performed under the BS model and the VG model, and the convergence rates of different QMC methods are also compared. In particular, we show how to use the SIS method for a more general model combining with dimension reduction methods in Section 6. Conclusions are given in Section 7. In the supplementary material, we present the numerical results of effective dimension related characteristics for the examples in Section 6, and report the performance of the proposed method on pricing lookback options under the BS model. 4

ACCEPTED MANUSCRIPT

2. Quasi-Monte Carlo Simulation QMC methods are usually used to approximate integrals over the unit cube (0, 1)d Z I( f ) = f (u) du by the average

n

1X f (ui ), Iˆn ( f ) = n i=1

CR IP T

(0,1)d

(1)

M

AN US

where u1 , . . . , un ∈ (0, 1)d are (deterministic) low discrepancy points that are distributed more uniformly than random points used in the MC sampling. We refer to the monographs Niederreiter (1992) and Dick and Pillichshammer (2010) for details on various constructions of low discrepancy points. For QMC, it is hard to estimate the actual error |Iˆn ( f ) − I( f )|. Randomized QMC (RQMC) thus becomes a practical technique in error estimation. RQMC randomizes the points u1 , . . . , un and treat the random version of the quadrature Iˆn ( f ) in (1) as an RQMC estimate. Usually, the randomized points are uniformly distributed over (0, 1)d , and the low discrepancy property of the points is preserved under the randomization (see L’Ecuyer and Lemieux (2005) for a survey of various RQMC methods). This ensures that the RQMC estimate is unbiased, and thus its variance can be used to measure the efficiency of the RQMC method as in the MC sampling. The payoff functions of barrier options are defined over the whole d-dimensional space Rd . A necessary first step in applying QMC methods to a practical integral formulated over Rd is to transform the integral into an integral over the unit cube (0, 1)d . The transformation for barrier options depends on the model that we use.

ED

3. Conditional Sampling Method Revisited

PT

We review the conditional sampling method proposed by Glasserman and Staum (2001). Let S (t) denote the price of the underlying asset at time t. Denote S i = S (ti ) and S = (S 1 , . . . , S d )> , where 0 < t1 < · · · < td = T . Suppose that S 1 , . . . , S d have the Markov property. Consider the barrier call option in Example 1.

CE

Example 1. Up-and-out barrier call option with discounted payoff e−rT max {S d − K, 0}1{max {S i } < Hu }, i

(2)

AC

where 1{·} is an indicator function, r is the risk-free interest rate, Hu > K > 0, and Hu > S 0 , the current asset price.   The option’s price can be expressed as p = E q(S) under the risk-neutral measure with q(S) = g(S)1{max {S i } < Hu }, i

where g(S) = e−rT max{S d − K, 0} is a real-valued function over Rd . Define p j (s) := P(S j < Hu |S j−1 = s), 5

(3)

ACCEPTED MANUSCRIPT

as the survival probability conditional on the one-step prior state. For 1 ≤ i ≤ d, define Li :=

i Y

p j (S j−1 ),

j=1

CR IP T

which can be regarded as the “likelihood” of surviving i steps via this path. Then the estimator based on conditionally sampling the asset path such that all S i < Hu , i = 1, . . . , d, is given by qˆ = Ld g(S).

(4)

AN US

This sampling scheme guarantees nonoccurrence of knocking out, and contributes to variance reduction as no path is wasted. The “likelihood” Ld can be seen as an adjustment for bias. Glasserman and Staum (2001) proved that the estimator (4) is unbiased and has lower MC variance. The use of the conditional sampling method in QMC is straightforward. Despite the smoothness, the estimator (4) may have high effective dimension that could make QMC inefficient. 4. Sequential Importance Sampling as a Smoothing Method We generalize the form (3) as

q(X) = g(X)

d Y

1{ai < Xi < bi },

(5)

M

i=1

such that it is tailored to not only Example 1 but also the following knock-out barrier options, where ai , bi are constants satisfying ai < bi , and X = (X1 , . . . , Xd )> .

ED

Example 2. Down-and-out barrier call option with discounted payoff e−rT max {S d − K, 0}1{min {S i } > Hl },

PT

where Hl < S 0 .

i

(6)

CE

Example 3. Double-barrier call option with discounted payoff e−rT max {S d − K, 0}1{max {S i } < Hu }1{min {S i } > Hl }, i

i

(7)

where Hu > K > 0, and Hl < S 0 < Hu .

AC

Note that we use the notation X in (5) instead of S to stress that X may not be the same as S. In fact, {X1 , . . . , Xd } presents the order to generate the asset prices, i.e., a permutation of {S 1 , . . . , S d }. Also it is possible that ai = −∞ or bi = ∞, in which case the target function (5) corresponds to Examples 1 and 2, respectively. In pricing knock-out barrier options with the form (5), a lot of generated paths will end up with a knock-out if the probability of survival is low. This can make the variance among all paths quite large relative to the price. The situation is even more critical in QMC, since the discontinuities in the integrand could cause a huge loss of efficiency. To accelerate the convergence rate of QMC, a natural strategy is to smooth the target function (5). In this section, we present 6

ACCEPTED MANUSCRIPT

f (x1 , . . . , xd ) =

d Y i=1

fi (xi |x1:(i−1) ).

CR IP T

an IS-based approach for smoothing the function in (5) which, at the same time, provides an unbiased estimate with reduced variance.   To estimate the expectation p = E q(X) by QMC methods, one needs to generate X using low-discrepancy points (see Section 2). To this end, we assume that X can be generated via the multivariate inverse transformation proposed by Rosenblatt (1952), which can be used to sample from a non-uniform distribution provided that we can invert all the required conditional cumulative distribution functions (CDFs). More precisely, let F(x) be the joint distribution function of X. Further, let Fi (xi ) be the marginal distribution of Xi , and for i = 2, . . . , d, let Fi (xi |x1:(i−1) ) and fi (xi |x1:(i−1) ) be the conditional CDF and the conditional probability density function (PDF) of Xi given X1 , . . . , Xi−1 , respectively, where x1:(i−1) = (x1 , . . . , xi−1 )> . The joint PDF of X is given by (8)

AN US

The multivariate inverse transformation x = F −1 (u) = F −1 (u1 , . . . , ud ) is defined by     x1 = F1−1 (u1 ),    xi = F −1 (ui |x1:(i−1) ), i = 2, . . . , d. i

CE

PT

ED

M

If u = (u1 , . . . , ud )> ∼ U((0, 1)d ), then X = F −1 (u) has the joint distribution F(x) (see Rosenblatt, 1952). When the random variables X1 , . . . , Xd constitute a Markov chain, the conditional CDFs Fi−1 (ui |x1:(i−1) ) reduce to Fi−1 (ui |xi−1 ). Assume that the CDF Fi and the inverse CDF Fi−1 can be computed analytically or numerically1 . Motivated by Glasserman and Staum (2001), we are going to remove the indicator functions in (5) via a specific IS. To this end, our basic idea is to employ an IS to twist the probability distribution so that the probability of survival becomes one. Under the IS framework, we call the Q support Ω of the product of the indicator functions as the important region, i.e., Ω = di=1 (ai , bi ). In such a way, all observations of X generated from the IS distribution fall into the important region Ω with probability 1, and the indicator functions disappear. To make the idea more precise, we need to define the IS distribution. Let f˜(x1 , . . . , xd ) denote the joint importance density function with the property described above. Typically, f (·) is required to be absolutely continuous with respect to f˜(·), i.e., f (x) = 0 if f˜(x) = 0. In our context, however, only absolute continuity within the support Ω is required because of the presence of the indicator functions 1{ai < Xi < bi }, 1 ≤ i ≤ d. Let L(x) = f (x)/ f˜(x) be the likelihood ratio between f (x) and f˜(x), then we have Z Z ˜ X)L( ˜ X)], ˜ p= g(x) f (x) dx = g(x)L(x) f˜(x) dx = E[g( Ω

Ω

AC

where X˜ is a random vector that is governed by the IS distribution and thus satisfies ai < X˜ i < bi for all 1 ≤ i ≤ d with probability 1 and E˜ denotes the expectation taken under the IS distribution. In the context of financial engineering, X˜ can be regarded as an asset path that would always result in positive payoff. In our applications, g(·) and the likelihood ratio L(·) are usually smooth ˜ X) ˜ which is unbiased. funcions. Therefore, we arrive at a smoothed estimator g( X)L( 1 If the characteristic function or the PDF of the variable X conditional on X , . . . , X i 1 i−1 is only known explicitly, one may resort to some numerical inversion methods to approximate the conditional CDF Fi and its inverse CDF Fi−1 from the characteristic function or the PDF (see, e.g., Chen et al., 2012; H¨ormann and Leydold, 2003).

7

ACCEPTED MANUSCRIPT

where

f˜i (xi |x1:(i−1) ) =

(9)

CR IP T

For our purpose, we construct an IS distribution with a joint density of the form  d  Y     f˜i (xi |x1:(i−1) ), x ∈ Ω,  f˜(x) =   i=1    0, x < Ω, fi (xi |x1:(i−1) ) . Fi (bi |x1:(i−1) ) − Fi (ai |x1:(i−1) )

Recall that Fi (xi |x1:(i−1) ) is the conditional CDF of Xi given X1 , . . . , Xi−1 . In this case, the likelihood ratio L(x) becomes d Y i=1

[Fi (bi |x1:(i−1) ) − Fi (ai |x1:(i−1) )].

(10)

AN US

L(x) =

In simulation, we may use the inversion method sequentially. In the first step, generate a random sample u1 ∼ U((0, 1)), and let    x˜1 = F1−1 F1 (b1 ) − F1 (a1 ) u1 + F1 (a1 ) .

M

For i = 2, . . . , d, if x˜1 , . . . , x˜i−1 have already been generated and do not make x˜ fall out of Ω, set    x˜i = Fi−1 Fi (bi | x˜ 1:(i−1) ) − Fi (ai | x˜ 1:(i−1) ) ui + Fi (ai | x˜ 1:(i−1) ) x˜ 1:(i−1) ,

iid ˜ X) ˜ where ui ∼ U((0, 1)) and x˜ 1:(i−1) = ( x˜1 , . . . , x˜i−1 )> . By doing so, the smooth estimate g( X)L( d can be viewed as a function of u ∈ (0, 1) , say, H(u). We can estimate p by n

n

ED

1X ˜ 1X H(ui ) = g( Xi )L( X˜ i ). pˆ = Iˆn (H) = n i=1 n i=1

(11)

AC

CE

PT

In the context of QMC, it is also straightforward to simulate X˜ under the IS measure by using low discrepancy points ui ∈ (0, 1)d instead of random points. Using the above IS density, the d indicator functions in the right-hand side of (5) are removed sequentially. Also, it has the effect of smoothing which is crucial for QMC. Particularly, if the functions g(x), Fi (xi |x1:(i−1) ) and fi (xi |x1:(i−1) ) have unlimited smoothness, then H(u) also has unlimited smoothness. From this point of view, the SIS estimator is much smoother than the estimators proposed by Achtsis et al. (2013a,b) and He and Wang (2017). The later ones involve some max functions. That is why we call our method as the SIS smoothing method. The smoothness enables us to use some dimension reduction techniques. We should clarify that in this paper SIS is not the IS technique used in the context of particle filters, although they bear some similarities. In the MC context, the SIS method mainly contributes to variance reduction. The next theorem indicates the unbiasedness and the effectiveness of the proposed estimator.

Theorem 1. Suppose that q(X) is the payoff given in (5) and X has the PDF f (x) given by (8). Let X˜ be a random vector having a joint PDF f˜(x) given by (9). Then ˜ X)L( ˜ X)] ˜ = E q(X) , E[g( (12) 8

ACCEPTED MANUSCRIPT

where E˜ denotes the expectation taken under the IS measure and L(x) = f (x)/ f˜(x) is the likelihood ratio. Moreover, we have g X)L( ˜ X)] ˜ ≤ [F1 (b1 ) − F1 (a1 )]Var q(X) , Var[g(

(13)

CR IP T

g denotes the variance taken under the IS measure. If Xi ’s are independent, then we where Var have d Y   g X)L( ˜ X)] ˜ ≤ [Fi (bi ) − Fi (ai )]Var q(X) . Var[g( i=1

AN US

Proof. For the first part, it is easy to see that f (x) is absolutely continuous with respect to f˜(x) Q within the important region Ω = di=1 (ai , bi ). Hence, the unbiasedness is guaranteed by (9). Now we prove the second part. Let L1 = F1 (b1 ) − F1 (a1 ). We first prove that h i ˜ X) ˜ 2 L( X) ˜ 2 ] ≤ L1 E g(X)2 1{X ∈ Ω} . E[g(

Based on some basic algebra, we can find that

h i ˜ X) ˜ 2 L( X) ˜ 2 ] = E g(X)2 L(X)1{X ∈ Ω} E[g( h i ≤ L1 E g(X)2 1{X ∈ Ω} ,

(14)

where we use the fact from (10) that L(x) ≤ L1 ≤ 1 for all x ∈ Ω. Together with (14), we have

M

˜ X) ˜ X)L( g X)L( ˜ X)] ˜ = E[g( ˜ 2 L( X) ˜ 2 ] − (E[g( ˜ X)]) ˜ 2 Var[g( h i   ≤ L1 E g(X)2 1{X ∈ Ω} − L1 (E g(X)1{X ∈ Ω} )2   = L1 Var g(X)1{X ∈ Ω} . Qd

ED

If Xi ’s are independent, then L(X) = proved similarly.

i=1 [F i (bi )

− Fi (ai )] is a constant. The third part can be

CE

PT

Theorem 1 shows that the MC variance is reduced by at least a factor 1/[F1 (b1 ) − F1 (a1 )], compared to the original one. The factor is defined as the guaranteed efficiency gain (GEG), the least possible efficiency improvement that can be assured by the SIS method in the MC setting, i.e., RGEG := 1/[F1 (b1 ) − F1 (a1 )]. (15)

AC

Clearly, RGEG > 1 confirms that the SIS method indeed reduces the variance. The actual efficiency improvement can be much larger than RGEG . Moreover, RGEG depends on the probability of samples falling out of the important region in the first step under the original probability measure. This insight raises the question of how to optimally generate S if S i can be generated in any order. That is the topic of the next section.

Remark 1. In fact, the proposed SIS method is by no means limited to the application in pricing knock-out barrier options. Consider the down-and-in barrier option whose discounted payoff is given by e−rT max {K − S d , 0}1{min {S i } < Hl }. i

9

ACCEPTED MANUSCRIPT

The critical step of using SIS is to write the payoff function in the form (5). After some rearrangements, the payoff can be decomposed into two parts, i.e., e−rT max {K − S d , 0} − e−rT (K − S d )1{Hl ≤ S d ≤ K}

d−1 Y i=1

1{S i ≥ Hl }.

CR IP T

The first (continuous) part is actually the payoff of a European put option that may have an analytical solution or can be handled easily using QMC. The second (discontinuous) part has exactly the form (5), and thus can be handled by the SIS method. Similarly, the SIS method can also be applied to pricing lookback options; see Appendix B in the supplementary material for details.

AN US

Remark 2. Genz and Bretz (2009) introduced the separation-of-variables method for numerical integration, which employed the similar idea, though focused on a different problem. They dealt with the numerical computation of multivariate normal and t distributions, while in this paper, we cope with derivative pricing problems, whose distributions are not confined to the above two distributions. Moreover, the initiative here is to smooth the integrand to achieve better efficiency in QMC integration. 5. Maximizing Variance Reduction by Choosing a Good PGM

5.1. Black-Scholes Model

M

We are interested in the problem of how to maximize the GEG of the SIS method for the barrier options in Examples 1-3 under the BS model and the VG model by choosing a proper path generation order.

ED

We assume that under the risk neutral measure, the underlying asset price follows a geometric Brownian motion dS (t) = rS (t) dt + σBS S (t) dB(t), (16)

PT

where σBS is the volatility and B(t) is the standard Brownian motion. The explicit solution of (16) is S (t) = S (0) exp{ωt + σBS B(t)}, where ω = r − σ2BS /2. Therefore, we have S i := S (tt ) = S (0) exp{ωt + σBS B(ti )}, i = 1, . . . , d,

(17)

AC

CE

where 0 < t1 < · · · < td = T , and T is the expiration time. There are many constructions of the standard Brownian motion, among which the simplest one is the standard (ST) construction (i.e., the random walk construction). The ST construction generates the standard Brownian motion sequentially: given B(0) = 0, √ (18) B(ti ) = B(ti−1 ) + ti − ti−1 zi , zi ∼ N(0, 1), i = 1, . . . , d.

In fact, we may generate the Brownian motion in any order, provided that at each time step we sample from the correct conditional distribution given the values already generated. Brownian bridge (BB) construction is an alternative method widely used in practice (Caflisch et al., 1997; Moskowitz and Caflisch, 1996). Under the BB construction, one generates the final value B(td ) firstly, then generates the remaining values conditional on the previously generated values. 10

ACCEPTED MANUSCRIPT

Specifically, suppose 0 < u < s < t, the BB construction generates B(s) conditioning on B(u) and B(t) based on the following formula r (t − s)(s − u) (t − s)B(u) + (s − u)B(t) B(s) = + z, z ∼ N(0, 1). (19) t−u t−u

AN US

CR IP T

In this paper, we default that the BB construction generates the Brownian motion in the order: B(td ), B(tbd/2c ), B(tbd/4c ), B(tb3d/4c ) and so on. It means that after generating the final value B(td ), we generate B(tbd/2c ) conditional on the values of B(0) and B(td ), and proceed by progressively filling in the intermediate values. In the simulation framework, it suffices to generate the Brownian motion B(ti ), which can be performed in any order as mentioned. According to (17), the asset prices S i are then generated in the same order of generating the Brownian motion. Generally, assume that the Brownian motion and the asset prices are generated in the time order {τ1 , . . . , τd }, which is a permutation of {t1 , . . . , td }. Let Xi = S (τi ) for i = 1, . . . , d and τ0 = 0. Taking Example 3 as an illustrative case, the function (7) can be rewritten as the form (5), i.e., q(X) = g(X)

d Y i=1

1{Hl < Xi < Hu },

(20)

ED

M

where g(X) = e−rT max {Xi∗ − K, 0} with Xi∗ = S d for some index i∗ (depending on the specific generation order). To apply the SIS method, it suffices to know the conditional CDF Fi (·|X1:(i−1) ) and its inversion Fi−1 . Note that   Fi (H|X1:(i−1) ) = P Xi ≤ H X1:(i−1)   = P B(τi ) ≤ γ(τi , H) B(τ1 ), . . . , B(τi−1 ) ,

CE

PT

where γ(τi , H) = [log(H/S 0 )−ωτi ]/σBS . If τi > Mi := max{τ0 , . . . , τi−1 }, we generate B(τi ) using the sequential construction based on B(Mi ). By replacing B(τi ) in the above formula according to (18), we get   p Fi (H|X1:(i−1) ) = P B(Mi ) + τi − Mi z ≤ γ(τi , H)  γ(τ , H) − B(M )  i i = Φ , √ τ i − Mi

AC

where Φ(·) is the CDF of the standard normal distribution. Otherwise, if there exist two time points τ j and τk in the set {τ0 , τ1 , . . . , τi−1 } satisfying τ j < τi < τk , then by applying (19), B(τi ) can be generated by the BB construction, and we have   Fi (H|X1:(i−1) ) = P α(τi , τ j , τk ) + β(τi , τ j , τk )z ≤ γ(τi , H) B(τ j ), B(τk )    = Φ γ(τi , H) − α(τi , τ j , τk ) β(τi , τ j , τk ) , where

α(τi , τ j , τk ) =

(τk − τi )B(τ j ) + (τi − τ j )B(τk ) , τk − τ j 11

ACCEPTED MANUSCRIPT

and β(τi , τ j , τk ) =

s

(τk − τi )(τi − τ j ) . τk − τ j

CR IP T

It is straightforward to compute the inverse CDF Fi−1 . Thus the SIS estimator can be obtained as constructed in Section 4. For Example 1, the SIS estimator can be obtained by taking Hl = −∞. Similarly, taking Hu = ∞ gives the estimator for Example 2. Remark 3. If the asset prices are generated sequentially by setting τi = ti for all i = 1, . . . , d (i.e., in every step the Brownian motion is generated using the ST construction), then our method reduces to the conditioning on one-step survival method proposed by Glasserman and Staum (2001), revisited in Section 3. We include the conditioning on one-step survival method for comparison in our numerical experiments. Lindset and Persson (2016) used non-equal grid sizes in generating the asset prices. Our method could be fitted into this setting.

AN US

We have showed in Theorem 1 that the SIS method brings a variance reduction regardless of the order in which Xi ’s are generated. However, the order of generating the asset prices may have an important impact on the performance of the SIS method. The next theorem shows that RGEG defined in (15) depends heavily on the time τ1 at which the asset price is firstly generated.

Theorem 2. Assume that we generate the Brownian motion and the asset prices in the order of times τ1 , . . . , τd . Suppose that the function q(X) is given by (20). Then the associated GEG of the SIS method is given by

M

! !#−1 " γ(τ1 , Hl ) γ(τ1 , Hu ) −Φ . RGEG (τ1 ) = Φ √ √ τ1 τ1

(21)

 √  F1 (b1 ) = P(S (τ1 ) ≤ Hu ) = P(B(τ1 ) ≤ γ(τ1 , Hu )) = Φ γ(τ1 , Hu )/ τ1 ,

PT

and

ED

Proof. To translate into the form (5), we have ai = Hl and bi = Hu . According to (17),  √  F1 (a1 ) = P(S (τ1 ) ≤ Hl ) = P(B(τ1 ) ≤ γ(τ1 , Hl )) = Φ γ(τ1 , Hl )/ τ1 , where we use B(τ1 ) ∼ N(0, τ1 ). It follows from (15) that

CE

RGEG

" ! !#−1 1 γ(τ1 , Hu ) γ(τ1 , Hl ) = = Φ −Φ , √ √ F1 (b1 ) − F1 (a1 ) τ1 τ1

which explicitly depends on τ1 , the time at which the asset price is firstly generated.

AC

Based on Theorem 2, one may choose the time τ1 to maximize RGEG for the purpose of improving efficiency. Alternatively, the payoff function (7) can be rewritten in another form Y q(X) = g˜ (X)1{max {K, Hl } < Xi∗ < Hu } 1{Hl < Xi < Hu }, (22) i,i∗

where g˜ (X) = e−rT (Xi∗ − K) and again Xi∗ = S d . Note that g˜ (X) is smoother than g(X) given in (20). Similar to Theorem 2, we have the following result for the function (22). 12

ACCEPTED MANUSCRIPT

CR IP T

Theorem 3. Assume that the Brownian motion and the asset prices are generated in the order of times τ1 , . . . , τd . Suppose that the function q(X) is given by (22). Then the associated GEG of the SIS method is given by     −1     , τ1 < T, Φ γ(τ√1 τ,H1 u ) − Φ γ(τ√1τ,H1 l )     −1   (23) R˜ GEG (τ1 ) =    γ(T,max {K,Hl }) γ(T,H )   √ , τ = T.  Φ √T u − Φ 1 T Corollary 1. Suppose that RGEG is given by (21) and R˜ GEG is given by (23), then we have R˜ GEG (τ1 ) ≥ RGEG (τ1 ).

For an arbitrary τ1 ∈ (0, T ], we have limHu −Hl →0 RGEG (τ1 ) = limHu −Hl →0 R˜ GEG (τ1 ) = ∞.

AN US

Corollary 1 shows the advantage of taking the (implied) indicator function 1{S d > K} into consideration, besides the indicator functions for the barrier. As a result, we propose that the first step τ1 is determined by maximizing the function R˜ GEG in numerical experiments. Generally, one could simply calculate the function values R˜ GEG (ti ) for i = 1, . . . , d and let τ1 = tk if R˜ GEG (tk ) ≥ R˜ GEG (ti ) for i , k. The computation time for this is negligible. Furthermore, the level of GEG depends on the gap between Hu and Hl and between Hu and max {K, Hl }. The closer Hu and Hl are, the greater GEG is, and thus the greater efficiency improvement could be expected. 5.2. Models Based on Subordinated Brownian Motion

ED

M

In this subsection we extend the results to the asset price model based on subordinated Brownian motion, which is an important class of models in derivatives pricing. A subordinated Brownian motion is defined by W(t) = B(L(t); θ, σ),

CE

PT

where the subordinator L(t) is an increasing L´evy process with independent and stationary increments, and B(t) := B(t; θ, σ) is a Brownian motion with drift parameter θ and variance parameter σ. Both processes have initial zero value L(0) = 0 and B(0) = 0. Examples of subordinated Brownian motion are the VG process and the normal inverse Gaussian (NIG) process (see Cont and Tankov, 2004). Assume the asset price process can be explicitly expressed as S (t) = S (0) exp{ωt + W(t)},

(24)

AC

where ω is a constant. Note that the Brownian motion B(t) is independent of the incremental process L(t). For t ≥ 0 and δ > 0, the Brownian motion’s increments follow normal distribution, i.e., B(t + δ) − B(t) ∼ N(θδ, σ2 δ).

Conditional on the increment L(t + δ) − L(t), the increment B(L(t + δ)) − B(L(t)) has a normal distribution with mean θ[L(t + δ) − L(t)] and variance σ2 [L(t + δ) − L(t)]. In simulation, we could generate the entire random time path L(t1 ), . . . , L(td ) firstly, and then simulate the asset prices at these times. To apply the SIS smoothing method, we can regard the entire random time path L(t1 ), . . . , L(td ) as known quantities, and calculate the conditional CDF Fi (xi |x1:(i−1) ) given X1 , . . . , Xi−1 . 13

ACCEPTED MANUSCRIPT

Similar to the case of the BS model, we could define the conditional GEG as the GEG of the SIS method conditional on the knowledge of the subordinator process values. By calculating the conditional probability of X1 falling between the corresponding a1 and b1 , we obtain the following theorem.

CR IP T

Theorem 4. Assume that the Brownian motion and the asset prices are generated in the order of times τ1 , . . . , τd . Suppose that the function q(X) is given by (20). Then for the asset price model based on the subordinated Brownian motion, the conditional GEG of the SIS method is given by ! !#−1 " log(Hl /S 0 ) − ωτ1 − θL(τ1 ) log(Hu /S 0 ) − ωτ1 − θL(τ1 ) −Φ . RCGEG (τ1 ) = Φ √ √ σ L(τ1 ) σ L(τ1 )

(25)

PT

ED

M

AN US

Note that RCGEG (τ1 ) in (25) depends on the random variate L(τ1 ). In practice, we substitute L(τ1 ) in (25) with its expectation value E [L(τ1 )] to estimate the expected conditional GEG R˜ CGEG (τ1 ). In this way, we circumvent the random factor within the expression of RCGEG . Now, one can determine the first generation step τ1 by maximizing R˜ CGEG (τ1 ). Meanwhile, we can obtain the similar results as Theorem 3 and Corollary 1. In other words, there is further improvement in conditional GEG if we write the payoff function in the form (22). As an illustrative example, we consider the VG model. The VG model is a Brownian motion with a random time change that follows a stationary gamma process. Specifically, the subordinator L(t) is a gamma process G(t) := G(t; µ, ν) with drift µ = 1 and volatility ν > 0 (see Madan et al., 1998). One of the advantages of this model is that it permits more flexibility in modeling skewness and kurtosis relative to Brownian motion. Empirical evidence has been provided that this model gives a better fit to market option prices than the BS model (Avramidis and L’Ecuyer, 2006). Although there exists efficient techniques such as the fast Fourier transformation for option pricing under the VG model (Carr and Madan, 1999), it is not sufficiently well developed for path-dependent options such as the barrier option. In this context, MC and QMC methods can be simply implemented, but further efficiency improvement is needed. Under the risk-neutral measure, the asset price based on the VG process W(t) is of the form (24) with ω = log(1 − θν − σ2VG ν/2)/ν + r − c, where σ = σVG , r is the risk-free interest rate, and c is the asset’s continuously compounded dividend yield (see Madan et al., 1998). Based on Theorem 4, by substituting the subordinator value at t = τ1 in (25) with the expectation value E [G(τ1 )] = τ1 , we obtain the expected conditional GEG under the VG model (26)

CE

" ! !#−1 log(Hu /S 0 ) − (ω + θ)τ1 log(Hl /S 0 ) − (ω + θ)τ1 R˜ CGEG (τ1 ) = Φ −Φ . √ √ σVG τ1 σVG τ1

AC

The function R˜ CGEG (τ1 ) looks very like the GEG function under the BS model (21). When applying the SIS method, one can optimally choose the first generation step τ1 which gives the maximal R˜ CGEG (τ1 ) value. We present the numerical results of pricing various types of barrier options under the VG model in Section 6.2. Remark 4. To fit into QMC framework, one needs to generate the subordinator process at the discrete times L(t1 ), . . . , L(td ) via low discrepancy points. To keep the smoothing effect, it is desirable to choose the associated mapping function from the unit cube to the asset prices be smooth. This is the case for the Gamma processes (subordinator of the VG process) and the inverse Gaussian processes (subordinator of the NIG process). The increments of both processes follow a distribution with continuous density (see Cont and Tankov, 2004), which can be sampled 14

ACCEPTED MANUSCRIPT

by the inverse transform method. This method maps the input u ∈ (0, 1) monotonically and continuously to the desired random variable (see Glasserman, 2004). However, the requirement of smooth mapping is not generally guaranteed for other subordinators. 6. Numerical Experiments

CR IP T

This section carries out numerical experiments to examine the effectiveness of the proposed SIS smoothing method on the BS model and the VG model, and to compare different PGMs empirically. We use a variety of QMC simulation schemes to estimate prices of barrier options in Examples 1-3: • Up-and-out barrier call option: q(S) = e−rT max {S d − K, 0}1{maxi {S i } < Hu };

• Down-and-out barrier call option: q(S) = e−rT max {S d − K, 0}1{mini {S i } > Hl };

AN US

• Double-barrier call option: q(S) = e−rT max {S d − K, 0}1{maxi {S i } < Hu }1{mini {S i } > Hl }.

The implementation procedure for the BS model is outlined in Algorithm 1. The implementation for the VG model is similar, but with a slight modification that the SIS sampling is performed after one process is generated. We omit the details for sake of saving space. In our experiments, we use a randomized version of Sobol points proposed by Matouˇsek (1998). We choose the sample size n = 214 and the number of independent replications R = 100 required in Algorithm 1.

CE

PT

ED

M

Algorithm 1 Pricing barrier options using the SIS method with QMC under the BS model. ˜ GEG (τ1 ). 1: Determine the first step τ1 by maximizing the GEG R 2: Choose a subsequent generation order τ2 , . . . , τd , such as ST, BB, and RV. 3: Choose the number of independent replications R and the size of samples for each replication n. Let X = (S τ1 , . . . , S τd ), and let u1 , . . . , un be the deterministic QMC points (e.g, Sobol points). 4: for i = 1, . . . , R do (i) 5: Generate the RQMC points u(i) 1 , . . . , un by randomizing the deterministic QMC points. 1 Pn 6: Compute the SIS estimate pˆ i = n j=1 H(u(i) j ) with H given in (11). 7: end for 8: Deliver the total computational time tPGM over the R replications, the average of the R estiP PR 1 mates as the final estimate p¯ = R1 Ri=1 pˆ i , and its empirical variance σ ˆ 2PGM = R(R−1) ˆi − i=1 ( p   2 p) ¯ , which is an unbiased estimate of Var p¯ .

AC

We take the efficiency of the crude QMC method (which generates asset price paths with the ST construction) as the benchmark, and use the variance reduction factor (VRF) and the efficiency factor (EF) to demonstrate the performance of each QMC simulation scheme of interest compared to the benchmark. Specifically, let σ ˆ 2PGM be the empirical variance of using some particular QMC scheme of interest and tPGM the computational time as in Step 8 of Algorithm 1, and similarly define σ ˆ 2crude and tcrude for the crude QMC. The VRF and the EF are then defined as VRF := 15

σ ˆ 2crude σ ˆ 2PGM

ACCEPTED MANUSCRIPT

and EF :=

σ ˆ 2crude tcrude σ ˆ 2PGM tPGM

= VRF ×

tcrude , tPGM

(27)

CR IP T

respectively. The EF takes the computational time into account. The larger is the EF, the more effective is the QMC scheme of interest. The standard deviation (SD) σ ˆ crude and the computational time tcrude for the crude QMC method are reported in tables for comparison. All experiments are conducted using MATLAB on Intel Xeon E5-2650 v4 2.20GHz CPU with 12 cores and 72GB RAM. Furthermore, to demonstrate how various PGMs behave with the increasing number of samples, we present the convergence graphs for the crude MC and all the QMC-based simulation schemes. To understand the source of huge diversity in performance under different schemes in the QMC framework, we intentionally compute some effective dimension related characteristics, which are deferred to the supplementary material. 6.1. Examples under Black-Scholes Model

AN US

Recall that the discretized underlying process is given by S i = S (0) exp{ωt + σBS B(ti )}, where ω = r − σ2BS /2. Note that there is no discretization error for the BS model. We compare performance of six different QMC PGMs including the crude QMC method in pricing knock-out options under the BS model: • GS: the conditioning sampling method proposed by Glasserman and Staum (2001); • crude SIS: the SIS method with the sequential order τi = ti ;

M

• SIS+ST: the SIS method with τ1 optimally chosen by maximizing the GEG R˜ GEG (τ1 ) and the following d − 1 steps in the sequential order; • SIS+RV: similar to SIS+ST, but the last d − 1 steps are in the reverse order;

ED

• SIS+BB: similar to SIS+ST, but the last d − 1 steps are in the BB order.

AC

CE

PT

All the SIS schemes take the implied indicator function 1{S d > K} into consideration as illustrated in (22). The crude SIS does not take the GEG optimization procedure into consideration. To investigate the influence of the dimension d, the volatility level σBS and the gap between the barrier and the strike price K, we intentionally design six parameter settings outlined in Table 1 with S 0 = 100, r = 0.01 and T = 0.25 fixed. Note that the lower barrier Hl is not required in Example 1, while the higher barrier Hu is not required in Example 2. Tables 24 report the numerical results on various QMC methods for Examples 1-3, respectively. We only present the computational time for the benchmark scheme – the crude QMC. For other schemes, the computational time can be recovered by the formula (27). We observe that for all the SIS methods and the GS method, VRF is very close to EF. This implies that the computational time for these SIS schemes is close to that of the crude QMC. In other words, using the SIS method does not bring (significant) extra computation burden. Based on the numerical results, we conclude the following. • There is consistent significant advantage of the SIS+BB scheme over other PGMs for the six cases of parameters. Additionally, the improvement achieved by the SIS+ST and SIS+RV methods is also remarkable. The efficiency of the proposed SIS method with optimal first generation step is therefore justified. Under the BB construction, the smoothed 16

ACCEPTED MANUSCRIPT

integrand has the smallest effective dimension (see Appendix A in the supplementary material for details). As a result, SIS+BB has not only the smoothing effect, but also the dimension reduction effect, both beneficial for QMC.

CR IP T

• Compared to the crude SIS scheme, the GS scheme proposed by Glasserman and Staum (2001) targets the payoff of the form (20), which does not take the non-smooth component max {S d − K, 0} into consideration. Corollary 1 shows that the crude SIS scheme yields a GEG no less than that of the GS scheme. That is why the crude SIS scheme performs better than the GS scheme. However, the advantage is minor because the crude SIS scheme uses the sequential order, without optimizing the first generation step.

AN US

• For the example of pricing double-barrier option, large VRFs and EFs are observed in Table 4 for Case (1). We should note that for this case, the probability of knock-out is high, making the crude sampling scheme less effective. This can be verified from the standard deviation, which is much larger than the estimated price for the crude QMC. From this point of view, the proposed SIS method has the power of handling extreme cases (rare events).

ED

M

Figure 1 presents the convergence graphs for the crude MC and all the QMC schemes listed above. We find that the SIS+BB scheme results in lowest standard deviation for all the sample sizes n = 2i , i = 0, . . . , 14. Additionally, the standard deviation of the SIS+BB scheme decays fastest as the sample size increases. The standard deviations of the crude QMC and the crude MC are very close, implying that a naive use of QMC cannot reduce the MC variance significantly. Both the SIS+RV and SIS+ST schemes perform better than the crude SIS and GS schemes, implying the advantage of choosing the first step τ1 optimally. All SIS methods perform much better than the crude sampling schemes, implying the good effect of smoothing. Reducing the effective dimension can further bring orders of magnitude reduction in integration error, as observed in SIS+BB. Table 1: The parameter settings for the BS model

d 64 64 64

σBS 0.5 0.3 0.5

K 95 100 105

Hl 94 92 90

Hu 106 108 110

Case (4) (5) (6)

d 128 128 128

CE

PT

Case (1) (2) (3)

σBS 0.3 0.5 0.3

K 95 100 105

Hl 94 92 90

Hu 106 108 110

6.2. Examples under Variance Gamma Model

AC

In this subsection, we present the numerical results of the VG model. Besides the case we described in Section 5.2 which uses the discretized underlying process S i = S (0) exp{ωti + B(G(ti ), θ, σVG )} (we refer it as Brownian gamma (BG) representation of the VG model), we also include another representation of the VG process for comparison. An alternative representation of the VG process is the difference between two independent gamma processes (see Madan et al., 1998) W(t) = Γ+ (t) − Γ− (t), 17

ACCEPTED MANUSCRIPT

Table 2: Numerical results for up-and-out barrier option under the BS model

Case (1) (2) (3) (4) (5) (6)

crude SD 5.44e-3 6.51e-3 1.85e-3 1.07e-2 3.39e-3 2.86e-3

Time 8 8 9 18 19 18

Price 0.1153 0.1830 0.0201 0.3374 0.0485 0.0496

SIS+ST VRF 56 75 167 25 39 48

EF 50 65 149 23 38 45

GS VRF 2.5 2.4 3.0 2.0 2.1 1.3

EF 2.0 1.9 2.4 1.7 1.8 1.1

SIS+RV Price VRF 0.1152 82 0.1830 131 0.0201 452 0.3369 30 0.0484 80 0.0496 141

EF 73 115 404 27 77 132

Price 0.1153 0.1832 0.0203 0.3376 0.0488 0.0496

crude SIS Price VRF 0.1152 3.1 0.1832 2.7 0.0202 6.3 0.3376 2.0 0.0487 2.5 0.0495 1.5

EF 2.5 2.2 5.0 1.7 2.2 1.3

SIS+BB VRF 1785 2374 6393 455 1520 1551

EF 1600 2129 5805 425 1485 1405

CR IP T

(1) (2) (3) (4) (5) (6)

Price 0.1150 0.1827 0.0200 0.3380 0.0491 0.0496

Price 0.1153 0.1830 0.0201 0.3370 0.0484 0.0496

AN US

Case

M

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes GS (Glasserman and Staum (2001)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.01, T = 0.25, n = 214 and R = 100, and the other parameters for the six cases are given in Table 1. The computational time is reported in seconds.

Table 3: Numerical results for down-and-out barrier option under the BS model

CE

Case

AC

(1) (2) (3) (4) (5) (6)

crude SD 6.88e-2 2.47e-2 4.60e-2 3.76e-2 5.94e-2 3.50e-2

Time 8 8 8 18 18 18

SIS+ST VRF 42 29 54 9 15 130

EF 37 25 46 8 14 114

ED

(1) (2) (3) (4) (5) (6)

Price 7.1763 5.3416 6.4914 6.2621 6.8750 3.8334

PT

Case

Price 7.1827 5.3414 6.4951 6.2666 6.8812 3.8352

GS VRF 3.7 1.6 1.7 1.5 1.6 1.1

EF 3.0 1.3 1.4 1.2 1.3 0.9

crude SIS Price VRF 7.1824 3.7 5.3400 1.6 6.4914 1.7 6.2640 1.5 6.8772 1.6 3.8334 1.1

EF 3.0 1.3 1.4 1.2 1.3 0.9

SIS+RV Price VRF 7.1805 13 5.3412 17 6.4940 25 6.2655 6 6.8800 11 3.8351 74

EF 12 14 22 5 10 65

SIS+BB Price VRF 7.1831 258 5.3414 200 6.4956 307 6.2653 52 6.8804 86 3.8350 581

EF 225 170 263 45 77 512

Price 7.1824 5.3401 6.4911 6.2640 6.8772 3.8335

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes GS (Glasserman and Staum (2001)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.01, T = 0.25, n = 214 and R = 100, and the other parameters for the six cases are given in Table 1. The computational time is reported in seconds.

18

ACCEPTED MANUSCRIPT

Table 4: Numerical results for double-barrier option under the BS model

Case (1) (2) (3) (4) (5) (6)

crude SD 8.67e-5 4.23e-3 5.54e-4 2.38e-3 2.73e-4 2.17e-3

Time 8 8 8 18 18 18

Price 2.10e-5 6.50e-2 1.74e-3 1.44e-2 2.49e-4 3.17e-2

GS VRF 34808.5 6.5 39.7 14.3 135.1 2.3

EF 22019.6 3.8 24.0 8.7 83.7 1.4

Price 2.10e-5 6.50e-2 1.74e-3 1.44e-2 2.49e-4 3.16e-2

crude SIS VRF EF 36273.9 22932.9 7.0 4.2 61.6 36.9 14.3 8.3 153.3 95.3 2.8 1.8

Price 2.10e-5 6.50e-2 1.74e-3 1.45e-2 2.50e-4 3.16e-2

SIS+ST VRF 72674 65 483 36 646 27

EF 53422 45 337 25 463 19

Price 2.10e-5 6.49e-2 1.73e-3 1.44e-2 2.51e-4 3.16e-2

SIS+RV VRF 75810 67 614 27 802 56

EF 55913 46 426 19 575 41

Price 2.10e-5 6.49e-2 1.73e-3 1.45e-2 2.52e-4 3.16e-2

SIS+BB VRF 586917 1679 6976 698 9252 791

CR IP T

(1) (2) (3) (4) (5) (6)

Price 1.83e-5 6.52e-2 1.70e-3 1.44e-2 2.74e-4 3.15e-2

AN US

Case

ED

M

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes GS (Glasserman and Staum (2001)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.01, T = 0.25, n = 214 and R = 100, and the other parameters for the six cases are given in Table 1. The computational time is reported in seconds.

Double-barrier option with parameter (2)

10 0

CE

10 0

PT

Figure 1: The convergence rate for the double-barrier option under the BS model with parameters settings (2) and (6) given in Table 1. The sample sizes are n = 2i , i = 0, . . . , 14. The SDs are estimated based on R = 100 independent replications.

Double-barrier option with parameter (6)

10 -1

Standard Deviation (SD)

10 -2

AC

Standard Deviation (SD)

10 -1

10 -3

10 -4 10 0

crude MC crude QMC GS crude SIS SIS+ST SIS+RV SIS+BB

10 1

10 2

10 3

10 4

10 5

Number of Samples

10 -2

10 -3

10

-4

10 -5 10 0

crude MC crude QMC GS crude SIS SIS+ST SIS+RV SIS+BB

10 1

10 2

10 3

Number of Samples

19

10 4

10 5

EF 437355 1167 4900 501 6725 580

ACCEPTED MANUSCRIPT

CR IP T

where Γ+ (t) and Γ− (t) are qindependent gamma processes with parameters (µ1 , ν1 ) and (µ2 , ν2 ), respectively, with µ1 = ( θ2 + 2σ2VG /ν + θ)/2, µ2 = µ1 − θ, ν1 = µ21 ν and ν2 = µ22 ν. We refer to this representation of the VG model as the difference of gammas (DG). To apply the SIS method, we generate one of the gamma process firstly, and then we use the SIS smoothing method for the other gamma process conditional on the previous generated one. The conditional distribution under the SIS smoothing method is again calculated conditional on all the previously generated asset prices and the entire path of the firstly generated gamma process. The conditional GEG can be derived and maximized to determine the optimal first generation step τ1 . Generally speaking, one could generate either gamma process firstly. However, for the upand-out barrier option, if the gamma process {Γ− (t)} is generated firstly, the likelihood ratio L(x) would be calculated as d   Y L(x) = P Γ+ (ti ) < Γ− (ti ) + log(Hu /S 0 ) − ωti Γ− (ti ) , i=1

−

AN US

which can be zero as Γ (ti ) + log(H/S 0 ) − ωti can be negative and the values of a gamma process must be positive. Zero likelihood ratio would result in zero output, and would reduce the efficiency in either MC setting or QMC setting, as samples are wasted. On the other hand, if the gamma process {Γ+ (t)} is generated firstly, the likelihood ratio L(x) would be calculated as d   Y P Γ− (ti ) > Γ+ (ti ) − log(Hu /S 0 ) + ωti Γ+ (ti ) , L(x) = i=1

CE

PT

ED

M

which would always be positive. Based on this consideration, we generate {Γ+ (t)} firstly, and use the SIS smoothing method for the gamma process {Γ− (t)} with the knowledge of {Γ+ (t)} for the up-and-out barrier option. For the same reason, one may generate the gamma process {Γ− (t)} firstly for the down-and-out barrier option. However, when pricing the knock-out option with double barriers with the SIS method, the problem of zero likelihood ratio cannot be completely avoided whichever gamma process is generated firstly. In our numerical experiments, we generate {Γ+ (t)} firstly for up-and-out and double-barrier options, while we generate {Γ− (t)} firstly for the down-and-out option. We next present numerical results of the VG model with the BG and DG representations. Note that there is no discretization error for either representation. Again in application of the SIS method, we optimize the first generation step, and consider three alternative generation orders for the subsequent steps. As the efficiency of gamma bridge sampling has been shown in Avramidis and L’Ecuyer (2006), we combine it with the SIS method by implementing it on the firstly generated gamma process. We also take their original sampling scheme into comparison. The QMC schemes are outlined as follows.

AC

• AL: the sampling method proposed by Avramidis and L’Ecuyer (2006), where both processes are generated using the bridge method; • crude SIS: the firstly generated gamma process is generated using the gamma bridge, and the other process (Bronwian motion under the BG representation, gamma process under the DG representation) is sequentially generated with SIS; • SIS+ST: the firstly generated gamma process is generated using the gamma bridge, and the other process is generated using SIS with the first generation step τ1 optimally chosen by maximizing the CGEG R˜ CGEG (τ1 ) and the following d − 1 steps in the sequential order; 20

ACCEPTED MANUSCRIPT

• SIS+RV: similar to SIS+ST, but the last d − 1 steps are in the reverse order; • SIS+BB: similar to SIS+ST, but the last d − 1 steps are in the BB order.

CR IP T

The parameter settings are outlined in Table 5 with S 0 = 100, r = 0.0548, c = 0, θ = −0.2859, σVG = 0.1927, ν = 0.2505 and T = 0.40504 fixed. Tables 6–8 present the numerical results on the QMC schemes for Examples 1–3, respectively. We note that the nominal dimension is 2d since we need 2d random variables to generate the two processes. We observe the following.

• Overall speaking, the effectiveness of the SIS method under the VG model is not as striking as under the BS model, and the rank between the two generation orders (BB and RV) varies with different examples and parameter settings. Nevertheless, the SIS method with optimal first generation step always exhibits superiority over the AL scheme proposed by Avramidis and L’Ecuyer (2006) and the crude SIS scheme without optimizing the first generation step.

AN US

• For the up-and-out barrier option, the SIS+RV scheme always performs best with both representations. In pricing the down-and-out barrier option, the SIS+ST scheme performs best with the BG representation, while the three generation orders lead to a similar performance with the DG representation. As for the option with double barriers, the SIS+RV and SIS+BB schemes demonstrate consistent substantial advantages over the SIS+ST scheme.

M

• For the up-and-out barrier option and the double-barrier option, the SIS+RV and SIS+BB schemes seem to work better with the DG representation than with the BG representation. While for the down-and-out barrier option, the SIS methods with optimal first generation order always work better on the BG representation.

ED

• When using the SIS method without carefully choosing the first generation step, the crude SIS scheme results in very limited improvement over the crude QMC, and hardly outperforms the AL scheme proposed by Avramidis and L’Ecuyer (2006). This emphasizes the significance of optimizing the first generation step.

CE

PT

• It seems that the number of time steps d has little effect on the best SIS method. One possible reason is that the effective dimension of the smoothed integrand is reduced greatly, as confirmed in the BS model.

d 16

K 95

Hl 90

Hu 105

Case (2)

d 64

K 100

Hl 96

Hu 104

Case (3)

d 128

K 105

Hl 95

Hu 110

AC

Case (1)

Table 5: The parameter settings for the VG model

6.3. Extension to general models

As an extension, we demonstrate how to apply the SIS method to the general asset pricing model. To this end, we consider a scalar stochastic differential equation (SDE) with general drift a(S (t), t) and volatility b(S (t), t), i.e., dS (t) = a(S (t), t) dt + b(S (t), t) dB(t), 0 ≤ t ≤ T, 21

(28)

ACCEPTED MANUSCRIPT

Table 6: Numerical results for up-and-out barrier option under the VG model

(2) (3)

Case (1) (2) (3)

BG DG BG DG BG DG

crude SD 1.34e-2 1.12e-2 6.56e-3 5.80e-3 3.05e-3 2.89e-3

Time 34 60 97 167 153 258

Price 0.6145 0.6155 0.2549 0.2550 0.0660 0.0663

AL VRF 4.4 3.1 4.3 4.1 2.2 3.3

EF 1.9 1.3 0.8 0.7 0.3 0.4

crude SIS Price VRF 0.6137 2.3 0.6148 1.1 0.2547 1.4 0.2556 0.9 0.0662 2.1 0.0661 1.1

EF 1.0 0.6 0.3 0.3 0.3 0.2

Price 0.6148 0.6154 0.2551 0.2552 0.0660 0.0661

SIS+ST VRF 42 26 43 37 45 21

EF 19 11 8 6 6 2

Price 0.6151 0.6151 0.2552 0.2553 0.0661 0.0661

SIS+RV VRF 696 8139 341 21275 446 8667

EF 309 2963 65 3051 56 843

SIS+BB Price VRF 0.6150 144 0.6151 253 0.2552 121 0.2552 179 0.0660 125 0.0660 333

EF 63 102 23 30 16 37

CR IP T

(1)

BG DG BG DG BG DG

Price 0.6154 0.6170 0.2552 0.2545 0.0659 0.0657

AN US

Case

M

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes AL (Avramidis and L’Ecuyer (2006)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.0548, c = 0, θ = −0.2859, σVG = 0.1927, ν = 0.2505, T = 0.40504, n = 214 and R = 100, and the other parameters for the three cases are given in Table 5. The computational time is reported in seconds.

PT

ED

where both functions a and b map R × [0, +∞) onto R and S (0) = S 0 . If the coefficients a and b satisfy the local Lipschitz and linear growth conditions, the scalar SDE has a unique solution. The Lipschitz condition can be relaxed considerably. Particularly, Yamada and Watanabe (1971) established the uniqueness under the following conditions R 1. there exists a strictly increasing function ρ(u) on [0, ∞) such that ρ(0) = 0, 0+ ρ−2 (u)du = ∞ and |b(x, t) − b(y, t)| ≤ ρ(|x − y|) for all t ∈ [0, T ] and all x, y ∈ R; 2. there exists an increasing and concave function κ(u) defined on [0, ∞) such that κ(0) = 0, R −1 κ (u)du = ∞ and |a(x, t) − a(y, t)| ≤ κ(|x − y|) for all t ∈ [0, T ] and all x, y ∈ R. 0+

AC

CE

We refer to Ikeda and Watanabe (1989) for the detailed discussions. We take the Euler discretization to the SDE (28), i.e., the asset price at time t is approximated by Sˆ (t) = Sˆ (s) + a(Sˆ (s), s)(t − s) + b(Sˆ (s), s)(B(t) − B(s)), 0 ≤ s < t ≤ T, where B(t) − B(s) ∼ N(0, t − s). The approximated asset price Sˆ (t) is explicitly determined by the approximated asset price at some previous time s and a normal variable. Note that the relationship between Sˆ (t) and the normal variable is linear, which ensures viability and simplicity of working out the conditional probability     H−Sˆ (s)+a(Sˆ (s),s)(t−s)  √  Φ , b(Sˆ (s), s) > 0,   b(Sˆ (s),s) t−s  P(Sˆ (t) ≤ H|Sˆ (s)) =   H−Sˆ (s)+a(Sˆ (s),s)(t−s)   √ , b(Sˆ (s), s) < 0. 1 − Φ ˆ b(S (s),s) t−s 22

ACCEPTED MANUSCRIPT

Table 7: Numerical results for down-and-out barrier option under the VG model

(2) (3)

Case (1) (2) (3)

BG DG BG DG BG DG

crude SD 3.36e-2 1.41e-2 3.00e-2 2.94e-2 4.58e-2 3.75e-2

Time 34 61 98 165 158 263

Price 9.9065 9.9061 3.9351 3.9350 6.1605 6.1599

SIS+ST VRF 2242 16 1434 64 347 14

EF 986 7 279 11 46 2

AL VRF 30.2 6.3 20.4 28.7 14.8 10.8

EF 12.7 2.5 3.9 4.8 2.0 1.3

crude SIS Price VRF 9.9032 2.1 9.9060 4.7 3.9395 1.6 3.9332 7.4 6.1610 1.6 6.1595 7.0

EF 0.9 2.7 0.3 2.1 0.2 2.1

SIS+RV Price VRF 9.9066 52 9.9064 21 3.9350 36 3.9350 62 6.1599 25 6.1603 24

EF 23 8 7 9 3 3

SIS+BB Price VRF 9.9065 170 9.9062 17 3.9354 171 3.9350 96 6.1602 94 6.1602 26

EF 74 7 33 16 12 3

Price 9.9067 9.9073 3.9347 3.9355 6.1588 6.1602

CR IP T

(1)

BG DG BG DG BG DG

Price 9.9026 9.9040 3.9316 3.9367 6.1533 6.1608

AN US

Case

M

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes AL (Avramidis and L’Ecuyer (2006)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.0548, c = 0, θ = −0.2859, σVG = 0.1927, ν = 0.2505, T = 0.40504, n = 214 and R = 100, and the other parameters for the three cases are given in Table 5. The computational time is reported in seconds.

AC

CE

PT

ED

Based on this, one can easily obtain the likelihood ratio and then the IS-based estimator (11). Note that the Euler discretization is always performed in the forward order. Consequently, the asset prices are always generated in the forward order. This leaves no room for optimizing the first generation step in such a model as in the BS model and the VG model. As observed in Sections 6.1 and 6.2, using SIS alone in QMC does not necessarily lead to high efficiency. The SIS combining with a good choice of the first generation step has the effect of smoothing and dimension reduction, leading to large efficiency improvement. Choosing a good generation order in the BS model and the VG model is a way of dimension reduction (see Appendix A in the supplementary material). Motivated by this, instead of choosing a good generation order, we suggest to use some dimension reduction techniques to reduce the effective dimension of the crude SIS estimator. There are some dimension reduction techniques in the literature, such as the LT method (Imai and Tan, 2006), the clustering algorithm based QR method (Weng et al., 2017), and the gradient PCA (GPCA) method (Xiao and Wang, 2017). These methods are generally designed for smooth functions, not for discontinuous functions. Thanks to the smoothing effect of the SIS method, we can combine the SIS method with a dimension reduction technique. We take the CEV model as an example to illustrate the effectiveness of the crude SIS method and the SIS method combined with dimension reduction techniques. The CEV model is considered as an important alternatives to the lognormal model for an asset price (see Glasserman, 2004). In this framework, the asset price S (t) is assumed to follow the SDE dS (t) = rS (t) dt + σCEV S (t)β/2 dB(t), where r, σCEV and β > 0 are constants, and B(t) is the standard Brownian motion. In particular, 23

ACCEPTED MANUSCRIPT

Table 8: Numerical results for double-barrier option under the VG model

(2) (3)

Case (1) (2) (3)

BG DG BG DG BG DG

crude SD 1.20e-2 1.13e-2 6.19e-3 5.70e-3 1.47e-3 1.71e-3

Time 33 59 96 166 149 255

Price 0.5193 0.5197 0.2172 0.2173 0.0135 0.0135

SIS+ST VRF 36 30 41 48 52 26

EF 15 12 8 8 7 3

AL VRF 2.9 3.3 3.5 3.5 1.4 2.6

EF 1.2 1.4 0.7 0.6 0.2 0.3

crude SIS Price VRF 0.5187 4.7 0.5199 3.6 0.2170 2.0 0.2166 2.9 0.0135 5.9 0.0135 6.0

EF 2.0 1.7 0.4 0.8 0.7 1.1

SIS+RV Price VRF 0.5197 84 0.5197 187 0.2173 71 0.2173 198 0.0135 67 0.0135 108

EF 36 67 13 28 8 10

SIS+BB Price VRF 0.5195 92 0.5196 275 0.2173 89 0.2172 153 0.0135 63 0.0135 156

EF 39 110 17 26 8 17

Price 0.5191 0.5207 0.2170 0.2172 0.0135 0.0136

CR IP T

(1)

BG DG BG DG BG DG

Price 0.5195 0.5206 0.2175 0.2167 0.0136 0.0134

AN US

Case

M

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of schemes AL (Avramidis and L’Ecuyer (2006)), crude SIS, and SIS combined with different generation orders (ST, RV, BB) relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.0548, c = 0, θ = −0.2859, σVG = 0.1927, ν = 0.2505, T = 0.40504, n = 214 and R = 100, and the other parameters for the three cases are given in Table 5. The computational time is reported in seconds.

AC

CE

PT

ED

the special case β = 2 is exactly the BS model. Some empirical studies have found that β < 2 gives a better fit to stock price data (see Glasserman, 2004), so we are only interested in the case β < 2. In the numerical examples, we take β = 1.5. To facilitate comparison between the results of the BS model and the CEV model, we still use the same parameter setting in BS model, except the parameter σCEV which is adjusted to satisfy σBS S 0 = σCEV S 0β/2 . Tables 9 and 10 show the numerical results of pricing the up-and-out barrier option and the down-and-out barrier option under the CEV model using the general dimension reduction techniques LT and GPCA, the crude SIS method, the SIS method combined with LT (SIS+LT), and the SIS method combined with GPCA (SIS+GPCA). Note that when directly using LT and GPCA for the unsmoothed payoff functions, the methods are only applied to the smooth structure S d − K within the payoff function (please refer to the very beginning of Section 6 for specific payoff functions), while when combining them with the SIS method, the methods are applied to the whole smoothed estimators. The observations are as follows. • The schemes that combine the SIS method with the general dimension reduction methods (SIS+LT and SIS+GPCA) always show absolute advantage over other simulation schemes. Although the EFs are usually half cut compared to the VRFs due to the relatively higher complexity, they are still much higher than those with other schemes. By using SIS, the estimators are smoothed and thus these dimension reduction methods (LT and GPCA) can play a better role. • The crude SIS method without dimension reduction brings small variance reduction under 24

ACCEPTED MANUSCRIPT

the CEV model with the Euler discretization. It shows comparable effectiveness on the CEV model and the BS model.

CR IP T

• When we directly use the general dimension reduction methods (LT and GPCA) without SIS, the corresponding EF is moderately higher than that of the crude QMC, but it is never better than that of combining dimension reduction methods with SIS. We thus do not recommend to use these dimension reduction methods in this way, as these methods are designed for smooth functions.

• The LT method and the GPCA method consume similar computational cost. Regarding their effectiveness in improving QMC efficiency, without SIS the LT method exhibits better performance than GPCA. When combining with SIS, GPCA usually gives higher VRF and EF than LT.

AN US

Remark 5. The general dimension reduction techniques (such as LT and GPCA) can also be applied to the crude SIS estimator under the BS model and the VG model. However, these methods are much more time consuming than just choosing a good generation order as we did in Section 6.1 and 6.2. We thus focused there on combining SIS with approaches of choosing the optimal generation order.

Table 9: Numerical results for up-and-out barrier option under the CEV model

Case

Price 0.1169 0.1896 0.0211 0.3428 0.0501 0.0522

LT VRF 3.4 3.8 3.3 3.5 2.0 3.9

EF 2.1 2.3 2.0 2.1 1.2 2.4

Price 0.1171 0.1891 0.0212 0.3420 0.0503 0.0520

SIS+LT VRF 21.5 24.1 54.8 10.0 8.8 18.7

EF 10.7 12.0 27.3 5.1 4.4 9.5

M

crude SIS VRF 2.7 1.8 4.8 1.6 2.0 2.3

CE

(1) (2) (3) (4) (5) (6)

Price 0.1171 0.1886 0.0210 0.3419 0.0501 0.0520

crude QMC SD Time 6.56e-3 12 6.27e-3 12 2.00e-3 12 9.66e-3 25 3.27e-3 26 2.85e-3 26

ED

(1) (2) (3) (4) (5) (6)

Price 0.1174 0.1884 0.0207 0.3412 0.0501 0.0515

PT

Case

EF 2.0 1.3 3.6 1.2 1.5 1.8

Price 0.1174 0.1893 0.0214 0.3418 0.0501 0.0524

GPCA VRF 1.0 1.2 1.2 1.9 1.3 1.8

EF 0.6 0.8 0.7 1.2 0.8 1.1

SIS+GPCA Price VRF EF 0.1170 49.7 25.2 0.1892 49.0 24.5 0.0212 48.0 23.7 0.3419 33.0 16.7 0.0501 17.2 8.7 0.0521 28.4 14.5

AC

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of general dimension reduction techniques (LT, GPCA), the crude SIS scheme and SIS combined with these general dimension reduction techniques relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.01, β = 1.5, T = 0.25, 1−β/2 n = 214 and R = 100, and the other parameters for the six cases are given in Table 1 with σCEV = σBS S 0 . The computation time is reported in seconds.

25

CR IP T

ACCEPTED MANUSCRIPT

Table 10: Numerical results for down-and-out barrier option under the CEV model

(1) (2) (3) (4) (5) (6)

Price 7.2151 5.3237 6.4189 6.2749 6.8573 3.7936

crude SIS VRF 1.8 1.2 1.2 1.5 1.5 1.0

EF 1.3 0.9 0.9 1.2 1.1 0.8

Price 7.2229 5.3257 6.4241 6.2694 6.8529 3.7920

LT VRF 3.7 21.8 16.7 5.2 5.4 30.7

EF 2.2 13.3 10.2 3.1 3.3 18.6

Price 7.2178 5.3262 6.4252 6.2719 6.8540 3.7920

SIS+LT VRF 43.2 82.9 105.5 29.8 20.1 79.6

EF 21.5 40.8 51.3 14.9 10.0 39.0

Price 7.2233 5.3247 6.4225 6.2679 6.8537 3.7925

AN US

Case

crude QMC SD Time 7.00e-2 12 5.23e-2 12 8.61e-2 12 4.83e-2 25 7.55e-2 26 3.80e-2 26

M

(1) (2) (3) (4) (5) (6)

Price 7.2158 5.3241 6.4207 6.2737 6.8628 3.7938

ED

Case

GPCA VRF 1.7 12.2 9.5 1.6 2.3 8.9

EF 1.1 7.5 5.8 1.0 1.4 5.4

SIS+GPCA Price VRF EF 7.2179 65.9 32.8 5.3261 220.3 108.6 6.4236 243.1 110.2 6.2725 43.3 21.7 6.8526 51.8 25.8 3.7927 158.5 78.0

AC

CE

PT

Note. The standard deviation (SD) and the computational time of the crude QMC scheme, the variance reduction factor (VRF) and the efficiency gain (EF) of general dimension reduction techniques (LT, GPCA), the crude SIS scheme and SIS combined with these general dimension reduction techniques relative to the crude QMC, and the price estimation under all schemes are reported. The fixed parameters are S 0 = 100, r = 0.01, β = 1.5, T = 0.25, 1−β/2 n = 214 and R = 100, and the other parameters for the six cases are given in Table 1 with σCEV = σBS S 0 . The computation time is reported in seconds.

26

ACCEPTED MANUSCRIPT

7. Concluding Remarks

AN US

CR IP T

In this paper, we made an attempt to deal with discrete barrier option pricing, where the discontinuous structures in the payoff functions become obstacles for using QMC. We developed an IS-based smoothing method which can remove the discontinuities of the payoff function to improve the efficiency of QMC. The proposed SIS method includes the conditional on one-step survival method proposed by Glasserman and Staum (2001) as a special case. We went a step further and found that the path generation order has a substantial influence on the efficiency. Hence we made an attempt to minimize the variance contributed at the first step. We focused on the barrier option pricing problems under the BS model and the class of asset price models based on subordinated Brownian motion. We designed a PGM which combines the IS-based estimator with an optimal first generation step. With a proper generation order in subsequent steps, the effective dimension of the SIS estimate can be reduced greatly. Numerical experiments have demonstrated the great power of the proposed method. For future work, we plan to study further a general rule to determine the whole generation order, beyond the first step. It is also interesting to extend the SIS method to general models. As an illustrative extension, we showed that the SIS method can be applied to the discretized process (using Euler scheme). However, we can not apply the proposed combined procedure due to the prescribed generation order of a discretized process. In this situation, an alternative scheme is to combine the SIS method with general dimension reduction techniques. Our numerical results showed that general dimension reduction techniques help in efficiency improvement of the IS-based estimator.

M

Acknowledgments

PT

References

ED

The authors gratefully thank the editor and anonymous reviewers, whose comments on the paper greatly improved its quality. Zhijian He was supported by the National Science Foundation of China (No. 71601189) and the Research Funds of South China University of Technology (No. K5180420). Fei Xie and Xiaoqun Wang was supported by the National Science Foundation of China (No. 71471100) and the National Key R&D Program of China (No. 2016QY02D0301).

AC

CE

Achtsis, N., Cools, R., Nuyens, D., 2013a. Conditional sampling for barrier option pricing under the Heston model. In: Monte Carlo and Quasi-Monte Carlo Methods 2012. Springer, pp. 253–269. Achtsis, N., Cools, R., Nuyens, D., 2013b. Conditional sampling for barrier option pricing under the LT method. SIAM J. Financ. Math. 4 (1), 327–352. Avramidis, A. N., L’Ecuyer, P., 2006. Efficient Monte Carlo and quasi-Monte Carlo option pricing under the variance gamma model. Manage. Sci. 52 (12), 1930–1944. Birge, J., Linetsky, V., 2007. Discrete barrier and lookback options. Handb. Oper. Res. Manage. Sci. 15, 343. Broadie, M., Glasserman, P., Kou, S., 1997. A continuity correction for discrete barrier options. Math. Finance 7 (4), 325–349. Caflisch, R. E., Morokoff, W. J., Owen, A. B., 1997. Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension. J. Comput. Financ. 1 (1), 27–46. Carr, P., Madan, D., 1999. Option valuation using the fast Fourier transform. J. Comput. Financ. 2 (4), 61–73. Chen, Z., Feng, L., Lin, X., 2012. Simulating L´evy processes from their characteristic functions and financial applications. ACM Trans. Model. Comput. Simul. 22 (3), 14:1–14:26. Cont, R., Tankov, P., 2004. Financial Modelling with Jump Processes. Vol. 2. CRC Press, Boca Raton. Dai, T.-S., Lyuu, Y.-D., 2010. The bino-trinomial tree: A simple model for efficient and accurate option pricing. J. Deriv. 17 (4), 7–24. Derman, E., Kani, I., 1996. The ins and outs of barrier options: Part 1. Deriv. Quarterly 2, 55–67.

27

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Derman, E., Kani, I., 1997. The ins and outs of barrier options: Part 2. Deriv. Quarterly 3, 73–80. Dia, E. H. A., Lamberton, D., 2011. Continuity correction for barrier options in jump-diffusion models. SIAM J. Financ. Math. 2 (1), 866–900. Dick, J., Pillichshammer, F., 2010. Digital Nets and Sequences: Discrepancy Theory and Quasi–Monte Carlo Integration. Cambridge University Press. Duan, J.-C., Dudley, E., Gauthier, G., Simonato, J.-G., 2003. Pricing discretely monitored barrier options by a Markov chain. J. Deriv. 10 (4), 9–31. Fang, F., Oosterlee, C. W., 2009. Pricing early-exercise and discrete barrier options by fourier-cosine series expansions. Numer. Math. 114 (1), 27–62. Feng, L., Linetsky, V., 2008. Pricing discretely monitored barrier options and defaultable bonds in L´evy process models: a fast Hilbert transform approach. Math. Financ. 18 (3), 337–384. Fuh, C.-D., Luo, S.-F., Yen, J.-F., 2013. Pricing discrete path-dependent options under a double exponential jump– diffusion model. J. Bank Financ. 37 (8), 2702–2713. Fusai, G., Abrahams, I. D., Sgarra, C., 2006. An exact analytical solution for discrete barrier options. Financ. Stoch. 10 (1), 1–26. Fusai, G., Germano, G., Marazzina, D., 2016. Spitzer identity, wiener-hopf factorization and pricing of discretely monitored exotic options. Eur. J. Oper. Res. 251 (1), 124–134. Genz, A., Bretz, F., 2009. Computation of Multivariate Normal and t Probabilities. Vol. 195. Springer Science & Business Media. Glasserman, P., 2004. Monte Carlo Methods in Financial Engineering. Springer. Glasserman, P., Staum, J., 2001. Conditioning on one-step survival for barrier option simulations. Oper. Res. 49 (6), 923–937. Golbabai, A., Ballestra, L., Ahmadian, D., 2014. A highly accurate finite element method to price discrete double barrier options. Comput. Econ. 44 (2), 153–173. Green, R., Fusai, G., Abrahams, I. D., 2010. The wiener–hopf technique and discretely monitored path-dependent option pricing. Math. Financ. 20 (2), 259–288. He, Z., Wang, X., 2014. Good path generation methods in quasi-Monte Carlo for pricing financial derivatives. SIAM J. Sci. Comput. 36 (2), B171–B197. He, Z., Wang, X., 2015. On the convergence rate of randomized quasi–Monte Carlo for discontinuous functions. SIAM J. Numer. Anal. 53 (5), 2488–2503. He, Z., Wang, X., 2017. Dimension reduction and smoothing in quasi–Monte Carlo for financial engineering. Preprint, arXiv:1709.02577. Heynen, R. C., Kat, H. M., 1994a. Crossing barriers. RISK 7 (June), 46–49. Heynen, R. C., Kat, H. M., 1994b. Partial barrier options. J. Financ. Engineering 3 (September/December), 253–274. H¨ormann, W., Leydold, J., 2003. Continuous random variate generation by fast numerical inversion. ACM Trans. Model. Comput. Simul. 13 (4), 347–362. Ikeda, N., Watanabe, S., 1989. Stochastic Differential Equations and Diffusion Processes. North-Holland. Imai, J., Tan, K. S., 2006. A general dimension reduction technique for derivative pricing. J. Comput. Financ. 10 (2), 129–155. Imai, J., Tan, K. S., 2014. Pricing derivative securities using integrated quasi–Monte Carlo methods with dimension reduction and discontinuity realignment. SIAM J. Sci. Comput. 36 (5), A2101–A2121. L’Ecuyer, P., 2009. Quasi-Monte Carlo methods with applications in finance. Financ. Stoch. 13 (3), 307–349. L’Ecuyer, P., Lemieux, C., 2005. Recent advances in randomized quasi-Monte Carlo methods. In: Dror, M., L’Ecuyer, P., Szidarovszky, F. (Eds.), Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications. Kluwer Academic Publishers, pp. 419–474. Lian, G., Zhu, S.-P., Elliott, R. J., Cui, Z., 2017. Semi-analytical valuation for discrete barrier options under timedependent L´evy processes. J. Bank Financ. 75, 167–183. Lindset, S., Persson, S.-A., 2016. A stochastic mesh size simulation algorithm for pricing barrier options in a jumpdiffusion model. J. Appl. Oper. Res. 8 (1), 15–25. Madan, D. B., Carr, P. P., Chang, E. C., 1998. The variance gamma process and option pricing. Eur. Financ. Rev. 2 (1), 79–105. Matouˇsek, J., 1998. On the L2 -discrepancy for anchored boxes. J. Complexity 14 (4), 527–556. Merton, R. C., 1973. Theory of rational option pricing. Bell J. Econ. Manage. Sci, 141–183. Moskowitz, B., Caflisch, R. E., 1996. Smoothness and dimension reduction in quasi-Monte Carlo methods. Math. Comput. Model. 23 (8), 37–54. Niederreiter, H., 1992. Random Number Generation and Quasi-Monte Carlo Methods. SIAM, Philadelphia. Rosenblatt, M., 1952. Remarks on a multivariate transformation. Ann. Math. Stat. 23 (3), 470–472. Wang, X., Tan, K. S., 2013. Pricing and hedging with discontinuous functions: Quasi-Monte Carlo methods and dimension reduction. Manage. Sci. 59 (2), 376–389.

28

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Weng, C., Wang, X., He, Z., 2017. Efficient computation of option prices and Greeks by quasi-Monte Carlo method with smoothing and dimension reduction. SIAM J. Sci. Comput. 39 (2), B298–B322. Xiao, Y., Wang, X., 2017. Enhancing quasi-Monte Carlo simulation by minimizing effective dimension for derivative pricing. Comp. Econ., To appear. Yamada, T., Watanabe, S., 1971. On the uniqueness of solutions of stochastic differential equations. Journal of Mathematics of Kyoto University 11 (1), 155–167. Zhu, Z., De Hoog, F., 2010. A fully coupled solution algorithm for pricing options with complex barrier structures. J. Deriv. 18 (1), 9–17.

29