Yijie Peng

Michael C Fu

Jiaqiao Hu

Pierre L'ecuyer

Bruno Tuffin

Michael C Fu The

Robert H Smith

Variance Reduction for Generalized Likelihood Ratio Method in Quantile Sensitivity Estimation

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Quantile, also known as value-at-risk (VaR), is an important risk measure. Quantile sensitivity estimation has been actively studied in simulation due to its centrality in risk management. A seminal work is [START_REF] Hong | Estimating Quantile Sensitivities[END_REF], which proposes a conditional infinitesimal perturbation analysis (IPA). [START_REF] Jiang | Technical Note -On Estimating Quantile Sensitivities Via Infinitesimal Perturbation Analysis[END_REF] avoid the need of computing the conditional expectation for some special cases. [START_REF] Liu | Kernel Estimation of Quantile Sensitivities[END_REF] propose a kernel-based estimator for estimating the quantile sensitivity, which can be extended to estimate sensitivities of conditional VaR (CVaR) [START_REF] Hong | Simulating Sensitivities of Conditional Value At Risk[END_REF]. A conditional Monte Carlo (CMC) method is applied to estimate quantile sensitivity in [START_REF] Fu | Conditional Monte Carlo Estimation of Quantile Sensitivities[END_REF]. [START_REF] Heidergott | A Measure-valued Differentiation Approach to Sensitivity Analysis of Quantiles[END_REF] propose a measure-valued differentiation (MVD) estimator for quantile sensitivity. CMC and MVD for quantile sensitivity rely on estimating two distribution sensitivities, i.e., the derivatives of the distribution function with respect to both argument and parameters in the stochastic model.

The difficulty in distribution sensitivity estimation lies in the discontinuity introduced by the indicator function. Recently, [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF] propose a generalized likelihood ratio (GLR) method to handle discontinuities for a wide scope of sample performance in sensitivity analysis. [START_REF] Peng | Maximum Likelihood Estimation By Monte Carlo Simulation: Towards Data-Driven Stochastic Modeling[END_REF] develop GLR estimators for any distribution sensitivity and use it to calculate maximum likelihood estimation for complex stochastic models. [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF] derive a GLR estimator for stochastic models with uniform random numbers as input, relaxing a requirement in [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF] that the tails of the input distribution of the stochastic model go smoothly to zero fast enough. [START_REF] Peng | On the asymptotic analysis of quantile sensitivity estimation by Monte Carlo simulation[END_REF] apply GLR to estimate quantile sensitivity, and [START_REF] Peng | Estimating quantile sensitivity for financial models with correlations and jumps[END_REF] extend the GLR quantile sensitivity estimator to stochastic model with correlations and jumps. [START_REF] Glynn | Computing Sensitivities for Distortion Risk Measures[END_REF] use the GLR estimators for distribution sensitivities to estimate derivatives of distortion risk measures, which cover VaR and CVaR as special cases.

