Surveys in Operations Research and Management Science 17 (2012) 38–59

Contents lists available at SciVerse ScienceDirect

Surveys in Operations Research and Management Sciencejournal homepage: www.elsevier.com/locate/sorms

Review

State-dependent importance sampling for rare-event simulation: An overview
and recent advancesJose Blanchet a,∗ , Henry Lam baColumbia University, United Statesb

Boston University, United Statesarticleinfo

Article history:
Received 7 February 2011
Received in revised form
19 September 2011
Accepted 20 September 2011

abstract
This paper surveys recent techniques that have been developed for rare-event analysis of stochastic
systems via simulation. We review standard (state-independent) techniques that take advantage of large
deviations results for the design of efficient importance sampling estimators. Classical examples andcounter-examples are discussed to illustrate the reach and limitations of the state-independent approach.
Then we move to state-dependent techniques. These techniques can be applied to both light and heavytailed systems and are based on subsolutions (see e.g. Dupuis and Wang (2004) [5], Dupuis and Wang(2007) [6], Dupuis and Wang (2009) [80], Dupuis et al. (2007) [7]) and Lyapunov bounds (Blanchet and
Glynn (2008) [9], Blanchet et al. (2007) [11], Blanchet (2009) [12]). We briefly review the ideas behindthese techniques, and provide several examples in which they are applicable.
© 2011 Published by Elsevier B.V.

Contents
1.
2.

3.

4.

5.

Introduction........................................................................................................................................................................................................................ 38
Notions on efficiency and importance sampling techniques .......................................................................................................................................... 40
2.1.
Notions of efficiency .............................................................................................................................................................................................. 40
2.2.
Basic properties of importance sampling estimators .......................................................................................................................................... 41
On large deviations techniques and importance sampling ............................................................................................................................................. 41
3.1.
Large deviations results for random walks .......................................................................................................................................................... 41
3.2.
Applications of large deviations to rare-event simulation.................................................................................................................................. 42
3.3.
Notes on state-independent importance sampling and large deviations .......................................................................................................... 44
State-dependent techniques for light-tailed systems ..................................................................................................................................................... 44
4.1.
Solutions to a deterministic control problem and efficient importance sampling ........................................................................................... 44
4.2.
Subsolutions to a deterministic control problem and efficient importance sampling ..................................................................................... 46
4.3.
Notes on state-dependent importance sampling for light-tailed systems ........................................................................................................ 50
Large deviations and importance sampling for heavy-tailed random walks ................................................................................................................. 50
5.1.
Classical simulation techniques for heavy-tailed random walks ....................................................................................................................... 50
5.2.
Large deviations for heavy-tailed random walks ................................................................................................................................................ 52
5.3.
Conditional distributions and implications for simulation ................................................................................................................................. 53
5.4.
State-dependent importance sampling techniques for heavy-tailed random walks........................................................................................ 54
5.5.
Notes on importance sampling with heavy tails ................................................................................................................................................. 57
Acknowledgments ............................................................................................................................................................................................................. 57
References........................................................................................................................................................................................................................... 57

1. Introduction
In recent years, there has been a substantial amount of research
related to the efficient design of state-dependent importance∗ Corresponding author.
E-mail address: jose.blanchet@gmail.com (J. Blanchet).
1876-7354/$ – see front matter © 2011 Published by Elsevier B.V.
doi:10.1016/j.sorms.2011.09.002

sampling estimators. This research has been motivated by a seriesof examples and counter-examples (see [1,2]) related to the useof large deviations principles in the design of efficient importancesampling estimators for light-tailed systems. The development ofefficient importance sampling estimators for heavy-tailed systems,
which as explained in [3] must often be state-dependent, hasalso driven the research developments that are the focus of thissurvey (the paper [4] discusses additional challenges that arise

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

in the context of efficient rare-event simulation for heavy-tailed
systems).
In order to systematically address the construction of statedependent importance sampling estimators that can be rigorously
shown to be efficient in large deviations environments for lighttailed systems, Dupuis, Wang, and their students have developed
a method based on control theory and the use of subsolutionsof an associated Isaacs equation for constructing and testing
asymptotically optimal importance sampling estimators [5–8].
A related approach, based on the construction of Lyapunov
inequalities, has been used in [9–11] for the construction andanalysis of state-dependent importance sampling estimators forheavy-tailed systems. The Lyapunov method has also been used forcounting bipartite graphs with a given degree sequence in [12].
In this paper we summarize the ideas behind the developmentof state-dependent importance sampling estimators and illustrate
their applications in some simple examples. However, in order to
put these recent ideas in context we will review some of the basic
notions and techniques that motivated the development of statedependent samplers and their efficiency analysis.
We shall start by discussing standard notions of statistical
efficiency that are used in order to rigorously test the optimalityof a simulation estimator. This will be done in Section 2. As
we shall see, when the cost per replication of the estimator isappropriately controlled, these notions guarantee a substantial
improvement, in terms of running time, relative to crude
Monte Carlo evaluation—which typically offers a natural generalbenchmark. The introduction of these notions of efficiency will
naturally lead to the question of how to actually construct
importance sampling estimators that can be rigorously shown tobe 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. It is very
important to emphasize, as we shall see, that importance samplingallows to efficiently estimate conditional expectations given theevent of interest.
The zero-variance change-of-measure turns out to be theconditional 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 departing point in
the design of efficient importance sampling estimators for rare
events, especially in the setting of light-tailed systems, which
will be our main focus in the first part of the paper and which
encompasses most classical models such as Jackson networks.
Regardless of whether one is studying systems with light or heavy
tails, it is extremely important, however, to note that once anefficient importance sampling estimator is in place, then one
obtains a computational mechanism (via Monte Carlo) that allows
to enhance a given asymptotic result by providing an estimate(obtained by averaging iid replications) whose accuracy can be
improved to any desired level of relative precision by increasing
the sample size. In turn, the guaranteed efficiency will typicallyensure that the number of replications will remain under control
as the event of interest becomes more and more rare according to
the underlying large deviations parameter.
In Section 3 we will review some of the earlier ideas in theliterature showing the basic connection between large deviations
theory and the design of efficient importance sampling estimators.
As an illustrative example, we will study the sample-path large
deviations of random walk, how it suggests a conditional path andhow it leads to an efficient state-independent exponential tilting
scheme for a first passage problem.
The successful application of large deviations theory inthe design of efficient importance sampling estimators is notstraightforward and this feature turned out to be the norm

39

rather than the exception, especially in Operations Research
applications such as the analysis of queueing networks. This might
be somewhat surprising given the fact that large deviations theory
allows to characterize the asymptotic conditional distribution (i.e.
the zero-variance change-of-measure) leading to the occurrenceof the rare event. The problem, as we shall see, arises from the factthat such characterization is done in fluid scale (i.e. in the spatial
and temporal scales corresponding to the Law of Large Numbers)
and there are many equivalent changes-of-measure that possess
the same Law of Large Numbers limit and that coincide with thatof the zero-variance change-of-measure. In Section 3 we will firstillustrate with an example the successful application of direct large
deviations results and then study in detail counter-examples thatflesh-out precisely the issues that can go wrong in the variance
control of importance sampling estimators.
Our discussion in Section 3 then will lead us to introduce
the approach proposed by Dupuis and Wang, which formulates
the design of state-dependent importance sampling estimators interms 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 4, the stochastic control problem becomesa deterministic control problem. The simplification comes fromthe fact that only a subsolution is needed for the design of an
efficient sampling scheme. A trivial supersolution is given, thanks
to Jensen’s inequality, by twice the underlying large deviations
rate function. The whole point of Dupuis and Wang’s approach is
to construct a subsolution that matches the trivial supersolution
at one point, namely, the initial condition of the system underconsideration. In many problems of interest it turns out that sucha subsolution is piecewise affine and therefore relatively easy tomanipulate.
We will then move to the analysis of heavy-tailed systems,
which will be introduced in Section 5. These types of systems
arise in the analysis of finance and insurance models. There arefundamental qualitative differences between light and heavytailed systems when it comes to the analysis of rare 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 fromheavy-tailed building blocks. Light-tailed large deviations occur
gradually, while heavy-tailed large deviations occur suddenly andthis feature manifests itself in the design of efficient importancesampling estimators for heavy-tailed systems. This difference isillustrated by, for example, the problem of estimating the firstpassage probability for a random walk i.e. P (τn < ∞), whereτn = inf{k ≥ 0: Sk > n}, Sk = X1 + · · · + Xk and Xi ’s are iid
increments with negative mean. If Xi has exponential moments i.e.
E exp(θ Xi ) < ∞ for any θ ∈ R, then the most likely way for the
first passage to occur is by small gradual contribution from eachincrement to drift upwards. This translates into an exponential
tilting for each increment until τn is hit. On the other hand, if
Xi is regularly varying (i.e. basically having a tail distributionfunction with power-law decay) then the large deviations event
occurs typically by a big jump in one of the increments, whereas
all the other increments behave like the nominal (unconditional)
distribution (i.e. with the same unconditional negative meanunder fluid scale). One efficient importance sampling scheme inthis case will be a mixture algorithm that recognizes a suitablychosen probability for a sudden jump to occur. This first passage
problem will be visited in Example 1 (for the light-tailed case) and
Example 5 (for the regularly varying case and beyond).
Throughout the paper we use Landau’s notation for theasymptotic behavior of functions. In particular, given two nonnegative functions (f (N ): N ≥ 1) and (g (N ): N ≥ 1) we say thatf (N ) = O (g (N )) if there exists a constant c ∈ (0, ∞) independentof N such that f (N ) ≤ cg (N ) for N ≥ 1. We write f (N ) = Ω (g (N ))if f (N ) ≥ c ′ g (N ) for some constant c ′ ∈ (0, ∞) independent

40

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

of N. We write f (N ) = Θ (g (N )) if both f (N ) = O (g (N )) andf (N ) = Ω (g (N )). Finally we say that f (N ) = o (g (N )) as N → ∞if f (N )/g (N ) → 0 as N → ∞.
We emphasize that this paper aims to introduce readers the
key concepts and ideas in rare-event problems. Throughout the
exposition we may drop some mathematical details and highlight
only the central arguments, pointing readers to other references.
The rest of the paper is organized as follows. In Section 2 weintroduce the basic notions, criteria and guidelines for rare-event
simulation estimators. In Section 3, we discuss the classical large
deviations theory for light-tailed systems from a sample pathperspective, and how it leads to the idea of exponential tilting
in importance sampling. There we will also point out someissues of this classical idea through examples. Section 4 willthen be devoted to the recent developments of state-dependent
techniques that aim to fix these issues, using a subsolution
approach as discussed briefly above. Lastly, Section 5 focuses onheavy-tailed systems, including large deviations results, classicalestimators using conditional Monte Carlo and hazard rate twisting,and the recent work on state-dependent importance samplers via
weak convergence and Lyapunov technique.

1/2

squared error equals Var(αnc (N ))1/2 = αn (1 − αn )1/2 /N 1/2 .
In turn, we obtain a relative mean squared error equal to thecoefficient of variation of 
αnc (N ), namely,CV (αnc (N )) :=The central question that we study consists in estimating via
simulation P (A), where A is some rare event of interest, which
implies that P (A) ≈ 0. For instance, the event A might correspondto the probability of bankruptcy within a given time horizon ofa financial or insurance company. A related question consistsin computing conditional expectations of certain quantities ofinterest given that A occurs. For instance, what is the mean valueof the deficit at the time of ruin or, in the case of companies thatpay dividends to their shareholders, what is the distribution of the
net present value of the dividends paid up to the time of ruin. In
this section we discuss notions that allow to quantify the efficiencyof rare event simulation estimators and review basic properties ofimportance sampling estimators for rare events.

as n ↗ ∞.

Nαnc (N ) =1 −N j =1In (j).

(1)

The super-index ‘‘c’’ that appears in 
αnc (N ) corresponds to crudecMonte Carlo. Since 
αn (N ) is an unbiased estimator of αn , its mean.

Definition 1. An estimator Rn is said to be strongly efficient if
ER2n = O αn2

In order to design algorithms that achieve high accuracy inrelative terms and that are designed to perform well in rare
event settings, a typical approach is to embed the question ofinterest in a suitable asymptotic regime, often motivated by anapplied standpoint. For example, in the setting of bankruptcy of an
insurance company, it is often the case that the capital or reserveof the company is large relative to individual claim sizes, and thisis basically why one might expect the probability of bankruptcy,
P (A), to be small. So, one might introduce a parameter, say n,
for instance as the initial capital in the insurance example, andconsider αn = P (An ) as a function of the parameter n which isoften referred to as the ‘‘rarity parameter’’. From now on we willassume that underlying the rare event there is a rarity parametern 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 crude(or naive) Monte Carlo, which provides an obvious alternative
and, therefore, a general benchmark. Note that crude Monte Carlo
involves simulating N iid Bernoulli’s, In (1), . . . , In (N ), with successparameter αn , thereby producing an estimator of the form

1/2

αn N 1/2

In order to control the relative error of 
αnc (N ) (i.e. CV (α̂nc (N ))remains bounded as n ↗ ∞), one needs N = Ω (1/αn ) Bernoulli’s.
The problem with the estimator 
αnc (N ) is that each of theunderlying replications, the In (j)’s, has a very big variance inrelative terms. The number of replications N must grow at the samerate as CV (In (j))2 = (1 − αn )/αn . The same phenomenon occurs inany estimator obtained out of averaging iid replications, as in (1).
The objective is then to design an estimator Rn with a controlledmean 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 . Motivatedby these considerations we have the following definitions.

2. Notions on efficiency and importance sampling techniques

2.1. Notions of efficiency(1 − αn )1/2





Definition 2. Rn is said to be weakly efficient or asymptotically
optimal if for each ε > 0 we have that
ER2n = O αn2−ε

as n ↗ ∞.
Definition 3. We say that Rn has polynomial complexity of order at
most l ≥ 0 if
ER2n = O αn2 (log (1/αn ))2l ,





(2)

or in other words, if its relative mean squared error (or coefficientof variation) grows at most at a polynomial rate with degree l inlog (1/αn ). The estimator is strongly efficient if l = 0 in (2).
In most cases the analysis of importance sampling estimators(especially state-dependent samplers in light-tailed cases) concludes only weak efficiency of the estimators. Another class of importance sampling changes-of-measure that have been proposed
recently exhibits asymptotically negligible relative mean squared
error, defined next (see [13]):
Definition 4. The estimator Rn exhibits asymptotically negligiblerelative error if given any ε > 0 one has that
lim

n→∞

ER2n

αn2

≤ 1 + ε,or in other words, if its relative mean squared error (or coefficientof variation) is less than ε .
In order to see the gain in efficiency, suppose that an estimator
Rn has polynomial complexity of order at most l and consider∑Nα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 log (1/αn )2l ] replications are required toconclude that 
αn (N ) is ε -close to αn in relative terms with atleast 1 − δ confidence (i.e. P (|α̃n (N ) − αn | > εαn ) ≤ δ ).
Observe that in some sense, assuming that the cost per replicationof Rn (j) is comparable to that of In (j) (in the setting of crude Monte
Carlo), polynomially efficient estimators provide an exponential gainin efficiency relative to crude Monte Carlo.The paper of [14] provides

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

further study and examples related to each of the aforementioned
classes of estimators.
Besides variance control, it is worth pointing out othermeasures for importance sampling estimators, and we referinterested readers to the following references. One criterion
is the robustness of higher moment estimates as the rarity
parameter goes to zero, as discussed in [15]. In the special caseof second moment, [16,17] introduce the notion of bounded
normal approximation, suggested by Berry–Esseen theorem,in the context of highly reliable Markovian systems. Related
diagnostic robustness check of confidence interval coverage isdiscussed in [18] with so-called coverage functions. Secondly,
recall that our criteria in this section all presume that the costsof replication are comparable between importance sampler andcrude Monte Carlo. This might not always be true, i.e. the cost ofgenerating each sample from importance sampling may increase
disproportionately with the rarity parameter. The paper [19]discusses this issue and suggests the use of so-called relative
efficiency. The paper [20] justifies theoretically how to include
computational effort per sample rightly into efficiency measures.
Finally, [21] provides a survey of the above robustness issues inrare-event estimators.
2.2. Basic properties of importance sampling estimators
Importance sampling is a variance reduction technique thatis often applied to design an estimator, Rn , with efficiency
characteristics such as those described by Definitions 1–4. We
shall review basic properties of importance sampling next; for amore detailed discussion see for example [22–24]. The basic idea
is to introduce a probability measure P (·) such that the likelihood
ratio or Radon–Nikodym derivative between the nominal (original)
probability measure, P (·), and P (·) is well defined on the event Anof interest. Then we can simply letdP

