Proceedings of the 2011 Winter Simulation Conference
S. Jain, R. R. Creasey, J. Himmelspach, K. P. White, and M. Fu, eds.

RARE EVENT SIMULATION TECHNIQUESJose Blanchet

Henry Lam

Columbia University
500 W 120th Street
New York, NY 10027, USA

Boston University
111 Cummington Street
Boston, MA 02215, USA

ABSTRACT
We discuss rare event simulation techniques based on state-dependent importance sampling. Classicalexamples and counter-examples are shown to illustrate the reach and limitations of the state-independent
approach. State-dependent techniques are helpful to deal with these limitations. These techniques can beapplied to both light and heavy tailed systems and often are based on subsolutions to an associated Isaacs
equation and on Lyapunov bounds.
1

INTRODUCTION

We study recent techniques that have been developed for the design of provably efficient rare event simulation
estimators. We concentrate mostly on importance sampling techniques. Our discussion here is based on
the survey paper Blanchet and Lam (2011).
We mention that besides importance sampling, another powerful set of tools for variance reduction
in rare event problems are based on splitting. Interested readers can consult L’Ecuyer et al. (2009) and
Chapter VI Section 9 of Asmussen and Glynn (2007) for basic notions on splitting.
The types of techniques that we consider in this paper focus on the interplay between large deviations
results and their applications to provably efficient importance sampling estimators. Early references in thisspirit include Siegmund (1976) and Sadowsky (1991). Applications of these techniques in settings such as
insurance and queueing include Asmussen (1985), Collamore (2002), Asmussen and Rubinstein (1995),
and Parekh and Walrand (1989). We refer the reader to Juneja and Shahabuddin (2006), Bucklew (1990)
and Chapter VI of Asmussen and Glynn (2007) for a wider scope on importance sampling for rare event
simulation.
A series of examples and counter-examples (Glasserman and Kou 1995, Glasserman and Wang 1997)
related to the use of large deviations principles in the design of efficient importance sampling estimatorsfor light-tailed systems has motivated research on state-dependent techniques. In addition, as explained in
Bassamboo et al. (2006), efficient importance sampling estimators in the setting of heavy tails often must
be state-dependent (see also Asmussen et al. 2000 for more on the challenges that arise in the context ofefficient rare-event simulation for heavy-tailed systems).
In order to systematically address the construction of provably efficient state-dependent importancesampling estimators for light-tailed systems P. Dupuis, H. Wang and their students have developed a method
based on control theory and the use of subsolutions of an associated Isaacs equation for constructing andtesting importance sampling estimators (Dupuis and Wang 2004; Dupuis and Wang 2007; Dupuis et al.
2007; Dupuis et al. 2007). A related approach, based on the construction of Lyapunov inequalities has
been also developed (Blanchet and Glynn 2008; Blanchet and Liu 2008; Blanchet et al. 2007) for theconstruction and analysis of state-dependent importance sampling estimators for heavy-tailed systems.
In this paper we review the development of state-dependent importance sampling estimators and illustrate
their application. We start by discussing standard notions of statistical efficiency that are used in order torigorously test the optimality of a simulation estimator. This will be done in Section 2. The introduction of978-1-4577-2109-0/11/$26.00 ©2011 IEEE

146

Blanchet and Lamthese notions of efficiency leads to the question of how to construct importance sampling estimators thatcan be rigorously shown to be efficient. The zero-variance change-of-measure, while virtually impossible
to implement in most situations of interest, provides a good guidance in order to construct such estimators.
The zero-variance change-of-measure turns out to be the conditional distribution of the underlying
process given the rare event of interest. So, it is natural to leverage-off the available asymptotic theory
for rare events in order to understand such a conditional distribution. This is basically the starting pointin the design of efficient importance sampling estimators for rare events. The idea is to take advantageof the available large deviations theory, which provides insight into the occurrence of rare events. In thelight-tailed case, the theory of large deviations is comprehensively documented in Dembo and Zeitouni(1998). References in connection to large deviations for heavy-tailed problems (mainly related to random
walk settings) can be found in Embrechts et al. (1997) and Chapter IX of Asmussen (2000).
The successful application of large deviations theory in the design of efficient importance sampling
estimators is not straightforward and this feature turned out to be the norm rather than the exception,
especially in applications such as the analysis of queueing networks (see Glasserman and Kou 1995). To
avoid such poor variance behavior, Dupuis and Wang formulate the design of state-dependent importancesampling estimators in terms of a stochastic control problem in which the value function is the second
moment of the estimator. Taking an asymptotic limit, as we shall see in Section 3, the stochastic control
problem becomes a deterministic control problem. The application of Dupuis and Wang’s method is
simplified thanks to the fact that only a subsolution to the value function is needed for the design of an
efficient sampling scheme. A trivial supersolution is given, thanks to Jensen’s inequality, by twice theunderlying large deviations rate function. The whole point of Dupuis and Wang’s approach is to constructa subsolution that matches the trivial supersolution at one point, namely, the initial condition of the system
under consideration. In many problems of interest it turns out that such a subsolution is piecewise affine
and therefore relatively easy to manipulate.
We will then move to the analysis of heavy-tailed systems, which will be introduced in our last section,
namely Section 4. These types of systems arise in the analysis of finance and insurance models. There arefundamental qualitative differences between light and heavy tailed systems when it comes to the analysis ofrare events. This is why it is incorrect to use a model with light-tailed components to study large deviations
properties of systems that are built from heavy-tailed building blocks. Light-tailed large deviations occur
gradually, while heavy-tailed large deviations occur suddenly and this feature manifests itself in the designof efficient importance sampling estimators for heavy-tailed systems.
We emphasize that this paper aims to introduce readers the key concepts and ideas in rare-event
problems. Throughout the exposition we will drop many of the mathematical details and highlight only
the central arguments, pointing readers to other references. Finally, we briefly mention the existence ofsome empirical approaches aimed at selecting good importance sampling distributions. We do not discussthese methods in detail but rather refer the reader to some related references including: adaptive Monte
Carlo (Kollman, Baggerly, Cox, and Picard 1999) and the cross-entropy method (Rubinstein and Kroese
2004). In this vein, we mention the recent work of Chan and Lai (2010) who considers an interesting
importance sampling approach for Markov-modulated problems that avoids the need for computing theassociated eigenfunctions.
2