CMC methods can reduce the variance and smooth the performance function in simulation by conditioning on certain events or random variables and then integrating out the remaining randomness [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. [START_REF] Fu | Conditional Monte Carlo Estimation of Quantile Sensitivities[END_REF] uses CMC to smooth the discontinuity introduced by the indicator function, after which IPA is applied to differentiate the conditional expectation. GLR does not need smoothing to obtain an unbiased distribution sensitivity estimator, but CMC can be applied afterward to reduce the variance for GLR.

As other variance reduction techniques, RQMC methods replace the vectors of uniform random numbers that drive independent simulation runs by dependent vectors of uniform random numbers that cover the space more evenly. When estimating an expectation, they can provide an unbiased estimator with its variance converging to zero at a rate that is faster than with Monte Carlo [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Method[END_REF][START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. Such a faster rate can be proved when the estimator inside the expectation is sufficiently smooth as a function of the underlying uniform random numbers. When the estimator is not smooth (e.g., discontinuous), the convergence rate may not be improved, but RQMC could still reduce the variance by a constant factor. We show, through simulation experiments, that the variance of the quantile sensitivity estimator by GLR can be significantly reduced by appropriately combining the method with CMC and RQMC. Similar use of CMC and RQMC for reducing the variance of quantile estimation can be found in [START_REF] Nakayama | Quantile Estimation Via a Combination of Conditional Monte Carlo and Randomized Quasi-Monte Carlo[END_REF].

The rest of the paper is organized as follows. Section 2 introduces the quantile sensitivity estimation problem and the GLR method. Variance reduction by CMC and RQMC is discussed in Section 3. Section 4 provide applications and simulation experiments. The last section offers conclusions.

QUANTILE SENSITIVITY ESTIMATION

For 0 ≤ α ≤ 1, the α-quantile (also known as value-at-risk) of an output random variable Y θ = h(X; θ ) with input random variables X = (X 1 , . . . , X n ) and cumulative distribution function (cdf) F(•; θ ) is defined as q α (θ ) := arg inf{y : F(y; θ ) ≥ α}.

(1)

When F(•; θ ) is continuous, q α (θ ) = F -1 (α; θ ). Our goal is to estimate the derivative of the α-quantile with respect to parameter θ , dq α (θ )/dθ . Assume h(X; θ ) is a continuous r.v. with a positive and continuous density f (y; θ ) on (q α (θ )ε, q α (θ ) + ε), ε > 0. Using the formula for the derivative of an implicit function, we have [START_REF] Fu | Conditional Monte Carlo Estimation of Quantile Sensitivities[END_REF])

d dθ q α (θ ) = - ∂ F(y; θ ) ∂ θ y=q α (θ ) f (q α (θ ); θ ). Notice that ∂ F(y; θ ) ∂ θ = ∂ E [1{h(X; θ ) ≤ y}] ∂ θ , f (y; θ ) = ∂ E [1{h(X; θ ) ≤ y}] ∂ y .
From [START_REF] Peng | Maximum Likelihood Estimation By Monte Carlo Simulation: Towards Data-Driven Stochastic Modeling[END_REF], we have the GLR estimators for these distribution sensitivities: G 1,i (X; y, θ ) := 1{h(X; θ ) ≤ y}ψ 1,i (X; θ ) and G 2,i (X; y, θ ) := 1{h(X; θ ) ≤ y}ψ 2,i (X; θ ), i = 1, . . . , n, such that

∂ F(y; θ ) ∂ θ = E [G 1,i (X; y, θ )] , f (y; θ ) = E [G 2,i (X; y, θ )] ,
where

ψ 1,i (x; θ ) = ∂ log f X (x; θ ) ∂ θ - ∂ h(x; θ ) ∂ x i -1 ∂ h(x; θ ) ∂ θ ∂ log f X (x; θ ) ∂ x i - ∂ 2 h(x; θ ) ∂ x 2 i ∂ h(x; θ ) ∂ x i -1 + ∂ 2 h(x; θ ) ∂ θ ∂ x i , ψ 2,i (x; θ ) = ∂ h(x; θ ) ∂ x i -1 ∂ log f X (x; θ ) ∂ x i - ∂ 2 h(x; θ ) ∂ x 2 i ∂ h(x; θ ) ∂ x i -1
, with f X (•; θ ) being the (joint) density of the vector of input random variables X. From [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF], we have the GLR estimators for distribution sensitivities when X is a vector of uniform random numbers U = (U 1 , . . . ,U n ):

G 1,i (U; y, θ ) := 1{h(U i ; θ ) -y ≤ 0}r i (U i ; θ ) -1{h(U i ; θ ) -y ≤ 0}r i (U i ; θ ) + 1{h(U; θ ) -y ≤ 0}d(U; θ ),
where

r i (u; θ ) = ∂ h(u; θ ) ∂ u i -1 ∂ h(u; θ ) ∂ θ , d(u; θ ) = ∂ h(u; θ ) ∂ u i -1 ∂ h(u; θ ) ∂ u i -1 ∂ h(u; θ ) ∂ θ ∂ 2 h(u; θ ) ∂ u 2 i - ∂ 2 h(u; θ ) ∂ u i ∂ θ , and
G 2,i (U; z, θ ) := 1{h(U i ; θ ) -y ≤ 0}r i (U i ; θ ) -1{h(U i ; θ ) -y ≤ 0}r i (U i ; θ ) + 1{h(U; θ ) -y ≤ 0} d(U; θ ),
where

U i := (U 1 , . . . , 1 - ith element , . . . ,U n ), U i := (U 1 , . . . , 0 + ith element , . . . ,U n ), and ri (u; θ ) = - ∂ h(u; θ ) ∂ u i -1 d(u; θ ) = - ∂ h(u; θ ) ∂ u i -2 ∂ 2 h(u; θ ) ∂ u 2 i .
Let X ( j) , j = 1, . . . , m, be i.i.d. realizations of X, and Fm (•) be the empirical distribution of h(X ( j) ; θ ), j = 1, . . . , m. The empirical α-quantile F-1 m (α) is the inverse of the empirical cdf evaluated at α, defined as in (1) with F replaced by F-1 m . Then the quantile sensitivity can be estimated by

- ∑ m j=1 1{h(X ( j) ; θ ) ≤ y}ψ 1,i (X ( j) ; θ ) ∑ m j=1 1{h(X ( j) ; θ ) ≤ y}ψ 2,i (X ( j) ; θ ) y= F-1 m (α)
.

(2) This is a random function of a ratio of expectations by a ratio of averages evaluated at a random estimate.

Even if each term in the ratio ( 2) is an unbiased estimator of the corresponding expectation, the ratio (2) is not an unbiased estimator of the ratio of expectations. There is also another source of bias for the sensitivity: the quantile estimator itself is also biased. Also, the two terms in (2) and quantile estimators are not independent because they are estimated by a same batch of inputs and outputs: X ( j) and h(X ( j) ; θ ), j = 1, . . . , m. The strong consistency and a central limit theorem for (2) have been established in [START_REF] Peng | On the asymptotic analysis of quantile sensitivity estimation by Monte Carlo simulation[END_REF] using an empirical process theory. [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF] consider sensitivity analysis for more general stochastic models:

∂ ∂ θ E[ϕ(g(X; θ ))], (3)
where ϕ : R n → R is a measurable function not necessarily continuous, g(x; θ ) = (g 1 (x; θ ), . . . , g n (x; θ )) T is a vector of functions with sufficient smoothness for x ∈ R n , and X = (X 1 , . . . , X n ) is a vector of input random variables with a joint density f X (x; θ ) supported on Ω ⊆ R n . The Jacobian of g(•; θ ) is

J g (x; θ ) :=       ∂ g 1 (x;θ ) ∂ x 1 ∂ g 1 (x;θ ) ∂ x 2 • • • ∂ g 1 (x;θ ) ∂ x n ∂ g 2 (x;θ ) ∂ x 1 ∂ g 2 (x;θ ) ∂ x 2 • • • ∂ g 2 (x;θ ) ∂ x n . . . . . . . . . . . . ∂ g n (x;θ ) ∂ x 1 ∂ g n (x;θ ) ∂ x 2 • • • ∂ g n (x;θ ) ∂ x n      
.

The Jacobian matrix is required to be invertible almost everywhere, which justifies local invertibility of function g by the implicit function theory. The conditions for justifying unbiasedness are referred to [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF] and [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF]. Let e i be the i-th unit vector. We define

d(x; θ ) := n ∑ i=1 e T i J -1 g (x; θ ) (∂ x i J g (x; θ )) J -1 g (x; θ )∂ θ g(x; θ ) -trace J -1 g (x; θ )∂ θ J g (x; θ ) -J -1 g (x; θ )∂ θ g(x; θ ) T ∇ x log f X (x; θ ),
where and ∂ y M(y) is the matrix obtained by differentiating M with respect to y element-wise. From [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF], when the tails of f X (•; θ ) go smoothly to zero fast enough, then, under certain regularity conditions, it can be shown that an unbiased GLR estimator for ( 3) is given by

G(X; θ ) := ϕ(g(X; θ ))w(X; θ ) with w(x; θ ) := ∂ log f X (x; θ ) ∂ θ + d(x; θ ). (4)
As a particular stochastic activity network (SAN) example, the output could be the maximum of the durations of activities on different paths. By the GLR method under general framework (4), we can consider quantile sensitivity for max i=1,...,n g i (X; θ ), which has a distribution function expressed as [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF] consider the case when X is a vector of uniform random numbers U = (U 1 , . . . ,U n ) such that Ω = (0, 1) n and ∂ log f X (x; θ )/∂ θ = 0. Then under certain regularity conditions, we have the following unbiased GLR estimator for (3):

E 1 max i=1,...,n g i (X; θ ) ≤ y = E n ∏ i=1 1 {g i (X; θ ) -y ≤ 0} , with ϕ(y) = ∏ n i=1 1{y i ≤ 0}.
G(U; θ ) := n ∑ i=1 ϕ(g(U i ; θ ))r i (U i ; θ ) -ϕ(g(U i ; θ ))r i (U i ; θ ) + ϕ(g(U; θ ))d(U; θ ), (5)
where

r i (x; θ ) := J -1 g (x; θ ) ∂ θ g(x; θ ) T e i , i = 1, . . . , n.

VARIANCE REDUCTION

CMC methods can reduce the variance and smooth the performance function in simulation by conditioning on certain events or random variables and then integrating out the remaining randomness [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]). For an estimator H(Z), we have

E[H(Z)] = E[ H( Z)],
where

H( Z) := E[H(Z)| Z],
with Z being a subset of input random variables in Z. The variance reduction for the conditional estimator H( Z) can be seen from the following variance decomposition formula:

Var(H(Z)) = Var(E[H(Z)| Z]) + E[Var(H(Z)| Z)] ≥ Var( H( Z)).
Typically, H( Z) is smoother than H(Z), due to the integration taken in the conditional expectation. Quasi-Monte Carlo (QMC) refers to a class of deterministic numerical integration methods in which the integrand is evaluated at a fixed set of m points, and the average is used as an approximation. One limitation of the method is that it is very hard to estimate the approximation error in practice. RQMC takes the QMC points and randomizes them in a way that each point has the uniform distribution over (0, 1) n , so that each randomized point represents a proper realization of U while the set of m points still covers the unit hypercube (0, 1) n more uniformly than typical independent random points (so the points are not independent) (L'Ecuyer 2018).

In general, for a given function p, RQMC estimates the integral µ = (0,1) n p(u)du by the average

μm := 1 m m ∑ j=1 p(U ( j) ),
where U (1) , . . . ,U (m) form an RQMC point set. The most common types of QMC point set constructions are lattice rules, polynomial lattices rules, and digital nets [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] L'ecuyer | Quasi-Monte Carlo methods with applications in finance[END_REF].

For lattice rules, an appropriate randomization is a random shift modulo 1, which adds a single uniform random point to all the lattice points, and retains the shifted points that are in (0, 1) n as the m RQMC points. This randomization preserves the lattice structure, and there are explicit expressions for Var[ μm ] in terms of the Fourier coefficients of h, and computable bounds on this variance for certain classes of smooth functions (L'Ecuyer and Lemieux 2000; L'Ecuyer and Munger 2012; L'Ecuyer and Munger 2016). When the mixed derivatives of p are sufficiently smooth, the variance can converge at a faster rate than O(m -1 ), sometimes nearly O(m -2 ) and even faster in some cases. When p is not smooth (e.g., discontinuous), these convergence rate results do not apply, although weaker results do apply [START_REF] He | On the Convergence Rate of Randomized Quasi-Monte Carlo for Discontinuous Functions[END_REF], and even when the convergence rate is not improved, the variance is often reduced by a constant factor. For polynomial lattices rules and digital nets in general, which include Sobol' points, the random shift does not preserve the structure and net properties, but other appropriate randomizations do, including nested uniform scrambling, some affine scrambles, and random digital shifts. Variance bounds and convergence rate results are also available for these rules [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF][START_REF] L'ecuyer | Quasi-Monte Carlo methods with applications in finance[END_REF]. Our model formulation (5) in terms of a function of independent U (0, 1) random variables makes it an obvious candidate for the application of RQMC, which is designed exactly for this type of formulation. For formulation (4), suppose X can be generated by Γ(U), and RQMC can be applied. For example, when X 1 , . . . , X n are independent random variables with marginal distribution functions F 1 , . . . , F n , then they can be generated by

X i = F -1 i (U i ), i = 1, . . . , n.
For the α-quantile sensitivity estimation, the cdf F(•; θ ) of p(U) can be estimated by its empirical RQMC counterpart

F m (y) := 1 m m ∑ j=1 1{p(U ( j) ) ≤ y},
where {U (1) , . . . ,U (m) } is an RQMC point set, and the quantile q α (θ ) can be estimated by the pseudo-inverse F -1 m (α). We can also use RQMC points in the sensitivity estimate (2).

NUMERICAL EXPERIMENTS

We apply the proposed method to two applications. GLR is compared with the finite difference method using common random numbers (FDC(δ )), where δ is the perturbation size. The GLR method together with CMC is called conditional GLR (CGLR), and CGLR together with RQMC is denoted as CGLRQ.

For RQMC, we use the Sobol sequence scrambled by the algorithm of [START_REF] Matousek | On the L 2 -discrepancy for anchored boxes[END_REF] in Matlab. In the two examples, the sample size for estimating quantile sensitivities are set as m = 2 15 for the standard Monte Carlo and RQMC estimators, and the quantile senstivities and standard errors of the estimators are estimated by 10 4 independent experiments.

TOY EXAMPLE

To illustrate, we estimate quantile sensitivity of a simple stochastic model θ X 1 + X 2 +U, where X 1 , X 2 are standard normal random variables with the common cdf Φ(•), and U is an independent uniform random variable U (0, 1). In the stochastic model for deriving GLR estimators, the input random variable could either be taken as the normal random variables X 1 , X 2 or uniform random number U. If we choose U to be the input random variable, then h(u; θ ,

x 1 , x 2 ) = θ x 1 + x 2 + u, ∂ h(u; θ , x 1 , x 2 )/∂ u = 1, ∂ h(u; θ , x 1 , x 2 )/∂ θ = x 1 ,
and other derivatives in (5) are zeros, so GLR estimators for distribution sensitivities are given by

G 1,1 (U; θ , X 1 , X 2 ) = 1{θ X 1 + X 2 + 1 ≤ y}X 1 -1{θ X 1 + X 2 ≤ y}X 1 , G 2,1 (U; θ , X 1 , X 2 ) = 1{θ X 1 + X 2 ≤ y} -1{θ X 1 + X 2 + 1 ≤ y}.
By conditioning on X 1 , we have

E[ G 1,1 (U; y, θ , X 1 , X 2 )] = E E 1{θ X 1 + X 2 + 1 ≤ y}X 1 -1{θ X 1 + X 2 ≤ y}X 1 X 1 = E[(Φ(y -θ X 1 -1) -Φ(y -θ X 1 ))X 1 ], E[ G 2,1 (U; y, θ , X 1 , X 2 )] = E E 1{θ X 1 + X 2 ≤ y} -1{θ X 1 + X 2 + 1 ≤ y} X 1 = E[Φ(y -θ X 1 ) -Φ(y -θ X 1 -1)].
We can take the expressions inside the expectations on the right-hand sides of the above equations as conditional GLR estimators. From this derivation, three quantile sensitivity estimators of the form (2) using GLR, CGLR, and CGLRQ to estimate distribution sensitivities are referred to as GLR-1, CGLR-1, and CGLRQ-1, respectively. If we choose X 1 to be the input random variable, then h(x 1 , x

2 ; θ , u) = θ x 1 +x 2 +u, ∂ h(x 1 , x 2 ; θ , u)/∂ x 1 = θ , ∂ h(x 1 , x 2 ; θ , u)/∂ θ = x 1 , ∂ 2 h(x 1 , x 2 ; θ , u)/∂ θ ∂ x 1 = 1, ∇ x log f X (x 1 , x 2 ) = -(x 1 , x 2 )
, and the other derivatives in (4) are zeros, so GLR estimators for distribution sensitivities can be given by

G 1,1 (X 1 , X 2 ; y, θ ,U) = 1 θ 1{θ X 1 + X 2 +U ≤ y}(X 2 1 -1), G 2,1 (X 1 , X 2 ; y, θ ,U) = - 1 θ 1{θ X 1 + X 2 +U ≤ y}X 1
By conditioning on X 1 ,U, we have

E[G 1,1 (X 1 , X 2 ; y, θ ,U)] = E E 1 θ 1{θ X 1 + X 2 +U ≤ y}X 2 1 X 1 ,U = 1 θ E[Φ(y -θ X 1 -U)X 2 1 ], E[G 2,1 (X 1 , X 2 ; y, θ ,U)] = E E - 1 θ 1{θ X 1 + X 2 +U ≤ y}X 1 X 1 ,U = - 1 θ E[Φ(y -θ X 1 -U)X 1 ].
We can take the expressions inside the expectations on the right-hand sides of the above equations as conditional GLR estimators. From this derivation, three quantile sensitivity estimators of the form (2) using Table 1:

α = 0.1 α = 0.3 α = 0.5 α = 0.7 α = 0.9 GLR-1
-0.888 ± 1.4 × 10 -4 -0.363 ± 1.0 × 10 -4 0 ± 1.0 × 10 -4 0.363 ± 1.0 × 10 -4 0.887 ± 1.4 × 10 -4 CGLR-1 -0.888 ± 1.0 × 10 -4 -0.363 ± 8 × 10 -5 0 ± 7 × 10 -5 0.364 ± 8 × 10 -5 0.888 ± 1.0 × 10 -4 CGLRQ-1 -0.888 ± 4.7 × 10 -5 -0.363 ± 4.1 × 10 -5 0 ± 3.9 × 10 -5 0.363 ± 4.1 × 10 -5 0.888 ± 4.7 × 10 -5 GLR-2 -0.888 ± 2.3 × 10 -4 -0.363 ± 1.8 × 10 -4 0 ± 2.1 × 10 -4 0.363 ± 2.9 × 10 -4 0.890 ± 7.3 × 10 -4 CGLR-2 -0.888 ± 2.0 × 10 -4 -0.363 ± 1.8 × 10 -4 0 ± 2.0 × 10 -4 0.364 ± 2.8 × 10 -4 0.889 ± 7.2 × 10 -4 CGLRQ-2 -0.888 ± 4.9 × 10 -5 -0.363 ± 4.1 × 10 -5 0 ± 3.9 × 10 -5 0.363 ± 4.1 × 10 -5 0.888 ± 4.9 × 10 -5 FDC(0.01) -0.888 ± 1.2 × 10 -3 -0.364 ± 8.5 × 10 -4 0 ± 7.9 × 10 -4 0.365 ± 8.6 × 10 -4 0.890 ± 1.2 × 10 -3 FDC(0.1) -0.911 ± 3.9 × 10 -4 -0.373 ± 2.7 × 10 -4 0 ± 2.6 × 10 -4 0.372 ± 2.8 × 10 -4 0.910 ± 3.9 × 10 -4

Figure 1: A SAN with seven activities.

GLR, CGLR, and CGLRQ to estimate distribution sensitivities are referred to as GLR-2, CGLR-2, and CGLRQ-2, respectively.

For our experiments, we took θ = 1 and five probability levels: α = 0.1, α = 0.3, α = 0.5, α = 0.7, and α = 0.9. From Table 1, we can see that the quantile sensitivity estimates by GLR-1, CGLR-1, CGLRQ-1, GLR-2, CGLR-2, CGLRQ-2, and FDC(0.01) are statistically indistinguishable, whereas the quantile sensitivity estimates by FDC(0.1) lie outside of the 90% confidence intervals of other methods. GLR leads to smaller variance than FDC, the variances of GLR-1, CGLR-1, CGLRQ-1 are smaller than GLR-2, CGLR-2, CGLRQ-2, respectively, and the CMC and RQMC lead to significant variance reduction.

STOCHASTIC ACTIVITY NETWORK (SAN)

We estimate quantile sensitivity for a simple SAN. There are three different paths representing the tasks to reach the final stage of a project, i.e., (1, 4, 6), (2, 5, 6), (1,3,5,6). The completion time for the entire project is max(Y 1 +Y 4 +Y 6 ,Y 2 +Y 5 +Y 6 ,Y 1 +Y 3 +Y 5 +Y 6 ), and the sample performance for the cdf of the completion time is

1 {max(Y 1 +Y 4 +Y 6 ,Y 2 +Y 5 +Y 6 ,Y 1 +Y 3 +Y 5 +Y 6 ) ≤ y} .
We assume that the first three activities follow independent exponential distributions:

Y i = -1 λ i log(U i ), i = 1, 2,
3, and the other three activities follow independent log-normal distributions, Y i = exp (µ i + σ i X i ), i = 4, 5, 6, where X i 's are standard normal random variables. Let λ 1 = θ . The distribution function F(y; θ ) of the completion time Y is

F(y; θ ) = E[1 {Y 1 + max(Y 4 ,Y 3 +Y 5 ) +Y 6 -y ≤ 0} 1 {Y 2 +Y 5 +Y 6 -y ≤ 0}]
If we keep the original indicator above, then the image of g is in R 3 and we need three inputs. With the rewriting, the image of g is in R 2 and we select two inputs in a way that the Jacobian is invertible. To estimate the two distribution sensitivities, we can view (U 1 ,U 2 ) as the input random variables in the stochastic model of [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF], and we have

g 1 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 ) = - 1 λ 1 logU 1 + max(Y 4 ,Y 3 +Y 5 ) +Y 6 -z, g 2 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 ) = - 1 λ 2 logU 2 +Y 5 +Y 6 -z, ∂ θ g(u 1 , u 2 ; θ , y, u 3 , x 4 , x 5 , x 6 ) = log u 1 /λ 2 1 , 0 T , ∂ y g(u 1 , u 2 ; θ , y, u 3 , x 4 , x 5 , x 6 ) = -1, 1 T .
The Jacobian matrix and its inverse are

J g (u 1 , u 2 ; θ , y) = - 1 λ 1 u 1 0 0 1 λ 2 u 2 , J -1 g (u 1 , u 2 ; θ , y) = - λ 1 u 1 0 0 λ 2 u 2 , ∂ u 1 J g (u 1 , u 2 ; θ , y) = 1 λ 1 u 2 1 0 0 0 , ∂ u 2 J g (u 1 , u 2 ; θ , y) = 0 0 0 1 λ 2 u 2 2 , ∂ θ J g (u 1 , u 2 ; θ , y) = 1 λ 2 1 u 1 0 0 0 and r 1 (u 1 , u 2 ; θ ) = 1 λ 1 u 1 log u 1 , r 2 (u 1 , u 2 ; θ ) = 0, d(u 1 , u 2 ; θ ) = 1 λ 1 (log u 1 + 1), r1 (u 1 , u 2 ; θ ) = λ 1 u 1 , r2 (u 1 , u 2 ; θ ) = λ 2 u 2 , d(u 1 , u 2 ; θ ) = -λ 1 -λ 2
Then the GLR estimator of ∂ F(y; θ )/∂ θ is given by

G 1 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 ) = 1 λ 1 -Y 1 1 {max(Y 1 +Y 4 +Y 6 ,Y 2 +Y 5 +Y 6 ,Y 1 +Y 3 +Y 5 +Y 6 ) ≤ y} ,
and the GLR estimator of ∂ F(y; θ )/∂ y is given by