Rn (ω) =dP(ω) I (ω ∈ An ) .

∫

ERn =

dP
An d PE( ) =

dP

An

d
P

∫
=

dP
An d P(ω)

1

log P (Yn (·) ∈ C )≤ − inf{I (x(·)) : x (·) ∈ AC [0, T ] ∩ C , x (0) = 0},

and for any open set O
lim

1

log P (Yn (·) ∈ O)

∫ T

J (ẋ(s)) ds

and
(3)


P (·) = Pn∗ (·) := P (· | An ) ,
then
dPn∗

Theorem 1. Consider the continuous approximation Yn (t ) = S⌊nt ⌋/n +(t −⌊nt ⌋ /n)X⌊nt ⌋+1 obtained by joining the values of the random
walk at the lattice points and by scaling space by 1/n and time by n.
Consider the space C [0, T ] for any T > 0, with the topology generated
by the uniform norm and let us write AC [0, T ] to denote the space ofabsolutely continuous function in [0, T ]. Then, for any closed set C

0

In principle one can design an importance sampling estimator withzero variance. Indeed, if we letdP

We shall take advantage of the following fundamental result inthe theory of large deviations, known as Mogulskii’s theorem (seefor instance [25, p. 176 Theorem 5.1.2]).

I (x(·)) =R∗n :=

3.1. Large deviations results for random walks

where

d
P (ω)

(ω) dP (ω) = E (Rn ).

where the Xi ’s are iid random variables. We assume that ψ (θ ) =log E exp (θ Xi ) is finite in a neighborhood of the origin over thereal axis. We use ψ̇ (·) to denote the derivative of ψ (·). Given any
parameter a ∈ R, we will also assume that a solution θa to a root
equation such as ψ̇ (θa ) = a exists whenever we require such asolution. The role of the parameter a will depend on the context ofthe examples and applications that we shall discuss.

≥ − inf{I (x(·)) : x (·) ∈ AC [0, T ] ∩ O, x (0) = 0},

so Rn is an unbiased estimator, andR2nS k = X1 + · · · + Xk ,n→∞ n

P (ω) = P (An ) = αn ,(ω) d

2

In order to explain the interplay between large deviations
techniques and the design of importance sampling estimators
we shall concentrate on models whose fundamental building
blocks are given by random walks. Many questions of interest
in Operations Research can be posed in terms of sample-path
properties of random walks; we will see some examples in thesequel. In order to guide our discussion we shall concentrate on
several specific instances underlying a generic random walk S =(Sk : k ≥ 1), defined via

lim



∫ 

3. On large deviations techniques and importance samplingn→∞ n

As usual, we use ω to denote the underlying random outcome
which is simulated according to the probability measure P (·).
We use E (·) to denote the expectation operator associated to theprobability measure P (·). Note that41(ω) I (ω ∈ An ) = P (An ) ,which is clearly an exact estimate of αn = P (An ). Of course,
assuming that one has access to P ∗ (·) defeats the point ofusing simulation to estimate P (An ). However, the underlying
lesson behind the characterization of the zero-variance changeof-measure as the conditional distribution given the rare event ofinterest is that one can use asymptotic theory to describe Pn∗ (·)
as n ↗ ∞. As we shall review in the next section, this is often
the starting point in the design of efficient importance sampling
estimators.

J (z ) = sup[θ z − ψ (θ)].
θ∈R

Remark 1. Given a set A we often write
I (A) = inf{I (x(·)) : x (·) ∈ AC [0, T ] ∩ A, x (0) = 0}.
The previous result provides rigorous support for the heuristic
approximation P (Yn (·) ≈ y) ≈ exp (−nI (y)). In most applications
in Operations Research, and in particular in Examples 1 and 2 to be
discussed below, we are interested in a set A such that
P (Yn (·) ∈ A) = exp (−nI (A) + o(n)) .

(4)

That is, sets for which the limits in the previous theorem coincide
when applied to the interior and the closure of A.
Now, the solution x∗ (·) to the calculus of variations problem
defining I (A) is often unique and describes the law of large
numbers of Yn (·) given that Yn (·) ∈ A and in this sense we

42

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

say that x∗ (·) approximates the conditional distribution of Yn (·)
given Yn (·) ∈ A. This is easy to see using the earlier heuristic
approximation (a rigorous verification is also almost direct and canbe found in [26, p. 61 Lemma 4.2]). Indeed, note that if x∗ ̸= y ∈ A,
then
P (Yn (·) ≈ y|Yn (·) ∈ A) ≈ exp (−nI (y) + nI (A))

 
= exp −nI (y) + nI x∗ −→ 0
as n → ∞, since x∗ is the unique optimizer of I (A).
Given that x∗ approximates the conditional distribution ofYn (·) given the event of interest it makes sense to bias thepaths to track x∗ . In order to do this one uses a very particularimportance sampling mechanism called exponential tilting, which
is implemented increment-by-increment in the simulation viaPθk (Xk ∈ dx) = exp(θk x − ψ(θk ))P (Xk ∈ dx).

(5)

If Xk has a density with respect to the Lebesgue measure, sayf (·), then P (Xk ∈ dx) = f (x) dx. The notation adopted in theprevious display is simply a mechanism that allows us to consider
general distributions (including discrete or mixtures of discreteand continuous distributions). We will use Eθk (·) to denote theexpectation operator associated to Pθk (·). The parameter θk might
depend on (X1 , . . . , Xk−1 ) but it is not allowed to depend on future
observations. For example, if under P (·), Xk is standard Gaussian,then under Pθk (·), Xk is normal with mean θk and unit variance.
On the other hand, if under P (·), Xk is exponentially distributedwith unit mean, then under Pθk (·) (for θk < 1), Xk is exponentiallydistributed with mean 1/(1 − θk ).
In order to use exponential tilting to track the optimal path x∗ (·)
it is natural to select for the k-th increment a value θk such thatEθk Xk = ẋ∗ (k/n).

(6)

Since

= ψ(η + θk ) − ψ(θk ),we conclude that Eθk Xk = ψ̇(θk ) and therefore Eq. (6) is equivalent
to(7)Note that selecting θk via (7) gives rise to a state-independent
importance sampling estimator because the simulation of the k-th
increment does not depend on the position Sk of the random walk.
A natural question that might arise in the reader’s mind is why
shall we use exponential tilting as a biasing mechanism? One could
think, for instance, to simply shift the means of the increments.
This, as we saw earlier, is equivalent to exponential tilting in theGaussian case, but it is not so in the case of exponential randomvariables. Moreover, in the later case, the change-of-measureinduced by shifting might not be admissible because one could
violate absolute continuity (the supports of the distributions might
not even match) and the likelihood ratio might not be even welldefined. The fundamental reason behind the use of exponential
tilting arises from entropy considerations. Recall that the entropy
between two probability measures P (·) and P (·) is defined via  dPH (P ‖ P ) = E log.
dP

3.2. Applications of large deviations to rare-event simulation
We are now ready to apply the approach discussed in theprevious subsection to a couple of examples. Our goal is toinvestigate the extent to which the tracking of the optimal path
via state-independent exponential tilting allows us to obtaina provably efficient importance sampling estimator. Example 1
illustrates a situation in which weak efficiency is obtained and, as
we shall see, Example 2 shows an instance in which no efficiency
is guaranteed.
The first example is classical in insurance and is of significant
historical importance in rare-event simulation because it wasamong the first problems for which an explicit connection to large
deviations theory was used in order to construct and analyze
efficient importance sampling estimators (see [29]).
Example 1 (Using Large Deviations For Efficient Importance Sampling). Suppose that S0 = 0 and let τn = inf{k ≥ 0: Sk > n}.
Assume that EXk < 0 and consider the problem of computing

αn = P (τn < ∞)as n ↗ ∞. In order to motivate this example note that the surplus
process of a insurance company at time t > 0 in continuous time
can be represented via the process

log Eθk exp(ηXk ) = log E exp((η + θk )Xk − ψ(θk ))ψ̇(θk ) = ẋ∗ (k/n).

is precisely given by 
P (·) = Pθk (·). A proof of this fact can be foundin [27] p. 409 and 421. The reader might wish to verify this fact atleast in the case of distributions with finite support, in which caseis an elementary exercise in convex optimization.
Another reason that motivates the use of exponentially tilted
distributions is that one can sometimes show that the distributionof (X1 , . . . , Xk ) given a large deviations event of interest converges
to the optimal exponential tilting (i.e. selecting each tilting
parameter according to (7)); see for instance [23] Chapter VI
Section 5 and [28].

(8)

Y (t ) = n + pt −N (t )
−

Zj ,

(9)

j =1

where N (·) represents a renewal process that counts the numberof arrivals up to time t, p > 0 is a premium rate paid continuously
in time, n := Y (0) is the initial surplus, and the Zj ’s represent iid
claim sizes—which are assumed to be positive and independent of
the process N (·). Eventual ruin, namely the event that Y (·) hits
a level below zero at some point into the future, can only occur at
arrival times. Therefore, if Ak is the arrival of the k-th claim into thesystem then eventual ruin occurs if and only if there exists k ≥ 1such that Y (Ak ) < 0. In particular, if we write Ak = T1 + · · · + Tk ,where the Tj ’s are iid inter-arrival times, then ruin happens if andonly if τn < ∞, where Xj = Zj − pTj .
We now compute the optimal path introduced in Theorem 1,
which involves solving the calculus of variation problem. Recall
that we are assuming EXk = ψ̇(0) < 0. We also will assume
that there exists θ ∗ > 0 such that ψ(θ ∗ ) = 0 and ψ̇(θ ∗ ) < ∞.
We are naturally assuming that Xk is non-deterministic (otherwise
αn = 0), therefore ψ (·) is strictly convex and θ ∗ is unique (see[30, p. 352]). Moreover, also by convexity we have that ψ̇(θ ∗ ) > 0.
We need to find x∗ (·) that solves the problem

∫ ∞



Note that H (P ‖ P ) ≥ 0 and H (P ‖ P ) = 0 if and only if
P = P, so H (·) serves as a measure of ‘‘discrepancy’’ betweenmin

two probability measures. Exponential tilting arises as a natural
biasing mechanism because, in connection to (6), the solution tothe optimization problem

Since I (z ) ≥ 0 it is evident that the optimal path must be such that
I (ẋ∗ (t )) = 0 for t > T ∗ , where T ∗ = inf{t ≥ 0: x∗ (t ) ≥ 1}. This
forces ẋ∗ (t ) = ψ̇ (0) for t > T ∗ . Therefore, we need to minimize

min{H (P ‖ P ): 
P (·) such that 
EXk = ẋ∗ (k/n)}

I (ẋ(s))ds: x(0) = 0, x(t ) ≥ 1 for some t > 0 .

0

T
0

I (ẋ(s))ds, where T = inf{t ≥ 0: x (t ) ≥ 1}. Now, the function

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

I (·) is convex because it is the Legendre transform of the convex
function ψ (·) (see for example, [31, Part III, Section 12]). So, by
convexity of I (·) and by Jensen’s inequality we have that
T

1

∫ T

T

0

I (ẋ(s))ds ≥ TI

 ∫ T
1

T

ẋ(s)ds

= TI (1/T ).

(10)

0

Optimizing the right hand side of (10) over T we then must selectT ∗ in order to satisfy the first order conditions of optimality,
namely,− T∗

dI (1/T ∗ )
dT ∗



= I 1/T ∗ .

(11)

Now we will use the previous equation to characterize theoptimal path. Note that
I (z ) = z θz − ψ(θz ),11

θ1/T ∗ = ∗ θ1/T ∗ − ψ(θ1/T ∗ ) ⇒ ψ(θ1/T ∗ ) = 0.
T∗
T
The previous equation implies that θ1/T ∗ = θ and therefore thesecalculations suggest that x∗ (t ) = t ψ̇(θ ∗ ) for t ≤ T ∗ = 1/ψ̇(θ ∗ ). To
complete the argument we must show that the lower bound in (10)is achieved by our selection, but this is clearly the case given that
our candidate path is linear up to T ∗ and therefore its derivative isconstant.
The corresponding importance sampling algorithm then consists in simulating each increment Xk according to the exponentially tilted distribution induced by the tilting parameter θ ∗ . Of
course, one stops the simulation until time τn which is finite withprobability one because under the importance sampling distribution Eθ ∗ Xk = ψ̇(θ∗ ) > 0. A single replication of the importancesampling estimator is then given by
∗

Rn :=

τn
∏
dP
k=1

dPθ ∗

Example 2 (Lack of Efficiency from Tracking the Optimal Path in aState-Independent Way). Given a < 0 < b and S0 = 0 we areinterested in efficient estimation via simulation of the probability

αn = P [ min Sk ≤ an, Sn > nb].

(12)

0≤k≤n

Such probability can be interpreted as the price of a digital knockin option in the absence of interest rates. That is, the price of an
option that upon exercise pays one unit at time n if and only if
both Sn > nb and, in addition, the random walk hits a level lessor equal than an before time n. (See [32] for more applications ofimportance sampling to option pricing.)
For concreteness, in order to facilitate the discussion, let us
assume that Xk is standard Gaussian. The calculus of variations
problem that must be solved takes the form

 ∫ 1

where ψ̇ (θz ) = z. This implies that dI (z )/dz = θz and therefore(11) is equivalent to

43

inf

1
2



ẋ (s) ds: x(0) = 0, min x (u) ≤ a, x (1) ≥ b .
2

0≤u≤1

0

(13)

The solution x∗ attaining the infimum in (13) can be computed
explicitly asx∗ (t ) = s1 tI (t ≤ t0 ) + [s2 (t − t0 ) − a]I (t0 < t ≤ 1) ,where
s1 = a/t0 ,s2 = (b − a) / (1 − t0 ) ,t0 = −a/(b − 2a).

Verifying that this indeed is the optimal path is done using Jensen’s
inequality and convexity arguments similar in spirit to those usedin Example 1. In turn, we obtain thatαn = exp(−n(b − 2a)2 /2 + o(n))(14)

as n → ∞ by evaluating I (x ). Alternatively, given that we areworking with Gaussian increments, we can replace the random
walk by Brownian motion B (·) and then, using the reflection
principle, we can conclude that∗P ( min B (t ) ≤ an, B(n) > bn) = P (B(n) ≤ − (b − 2a) n)(Xk )I (τn < ∞) = exp(−θ ∗ Sτn ).

We now proceed to verify weak efficiency. We first need toobtain the asymptotics for αn . One can use the previous optimality
argument for x∗ combined with Theorem 1 to conclude thatαn = exp(−nI (x∗ ) + o(n)) = exp(−nθ ∗ + o(n))
as n → ∞ (in particular this corresponds to the limit (4)). On theother hand, using the expression for the second moment given in(3) we conclude that
ERn = E [exp(−θ ∗ Sτn )I (τn < ∞)] ≤ exp(−nθ ∗ )P (τn < ∞)= exp(−2nθ ∗ + o(n)).
The previous expression therefore implies that if there exists θ ∗ >0 such that ψ(θ ∗ ) = 0 and ψ̇(θ ∗ ) ∈ (0, ∞), then Rn is a weakly
efficient estimator. Now, in fact, assuming only a mild regularly
condition, for instance if the Xk ’s have a continuous distribution,
it turns out that one can apply renewal theory to the incrementsof the strictly increasing ladder heights (i.e. the increments of thesequence of new records of the random walk) to conclude that Rn
is in fact a strongly efficient estimator (see [30, p. 375 Theorem7.1]). If we have that ψ̇(θ ∗ ) = ∞, then strong efficiency is violated
but one can obtain weak efficiency or polynomial complexity
depending on the tail of the distribution of Xk .
We now move on to study Example 2. We have selected thisexample because it illustrates that direct application of stateindependent importance sampling, guided by direct, yet sensible,
large deviations principles, might not yield efficient estimators.

0≤t ≤n= exp(−n(b − 2a)2 /2 + o(n))
as n → ∞, obtaining the asymptotic rate indicated earlier for αn .
Since we are assuming Gaussian increments, the tilting
parameter for the increment Xk , namely θk∗ , matches ẋ∗ (k/n) andthis yields

