Proceedings of the 2022 Winter Simulation Conference
B. Feng, G. Pedrielli, Y. Peng, S. Shashaani, E. Song, C.G. Corlu, L.H. Lee, E.P. Chew, T. Roeder, andP. Lendermann, eds.
IMPORTANCE SAMPLING FOR RARE-EVENT GRADIENT ESTIMATION
Yuanlu Bai
Shengyi He
Henry Lam

Guangxin Jiang
School of Management
Harbin Institute of Technology
No. 92 Xidazhi Street
Harbin, Heilongjiang 150001, CHINA

Department of Industrial Engineering& Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027, USA

Michael C. Fu
Robert H. Smith School of Business andInstitute for Systems Research
University of MarylandCollege Park, MD 20742, USA

ABSTRACT
Importance sampling (IS) is a powerful tool for rare-event estimation. However, in many settings, we needto estimate not only the performance expectation but also its gradient. In this paper, we build a bridgefrom the IS for rare-event estimation to gradient estimation. We establish that, for a class of problems, anefficient IS sampler for estimating the probability of the underlying rare event is also efficient for estimating
gradients of expectations over the same rare-event set. We show that both the infinitesimal perturbation
analysis and the likelihood ratio estimators can be studied under the proposed framework. We use twonumerical examples to validate our findings.
1

INTRODUCTION

Gradient estimation for the performance of stochastic systems is an important topic in simulation analysis,
which naturally appears in many problems of interest. A prime example is to solve stochastic optimization
problems, such as in simulation-based optimization or empirical risk minimization in machine learning. In
these contexts, gradient estimation is used in each iteration of the gradient descent or stochastic approximation
algorithm. In addition to optimization, gradient estimation arises in model sensitivity analysis, hedging
financial derivatives, among other applications.
Many gradient estimation methods have been proposed (see Fu (2006) for a thorough introduction).
If we do not have any information about the system, i.e., the model is “black-box” or that we are ina “zeroth-order” setting, finite-difference approximation is commonly used to estimate gradients. Thisapproach perturbs the value of each component of the parameter of interest (see Glynn (1989), Fox and
Glynn (1989), and Section 5.2.1 in Fu (2006)). If we know some information about the simulation model
structure or the input probability distribution, then we can devise unbiased or “first-order” estimators. In
particular, infinitesimal perturbation analysis (IPA) or pathwise differentiation, and the likelihood ratio

Bai, He, Lam, Jiang, and Fu(LR) or the score function (SF) method, and their variants, have been commonly used. These methods
operate by interchanging the derivative and expectation operations. IPA, introduced by Ho and Cao (1983),
takes derivative of the sample performance directly, and requires the continuity of the sample performance(Glasserman (1988)). To handle discontinuities, smoothed perturbation analysis is proposed to smooth thesample performance by taking a conditional expectation (see Fu and Hu (1997)). The LR or SF method,
proposed by Reiman and Weiss (1989), Rubinstein (1989) and Glynn (1990), takes derivatives of parametersinside the input distribution and requires the continuity therein (L’Ecuyer (1990)). While the LR estimator
does not require smoothness of the sample performance, it typically have larger variances than those of
IPA. Improvements and complements of the IPA and the LR method include kernel methods (Liu andHong (2011)), weak derivatives (Pflug (1989); Heidergott et al. (2010)), and the generalized likelihood
ratio method (Wang et al. (2012); Peng et al. (2018)).
In this paper, we study gradient estimation where the performance is evaluated on a rare-event region.
This problem arises in optimization or sensitivity analysis over tail-related quantities, which appears whentraining risk-averse predictive models or conducting performance analyses on extremal risk measures (Heet al. (2022)). In these settings, estimation based on crude Monte Carlo (MC), even just for expectation
evaluations, would be inefficient, because a prohibitively large sample size is required to ensure the rare
event is observed sufficiently adequately. Statistically, this inefficiency manifests as a huge estimation
variance. To address this challenge, variance reduction techniques such as importance sampling (IS), are
known to be powerful in enhancing efficiency. The idea of IS is to draw samples from an alternative
distribution that hits the rare event more frequently and adjusts the estimator with a likelihood ratio toretain unbiasedness. When the alternative distribution is properly chosen, IS can significantly reduce thevariance of the estimation. Selecting a good alternative distribution has been the key of IS and has been
extensively studied (see, e.g., the surveys Bucklew (2004), Juneja and Shahabuddin (2006), Blanchet andLam (2012)).
Despite a vast literature on the application of IS to rare-event estimation problems, there are few
studies on using IS for gradient estimation problems that involve rare events which, as mentioned above,
nonetheless arises in risk-averse model training and performance analysis. In this paper, we study rare-event
gradient estimation by looking at a class of problems that generically arise in this context. More concretely,
we consider estimating the expectation of an objective function that is piecewise polynomial, evaluatedon a rare event that is defined as the excursion of a piecewise linear function applied on a random object.
Under suitable conditions, this objective function can result from both the IPA and the LR or SF methods
in estimating the gradients of rare-event probabilities or related risk quantities. We show that in the specificcase of Gaussian inputs, the IS based on Gaussian mixture designed for the estimation of the underlying
rare-event probability is also efficient for the estimation of the more general objective function. Thissuggests that once we obtain an efficient sampler for the rare-event probability, the same sampler could be
used for a wide class of gradient estimation problems on the same rare-event set. Our result thus facilitates
the design of rare-event gradient estimators by identifying a readily obtainable efficient IS. Although our
results only focus on Gaussian inputs, similar theories are conjectured to hold for more general light-tailed
distributions which we relegate to a full journal version of this work.
We briefly review other studies on variance reduction techniques for gradient estimation. The most
relevant work we have seen is Nakayama (1995), which studies the gradients on the performances of highly
reliable Markovian systems. It is shown that when applying a type of IS called “balanced failure biasing”,
all gradients can be estimated as accurately as the performance measure itself. Our study can be viewed
as a generalization of this observation. Besides Nakayama (1995), control variates have been applied instochastic gradient descent algorithms to reduce the variance of the gradient estimate and improve theconvergence rates of the algorithms (see Wang et al. (2013) and Driggs et al. (2022)). For empirical risk
minimization problems, stochastic variance reduced gradient (SVRG) was proposed, which progressively
reduces the variance using the averaged gradient at an estimated optimal solution. To save computationalcost, the estimated optimal solution is updated only periodically (Johnson and Zhang (2013); Xiao andBai, He, Lam, Jiang, and Fu
Zhang (2014)). Another work on variance reduction in empirical risk minimization is Zhao and Zhang(2015), which uses IS to select the data when running stochastic gradient descent. However, as mentioned
before, there are few works on IS to reduce the variance of the gradient estimator with rare events, eventhough IS has been extensively studied in estimating a rare-event probability or its related expectation.
The remainder of this paper is organized as follows: Section 2 introduces our problem class and the
proposed IS for gradient estimation, and shows the efficiency of the proposed IS. In Section 3, we use twonumerical examples to show the efficiency of the mixture IS in the estimation of both rare-event expectationand gradient. Finally, Section 4 provides concluding remarks.
2