G 2 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 ) = λ 1 1 {max (Y 4 +Y 6 ,Y 2 +Y 5 +Y 6 ,Y 3 +Y 5 +Y 6 ) ≤ y} + λ 2 1 {max (Y 1 +Y 4 +Y 6 ,Y 5 +Y 6 ,Y 3 +Y 5 +Y 6 ) ≤ y} -(λ 1 + λ 2 )1 {max(Y 1 +Y 4 +Y 6 ,Y 2 +Y 5 +Y 6 ,Y 1 +Y 3 +Y 5 +Y 6 ) ≤ y} .
Furthermore, we have

E[G 1 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 )] = E[E[G 1 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 )|Y 1 , . . . ,Y 5 ]] =E 1 λ 1 -Y 1 Φ 1 σ 6 log[(y -Y ) + ] -µ 6 ,
and

E[G 2 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 )] = E[E[G 2 (U 1 ,U 2 ; θ , y,U 3 , X 4 , X 5 , X 6 )|Y 1 , . . . ,Y 5 ]] =λ 1 E Φ 1 σ 6 log[(y -max(Y 4 ,Y 2 +Y 5 ,Y 3 +Y 5 )) + ] -µ 6 + λ 2 E Φ 1 σ 6 log[(y -max(Y 1 +Y 4 ,Y 5 ,Y 1 +Y 3 +Y 5 )) + ] -µ 6 -(λ 1 + λ 2 )E Φ 1 σ 6 log[(y -Y ) + ] -µ 6 ,
where Y = max(Y 1 +Y 4 ,Y 2 +Y 5 ,Y 1 +Y 3 +Y 5 ). We can take the expressions inside the expectations on the right-hand sides of the above equation as conditional GLR estimators. From this derivation, the three quantile sensitivity estimators of the form (2) using GLR, CGLR, and CGLRQ to estimate distribution sensitivities are referred to as GLR-1, CGLR-1, and CGLRQ-1, respectively. On the other hand, we have

F(y; θ ) = E[1 {Y 4 +Y 1 +Y 6 -y ≤ 0} 1 {Y 5 + max(Y 2 ,Y 1 +Y 3 ) +Y 6 -y ≤ 0}],
and we can also let (X 4 , X 5 ) be the input random variables in the stochastic model of [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF], leading to

g 1 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 ) = exp (µ 4 + σ 4 X 4 ) +Y 1 +Y 6 -y, g 2 (X 4 , X ; θ , y,U 1 ,U 2 ,U 3 , X 6 ) = exp (µ 5 + σ 5 X 5 ) + max(Y 2 ,Y 1 +Y 3 ) +Y 6 -y, ∂ θ g(x 4 , x 5 ; θ , y, u 1 , u 2 , u 3 , x 6 ) = log u 1 /λ 2 1 1, 1{log u 1 /λ 1 + log u 3 /λ 3 < log u 2 /λ 2 } T , ∂ y g(x 4 , x 5 ; θ , y, u 1 , u 2 , u 3 , x 6 ) = -1, 1 T , ∇ log f (X 4 ,X 5 ) (x 4 , x 5 ) = -x 4 , x 5 T .
The Jacobian matrix and its inverse are

J g (x 4 , x 5 ) = σ 4 e µ 4 +σ 4 x 4 0 0 σ 5 e µ 5 +σ 5 x 5 , J -1 g (x 4 , x 5 ) = 1 σ 4 e -µ 4 -σ 4 x 4 0 0 1 σ 5 e -µ 5 -σ 5 x 5 , ∂ x 4 J g (x 4 , x 5 ) =
σ 2 4 e µ 4 +σ 4 x 4 0 0 0 , ∂ x 5 J g (x 4 , x 5 ) = 0 0 0 σ 2 5 e µ 5 +σ 5 x 5 , and

d(x 4 , x 5 ; θ ) = log u 1 λ 2 1 1 + x 4 σ 4 e -µ 4 -σ 4 x 4 + 1 log u 1 λ 1 + log u 3 λ 3 < log u 2 λ 2 1 + x 5 σ 5 e -µ 5 -σ 5 x 5 , d(x 4 , x 5 ; θ ) = -1 + x 4 σ 4 e -µ 4 -σ 4 x 4 -1 + x 5
σ 5 e -µ 5 -σ 5 x 5 .

Then the GLR estimator of ∂ F(y; θ )/∂ θ is given by

G 1 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 ) = -1 {Y ≤ z} Y 1 λ 1 1 + X 4 σ 4 1 Y 4 + 1 Y 1 +Y 3 > Y 2 1 + X 5 σ 5 1 Y 5 . the GLR estimator of ∂ F(y; θ )/∂ y is given by G 2 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 ) = -1 {Y ≤ z} 1 + X 4 σ 4 1 Y 4 + 1 + X 5 σ 5 1 Y 5 . Furthermore, we have E[G 1 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 )] = E[E[G 1 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 )|Y 1 , . . . ,Y 5 ]] = -E Y 1 λ 1 1 + X 4 σ 4 1 Y 4 + 1 Y 1 +Y 3 > Y 2 1 + X 5 σ 5 1 Y 5 Φ 1 σ 6 log[(y -Y ) + ] -µ 6 .
Table 2:

α = 0.1 α = 0.3 α = 0.5 α = 0.7 α = 0.9 GLR-1
-0.456 ± 1.0 × 10 -4 -0.697 ± 1.3 × 10 -4 -0.907 ± 2.0 × 10 -4 -1.136 ± 3.9 × 10 -4 -1.374 ± 1.6 × 10 -3 CGLR-1 -0.456 ± 6.8 × 10 -5 -0.696 ± 1.0 × 10 -4 -0.907 ± 1.7 × 10 -4 -1.136 ± 3.6 × 10 -4 -1.374 ± 1.5 × 10 -3 CGLRQ-1 -0.456 ± 5.0 × 10 -5 -0.696 ± 5.6 × 10 -5 -0.907 ± 6.8 × 10 -5 -1.135 ± 1.1 × 10 -4 -1.377 ± 3.1 × 10 -4 GLR-2 -0.456 ± 1.0 × 10 -4 -0.700 ± 6.7 × 10 -4 -0.911 ± 1.1 × 10 -3 -1.163 ± 2.9 × 10 -3 -0.094 ± 7.4 × 10 -1 CGLR-2 -0.456 ± 3.0 × 10 -4 -0.698 ± 5.4 × 10 -4 -0.913 ± 1.0 × 10 -3 -1.161 ± 2.5 × 10 -3 -1.026 ± 2.5 × 10 -1 CGLRQ-2 -0.456 ± 2.6 × 10 -4 -0.697 ± 4.4 × 10 -4 -0.910 ± 7.9 × 10 -4 -1.145 ± 1.7 × 10 -3 -1.957 ± 2.0 × 10 -1 FDC(0.01) -0.453 ± 9.5 × 10 -4 -0.690 ± 9.8 × 10 -4 -0.897 ± 1.2 × 10 -3 -1.123 ± 1.6 × 10 -3 -1.366 ± 3.3 × 10 -3 FDC(0.1) -0.423 ± 2.9 × 10 -4 -0.639 ± 3.0 × 10 -4 -0.826 ± 3.6 × 10 -4 -1.024 ± 5.0 × 10 -4 -1.223 ± 1.0 × 10 -3

and

E[G 2 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 )] = E[E[G 2 (X 4 , X 5 ; θ , y,U 1 ,U 2 ,U 3 , X 6 )|Y 1 , . . . ,Y 5 ]] = -E 1 Y 4 + X 4 σ 4 Y 4 + 1 Y 5 + X 5 σ 5 Y 5 Φ 1 σ 6 log[(y -Y ) + ] -µ 6 .
We can take the expressions inside the expectations on the right-hand sides of the above equations as conditional GLR estimators. From this derivation, the three quantile sensitivity estimators of the form (2) using GLR, CGLR, and CGLRQ to estimate distribution sensitivities are referred to as GLR-2, CGLR-2, and CGLRQ-2, respectively.

In the experiment, we set λ i = 1, µ i = 0, σ i = 1, i = 1, 2, 3, and test five probability levels at α = 0.1, α = 0.3, α = 0.5, α = 0.7, and α = 0.9. From Table 2, we can see that the quantile sensitivity estimates by GLR-1, CGLR-1, CGLRQ-1, GLR-2, CGLR-2, CGLRQ-2 are statistically indistinguishable, FDC(0.01) appears to be slightly biased, and FDC(0.1) exhibits a significant bias compared to the GLR method. The variances of GLR-1, CGLR-1, CGLRQ-1 are smaller than GLR-2, CGLR-2, CGLRQ-2, respectively. The variances of all quantile sensitivity estimators go up as the probability levels grow up, and in particular, the variances of GLR-2, CGLR-2, CGLRQ-2 go up dramatically when the probability levels grow up from α = 0.7 to α = 0.9. The variances of GLR-1, CGLR-1, CGLRQ-1 are smaller than those of FDC(0.1) and FDC(0.01), and the CMC and RQMC lead to significant variance reduction.

CONCLUSIONS

We combine GLR methods in [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF] and [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF] with CMC and RQMC to reduce the variance of quantile sensitivity estimation, and apply the proposed methods to a toy example and a SAN example. In both examples, the GLR method in [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF] leads to a better performance than the GLR method in [START_REF] Peng | A New Unbiased Stochastic Derivative Estimator for Discontinuous Sample Performances with Structural Parameters[END_REF], and the variance reduction by CMC and RQMC is significant.