θk∗ = − (b − 2a) I (−na/(b − 2a) > k)(b − a) I (−na/(b − 2a) ≤ k)+
1 − (−a) /(b − 2a)= − (b − 2a) I (−na/(b − 2a) > k)+ (b − 2a) I (−na/(b − 2a) ≤ k) .
The importance sampling estimator then takes the form
Rn = exp(−[θ0∗ Sk0 − k0 θ0∗2 /2] − θn∗ (Sn − Sk0 ))× exp((n − k0 )θn∗2 /2)I ( min Sk ≤ an, Sn > bn),0≤k≤nwhere k0 = ⌊−an/ (b − 2a)⌋.
We now show that Rn might not be efficient. We write P (·)
to denote the probability measure induced by sampling the Xk ’s
independently according to the exponential tilting parameter θk∗ .
First we develop a lower bound for 
E [R2n ] by constructing a suitable
set under which the path deviates a distance of size ε > 0 (insupremum norm) from the optimal trajectory. In order to simplify
the notation let θ ∗ = (b − 2a) and then note thatR2n = exp(−nθ ∗ ) exp 2θ ∗ Sk0 − Sn − Sk0 + nθ ∗2

× I ( min Sk ≤ an, Sn > bn).
0≤k≤n







44

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

We write An = {min0≤k≤n Sk ≤ an, Sn > bn} and introduce thesets
Bε (n) = An ∩ Sk0 + nθ ∗ ≥ nε ,Cε (n) = Bε (n) ∩ Sk0 − Sn − Sk0 + nθ ∗ ≥ nε .









Then,






 



E [R2n ] = exp −nθ ∗2 
E exp 2θ ∗ Sk0 − Sn − Sk0 + nθ ∗ ; An





 


E exp 2θ ∗ Sk0 − Sn − Sk0 + nθ ∗ ; Bε (n)≥ exp −nθ ∗2 ≥ exp −nθ ∗2 exp 4θ ∗ ε n 
P (Cε (n)) .
Next, because Cε (n) involves a deviation of size ε in supremum
norm relative to the optimal trajectory, it is not difficult to verify
that
lim

log 
P (Cε (n))
n

n→∞

 
= −2θ ∗ ε + O ε2 .

So,
lim

1

n→∞ n

log 
E [R2n ] ≥ −θ ∗2 + 2θ ∗ ε + O ε 2 ≥ −θ ∗2 + δ,

 

(15)

 
if one picks ε small enough such that 2θ ∗ ε + O ε 2 ≥ δ > 0 forsome δ > 0.
Combining (14) together with (15) we conclude that Rn is not ingeneral even weakly efficient.
3.3. Notes on state-independent importance sampling and large
deviations
Theory of large deviations from a function space and sample
path perspective is comprehensively documented in [25]. On
more specific application areas, [33] discusses large deviations inqueueing and computer science, while [34] focuses on applications
in statistical mechanics. Ref. [35] provides an alternative approach,based on control theory and weak convergence analysis, to develop
large deviations results.
The types of importance samplers discussed in Examples 1 and2 originate from the work of [29] in the context of sequential
analysis. Another early reference in this spirit is [36]. The paper [37]develops and proves several general results on importancesamplers for Markov chains and semi-Markov processes. Levelcrossing probabilities on random walks, as discussed in Examples 1
and 2, have important connections to insurance and performance
evaluation of queueing systems. For the former, readers arereferred to [38,39] and Chapter X of [40]. The paper [41] studies
the extension to Markov-modulated walks on general state space.
For queueing, [28] discusses exponential change-of-measure inseveral applications, and points out a duality link to steady-state
estimation via level-crossing problems. Other related papers onchange-of-measure from a random walk perspective applied toqueueing models include [42–44].
The fact that reasonable importance samplers based on large
deviations results can lead to efficiency deterioration was firstpointed out in the 90’s. The paper [1] studies the problem in thecontext of tandem queues, which is further analyzed in [45]. The
issue is elaborated and strengthened in several examples in [2].
Ref. [46] and Chapter 5 of [47] discuss such inefficiency from socalled dominating point geometry in a general setting, while [48]extends the conditions using Varadhan’s lemma. The recent workof [6] aims to remedy such problems in generality, drawing upon
the theory of subsolution in Isaacs equation (this will be our main
topic of discussion in the next section).
For general review on rare-event simulation and exponentialchange-of-measure, readers are referred to [49,47] and Chapter
VI of [23]. The work of [50] provides a comprehensive discussion

on the active research directions and list of references. For easeof explanation of recent advances and as supplement to [50],in this paper we focus mainly on random walk problems, andwe emphasize that the scope of applications for importancesampling is much wider. For more specific references on general
queueing applications, readers are referred to [51,52]. In thecontext of infinite and many-server systems especially, recent
work include [53–55]. Another related application area is thereliability estimation of highly dependable Markovian systems, ofwhich [51] again provides an excellent list of earlier references.
More recent documentation is in [56]. Other applications ofimportance sampling are in finance (see for example [32]), physics(see for example [57]), systems biology (see [58]), combinatorics
and counting (see [59,60,12]), among others.
We also want to mention that besides importance sampling,
another powerful tool for variance reduction in rare-event problems is known as the splitting algorithm. Rather than modifying
the underlying probability distribution, splitting algorithm divides
the simulation trajectory space into levels (more precisely a nested
sequence of subsets), with the hope that the success rate of passing through each level stands high even though the target probability is rare. Sample trajectories then split or branch out at eachlevel to enhance arrivals at the next level, and finally up to thetarget set of interest. Hence the weights of the samples are basedon the calculation of the associated branching process instead ofa change-of-measure. Interested readers can go to [61] and Chapter VI Section 9 of [23]. The use of large deviations in analyzing
splitting algorithms is studied in [62,63]. After that, [64] proposes a
subsolution approach which is also closely related to our next section, while recently [65] discusses such algorithms in counting andcomputer science. Related methodology for the case of a fixed large
deviations parameter, based on Feynman–Kac flows which accommodate most of the particle methods used in practice, is studied
in [66].
Finally, besides drawing upon the theory of large deviations,
importance samplers, with the goal of minimizing estimation
variance, can also be implemented from an empirical perspective.
These are schemes whose parameters are unknown and are
updated from past samples in the simulation (note that theset of parameters can be huge, such as all the entries ofa transition matrix in a Markov chain). As opposed to the
optimal exponential change-of-measure in Examples 1 and 2 (andalso the state-dependent samplers in the next sections) whose
exponential parameters can be solved by looking at the associatedoptimization or variational problems, these empirical Monte
Carlo schemes do not have analytical solution to the change-ofmeasure parameters, and hence empirical measurement is needed.
Under appropriate assumptions, they can be proved to reduce
variance or lead to exponential convergence of the estimate. Some
important approaches in such setting are adaptive Monte Carlo(see e.g. [67,68]), the cross-entropy method (see [69–72]), andstochastic approximation based method (see [73,74]). Moreover,recent work of [75] considers an interesting importance samplingwith resampling approach for Markov-modulated problems thatavoids the need for computing the associated eigenfunctions.
4. State-dependent techniques for light-tailed systems
4.1. Solutions to a deterministic control problem and efficientimportance sampling
In order to deal with the types of issues illustrated by Example 2
above [5] proposed a control theoretic approach for the problem
of designing efficient importance sampling estimators. We will
formally characterize the solution to a deterministic control
problem proposed by Dupuis and Wang [5] in the setting ofmultidimensional random walks.

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

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(⟨Xk , θ⟩) isfinite everywhere. We then let Sk = X1 + · · · + Xk and consider ascaling similar to that introduced in Mogulskii’s theorem. We let
Yn (t ) = y + S⌊nt ⌋ /n,
for some y ∈ Rd . Let B be a closed set with non-empty interior
and define τB (n) = inf{t > 0: Yn (t ) ∈ B}. We assume that B isattainable in the sense that P (τB (n) < ∞|Yn (0) = y) > 0 for eachy ∈ Rd .
Consider the question of computing via simulation for two
disjoint sets B and C ; in addition to the regularity imposed on B
we will impose an additional large deviations requirement below.
We are interested in computing

αn (y) = Py (τB (n) < τC (n), τB (n) < ∞):= P (τB (n) < τC (n), τB (n) < ∞|Yn (0) = y).

∫ t

θ

First order optimality conditions imply that at the optimal value∇ψ(θ ∗ (y)) = ∇ψ(−∇ H (y) − θ ∗ (y)),

nτB (n)⟨θk , Xk ⟩ψ(θk ) I (τB (n) < τC (n), τB (n) < ∞),

(16)

k=1

where θk is adapted to the sequence Xj : j ≥ 1 in the sense that
θk is allowed to depend only on X1 , . . . , Xk−1 . The selection of
exponential tiltings as the family of controls is motivated by the
discussion involving the optimization problem (8).
Let Vn (y) denote the optimal second moment of Rn with y being
the initial position of the random walk. It follows immediatelythat the HJB (Hamilton–Jacobi–Bellman) equation corresponding
to finding the optimal policy (i.e. the optimal sequence of θk ’s) tominimize the second moment of Rn takes the formVn (y) = inf Eθ [exp(−2⟨θ , X ⟩ + 2ψ(θ ))Vn (y + X /n)]θ= inf E [exp(−⟨θ , X ⟩ + ψ(θ ))Vn (y + X /n)]θ(17)

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)Vn (y) = E exp −−⟨θk∗ , Xk ⟩k=1nτB (n)+−ψ(θk∗ ) I (τB (n) < τC (n), τB (n) < ∞)k=1nτB (n)≤ E exp −−⟨θk , Xk ⟩k=1nτB (n)+−k=1

which yields θ ∗ (y) = −∇ H (y)/2 and therefore we conclude thatEq. (19) can be expressed as2ψ(−∇ H (y)/2) = 0,

(20)

subject to the boundary conditions (inherited from (17)) equal to
H (y) = 0 for y ∈ B and H (y) = ∞ for y ∈ C . Eq. (20) is an instanceof a so-called Isaacs equation.
On the other hand, we have assumed

αn (y) = exp(−nIB,C (y) + o(n))

k=1

+

(19)

θ ∗ (y) one has

See for instance the text of [35, p. 66] or [25, p. 176].
We are interested in choosing a suitable importance sampling
estimator with optimal performance in the sense of minimizing its
second moment. The class of policies that we consider are based on
exponential tilting so that the estimator ultimately takes the form

−

(18)

θ

= min[ψ(θ ) + ψ(−∇ H (y) − θ )].



nτB (n)θ+ ψ(θ ) − nH (y) − ⟨∇ H (y), X ⟩)].