IMPORTANCE SAMPLING FOR GRADIENT ESTIMATIONAs we will show momentarily, a class of gradients derived using either IPA or LR can be formulated as
E [ f (X)1{g(X) ≥ γ}] ,(1)where {g(X) ≥ γ} represents a rare event. We impose the following assumption regarding the underlying
distribution, the function f and the structure of the rare event.
Assumption 1 X ∼ N(µ, Σ) for some µ ∈ Rd and Σ ∈ Rd×d . f (·) is a piecewise polynomial and g(·)
is a piecewise linear function such that we can split {g(x) ≥ γ} into n mutually exclusive polyhedrons{A(i) x + b(i) ≥ t (i) }(i = 1, · · · , n) and f (x) = fi (x) over polyhedron {A(i) x + b(i) ≥ t (i) }, where f (i) (·) isa polynomial function. Moreover, for any i = 1, 2, . . . , n, {A(i) x + b(i) ≥ t (i) } is full-dimensional, i.e.,P({A(i) X + b(i) ≥ t (i) }) > 0 where X ∼ N(µ, Σ).
We show that a class of IPA and LR can be formulated as (1) such that Assumption 1 holds.
Example 1 (IPA) Suppose we are interested in the gradient of E [ f (X, θ )1{g(X) ≥ γ}] w.r.t. θ , where
X ∼ N(µ, Σ) and f (x, θ ) is a piecewise polynomial function of x whose coefficients are determined by θ .
Then we have that∂∂
E [ f (X, θ )1{g(X) ≥ γ}] = E
f (X, θ )1{g(X) ≥ γ} ,∂θ∂θwhich has the form of (1).
Example 2 (LR) Suppose we are interested in the gradient ∂∂θ Eθ [ f (X)1{g(X) ≥ γ}] where X ∼ N(µ(θ ), Σ).
Using LR and noting that∂∂log pθ (X) =µ(θ )> Σ−1 (X − µ(θ )),∂θ∂θwhere pθ (X) is the density of N(µ(θ ), Σ), we have



∂
∂
> −1
Eθ [ f (X)1{g(X) ≥ γ}] = Eθµ(θ ) Σ (X − µ(θ )) f (X)1{g(X) ≥ γ} .
∂θ
∂θ


If f is a piecewise polynomial, then the function ∂∂θ µ(θ )> Σ−1 (X − µ(θ )) f (X) is also a piecewise
polynomial. Therefore, the LR estimator can also be written as the form of (1).
For the estimation of the probability P(g(X) ≥ γ), an efficient sampler is given in Bai et al. (2022),
which is described as follows. For each i = 1, 2, . . . , n, let xi∗ denote the minimizer of (x − µ)> Σ−1 (x − µ)over the polyhedron {A(i) X + b(i) ≥ t (i) }. The proposed sampler is a Gaussian mixture whose density isp(x) =

1 n∑ pi (x)
n i=1

(2)

Bai, He, Lam, Jiang, and Fu
where pi (x) is the density of N (xi∗ , Σ). The associated likelihood ratio is given by
T −1 (x−µ)1

L(x) =

ne− 2 (x−µ) Σ1∗ T −1 (x−x∗ )1

e− 2 (x−x1 ) Σ

1∗ T −1 (x−x∗ )n+ · · · + e− 2 (x−xn ) Σ.

We propose to use the same sampler for the problem (1). Therefore, the proposed estimator for (1) isgiven byf (X̃)L(X̃)1{g(X̃) ≥ γ},
where X̃ ∼ p(·). We will show that the proposed estimator is efficient. Typically, for rare-event estimation
problems, the efficiency is defined relative to the mean of the estimator. However, since we are looking at
gradient estimation problems, usually we do not have f > 0 in (1), so it might happen that the mean isextremely small due to the cancellation between the positive parts of f and the negative parts of f . For
this reason, we use the rare-event probability as the benchmark. This benchmark is reasonable, although
its usefulness would depend on the nature of the downstream tasks in consideration. More precisely, we
have the following theorem.
Theorem 1 Suppose that Assumption 1 holds. Let X̃ be a sample drawn from p(·) as defined in (2) andX ∼ N(µ, Σ). Then we have that
E f (X̃)2 L2 (X̃)1{g(X̃) ≥ γ}
(3)
(P(g(X) ≥ γ))2
does not grow exponentially as γ → ∞.
Proof.
Without loss of generality, throughout the proof we suppose that µ = 0. Under Assumption 1,
the expectation in (1) can be written as
n