However, the variance reduction by RQMC for estimator (2) after smoothing the GLR estimators for two distribution sensitivities by CMC are not as substantial as that for distribution sensitivity estimators themselves (see [START_REF] Peng | Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs[END_REF]. The reason may be due to that the quantile estimator by order statistics is discontinuous with respect to the input random variables. This discontinuity issue can be treated by smoothing the empirical distribution function by certain kernel method, which could be a future work.

ACKNOWLEDGMENTS

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grants 71901003 and 72022001, by the Air Force Office of Scientific Research under Grant FA95502010211, and by Discovery Grant RGPIN-2018-05795 from NSERC-Canada.

AUTHOR BIOGRAPHIES

YIJIE PENG is an Assistant Professor in Guanghua School of Management, with affiliate faculty appointments in the Institute of Artificial Intelligence and National Institute of Health Data Science, all at Peking University. His research interests include stochastic modeling and analysis, simulation optimization, machine learning, data analytics, and healthcare. He is a member of INFORMS and IEEE, and serves


Yijie Peng

Michael C Fu

Pierre L'ecuyer

Bruno Tuffin

Variance Reduction for Generalized Likelihood Ratio Method By Conditional Monte Carlo and Randomized Quasi-Monte Carlo

The generalized likelihood ratio (GLR) method is a recently introduced gradient estimation method for handling discontinuities for a wide scope of sample performances. We put the GLR methods from previous work into a single framework, simplify regularity conditions for justifying unbiasedness of GLR, and relax some of those conditions that are difficult to verify in practice. Moreover, we combine GLR with conditional Monte Carlo methods and randomized quasi-Monte Carlo methods to reduce the variance. Numerical experiments show that the variance reduction could be significant in various applications.

Introduction

Simulation is widely used for analyzing and optimizing complex stochastic systems [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. Specifically, we first generate simple input random variables X = (X1, . . . , Xn), and then simulate a complex output random variable Y by evaluating the performance function ϑ(•) at the input random variables, i.e., Y = ϑ(X), which is also known as the sample performance. Depending on the application, ϑ(•) could take complicated functional forms. In stochastic activity networks used in project management, the performance function is nonlinear, because we are interested in the duration of the critical path, i.e., the activity path taking the longest duration; in a queueing system, the performance function can be estimated using a recursive equation when we are interested in waiting times of customers following certain dynamic mechanism; in financial options, the performance may be discontinuous due to thresholds determining whether the option is in the money or not.

Stochastic gradient estimation has been an important topic studied in simulation for decades [START_REF] Fu | Stochastic gradient estimation[END_REF], because it can be used for sensitivity analysis and optimization in complex stochastic models. Specifically, suppose the input random variables X(θ) and performance function ϑ(•; θ) are parameterized by θ, and we aim to estimate the derivative of the expectation of the output random variable, i.e, ∂E[Y (θ)]/∂θ. When θ is multi-dimensional, the gradient is a vector of derivatives with respect to parameters in each dimension. By the definition of derivative, ∂E[Y (θ)]/∂θ = lim δ→0 E[Y (θ + δ) -Y (θ)]/δ, and the finite difference (FD) estimator, (Y (θ + δ) -Y (θ))/δ, is straightforward and always implementable, and it treats the simulation model as a black-box. However, FD requires simulating one more sample path for each dimension of the parameter, and it suffers from a bias-variance tradeoff, i.e., the balance between choosing large δ to reduce the variance and using small values of δ to make the bias low. These undesirable properties are especially apparent when the dimension of the parameters is high.

In contrast, single-run stochastic gradient estimation techniques may provide unbiased estimates for the gradient with respect to all parameters simultaneously in the process of estimating the sample performance. Infinitesimal perturbation analysis (IPA) and the likelihood ratio (LR) method are two classic single-run unbiased stochastic gradient estimation techniques [START_REF] L'ecuyer | A unified view of the IPA, SF, and LR gradient estimation techniques[END_REF], [START_REF] Ho | Discrete Event Dynamic Systems and Perturbation Analysis[END_REF], [START_REF] Rubinstein | Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method[END_REF]. A measuredvalue differentiation can be viewed as a technique to reduce the variance of LR at a cost of extra simulations [START_REF] Heidergott | Weak differentiability of product measures[END_REF]. IPA is a pathwise derivative estimate obtained from The interchange of derivative and expectation is usually jus-tified by the dominated convergence theorem, which typically requires continuity and almost sure differentiability of the sample performance Y (θ) = ϑ(X; θ) with respect to θ.

∂E[Y (θ)] ∂θ = E ∂Y (
A detailed theoretical discussion on IPA can be found in [START_REF] Glasserman | Gradient Estimation via Perturbation Analysis[END_REF]. The LR estimator on the other hand can be obtained from

∂E[Y (θ)] ∂θ = ∂ ∂θ ϑ(x)fX (x; θ)dx = ϑ(x) ∂fX (x; θ) ∂θ dx = E ϑ(X) ∂ log fX (x; θ) ∂θ ,
where fX (•; θ) is the joint density of the input random variables. To justify unbiasedness of the estimator, LR does not require continuity for the sample performance, but it can only estimate derivatives with respect to a distributional parameter, i.e., a parameter in the distribution of input random variables, but not a structural parameter, i.e., a parameter directly appearing in function ϑ(•; θ). Detailed theoretical discussion on LR can be found in [START_REF] L'ecuyer | Note: On the interchange of derivative and expectation for likelihood ratio derivative estimators[END_REF]. A common belief in the simulation literature is that IPA usually has smaller variance than LR when they both apply. This may not be true in general (see, e.g., [START_REF] Cui | On the variance of single-run unbiased stochastic derivative estimators[END_REF] for a counterexample), but has been substantiated in [START_REF] Cui | On the variance of single-run unbiased stochastic derivative estimators[END_REF] under a sufficient condition. A hybrid of IPA and LR (IPA-LR) in [START_REF] L'ecuyer | A unified view of the IPA, SF, and LR gradient estimation techniques[END_REF] can be obtained by

∂E[Y (θ)] ∂θ = ∂ ∂θ ϑ(x; θ)fX (x; θ)dx =E ∂ϑ(x; θ) ∂θ x=X + ϑ(X; θ) ∂ log fX (x; θ) ∂θ .
Structural parameter is allowed for IPA-LR, but it requires differentibility of function ϑ(•; θ) with respect to the structural parameter.

The single-run unbiased stochastic gradient estimation techniques are considered in traditional backgrounds including discrete event systems [2], and risk management in financial engineering [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF], [START_REF] Hong | Estimating quantile sensitivities[END_REF], [START_REF] Hong | Simulating sensitivities of conditional value at risk[END_REF], [START_REF] Liu | Kernel estimation of the Greeks for options with discontinuous payoffs[END_REF], [START_REF] Wang | A new stochastic derivative estimator for discontinuous payoff functions with application to financial derivatives[END_REF], [3], [START_REF] Hong | Estimating sensitivities of portfolio credit risk using Monte Carlo[END_REF], [START_REF] Jiang | Technical note -On estimating quantile sensitivities via infinitesimal perturbation analysis[END_REF], [START_REF] Heidergott | A[END_REF], [START_REF] Glynn | Computing sensitivities for distortion risk measures[END_REF]. Recently, this topic has attracted attention in artificial intelligence and machine learning [START_REF] Mohamed | Monte Carlo gradient estimation in machine learning[END_REF], [START_REF] Peng | A new likelihood ratio method for training artificial neural networks[END_REF], where the dimension of the parameter is typically extremely high, so that the single-run and unbiasedness properties are particularly helpful in gradient-based optimization. Backpropagation (BP), the most popular gradient estimation technique for training artificial neural networks (ANN), is shown in [START_REF] Peng | A new likelihood ratio method for training artificial neural networks[END_REF] to be pathwise equivalent to IPA, and backward propagation of errors can reduce computational complexity. In addition, [START_REF] Peng | A new likelihood ratio method for training artificial neural networks[END_REF] shows that an LR-based method can significantly improve the robustness of ANN.

In simulation, discontinuities often appear in the sample performance with respect to the structural parameter, so that IPA and LR fail to be unbiased. For example, the derivative of the distribution function ∂F (z)/∂z = ∂E[1{Y ≤ z}]/∂z has a discontinuous sample performance with respect to z and the derivative and expectation cannot be interchanged. Smoothed perturbation analysis (SPA) and push-out LR can be used to address the discontinuity issue [START_REF] Fu | Conditional Monte Carlo: Gradient Estimation and Optimization Applications[END_REF], [START_REF] Rubinstein | Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method[END_REF], but SPA has to choose what to be conditioned on and push-out LR requires an explicit transformation to push the parameter into the density of input distribution, which are problem dependent. Peng et al. [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] propose a generalized LR (GLR) method to systematically treat sensitivity analysis for a large class of sample performances with discontinuities. In contrast to SPA and push-out LR, the GLR estimator has an analytical form without the need for conditioning and transformation. In [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], the unbiasedness of GLR is justified by a set of conditions, including that the tails of the input distribution go smoothly to zero fast enough, which excludes exponentially distributed random variables and uniform random numbers, for example. This smoothness requirement is relaxed in [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF], where the inputs of the stochastic model are uniform random numbers, which are the basic building blocks for generating other random variables. GLR is shown to generalize both LR and pushout LR.

In this work, we put the GLR methods in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] and [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF] under a common framework. Then, by adopting a similar technique developed in [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF], we simplify the assumptions required for justifying the unbiasedness of GLR in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] and further relax some of the conditions that are difficulty to verify in practice. Although the GLR method has broad applicability, previous work indicates that it may inherit some of the undesirable variance properties of the LR method. We address this issue by combining GLR with conditional Monte Carlo (CMC) methods and randomized quasi-Monte Carlo (RQMC) methods.

CMC methods can reduce the variance and smooth the performance function in simulation by conditioning on certain events or random variables and then integrating out the remaining randomness [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. For an estimator H(Z), we have

E[H(Z)] = E[ H( Z)], (1)
where

H( Z) := E[H(Z)| Z]
with Z being a part of input random variables in Z. The variance reduction for the conditional estimator H( Z) can be seen from the following variance decomposition formula:

Var(H(Z)) = Var(E[H(Z)| Z]) + E[Var(H(Z)| Z)] ≥ Var( H( Z)).
Typically, H( Z) is smoother than H(Z), due to the integration taken in the conditional expectation. SPA uses CMC to smooth the sample performance, after which IPA is applied to differentiate the conditional expectation. GLR does not need smoothing to obtain an unbiased derivative estimator, but CMC can be applied afterward to reduce the variance for GLR. The relation between applications of CMC in SPA and GLR is illustrated in Figure 1. RQMC methods replace the vectors of uniform random numbers that drive independent simulation runs by dependent vectors of uniform random numbers that cover the space more evenly. When estimating an expectation, they can provide an unbiased estimator with a variance that converges to zero at a faster rate than with Monte Carlo [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Method[END_REF],

Figure 1: Applications of CMC to obtain SPA and conditional GLR. [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF]. Such a faster rate can be proved when the estimator inside the expectation is sufficiently smooth as a function of the underlying uniform random numbers. When the estimator is not smooth (e.g., discontinuous), the convergence rate may not be improved, but RQMC could still reduce the variance by a constant factor. We show, through numerical experiments, that the variance of the GLR estimator can be significantly reduced by appropriately combining GLR with CMC and RQMC. Similar type of combination of CMC and RQMC for reducing the variance of quantile estimation can be found in [START_REF] Nakayama | Quantile estimation via a combination of conditional monte carlo and randomized quasi-monte carlo[END_REF].

The rest of the paper is organized as follows. Section 2 introduces the GLR method. Variance reduction for GLR by CMC and RQMC is discussed in Section 3. Section 4 exemplifies the applications of the method to a stochastic activity network, a single-server queue, and a barrier option, with numerical experiments on them presented in Section 5. The last section offers conclusions.

Generalized Likelihood Ratio Method

We consider the sensitivity analysis problem of estimating

∂ ∂θ E[ϕ(g(X; θ))], (2)
where ϕ : R n → R is a measurable function not necessarily continuous, g(x; θ) = (g1(x; θ), . . . , gn(x; θ)) T is a vector of functions with sufficient smoothness for x ∈ R n , and X = (X1, . . . , Xn) is a vector of input random variables with a joint density fX (x; θ) supported on Ω ⊆ R n . From [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], the dimension of g is required to be smaller or equal to the dimension of x in deriving a GLR estimator. For simplicity of theoretical discussion, we require them to be the same. In application, X does not have to include all the input random numbers given to a simulation model. Instead, X can be only a subset of input random variables in simulation and we can condition on the remaining input random variables outside of X when deriving the GLR estimators. This leaves us freedom to select input random variables in X for deriving GLR estimators as long as the conditions to ensure unbiasedness introduced later can be satisfied, and some choices of X could lead to simpler forms and more desirable variance properties for GLR. The Jacobian of g(•; θ) is

Jg(x; θ) :=       ∂g 1 (x;θ) ∂x 1 ∂g 1 (x;θ) ∂x 2 • • • ∂g 1 (x;θ) ∂xn ∂g 2 (x;θ) ∂x 1 ∂g 2 (x;θ) ∂x 2 • • • ∂g 2 (x;
• • • ∂gn(x;θ) ∂xn      
.

In this work, we require the function g(.; θ) to be one-to-one (invertible) so that the Jacobian will be invertible. That is, g must represent a valid multivariate change of variable. However, we only need the existence of the inversion or change of variable for g rather than the explicit form. This requirement can be relaxed to the extent that that g and its the Jacobian are only invertible locally [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF]. Let ei be the i-th unit vector. We define

d(x; θ) := -trace J -1 g (x; θ)∂ θ Jg(x; θ) + n i=1 e T i J -1 g (x; θ) (∂x i Jg(x; θ)) J -1 g (x; θ)∂ θ g(x; θ) -J -1 g (x; θ)∂ θ g(x; θ) T ∇x log fX (x; θ).
From [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], when the tails of fX (•; θ) go smoothly to zero fast enough, then, under certain regularity conditions, it can be shown that an unbiased GLR estimator for (2) is given by G(X; θ) := ϕ(g(X; θ))w(X; θ),

w(x; θ) := ∂ log fX (x; θ) ∂θ + d(x; θ), (3)
where ∂yh(y) is the matrix obtained by differentiating h with respect to y element-wise.

Peng et al. [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF] consider the case when X is a vector of uniform random numbers U = (U1, . . . , Un) such that Ω = (0, 1) n and ∂ log fX (x; θ)/∂θ = 0. Then under certain regularity conditions, we have the following unbiased GLR estimator for (2):

G(U ; θ) := ϕ(g(U ; θ))d(U ; θ) + n i=1 ϕ(g(U i; θ))ri(U i; θ) -ϕ(g(U i ; θ))ri(U i ; θ) (4)
where

U i := (U1, . . . ,
:= J -1 g (x; θ) ∂ θ g(x; θ) T ei, i = 1, . . . , n.
In a special case when the input random variables are independent and non-parameterized with the tails of the distribution going smoothly to zero fast enough, from [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF], the GLR estimator in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] and the GLR estimator in [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF] coincide even though the input random variables or sample point (in probability space) are interpreted differently in two GLR methods. This is in contrast to the unified IPA and LR framework in [START_REF] L'ecuyer | A unified view of the IPA, SF, and LR gradient estimation techniques[END_REF], where two different interpretations of the input random variables or sample point lead to two distinctive estimators, i.e, the IPA and LR estimators.

Example 1. We use a simple density estimation problem to illustrate how to apply the two versions of the GLR method: We want to estimate

∂E[1{X + U ≤ z}] ∂z ,
where X is a standard normal random variable, and U is a uniform random variable U (0, 1) that is independent of X.

Although there are two input random variables X and U , we select one and condition on the other to apply GLR. In this example,

ϕ(•) = 1{• ≤ 0}. This gives n = 1 in the preceding development. If we select X, then g(x; z, u) = x + u -z, ∂g(x; z, u)/∂x = 1, ∂g(x; z, u)/∂z = -1, ∇x log fX (x) = -x
, and the other derivatives in (3) are zeros, so the GLR estimator can be given by G(X; z, U ) = -1{X + U ≤ z}X.

If we select U , then g(u; z, x) = x+u-z, ∂g(u; z, x)/∂u = 1, ∂g(u; z, x)/∂z = -1, and other derivatives in (4) are zeros, so the GLR estimator can be given by G(U ; z, X) = 1{X ≤ z} -1{X + 1 ≤ z} = 1{z -1 < X ≤ z}, which coincides with an SPA estimator derived in [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF] for the same problem. This example falls under the umbrella of distribution sensitivity estimation studied in both [START_REF] Peng | Maximum likelihood estimation by Monte Carlo simulation: Towards data-driven stochastic modeling[END_REF] and [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF], which estimates the derivatives of the distribution function

F (z; θ) = E[1{g(X; θ) ≤ z}]
with respect to both θ and z, with g(•; θ) being a function with sufficient smoothness.

As a particular stochastic activity network (SAN) example, the output could be the maximum of the durations of activities on different paths. We then consider the following distribution function:

E 1 max i=1,...,n gi(X; θ) ≤ z = E n i=1 1 {gi(X; θ) -z ≤ 0} ,
and the distribution sensitivities can be estimated by the GLR method using general frameworks (3) and ( 4) with

ϕ(y) = n i=1 1{yi ≤ 0}.
For the SAN example, gi would be the duration of activities on the i-th path, which will be analyzed in Section 4.

Example 2. For 0 ≤ α ≤ 1, the α-quantile (also known as value-at-risk) of a random variable g(U ; θ) with cdf F (•; θ) is defined as

qα(θ) := arg min{z : F (z; θ) ≥ α}. When F (•; θ) is continuous, qα(θ) = F -1 (α; θ). Let U (j) , j = 1, . . . , m, be i.i.d. realizations of U ∼ U (0, 1) d ,
and Fm(•) be the empirical distribution of g(U (j) ; θ), j = 1, . . . , m. The empirical α-quantile F -1 m (α) is the inverse of the empirical distribution evaluated at α. This empirical quantile satisfies the following central limit theorem [START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF]:

√ m F -1 m (α) -qα(θ) d → N 0, α(1 -α) f (qα(θ); θ) . ( 5
)
We can estimate the asymptotic variance by

mα(1 -α) m j=1 G (j) (z, θ)| z= F -1 m (α)
, where G (j) (z, θ) is a GLR estimator for estimating the density using the same realizations of uniform random variables U (j) as in quantile estimator F -1 m (α). From [START_REF] Fu | Conditional Monte Carlo estimation of quantile sensitivities[END_REF], we can also apply distribution sensitivity estimators to estimate the quantile sensitivity

dqα(θ) dθ = - ∂F (z; θ) ∂θ z=qα(θ) f (qα(θ); θ),
where the numerator and denominator can be estimated by GLR.

We then simplify the regularity conditions in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] and relax some of them that are difficult to verify in practice by adapting a similar technique as developed in [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF]. In particular, we introduce the following conditions to justify the unbiasedness of the GLR estimation in the form (3).

(A.1) The function g(•; θ) is invertible and the density fX (x; θ) is supported on R n ; g(x; θ) is twice continuously differentiable and fX (x; θ) is continuously differentiable with respect to (x, θ) ∈ R n × Θ, where Θ is a compact neighborhood of parameter θ.

(A.

2) The following limiting condition holds: lim

x i →±∞ R n sup θ∈Θ ri(x; θ)fX (x; θ) j =i dxj = 0, i = 1, . . . , n. (A.
3) The following integrability condition holds:

R n sup θ∈Θ ϕ(g(x; θ))w(x; θ)fX (x; θ) dx < ∞.
Remark 1. Condition (A.1) requires that the function g and density fX be sufficiently smooth. Condition (A.2) requires that the tails of fX (•; θ) go to zero fast enough. In contrast to the conditions in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] for establishing unbiasedness of the GLR estimator, conditions (A.1)-(A.3) avoid certain integrability condition imposed on some intermediate smoothed function of ϕ in the proof, which could be harder to verify in practice.

Theorem 1. Under conditions (A.1)-(A.3), ∂E[ϕ(g(X; θ))] ∂θ = E[ϕ(g(X; θ))w(X; θ)],
where w(•) is defined by (3).

The proof the theorem can be found in the appendix. The idea is to first smooth ϕ(•), which may not be continuous, then apply integration by parts to move differentiation from ϕ(•) to other smoother terms, and finally take limit to establish the unbiasedness of the final GLR estimator.

Variance Reduction

In this section, we discuss how to apply CMC and RQMC for reducing the variance of the GLR estimators.

Conditional Monte Carlo Method

The GLR estimator can be combined with CMC for reducing its variance. This will also be a key transformation to improve RQMC accuracy. We want to find an appropriate Z to condition on such that the conditional expectation H( Z) in (1) becomes smoother and can be computed efficiently. To illustrate, we consider a special case of function ϕ(•):

ϕ(y) = n i=1 1{yi ∈ Yi}h(y),
where h(•) is a smooth function of y = (y1, . . . , ym), and Yi is a set. This special case of ϕ(•) covers the examples of distribution sensitivities discussed later. Then the GLR estimator (3) can be written as n i=1 1{gi(X; θ) ∈ Yi}W (X; θ), where W (x; θ) := h(g(x; θ))w(x; θ), and the GLR estimator (4) can be written as

n i=1 n j=1 1{gj(U i; θ) ∈ Yi}Ri(U i; θ) - n j=1 1{gj(U i ; θ) ∈ Yi}Ri(U i ; θ) + n i=1 1{gi(Ui; θ) ∈ Yi}D(U ; θ), where Ri(u; θ) := ri(u; θ)h(g(u; θ)), D(u; θ) := d(u; θ)h(g(u; θ)).
Let Z be a vector of generic random variables, which could either be X or U , and Q(•) be a generic function which could either be W (•) or D(•). To smooth the sample performance, suppose there exists i * ∈ {1, . . . , n} such that for i = 1, . . . , n,

1{gi(Z; θ) ∈ Yi} = 1{Zi * ∈ Zi(Z-i * ; θ)},
where z-j = (z1, . . . , zj-1, zj+1, . . . , zn), and Zi(•) is a set depending on the argument. Here we only integrate out one random variable, and there are discussions on how to integrate more than one random variable for some applications in [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF], which could lead to smaller variance. We have

E n i=1 1{gi(Z; θ) ∈ Yi}Q(Z; θ) = E [1{Zi * ∈ ∩ n i=1 Zi(Z-i * )}Q(Z; θ)] = E 1{zi * ∈ Zi(Z-i * ; θ)}Q(zi * , Z-i * )fi * (zi * )dzi * = E P Zi * ∈ ∩ n i=1 Zi(Z-i * ) Z-i * ×E Q(Z; θ) Zi * ∈ ∩ n i=1 Zi(Z-i * ), Z-i * .
If, in addition, Q(Z; θ) is independent of Zi * , then the equation above simplifies to

E P Zi * ∈ ∩ n i=1 Zi(Z-i * ) Z-i * Q(Z; θ) .
When Yi is an interval (ai, bi) and gi(•) is strictly increasing with zi * , we have

Zi(z-i * ; θ) = (g -1 i (ai; z-i * ), g -1 i (bi; z-i * )).
For the GLR estimator -1{X + U ≤ z}X in Example 1 of Section 2 for estimating the density, a conditional GLR estimator (z -X)1{z -1 < X < z} can be obtained by

E[-1{X + U ≤ z}X] = E[E[-1{U ≤ z -X}X|X]] = -E[(z -X)1{z -1 < X < z}].

Randomized

Quasi-Monte Carlo Method

Quasi-Monte Carlo (QMC) refers to a class of deterministic numerical integration methods in which the integrand is evaluated at a fixed set of m points, and the average is used as an approximation. One limitation of the method is that it is very hard to estimate the approximation error in practice. RQMC takes the QMC points and randomizes them in a way that each point has the uniform distribution over (0, 1) n , so that each randomized point represents a proper realization of U while the set of m points still covers the unit hypercube (0, 1) n more uniformly than typical independent random points (so the points are not independent) [START_REF] L'ecuyer | Randomized quasi-Monte Carlo: An introduction for practitioners[END_REF].

In general, for a given function h, RQMC estimates the integral µ = (0,1) n h(u)du by the average

μm := 1 m m j=1 h(U (j) ),
where U (1) , . . . , U (m) form an RQMC point set. The most common types of QMC point set constructions are lattice rules, polynomial lattices rules, and digital nets [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF], [START_REF] L'ecuyer | Quasi-Monte Carlo methods with applications in finance[END_REF].

For lattice rules, an appropriate randomization is a random shift modulo 1, which adds a single uniform random point to all the lattice points, and retains the shifted points that are in (0, 1) n as the m RQMC points. This randomization preserves the lattice structure, and there are explicit expressions for Var[μm] in terms of the Fourier coefficients of h, and computable bounds on this variance for certain classes of smooth functions [START_REF] L'ecuyer | Variance reduction via lattice rules[END_REF], [START_REF] L'ecuyer | Randomized quasi-Monte Carlo: An introduction for practitioners[END_REF], [START_REF] L'ecuyer | Algorithm 958: Lattice builder: A general software tool for constructing rank-1 lattice rules[END_REF]. When the mixed derivatives of h are sufficiently smooth, the variance can converge at a faster rate than O(m -1 ), sometimes nearly O(m -2 ) and even faster in some cases. When h is not smooth (e.g., discontinuous), these convergence rate results do not apply, although weaker results do apply [START_REF] He | On the convergence rate of randomized quasi-Monte Carlo for discontinuous functions[END_REF], and even when the convergence rate is not improved, the variance is often reduced by a constant factor. For polynomial lattices rules and digital nets in general, which include Sobol' points, the random shift does not preserve the structure and net properties, but other appropriate randomizations do, including nested uniform scrambling, some affine scrambles, and random digital shifts. Variance bounds and convergence rate results are available for these, as well [START_REF] Dick | Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration[END_REF], [START_REF] L'ecuyer | Quasi-Monte Carlo methods with applications in finance[END_REF].

We now discuss how to combine GLR with RQMC. Our model formulation (4) in terms of a function of independent U (0, 1) random variables makes it an obvious candidate for the application of RQMC, which is designed exactly for this type of formulation. In our setting, we can apply RQMC by taking h as the GLR derivative estimator G( • ; θ) to obtain the RQMC estimator

Ḡm := 1 m m j=1 G U (j) ; θ ,
in which {U (1) , . . . , U (m) } is an RQMC point set. For formulation (3), suppose X can be generated by Γ(U ), and an RQMC estimator can be obtained by

Ḡm := 1 m m j=1
G Γ(U (j) ); θ .

For example, when X1, . . . , Xn are independent random variables with marginal distribution functions F1, . . . , Fn, then they can be generated by Xi = F -1 i (Ui), i = 1, . . . , n. Note that with RQMC, the terms in the sum of Ḡm are not independent, so one cannot estimate the variance of Ḡm through a straightforward application of the sample variance, as in standard Monte Carlo. To estimate the RQMC variance, one can make, say l, independent randomizations of the QMC point set, to obtain l independent replicates of Ḡm, and compute the sample variance of these l replicates [START_REF] L'ecuyer | Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications[END_REF]. This could be used to compute a confidence interval on the true derivative, although one must be careful, because the distribution of Ḡm does not always converge to a normal distribution when m → ∞ with RQMC [START_REF] L'ecuyer | On the distribution of integration error by randomly-shifted lattice rules[END_REF]. The overall RQMC estimator in this case will be

Ḡm,l := 1 l l =1 Ḡ( ) m ,
where Ḡ( ) m is the th independent replicate of Ḡm, and Var[ Ḡm,l ] is estimated by

1 (l -1)l l =1 Ḡ( ) m -Ḡm,l 2 .
For the α-quantile sensitivity estimation discussed in Example 2 of Section 2, the cdf F (•; θ) of g(U ; θ) can be estimated by its empirical RQMC counterpart

Fm(z) := 1 m m j=1 1{g(U (j) ; θ) ≤ z},
where {U (1) , . . . , U (m) } is an RQMC point set, and the quantile qα(θ) can be estimated by the pseudo-inverse F -1 m (α). With l independent randomizations of the RQMC points, we can average the l independent randomizations of this empirical RQMC cdf,

F m, (z) := 1 m m j=1 1{g U (j, ) ; θ ≤ z}, Fm,l (z) := 1 l l =1 F m, (z),
and estimate qα(θ) by F-1 m,l (α), which is a consistent estimator and satisfies the central limit theorem:

√ l F-1 m,l (α) -qα(θ) d → N 0, Var (F m, (qα(θ))) f 2 (qα(θ)) .
The asymptotic variance can be estimated by

1 l-1 l =1 (F m, (z) -Fm,l (z)) 2
Ḡm,l (z) where Ḡm,l (z) is a GLR estimator for density using RQMC. Kaplan et al. [START_REF] Kaplan | Randomized quasi-Monte Carlo for quantile estimation[END_REF] point out that using

2 z= F -1 m,l (α)
1 l l =1 F -1 m, (α) instead of F-1
m,l (α) would not provide a consistent estimator of qα(θ) as l → ∞. Similarly, we can also estimate the quantile sensitivity by GLR using RQMC.

Applications

In this section, we discuss the applications of GLR and CMC to estimating distribution sensitivities for a SAN, single-server queue, and barrier option. Distribution sensitivities can also be applied to calibrate parameters and optimize quantile-based risk measures [START_REF] Peng | Maximum likelihood estimation by Monte Carlo simulation: Towards data-driven stochastic modeling[END_REF], [START_REF] Hu | A stochastic approximation method for simulation-based quantile optimization[END_REF].

Stochastic Activity Network

We estimate a simple SAN studied in [START_REF] Heidergott | Sensitivity estimation for Gaussian systems[END_REF], which is depicted in Figure 2. Distribution sensitivity estimation by SPA for a SAN with a different structure can be found in [START_REF] Fu | Sensitivity analysis in monte carlo simulation of stochastic activity networks[END_REF] and [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF]. As in [START_REF] Heidergott | Sensitivity estimation for Gaussian systems[END_REF], the duration of the seventh activity is assumed to be deterministic. There are four different paths representing the tasks to reach the final stage of a project, i.e., [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF][START_REF] Cui | On the variance of single-run unbiased stochastic derivative estimators[END_REF][START_REF] Fu | Sensitivity analysis in monte carlo simulation of stochastic activity networks[END_REF], (2, 5, 6), (1, 3, 5, 6), [START_REF] Fu | Stochastic gradient estimation[END_REF].

The completion time for the entire project is max(Y1 + Y4 + Y6, Y2 + Y5 + Y6, Y1 + Y3 + Y5 + Y6, Y7), and the sample performance for the distribution function of completion time is

Y = 1 {max(Y1 + Y4 + Y6, Y2 + Y5 + Y6, Y1 + Y3 + Y5 + Y6, Y7) ≤ z} =1 {Y1 + Y4 + Y6 ≤ z} 1 {Y2 + Y5 + Y6 ≤ z} × 1 {Y1 + Y3 + Y5 + Y6 ≤ z} 1{Y7 ≤ z}.
Unlike [START_REF] Heidergott | Sensitivity estimation for Gaussian systems[END_REF], where the durations of the first six activities all follow independent normal distributions, we assume that the first three activities follow independent exponential distributions: Yi = -1 λ i log(Ui), i = 1, 2, 3, and the other three activities follow independent log-normal distributions, Yi = exp (µi + σiXi), i = 4, 5, 6. We note that

1{Y1 + Y4 + Y6 ≤ z}1{Y2 + Y5 + Y6 ≤ z} = 1{Y1 + max(Y4, Y3 + Y5) + Y6 ≤ z} and 1{Y2 + Y5 + Y6 ≤ z}1{Y1 + Y3 + Y5 + Y6, Y7 ≤ z} = 1{Y5 + max(Y2, Y1 + Y3) + Y6 ≤ z}.
For z ≤ Y7, the distribution function F (z) of the completion time Y is equal to zero, and for z > Y7, the distribution function is

F (z) = E[1 {Y1 + max(Y4, Y3 + Y5) + Y6 -z ≤ 0} × 1 {Y2 + Y5 + Y6 -z ≤ 0}] = E[1 {Y4 + Y1 + Y6 -z ≤ 0} × 1 {Y5 + max(Y2, Y1 + Y3) + Y6 -z ≤ 0}].
To estimate the density f (z) = ∂ ∂z F (z), we can view (U1, U2) as the input random variables in the stochastic model ϕ(g(U ; θ)), and we have

g1(U1, U2; z) = - 1 λ1 log U1 + max(Y4, Y3 + Y5) + Y6 -z, g2(U1, U2; z) = - 1 λ2 log U2 + Y5 + Y6 -z, ∂zg(u1, u2; z) = -1, 1 T .
The Jacobian matrix and its inverse are Jg(u1, u2; z) = -

1 λ 1 u 1 0 0 1 λ 2 u 2 , J -1 g (u; z) = - λ1u1 0 0 λ2u2 , and
r1(u1, u2; z) = λ1u1, r2(u1, u2; z) = λ2u2, d(u1, u2; z) = -λ1 -λ2.
We emphasize that the input random variables in the stochastic model ϕ(g(U ; θ)) must be selected carefully in a way that the Jacobian of g is invertible. Also, in the example, the three paths are merged into two by some transformation, so we end up with a 2×2 Jacobian. But we could also have kept the three paths and selected (U1, U2, U3) as the input, which is more direct. So there are many valid choices for what we take as "input" when applying the method. For the SAN, the authors of [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF] select a minimal cut and compute the density on the completion time conditional on the other links (not in the cut) as a density estimator. For GLR, we can take the lengths of the links of the minimal cut as the input random variables. For z > Y7, the GLR estimator of f (z) is given by

G(U1, U2; z) = λ11 {max (Y4 + Y6, Y2 + Y5 + Y6, Y3 + Y5 + Y6) ≤ z} + λ21 {max (Y1 + Y4 + Y6, Y5 + Y6, Y3 + Y5 + Y6) ≤ z} -(λ1 + λ2)1 Y ≤ z ,
where Y = max(Y1 +Y4, Y2 +Y5, Y1 +Y3 +Y5). Furthermore, we have

E[G(U 1 , U 2 ; z)] = E[E[G(U 1 , U 2 ; z)|Y 1 , . . . , Y 5 ]] =λ 1 E Φ 1 σ 6 log[(z -max(Y 4 , Y 2 + Y 5 , Y 3 + Y 5 )) + ] -µ 6 +λ 2 E Φ 1 σ 6 log[(z -max(Y 1 + Y 4 , Y 5 , Y 1 + Y 3 + Y 5 )) + ] -µ 6 -(λ 1 + λ 2 )E Φ 1 σ 6 log[(z -Y ) + ] -µ 6 .
On the other hand, we can also let (X4, X5) be the input random variables in the stochastic model ϕ(g(X; θ)), leading to g1(X4, X5; z) = exp (µ4 + σ4X4) + Y1 + Y6 -z, g2(X4, X5; z) = exp (µ5 + σ5X5)

+ max(Y2, Y1 + Y3) + Y6 -z, ∂zg(x4, x5; z) = -1, 1 T , ∇ log f (X 4 ,X 5 ) (x4, x5) = -x4, x5 T .
The Jacobian matrix and its inverse are

Jg(x4, x5; z) = σ4e µ 4 +σ 4 x 4 0 0 σ5e µ 5 +σ 5 x 5 , J -1 g (x4, x5; z) = 1 σ 4 e -µ 4 -σ 4 x 4 0 0 1 σ 5 e -µ 5 -σ 5 x 5
, and

d(x4, x5; z) = -1 + x4 σ4 e -µ 4 -σ 4 x 4 -1 + x5 σ5 e -µ 5 -σ 5 x 5 .
Again, the input random variables in the stochastic model ϕ(g(X; θ)) must be selected carefully in a way that the Jacobian of g is invertible. For z > Y7, the GLR estimator of f (z) is given by

G(X4, X5; z) = -1 {Y ≤ z} 1 + X4 σ4 1 Y4 + 1 + X5 σ5 1 Y5 .
Furthermore, we have

E[G(X4, X5; z)] = E[E[G(X4, X5; z)|Y1, . . . , Y5]] = -E 1 Y4 + X4 σ4Y4 + 1 Y5 + X5 σ5Y5 ×Φ 1 σ6 log[(z -Y ) + ] -µ6 .
We can take the expression inside the expectation as a conditional GLR estimator.

Single-Server Queue

We consider distribution sensitivity estimation for the waiting time of the customers in a single-server first-come-firstserved queue depicted in Figure 3. Density estimation for this single-server queue by SPA and RQMC can be found in [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF]. When the i-th customer arrives, he/she may need to wait if the system time (waiting time plus service time), for the (i -1)-th customer is longer than the interarrival time between the i-th customer and (i -1)-th customer. Otherwise, the waiting time of the i-th customer is zero, i.e., the waiting time of customers follows the Lindley equation:

Wi = max{0, Wi-1 + Si-1 -Ai}, i ≥ 2,
where Wi and Si are the waiting time and service time of the i-th customer, and Ai is the interarrival time between the i-th customer and the (i -1)-th customer. Suppose Si = exp (µ1 + σ1Xi) and Ai = exp (µ2 + σ2X i ) follow independent log-normal distributions, where Xi and X i follow the standard normal distribution for i ≥ 1. We can also represent the service time as Si = exp µ1 + σ1Φ -1 (Ui) so that uniform random numbers could be viewed as input random variables for deriving the GLR estimator. Let µ1 = θ.

The density f (z; θ) of Wn over (0, ∞) can be written as the derivative of the distribution function F (z; θ) of Wn:

∂F (z; θ) ∂z = ∂ ∂z E[1{Wn ≤ z}] = ∂ ∂z E[1{Wn-1 + Sn-1 -An ≤ z}].
If we let

g(Xn-1, z) = Wn-1 + exp (µ1 + σ1Xn-1) -An -z, then ∂g(xn-1, z) ∂xn-1 = σ1 exp (µ1 + σ1xn-1) , ∂ 2 g(xn-1, z) ∂x 2 n-1 = σ 2 1 exp (µ1 + σ1xn-1) , ∂ log fX n-1 (xn-1) ∂xn-1 = -xn-1,
and from [START_REF] Peng | Maximum likelihood estimation by Monte Carlo simulation: Towards data-driven stochastic modeling[END_REF], the GLR estimator for the density is

G1(Xn-1; z) = -1 {Wn ≤ z} Xn-1 + σ1 Sn-1σ1 , z > 0.
To estimate the derivative with respect to θ, let gi(Xi; θ) = θ + σ1Xi, i = 1, . . . , n -1.

Then from [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], the GLR estimator for ∂F (z; θ)/∂θ is given by

G2(X; θ) = 1 σ1 1 {Wn ≤ z} n-1 i=1
Xi, z > 0, which coincides with the classic LR estimator. The GLR estimators in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF] and [START_REF] Peng | Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[END_REF] also coincide for this example, because the input random variables are independent and non-parameterized, with the tails of the distribution going smoothly to zero fast enough (standard normal in this example). From [START_REF] Peng | Maximum likelihood estimation by Monte Carlo simulation: Towards data-driven stochastic modeling[END_REF], higher-order distribution sensitivities can also be obtained by the GLR method, e.g.,

∂ 2 F (z; θ) ∂θ∂z = E 1 {Wn ≤ z} 1 - n-1 i=1 Xi (Xn-1 + σ1) Sn-1σ1 .
By noticing that

{Wn ≤ z} = X n ≥ 1 σ2 log (Wn-1 + Sn-1 -z) + -µ2 ,
we have

E [G1(Xn-1; z)] = -E E 1 {Wn ≤ z} Xn-1 + σ1 Sn-1σ 2 1
Wn-1, Xn-1

= -E 1 -Φ 1 σ2 log (Wn-1 + Sn-1 -z) + -µ2 × Xn-1 + σ1 Sn-1σ 2 1 ,
which offers a conditional GLR estimator for the density, and

E [G2(X; θ)] = E E 1 σ1 1 {Wn ≤ z} n-1 i=1 Xi Wn-1, Xn-1, . . . , X1 = 1 σ1 E 1 -Φ 1 σ2 log (Wn-1 + Sn-1 -z) + -µ2 × n-1 i=1
Xi , which offers a conditional GLR estimator for the distribution sensitivity with respect to θ. Furthermore,

∂ 2 F (z; θ) ∂z∂θ = E E 1 {Wn ≤ z} 1 - n-1 i=1 Xi (Xn-1 + σ1) Sn-1σ1
Wn-1, Xn-1, . . . , X1

= E 1 -Φ 1 σ2 log (Wn-1 + Sn-1 -z) + -µ2 × 1 - n-1 i=1 Xi (Xn-1 + σ1) Sn-1σ1 .
The expression inside the expectation offers a conditional GLR estimator for the second-order distribution sensitivity.

Barrier Option

We consider the distribution sensitivity estimation for the payoff of a barrier option. An up-and-out (knockout) barrier option is worthless if the path of the underlying asset exceeds a barrier L. The event when the barrier option stays "alive" is {maxi=1,...,n St i < L}, where St ∈ R is the underlying asset price at time t. In Figure 4, the price curve on the top breaches the barrier prior to expiration T , so that its payoff is zero. Suppose St = S0 exp{(r -σ 2 /2)t + σBt + N (t) j=1 Jj} follows a geometric jump-diffusion process, where S0 is the initial underlying asset price, r is the interest rate, σ is the implied volatility, {N (t)} is a counting process, and Ji, i ∈ Z + , are the jump sizes. Let S0 = θ. Assume discrete monitoring, and ∆ is the step size of the discrete monitoring points ti = i∆, i = 1, . . . , n, T = n∆. The barrier option would have a positive payoff if it stays "alive" and at T the price of the underlying asset is above the strike price K. For a European barrier option, the payoff is

e -rn∆ (Sn -K) n-1 i=1 1{Si ≤ L}1{K < Sn < L}, where Si = θ exp i r - σ 2 2 ∆ + σ √ ∆ i j=1 Xj + N i =1 J i, ,
where Xi = (Bi∆ -B (i-1)∆ )/ √ ∆, i = 1, . . . , n, which are i.i.d. standard normal random variables, Ni = N (i∆) -N (i(∆ -1)), and J i, is the -th jump in the ith period.

For 0 ≤ z ≤ e -rn∆ (L -K), we can represent the distribution function by

F (z; θ) = E n-1 i=1 1{gi(X; θ) ≤ 0}1{gn(X; θ, z) ≤ 0} + n i=1 i-1 j=1 1{gj(X; θ) ≤ 0}1{gi(X; θ) > 0} = 1 -E n-1 i=1 1{gi(X; θ) ≤ 0} (1 {gn(X; θ) ≤ 0} -1 {gn(X; θ, z) ≤ 0})] , where gi(X; θ) = log θ + i r - σ 2 2 ∆ + σ √ ∆ i j=1 Xj + i j=1 N j =1 J j, -log L, i = 1, . . . , n,
gn(X; θ, z) = gn(X; θ) + log L e rn∆ z + K .