NOTIONS ON EFFICIENCY AND IMPORTANCE SAMPLING

2.1 Notions of Efficiency
In order to design algorithms that achieve high accuracy in relative terms and that perform well in rare
event settings it is useful to embed the question of interest in a suitable asymptotic regime, often motivatedby an applied standpoint. For example, in the setting of bankruptcy of an insurance company, it is often
the case that the capital or reserve of the company is large relative to individual claim sizes, and this isbasically why one might expect the probability of bankruptcy, P (A), to be small. So, one might introduce

147

Blanchet and Lama parameter, say n, for instance the initial capital in the insurance example, and consider αn = P (An ) asa function of the parameter n which is often referred to as the “rarity parameter”. From now on we willassume that underlying the rare event there is a rarity parameter n such that αn −→ 0 as n −→ ∞.
The idea is then to design a simulation estimator whose error is controlled in relative terms as n % ∞compared to naive (or crude) Monte Carlo, which provides an obvious alternative and, therefore, a generalbenchmark. Note that naive Monte Carlo involves simulating N (independent and identically distributed)
iid Bernoullies, In (1), ..., In (N ), with success parameter αn , thereby producing an estimate of the formPN
αbnc (N ) =The super-index “c” that appears in α
bnc (N ) corresponds to crude Monte
j=1 In (j)/N .
Carlo. Since α
bnc (N ) is an unbiased estimator of αn , its mean squared error equals V ar(bαnc (N ))1/2 =
1/2
αn (1 − αn )1/2 /N 1/2 . In turn, we obtain a relative mean squared error equal to the coefficient of variation
1/2of α
bnc (N ), namely, CV (bαnc (N )) = (1 − αn )1/2 /(αn N 1/2 ). In order to control the relative error of αbnc (N )one needs N = Ω(1/αn ) Bernoullies.
The problem with the estimator α
bnc (N ) is that each of the underlying replications, the In (j)’s, havea very big variance in relative terms. The number of replications N must grow at the same rate as
1/2
CV (In (j)) = (1 − αn )1/2 /αn . The same phenomenon occurs in any estimator obtained out of averaging
iid replications, as in the definition of αbnc (N ). The objective is then to design an estimator Rn with acontrolled mean squared error. We concentrate only on unbiased estimators, so controlling the behavior of
the mean squared error boils down to controlling the behavior of the second moment. The overall estimate
is then obtained by averaging N iid copies of Rn . Motivated by these considerations we have the following
definitions.

Definition 1 An estimator Rn is said to be strongly efficient if ERn2 = O αn2 as n % ∞.
Definition 2 Rn is said to be weakly efficient or asymptotically optimal if for each ε > 0 we have thatERn2 = O αn2−ε as n % ∞.
In most cases the analysis of importance sampling estimators (especially in light-tailed cases) concludes
only weak efficiency of the estimators.
In order to see the gain in efficiency, suppose that an estimator Rn is strongly efficient and considerPαen (N ) = N
i=1 Rn (i) /N where the Rn (i)’s are independent replications of Rn . It follows from Chebyshev’s inequality that at most N = O(ε−2 δ −1 ) (uniformly in n) replications are required to conclude thatαen (N ) is ε-close to αn in relative terms with at least 1 − δ confidence.
Besides variance control, it is worth pointing out the other measures of good importance sampling
estimators, and we refer interested readers to L’Ecuyer et al. (2010) for these considerations (including
work-normalized measures to account for cost-per replication).
2.2 Basic Properties of Importance Sampling Estimators
Importance sampling is a variance reduction technique that is often applied to design an estimator, Rn , withefficiency characteristics such as those described by Definitions 1 and 2. We shall review basic properties ofimportance sampling next; for a more detailed discussion see, for example, Rubinstein and Kroese (2009)
and Asmussen and Glynn (2007). The basic idea is to introduce a probability measure Pe (·) such that thelikelihood ratio or Radon-Nikodym derivative between the nominal (original) probability measure, P (·),
and Pe (·) is well defined on the event An of interest. Then we can simply let
Rn =dP(ω) I (ω ∈ An ) ,dPe

as usual, we use ω to denote the underlying random outcome which is simulated according to the probabilitye (·) to denote the expectation operator associated to the probability measure Pe (·).
measure Pe (·). We use E

148

Blanchet and Lam
Note that

Z
e n=
ER

dP
(ω) dPe (ω) = P (An ) = αn ,e
An dP