h
i
(i)
(i)
(i)
E
f
(X)1{A
X
+
b
≥
t
}
.
∑ i

i=1

For each i = 1, 2, . . . , n, let N(xi∗ ) = {kx − xi∗ k2 ≤ ε} denote a neighborhood of xi∗ . Since {A(i) x + b(i) ≥ t (i) }is full-dimensional, by introducing redundant constraints, we may find a full-rank B(i) and a vector s(i)such that {A(i) x + b(i) ≥ t (i) } ⊂ {B(i) x ≥ s(i) } and when ε is sufficiently small,N(xi∗ ) ∩ {A(i) x + b(i) ≥ t (i) } = N(xi∗ ) ∩ {B(i) x ≥ s(i) }.
See Lemma 7.2 and Lemma 7.1 of Bai et al. (2022) for more details on the construction of B(i) and s(i) .
From Lemma 7.2 of Bai et al. (2022), we have thatpoly

1 ∗> −1 ∗xi

P({A(i) X + b(i) ≥ t (i) }) ∼ e− 2 xi Σ

(4)

poly

where f1 (γ) ∼ f2 (γ) if f1 (γ)/ f2 (γ) changes at most polynomially in γ.
The second moment of the estimator can be written asn

E[( f (X̃))2 1{g(X̃) ≥ γ}L2 (X̃)] = E[( f (X))2 1{(g(X) ≥ γ}L(X)] = ∑ E[( fi (X))2 L(X)1{A(i) X + b(i) ≥ t (i) }].
i=1

For each i, we are going to give an upper bound of E[( fi (X))2 L(X)1{A(i) X + b(i)} ≥ t (i) ]. We knowthat

1 T −1x

L(x) ≤

ne− 2 x Σ1
∗ T −1∗e− 2 (x−xi ) Σ (x−xi )1 ∗> −1 ∗xi −xT Σ−1 xi∗= ne 2 xi Σ.

Bai, He, Lam, Jiang, and Fu
We omit the index i in A(i) , B(i) , b(i) , s(i) , xi∗ from now on for convenience. Let φΣ (x) denote the density ofN(0, Σ). Using the above inequality, we have that2

2

E[( fi (X)) L(X)1{AX + b ≥ t}] ≤ E[( fi (X)) L(X)1{BX ≥ s}] =≤Z

1 ∗> −1 ∗Σ x −xT Σ−1 x∗

ne 2 x{Bx≥s}

Z{Bx≥s}( fi (X))2 L(x)φΣ (x)dx( fi (x))2 φΣ (x)dx

 d2 Z
1∗ T −1∗∗ T −1 ∗T −1 ∗
1=n( fi (x))2 e− 2 (x−x ) Σ (x−x )−(x−x ) Σ x −x Σ x dx
2π|Σ|{Bx≥s}
 d2
Z1∗> −1 ∗
∗ T −1∗∗ T −1 ∗1e−x Σ x( fi (x))2 e− 2 (x−x ) Σ (x−x )−2(x−x ) Σ x dx
=n2π|Σ|{Bx≥s}
 d2
Z1∗> −1 ∗
∗ T −1∗1e−x Σ x( fi (x))2 e− 2 (x−x ) Σ (x−x ) dx≤n2π|Σ|{Bx≥s}
 d2
Z1∗> −1 ∗T11≤ne−x Σ x( fi (x∗ + B−1 y))2 e− 2 λmin y y dy.
∗|B| 2π|Σ|{y≥s−Bx }

Here we used the substitution x = x∗ + B−1 y in the last step and λmin is the minimum eigenvalue ofR1T(B−1 )> Σ−1 B−1 . Since fi is a polynomial, {y≥s−Bx∗ } ( fi (x∗ + B−1 y))2 e− 2 λmin y y dy can be bounded fromabove by a polynomial of absolute values of components of x∗ . Therefore, combining the above inequality
with (4), we get the desired result.
In Bai et al. (2022), the same sampler is shown to be relatively efficient for the estimation of theprobability of the rare event {g(X) ≥ γ}. Combining this with Theorem 1, we observe that in this case, anefficient sampler for the underlying rare-event probability is also efficient for the estimation of gradients of awide class of expectations defined on the same rare-event set that satisfies Assumption 1. A generalizationof this observation to settings beyond Assumption 1 (e.g., models with general light-tailed underlying
distributions and general performance functions) is left as future work. Another direction of our future
work is to give a tighter bound for the growth of (3) as γ grows in terms of a polynomial whose degreeand coefficients are determined by model parameters.
In many applications, we want the gradient of expectations of positive functions. An example isexpected loss which is typically positive over the entire region in consideration. In this case, with Theorem
1, we can also show relative efficiency using the expectation before taking derivative as the benchmark.
Corollary 2 (IPA) In Example 1, if f (x, θ ) > 0 for all x ∈ Rd and Assumption 1 holds when the problem
is formulated as (1), then we have


2
∂
2
E ∂ θ f (X̃, θ ) L (X̃)1{g(X̃) ≥ γ}2
E f (X̃, θ )L(X̃)1{g(X̃) ≥ γ}
does not grow exponentially as γ → ∞.
Corollary 3 (LR) In Example 2, suppose that f (x) > 0 for all x ∈ Rd and Assumption 1 holds when the
problem is formulated as (1). Let X̃ be a sample drawn from p(·) as defined in (2) with µ = µ(θ ), thenwe have


2
∂
>
−1
2
E ∂ θ µ(θ ) Σ (X̃ − µ(θ )) f (X̃) L (X̃)1{g(X̃) ≥ γ}2
E f (X̃)L(X̃)1{g(X̃) ≥ γ}
does not grow exponentially as γ → ∞.

Bai, He, Lam, Jiang, and Fu
3

NUMERICAL EXPERIMENTS

In this section, we consider two toy examples and a finance example (pricing a basket option and its delta) todemonstrate the effectiveness of the IS approach in both the performance function and gradient estimation.
3.1 Toy Examples
In the first toy example, we verify the effectiveness of IS for IPA. Suppose that X = (X1 , X2 )> ∼ N(0, I2 )where I2 is the 2 × 2 identity matrix. The performance function isf (X, θ ) = θ 2 X12 + X22 + 0.01,
and the constraint function isg(X) = max {2X1 + X2 , X2 + 2X2 } .

(5)

We would like to estimate E[ f (X, θ )1{g(X) ≥ γ}]. Specifically, we set γ = 7 and θ = 2, and use the IPA
estimator given by 2θ X12 1{g(X) ≥ γ}.
According to (2) and (5), we obtain x1∗ = (1.4, 2.8)> and x2∗ = (2.8, 1.4)> . Therefore, the IS distribution
is a mixture of N(x1∗ , I2 ) with probability 50% and N(x2∗ , I2 ) with probability 50%. For sample size from1000 to 500000, Table 1 shows that the proposed mixture IS works well both for the performance andTable 1: Means, variances, and variance reduction ratios of performance and gradient estimation (IPA) for
MC and IS in the first toy example (true performance and gradient are 4.72E-02 and 3.78E-02, respectively)
Size

IS
mean

1000
2000
5000
10000
20000
50000
1000002000005000004.70E-02
4.72E-02
4.72E-02
4.71E-02
4.72E-02
4.72E-02
4.72E-02
4.72E-02
4.72E-02

PerformanceMCvariancemean
variance

ratio

mean

variancemean

variance

ratio

8.31E-06
3.94E-06
2.01E-06
8.28E-07
4.51E-07
1.44E-07
7.72E-08
5.33E-08
1.61E-08

231.3
231.7
166.4
219.8
168.0
274.7
236.5
180.2
227.8

3.77E-02
3.78E-02
3.77E-02
3.77E-02
3.77E-02
3.78E-02
3.77E-02
3.78E-02
3.78E-02

7.04E-06
3.53E-06
1.74E-06
7.40E-07
3.88E-07
1.30E-07
7.14E-08
4.55E-08
1.38E-08

3.88E-02
4.02E-02
3.89E-02
3.66E-02
3.79E-02
3.79E-02
3.79E-02
3.75E-02
3.78E-02

1.51E-03
7.36E-04
2.84E-04
1.48E-04
6.35E-05
3.43E-05
1.45E-05
8.11E-06
2.89E-06

214.1
208.8
163.4
200.0
163.5
263.2
203.5
178.1
209.2

4.89E-02
5.03E-02
4.82E-02
4.60E-02
4.73E-02
4.71E-02
4.73E-02
4.69E-02
4.72E-02

1.92E-03
9.14E-04
3.35E-04
1.82E-04
7.58E-05
3.97E-05
1.83E-05
9.61E-06
3.66E-06

Gradient
ISMC

gradient estimation (via IPA). For example, when the sample size is 500000, the variance reduction ratios
for performance and gradient estimation are 227.8 and 209.2, respectively. When the sample size is small,e.g., less than 10000, the IS for both performance and gradient estimation has a good accuracy, whileclassical MC does not.
Furthermore, we change the rare-event parameter γ from 4 to 11 with sample size 100000, and calculate
the efficiency ratios of performance estimation and gradient estimation, which are given by


∂
Var
f
(X,
θ
)1{g(X)≥γ}∂θVar ( f (X, θ )1{g(X) ≥ γ})
and
,
(6)
2
(E[ f (X, θ )1{g(X) ≥ γ}])(E[ f (X, θ )1{g(X) ≥ γ}])2respectively, for classical MC, and



∂
Var
f
(
X̃,
θ
)L(
X̃)1{g(X̃)≥γ}
Var f (X̃, θ )L(X̃)1{g(X̃) ≥ γ}∂θ
2 and2 ,E[ f (X̃, θ )L(X̃)1{g(X̃) ≥ γ}]
E[ f (X̃, θ )L(X̃)1{g(X̃) ≥ γ}](7)respectively, for IS. The results are shown in Table 2, which indicates that, as γ increases, the efficiency
ratios of performance and gradient estimation will blow up for MC, while they grow very slowly for IS.

Bai, He, Lam, Jiang, and Fu
Table 2: Efficiency ratios in the first toy example (sample size 105 )γ
Objective
Gradient

IS
MC
IS
MC

4

5

6

7

8

9

10

11

2.31E+00
2.08E+01
2.18E+00
1.79E+01

2.53E+00
5.56E+01
2.37E+00
4.73E+01

3.19E+00
2.07E+02
2.80E+00
1.64E+02

3.88E+00
8.18E+02
3.54E+00
6.46E+02

3.72E+00
3.63E+03
3.08E+00
2.95E+03

4.91E+00
1.97E+04
4.14E+00
1.66E+04

5.10E+00
1.35E+05
4.26E+00
1.08E+05

6.32E+00
9.72E+05
5.17E+00
7.61E+05

In the second toy example, we verify the effectiveness of IS for LR. Suppose that the performance
function isf (X) = X12 + X22 + 0.01,
where X = (X1 , X2 ) ∼ N([θ , 0]> , I2 ). The constraint function is the same as (5). Specifically, we set γ = 7
and θ = 0, and the LR estimator is given by f (X)(X1 − θ )1{g(X) ≥ γ}.
Similar to the first toy example, we vary the sample size from 1000 to 500000 with γ = 7 and vary therare-event parameter γ from 4 to 11 with sample size 100000, respectively, to obtain Tables 3 and 4. Fromthese tables, we observe the similar results as in the first toy example: The mixture IS works well for boththe performance and gradient estimation (via LR). As the rare-event parameter grows, the efficiency ratiosof performance and gradient estimation for IS grow very slowly, while they blow up for classical MC.
Table 3: Means, variances, and variance reduction ratios of performance and gradient estimation (LR) for
MC and IS in the second toy example (true objective and gradient is 1.89E-02 and 4.25E-02, respectively)
Objective
Size

IS

1000
2000
5000
10000
20000
50000
100000200000500000

Gradient
MC

IS

MC

mean

variancemean

variance

ratio

mean

variancemean

variance

ratio

1.89E-02
1.89E-02
1.89E-02
1.89E-02
1.89E-02
1.89E-02
1.89E-02
1.89E-02
1.89E-02

1.10E-06
5.07E-07
2.23E-07
1.17E-07
5.29E-08
2.17E-08
1.29E-08
5.43E-09
1.91E-09

1.76E-02
1.81E-02
1.88E-02
1.84E-02
1.92E-02
1.88E-02
1.90E-02
1.89E-02
1.88E-02

2.28E-04
9.74E-05
5.22E-05
2.49E-05
1.14E-05
4.46E-06
2.32E-06
1.40E-06
4.31E-07

207.9
192.2
234.1
213.0
214.6
206.1
180.6
257.7
226.1

4.25E-02
4.25E-02
4.26E-02
4.25E-02
4.25E-02
4.25E-02
4.25E-02
4.25E-02
4.25E-02

7.13E-06
3.29E-06
1.44E-06
7.19E-07
3.70E-07
1.39E-07
7.88E-08
3.85E-08
1.24E-08

3.97E-02
4.16E-02
4.24E-02
4.13E-02
4.35E-02
4.21E-02
4.25E-02
4.26E-02
4.25E-02

1.33E-03
6.82E-04
3.09E-04
1.72E-04
8.48E-05
2.90E-05
1.60E-05
9.30E-06
3.46E-06

186.5
207.7
214.5
238.5
229.2
209.3
203.0
241.7
278.7

Table 4: Efficiency ratios in the second toy example (sample size 105 )γ
Objective
Gradient

IS
MC
IS
MC

4

5

6

7

8

9

10

11

1.62E+00
1.82E+01
5.12E+00
8.17E+01

1.94E+00
5.13E+01
8.90E+00
2.42E+02

2.51E+00
1.57E+02
1.30E+00
7.50E+02

2.85E+00
7.54E+02
1.86E+01
4.98E+03

3.87E+00
3.17E+03
2.97E+01
2.58E+04

4.68E+00
1.57E+04
4.43E+01
1.69E+05

4.06E+00
1.52E+05
4.16E+01
1.67E+06

6.29E+00
1.29E+06
7.90E+01
1.92E+07

3.2 Basket Option
In this example, we estimate the price of a basket option and its delta, i.e., the partial derivative of theoption price w.r.t. the underlying asset prices. Suppose that there is a basket option with two underlying
assets, which are driven by geometric Brownian motions S1 (t) and S2 (t), i.e.,




√
1 2
1 2
Si (T ) = Si (0) exp (r − σi )T + σiWi (T ) = Si (0) exp (r − σi )T + σi T Xi ,(8)22

Bai, He, Lam, Jiang, and Fu
where Si (T ) is the price of underlying asset i at time T , r is the risk-free interest rate, σi is the volatility of asset i,
and Wi (t), i = 1, 2 are standard Brownian motions, which are mutually independent, i.e., (X1 , X2 )> ∼ N(0, I2 ).
Let X = (X1 , X2 )> and S(T ) = (S1 (T ), S2 (T ))> .
We consider a geometric-average basket option with strike price K, i.e., the payoff function is


 
! 12
! 21
 2

2


V (S(T ), K) =(9)
∏ Si (T ) − K 1  ∏ Si (T ) ≥ K  .
i=1
i=1
1/2
The performance function of interest is E[V (S(T ), K)]. When K is large, { ∏2i=1 Si (T )≥ K} is arare-event set, and can be reformulated as
! 21
! 12

2
2
√
1 2
=
∏ Si (T )∏ Si (0) exp (r − 2 σi )T + σi T Xi
i=1
i=1



2 +σ2
√
1
σ
1
1
2
r−
T (σ1 X1 + σ2 X2 ) ≥ K.
T+
= (S1 (0)S2 (0)) 2 exp
4
2Therefore, setting


1√ >
1
σ12 + σ22
g(X) ,T σ X + (log S1 (0) + log S2 (0)) + r −T,2
2
4
where σ = (σ1 , σ2 )> , the rare event can be written as {g(X) ≥ log K}. Note that g(X) is a linear function.
The objective function in (9) is given in terms of X as
! 12



2
1
σ12 + σ22
1√
>
2
T+
(10)
T (σ X) − K.
∏ Si (T ) − K = (S1 (0)S2 (0)) exp r − 4
2
i=1
We note that the above is exponential in X and is thus heavy-tailed, which does not satisfy Assumption
1. But when K is not too large, we expect that a Taylor series expansion could approximate the original
objective function reasonably well. This approximation is heuristic, but from the numerical examples as
we will show, the approximated objective is quite close to the original quantity. More specifically, let 
β , (S1 (0)S2 (0))1/2 exp r − (σ12 + σ22 )/4 T . Then the payoff function in (9) can be approximated by al-th order Taylor series expansion:
! 12



2
1
σ12 + σ22
1√
2 exp
S
(T
)
−
K
=
(S
(0)S
(0))
r
−
T
+
T
(σ
X
+
σ
X
)
−K
1
2
1 1
2 2
∏ i
4
2
i=1


1√ >
= β exp
Tσ X −K21√ >1 l > l1>2
≈ β −K +β
T σ X + 2 T (σ X) + · · · + l T 2 (σ X)
2
2 2!
2 l!
, f (X).
(11)√
Furthermore, letting A = T σ /2, b = (log S1 (0) + log S2 (0))/2 + r − (σ12 + σ22 )/4 T , and t = log K, thenthe price of the basket option (without discount factor) can be approximated by

 

! 21
! 12
 2

h
i
2
>

E  ∏ Si (T ) − K  1
S
(T
)
≥
K
≈
E
f
(X)1{A
X
+
b
≥
t}
.
i
 ∏

i=1
i=1

Bai, He, Lam, Jiang, and Fu
Note that the RHS above has the form of (1), which satisfies Assumption 1. Therefore, the IS distribution
is chosen as N(x∗ , I2 ), where x∗ is the minimizer of kxk2 over polyhedron {AX + b ≥ t} given by
x∗ =

t −b
A.
kAk2

(12)

Using this IS distribution, we would sample X̃ = (X̃1 , X̃2 )> ∼ N(x∗ , I2 ). Then we can compute S̃i (T ),which is the counterpart of Si (T ) driven by the normal random variable X̃i based on (8). Let L(X̃) be thelikelihood ratio. We have thatE[V (S(T ), K)] ≈ E[ f (X)1{g(X) ≥} log K] = E[ f (X̃)1 g(X̃) ≥ log K L(X̃)].
Gradient Estimation
For gradient estimation, we use the IPA to derive the delta estimator. We consider both the original objectivefunction given in (9) and the approximated function given in (11). For the original objective function, the
IPA estimator of the delta with respect to Si (0) is
! 21

2

1
2Si (0)∏ Si (T )1 {g(X) ≥ log K} ,i=1

and the IPA estimator with IS is1
2Si (0)2! 21∏ S̃i (T )
1 g(X̃) ≥ log K L(X̃).

i=1

For the approximated function (11), the IPA estimator for classical MC is1( f (X) + K) 1 {g(X) ≥ log K} ,2Si (0)and the IPA estimator with IS is 
1f (X̃) + K 1 g(X̃) ≥ log K L(X̃).
2Si (0)
Simulation Results
We investigate the variance reduction effects both for the approximated objective function (11) that satisfies
the polynomial setting in the paper and the original exponential objective function (10) that does notsatisfy the piece-wise polynomial setting. Let S1 (0) = 50, S2 (0) = 60, σ1 = 0.2, σ2 = 0.15, T = 0.1, andK = 60. For the approximated objective function (11), we set l = 4, i.e., we use a fourth-order Taylor
series expansion (a fourth-order polynomial) to approximate the true objective function. We vary the total
simulation budget from 1000 to 500000, and obtain the following tables.
From Tables 5 and 6, we observe the following: (i) The fourth-order polynomial is quite accurate inapproximating the original exponential objective function and also works in the delta estimation. In Table
5, the price means of IS for the fourth-order polynomial approximation are all about 8.90E-03, which is
the true price according to the Black-Scholes formula. And the price means of MC for the fourth-orderpolynomial approximation are also about 8.90E-03 when the sample size is larger than 100000. Similar
results can be found in Table 6. The means of the delta derived by IPA based on the fourth-order polynomial
are close to 6.59E-03, the true delta with respect to S1 (0). (ii) Even though the IS parameter (12) is designed

Bai, He, Lam, Jiang, and Fu
Table 5: Means, variances, and variance reduction ratios of option prices for MC and IS with the fourth-orderpolynomial approximation and the true objective function (true price is 8.90E-03)
Fourth-order polynomial approximation
IS
MC
mean
variancemean
variance

Size
1000
2000
5000
10000
20000
50000
1000002000005000008.93E-03
8.90E-03
8.91E-03
8.90E-03
8.89E-03
8.90E-03
8.90E-03
8.90E-03
8.90E-03

1.15E-07
5.82E-08
2.13E-08
1.26E-08
5.93E-09
2.03E-09
1.22E-09
5.05E-10
2.01E-10

8.54E-03
8.88E-03
8.94E-03
8.93E-03
8.80E-03
8.84E-03
8.92E-03
8.91E-03
8.91E-03

1.26E-05
5.76E-06
3.25E-06
1.08E-06
6.70E-07
2.65E-07
1.36E-07
5.91E-08
2.73E-08

IS
ratio

mean

109.4
98.8
152.5
85.6
113.0
130.6
111.3
117.0
135.7

8.88E-03
8.90E-03
8.89E-03
8.90E-03
8.90E-03
8.90E-03
8.90E-03
8.90E-03
8.90E-03

True objective function
MCvariancemean
variance

ratio

1.07E-07
4.51E-08
2.16E-08
1.02E-08
5.90E-09
1.94E-09
1.26E-09
4.82E-10
2.25E-10

145.7
184.8
158.1
113.7
111.3
150.3
114.8
152.5
119.2

9.20E-03
9.12E-03
8.92E-03
8.87E-03
8.87E-03
8.97E-03
8.92E-03
8.92E-03
8.92E-03

1.55E-05
8.33E-06
3.42E-06
1.16E-06
6.57E-07
2.92E-07
1.45E-07
7.35E-08
2.68E-08

Table 6: Means, variances, and variance reduction ratios of deltas with respect to S1 (0) for MC and IS
with the fourth-order polynomial approximation and the true objective function (true delta with respect toS1 (0) is 6.59E-03)
Delta for polynomial approximation
IS
MC
mean
variancemean
variance

Size
1000
2000
5000
10000
20000
50000
1000002000005000006.55E-03
6.59E-03
6.58E-03
6.59E-03
6.58E-03
6.59E-03
6.59E-03
6.60E-03
6.59E-03

1.16E-07
5.88E-08
2.36E-08
1.28E-08
5.91E-09
2.08E-09
1.24E-09
4.49E-10
2.30E-10

6.41E-03
6.54E-03
6.68E-03
6.59E-03
6.56E-03
6.58E-03
6.61E-03
6.60E-03
6.60E-03

3.46E-06
1.72E-06
9.35E-07
3.68E-07
1.72E-07
7.40E-08
4.07E-08
2.11E-08
8.38E-09

ratio
29.8
29.3
39.6
28.8
29.1
35.5
32.7
47.1
36.4

Delta for true objective function
IS
MC
mean
variancemean
variance
6.61E-03
6.57E-03
6.59E-03
6.59E-03
6.59E-03
6.59E-03
6.59E-03
6.59E-03
6.59E-03

1.08E-07
5.63E-08
2.24E-08
1.27E-08
6.70E-09
2.03E-09
1.13E-09
6.04E-10
2.39E-10

6.81E-03
6.64E-03
6.58E-03
6.56E-03
6.57E-03
6.63E-03
6.59E-03
6.61E-03
6.59E-03

3.66E-06
2.21E-06
9.52E-07
3.62E-07
1.84E-07
7.64E-08
4.04E-08
1.96E-08
7.09E-09

ratio
33.8
39.3
42.6
28.6
27.5
37.7
35.8
32.4
29.6

for the polynomial function, it also works for the original exponential objective function. For example,
in Table 5, when the sample size is 200000, the variance reduction ratio for the fourth-order polynomial
approximation is 117.0, while the variance reduction ratio for the true objective function is 152.5. (iii)
The proposed IS parameter (12) works well for both the price estimation (objective estimation) and delta
estimation (gradient estimation). As seen in the tables, when the sample size is small (1000 to 20000), theIS method has good accuracy but the MC does not. The variance reduction ratios for the price estimationand delta estimation are about 100 and 30, respectively.
Next, we change K from 58 to 70, so that {g(X) ≥ log K} becomes rarer, and we use the efficiency ratios
defined in (6) and (7) to measure the efficiency of the price and delta estimation. Let the sample size be 107 .
Then we obtain Tables 7 and 8. From the tables, we can see that for both the approximated fourth-order
Table 7: Efficiency ratios for IS and MC in price estimation (sample size 107 )KISMCApproximatedTrue objective
ApproximatedTrue objective

58

60

62

64

66

68

70

1.14E+00
1.14E+00
2.60E+01
2.60E+01

1.27E+00
1.27E+00
1.73E+02
1.73E+02

1.71E+00
1.84E+00
2.11E+03
2.14E+03

2.20E+00
1.94E+00
4.62E+04
4.63E+04

2.28E+00
2.28E+00
1.62E+06
1.62E+06

2.99E+00
3.04E+00
8.42E+07
8.71E+07

3.63E+00
3.62E+00
6.38E+10
6.37E+10Bai, He, Lam, Jiang, and Fu
Table 8: Efficiency ratios for IS and MC in delta estimation (sample size 107 )K
IS
MCapproximatedTrue objective
approximatedTrue objective

58

60

62

64

66

68

70

5.49E-01
5.49E-01
4.74E+00
4.74E+00

1.18E+00
1.18E+00
4.94E+01
4.94E+01

2.74E+00
2.82E+00
9.66E+02
1.16E+03

5.04E+00
4.67E+00
2.90E+04
2.68E+04

7.82E+00
7.81E+00
1.27E+06
1.27E+06

1.27E+01
1.36E+01
8.51E+07
8.11E+07

1.83E+01
1.83E+01
1.64E+10
1.63E+10

polynomial function and the original exponential objective function, as the rare-event parameter becomeslarger, the efficiency ratios when using IS do not grow significantly, whereas the efficiency ratios growvery fast when using classical MC.
4CONCLUSION

This paper studies importance sampling for stochastic gradient estimation involving rare events. We develop
theoretical results to show that, for a class of problems, an efficient IS mixture sampler for estimating
the probability of the underlying rare event will also be efficient for gradient estimations of a wide classof performance measures defined on the same rare-event set. Several simulation examples, including afinancial derivatives example, are provided and support the theory.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1834710, IIS-1849280 and IIS-2123684, from the Air Force Office of Scientific Research under grantFA95502010211, and from the National Natural Science Foundation of China under grants 72171060 and
71801148.
REFERENCES
Bai, Y., Z. Huang, H. Lam, and D. Zhao. 2022. “Rare-Event Simulation for Neural Network and Random Forest Predictors”.
ACM Transactions on Modeling and Computer Simulation Forthcoming.
Blanchet, J. H., and H. Lam. 2012. “State-dependent importance sampling for rare-event simulation: An overview and recent
advances”. Surveys in Operations Research and Management Science 17(1):38–59.
Bucklew, J. A. 2004. Introduction to Rare Event Simulation. New York, NY: Springer Science & Business Media.
Driggs, D., M. J. Ehrhardt, and C.-B. Schönlieb. 2022. “Accelerating variance-reduced stochastic gradient methods”. Mathematical
Programming 191(2):671–715.
Fox, B. L., and P. W. Glynn. 1989. “Replication Schemes For Limiting Expectations”. Probability in the Engineering andInformational Sciences 3(3):299–318.
Fu, M. C. 2006. “Gradient Estimation”. Handbooks in Operations Research and Management Science 13:575–616.
Fu, M. C., and J. Q. Hu. 1997. Conditional Monte Carlo: Gradient Estimation and Optimization Applications. Norwell, MA:
Kluwer Academic Publishers.
Glasserman, P. 1988. “Performance continuity and differentiability in Monte Carlo optimization”. In 1988 Winter Simulation
Conference Proceedings, 518–524. IEEE.
Glynn, P. 1989. “Optimization of Stochastic Systems via Simulation”. In 1989 Winter Simulation Conference Proceedings,edited by E. A. MacNair, K. J. Musselman, and P. Heidelberger, 90–105. Piscataway, New Jersey: Institute of Electricaland Electronics Engineers, Inc.
Glynn, P. W. 1990. “Likelihood ratio gradient estimation for stochastic systems”. Communications of the ACM 33(10):75–84.
He, S., G. Jiang, H. Lam, and M. C. Fu. 2022. “Adaptive importance sampling for efficient stochastic root finding and quantile
estimation”. working paper, https://arxiv.org/abs/2102.10631.
Heidergott, B., F. Vazquez-Abad, G. Pflug, and T. Farenhorst-Yuan. 2010. “Gradient estimation for discrete-event systems bymeasure-valued differentiation”. ACM Transactions on Modeling and Computer Simulation 20(1):No. 5.
Ho, Y. C., and X. Cao. 1983. “Perturbation Analysis and Optimization of Queueing Networks”. Journal of Optimization Theory
and Applications 40:559–582.
Johnson, R., and T. Zhang. 2013. “Accelerating stochastic gradient descent using predictive variance reduction”. In Advancesin Neural Information Processing Systems, 315–323. Red Hook, NY, USA: Curran Associates Inc.

Bai, He, Lam, Jiang, and Fu
Juneja, S., and P. Shahabuddin. 2006. “Rare-Event Simulation Techniques: An Introduction and Recent Advances”. In Handbooksin Operations Research and Management Science, edited by S. G. Henderson and B. L. Nelson, Volume 13, Chapter 11,
291 – 350. Amsterdam, The Netherlands: Elsevier.
L’Ecuyer, P. 1990. “A unified view of the IPA, SF, and LR gradient estimation techniques”. Management Science 36(11):1364–1383.
Liu, G., and L. J. Hong. 2011. “Kernel Estimation of the Greeks for Options with Discontinuous Payoffs”. Operations
Research 59:96–108.
Nakayama, M. K. 1995. “Asymptotics of Likelihood Ratio Derivative Estimators in Simulations of Highly Reliable Markovian
Systems”. Management Science 41(3):524–554.
Peng, Y. J., M. C. Fu, J. Q. Hu, and B. Heidergott. 2018. “A new unbiased stochastic derivative estimator for discontinuous
sample performances with structural parameters”. Operations Research 66(2):487–499.
Pflug, G. C. 1989. “Sampling derivatives of probabilities”. Computing 42:315–328.
Reiman, M. I., and A. Weiss. 1989. “Sensitivity analysis for simulations via likelihood ratios”. Operations Research 37(5):830–844.
Rubinstein, R. Y. 1989. “Sensitivity analysis of computer simulation models via the score efficient”. Operations Research 37:72–81.
Wang, C., X. Chen, A. J. Smola, and E. P. Xing. 2013. “Variance Reduction for Stochastic Gradient Optimization”. In Advancesin Neural Information Processing Systems, edited by C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger,
Volume 26: Curran Associates, Inc.
Wang, Y., M. C. Fu, and S. I. Marcus. 2012. “A new stochastic derivative estimator for discontinuous payoff functions withapplication to financial derivatives”. Operations Research 60(2):447–460.
Xiao, L., and T. Zhang. 2014. “A Proximal Stochastic Gradient Method with Progressive Variance Reduction”. SIAM Journal
on Optimization 24(4):2057–2075.
Zhao, P., and T. Zhang. 2015, 07–09 Jul. “Stochastic Optimization with Importance Sampling for Regularized Loss Minimization”.
In Proceedings of the 32nd International Conference on Machine Learning, edited by F. Bach and D. Blei, Volume 37,
1–9. Lille, France: PMLR.

AUTHOR BIOGRAPHIESYUANLU BAI is a PhD candidate in the Department of Industrial Engineering and Operations Research at Columbia University.
She received an M.S. in operations research from Columbia University, and a B.S. in statistics as well as a B. Ec in economics
from Peking University. Her research interest lies in uncertainty quantification and rare-event simulation. Her email address isyb2436@columbia.edu.
SHENGYI HE is a PhD student in the Department of Industrial Engineering and Operations Research at Columbia University.
He received his B.S. degree in statistics from Peking University in 2019. His research interests include variance reduction anduncertainty quantification via stochastic and robust optimization. His email address is sh3972@columbia.edu.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. He received his Ph.D. degree in statistics from Harvard University in 2011, and was on the faculty of Boston
University and the University of Michigan before joining Columbia in 2017. His research focuses on Monte Carlo simulation,
uncertainty quantification, risk analysis, and stochastic and robust optimization. His email address is henry.lam@columbia.edu.
GUANGXIN JIANG is a professor in the School of Management at Harbin Institute of Technology. He earned his Ph.D.
in Applied Mathematics from Tongji University. His research interests include stochastic models and simulation, machinelearning, financial engineering and risk management, FinTech, etc. His email address is gxjiang@hit.edu.cn.
MICHAEL C. FU holds the Smith Chair of Management Science in the Robert H. Smith School of Business, with a joint
appointment in the Institute for Systems Research and affiliate faculty appointment in the Department of Electrical and Computer
Engineering, all at the University of Maryland, College Park. His research interests include stochastic gradient estimation,
simulation optimization, and applied probability. He served as WSC2011 Program Chair and received the INFORMS Simulation
Society’s Distinguished Service Award in 2018. He is a Fellow of INFORMS and IEEE. His e-mail addresses is mfu@umd.edu.