From [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], the GLR estimators for f (z; θ) and ∂F (z; θ)/∂θ are given respectively by

G1(X; z) = - n-1 i=1 1{gi(X; θ) ≤ 0}1 {gn(X; θ, z) ≤ 0} × e rn∆ Xn σ √ ∆(e rn∆ z + K) , and
G2(X; θ) = - n-1 i=1 1{gi(X; θ) ≤ 0} × (1 {gn(X; θ) ≤ 0} -1 {gn(X; θ, z) ≤ 0}) X1 θσ √ ∆ . Define T (l) i = 1 σ √ ∆   log(L/θ) -i r - σ 2 2 ∆ - i j=1 N j =1 J j,   - i j=1,j =l Xj, i = 1, . . . , n, T (l) n = 1 σ √ ∆ log[(e rm∆ z + K)/θ] -n r - σ 2 2 ∆ - n j=1 N j =1 J j,   - n j=1,j =l Xj.
We have

E[E[G 1 (X; z)|X 2 , . . . , Xn]] = - e rn∆ σ √ ∆(e rn∆ z + K) ×E E 1 X 1 ≤ min T (1)
n , min i=1,...,n-1 T

(1) i Xn X 2 , . . . , Xn = -

e rn∆ σ √
∆(e rn∆ z+K)