so Rn is an unbiased estimator, and
2
Z 
Z
dP
dP
e n2 ) =E(R(ω) dPe (ω) =(ω) dP (ω) = E(Rn ).
e
e
dP
An
An dP

(1)

In principle one can design an importance sampling estimator with zero variance. Indeed, if we let
Pe (·) = Pn∗ (·) := P (· | An ), thenRn∗ :=dP(ω) I (ω ∈ An ) = P (An ) ,dPn∗which is clearly an exact estimate of αn = P (An ). Of course, assuming that one has access to P ∗ (·)
defeats the point of using simulation to estimate P (An ). However, the underlying lesson behind thecharacterization of the zero-variance change-of-measure as the conditional distribution given the rare eventof interest is that one can use asymptotic theory to describe Pn∗ (·) as n % ∞.
3

STATE-DEPENDENT TECHNIQUES FOR LIGHT-TAILED SYSTEMS

3.1 Solutions to a Deterministic Control Problem and Efficient Importance Sampling
As we indicated in the Introduction, Dupuis and Wang introduced a technique to guarantee efficiency (inthe sense of Definition 2) of importance sampling estimators based on control theoretic considerations.
We will formally characterize the solution to a deterministic control problem proposed by Dupuis andWang (2004). In this section we assume that (Xk : k ≥ 1) is a sequence of iid rv’s taking values in Rd
and such that ψ(θ) = log E exp(hXk , θi) is finite everywhere, where h·, ·i is defined as the inner productunder Euclidean norm. We then let Sk = X1 + ... + Xk and consider, for some y ∈ Rd ,Yn (t) = y + Sbntc /n.

(2)

Let B be a closed set with non-empty interior and define τB (n) = inf{t > 0 : Yn (t) ∈ B}. We
assume that B is attainable in the sense that P (τB (n) < ∞|Yn (0) = y) > 0 for each y ∈ Rd .
Consider two disjoint sets B and C. We are interested in computing
αn (y) = Py (τB (n) < τC (n), τB (n) < ∞) := P (τB (n) < τC (n), τB (n) < ∞|Yn (0) = y).
We assume that B and C are regular in the sense that limn→∞ n−1 log αn (y) = −IB,C (y), where
IB,C (y) is the associated large deviations rate function (which we do not need to specify but can be
consulted in Dembo and Zeitouni 1998, p. 176).
We are interested in choosing a suitable importance sampling estimator with optimal performance inthe sense of minimizing its second moment. The class of policies that we consider are based on exponential
tilting. In particular, the increment Xk is ultimately simulated according to the distributionPθk (Xk ∈ dx) = exp(hθk , Xk i − ψ(θk ))P (Xk ∈ dx).

(3)

If Xk has a density with respect to the Lebesgue measure, say f (·), then P (Xk ∈ dx) = f (x) dx. The
notation adopted in (3) is simply a mechanism that allows us to consider general distributions (including
discrete or mixtures of discrete and continuous distributions). We will use Eθk (·) to denote the expectation
operator associated to Pθk (·). The parameter θk might depend on (X1 , ..., Xk−1 ) but it is not allowed todepend on future observations. We say that θk is non-anticipating. As examples of exponentially tilted
149

Blanchet and Lam
distributions consider the following: If Xk is a one dimensional standard Gaussian under P (·), then, underPθk (·), Xk is normal with mean θk and unit variance; on the other hand, if under Xk is exponentiallydistributed with unit mean under P (·), then under Pθk (·) (for θk < 1), Xk is exponentially distributed withmean 1/(1 − θk ).
One reason that motivates the use of exponentially tilted distributions is that one can sometimes showthat the distribution of (X1 , ..., Xk ) given a large deviations event of interest converges to a suitable optimal
exponential tilting as n → ∞ (see, for instance, Asmussen and Glynn 2007, Section VI.5).
Using exponential tilting the estimator ultimately takes the form
Rn = exp(−nτB (n)nτB (n)XXhθk , Xk i +k=1

ψ(θk ))I(τB (n) < τC (n), τB (n) < ∞).

(4)

k=1

It follows immediately that the HJB (Hamilton-Jacobi-Bellman) equation corresponding to finding the
optimal policy (i.e., the optimal sequence of θk ’s) to minimize the second moment of Rn takes the form
Vn (y) = inf Eθ [exp(−2hθ, Xi + 2ψ(θ))Vn (y + X/n)]θ= inf E[exp(−hθ, Xi + ψ(θ))Vn (y + X/n)]θ(5)

for y ∈
/ B ∪ C and subject to the boundary condition that Vn (y) = 1 for y ∈ B and Vn (y) = 0 if y ∈ C.
If there is an optimal policy (θk∗ : k ≥ 1) generated by the previous HJB equation, then generally it would
be the case thatnτB (n)nτB (n)X

X

Vn (y) = E exp(−hθk∗ , Xk i +

k=1

k=1

nτB (n)nτB (n)

XX≤ E exp(−k=1

hθk , Xk i +ψ(θk∗ ))I(τB (n) < τC (n), τB (n) < ∞)ψ(θk ))I(τB (n) < τC (n), τB (n) < ∞),k=1

for any non-anticipating policy (θk : k ≥ 1).
The large deviations scaling suggests writing Vn (y) = exp(−nHn (y)) and thus we should expect
Hn (y) → H(y) as n → ∞ for some function H (y). If we proceed using this postulated limit in theprevious HJB equation, after taking logarithms of (5), formally and without being careful about underlying
smoothness assumptions and errors incurred, we arrive at the approximation−nH(y) ≈ min log E[exp(−hθ, Xi + ψ(θ) − nH(y + X/n))]θ(6)≈ min log E[exp(−hθ, Xi + ψ(θ) − nH(y) − h∇H(y), Xi)].
θ

Equivalently, we have that0 ≈ min log E[exp(−hθ, Xi + ψ(θ) − h∇H(y), Xi)]θ= min log exp(ψ(θ) + ψ(−∇H(y) − θ)) = min[ψ(θ) + ψ(−∇H(y) − θ)].
θ

θ

(7)

First order optimality conditions imply that at the optimal value θ∗ (y) one has ∇ψ(θ∗ (y)) =∇ψ(−∇H(y) − θ∗ (y)), which yields θ∗ (y) = −∇H(y)/2 and therefore we conclude that equation(7) can be expressed as2ψ(−∇H(y)/2) = 0,
(8)

150

Blanchet and Lam
subject to the boundary conditions (inherited from (5)) equal to H(y) = 0 for y ∈ B and H(y) = ∞ fory ∈ C. Equation (8) corresponds to a so-called Isaacs equation. The name originates from a connectionthat can be established to a solution of a two-player differential game, see Dupuis and Wang (2004).
On the other hand, we have assumed αn (y) = exp(−nIB,C (y) + o(n)) as n → ∞ and, by conditioning
on the first increment of the random walk, we also have thatαn (y) = Eαn (y + X/n)

(9)subject to the constraints that αn (y) = 1 for y ∈ B and αn (y) = 0 if y ∈ C. Proceeding to analyze
equality (9) formally as we did for the discussion leading to (8), we conclude thatexp(−nIB,C (y) + o(n)) = E exp(−nIB,C (y + X/n)) ≈ E exp(−nIB,C (y) − h∇IB,C (y), Xi).
Equivalently, taking logarithms we arrive at ψ(−∇IB,C (y)) = 0, subject to the boundary conditions
implied by αn (y), namely, IB,C (y) = 0 for y ∈ B and IB,C (y) = ∞ if y ∈ C. These considerations,
together with our analysis leading to (8) yield that H(y) = 2IB,C (y) and therefore, under suitable regularity
conditions we obtain that applying importance sampling with exponential tilting given by the tilting parameterθ∗ (y) = −∇IB,C (y) allows us to obtain an asymptotically optimal estimator (recall Definition 2). Sufficientconditions required to rigorously substantiate this result are given in Dupuis and Wang (2004), at least in
the case of large deviations problems for Sn /n.
The problem with (8) is that it is rarely applicable in most situations of interest, especially because it isoften the case that there is no classical (smooth) solution to the Isaacs equation (8). (It is apparent from ourheuristic derivation that smoothness appears to be an important feature.) The function H(y) = 2IB,C (y)typically provides a solution only in a weak sense. Fortunately, as we shall see next, only a subsolution isrequired, as long as it coincides with the value function at the initial condition.
3.2 Subsolutions to a Deterministic Control Problem and Efficient Importance Sampling
We take as starting point the HJB equation introduced in (5). The next lemma provides an inequality whose
solution gives an upper bound for the value function in (5); for a proof see Blanchet and Lam (2011).
Lemma 1 Suppose that one finds a non-negative function Un (·) such that
Un (y) ≥ E[exp(−hθ, Xi + ψ(θ))Un (y + X/n)] ≥ inf E[exp(−hθ, Xi + ψ(θ))Un (y + X/n)],θ

where y ∈/ B ∪ C and subject to the boundary condition that Un (y) ≥ 1 for y ∈ B. Then,
Un (y) ≥ E[exp(−

τBXhθj , Xj i +

j=1

τB
X

ψ(θj ))I(τB < ∞, τB < τC )]j=1

for any non-anticipating policy (θj : j ≥ 1).
Using a similar heuristic argument as the one introduced in the previous section in (6) we write
Un (y) = exp(−nGn (y)) and postulate Gn (y) → G(y) as n → ∞ for some function G (y). If we proceedusing this postulated limit in the inequality obtained in the previous lemma, after taking logarithms, wearrive at−nG(y) ' min log E[exp(−hθ, Xi + ψ(θ) − nG(y + X/n))]θ≈ min log E[exp(−hθ, Xi + ψ(θ) − nG(y) − h∇G(y), Xi)].
θ

151

Blanchet and Lam
We then conclude, as in the analysis leading to (7) that0 ≥ min[ψ(θ) + ψ(−∇G(y) − θ)].

(10)

θ

First order optimality conditions imply that at the optimal value θ∗ (y) one has ∇ψ(θ∗ (y)) =∇ψ(−∇G(y) − θ∗ (y)), which yields θ∗ (y) = −∇G(y)/2 and therefore we conclude that inequality(10) can be expressed as
0 ≥ 2ψ(−∇G(y)/2),(11)
for y ∈
/ B ∪ C and subject to the boundary conditions (inherited from Lemma 1) equal to G(y) ≤ 0 fory ∈ B. In order to conclude asymptotic optimality for αn (y0 ) we must have G(y0 ) ≥ 2IB,C (y0 ). Sufficientconditions required to rigorously substantiate this subsolution equation are given in Theorem 8.1 of Dupuis
and Wang (2007).
We now apply the previous result to an example, which is a generalization of a one dimensional ruin
problem.
(j)

Example 1 (multidimensional first passage time problems) Assume that d = 2 and let us write Xi todenote the j-th component of the i-th increment; in particular, we have that j = 1, 2. Similarly, we use(i)
Yn (t) to denote the i-th component of the process Yn (·) at time t, recall equation 2. Assume that Xi
follows a bivariate Gaussian distribution with mean EXi = (−1, −1)T and variance V ar (Xi ) = I (thatis, the components of Xi are independent). Therefore, we have that ψ (θ) = ||θ||2 /2.
Define the set B (i) = {y = (y (1) , y (2) )T : y (i) ≥ 1} and let B = B (1) ∪ B (2). Consider the firstpassage time probability
αn (y) = P (τB (n) < ∞|Yn (0) = y),where τB (n) = inf{t ≥ 0 : Yn (t) ∈ B}. We will construct an efficient estimator for αn (y) based on (11).
(i)
Similarly, we let τB(i) (n) = inf{t > 1 : Yn (t) ≥ 1} for each i = 1, 2, and note that one has theelementary inequalities
max P (τB(i) (n) < ∞|Yn (0) = y) ≤ αn (y) ≤
i=1,2

2
X

P (τB(i) (n) < ∞|Yn (0) = y).

(12)

i=1

Evidently, P (τB(i) (n) < ∞|Yn (0) = y) is equivalent to a one dimensional first passage time probability.
It is well known (see Asmussen 2003 p. 352) that
P (τB(i) (n) < ∞|Yn (0) = y) = exp[−nθi∗ (1 − y (i) ) + o(n)],(13)

where θi∗ = 2. We then obtain from (12) that limn→∞ n−1 log αn (y) = − mini=1,2 {θi∗ (1 − y (i) )}. The
value of θi∗ corresponds to the conditional asymptotic mean of the i − th component if the random walkhits B through B(i).
To apply the subsolution approach we need to find a function G(·) satisfying the inequality (11) together
with the corresponding boundary condition that G(y) ≤ 0 for y ∈ B. Moreover, since in the end we areafter an asymptotic upper bound for the second moment of the importance sampling estimator and we wish
to achieve the same decay rate as αn (y)2 , it makes sense to use the union bound (12) and (13) to postulateexp(−nG(y)) = (2X

exp[−nθi∗ (1 − y (i) )])2 .

i=1

We now verify that our selection is indeed a valid subsolution. First, note that ∇G(y) = −2
where a1 = (1, 0), a2 = (0, 1) and
exp[−nθj∗ (1 − y (j) )])pj (y, n) = P2.
∗(i)i=1 exp[−nθi (1 − y )]
152