θJ (ẋ(s)) ds: x(0) = y, x(t ) ∈ B for some−θ≈ min log E [exp(−⟨θ , X ⟩= min log exp(ψ(θ ) + ψ(−∇ H (y) − θ ))t < ∞ & x (s) ̸∈ C if s < t .

Rn = exp −− nH (y) ≈ min log E [exp(−⟨θ , X ⟩ + ψ(θ ) − nH (y + X /n))]0 ≈ min log E [exp(−⟨θ , X ⟩ + ψ(θ ) − ⟨∇ H (y), X ⟩)]0

for any adapted policy (θk : k ≥ 1) (adapted in the sense describedearlier below display (16)).
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 the previous HJB equation, after taking logarithms, formally andwithout being careful about underlying smoothness assumptions
and errors incurred, we arrive at the approximation

Equivalently, we have that

We assume that B and C are regular in the sense that1lim
log αn (y) = −IB,C (y),
n→∞ n
where
IB,C (y) = inf

45

as n → ∞ and, by conditioning on the first increment of therandom walk, we also have thatαn (y) = E αn (y + X /n)subject to the constraints that αn (y) = 1 for y ∈ B and αn (y) = 0
if y ∈ C . Proceeding to analyze equality (21) formally as we did forthe discussion leading to (20), we conclude thatexp(−nIB,C (y) + o(n)) = E exp(−nIB,C (y + X /n))

≈ E exp(−nIB,C (y) − ⟨∇ IB,C (y), X ⟩).
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 (20) 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. Sufficientconditions required to rigorously substantiate this result are givenin the following theorem formulated by Dupuis and Wang [5] inthe case of large deviations problems for Sn /n.
Theorem 2 (Adaptation From [5] Theorem 2.1). Let Sn , Yn and ψ(θ )be defined as in the beginning of this section. Consider the problem of
estimating P (Sn /n ∈ A|S0 = 0) where A is a closed set that satisfies
inf I (y) = inf I (y)y∈A◦ψ(θk ) I (τB (n) < τC (n), τB (n) < ∞),

(21)

y∈A

and I (y) = supθ∈Rd [⟨θ , y⟩ − ψ(θ )]. Let us use ∇x for the gradientof a function with respect to the component x. Suppose that there

46

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

is a classical sense solution corresponding to the Isaacs equation
analogous to (20) given by

It follows that on {τB,C > k},∂t H (x, t ) = 2ψ (∇x H (x, t ) /2)

E [Mk+1 |(Sj : 1 ≤ j ≤ k)] = exp −subject to the boundary condition limt ↘0 H (x, t ) = ∞ × I (x ̸∈ A)(by convention we consider ∞× 0 = 0). The optimal state-dependent
exponential tilting scheme with k-step exponential parameter θk∗ =∇x H (Sk /n, k/n) generates the estimator




n
−
∗
∗
exp
(−⟨θk−1 , Xk ⟩ + ψ(θk−1 )) I (Sn /n ∈ A).



j =1



k

+

−

ψ(θj ) EUn (Sk+1 /n) exp(−⟨θ , Xk+1 ⟩ + ψ(θk+1 ))

j =1




k
k
−
−
≤ exp −
⟨θj , Xj ⟩ +ψ(θj ) Un (Sk /n).
j=1

k=1

Let Vn (0) be the second moment of this estimator and let Hn (0) be
such that Vn (0) = exp(−nHn (0)). We have thatlim Hn (0) = 2 inf I (y) = H (0) .

n→∞

k
−
⟨θj , Xj ⟩j =1

Evidently, on {τB,C ≤ k},
E [Mk+1 |(Sj : 1 ≤ j ≤ k)]


τB,C
τB,C
−
−
ψ(θj ) .
⟨θj , Xj ⟩ += Un (SτB,C /n) exp −y∈AThe problem with the previous result is that it is rarely
applicable in most situations of interest, especially because it isoften the case that there is no classical solution to the Isaacs
equation. Other connections are discussed in [76]. The function
H (y) = 2IB,C (y) typically provides a solution in a weak sense (aso-called viscosity subsolution). The extent to which this condition
can be relaxed and the type of optimality results that can beobtained beyond weak efficiency is not fully understood.

Therefore,
E [Mk+1 |(Sj : 1 ≤ j ≤ k)]




k
k
−
−
≤ exp −
⟨θj , Xj ⟩ +ψ(θj ) Un (Sk /n)I (τB,C > k)
j=1

j =1


4.2. Subsolutions to a deterministic control problem and efficientimportance sampling
We take as starting point the HJB equation introduced in (17).
The next lemma provides an inequality whose solution gives anupper bound for the value function in (17).
Lemma 1. Suppose that one finds a non-negative function U (·) suchthat
Un (y) ≥ E [exp(−⟨θ , X ⟩ + ψ(θ ))Un (y + X /n)]≥ inf E [exp(−⟨θ , X ⟩ + ψ(θ ))Un (y + X /n)],

j =1

j =1

+ Un (SτB,C /n) exp −τB,C−⟨θj , Xj ⟩
j =1

+

τB,C
−


ψ(θj ) I (τB,C ≤ k)

j=1

= Mk
and, if S0 = ny,

Un (y) ≥ E Un (SτB,C ∧k /n) exp −τ−B,C ∧kθwhere y ̸∈ B ∪ C and subject to the boundary condition that Un (y) ≥ 1
for y ∈ B. Then,

+

τB,C
−



ψ(θj ) I (τB,C < ∞, τB < τC )are defined as in the beginning of Section 4.1.
Proof. Recall τB,C = n × (τB (n) ∧ τC (n)) and fix any policy (θk : k ≥0) for which
Un (y) ≥ E [exp(−⟨θ , X ⟩ + ψ(θ ))Un (y + X /n)],and τB,C < ∞. Consider the process
Mk+1 = Un (SτB,C ∧(k+1) /n) exp −τB,C−∧(k+1)⟨θj , Xj ⟩

j =1

τB,C−
∧(k+1)

τB,C−

By Fatou’s lemma, we obtain thatj =1

We claim that (Mk : k ≥ 0) is a supermartingale. Note thatE [Mk+1 |(Sj : 1 ≤ j ≤ k)] = E [Mk+1 |(Sj : 1 ≤ j ≤ k)]I (τB,C > k)

+ E [Mk+1 |(Sj : 1 ≤ j ≤ k)]I (τB,C ≤ k).

Un (y) ≥ E Un (SτB,C /n) exp −τB,C−⟨θj , Xj ⟩
j=1

+

τB,C
−





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

j =1

Using the boundary condition on B and the fact that Un (y), being
non-negative, must satisfy that Un (y) ≥ 0 on C we obtain

Un (y) ≥ E exp −τB,C−⟨θj , Xj ⟩j =1ψ(θj ) .

ψ(θj ) I (τB,C ≤ k) .


for any adapted policy θj : j ≥ 1 , where τB,C = n(τB ∧ τC ) and τB , τC



τB,C
−
⟨θj , Xj ⟩

j =1



j =1

+

≥ E Un (SτB,C /n) exp −



j =1

+



ψ(θj )

j =1

j =1





j=1

τB,C
−
⟨θj , Xj ⟩Un (y) ≥ E exp −⟨θj , Xj ⟩ +τ−B,C ∧k

+

τB,C
−





ψ(θj ) I (τB,C < ∞, τB < τC ) . 

j =1

Using a similar formal argument as the one introduced in theprevious section in (18) we write Un (y) = exp(−nGn (y)) andJ. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

postulate Gn (y) → G(y) as n → ∞ for some function G (y). If
we proceed using this postulated limit in the inequality obtainedin the previous lemma, after taking logarithms, we arrive at−nG(y) ' min log E [exp(−⟨θ , X ⟩ + ψ(θ ) − nG(y + X /n))]θ

We can take advantage of Example 1 in order to obtain
the asymptotics for αn . In particular, define τi (n) = inf{t > 1:
⟨Y (t ), ai ⟩ ≥ 1} for each i ∈ {1, . . . , m} and note that one has theelementary inequalitiesmmax P (τi (n) < ∞|Yn (0) = y)≈ min log E [exp(−⟨θ , X ⟩ + ψ(θ ) − nG(y)

i =1

θ

− ⟨∇ G(y), X ⟩)].

≤ αn (y) ≤m−

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

(22)First order optimality conditions imply that at the optimal value∇ψ(θ ∗ (y)) = ∇ψ(−∇ G(y) − θ ∗ (y)),
which yields θ ∗ (y) = −∇ G(y)/2 and therefore we conclude that
inequality (22) can be expressed as(23)

for y ̸∈ B ∪ C and subject to the boundary conditions (inheritedfrom Lemma 1) equal to G(y) ≤ 0 for y ∈ B. Inequality (23) isthe corresponding subsolution to the Isaacs equation developed in(20). In order to conclude asymptotic optimality for αn (y0 ) (here
we let y0 as the initial position) we must have G(y0 ) ≥ 2IB,C (y0 ).
Sufficient conditions required to rigorously substantiate this result
are given in the following theorem.
Theorem 3 (Adaptation from [6] Theorem 8.1). Let Xk and Sk be
defined as usual, with ψ(θ ) = E exp(⟨θ , Xk ⟩) < ∞ for any θ ∈ Rd .
Suppose that G(y) satisfies (22) with the boundary condition G(y) ≤0 for y ∈ B. Let θ ∗ (y) = −∇ G(y)/2. Then, the state-dependent
sampler using k-th step exponential parameters θ ∗ (Sk /n) given by

τ
B, C
−

Rn = exp−⟨θ ∗ (Sk−1 /n), Xk ⟩k=1

+ ψ(θ (Sk−1 /n))
∗




n→∞

1
n

(24)

I (τB,C < ∞, τB < τC )Evidently, P (τi (n) < ∞|Yn (0) = y) is equivalent to the one dimensional first passage time probability discussed in Example 1. To
see this, one simply needs to work with the one dimensional random walk with iid increments (X1 (i), X2 (i), . . .) defined via Xk (i) =⟨Xk , ai ⟩. We then obtain that
P (τi (n) < ∞|Yn (0) = y) = exp[−nθi∗ (1 − ⟨y, ai ⟩) + o(n)]log E [R2n ] ≥ G(y0 )where S0 = y0 . Consequently, if G(y0 ) ≥ 2IB,C (y0 ), then Rn is an
asymptotically optimal estimator.
We now apply the previous result to a couple of examples.
Example 3 (Multidimensional First Passage Time Problems). Assume
that µ = EXi ̸= 0 and consider vectors a1 , . . . , am ∈ Rd suchthat ⟨µ, ai ⟩ < 0 for all i ∈ {1, . . . , m}, where ⟨x, y⟩ is used todenote the inner product between the vectors x, y ∈ Rd . Define
the event B = ∪mi=1 {y: ⟨y, ai ⟩ ≥ 1} and consider the first passage
time probability

αn (y) = P (τB (n) < ∞|Yn (0) = y),where τB (n) = inf{t ≥ 0: Yn (t ) ∈ B}. Suppose that for eachi ∈ {1, . . . , m} there exists θi∗ ̸= 0 such that ψ(ai θi∗ ) = 0 and⟨ai , ∇ψ(ai θi∗ )⟩ < ∞. Just as we explained in Example 1, since
EaTi Xk < 0 we must have that θi∗ > 0 and also that⟨ai , ∇ψ(ai θi∗ )⟩ ∈ (0, ∞). This problem is the natural generalization of Example 1 and can also be interpreted in terms of an insurance formulation in which one computes a ruin probability of
an insurance company with several lines of business and for which
one is allowed to borrow some resources from one line to anotherone. This formulation is discussed in [77].

(25)

as n → ∞ and therefore,
lim

1

n→∞ n

m

log αn (y) = − min θi∗ (1 − ⟨y, ai ⟩).
i =1

The natural state-independent technique, analogous to that of Example 1, consists in computing j such thatmθj∗ (1 − ⟨y, aj ⟩) = min θi∗ (1 − ⟨y, ai ⟩)i=1

and applying exponential tilting according to the tilting parameter aj θj∗ . Unfortunately, just as we showed in Example 2, it turns
out that this strategy is not guaranteed to yield an efficient estimator; such a counter-example can be constructed using a two dimensional Gaussian random walk using a lower bound analogous
to that developed for Example 2; the details are left as an exercise
to the reader.
We now use the subsolution approach to construct an importance sampling change-of-measure for this problem. Ultimately
we need to find a function G(·) satisfying the inequality (23) together with the corresponding boundary condition that G(y) ≤ 0
for y ∈ B. Moreover, since in the end we are after 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 (24) and (25) to postulate

has second moment satisfying
lim inf −P (τi (n) < ∞|Yn (0) = y).

i=1θ ∗ (y) one has

0 ≥ 2ψ(−∇ G(y)/2),47
exp(−nG(y)) =m−2exp[−nθi (1 − ⟨y, ai ⟩)]
∗

.

i =1

We now verify that our selection is indeed a valid subsolution. First,
note that∇ G(y) = −2m−θj∗ aj pj (y, n),

j =1

where
pj (y, n) =exp[−nθi∗ (1 − ⟨y, ai ⟩)]m∑.

exp[−nθi (1 − ⟨y, ai ⟩)]∗i =1Note that for each y, (p1 (y, n), . . . , pm (y, n)) can be interpreted
as a probability vector. Moreover, observe that as n → ∞, theprobability vector degenerates and concentrates its mass on theindexes corresponding to the tilting parameters θj∗ that drive themost likely way in which the random walk hits the target set. Now,
using Jensen’s inequality and convexity of ψ (·) we conclude thatψ(−∇ G(y)/2) = ψm−
θj aj pj (y, n)
∗

j =1

≤

m
−
j =1

pj (y, n)ψ(θj∗ aj ) = 0.

48

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

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

(26)

j =1

for k/n < τB (n) is asymptotically optimal.
Remark 2. Note that we have restricted ourselves to importancesampling estimators of exponential tilting form. Nevertheless, theform of (26), which is expressed in terms of a mixture, suggests
considering mixtures of exponential tiltings, or in other words,
distributions of the form
Q (Xk ∈ dy) =m−pj (Sk−1 /n, n) exp(θj∗ ay)P (Xk ∈ dy).

(27)

j=1

Dupuis and Wang [6] also study these types of ‘‘generalized’’
controls which are asymptotically equivalent (in the sense ofachieving also weak efficiency) to those in (26), see also [7]. The
paper [78] shows that using Lyapunov inequalities the sampler(27) actually yields strong efficiency.
The subsolution approach is better appreciated in situations inwhich the system exhibits piecewise random walk dynamics, as
it happens in the setting of Jackson networks. The next example
illustrates this situation in the context of a two dimensionalrandom walk.
Example 4 (A First-Passage Time Problem With Discontinuous
Statistics). Let (Zk : k ≥ 1) be a sequence of iid two dimensional
standard Gaussian random variables (i.e. having zero mean andidentity covariance matrix). Let µ1 = (cos(9π /8), sin(9π /8))Tand µ2 = (cos(11π /8), sin(11π /8))T . We let W = (Wk : k ≥ 0)
be a two dimensional process defined via
Wk+1 = Wk + Zk+1 + I (Wk (1) < Wk (2))µ1 + I (Wk (1)≥ Wk (2))µ2 ,with W0 = ny, where y = (y(1), y(2))T is such that y(1) < 1 andy(2) < 1. Define Yn (t ) = W⌊nt ⌋ /n, B = {x = (x(1), x(2))T : x(1) ≥1 or x(2) ≥ 1} and consider the problem of efficiently estimating
via simulation

αn (0) = P (τB (n) < ∞|Yn (0) = 0),
where τB (n) = inf{t > 0: Yn (t ) ∈ B}.
Although the process (Wk : k ≥ 0) is not a random walk,note that it behaves like a Gaussian random walk above andbelow the 45° line (i.e. the line y(1) = y(2)); only the driftchanges in each of these two regions. We have selected the drift
parameters to simplify our analysis and flesh out the main issue
arising in constructing the subsolution in the current setting. A
picture illustrating the drift of the process in the whole space isgiven next (see Fig. 1).
In this case it is not difficult to convince ourselves that if y(1) <y(2) < 1, then
lim

1

n→∞ n

log αn (y) = 2µ1 (2)(1 − y(2)).

and if 1 > y(1) ≥ y(2), then
lim

1

n→∞ n

log αn (y) = 2µ2 (1)(1 − y(1)).

Fig. 1. Vector field of fluid paths.

The previous limits correspond to the first passage time problem ofa one dimensional random walk with Gaussian increments having
negative drift equal to µ1 (2) = sin(9π /8) = µ2 (1).
A similar development to that leading to (23) leads to thesubsolution inequality for G (·) given by2ψ(y, −∇ G(y)/2) ≤ 0,
where, for y = (y(1), y(2))T ,ψ(y, θ ) = log E exp (⟨Z1 , θ⟩
+ I (y(1) ≥ y(2)) ⟨µ1 , θ⟩ + I (y(1) < y(2)) ⟨µ2 , θ⟩)= ‖θ ‖2 /2 + I (y(1) ≥ y(2)) ⟨µ1 , θ⟩
+ I (y(1) < y(2)) ⟨µ2 , θ⟩ ,and the corresponding boundary condition G(y) ≤ 0 for y ∈ B.
Now, let us first attempt a construction of the subsolution
similar to that given in Example 3, namely by considering thefunction G(·) such that2G(y) = − log(exp[2nµ1 (2)(1 − y(2))]
n+ exp[2nµ2 (1)(1 − y(1))]).
Note that in this case we have that− ∇G(y)/2 = 
p1 (y, n)θ1∗ +
p2 (y, n)θ2∗ ,

(28)

where


p1 (y, n) =exp[2nµ1 (2)(1 − y(2))]
exp[2nµ1 (2)(1 − y(2))] + exp[2nµ2 (1)(1 − y(1))];θ1∗ = −(0, 2µ1 (2)),
andp2 (y, n) =exp[2nµ2 (1)(1 − y(1))]
exp[2nµ1 (2)(1 − y(2))] + exp[2nµ2 (1)(1 − y(1))];θ2∗ = −(2µ2 (1), 0).
Observe that

ψ(y, −∇G(y)/2) ≤ 
p1 (y, n)ψ(y, θ1∗ ) +p2 (y, n)ψ(y, θ2∗ )≤ 2
p1 (y, n)I (y(1) ≥ y(2))+ 2p2 (y, n)I (y(2) > y(1)).
Therefore, if |y(1) − y(2)| ≥ δ for some δ > 0 thenψ(y, −∇G(y)/2) ≤ 2 exp(2δ nµ1 (2)).

(29)

So, while the subsolution inequality is not satisfied, the size of theviolation converges to zero exponentially fast as n → ∞. On theother hand, we can easily deduce that

ψ(y, θ1∗ )
p1 (y, n) + ψ(y, θ2∗ )
p2 (y, n) > δ > 0
on a region of size O(1/n) around the 45° line, namely the liney(1) = y(2). While this region might look small, note that in the

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

49

Fig. 2. Plot of the function max −4µ1 (2)y(1), −4µ1 (2)y(2), 2 θ ′ , y + ∆ .









Fig. 4. Plot of the level curves {θ : ψ (θ, y) = 0} for y = (y (1) , y (2)) withy (1) < y(2) and y(1) > y(2).

we indicated before, we also wish to select ψ(θ0∗ , y) ≤ 0. So, it isuseful to plot the functions ψ(·, y) for y(2) > y(1) and y(2) < y(1)respectively, as we do next.
F0 (y) = −4µ1 (2) max(y(1), y(2)).
Fig. 3. Plot of the function 

original scaling of the process W the region is of order O(1), so theprocess might spend a significant amount of time (in fact a randomamount of time of order O(1)) in a region where the likelihood
ratio of the importance sampling distribution is not appropriately
controlled.
The main difficulty with G(y) arises precisely because of the
discontinuous dynamics of the process W right on the line y(1) =y(2). In order to deal with this problem we note thatlim G(y) = −2 max {⟨θ1 , y⟩ + µ2 (1), ⟨θ2 , y⟩ + µ1 (2)}n→∞

= 2µ2 (1) max{y(2), y(1)} − 2µ2 (1).
(30)
We will need to slightly modify the function G(·) right on the set
y(1) = y(2) and in such a way that one can obtain a smoothfunction G(·) whose gradient can be represented as a weighted

The intersection of the areas enclosed by both circles shown in
Fig. 4 above correspond to the points θ for which ψ(θ , y) ≤ 0,for each y ∈ R2 . A natural selection for θ0∗ (although clearly not the
only one) is θ0∗ = (−µ2 (1), −µ2 (1))T . Clearly, since −µ2 (1) ∈(0, 1) we have that ψ(θ0∗ , y) < 0 for each y ∈ R2 and such a
selection indeed corresponds to the supporting hyperplane thatmakes the plot of Fig. 2 symmetric around the line y(1) = y(2).
We now have all the necessary elements required to define our
subsolution. In particular, by virtue of our considerations leading
to the definition of 
F0 (·) and our selection of θ0∗ , we let
G(y) = −2nlog(exp[n(f1 (y) + µ2 (1))] + exp[n(f2 (y) + µ2 (1))]+ exp[n(f0 (y) + µ2 (1))]).
The importance sampling scheme is generated by using as
exponential tilting parameter at time k:sum, just as in (28). In order to do this, definef1 (y) = 2 θ1∗ , y ,
f2 (y) = 2 θ2∗ , y ,2−

k−1 /n) =θk (Wk−1 /n, n),θj∗ pj (Wj =0

and define the convex functionF0 (y) = max{
f1 (y),
f2 (y)} = −4µ2 (1) max{y(2), y(1)},which is just a translation to the origin of the limit obtained in (30).
The plot of 
F0 (·) is given in Fig. 3.
The idea is to introduce an affine function 
f0 (y) = 2 θ0∗ , y +∆, for some ∆ > 0, so that we can define corresponding
probabilities p0 (y, n), p1 (y, n) and p2 (y, n), associated to θ0∗ , θ1∗and θ2∗ , respectively, with the property that p1 (y, n), p2 (y, n) =o(p0 (y, n) exp(−ε n)) for some ε := ε(∆) > 0 uniformly over ysuch that |y (1) − y (2)| ≤ δ := δ(∆) and such that ψ(y, θ0∗ ) ≤0 for all y. We also wish the behavior of p0 (y, n), p2 (y, n) =o(p1 (y, n) exp(−ε n)) uniformly over y(2) > y(1) + δ and similarly
p0 (y, n), p1 (y, n) = o(p2 (y, n) exp(−ε n)) uniformly over y(1) >y(2) + δ . If in addition, we select ψ(θ0 , y) ≤ 0, just as we obtainedin the analysis (29), we will be violating the subsolution only by an
exponentially small amount in n, which as it turns out, is enoughto conclude asymptotic optimality.
A natural approach that comes to mind to construct such afunction 
f0 (·) is by considering a supporting hyperplane to thefunction 
F0 (·) right along the axis y(1) = y(2). Let us denotesuch supporting hyperplane by f (y) = 2 ⟨η, y⟩. By shifting such
hyperplane an amount ∆ we obtain a picture such as the one
shown in Fig. 2.
All subgradients corresponding to supporting hyperplanes strictlyalong y(1) = y(2) are potential candidates for 2θ0∗ . However, ask−1 /n, n) is as given in Box I below. Under thewhere pj (Wimportance sampling distribution the process evolves according to
the Markov chaink+1 = W
k + 
W
Zk+1 + θk (Wk /n),0 = 0,W

where (Zk : k ≥ 1) is a sequence of iid two dimensional standard
Gaussian random variables. The importance sampling estimator,
namely,nτB (n)−Rn = exp −k−1 /n)(Wk − Wk−1 )θk (Wk=1nτB (n)+−k−1 /n, θk (Wk−1 /n)) I (τB (n) < ∞),
ψ(W

k =1

satisfieslim

log ERnn→∞ log

αn (0)≥ 2 + ∆/µ2 (1).

One can select ∆ := ∆n −→ 0 at a suitable speed inorder to conclude asymptotic optimality. These types of issues
are discussed in the next subsection on remarks and further
considerations on state-dependent importance sampling for lighttailed systems.

50

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

k−1 /n, n) =pj (Wk−1 /n) + µ2 (1))]
exp[n(
fj (W
,



k−1 /n) + µ2 (1))] + exp[n(k−1 /n) + µ2 (1))]
exp[n(f1 (Wk−1 /n) + µ2 (1))] + exp[n(f2 (W
f0 (Wwhenever k/n < τB (n).
Box I.

4.3. Notes on state-dependent importance sampling for light-tailed
systems
The discussion so far in this section is taken mainly from a seriesof papers by [5,79,6]. It is important to note, however, that wehave restricted ourselves to exponential changes-of-measure. They
formulate optimal state-dependent twisting as an HJB equationthat leads to a differential game and a corresponding Isaacs
equation in [5]. The case for uniformly recurrent Markov chain
is studied in [79], and the general method of finding subsolution
in a variety of examples is shown in [6]. Moreover, [7] studies
the subsolution approach to the tandem queue problem spottedby [1], as well as some generalizations. Importance sampling forJackson network using similar approach is done in [80]. The use ofviscosity solution in large deviations analysis of queueing networks
can be seen in [81]. Recently, rare event simulation estimatorsfor measure-valued light-tailed systems are developed in [55,82].
Developing Isaacs equations in these settings appears to be a very
challenging task.
Other related approaches and general conditions for asymptotically optimal importance samplers are studied in [83,46,84,
48]. Sharper analysis (beyond weak efficiency) of state-dependent(subsolution based) importance sampling estimators is given
in [55].

as n ↗ ∞. Moreover, for each Borel set A ⊂ Rd , define 
Pn (·) via


Pn ((X1 , . . . , Xm ) ∈ A) =5.1. Classical simulation techniques for heavy-tailed random walks
Let us start by collecting some useful definitions related toheavy-tailed systems (see, for instance, [96]).

P ((X1 , . . . , Xm ) ∈ A|Xj > n)

j =1

1
m

,

then
sup |P ((X1 , . . . , Xm ) ∈ A|Sm > n) − Pn ((X1 , . . . , Xm ) ∈ A)|
A−→ 0
as n ↗ ∞.
Proof. The first part follows by definition of subexponentiality.
Indeed, it follows by Bonferroni inequalities that
P (Sm > x) ∼ mF (x) ∼ P

max Xj > x

as x ↗ ∞.

(32)

1≤j≤m

For the second part, first note that
P ((X1 , . . . , Xm ) ∈ A|Sm > n)
= P (X1 , . . . , Xm ) ∈ A| max Xj > n
1≤j≤m


P max Xj > n
1≤j≤m
×P (S > n) m
P (X1 , . . . , Xm ) ∈ A, max Xj ≤ n, Sm > n
1≤j≤m
+
.
P (Sm > n)

5. Large deviations and importance sampling for heavy-tailed
random walks
We now move to the study of Monte Carlo methods for
random walks with heavy-tailed increments. Large and moderate
deviations for heavy-tailed random walks are studied in theclassical papers of [85–90]. More recent work on first passage andother related quantities include [91–94]. The paper [95] studies
conditional limit theorems and their applications to insurance risk
analysis.
Let us start by discussing classical simulation techniques in thissetting. This discussion will allow us to explain the basic principles
that govern the occurrence of rare events when heavy tails arepresent and it will also allow us to set the stage for state-dependent
techniques.

m−

Noting thatP
(X1 , . . . , Xm ) ∈ A, max Xj ≤ n, Sm > n
1≤j≤mP (Sm > n)P≤

max Xj ≤ n, Sm > n

1≤j≤mP (Sm > n)→0

as n → ∞ and using (32) we conclude that
P ((X1 , . . . , Xm ) ∈ A|Sm > n)
− P (X1 , . . . , Xm ) ∈ A| max Xj > n  → 0
1≤j≤m

uniformly over all Borel sets A as n → ∞. Similarly, it follows that
Definition 5. A non-negative random variable X is said to besubexponential if given m independent copies X1 , X2 , . . . , Xm of X
we have P (X1 + X2 + · · · + Xm > x) ∼ mP (X > x) as x ↗ ∞.
The subexponential property, as we shall see, provides insight
into 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 subexponential rv’s.

P

max Xj > n|Sm > n

1≤j≤m


−→ 1

(31)

 


(X1 , . . . , Xm ) ∈ A| max Xj > n
1≤j≤m
A



−
Pn ((X1 , . . . , Xm ) ∈ A)  −→ 0sup P

and therefore we conclude our result.

The previous proposition illustrates the so-called ‘‘catastrophe
principle’’ behind extreme behavior of heavy-tailed systems,
which postulates that large deviations are caused by extremes
in one or few components. In this case, the sum is large

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

because one single component, namely the maximum, is large. In
previous sections we saw that a convenient change-of-measure
for importance sampling is suggested by studying the asymptotic
conditional distribution of the underlying process given the rare
event of interest. Now, consider applying the idea behind suchconditional description to a simple example. Suppose that weare interested in estimating P (Sm > n) efficiently as n ↗ ∞.
Proposition 1 indicates that an asymptotic conditional descriptionof the Xj ’s given Sm > n assigns zero mass to sample paths forwhich all the random variables are less than n. As a consequence,
the natural asymptotic description of the Xj ’s given Sm > n is
singular with respect to the nominal (original) distribution andtherefore, contrary to the light-tailed case, a direct importancesampling approach is not feasible. This feature was observed by [4],which provides an extended discussion of the difficulties that are
inherent to the design of efficient simulation estimators for heavytailed systems.
Perhaps the most relevant special subclass of subexponential
random variables is given by those whose right tail is regularly
varying, which we define next.
Definition 6. 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 any β > 0. Similarly,
X has a regularly varying left tail with index α > 0 if −X has aregularly 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 ).
In addition to regularly varying distributions, subexponential
random variables include models that incorporate Weibull-type
tails with index γ ∈ (0, 1), which satisfy F (x) ∼ c exp (−β xγ )for c , β > 0 as x ↗ ∞. Other popular special models include lognormal, log-gamma and t distributions (the last two models are
particular cases of regular variation).
We summarize a few important properties of regularly varying
distributions that will be used in our future development; see for
instance, [97] or [96].
Theorem 4. Suppose that X possesses a regularly varying right tailwith index α > 0; in particular, P (X > t ) := F (t ) = t −α L(t ). Then,(i) (Karamata’s theorem) If α > 1, then∫ ∞F (s)ds =t −α+1 L(t )α−1t(1 + o(1))(ii) (Breiman’s theorem) Let Z be positive random variable such that
EZ α+ε < ∞ for some ε > 0, then
P (ZX > t ) = t

The first provably efficient rare-event simulation algorithms forheavy-tailed random walk problems focused on the problem of
estimating

αn = P (Sm > n),(33)

efficiently as n → ∞. This calculation is motivated by theruin problem explained in Example 2 in the case of Poisson
arrivals. In this case, there exists a representation, known as
Pollaczek–Khintchine’s formula (see [30, p. 237]) which allows
to express the probability of eventual ruin, given that the initialreserve of the company is n, as P (SM > n), where M is asuitably defined geometric random variable independent of an
appropriately defined random walk (Sk : k ≥ 1).
The first weakly efficient algorithm for estimating (33) was
designed by [98] for regularly varying increments and it was basedon a conditional Monte Carlo idea (see also [4]). They note that
it is straightforward to have direct access to the distribution ofX(m) := max1≤j≤m Xj . Using the fact that the event {Sm > n} belongs
to the single jump domain as n → ∞ it is natural to condition on
the first m − 1 order statistics X(1) , . . . , X(m−1) corresponding to
the Xj ’s (j ≤ m) (i.e. integrating out the contribution of X(m) whichis the most relevant one). Following this idea, [98] provided thefollowing result.
Proposition 2. Assume that X1 has a density and a regularly varying
right tail with index α > 0. Define S(m−1) = X(1) + · · · + X(m−1) andset
Z0 (n) = P Sm > n|X(1) , . . . ., X(m−1)



=

F





n − S(m−1) ∨ X(m−1)



F X(m−1)αL(t )EZ (1 + o(1))

P (X − nx > nt |X > nx) →Proof. Note that
EZ0 (n)2 = E [Z0 (n)2 ; X(m−1) ≤ n/2]

+ E [Z0 (n)2 ; X(m−1) > n/2].
Now, let us denote the density of X1 by f (·). It follows that thedensity of X(m−1) , which we shall denote by f(m−1) (·) satisfies
fm−1 (x) = m (m − 1) F (x)m−2 F (x) f (x) .
Therefore,
E (Z0 (n)2 ; X(m−1) > n/2) ≤ P X(m−1) > n/2



as n ↗ ∞. In addition,
E (Z0 (n)2 ; X(m−1) ≤ n/2) =∫ n/m

E (Z0 (n)2 ; X(m−1) ∈ dx)

0

∫ n/2

1

(1 + t /x)−αn/m(iv) (Long tails) For any x ∈ R−→ 1

E (Z0 (n)2 ; X(m−1) ∈ dx).

Observe that on X(m−1) < n/m we have S(m−1) < (m − 1) n/m andtherefore (by Karamata’s theorem)∫ n/m

as n → ∞.

as t → ∞.

.

+(iii) (Pareto conditional excess tail) For any x > 0 we have that

P (X > t )Then, Z0 (n) is a logarithmically efficient estimator for P (Sm > n) as
n ↗ ∞.

as t → ∞.

P (X > t + x)= O(F (n)2 )as t → ∞.

−α51

E (Z0 (n)2 ; X(m−1) ∈ dx) ≤

0

∫ n/m

m2 F (n/m)2 f (x)

0

F (x)

dx

= O(F (n)2 log (n)).
The last integral piece (from n/m to n/2) is handled as in theprevious display. Logarithmic efficiency then follows as a result of
the previous estimates. 

52

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

Regular variation has been used extensively in the previous
result, so it may not be surprising that the estimator Z0 (n), in
the presence of other types of heavy-tailed assumptions such as
Weibull tails, fails to be logarithmic efficient. Improved Monte
Carlo estimators have been recently proposed by [99]. In particular,
they observe that in order to reduce some uncertainty, one can also
consider the index corresponding to the largest jump and note that
P (Sm > n) =

m
−
 

E P Sm > n, X(m) = Xj |X1 ,of Xj . In particular, if T is exponentially distributed with mean 1,
then Λ−1 (T ) is a copy of Xj and therefore, for appropriate θ , wecan define hazard rate tilting densities, fθ (·), viaf θ ( x) =exp (θ Λ (x)) f (x)

The corresponding hazard rate importance sampling estimator for
P (Sm > n) takes the form



j =1


. . . , Xj−1 , Xj+1 , . . . , Xm = mE P Sm > n, X(m) = Xm |X1 , . . . , Xm−1 .

  = exp (θ Λ (x)) f (x) (1 − θ) .

E exp θ Λ Xj

exp −Z 2 ( n) =m∑
θ Λ Xj 

j =1

I (Sm > n) .

(1 − θ )mTherefore, a natural conditional Monte Carlo estimator that one canconsider isUsing Eθ (·) to denote the expectation operator assuming that the
Xj ’s are iid with density fθ (·) we obtain

Z1 (n) = mP Sm > n, X(m) = Xm |X1 , . . . , Xm−1

Eθ (Z2 (n)2 ) = (1 − θ)−2m Eθ









= mF (n − Sm−1 ) ∨ max(Xj : j ≤ m − 1) .
The paper [99] proves the following result (the analysis is similar
to that given in the proof of Proposition 2 and therefore the details
are omitted).
Proposition 3. Assume that X1 has a density and regularly varying
right tails. Then, the estimator Z1 (n) is strongly efficient as n ↗ ∞.
In addition, if X1 has Weibull-type tails with index γ ∈ (0, .58) thenZ1 (n) is logarithmically efficient.
The previous proposition indicates that Z1 (n) can only beguaranteed to be efficient if γ is sufficiently small. In the designof efficient estimators in the heavy-tailed setting it often occursthat models with less heavy tails than, say regularly varying, tend
to require more and more information of all the components (notonly the largest one). It is partly due to this feature that the
majority of the efficient estimators developed for heavy-tailed
systems assume special characteristics (such as regular variation
or Weibullian-type tails). Notable exceptions are the algorithms byBlanchet and Glynn [9], Blanchet and Li [123] which are designed
to work basically for all subexponential distributions. Additional
conditional Monte Carlo algorithms have been proposed for thetransient distribution of an M /G/1 queue by [100]. Their estimatoris proved to be strongly efficient for regularly varying input.
The paper [4] discusses the difficulties of applying importancesampling to heavy-tailed problems—an important problem concerns the singularity issue that we mentioned in the paragraph after the proof of Proposition 1. Nevertheless, they also study ideasthat give rise to provably efficient algorithms. For instance, theynote that if the Xj ’s have a density and regularly varying tails withindex α > 0, then a logarithmically efficient importance sampling
scheme for estimating P (Sm > n) is obtained by sampling the Xj ’s
in an iid fashion according to the tail distribution 
P Xj > x =1/ log (e + x) for x ≥ 0. Logarithmic efficiency for this estimator
is very easy to obtain by means of Karamata’s theorem (see, forinstance, [23, p. 176]). This selection of importance sampling distribution biases the increments to induce very heavy-tailed distributions and therefore over-samples such sample paths for which
several (and not only the maximum) components contribute to theoccurrence of the rare event. However, this importance sampler,although logarithmically efficient, does not seem to perform wellin practice as reported by [23, p. 176].
Another importance sampling approach was suggested by [101]and is based on applying exponential tilting type ideas via thehazard rate corresponding to the Xj ’s. A basic observation behind
Juneja and Shahabuddin’s hazard rate tilting approachis thefact that if the Xj ’s have a positive density, then P Xj > t =exp (−Λ (t )), where Λ (t ) =t0

λ (s) ds and λ (·) is the hazard rate

× exp −2m−θ Λ Xj ; Sm > n .

j=1

Assuming that Λ (·) is a concave function (as is the case forstandard Pareto and Weibull random variables), then we obtainthat if Sm > n then
m
−
j =1


Λ Xj ≥ Λ

 

m
−


Xj

= Λ (Sm ) ≥ Λ (n)

(34)

j =1

and therefore, if θ ≥ 0, we obtain
Eθ (Z2 (n)2 ) ≤ (1 − θ)−2m exp (−2θ Λ (n)) .
Taking θ = 1 − η/Λ (n) for some η > 0 yields the following result.
Proposition 4. If the Xj ’s are non-negative and have a concave
cumulative hazard rate function Λ (·) then Z2 (n) is a logarithmicallyefficient estimator for P (Sm > n) as n ↗ ∞.
The previous hazard rate tilting strategy has been improved forthe case in which the number of increments follows a geometric
distribution as in the M /G/1 setting. Such approaches involve
suitable translation of the function Λ (·) applied when tilting. Thisadjustment is convenient in order to apply the concavity argument
given in (34) at a more convenient location relative to the rare
event Sm > n (see [50] and references therein). The hazard ratetilting idea has inspired further studies in the field because it tries
to develop rare-event simulation methodology through a structurethat resembles that of light-tailed input systems (via exponential
tiltings). Nevertheless, virtually all estimators that take advantageof this idea utilize the random walk structure substantially andsome sort of subadditivity argument on Λ (·) as we did in (34).
5.2. Large deviations for heavy-tailed random walks
As we illustrated by means of the catastrophe principle, thelarge deviations theory for heavy-tailed random walks has acompletely different qualitative characteristics than its light-tailed
counter-part. Intuitively, one can explain this because in lighttailed situations large deviations occur, as we have seen, becauseof a ‘‘conspiratorial effect’’ among all the increments; each ofthem follows a specific biased (exponentially tilted) distributionthat tracks the ‘‘optimal’’ trajectory. In the heavy-tailed setting, incontrast, large deviations occur because only one or few jumps
contribute to the occurrence of the rare event; the rest of thejumps evolve according to the original / nominal distribution. Thisintuition has been made rigorous in considerable generality onlyin situations in which only one jump contributes to the occurrenceof the rare event in the regularly varying setting [102]. (We shall

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

refer to this setting as ‘‘the single jump domain’’.) A few examples
of interest in which large deviations happen by the occurrenceof more than one jump (‘‘multiple jump domain’’) have been
rigorously studied mostly in the setting of queueing systems (see[103–105]).
Throughout the rest of our technical discussion we shall concentrate on regularly varying distributions. We will provide further
remarks on other types of heavy-tailed increment distributions at
the end of the section.
In spite of not having a large deviations theory that is ascomprehensive as that of light-tailed systems, in the heavy-tailed
case there is a powerful heuristic that allows to ‘‘guess’’ a plausible
form for large deviations results and the associated ‘‘optimal path’’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. Since, as we shall
explain, the methodology for constructing the algorithms can beused to rigorously justify (up to a constant) the heuristics, we
content ourselves with explaining how to apply the fluid heuristic
in a few examples and point to relevant references when available.
In order to guide our discussion we shall revisit a couple ofexamples studied earlier in the light tailed case.
Example 5 (First Passage Time Probabilities For Heavy-Tailed Random Walk). We consider once again the surplus process introducedin Example 1. In this case, however, we shall assume that the claim
sizes, namely the Zj ’s, are regularly varying with index α > 1.
We define Xj = Zj − pTj , where Tj (the j-th inter-arrival time),
is exponentially distributed with unit mean. Assume that EXk =−µ < 0. Suppose that S0 = y, let τn = inf{k ≥ 0 : Sk > n} andconsider the problem of approximating the ruin probability

αn (y) = P (τn < ∞|S0 = y)
as n ↗ ∞.
Because the service times are heavy-tailed we expect that only
a few jumps would contribute to the rare event. By the geometryof the problem, however, it is natural to expect that only one jump
will contribute to ruin and thus we are in the single jump domain.
The fluid heuristic consists in substituting the random walk by
its fluid or Law of Large Numbers behavior prior to the big jump.
The associated fluid path in this case is S⌊nt ⌋ /n ≈ y(t ) = y − µt
and therefore the fluid heuristic suggests

αn (y) ∼∞−P (Xk+1 + y(k) > n)

k=0

∫ ∞
∼

The associated fluid path in this case is Sk /n ≈ y(k) := 0 andtherefore, given the geometry of the problem, it is not difficult toconvince ourselves that in the current situation the most likely
path leading to the rare event involves two jumps, therefore, we
are in the multiple jump domain.
The fluid heuristic then suggests

−

αn ∼

P (Xk1 + y(k1 − 1) < an)1≤k1 <k2 ≤n× P (Xk2 + y(k2 ) − y(k1 ) + Xk1 > nb|Xk1+ y(k1 − 1) < an)n(n − 1)∼
P (X1 < an)P (X2 + X1 > nb|X1 < an).
2

Once again, note that we replace the path in-between jumps by
the fluid trajectory and neglect the contributions of more than two
jumps. We can further simplify the previous expression using thefact that for each y < 0,
P (X1 < yn + an|X1 < an) → P (Za < y + a) :=1(y/a + 1)v,

as n → ∞. The previous limit indicates that the distribution of
X1 /n given that X1 /n < a can be well approximated by Za ≤ awith the above distribution. This observation, combined with theheuristic approximation developed earlier in turn suggests

αn ∼n(n − 1)
2

P (X1 < an)P (X2 > n(b − Za )),where Za is independent of X2 . By the symmetry of the tdistribution it follows that P (X1 < an) = P (X1 > |a| n). Further,since a t distribution with v degrees of freedom is regularly varying
with index v , it follows that P (X1 > |a| n) = |a|−v P (X1 >n)(1 + o(1)) as n → ∞. In this case, evidently, we have that
0 < 1/(b − Za ) < 1/(b − a) and therefore we can apply Breiman’s
theorem to conclude thatP (X2 > n(b − Za )) = E (b − Za )−v P (X2 > n)(1 + o(1)),
as n → ∞. We then conclude the following plausible asymptotic
approximation,

αn ∼ n2 P (X1 > n)2

|a|−v
2

E (b − Za )−v = Θ (n2 P (X1 > n)2 ),as n → ∞. We have not been able to locate a reference for thisasymptotic result in the literature, but it is not too difficult to show
rigorously its validity.

P (X > n + µt − y) dt
5.3. Conditional distributions and implications for simulation0

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 obtain

αn (y) ∼

53

1

∫ ∞

µ

n−y

P (X > u) dy.

as n ↗ ∞. This approximation turns out to be correct even beyond
the assumption of regularly varying service times (see [106–109]).
We now discuss an example in the multiple jump domain.
Example 6 (A Multiple Jump Large Deviations Event). We revisit
Example 2, now in the context of heavy-tailed random variables. In
particular, we assume that (Xk : k ≥ 1) is a sequence of iid randomvariables t distributed with v > 1 degrees of freedom and we areinterested in approximating the probability

αn = P [ min Sk ≤ an, Sn > nb]

As we discussed in Section 2, it is useful to understand theconditional distribution of the random walk given the occurrenceof the rare event of interest because such conditional distribution
provides the optimal change-of-measure (in terms of varianceminimization). A natural starting point, as we followed in the lighttailed setting 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 importancesampling distribution. Let us pursue this idea in the context of
Example 5.
Once again, we use the fluid heuristic to obtain an asymptotic
result for the conditional distribution of (X1 , . . . , Xk ) for any fixedk > 0 given that τn < ∞ and that S0 = 0 as n → ∞. Given+−∞ < x−
j < xj < ∞ for j ∈ {1, . . . , k} we note that
+
−
+
P x−
1 ≤ X1 ≤ x1 , . . . , xk ≤ Xk ≤ xk |τn < ∞



0≤k≤n

as n → ∞, where Sk = X1 + · · · + Xk and a < 0 < b.

+
−
+
P x−
1 ≤ X1 ≤ x1 , . . . , xk ≤ Xk ≤ xk , τn < ∞



=




αn (0)

54

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

E [I (x1 ≤ X1 ≤ x1 , . . . , xk ≤ Xk ≤ x+k )αn (Sk )]
−

=

+

−

αn (0)

If we apply the fluid heuristic locally, that is by replacing Xk+1 by
its mean, namely EXk+1 = −µ, we arrive at.

Note that αn (y) is increasing as a function of y. Therefore, letting
−
−
+
+
+
s−
k = x1 + · · · + xk and sk = x1 + · · · + xk ,αn (s−k )αn (0)

+−+≤ P x−
1 ≤ X1 ≤ x1 , . . . , xk ≤ Xk ≤ xk |τn < ∞

+
−
+
P (x−
1 ≤ X 1 ≤ x1 , . . . , xk ≤ X k ≤ xk ) ×+−+≤ P (x−1 ≤ X1 ≤ x1 , . . . , xk ≤ Xk ≤ xk ) ×

αn (s+k ).
αn (0)(35)

As we have seen in Example 5,

αn (y) ∼ G(n − y) :=

1

∫ ∞

µ

n −y

≈

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) .
Now,
P (Xk+1 ∈ dy, τn = k + 1|τn ∈ (k, ∞), Sk = s)
P (Xk+1 ∈ dy) I (y > n − s)
P (Xk+1 > n − s|s)

Therefore, applying the previous local version of the fluid heuristic
into (37) and combining this with (36) we arrive at the nonrigorous approximation
P (Xk+1 ∈ dy|τn ∈ [k + 1, ∞), Sk = s)

P (X > t ) dt

as n − y → ∞. Now, as a consequence of Karamata’s theorem(see for example [97, p. 17]), G (·) is also regularly varying (i.e.
exhibits basically power-law tail behavior) and therefore it follows
easily that G (n − y) /G(n) −→ 1 as n → ∞ for each y. We thenconclude (taking the limit in (35)) that the asymptotic conditional
distribution of (X1 , . . . , Xk ) given that τn < ∞ remains unchanged
as n → ∞.
The fact that the asymptotic conditional distribution of(X1 , . . . , Xk ) does not provide a useful description that can be used
as an importance sampling distribution in the prelimit is one of
the main reasons that makes the design of efficient importancesampling estimators for heavy-tailed systems a challenging
problem. This observation lies at the center of the discussion in
the paper by [4], which provides a number of counter-examples,in the spirit of the previous analysis, illustrating the challenges
that arise in the design of efficient Monte Carlo methods in heavytailed settings. In the same vein, the papers of [3,110] show thatno state-independent importance sampling estimator can achieve
weak efficiency.
Let us now take a more direct approach to approximate theincrement Xk+1 given the current position Sk = s of the random
walk at time k and given that k < τn < ∞. We note (adopting thenotation introduced in (5)) that

= p∗n (s)αn (s + Xk+1 )αn (s − µ)≈αn (s)
αn (s)(n − s + µ)−α+1≈ [1 + µ/(n − s)]−α+1≈(n − s)−α+1µ(1 − α)
1
= 1−
+o
.
n−s
n−s

,

(36)

P (Xk+1 ∈ dy) I (y > n − s)
P (Xk+1 ≤ n − s)+p(n − s)P (Xk+1 ∈ dy) I (y ≤ n − s)
P (Xk+1 ≤ n − s)

[1 − p(n − s)].

(38)

It turns out that this approximation is indeed valid in total variation[111], and it can be used to derive other conditional limit theorems,
such as the conditional distribution of τn given that τn < ∞. For
instance, applying the fluid heuristic we have that if τn > nt, thenSk ≈ −µk whenever k ≤ nt and therefore
P (τn > nt |τn < ∞) ≈nt∏
(1 − p(n + µk))k=1

≈ exp −nt−µ(α − 1)n + µkk=1
 ∫ ntµ(α − 1)
≈ exp −ds
n + µs0=

1(1 + µt )α−1(39)

as n → ∞. This result and analogous joint approximation results
involving τn and other quantities of interest such as the overshootand undershoot at the time of ruin have been derived, using ladderheight theory, by [95]; related extensions in continuous time are
given in [112]. The heuristic approach taken to derive (39) andother conditional limit theorems is made rigorous in [111]. It is
important to note that if 1 < α < 2, then the asymptotic
conditional distribution of the normalized hitting time τn /n giventhat τn < ∞ has infinite mean.
5.4. State-dependent importance sampling techniques for heavytailed random walkswherep∗n (s) =P (Xk+1 > n − s)α n ( s).

Using the approximation in Example 5, we have that
P (Xk+1 > n − s)µ(α − 1)p∗n (s) ∼ p(n − s) := µ  ∞∼(n − s)
P
X
>
t
dt
(
)
k
+
1
n−s
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)
P (Xk+1 ≤ n − s)× P (Xk+1 ≤ n − s)αn (s + y).
αn (s)(37)

In the previous section we argued, using a local version of thefluid heuristic, that a family of importance sampling distributions
based on a simple mixture appears suitable to design provably
efficient estimators. The family suggested by (38) is a mixtureof two components, one involving the occurrence of the rare
event in question and a second component corresponding to thenominal/original distribution of the increment distribution. The
form of the mixture is consistent with the way in which rare eventsoccur in a heavy-tailed setting. This type of mixture family was
introduced by [113] in the context of estimating P (Sm > n) as
n → ∞ for fixed m. In this context they were able to design astrongly efficient estimator for P (Sm > n). In fact, their estimator
was shown to give a relative error of at most ε > 0 for any fixed ε >0. Their proof technique is based on a weak convergence analysis,using Proposition 1, which is facilitated substantially because mis held fixed and independent of n. We now quickly illustrate theapproach proposed by [113].

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

Example 7 (Weak Convergence Approach And State-Dependent
Importance Sampling). We shall estimate αn = P (Sm > n)efficiently as n → ∞ for fixed m. For simplicity we shall assume
that the Xk ’s are non-negative random variables with densityf (·) and tail distribution function F (·). The family of importancesampling distributions that we consider is the following. Given that
Sk−1 = s for 1 ≤ k < m, then the next increment is sampledaccording to the density
fk (x|s) =pk (s) f (x) I (x > a (n − s))(1 − pk (s)) f (x) I (x ≤ a (n − s)),F (a (n − s))(40)

P (Xk > n − s)≈

P (Sm > n|Sk−1 = s)1
m−k+1

P (Xk > n − s)(m − k + 1)F (n − s).

Rn (X1 , . . . , Xm )f (X1 )

×f (X2 )

Rn I (υn = i) ≤ F a(1 − a)m n 2m I (υn = i)



× ···

f1 (X1 |S0 )f2 (X2 |S1 )f ( Xm )×
I (Sm > n)fm (Xm |Sm−1 )∏ F (a (n − Sk−1 )) I (Xk > a (n − Sk−1 ))=pk (Sk−1 )k=1F (a (n − Sk−1 )) I (Xk ≤ a (n − Sk−1 ))+(1 − pk (Sk−1 ))× I (Sm−1 < n)F (a (n − Sm−1 )) + I (Sm−1 ≥ n) .
Let P (·) be the probability measure induced by the proposed
importance sampling distribution; we use E (·) to denote theassociated expectation operator. Then1
P (Sm > n)21

E [Rn (X1 , . . . , Xm ) ]P (Sm > n)



1

∏

min(1/(m − k + 1), (1 − 1/(m − k + 1)))
k=1.

Therefore, we have thatRn (X1 , . . . , Xm )P (Sm > n)≤

F (a(1 − a)m n)
P (Sm > n)

m−1

×

2

∏

min(1/(m − k + 1), (1 − 1/(m − k + 1)))
k=1

= O(1)

as n → ∞. We then conclude thanks to Proposition 1 that1
P (Sm > n)=

E [Rn (X1 , . . . , Xm )|Sm > n]m−
E [Rn (X1 , . . . , Xm )|Xj > n]mP (Sm > n)

j =1

(1 + o(1))

as n → ∞. Evaluating the previous expectation we conclude thatmP (Sm > n)

j =1

→

a−α
m2



m
−

1

j =1

1/(m − j + 1)

×

j −1
∏

1

k=1

(1 − 1/(m − k + 1))= a−α .

m−1=

Now let υn = inf{k ≥ 1 : Xk > a(n − Sk−1 )} and note that on
υn = i we have that Si−1 ≤ n(1 − (1 − a)m )m−
E [Rn (X1 , . . . , Xm )|Xj > n]

Motivated by this observation, we shall propose pk (s) = 1/(m −k + 1) for 1 ≤ k ≤ m − 1 if s ≤ n and pk (s) = 1 is s > n. We thenconclude that the importance sampling estimator, Rn , satisfies=≤ an + a(1 − a)n + a(1 − a)2 n + · · · + a(1 − a)k−1 n≤ n(1 − (1 − a)m ) < n.

m−1

pk (s) ≈ P (Sk > n|Sm > n, Sk−1 = s)≈Sk ≤ an + (1 − a)[an + (1 − a)Sk−2 ]×where pk (s) ∈ (0, 1) is selected appropriately and a ∈ (0, 1). The
parameter a, as we shall see, is introduced in order to facilitate adominated convergence argument. If k = m (i.e. Sm−1 = s), thenthe increment is sampled according to the law of Xm given that
Xm > n − s. Dupuis et al. [113] obtains a limiting control problem
that allows to optimally select the pk ’s. Note that pk (s) represents
the probability that the rare event is caused in the k-th step giventhat Sk−1 = s < n. Therefore, if k ≤ m − 1, we would like to select=which iterating yields

F (a (n − s))+55

We can select a ∈ (0, 1) arbitrarily close to one in order to get
an asymptotically negligible relative error. Note however, that it is
crucial to select a ∈ (0, 1) in order to ensure that Rn /P (Sm > n)remains bounded. This suggests that a ≈ 1 might not be a good
option for a pre-determined n.
As it can be appreciated in the previous development, thefact that m is bounded substantially simplifies the use of the
weak convergence approach. An alternative technique, which
can be directly applied to Examples 7 and 8 and leverages thefluid heuristic is the use of Lyapunov inequalities, which wereintroduced in [9] and we discuss in the next lemma. 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 satisfying∫r −1 (x, y)P (Sj ∈ dy|Sj−1 = x) = 1

2

E [Rn (X1 , . . . , Xm )|Sm > n].

We wish to use Proposition 1 to approximate the previous
expectation. In order to do this we first will argue that Rn /P (Sm >n) is bounded by some constant uniformly as n → ∞. Note that
there exists 1 ≤ k ≤ m such that Xk > a(n − Sk−1 ). Otherwise, if
for all 1 ≤ k ≤ m, we have that Sk − Sk−1 = Xk ≤ a(n − Sk−1 ) or
equivalently we have thatSk ≤ an + (1 − a)Sk−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
TA ∧TBR=∏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)(41)