E Φ min T (1)
n , min i=1,...,n-1 T

and

E[E[G2(X; θ)|X2, . . . , Xn]] = - 1 θσ √ ∆ E E 1 min T (1) n , min i=1,...,n-1 T (1) i < X1 ≤ min i=1,...,n T (1) i X1 X2, . . . , Xn = - 1 θσ √ ∆ E min j=1,...,n T (1) j min T (1)
n ,min i=1,...,n-1 T

(1) i xφ(x)dx .

To avoid integration in the estimator, we can also condition on other random variables, which gives a less smooth estimator than that derived above:

E[E[G 2 (X; θ)|X 1 , X 3 , . . . , Xn]] = - 1 θσ √ ∆ E E 1 X 1 ≤ T (1) 1 × 1 min T (2) n , min i=2,...,n-1 T (2) i < X 2 ≤ min i=2,...,n T (2) i ×X 1 X 1 , X 3 , . . . , Xn = - 1 θσ √ ∆ E 1 X 1 ≤ T (1) 1 X 1 × Φ min i=2,...,n-1 T (2) i -Φ min T (2) n , min i=2,...,n-1 T (2) i
.

Numerical Experiments

We report numerical experiments on the three applications discussed in the previous section. GLR is compared with the finite difference method (Y (θ + δ) -Y (θ))/δ using common random numbers (FDC(δ)) to generate Y (θ + δ) and Y (θ). The GLR method together with CMC is called conditional GLR (CGLR), and CGLR together with RQMC is denoted as CGLRQ. Both CGLR and CGLRQ are compared with SPA and SPA when combined with RQMC (SPAQ). For RQMC, we use the Sobol sequence scrambled by the algorithm of [START_REF] Matousek | On the l2-discrepancy for anchored boxes[END_REF] in Matlab. Derivations of SPA for the density estimation in three applications are similar to those in [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF] and can be found in the appendix. We set the sample size as m = 2 13 for the standard Monte Carlo and RQMC estimators, and their variances are estimated by l = 100 independent experiments.