P2

∗
i=1 θj aj pj (y),Blanchet and Lam
Note that for each y, (p1 (y, n), p2 (y, n)) can be interpreted as a probability vector. So, using Jensen’s
inequality and convexity of ψ (·) we conclude that
2
2
X
X
ψ(−∇G(y)/2) = ψ(θj∗ aj pj (y, n)) ≤pj (y, n)ψ(θj∗ aj ) = 0.
j=1

j=1

In the last equality we used that ψ(θj∗ aj ) = 0 for j = 1, 2. Now, if y ∈ B, then there exists i such that
1 − y (i) ≤ 0 and therefore evidently exp(−nG(y)) ≥ 0 or, equivalently, G(y) ≤ 0. We then conclude thatthe importance sampling scheme generated by using as exponential tiltingθk (Sk−1 /n) =2Xθj∗ aj pj (Sk−1 /n, n)

(14)

j=1

for k/n < τB (n) is asymptotically optimal; actually in this case the estimator can be shown to be strongly
efficient in the sense of Definition 1.
4

STATE-DEPENDENT TECHNIQUES FOR HEAVY-TAILED SYSTEMSAs mentioned earlier, various challenges of rare-event simulation with heavy tails are pointed out in
Asmussen et al. (2000) and Bassamboo et al. (2006). The former focuses on the singularity of theasymptotic conditional distribution with respect to the original measure and the difficulty of finding asuitable importance sampler, while suggesting provably efficient conditional Monte Carlo scheme. The
latter suggests the non-existence of provably efficient state-independent importance samplers for some firstpassage problems. Conditional Monte Carlo is also studied in Asmussen and Binswanger (1997) in thecontext of insurance, and further improved in Asmussen and Kroese (2006). A technique called hazard
rate twisting is studied in Juneja et al. (1999) and Juneja and Shahabuddin (2002).
Our goal here is to illustrate a technique, based on Lyapunov inequalities and appropriate mixture
samplers to estimate rare event probabilities. The Lyapunov approach was proposed in Blanchet and Glynn(2008) and it has subsequently applied to various settings (see, for instance, Blanchet et al. 2007; Blanchetand Liu 2008; Blanchet and Liu 2010b). The mixture distributions that we shall use were proposed in
Dupuis et al. (2006) in the context of power-law decaying tails; see also Blanchet and Liu (2010a) for
extensions beyond power-law decaying tails.
We first start with some basic results on heavy-tailed large deviations.
4.1 Large Deviations for Heavy-tailed Sums
For simplicity we concentrate on a special case of heavy-tailed random variables, namely, regularly varying
random variables, which posses power-law type tails; see Embrechts et al. (1997) for more information
on more general heavy-tailed models.
Definition 3 A random variable X with tail distribution F (·) = P (X > ·) = 1 − F (·) has a regularly
varying right tail with index α > 0 if F (βx) /F (x) −→ β −α as x % ∞ for all β > 0. Similarly, X hasa regularly varying left tail with index α > 0 if −X has a regularly varying right tail with index α > 0.
A function L (·) is slowly varying at infinity if L (βx) /L (x) −→ 1 as x → ∞ for each β > 0.
Therefore, if X has regularly varying right tail with index α > 0 we often write F (t) = t−α L (t).
The following result provides insight on how large deviations tend to occur in heavy-tailed models.
Proposition 1 Suppose that Sm = X1 + ... + Xm where the Xi ’s are iid non-negative regularly varying
rv’s. Then P (max1≤j≤m Xj > n|Sm > n) −→ 1 as n % ∞.
The previous proposition illustrates the so-called “catastrophe principle” behind extreme behaviorof heavy-tailed systems, which postulates that large deviations are caused by extremes in one or few
components. In this case, the sum is large because one single component, namely the largest one, is large.
153