for x ̸∈ A ∪ B and G(x) ≥ ρ for x ∈ A. Then, E [R|S0 = x] ≤ G(x)/ρ.

56

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

The proof of the previous result is completely analogous toLemma 1 and therefore is omitted. We now explain how toapply the previous result, combined with the fluid heuristics todesign efficient importance sampling estimators in the context of
Example 5.
Example 8 (Lyapunov Inequalities And State-Dependent ImportanceSampling). We shall apply importance sampling following theanalysis in the previous section. We distinguish two situations,first, if the current position s is sufficiently far away from theboundary n, say if n − s > ∆ for some ∆ > 0, then the intuition
developed in the previous section 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 importance sampling. Also, if we use theprevious family when s is sufficiently close to n, we might create
variance because the mixture sampling distribution is derived
using the large deviations regime, so turning off the sampling whens 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|Sk = s)
P (Xk+1 ∈ dy) I (y > a(n − s))=p(n − s)I (n − s > ∆)P (Xk+1 > a(n − s))
P (Xk+1 ∈ dy) I (y ≤ a(n − s))+
P (Xk+1 ≤ a(n − s))× [1 − p(n − s)]I (n − s > ∆)+ P (Xk+1 ∈ dy) I (n − s ≤ ∆).
The parameter a ∈ (0, 1) is selected in order to enforce uniform
integrability conditions, just as we noted in Example 7. Following
the intuition developed earlier we propose

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 our current position s is suchthat G(n − s) < 1.
In order to proceed with the verification of inequality (41) wehenceforth assume that G(n − s) < 1. In this case inequality (41)takes the form
J1 + J2 ≤ 1 ,H ( n − s) =∫ ∞P (X > t )dt .

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, andusing the form of p(n − s), we obtain that
J1 ≤P (X > a(n − s))J2 ≤θ κ H (n − s)

 ∫ ∞