Stochastic Activity Network

The parameters for the stochastic models described in Section 4.1 are set as λi = 1, i = 1, 2, 3, µj = 0 and σj = 1, j = 4, 5, 6. We estimate the density function ∂F (z)/∂z for z ∈ (0, 20). Figure 5(a) shows the density curve estimated by GLR, and Figure 5(b) presents the coverage rate curve of the 90% confidence intervals based on the asymptotic normality result given by ( 5) for α-quantiles with α = 0.1 × i, i = 1, . . . , 9, by m = 2 13 samples and l = 100 independent experiments. Figure 6(a) shows a variance comparison between the two GLR estimators in Section 4.1. GLR-1 is derived using a stochastic model with uniform random numbers as inputs, and GLR-2 is derived using a stochastic model with normal random variables as inputs. Since the variance of GLR-1 is much smaller than GLR-2, GLR-1 is used for this example in the rest of the experiments, and we simply denote the method by GLR. Figure 6(b) shows that the variance of GLR is much smaller than FDC(0.01) and FDC(0.1), and the variance of FDC increases significantly when the perturbation size becomes smaller. The reason FDC works poorly is that the sample performance is discontinuous. If CMC is first applied to smooth the sample performance, then conditional FDC can achieve a comparable variance to SPA [START_REF] L'ecuyer | On the convergence rates of IPA and FDC derivative estimators[END_REF]. Figure 7(a) shows a variance comparison between GLR, CGLR, and CGLRQ. We can see that both CMC and RQMC substantially reduce the variance of GLR. From Figure 7(b), RQMC also reduces the variance of SPA, and CGLRQ and SPAQ achieve comparable magnitudes of variance.

Single-Server Queue

The parameters for the stochastic models described in Section 4.2 are set as µi = 0 and σi = 1, j = 4, 5, 6, and n = 10. We estimate distribution sensitivities ∂F (z; θ)/∂z, ∂F (z; θ)/∂θ, and ∂ 2 F (z; θ)/∂z∂θ for z ∈ (0, 30). Figure 8 presents the three distribution sensitivity curves estimated by GLR.

From Figure 9(a), FDC(0.01) has a much larger variance than GLR for estimating the second-order distribution sensitivity ∂ 2 F (z; θ)/∂z∂θ, and FDC suffers from a bias- variance tradeoff. In this example, the variances of GLR, CGLR, and CGLRQ cannot be distinguished statistically. Derivation of SPA for estimating the second-order distribution sensitivity is not straightforward due to the discontinuities in the SPA estimator for the first-order derivative (see the appendix).

From Figure 10(a), we can see that GLR and CGLR achieve comparable performance, whereas the variance of CGLRQ is much smaller than GLR and CGLR. The peak value on the variance curve of CGLRQ is about 1/6 of that of GLR and CGLR. Figure 10(b) shows that SPA and SPAQ have comparable variance, and the variance of CGLRQ is smaller than SPA and SPAQ when z ≤ 15. The variance of CGLRQ becomes larger than SPA and SPAQ when z > 15, because the SPA and SPAQ estimators go to zero as z → ∞, whereas the CGLRQ estimator converges to a random variable with a constant variance as z → ∞.

Barrier Option

The parameters for the stochastic models described in Section 4.3 are set as T = 10, n = 10, σ = 1, r = 0.001, ∆ = 1, K = 100, L = 120, and S0 = 100. The counting process is assumed to be a Poisson process with intensity 1, and the jump size is assumed to follow a normal distribution with mean 0 and standard deviation 0.01. We estimate distribution sensitivities ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ for z ∈ (0, 30). estimate ∂F (z; θ)/∂θ, we apply the GLR estimator in Section 4.3 that does not involve integration. Figure 11 presents the two distribution sensitivity curves estimated by GLR. Figure 12 presents the variance comparison between GLR and FDC for two first-order distribution sensitivities. FDC(0.01) has a much larger variance than GLR. The variances of GLR and FDC(0.1) are comparable, but GLR is unbiased and FDC(0.1) is biased. The bias estimate of FDC(0.01) and FDC(0.1) can be found in the appendix. Figure 13 presents the variance comparison between GLR, CGLR, and CGLRQ for two first-order distribution sensitivities. The variance of CGLR is slightly smaller than that of GLR, whereas the variance of CGLRQ is substantially smaller than those of CGLR and GLR. The peak value on the variance curve of CGLRQ is about 1/3 of that of GLR and CGLR.

A Appendix

Proof of Theorem 1

Proof. Proof. As in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], define a sequence of bounded functions ϕL(x) = max{min{ϕ(x), L}, -L}, and then |ϕL(x)| ≤ ϕ(x) and limL→∞ ϕL(x) = ϕ(x). From Theorem 1 in [START_REF] Peng | A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[END_REF], there exists a sequence of bounded and smooth functions ϕ matrix. Then we have