Blanchet and Lam
In the heavy-tailed case there is a powerful heuristic that allows to “guess” a plausible form for large
deviations results and the associated “most likely paths” to the rare event. This heuristic, which is known
as the “fluid heuristic”, is based on elementary Law of Large Numbers analysis and basic facts of regularly
varying distributions. In order to guide our discussion we shall revisit a one dimensional analogue of theexample studied earlier in the light tailed case.
Example 2 (first passage time probabilities for heavy-tailed random walk)
Suppose that the Xj ’s
are regularly varying with index α > 1. Assume that EXk = −µ < 0 and let S0 = y. Defineτn = inf{k ≥ 0 : Sk > n} and consider the problem of approximating the ruin probability
αn (y) = P (τn < ∞|S0 = y)
as n % ∞.
The fluid heuristic consists in substituting the random walk by its fluid or Law of Large Numbers behavior
prior to the big jump that causes ruin. The associated fluid path in this case is Sbntc /n ≈ y(t) = y − µt
and therefore the fluid heuristic suggests
Z ∞∞X
αn (y) ∼P (Xk+1 + y(k) > n) ∼P (X > n + µt − y) dt0

k=0

as n % ∞, thereby neglecting the contribution of more than one jump and fluctuations beyond the Law
of Large Numbers. Letting u = n + µt − y we obtainZ
1 ∞P (X > u) dy.
αn (y) ∼µ n−y
as n % ∞. This approximation turns out to be correct even beyond the assumption of regularly varying
service times (see Pakes 1975; Veraverbeke 1977; Embrechts and Veraverbeke 1982).

4.2 Conditional Distributions and Implications For Simulation
As we discussed in Section 2, it is useful to understand the conditional distribution of the random walk given
the occurrence of the rare event of interest because such conditional distribution provides the optimal changeof-measure (in terms of variance minimization). A natural starting point would be to obtain a conditional
limit theorem for the Xk ’s given the rare event of interest. The asymptotic conditional distribution would
then be a natural candidate for an importance sampling distribution. One can pursue this idea in the context+
of Example 2. In such case, it is not difficult to verify that for any −∞ < x−
i < xi < ∞, with i = 1, ..., k
we have that




+
− +
k
lim P ∩ki=1 {Xi ∈ (x−
,
x
)}|τ
<
∞
=
P
∩
{X
∈
(x
,
x
)}
.
n
i
i=1
i
i
i
i
n→∞

Thus, concluding that the asymptotic conditional distribution of (X1 , ..., Xk ) given that τn < ∞ remains
unchanged as n → ∞. This is one of the main reasons that makes the design of efficient importancesampling estimators for heavy-tailed systems a challenging problem.
Let us now take a more direct approach to approximate the increment Xk+1 given the current position
Sk = s of the random walk at time k and given that k < τn < ∞. We note (adopting the notation introducedin (3)) that
P (Xk+1 ∈ dy|τn ∈ [k + 1, ∞), Sk = s)
= P (Xk+1 ∈ dy, τn = k + 1|τn ∈ [k + 1, ∞), Sk = s)
+ P (Xk+1 ∈ dy, τn > k + 1|τn ∈ [k + 1, ∞), Sk = s) .
154

Blanchet and Lam
Now,
P (Xk+1 ∈ dy, τn = k + 1|τn ∈ (k, ∞), Sk = s) = p∗n (s)P (Xk+1 ∈ dy) I(y > n − s),
P (Xk+1 > n − s|s)(15)