P (X > t )dt , 1

is satisfied, we would conclude strong efficiency of the estimator.
Note that the inequality (41) holds trivially if G(n − s) = 1.
Indeed, in this case, since G (x) ≤ 1 for every x ∈ R, the inequality(41) takes the form

.

∫ x+y

2H (t )Ḣ (t ) dt

x

= 2E [yH (x + yU )Ḣ (x + yU )],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 [X F (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 thatX

F (n − s − XU )H (n − s − XU )I (X ≤ a(n − s))H2(n − s)∼Xα−1( n − s)almost surely as n − s ↗ ∞ and also we have, because of regularvariation and Karamata’s theorem that there exists a constant K ∈(0, ∞) such that0

= min(κ EX + , 1) > 0.

H 2 (n − s)[1 − p(n − s)]

H 2 (x + y) − H 2 (x) =We now need to construct a valid Lyapunov function. We wish toprove strong efficiency, so it is natural to suggest

G(n − s) ≥ ρ := min κ

θ 2κ

E [H 2 (n − s − X ); X ≤ a(n − s)]n−s

for some κ > 0. The strategy is then to select (given a ∈ (0, 1))θ , κ and ∆ in order to satisfy the Lyapunov inequality. If we canshow that there exists a selection of these parameters that satisfies(41), then, given that the corresponding boundary condition in
Lemma 2, namelyp(n − s)Using the following Taylor development with remainder, whichcomes from the elementary formula using t = x + yu−

G(n − s) = min(κ H 2 (n − s), 1) = O(αn (s)2 )=Then we analyze the term J2 . Note that

θ P (X > a(n − s))p(n − s) =  ∞P (X > t )dt
n −s
for some θ > 0. In fact, we expect θ ≈ µ to be optimal in the senseof variance minimization.
For notational convenience define

(42)

where

|X |

F (n − s − XU )H (n − s − XU )I (X ≤ a(n − s))≤ |X |

H 2 ( n − s)
F ((n − s)(1 − a))H ((n − s) (1 − a))H2(n − s)≤K |X |n−s.

EG(n − s − X ) ≤ G(n − s) = 1.

Consequently, if κ (or equivalently ∆ = H −1 (1/κ 1/2 )) ischosen sufficiently large, we conclude, applied the dominatedconvergence theorem that for each ε > 0

Now observe that G(n − s) < 1 holds if and only if H (n − s) <1/κ 1/2 , which in turn holds if and only if n − s > H −1 (1/κ 1/2 ).

J2 ≤ 1 + 2EX(α − 1)(1 − ε)n−sJ. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59

as long as G(n − s) < 1. Combining our estimates for J1 and J2 andchoosing κ sufficiently large so that(α − 1)(1 + ε),n−s
(α − 1)(1 − ε)1
≥ 1 + θ a−α,
1 − p(n − s)n−sp(n − s) ≤ θ a−α

we arrive at(α − 1)(1 + ε) −α(α − 1)
a + 1 − 2µ
(1 − ε)J1 + J2 ≤θ κ(n − s)n−s(α − 1)(1 − ε)× 1 + θ a−αn−s(α − 1)(1 + ε)
(α − 1)+ 1 − 2µ(1 − ε)≤ a−αθ κ(n − s)n−s(α − 1)(1 − ε).
+ θ a−αn−sIn order to enforce (42) we then need to ensure that
a−α

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 selectthe parameters appropriately. This makes sense given that whenEX ≥ 0, αn (y) = 1 and it would be therefore impossible to boundthe second moment of the estimator by a function that decreases
to zero as n → ∞, as we are proposing.
Note that since ε > 0 is arbitrarily small, κ can be selected
arbitrarily large and a can be chosen arbitrarily close to one, notethat the feasible region for parameter selection for θ becomesθ < 2µ.
It turns out that by selecting θ = µ we can make the relative
error arbitrarily small. Selecting θ ≈ 0 implies a low frequencyof the jumps and therefore an increase in the termination time of
the algorithm. In fact, if θ is chosen too low, one can obtain (justas we saw in (39)) an infinite expected termination time. One canopt, however, to increase θ , as long as θ < 2µ in order to increase
the frequency of the big jumps and therefore make the termination
time of the algorithm shorter. Considering for instance an extreme
case, we saw that if α ∈ (1, 2) then the optimal change-of-measure(in the sense of variance minimization) has infinite expected termination time (that is, whenever α ∈ (1, 2), E (τn |τn < ∞) = ∞).
By choosing θ ≈ 2µ, one can show that the termination time of
the importance sampling estimator constructed by means of theprevious Lyapunov procedure, which maintains strong efficiency,
is finite as long as α ∈ (3/2, 2). These optimality results are discussed in [111].
5.5. Notes on importance sampling with heavy tails
Detailed expositions on large deviations asymptotics for heavytailed random variables can be found in [96,97] and Chapter IX of[40]. General discussion on importance samplers with heavy tails
is covered in Chapter VI of [23,50].
As mentioned earlier this section, the various challenges ofrare-event simulation with heavy tails are pointed out in [4,110].
The former focuses on the singularity of conditional measurewith respect to the original measure and the difficulty of finding
suitable importance sampler, while suggesting provably efficientconditional Monte Carlo scheme. The latter suggests the nonexistence of state-independent importance samplers for some firstpassage problems. Conditional Monte Carlo is also studied in [98]57

in the context of insurance, and further improved in [99]. It is also
used in [114] for sums of dependent random variables. The hazard
rate twisting method that we discussed in Section 5.1 appearsin [115,101,116], as well as [117] in the context of finance and [118]in single-server queues. Another version of this method is studiedin [119] using cross-entropy, while [120] suggests a closely relatedtransform method that essentially equates hazard rate twisting.
Finally, [121] considers a different type of twisting, by scaling
the variance of the underlying variables, for sums of correlated
lognormals.
Regarding state-dependent importance samplers, as we havediscussed, [113] suggests the weak convergence approach forregularly varying sums and the use of a two-mixture sampler
for each increment distribution. Two mixtures are no longer
enough for other heavy tails, such as Weibull. The appropriate
family of mixtures in such settings is developed in [111]. The
Lyapunov approach is proposed and applied in a variety of contexts
by [11] in single-server queues, [10,122] in regularly varying
random walks, [9] in random walk maximum and [123] in randomsum. The last two papers in particular aim to construct samplers
under the assumptions of general subexponential distributions.
The paper [124] further extends the method to the case of Markov
modulation. Mixture-type algorithm is also used in [125] forcomputing quantile and other risk measures that involve regularly
varying tails. Recently, [126] proposes sequential importancesampling technique with resampling that appears to exhibit
good empirical performance in terms of variance reduction.
Moreover, [127] considers an adaptive mechanism based on crossentropy methods to guide the selection of provably efficientimportance sampling estimators for heavy-tailed sums.
Acknowledgments
We are grateful to the referees for their useful comments which
helped to improve the exposition in this paper. In addition, support
from the NSF through the grants DMS-0806145, DMS-0846816 andDMS-1069064 is gratefully acknowledged.
References[1] P. Glasserman, S. Kou, Analysis of an importance sampling estimator for
tandem queues, ACM TOMACS 5 (1995) 22–42.
[2] P. Glasserman, Y. Wang, Counter examples in importance sampling for large
deviations probabilities, Annals of Applied Probability 7 (1997) 731–746.
[3] A. Bassamboo, S. Juneja, A. Zeevi, On the inefficiency of state-independent
importance sampling in the presence of heavy-tails, Operations Research
Letters 34 (2006) 521–531.
[4] S. Asmussen, K. Binswanger, B. Hojgaard, Rare events simulation for heavytailed distributions, Bernoulli 6 (2000) 303–322.
[5] P. Dupuis, H. Wang, Importance sampling, large deviations, and differential
games, Stochastic and Stochastic Reports 76 (2004) 481–508.
[6] P. Dupuis, H. Wang, Subsolutions of an Isaacs equation and efficient schemesof importance sampling, Mathematics of Operations Research 32 (2007)
723–757.
[7] P. Dupuis, A. Sezer, H. Wang, Dynamic importance sampling for queueing
networks, The Annals of Applied Probability 17 (2007) 1306–1346.
[8] P. Dupuis, K. Leder, H. Wang, Large deviations and importance sampling fora tandem network with slowdown, QUESTA 57 (2007) 71–83.
[9] J. Blanchet, P. Glynn, Efficient rare-event simulation for the maximum ofa heavy-tailed random walk, The Annals of Applied Probability 18 (2008)
1351–1378.
[10] J. Blanchet, J.C. Liu, State-dependent importance sampling for regularly varying random walks, Advances in Applied Probability 40 (2008) 1104–1128.
[11] J. Blanchet, P. Glynn, J.C. Liu, Fluid heuristics, Lyapunov bounds and efficientimportance sampling for a heavy-tailed G/G/1 queue, QUESTA 57 (2007)
99–113.
[12] J. Blanchet, Efficient importance sampling for binary contingency tables, The
Annals of Applied Probability 19 (2009) 949–982.
[13] S. Juneja, Estimating tail probabilities of heavy tailed distributions withasymptotically zero relative error, QUESTA 57 (2007) 115–127.
[14] P. L’Ecuyer, M. Mandjes, B. Tuffin, Importance sampling in rare event
simulation, in: G. Rubino, B. Tuffin (Eds.), Rare Event Simulation Using
Monte Carlo Methods, Wiley, West Sussex, United Kingdom, 2009, pp. 17–38.
(chapter 2).

58

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59[15] P. L’Ecuyer, J. Blanchet, B. Tuffin, P. Glynn, Asymptotic robustness of
estimators in rare-event simulation, ACM TOMACS 20 (2010) 6:1–6:41.
[16] B. Tuffin, Bounded normal approximation for the simulation of highly reliable
Markovian systems, Journal of Applied Probability 36 (1999) 974–986.
[17] B. Tuffin, On numerical problems in simulation of highly reliable Markovian
systems, in: Proceedings of the 1st International Conference on Quantitative
Evaluation of SysTems (QEST), IEEE Computer Society Press, Los Alamitos, CA,
2004, pp. 156–164.
[18] L.W. Schruben, A coverage function for interval estimators of simulation
response, Management Science 26 (1980) 18–27.
[19] H. Cancela, G. Rubino, B. Tuffin, New measures of robustness in rare event
simulation, in: Proceedings of the 2005 Winter Simulation Conference, IEEE
Press, 2005, pp. 519–527.
[20] P. Glynn, W. Whitt, The asymptotic efficiency of simulation estimators,
Operations Research 40 (1992) 505–520.
[21] P. Glynn, G. Rubino, B. Tuffin, Robustness properties and confidence interval
reliability issues, in: Rare Event Simulation Using Monte Carlo Methods, West
Sussex, United Kingdom, 2009, pp. 63–84. (chapter 4).
[22] R.Y. Rubinstein, D.P. Kroese, Simulation and the Monte Carlo Method, second
ed. John Wiley, 2009.
[23] S. Asmussen, P. Glynn, Stochastic Simulation: Algorithms and Analysis,
Springer-Verlag, New York, NY, USA, 2007.
[24] J. Liu, Monte Carlo Strategies in Scientific Computing, Springer, 2001.
[25] A. Dembo, O. Zeitouni, Large Deviations Techniques and Applications, second
ed. Springer, New York, 1998.
[26] A. Ganesh, N. O’Connell, D. Wischik, Big Queues, Springer-Verlag, Berlin
Heidelberg, Germany, 2004.
[27] T.M. Cover, J.A. Thomas, Elements of Information Theory, second ed. Wiley,
2006.
[28] S. Asmussen, R. Rubinstein, Steady state rare events simulation in queueing
models and its complexity properties, in: J. Dshalalow (Ed.), Advances inQueueing, CRC Press, 1995, pp. 429–461.
[29] D. Siegmund, Importance sampling in the Monte Carlo study of sequential
tests, The Annals of Statistics 3 (1976) 673–684.
[30] S. Asmussen, Applied Probability and Queues, Springer-Verlag, New York, NY,USA.
[31] R.T. Rockafellar, Convex Analysis, Princeton University Press, 1970.
[32] P. Glasserman, Monte Carlo Methods in Financial Engineering, SpringerVerlag, New York, 2004.
[33] A. Shwartz, A. Weiss, Large Deviations for Performance Analysis. Queues,
Communications, and Computing, Chapman & Hall, London, UK, 1995.
[34] R.S. Ellis, Entropy, Large Deviations and Statistical Mechanics, SpringerVerlag, Berlin Heidelberg, Germany, 1985.
[35] P. Dupuis, R.S. Ellis, A Weak Convergence Approach to the Theory of Large
Deviations, Wiley, New York, 1997.
[36] J.S. Sadowsky, Large deviations theory and efficient simulation of excessive
backlogs in a GI /GI /m queue, IEEE Transactions on Automatic Control 36(1991) 1383–1394.
[37] P. Glynn, D. Iglehart, Importance sampling for stochastic simulations, IEEE
Transactions on Automatic Control 35 (1989) 1367–1392.
[38] S. Asmussen, Conjugate processes and the simulation of ruin problems,
Stochastic Processes and their Applications 20 (1985) 213–229.
[39] T. Lehtonen, H. Nyrhinen, Simulating level-crossing probabilities by importance sampling, Advances in Applied Probability 24 (1992) 858–874.
[40] S. Asmussen, Ruin Probabilities, World Scientific, London, UK, 2000.
[41] J. Collamore, Importance sampling techniques for the multidimensional ruin
problem for general markov additive sequences of random vectors, The
Annals of Applied Probability 12 (2002) 382–421.
[42] D.P. Kroese, V.F. Nicola, Efficient estimation of overflow probabilities inqueues with breakdowns, Performance Evaluation 36–37 (1999) 471–484.
[43] S. Parekh, J. Walrand, Quick simulation of rare events in networks, IEEE
Transactions on Automatic Control 34 (1989) 54–66.
[44] C.S. Chang, P. Heidelberger, S. Juneja, P. Shahabuddin, Effective bandwidthand fast simulation of ATM intree networks, Performance Evaluation 20(1994) 45–65.
[45] P.T. de Boer, Analysis of state-independent importance-sampling measures
for the two-node tandem queue, ACM Transactions on Modeling andComputer Simulation 16 (2006) 225–250.
[46] J.S. Sadowsky, On Monte Carlo estimation of large deviations probabilities,
The Annals of Applied Probability 6 (1996) 399–422.
[47] J. Bucklew, Introduction to Rare Event Simulation, Springer-Verlag, New
York, NY, 2004.
[48] A.B. Dieker, M. Mandjes, On asymptotically efficient simulation of large
deviation probabilities, Advances in Applied Probability 37 (2005) 539–552.
[49] J. Bucklew, Large Deviation Techniques in Decision, Simulation andEstimation, Wiley, New York, NY, USA, 1990.
[50] S. Juneja, P. Shahabuddin, Rare event simulation techniques: an introduction
and recent advances, in: S.G. Henderson, B.L. Nelson (Eds.), Simulation,
Handbooks in Operations Research and Management Science, Elsevier,
Amsterdam, The Netherlands, 2006.
[51] P. Heidelberger, Fast simulation of rare events in queueing and reliability
models, ACM Transactions on Modeling and Computer Simulation (TOMACS)
5 (1995) 43–85.
[52] J. Blanchet, M. Mandjes, Rare event simulation for queues. in: Rare Event
Simulation Using Monte Carlo Methods, West Sussex, United Kingdom, 2009,
pp. 87–124. (chapter 5).

[53] R. Szechtman, P. Glynn, Rare event simulation for infinite server queues. in:
E. Yucesan, C. H. Chen, J. L. Snowdon, and J. M. Charnes, Eds., Proceedings ofthe 2002 Winter Simulation Conference, pp. 416–423, 2002.
[54] A. Ridder, Importance sampling algorithms for first passage time probabilities in the infinite server queue, European Journal of Operational Research
199 (2009) 176–186.
[55] J. Blanchet, P. Glynn, H. Lam, Rare event simulation for a slotted time M /G/s
model, Queueing System 63 (2009) 33–57.
[56] G. Rubino, B. Tuffin, Markovian models for dependability analysis, in: G. Rubino, B. Tuffin (Eds.), Rare Event Simulation Using Monte Carlo Methods, Wiley, West Sussex, United Kingdom, 2009, pp. 125–144. (chapter 6).
[57] P. Grassberger, Go with the winners: a general Monte Carlo strategy,
Computer Physics Communications 147 (2002) 64–70.
[58] W. Sandmann, Rare event simulation methodologies in systems biology,
in: Rare Event Simulation Using Monte Carlo Methods, Wiley, West Sussex,
United Kingdom, 2009, pp. 243–266. (chapter 11).
[59] R.Y. Rubinstein, D.P. Kroese, The Cross-Entropy Method, Springer, New York,NY, 2004.
[60] M. Bayati, J. Kim, A. Saberi, A sequential algorithm for generating randomgraphs, in: Approximation, Randomization and combinatorial Optimization.
Algorithms and Techniques, Lecture Notes in Computer Science 4627 (2007)
326–340.
[61] P. L’Ecuyer, F. Le Gland, P. Lezaud, B. Tuffin, Splitting techniques, in: G. Rubino,
B. Tuffin (Eds.), Rare Event Simulation Using Monte Carlo Methods, Wiley,West Sussex, United Kingdom, 2009, pp. 39–62. (chapter 3).
[62] P. Glasserman, P. Heidelberger, P. Shahabuddin, T. Zajic, A large deviations
perspective on the effiency of multilevel splitting, IEEE Transactions on
Automated Control (1998) 1666–1679.
[63] P. Glasserman, P. Heidelberger, P. Shahabuddin, T. Zajic, Multilevel splitting
for estimating rare event probabilities, Operations Research 47 (1999)
585–600.
[64] T. Dean, P. Dupuis, Splitting for rare event simulation: a large deviation
approach to design and analysis, Stochastic Processes and their Applications
119 (2009) 562–587.
[65] R.Y. Rubinstein, Randomized algorithms with splitting: why the classic
randomized algorithms do not work and how to make them work,
Methodology and Computing in Applied Probability 12 (2010) 1–50.
[66] P. del Moral, Feynman–Kac Formulae: Genealogical and Interacting Particle
Systems with Applications, Springer-Verlag, New York, USA, 2004.
[67] C. Kollman, K. Baggerly, D. Cox, R. Picard, Adaptive importance sampling on
discrete Markov chains, Annals of Applied Probability 9 (1999) 391–412.
[68] P.Y. Desai, P.W. Glynn, A Markov chain perspective on adaptive Monte Carloalgorithms, Proceedings of 2001 Winter Simulation Conference 9 (2001)
391–412.
[69] P.T. De Boer, V.F. Nicola, R.Y. Rubinstein, Adaptive importance sampling
simulation of queueing networks, in: K. Kang, J.A. Joines, R.R. Barton,
P.A. Fishwick (Eds.), Proceedings of 2000 Winter Simulation Conference, IEEE
Press, Piscataway, New Jersey, 2000, pp. 646–655.
[70] P.T. de Boer, D. Kroese, S. Mannor, R. Rubinstein, A tutorial on the crossentropy method, Annals of Operations Research 134 (2005) 19–67.
[71] R.Y. Rubinstein, Optimization of computer simulation models with rare
events, European Journal of Operations Research 99 (1997) 89–112.
[72] R.Y. Rubinstein, Rare-event simulation via cross-entropy and importancesampling. Second Workshop on Rare Event Simulation, RESIM 99, pp. 1–17,
1999.
[73] T.P.I. Ahamed, V.S. Borkar, S. Juneja, Adaptive importance sampling techniquefor Markov chains using stochastic approximation, Operations Research 54(2006) 489–504.
[74] V.S. Borkar, S. Juneja, A.A. Kherani, Performance analysis conditioned on rare
events: an adaptive simulation scheme, Communications in Information 3(2004) 259–278.
[75] H.P. Chan, T.L. Lai, A sequential Monte Carlo approach to computing tailprobabilities in stochastic models, The Annals of Applied Probability, 2011
(forthcoming).
[76] E. Weinan, R. Weiqing, E. Vanden-Eijnden, Minimum action method for the
study of rare events, Communications on Pure & Applied Mathematics 52(2004) 637–656.
[77] H. Hult, F. Lindskog, Heavy-tailed insurance portfolios: buffer capital and ruinprobabilities, School of ORIE, Cornell University, Technical Report No. 1441,
2006.
[78] J. Blanchet, P. Glynn, K. Leder, On Lyapunov inequalities and subsolutions for
efficient importance sampling, 2010. Preprint,
http://www.columbia.edu/jb2814/papers/LyapSub7_Dec20_2009.pdf.
[79] P. Dupuis, H. Wang, Dynamic importance sampling for uniformly recurrent
Markov chains, The Annals of Applied Probability 15 (2005) 1–38.
[80] P. Dupuis, H. Wang, Importance sampling for Jackson networks, Queueing
Systems 62 (2009) 113–157.
[81] P. Dupuis, H. Ishii, H.M. Soner, A viscosity solution approach to the asymptotic
analysis of queueing systems, The Annals of Probability 18 (1990) 226–255.
[82] J. Blanchet, H. Lam, Rare-event simulation for many-server queues, Preprint.
[83] J.S. Sadowsky, J.A. Bucklew, On large deviations theory and asymptoticallyefficient Monte Carlo estimation, IEEE Transactions on Information Theory
36 (1990) 579–588.

J. Blanchet, H. Lam / Surveys in Operations Research and Management Science 17 (2012) 38–59[84] H.P. Chan, T.L. Lai, Efficient importance sampling for Monte Carlo evaluationof exceedance probabilities, The Annals of Applied Probability 17 (2007)
440–473.
[85] Y.V. Linnik, On the probability of large deviations for the sums of independent
variables, in: Proceedings of 4th Berkeley Symposium on Mathematical
Statistics and Probability, vol. 1, 1961, pp. 283–306.
[86] A.V. Nagaev, Integral limit theorems taking large deviations into accountwhen Cramer’s condition does not hold II, Theory of Probability and Its
Applications 14 (1969) 193–208.
[87] S. Nagaev, Large deviations of sums of independent random variables, Annalsof Probability 7 (1979) 745–789.
[88] L.V. Rozovskii, Probabilities of large deviations of sums of independent
random variables, Theory of Probability and Its Applications 34 (1990)
625–644.
[89] L.V. Rozovskii, Probabilities of large deviations on the whole axis, Theory ofProbability and its Applications 38 (1993) 53–79.
[90] A.A. Borovkov, K.A. Borovkov, On probabilities of large deviations for random
walks. I. Regularly varying distribution tails, Theory of Probability and its
Applications 49 (2001) 189–205.
[91] S. Asmussen, Subexponential asymptotics for stochastic processes: extremal
behavior, stationary distributions and first passage probabilities, The Annalsof Applied Probability 8 (1998) 354–374.
[92] S. Foss, Z. Palmowski, S. Zachary, The probability of exceeding a highboundary on a random time interval for a heavy-tailed random walk, Annalsof Applied Probability 15 (2005) 1936–1957.
[93] S. Foss, T. Konstantopoulos, S. Zachary, Discrete and continuous time modulated random walks with heavy-tailed increments, Journal of Theoretical
Probability 20 (2007) 581–612.
[94] S. Foss, S. Zachary, The maximum on a random time interval of a random walkwith long-tailed increments and negative drift, Annals of Applied Probability
13 (2003) 37–53.
[95] S. Asmussen, C. Kluppelberg, Large deviations results for subexponential
tails, with applications to insurance risk, Stochastic Processes and their
Applications 64 (1996) 103–125.
[96] P. Embrechts, C. Klpelberg, T. Mikosch, Modelling Extremal Events for
Insurance and Finance, Springer-Verlag, New York, 1997.
[97] S.I. Resnick, Extreme Values, Regular Variation, and Point Processes, New
York, 1987.
[98] S. Asmussen, K. Binswanger, Simulation of ruin probabilities for subexponential claims, ASTIN Bulletin 27 (1997) 297–318.
[99] S. Asmussen, D. Kroese, Improved algorithms for rare event simulation withheavy tails, Advances in Applied Probability 38 (2006) 545–558.
[100] L. Rojas-Nandayapa, S. Asmussen, Efficient simulation of finite horizon
problems in queueing and insurance risk, QUESTA 57 (2007) 85–97.
[101] S. Juneja, P. Shahabuddin, Simulating heavy-tailed processes using delayed
hazard rate twisting, ACM TOMACS 12 (2002) 94–118.
[102] H. Hult, F. Lindskog, T. Mikosch, G. Samordnitsky, Functional large deviations
for multivariate regularly varying random walks, The Annals of Applied
Probability 15 (2005) 2651–2680.
[103] A.P. Zwart, Queueing Systems with Heavy Tails, Ph.D. Dissertation, The
Netherlands, 2001.
[104] B. Zwart, S. Borst, M. Mandjes, Exact asymptotics for fluid queues fed bymultiple heavy-tailed on-off flows, Annals of Applied Probability 14 (2004)
903–957.
[105] S. Foss, D. Korshunov, Heavy tails in multi-server queue, QUESTA 52 (2006)
31–48.
[106] A.G. Pakes, On the tails of waiting-time distribution, Journal of Applied
Probability 12 (1975) 555–564.
[107] N. Veraverbeke, Asymptotic behavior of Wiener–Hopf factors of a random
walk, Stochastic Processes And Their Applications 5 (1977) 27–31.

59[108] P. Embrechts, N. Veraverbeke, Estimates for the probability of ruin withspecial emphasis on the possibility of large claims, Insurance: Mathematical
Economics 1 (1982) 55–72.
[109] J.W. Cohen, On the M /G/1 queue with heavy-tailed service time distributions, Technical Report PNA-R9702, CWI, Amsterdam, 1997.
[110] A. Bassamboo, S. Juneja, A. Zeevi, Importance sampling simulation in thepresence of heavy tails, in: M.E. Kuhl, N.M. Steiger, F.B. Armstrong, J.A. Joines(Eds.), Proceedings of 2005 Winter Simulation Conference, IEEE Press,
Piscataway, New Jersey, 2005, pp. 664–672.
[111] J. Blanchet, J.C. Liu, Efficient simulation and conditional functional limit
theorems for ruinous heavy-tailed random walks, 2010.
http://www.columbia.edu/jb2814/papers/ConcaveHazard9f.pdf (submitted
for publication).
[112] C. Kluppelberg, A. Kyprianou, R. Maller, Ruin probabilities and over-shootsfor general levy insurance risk processes, The Annals of Applied Probability
14 (2004) 1766–1801.
[113] P. Dupuis, K. Leder, H. Wang, Importance sampling for sums of randomvariables with regularly varying tails, ACM TOMACS 17 (2006).
[114] J. Blanchet, L. Rojas-Nandayapa, Efficient simulation of tail probabilities ofsums of dependent random variables, Journal of Applied Probability, Special
Vol. 48A (2011) 147–165.
[115] S. Juneja, P. Shahabuddin, A. Chandra, Simulating heavy-tailed processes
using delayed hazard-rate twisting, in: P.A. Farrington, H.B. Nembhard,
D.T. Sturrock, G.W. Evans (Eds.), Proceedings of 1999 Winter Simulation
Conference, IEEE Press, Piscataway, New Jersey, 1999, pp. 420–427.
[116] Z. Huang, P. Shahabuddin, A unified approach for finite dimensional, rareevent Monte Carlo simulation, in: R.G. Ignalls, M.D. Rossetti, J.S. Smith,
B.A. Peters (Eds.), Proceedings of 2004 Winter Simulation Conference, IEEE
Press, Piscataway, New Jersey, 2004, pp. 1616–1624.
[117] Z. Huang, P. Shahabuddin, Rare-event, heavy-tailed simulations using hazard function transformations, with applications to value-at-risk, in: S. Chick,
P.J. Sanchez, D. Ferrin, D. Morrice (Eds.), Proceedings of 2003 Winter Simulation Conference, IEEE Press, Piscataway, New Jersey, 2003, pp. 276–284.
[118] N.K. Boots, P. Shahabuddin, Simulating GI /GI /1 queues and insurance
processes with subexponential distributions, in: J.A. Joines, R.R. Barton,
K. Kang, P.A. Fishwick (Eds.), Proceedings of 2000 Winter Simulation
Conference, IEEE Press, Piscataway, New Jersey, 2000, pp. 656–665.
[119] S. Asmussen, D.P. Kroese, R.Y. Rubinstein, Heavy tails, importance sampling
and cross entropy, Stochastic Models 21 (2005) 57–76.
[120] D.P. Kroese, R.Y. Rubinstein, The transform likelihood ratio method for rareevent simulation with heavy tails, Queueing Systems 46 (2004) 317–351.
[121] S. Asmussen, J. Blanchet, S. Juneja, L.R. Nandayapa, Efficient simulation of tailprobabilities of sums of correlated lognormals, Annals of Operations Research
189 (2009) 5–23.
[122] J. Blanchet, J.C. Liu, Importance sampling in ruin problems for multidimensional regularly varying random walks, Journal of Applied Probability 47(2010) 301–322.
[123] J. Blanchet, C. Li, Efficient rare-event simulation for heavy-tailed compoundsums, ACM TOMACS Transactions in Modeling and Computer Simulation 21(2011) 1–10.
[124] J. Blanchet, H. Lam, Efficient state-dependent importance sampling for
Markov-modulated heavy-tailed random walk, Working Paper.
[125] H. Hult, J. Svensson, Efficient calculation of risk measures by importancesampling—the heavy tailed case. arXiv:0909.3335v1 [math.PR], 2009.
[126] H.P. Chan, S. Deng, T.L. Lai, Rare-event simulation of heavy-tailed random
walks by sequential importance sampling and resampling, Preprint.
[127] J. Blanchet, Y. Shi, Efficient rare event simulation for heavy-tailed systems
via cross entropy. S. Jain, R.R. Creasey, J. Himmelspach, K. P. White, and M.
Fu, (Eds.), in: Proceedings of 2011 Winter Simulation Conference. IEEE Press(in press).