∂ ∂θ [-M,M ] n ϕ ,L(g(x; θ))fX (x; θ)dx = [-M,M ] n ϕ ,L(g(x; θ)) ∂ log f (x; θ) ∂θ +∇yϕ ,L(y)|y=g(x;θ)∂θ g(x; θ) fX (x; θ)dx = [-M,M ] n ϕ ,L(g(x; θ)) ∂ log f (x; θ) ∂θ +∇xϕ ,L(g(x; θ))J -1 g (x; θ)∂ θ g(x; θ) fX (x; θ)dx.
The interchange of the differentiation and integration can be justified by the dominated convergence theory. Under condition (A.1), fX (x; θ) is continuously differentiable and g(x; θ) is twice continuously differentiable in R n ×Θ, so their function values and derivatives are bounded in a compact space [-M, M ] n × Θ. By the Gauss-Green Theorem,

[-M,M ] n ∇xϕ ,L(g(x; θ))J -1 g (x; θ)∂ θ g(x; θ) × fX (x; θ)dx1 • • • dxn = n i=1 [-M,M ] n-1 ϕ ,L(g(u; θ)) J -1 g (x; θ)∂ θ g(x; θ) T × eifX (x; θ) j =i dxj M x i =-M - [-M,M ] n ϕ ,L(g(x; θ)) × div J -1 g (x; θ)∂ θ g(x; θ)fX (x; θ) dx1 • • • dxn,
where for h(x) = (h1(x), . . . , hn(x)), By differentiating equation J -1 g (x; θ)Jg(x; θ) = I with respect to xi on both sides, we have

div(h(x)) := n i=1 ∂h(x) ∂xi . Then div J -1 g (x; θ)∂ θ g(x; θ)fX (x; θ) = n i=1 ∂ ∂xi e T i J -1 g (x; θ)∂ θ g(x; θ)fX (x; θ) = n i=1 e T i ∂x i J -1 g (x; θ) ∂ θ g(x; θ) +trace(J -1 g (x; θ)∂ θ Jg(x; θ)) fX (x; θ) + n i=1 J -1 g (x; θ)∂ θ g(x; θ) T ei ∂fX (x; θ) ∂xi .
0 =∂x i J -1 g (x; θ)Jg(x; θ) =∂x i J -1 g (x; θ)Jg(x; θ) + J -1 g (x; θ)∂x i Jg(x; θ),
which leads to

∂x i J -1 g (x; θ) = -J -1 g (x; θ) (∂x i Jg(x; θ)) J -1 g (x; θ).
Therefore, we have d(x; θ) = -div J -1 g (x; θ)∂ θ g(x; θ)fX (x; θ) /fX (x; θ). With the discussion above,

∂ ∂θ [-M,M ] n ϕ ,L(g(x; θ))fX (x; θ)dx = n i=1 [-M,M ] n-1 ϕ ,L(g(x; θ))ri(x; θ)fX (x; θ) j =i dxj M x i =-M + [-M,M ] n ϕ (g(x; θ))w(x; θ)fX (x; θ)dx. Under condition (A.2), lim →0 lim M →∞ sup θ∈Θ n i=1 [-M,M ] n-1 ϕ ,L(g(x; θ))ri(x; θ) × fX (x; θ) j =i dxj M x i =-M = 0.
By change of variables and Hölder's inequality, lim

→0 sup θ∈Θ [-M,M ] n (ϕ ,L(g(x; θ)) -ϕL(g(x; θ))) × w(x; θ)fX (x; θ)dx = lim →0 sup θ∈Θ S M (ϕ ,L(y) -ϕL(y))| det(Jg(x; θ))|w(x; θ) × fX (x; θ)| x=g -1 (y;θ) dy ≤ lim →0 ϕ ,L -ϕL p × sup θ∈Θ S M | det(Jg(x; θ))w(x; θ)fX (x; θ)| q x=g -1 (y;θ) dy 1/q = 0, (6)
where where the first term goes to zero because of ( 6) and the second term goes to zero because of [START_REF] Fu | Stochastic gradient estimation[END_REF]. From [START_REF] Rudin | Principles of Mathematical Analysis[END_REF],

SM := {y ∈ R n : y = g(x; θ), x ∈ [-M, M ] n }. With condition (A.3), lim M →∞ sup θ∈Θ [-M,M ] n ϕL(g(x; θ))w(x; θ)fX (x; θ)dx - (0,1) n ϕL(g(u; θ))w(x; θ)fX (x; θ)dx ≤ lim M →∞ (0,1) n \[-M,M ] n sup θ∈Θ ϕ(g(x; θ))w(x; θ)fX (x; θ) dx = 0. ( 7

SPA for SAN

There are many possible choices of conditioning random variables to derive SPA for density estimation, and we only use one possible choice. Certain conditions are needed to ensure unbiasedness for SPA. More general discussions of SPA for density estimation can be found in [START_REF] L'ecuyer | Monte Carlo and quasi-Monte Carlo density estimation via conditioning[END_REF]. The SPA estimator can be obtained by

∂ ∂z E[1{Y ≤ z}] = ∂ ∂z E[E[1{Y ≤ z}|Y1, . . . , Y5]] = ∂ ∂z E Φ 1 σ6 log[(z -Y ) + ] -µ6 = E 1{z ≥ Y } σ6(z -Y ) + φ 1 σ6 log[(z -Y ) + ] -µ6 .

SPA for Single-Sever Queue

The SPA estimator for the density is Figure 14 presents the variance comparison between CGLRQ, SPA, and SPAQ for estimating ∂F (z; θ)/∂θ in the barrier option example. We can see that SPA and SPAQ have a smaller variance than CGLRQ.

∂F (z; θ) ∂z = ∂ ∂z E[E[1{Wn ≤ z}|W n-1 , X n-1 , . . . , X 1 ]] = ∂ ∂z E 1 -Φ 1 σ 2 log (W n-1 + S n-1 -z) + -µ 2 = 1 σ 2 E 1{W n-1 + S n-1 ≥ z} (W n-1 + S n-1 -z) +
Figure 15 presents the bias of FDC(0.01) and FDC(0.1) for ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ in the barrier option example. The true distribution sensitivities are estimated by SPAQ with the sample size as m = 2 13 and l = 10 4 independent experiments. Bias of FDC are also estimated with the sample size as m = 2 13 and l = 10 4 independent experiments. From Figure 15, we can see that FDC(0.01) has a smaller bias than FDC(0.1). [3] Chen, N., Liu, Y., 2014. American option sensitivities estimation via a generalized infinitesimal perturbation analysis approach. Operations Research 62, 616-632. measure-valued differentiation approach to sensitivity analysis of quantiles. Mathematics of Operations Research 41, 293-317.

Figure 2 :

 2 Figure 2: A SAN with seven activities.

Figure 3 :

 3 Figure 3: A single-server queue.

Figure 4 :

 4 Figure 4: A European up-and-out barrier call option.

  Function curve of ∂F (z)/∂z.

Figure 5 :

 5 Figure 5: Density estimation by GLR and coverage rates of 90% confidence intervals for α-quantiles in the SAN example.

  Variance curves of GLR and FDC.

Figure 6 :

 6 Figure 6: Variances of density estimation by GLR and FDC in the SAN example.

  Variance curves of GLR and FDC.

Figure 7 :

 7 Figure 7: Variances of density estimation by GLR, CGLR, CGLRQ, SPA, and SPAQ in the SAN example.

  respect to z Distribution senstivities with respect to Second-order distribution sensitvities

Figure 8 :

 8 Figure 8:Estimation of distribution sensitivities ∂F (z; θ)/∂z, ∂F (z; θ)/∂θ, and ∂ 2 F (z; θ)/∂z∂θ by GLR in the queueing example.

  Variance curves of GLR and FDC. Variance curves of GLR, CGLR, and CGLRQ.

Figure 9 :

 9 Figure 9: Variances of estimation for ∂ 2 F (z; θ)/∂z∂θ by GLR, CGLR, CGLRQ and FDC in the queueing example.

  ,L(•) such that lim L→∞ ϕ ,L -ϕL p = 0, where p > 1, andh p := R n |h(x)| p dx 1/p . Since g(x; θ)is an invertible function vector, its Jacobian is an invertible Variance curves of GLR and FDC. Variance curves of GLR, CGLR, and CGLRQ.

Figure 10 :

 10 Figure 10: Variances of estimation for ∂F (z; θ)/∂θ by GLR, CGLR, CGLRQ, SPA and SPAQ in the queueing system example.

Figure 11 :

 11 Figure 11:Estimation of distribution sensitivities ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ by GLR in the barrier option example.

  Variance curves of GLR and FDC. Variance curves of GLR, CGLR, and CGLRQ.

Figure 12 :

 12 Figure 12: Variances of estimation for ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ by GLR and FDC in the barrier option example.

  Variance curves of GLR, CGLR, and CGLRQ. Variance curves of GLR, CGLR, and CGLRQ.

Figure 13 :

 13 Figure 13: Variances of estimation for ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ by GLR, CGLR, and CGLRQ in the barrier option example.

  [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]-δ] n ϕ ,L(g(u; θ))w(x; θ)fX (x; θ)dx -R n ϕL(g(x; θ))w(x; θ)fX (x; θ)dx≤ lim M →∞ lim →0 sup θ∈Θ [-M,M ] n ϕ ,L(g(x; θ))w(x; θ)fX (x; θ)dx -[-M,M ] n ϕL(g(x; θ))w(x; θ)fX (x; θ)dx + lim M →∞ sup θ∈Θ [-M,M ] n ϕL(g(x; θ))w(x; θ)fX (x; θ)dx -R nϕL(g(x; θ))w(x; θ)fX (x; θ)dx = 0,

d

  dθ lim →0 h (θ) = lim →0 h (θ) holds if h (θ) converges uniformly with respect to θ ∈ Θ as → 0. Therefore,∂ ∂θ R n ϕL(g(x; θ))fX (x; θ)dx = lim M →∞ lim →0 ∂ ∂θ [-M,M ] n ϕL(g(x; θ))fX (x; θ)dx = R n ϕL(g(u; θ))w(x; θ)fX (x; θ)du.With condition (A.3) and noticing that |ϕL(x) -ϕ(x)| ≤ ϕ(x),lim L→∞ sup θ∈Θ R n(ϕL(g(x; θ)) -ϕ(g(x; θ)))× w(x; θ)fX (x; θ)dx = 0.By this uniform convergence,∂ ∂θ R n ϕ(g(x; θ))fX (x; θ)dx = ∂ ∂θ R n lim L→∞ ϕL(g(x; θ))fX (x; θ)dx = lim L→∞ ∂ ∂θ R n ϕL(g(x; θ))fX (x; θ)dx = lim L→∞ R n ϕL(g(x; θ))w(x; θ)fX (x; θ)dx = R nϕ(g(x; θ))w(x; θ)fX (x; θ)dx, which proves the theorem.

×φ 1 σ 2 ≥ 1 . 1 T

 211 log (W n-1 + S n-1 -z) + -µ 2 , and the SPA estimator for ∂F (z; θ)/∂θ is ∂F (z;SPA for Barrier OptionTo derive the SPA estimators for f (z; θ) and ∂F (z; θ)/∂θ,F (z; θ) = 1 -E E n-1 i=1 1{gi(X; θ) ≤ 0} × (1 {gn(X; θ) ≤ 0} -1 {gn(X; θ, z) ≤ 0}) |X2, . . . , Xn]] rn∆ z + K)E φ min T(1) n , min i=1,...,n-

Figure 14 :

 14 Figure 14: Variances of estimation for ∂F (z; θ)/∂θ by CGLRQ, SPA and SPAQ in the barrier option example.

  Bias curves of FDC for ∂F (z; θ)/∂z. Bias curves of FDC for ∂F (z; θ)/∂θ.

Figure 15 :

 15 Figure 15: Bias of estimation for ∂F (z; θ)/∂z and ∂F (z; θ)/∂θ by FDC(0.01) and FDC(0.1) in the barrier option example.

Conclusions

In this paper, we have studied the GLR methods in [42] and [41] under a single framework, simplified regularity conditions for justifying unbiasedness of GLR in [42] and relaxed some of the conditions that are difficult to verify in practice by adapting a similar technique developed in [41]. Moreover, we have discussed how to combine GLR with CMC and RQMC to reduce the variance of the resultant estimators, and applied them to estimate the distribution sensitivities for a SAN, a single-server queue, and a barrier option. Numerical results show that CMC and RQMC may reduce the variance of GLR substantially, and CGLRQ sometimes achieves even smaller variance than SPAQ.

Acknowledgement

This document is the results of the research project funded by National Natural Science Foundation of China (NSFC) under Grants 71901003 and 72022001, by the Air Force Office of Scientific Research under Grant FA95502010211, by Discovery Grant RGPIN-2018-05795 from NSERC-Canada.