where p∗n (s) = P (Xk+1 > n − s) /αn (s). Using the approximation in Example 2, we have that
P (Xk+1 > n − s)µ(α − 1)p∗n (s) ∼ p(n − s) := µ R ∞∼(n − s)n−s P (Xk+1 > t) dt
as n − s → ∞. On the other hand, we have that
P (Xk+1 ∈ dy, τn > k + 1|τn ∈ [k + 1, ∞), Sk = s)
P (Xk+1 ∈ dy, τn ∈ [k + 2, ∞)|Sk = s) I(y ≤ n − s)=αn (s)P (Xk+1 ∈ dy) I(y ≤ n − s)αn (s + y)=× P (Xk+1 ≤ n − s).
P (Xk+1 ≤ n − s)αn (s)(16)

If we apply the fluid heuristic locally, that is by replacing Xk+1 by its mean, namely EXk+1 = −µ, wearrive at
αn (s − µ)(n − s + µ)−α+1µ(1 − α)1αn (s + Xk+1 )≈≈≈ [1 + µ/(n − s)]−α+1 = 1 −
+ o(
).
−α+1
αn (s)
αn (s)n−sn−s(n − s)Therefore, applying the previous local version of the fluid heuristic into (16) and combining this with (15)
we arrive at the non-rigorous approximation
P (Xk+1 ∈ dy) I(y > n − s)p(n − s)P (Xk+1 ≤ n − s)
P (Xk+1 ∈ dy) I(y ≤ n − s)+
[1 − p(n − s)].
P (Xk+1 ≤ n − s)

P (Xk+1 ∈ dy|τn ∈ [k + 1, ∞), Sk = s) ≈(17)

It turns out that this approximation is indeed valid in total variation (Blanchet and Liu 2010a) as n → ∞.
4.3 State-dependent Importance Sampling Techniques for Heavy-tailed Random Walks
We now discuss Lyapunov inequalities, which were introduced in Blanchet and Glynn (2008) in the context
of importance sampling, to show the efficiency of the family of samplers suggested by (17). This technique
is closely related and in some sense parallel to the subsolution approach explained earlier in the light-tailed
setting.
Lemma 2 Suppose that there exists r (x, y) > 0 satisfyingZ
r−1 (x, y)P (Sj ∈ dy|Sj−1 = x) = 1
for all x and j ∈ {1, 2, ...}. Then, we can define a Markov transition kernel via K(x, dy) = r−1 (x, y) P (Sj ∈dy|Sj−1 = x) and an importance sampling estimator of the form
R=

TAY
∧TB

r(Sj−1 , Sj )I(TA < TB ).

j=1

Suppose that there exists a non-negative function G (·) and a constant ρ ∈ (0, ∞) such thatE[r(x, S1 )G(S1 )|S0 = x] ≤ G(x)for x ∈/ A ∪ B and G(x) ≥ ρ for x ∈ A. Then, E [R|S0 = x] ≤ G(x)/ρ.
155

(18)

Blanchet and Lam
The proof of the previous result is completely analogous to Lemma 1 (see Blanchet and Lam 2011).
We now explain how to apply the previous result, combined with the fluid heuristics to design efficientimportance sampling estimators in the context of Example 2.
Example 3 (Lyapunov inequalities and state-dependent importance sampling) We shall apply importancesampling following the analysis in Section 4.2. We distinguish two situations; first, if the current positions is sufficiently far away from the boundary n, say if n − s > ∆ for some ∆ > 0, then the intuition
developed in Section 4.2 is accurate and then we apply importance sampling using the mixture family.
However, in the second situation, if s is sufficiently close to n we do not need to apply importancesampling. Also, if we use the previous family when s is sufficiently close to n, we might create variancebecause the mixture sampling distribution is derived using the large deviations regime, so turning off thesampling when s is sufficiently close to n is a sensible defense mechanism. We introduce a boundary layer
0 ≤ n − s ≤ ∆ for some ∆ > 0 and if s lies in this layer no importance sampling is applied.
Our family, then, takes the form
P (Xk+1 ∈ dy) I(y > a(n − s))p(n − s)I(n − s > ∆)Pe(Xk+1 ∈ dy|Sk = s) =P (Xk+1 > a(n − s))
P (Xk+1 ∈ dy) I(y ≤ a(n − s))[1 − p(n − s)]I(n − s > ∆)+P (Xk+1 ≤ a(n − s))+P (Xk+1 ∈ dy) I(n − s ≤ ∆) .
The parameter a ∈ (0, 1) is selected in order to enforce uniform integrability conditions, as we shall see.
Following the intuition developed earlier we proposeZ ∞p(n − s) = θP (X > a(n − s))/P (X > t)dt,
n−s

for some θ > 0. In fact, we expect θ ≈ µ to be optimal
R ∞ in the sense of variance minimization.
For notational convenience define H(n − s) = n−s P (X > t)dt. We now need to construct a valid
Lyapunov function. We wish to prove strong efficiency, so it is natural to suggest
G(n − s) = min(κH 2 (n − s), 1) = O(αn (s)2 )for some κ > 0. The strategy is then to select (given a ∈ (0, 1)) θ, κ and ∆ in order to satisfy the Lyapunov
inequality. If we can show that there exists a selection of these parameters that satisfies (18), then, giventhat the corresponding boundary condition in Lemma 2, namely
Z ∞P (X > t)dt, 1) = min(κEX + , 1) > 0
G(n − s) ≥ ρ := min(κ0is satisfied for s > n, we would conclude strong efficiency of the estimator.
Note that the inequality (18) holds trivially if G(n − s) = 1 is equivalent to n − s ≤ ∆. Indeed, in thiscase, since G (x) ≤ 1 for every x ∈ R, inequality (18) takes the form EG(n − s − X) ≤ G(n − s) = 1.
Now observe that G(n − s) < 1 holds if and only if H(n − s) < 1/κ1/2 , which in turn holds if andonly if n − s > H −1 (1/κ1/2 ). Therefore we can simply choose ∆ = H −1 (1/κ1/2 ) and simply select θand κ. In other words, we shall apply importance sampling using the mixture family if and only if ourcurrent position s is such that G(n − s) < 1. In order to proceed with the verification of inequality (18)
we henceforth assume that G(n − s) < 1. In this case inequality (18) takes the form
J1 + J2 ≤ 1,

156

(19)

Blanchet and Lam
where
E[G(n − s − X); X > a(n − s)]P (X > a(n − s)),κH 2 (n − s)p(n − s)E[G(n − s − X); X ≤ a(n − s)]P (X ≤ a(n − s))J2 =.
κH 2 (n − s)[1 − p(n − s)]J1 =We first bound the term J1 . Observe, noting that G(·) ≤ 1, and using the form of p(n − s), we obtainthat
P (X > a(n − s))p(n − s)J1 ≤=.
θκH(n − s)
θ2 κThen we analyze the term J2 . Note thatJ2 ≤E[H 2 (n − s − X); X ≤ a(n − s)].
H 2 (n − s)[1 − p(n − s)]Using the a Taylor expansion with remainder, expressed using elementary calculus with t = x + yu, we
obtainZx+y

H 2 (x + y) − H 2 (x) =2H(t)Ḣ (t) dt = 2E[yH(x + yU )Ḣ(x + yU )],x

where U is uniformly distributed in [0, 1]. Therefore, we have, letting x = (n − s) and y = −X,
E[H 2 (n − s − X)I(X ≤ a(n − s))]H 2 (n − s)
= P (X ≤ a(n − s)) − 2

E[X Ḣ(n − s − XU )H(n − s − XU ); X ≤ a(n − s)],
H 2 (n − s)where X and U are independent. Now,
E[X Ḣ(n − s − XU )H(n − s − XU ); X ≤ a(n − s)]
H 2 (n − s)
E[XF (n − s − XU )H(n − s − XU ); X ≤ a(n − s)]=.
H 2 (n − s)−We are assuming that G(n − s) < 1, or equivalently, that n − s > H −1 (1/κ1/2 ) = ∆. On the other hand,
we have that
F (n − s − XU )H(n − s − XU )I(X ≤ a(n − s))α−1X∼X
H 2 (n − s)(n − s)almost surely as n − s % ∞ and also we have, because of regular variation and Karamata’s theorem (seefor example Resnick 2004 p. 17) that there exists a constant K ∈ (0, ∞) such that
F (n − s − XU )H(n − s − XU )I(X ≤ a(n − s))
H 2 (n − s)
F ((n − s)(1 − a))H((n − s) (1 − a))
K |X|
≤ |X|
≤
.
2
H (n − s)
n−s

|X|

157

Blanchet and Lam
Consequently, if κ (or equivalently ∆ = H −1 (1/κ1/2 )) is chosen sufficiently large, we conclude, applying
the dominated convergence theorem (note here that a ∈ (0, 1) has been crucial) that for each ε > 0J2 ≤ 1 + 2EX(α − 1)(1 − ε)n−s

as long as G(n − s) < 1. Combining our estimates for J1 and J2 and choosing κ sufficiently large so thatp(n − s) ≤ θa−α

(α − 1)(1 + ε)1
(α − 1)(1 − ε),
≥ 1 + θa−α,n−s
1 − p(n − s)n−s

we arrive at(α − 1)(α − 1)(1 − ε)(α − 1)(1 + ε) −αa + (1 − 2µ
(1 − ε))(1 + θa−α)θκ(n − s)n−sn−s(α − 1)(1 + ε)
(α − 1)(α − 1)(1 − ε)≤ a−α+ (1 − 2µ(1 − ε) + θa−α).
θκ(n − s)n−sn−sJ1 + J2 ≤In order to enforce (19) we then need to ensure thata−α1+ε
− 2µ(1 − ε) + a−α θ(1 − ε) < 0,θκwhich clearly feasible given that µ > 0. Note that it is crucial that −µ = EX < 0, otherwise it would be
impossible to select the parameters appropriately. This makes sense given that when EX ≥ 0, αn (y) = 1
and it would be therefore impossible to bound the second moment of the estimator by a function that
decreases to zero as n → ∞, as we are proposing. See Blanchet and Liu (2010a) for further discussion on
the associated running time of the sampler.
ACKNOWLEDGEMENT
Support from the NSF foundation through the grants DMS-0806145, DMS-0846816 and DMS-1069064
is gratefully acknowledged.
REFERENCES
Asmussen, S. 1985. “Conjugate processes and the simulation of ruin problems”. Stochastic Processes andtheir Applications 20:213–229.
Asmussen, S. 2000. Ruin Probabilities. London, UK: World Scientific.
Asmussen, S. 2003. Applied Probability and Queues. New York, NY, USA: Springer-Verlag.
Asmussen, S., and K. Binswanger. 1997. “Simulation of Ruin probabilities for Subexponential Claims”.
ASTIN Bulletin 27:297–318.
Asmussen, S., K. Binswanger, and B. Hojgaard. 2000. “Rare events simulation for heavy-tailed distributions”.
Bernoulli 6:303–322.
Asmussen, S., and P. Glynn. 2007. Stochastic Simulation: Algorithms and Analysis. New York, NY, USA:
Springer-Verlag.
Asmussen, S., and D. Kroese. 2006. “Improved algorithms for rare event simulation with heavy tails”.
Advances in Applied Probability 38:545–558.
Asmussen, S., and R. Rubinstein. 1995. “Steady state rare events simulation in queueing models and its
complexity properties”. In Advances in Queueing, edited by J. Dshalalow, 429–461. CRC Press.
Bassamboo, A., S. Juneja, and A. Zeevi. 2006. “On the inefficiency of state-independent importancesampling in the presence of heavy-tails”. Operations Research Letters 34:521–531.

158

Blanchet and Lam
Blanchet, J., and P. Glynn. 2008. “Efficient rare-event simulation for the maximum of a heavy-tailed random
walk.”. Ann. of Appl. Probab. 18:1351–1378.
Blanchet, J., P. Glynn, and J. C. Liu. 2007. “Fluid heuristics, Lyapunov bounds and efficient importancesampling for a heavy-tailed G/G/1 queue”. QUESTA 57:99–113.
Blanchet, J., and H. Lam. 2011. “State-dependent importance sampling for rare-event simulation: An
overview and recent advances”. under review in Surveys in Operations Research and Management
Sciences.
Blanchet, J., and J. C. Liu. 2008. “State-dependent importance sampling for regularly varying random
walks”. Adv. in Appl. Probab. 40:1104–1128.
Blanchet,
J.,
and J. C. Liu. 2010a. “Efficient simulation and conditional functional
limit
theoremsforruinousheavy-tailed
randomwalks”.
submitted,
http://www.columbia.edu/ jb2814/papers/ConcaveHazard9f.pdf .
Blanchet, J., and J. C. Liu. 2010b. “Importance sampling in ruin problems for multidimensional regularly
varying random walks”. Journal of Applied Probability 47:301–322.
Bucklew, J. 1990. Large Deviation Techniques in Decision, Simulation and Estimation. New York, NY,
USA: Wiley.
Chan, H. P., and T. L. Lai. 2010. “A sequential Monte Carlo approach to computing tail probabilities instochastic models”. Ann. of Appl. Prob., forthcoming.
Collamore, J. 2002. “Importance sampling techniques for the multidimensional ruin problem for general
Markov additive sequences of random vectors”. Ann. Appl. Probab. 12:382–421.
Dembo, A., and O. Zeitouni. 1998. Large deviations techniques and applications, 2nd edition. New York,NY, USA: Springer.
Dupuis, P., K. Leder, and H. Wang. 2006. “Importance sampling for sums of random variables with regularly
varying tails.”. ACM TOMACS 17.
Dupuis, P., K. Leder, and H. Wang. 2007. “Large deviations and importance sampling for a tandem networkwith slowdown.”. QUESTA 57:71–83.
Dupuis, P., A. Sezer, and H. Wang. 2007. “Dynamic importance sampling for queueing networks.”. Ann.
Appl. Probab. 17:1306–1346.
Dupuis, P., and H. Wang. 2004. “Importance sampling, large deviations, and differential games”. Stoch.
Stoch. Reports 76:481–508.
Dupuis, P., and H. Wang. 2007. “Subsolutions of an Isaacs equation and efficient schemes of importancesampling”. Mathematics of Operations Research 32:723–757.
Embrechts, P., C. Klppelberg, and T. Mikosch. 1997. Modelling Extremal Events for Insurance and Finance.
New York: Springer-Verlag.
Embrechts, P., and N. Veraverbeke. 1982. “Estimates for the probability of ruin with special emphasis on
the possibility of large claims”. Insurance: Math. Economics 1:55–72.
Glasserman, P., and S. Kou. 1995. “Analysis of an importance sampling estimator for tandem queues”.
ACM Trans. Mod. Comp. Sim. 5:22–42.
Glasserman, P., and Y. Wang. 1997. “Counter examples in importance sampling for large deviations
probabilities”. Annals of Applied Probability 7:731–746.
Juneja, S., and P. Shahabuddin. 2002. “Simulating heavy tailed processes using delayed hazard rate twisting”.
ACM Trans. Model. Comput. Simul. 12 (2): 94–118.
Juneja, S., and P. Shahabuddin. 2006. “Rare event simulation techniques: An introduction and recent
advances”. Handbook on Simulation:291–350.
Juneja, S., P. Shahabuddin, and A. Chandra. 1999, December. “Simulating heavy tailed processes using
delayed hazard rate twisting”. In Proceedings of the 1999 Winter Simulation Conference, edited by
P. A. Farrington, H. B. Nembhard, D. T. Sturrock, and G. Evans, 420–427. Piscataway, New Jersey:
Institute of Electrical and Electronics Engineers, Inc.

159

Blanchet and Lam
Kollman, C., K. Baggerly, D. Cox, and R. Picard. 1999. “Adaptive importance sampling on discrete Markov
chains”. Annals of Applied Probability 9:391–412.
L’Ecuyer, P., J. Blanchet, B. Tuffin, and P. Glynn. 2010. “Asymptotic robustness of estimators in rare-event
simulation”. ACM TOMACS 20:6:1–6:41.
L’Ecuyer, P., F. L. Gland, P. Lezaud, and B. Tuffin. 2009. “Splitting techniques”. In Rare Event Simulation
Using Monte Carlo Methods, edited by G. Rubino and B. Tuffin, 39–62. West Sussex, United Kingdom:
Wiley. Chapter 3.
Pakes, A. G. 1975. “On the tails of waiting-time distribution”. J. Appl. Probab. 12:555–564.
Parekh, S., and J. Walrand. 1989. “Quick simulation of rare events in networks”. IEEE Trans. Automat.
Contr. 34:54–66.
Resnick, S. I. 2004. Extreme Values, Regular Variation, and Point Processes. New York, USA: Springer.
Rubinstein, R., and D. Kroese. 2004. The Cross-Entropy Method. New York, NY: Springer.
Rubinstein, R. Y., and D. P. Kroese. 2009. Simulation and the Monte Carlo Method, 2nd Edition. John
Wiley.
Sadowsky, J. S. 1991. “Large deviations theory and efficient simulation of excessive backlogs in a GI/GI/mqueue”. IEEE Trans. Autom. Control 36:1383–1394.
Siegmund, D. 1976. “Importance sampling in the Monte Carlo study of sequential tests”. Ann. Stat. 4:673–
684.
Veraverbeke, N. 1977. “Asymptotic behavior of Wiener-Hopf factors of a random walk”. Stochastic Process.
Appl. 5:27–31.
AUTHOR BIOGRAPHIESJOSE BLANCHET is a faculty member of the IEOR Department at Columbia University. Jose holds a
Ph.D. in Management Science and Engineering from Stanford University. Prior to joining Columbia he wasa faculty member in the Statistics Department at Harvard University. Jose is a recipient of the 2009 Best
Publication Award given by the INFORMS Applied Probability Society and of the 2010 Erlang Prize. He
also received a PECASE award sponsored by NSF in 2010. He worked as an analyst in Protego Financial
Advisors, a leading investment bank in Mexico. He has research interests in applied probability and Monte
Carlo methods. He serves in the editorial board of Advances in Applied Probability, Journal of Applied
Probability, Mathematics of Operations Research and QUESTA.
HENRY LAM is an Assistant Professor in the Department of Mathematics and Statistics at Boston
University. He recently graduated from Harvard University with a Ph.D. degree in statistics. His researchinterests lie in applied probability and Monte Carlo methods with applications in queueing, operations
management and insurance modeling.

160