Yijie Peng
email: pengyijie@pku.edu.cn

Michael C Fu The

Robert H Smith

Jiaqiao Hu
email: jqhu@ams.stonybrook.edu

Pierre L'ecuyer
email: lecuyer@iro.umontreal.ca

Bruno Tuffin
email: bruno.tuffin@inria.fr

Generalized Likelihood Ratio Method for Stochastic Models with Uniform Random Numbers As Inputs

Keywords: simulation, stochastic derivative estimation, discontinuous sample performance, uniform random numbers, generalized likelihood ratio method

) and Peng et al. (2020a) that precludes a direct formulation in terms of uniform random numbers as inputs. By overcoming this limitation, our new estimator greatly expands the applicability of the GLR method, which we demonstrate for several general classes of uniform input random numbers, including independent inverse transform random variates and dependent input random variables governed by an Archimedean copula. We show how the new derivative estimator works in specific settings such as density estimation, distribution sensitivity for quantiles, and sensitivity analysis for Markov chain stopping time problems, which we illustrate with applications to statistical quality control, stochastic activity networks, and credit risk derivatives. Numerical experiments substantiate broad applicability and flexibility in dealing with discontinuities in sample performance.

Introduction

Stochastic gradient estimation plays a central role in gradient-based optimization and sensitivity analysis [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. The finite difference (FD) method is easily implementable, but it must balance a bias-variance tradeoff and requires extra simulations. Infinitesimal perturbation analysis (IPA) and the likelihood ratio (LR) method are two well-established unbiased derivative estimation techniques [START_REF] Glynn | Likelihood ratio gradient estimation for stochastic systems[END_REF][START_REF] Ho | Discrete Event Dynamic Systems and Perturbation Analysis[END_REF][START_REF] Glasserman | Gradient Estimation via Perturbation Analysis[END_REF][START_REF] Rubinstein | Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method[END_REF][START_REF] Glynn | Likelihood ratio gradient estimation for stochastic recursions[END_REF]. IPA typically leads to lower variance than LR (L' Ecuyer, 1990;[START_REF] Cui | Technical note-On the variance of single-run unbiased stochastic derivative estimators[END_REF], and the weak derivative method reduces the variance of LR at the cost of performing extra simulations (Pflug, 1988, Heidergott and[START_REF] Heidergott | Weak differentiability of product measures[END_REF]. L 'Ecuyer (1990) provides a general framework unifying IPA and LR, under which the resulting estimator depends on the choice of what a sample point in a probability space represents, and could in particular be a hybrid between IPA and LR. See [START_REF] Fu | Stochastic gradient estimation[END_REF] for a recent review.

Traditional applications of stochastic gradient estimation are in discrete event dynamic systems (DEDS), including queueing systems [START_REF] Suri | Perturbation analysis gives strongly consistent sensitivity estimates for the M/G/1 queue[END_REF][START_REF] Fu | Second derivative sample path estimators for the GI/G/m queue[END_REF][START_REF] L'ecuyer | Stochastic optimization by simulation: Convergence proofs for the GI/G/1 queue in steady-state[END_REF], inventory management (Fu, 1994, Bashyam and[START_REF] Bashyam | Optimization of (s, S) inventory systems with random lead times and a service level constraint[END_REF], statistical quality control (Fu andHu, 1999, Fu et al., 2009b), maintenance systems [START_REF] Heidergott | Optimisation of a single-component maintenance system: a smoothed perturbation analysis approach[END_REF][START_REF] Heidergott | Gradient estimation for multicomponent maintenance systems with age-replacement policy[END_REF], and financial engineering and risk management, such as computing financial derivatives [START_REF] Fu | Sensitivity analysis for Monte Carlo simulation of option pricing[END_REF][START_REF] Broadie | Estimating security price derivatives using simulation[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][START_REF] Hong | Estimating sensitivities of portfolio credit risk using Monte Carlo[END_REF][START_REF] Chen | American option sensitivities estimation via a generalized infinitesimal perturbation analysis approach[END_REF][START_REF] Lei | Sensitivity analysis of portfolio credit derivatives by conditional Monte Carlo simulation[END_REF], value-at-risk (VaR) and conditional VaR (CVaR) (Hong, 2009[START_REF] Hong | Simulating sensitivities of conditional value at risk[END_REF], Fu et al., 2009a[START_REF] Jiang | Technical note -On estimating quantile sensitivities via infinitesimal perturbation analysis[END_REF][START_REF] Heidergott | A measure-valued differentiation approach to sensitivity analysis of quantiles[END_REF]. Recently, stochastic gradient estimation techniques have attracted attention in machine learning and artificial intelligence; see [START_REF] Mohamed | Monte Carlo gradient estimation in machine learning[END_REF] for a review paper written by a research team of Google's DeepMind. Peng et al. (2020b) show pathwise equivalence between IPA and backpropagation, and how the computational complexity for estimating the gradient is reduced by propagating the errors backwardly along the ANN. An LR-based method is then proposed to train ANNs, which can improve the robustness in classifying images under both adversarial attacks and natural noise corruptions.

IPA requires continuity in the sample performance, whereas LR does not directly apply for structural parameters (parameters directly appearing in the sample performance), which significantly limit their applicability. Smoothed perturbation analysis (SPA) deals with discontinuous sample performances by using a conditioning technique [START_REF] Gong | Smoothed perturbation analysis of discrete event dynamical systems[END_REF]Ho, 1987, Fu and[START_REF] Fu | Conditional Monte Carlo: Gradient Estimation and Optimization Applications[END_REF], but a good choice of conditioning is problem-dependent. Push-out LR addresses structural parameters by pushing the parameters out of the sample performance and into the density [START_REF] Rubinstein | Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method[END_REF], which can be achieved alternatively with the IPA-LR in L' Ecuyer (1990), but it requires an explicit transformation. Recently, Peng et al. (2018) proposed a generalized likelihood ratio (GLR) method that is capable of dealing with a large scope of discontinuous sample performances with structural parameters in a unified framework. The method extends the application domain of IPA and LR and does not require conditioning and transformation techniques tailored to specific problem structures.

The GLR method has the virtue of handling many applications in a uniform manner, and it has been used to deal with discontinuities in financial derivatives, statistical quality control, maintenance systems, and inventory systems [START_REF] Peng | On the regularity conditions and applications for generalized likelihood ratio method[END_REF], Peng et al., 2018). Distribution sensitivities, which mean the derivatives of the distribution function with respect to both the arguments and the parameters in the underlying stochastic model, lie at the center of many applications such as quantile sensitivity estimation, confidence interval construction for the quantile and quantile sensitivities, and statistical inference [START_REF] Peng | On the asymptotic analysis of quantile sensitivity estimation by Monte Carlo simulation[END_REF][START_REF] Lei | Applications of generalized likelihood ratio method to distribution sensitivities and steady-state simulation[END_REF]. Peng et al. (2020a) derive GLR estimators for any order of distribution sensitivities and apply them to maximum likelihood estimation for complex stochastic models without requiring analytical likelihoods. [START_REF] Glynn | Computing sensitivities for distortion risk measures[END_REF] apply the GLR method to estimate sensitivity of a distortion risk measure, which is a Lebesgue-Stieltjes integral of quantile sensitivities and includes Var and CVaR as special cases.

Although the existing GLR method has broad applicability, it requires that the density of the input distribution is known and that both tails of the density go down to zero smoothly and fast enough, which may not be satisfied in some applications, depending on what is interpreted as the input random variables. In this work, we relax this smoothness requirement and establish the unbiasedness of GLR gradient estimators for stochastic models whose inputs are uniform random numbers, which are the basic building blocks in generating other random variables. Unlike in Peng et al. (2018) where the surface integration part for the GLR estimator is zero, the surface integration part for the GLR estimator in the present work is not necessarily zero but can be estimated by simulation. If the surface integration part is zero, we are able to relax certain integrability conditions given in Peng et al. (2018) that are difficult to verify in practice.

We provide specific forms of the GLR estimators for two types of stochastic models and apply the GLR method to various problem settings, including distribution sensitivities, credit risk financial derivatives, and statistical quality control. The GLR estimator with independent input random variables generated from the inverse transform of uniform random numbers reduces to the classic LR estimator, which indicates that GLR is a generalization of LR from a different perspective than that of Peng et al. (2018). We also show how GLR can provide sensitivity estimators for models defined in terms of random vectors with given marginal distributions and whose dependence structures are specified by Archimedean copulas [START_REF] Nelsen | An Introduction to Copulas[END_REF]. The Gaussian copula has been widely used due to its simplicity [START_REF] Li | On default correlation: a copula approach[END_REF], and sensitivity analysis for portfolio credit risk derivatives with joint defaults governed by a Gaussian copula has been studied in [START_REF] Chen | Sensitivity estimates for portfolio credit derivatives using Monte Carlo[END_REF] using LR and SPA. However, the Gaussian copula was widely criticized after the 2008 financial crisis, because it underestimates the probability of joint defaults. Archimedean copulas [START_REF] Nelsen | An Introduction to Copulas[END_REF], not covered in [START_REF] Chen | Sensitivity estimates for portfolio credit derivatives using Monte Carlo[END_REF], are relatively easy to simulate and can better capture the asymmetric tail dependence structure of the joint default data [START_REF] Embrechts | Modelling dependence with copulas and applications to risk management[END_REF].

GLR can be used to estimate the distribution sensitivity functions, including the density and the quantile function together with confidence intervals, using a single batch of uniform random numbers. We provide numerical illustrations with various examples. Sample performances of control charts used in statistical quality control generally involve a stopping time defined by hitting a control limit (Fu andHu, 1999, Fu et al., 2009b). The IPA and LR methods do not apply to this problem due to the discontinuity with respect to the control limits. We formulate a stopping time problem with uniform random numbers as inputs in a Markov chain, and estimate its sensitivity by the GLR method. We also apply GLR to estimate distribution sensitivities for stochastic activity networks (SAN) with uniform random numbers as inputs. In practice, the duration of an activity in a SAN and the output sample in control charts may be supported on a compact space, so a distribution supported on the whole space, which was assumed in Peng et al. (2018) and Peng et al. (2020a), is unsuitable for modeling input distributions in these stochastic models.

The rest of the paper is organized as follows. Section 2 sets the framework. The GLR estimator is presented in Section 3 with the specific forms of the estimators for two types of models. Applications are given in Section 4. Numerical experiments can be found in Section 5. The last section concludes.

The technical proofs and additional numerical results can be found in the online appendix.

Problem Formulation

Consider a stochastic model of the following form:

ϕ(g(U ; θ)), (1)
where ϕ : R n → R is a measurable function (not necessarily continuous), g(•; θ) = (g 1 (•; θ), . . . , g n (•; θ)) is a vector of functions g i : (0, 1) n → R with certain smoothness properties to be made more precise shortly, and U = (U 1 , . . . , U n ) is a vector of i.i.d. U(0, 1) random variables (i.e., uniform over (0, 1)). For simplicity, we take θ as a scalar. When θ is a vector, each component of the gradient can be estimated separately using the method developed in this work. We consider the problem of estimating the following derivative:

∂E[ϕ(g(U ; θ))] ∂θ . (2)
A straightforward pathwise derivative estimator, i.e., IPA, obtained by directly interchanging derivative and expectation, may not apply because discontinuities in the sample performance of the stochastic model could be introduced by ϕ(•). In Peng et al. (2018), the stochastic model considered for the derivative estimation problem is ϕ(g(X; θ)), where the density of X = (X 1 , . . . , X n ) is assumed to be known and both tails go down to zero smoothly and fast enough. This assumption is not satisfied by discontinuous densities, such as the uniform and exponential distributions, and we want to address this limitation. Before deriving a GLR derivative estimator, we first introduce two examples to illustrate potential applications of the stochastic model (1).

Example 1 Independent Inputs Generated via the Inverse Transform Method. Suppose X = (X 1 , . . . , X n ) is a vector of independent random variables, where each X i has cumulative distribution function (cdf) F i (•; θ), i = 1, . . . , n, and is generated by (standard) inversion:

X i = F -1 i (U i ; θ), i = 1, . . . , n,
with i.i.d. U i ∼ U(0, 1). A stochastic model with i.i.d. U(0, 1) random numbers as input can be written as ϕ(g

(U ; θ)) = ϕ(F -1 1 (U 1 ; θ), . . . , F -1 n (U n ; θ)), where g(u; θ) = (F -1 1 (u 1 ; θ), . . . , F -1 n (u n ; θ)).
Example 2 Archimedean Copulas. Copulas are a general way of representing the dependence in a multivariate distribution. A copula is any multivariate cdf whose one-dimensional marginals are all U(0, 1). It can be defined by a function C(•; θ) : [0, 1] n → [0, 1] that satisfies certain conditions required for C to be a consistent cdf; see [START_REF] Nelsen | An Introduction to Copulas[END_REF]. For any given copula and arbitrary marginal distributions with continuous cdf's

F 1 (•), F 2 (•), ..., F n (•) with densities f i (•), i = 1, . . . , n,
one can define a multivariate distribution having exactly these marginals with joint cdf F X given by F X (x) = C(F 1 (x 1 ), F 2 (x 2 ), ..., F n (x n ); θ) for all x := (x 1 , . . . , x n ). [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF] shows that any multivariate distribution can be represented in this way. If C(•; θ) is absolutely continuous, the density of the joint distribution is

f X (x; θ) = c (F 1 (x 1 ), . . . , F n (x n ); θ) n i=1 f i (x i ), where c(v; θ) = ∂ n C(v; θ) ∂v 1 • • • ∂v n , v = (v 1 , . . . , v n )
, and the derivative is interpreted as a Radon-Nikodym derivative when C(•; θ) is not nth-order differentiable.

To generate X = (X 1 , . . . , X n ) from the joint cdf F X (•), generate V = (V 1 , . . . , V n ) from the copula and return X i = F -1 i (V i ) for each i. Generating V from the copula is not always obvious, but there are classes of copulas for which this can be easily done, one of them being the Archimedean copulas. This important family of copulas can model strong forms of tail dependence using a single parameter, which makes them convenient to use. An Archimedean copula C a is defined by

C a (v; θ) = ψ θ ψ [-1] θ (v 1 ) + . . . + ψ [-1] θ (v n ) ,
where the generator function ψ θ : [0, ∞) → [0, 1] is a strictly decreasing convex function such that lim x→∞ ψ θ (x) = 0, θ ∈ [0, ∞) is a parameter governing the strength of dependence, and ψ

[-1] θ is a pseudo-inverse defined by ψ [-1] θ (x) = 1{0 ≤ x ≤ ψ θ (0)}ψ θ (x)
, with the convention that ψ -1 θ (0) = inf{x : ψ θ (x) = 0}. Archimedean copulas are absolutely continuous, and their densities have the form:

c a (v; θ) = 1 0 ≤ n i=1 ψ -1 θ (v i ) ≤ ψ -1 θ (0) ∂ n ψ θ (x) ∂x n x= n i=1 ψ -1 θ (v i ) n i=1 ∂ψ -1 θ (v i ) ∂v i ,
assuming that the generator function ψ θ (•) is smooth [START_REF] Nelsen | An Introduction to Copulas[END_REF]. In general we do not have an analytical expression for c a (•; θ), and it can be discontinuous in θ. Therefore, the standard LR method typically does not apply to sensitivity analysis with Archimedean copulas. [START_REF] Marshall | Families of multivariate distributions[END_REF] propose the following simple algorithm to generate V from an Archimedean copula with generator function ψ θ (•):

(i) Generate a random variable Y θ from the distribution with Laplace transform ψ θ (•) (with at least one uniform random number as input).

(ii) For i = 1, . . . , n, let

V i = ψ θ (-(log U i )/Y θ ) with i.i.d. U i ∼ U(0, 1).
For a given Y θ , this gives a stochastic model with uniform random numbers U i as inputs:

ϕ(g(U ; θ)) = φ F -1 1 (ψ θ (-(log U 1 )/Y θ )) , . . . , F -1 n (ψ θ (-(log U n )/Y θ )) ,
where ϕ(v 1 , . . . , v n ) = φ(F -1 1 (v 1 ), . . . , F -1 1 (v n )) and g(u; θ) = (ψ θ (-(log u 1 )/Y θ ) , . . . , ψ θ (-(log u n )/Y θ )).

A Generalized Likelihood Ratio Method

In this section, we derive the GLR estimator for the derivative (2) of the expectation of stochastic model (1). The general theory for GLR is first derived, and then it is applied to the two examples in the previous section.

General Theory

Denote the Jacobian of g by

J g (u; θ) :=       ∂g 1 (u;θ) ∂u 1 ∂g 1 (u;θ) ∂u 2 • • • ∂g 1 (u;θ) ∂un ∂g 2 (u;θ) ∂u 1 ∂g 2 (u;θ) ∂u 2 • • • ∂g 2 (u;θ) ∂un . . . . . . . . . . . . ∂gn(u;θ) ∂u 1 ∂gn(u;θ) ∂u 2 • • • ∂gn(u;θ) ∂un      
, and

∂ θ g(u; θ) := ∂g 1 (u; θ) ∂θ , . . . , ∂g n (u; θ) ∂θ T ,
with the superscript T indicating vector transposition. In addition, we define two weight functions in the GLR estimator:

r i (u; θ) := J -1 g (u; θ) ∂ θ g(u; θ) T e i , i = 1, . . . , n, (3)
d(u; θ) := n i=1 e T i J -1 g (u; θ) (∂ u i J g (u; θ))J -1 g (u; θ)∂ θ g(u; θ) -trace(J -1 g (u; θ) ∂ θ J g (u; θ)), (4)
where e i is the ith unit column vector and ∂ z J g is the matrix obtained by differentiating J g with respect to z element-wise. Let x -and x + be limits taken from the left-hand side and right-hand side of x, respectively, and for a function h(•), denote h(x -) := lim x→x -h(x) and h(x

+ ) := lim x→x + h(x).
The following conditions are introduced to guarantee the unbiasedness of the proposed GLR derivative estimator.

(A.1) There exist a smooth function

ϕ (•) ∈ C ∞ and p > 1 such that lim →0 sup θ∈Θ (0,1) n |ϕ (g(u; θ)) -ϕ(g(u; θ))| p du = 0,
and if n > 1, for a fixed ε > 0 and any

u i ∈ (0, 1) \ [ε, 1 -ε], i = 1, . . . , n, lim →0 sup θ∈Θ (0,1) n-1 |ϕ (g(u; θ)) -ϕ(g(u; θ))| p du -i = 0,
where u -i := (u 1 , . . . , u i-1 , u i+1 , . . . , u n ), and if n = 1, for any u ∈ (0, 1)

\ [ε, 1 -ε], lim →0 sup θ∈Θ |ϕ (g(u; θ)) -ϕ(g(u; θ))| = 0. (A.
2) The Jacobian J g (u; θ) is invertible almost everywhere (a.e.), and the performance function g(u; θ) is twice continuously differentiable with respect to (u, θ) ∈ (0, 1) n × Θ, where Θ is a bounded neighborhood of the parameter θ of interest.

(A.

3) The following integrability conditions hold:

(0,1) n-1 sup θ∈Θ,u i ∈(0,1) ϕ(g(u; θ)) r i (u; θ) du -i < ∞, i = 1, . . . ,

n, and

(0,1) n sup θ∈Θ ϕ(g(u; θ)) d(u; θ) du < ∞.

(A.4) The function g(•; θ) is invertible, and lim

u i →1 - sup θ∈Θ,u -i ∈(0,1) n-1 |r i (u; θ)| = lim u i →0 + sup θ∈Θ,u -i ∈(0,1) n-1 |r i (u; θ)| = 0, i = 1, . . . , n.
Remark 1 Condition (A.1) can be checked in certain settings when ϕ (•) can be explicitly constructed; see Proposition 1. The invertibility of the Jacobian matrix in condition (A.2) justifies the local invertibility of function g(•; θ), whereas global invertibility of g(•; θ) in condition (A.4) is stronger, although much weaker than requiring an explicit inverse function for g(•; θ) in deriving the push-out LR estimator [START_REF] Rubinstein | Discrete Event Systems: Sensitivity Analysis and Stochastic Optimization by the Score Function Method[END_REF]. In general, it is difficult to find an explicit inverse function for a nonlinear function g(•; θ), but the existence of the inversion could be guaranteed by the inverse function theorem.

Unbiasedness of the new GLR estimator developed in this work is established under two sets of conditions in the following theorem.

Theorem 1 Under conditions (A.1) -(A.3) or (A.2) -(A.4), ∂E[ϕ(g(U ; θ))] ∂θ = E[G(U ; θ)], where (5)
G(U ; θ) := n i=1 ϕ(g(U i ; θ))r i (U i ; θ) -ϕ(g(U i ; θ))r i (U i ; θ) + ϕ(g(U ; θ))d(U ; θ), with U i := (U 1 , . . . , 1 - ith element , . . . , U n ), U i := (U 1 , . . . , 0 + ith element , .
. . , U n ), and r i (•) and d(•) defined by

(3) and (4), respectively.

Remark 2 The proof of the theorem can be found in the online Appendix A. Even if g(U i ; θ) = ∞ or g(U i ; θ) = ∞, the GLR estimator could still be well defined; see e.g. Section 3.1. The difference between this estimator and the GLR estimator in Peng et al. (2018) is that the surface integration part in Peng et al. (2018) is shown to be zero under certain conditions including that the tails of the input densities are required to go to zero smoothly and fast enough, whereas here the surface integration part is included in the estimator and can be estimated by simulation. In the case where the surface integration part becomes zero, we can prove the result without assuming (A.1), and we also avoid the integrability condition in Peng et al. (2018) on certain intermediate quantities (the smoothed function), which is difficult to verify in practice. The proof is obtained by first truncating the support of the input random variables to a compact set and then appropriately expanding it to the whole space.

We now examine the special case where g(u; θ) = (g 1 (u 1 ; θ), . . . , g n (u n ; θ)), which covers Examples 1 and 2. For n independent uniform random numbers, we can take g i (u i ; θ) = F -1 i (u i ; θ) for i = 1, . . . , n, while for dependent variables governed by an Archimedean copula, conditional on

Y θ = y, we have g i (u i ; θ) = ψ θ (-log u i /y) for i = 1, . . . , n.
In this special case, the Jacobian becomes

J g (u; θ) =       ∂g 1 (u 1 ;θ) ∂u 1 0 • • • 0 0 ∂g 2 (u 2 ;θ) ∂u 2 • • • 0 . . . . . . . . . . . . 0 0 • • • ∂gn(un;θ) ∂un       .
Then we have

r i (u; θ) = ∂g i (u i ; θ) ∂θ ∂g i (u i ; θ) ∂u i , i = 1, . . . , n, and
d(u; θ) = n i=1 ∂g i (u i ; θ) ∂θ ∂ 2 g i (u i ; θ) ∂u 2 i ∂g i (u i ; θ) ∂u i 2 - ∂ 2 g i (u i ; θ) ∂θ∂u i ∂g i (u i ; θ) ∂u i .
Moreover, condition (A.1) in Theorem 1 can be replaced by a set of a simpler assumptions when

the performance function ϕ(x) is a product of n indicators: ϕ(x) = n i=1 1{x i ≤ 0},
ϕ(g(U ; θ)) = n i=1 1{g i (U i ; θ) ≤ 0}. Condition (A.1) holds if for i = 1, . . . , n and a fixed ε > 0 in (6), inf θ∈Θ,u i ∈[ε,1-ε] ∂g i (u i ; θ) ∂u i > 0 and inf θ∈Θ,u i ∈(0,1)\[ε,1-ε] g i (u i ; θ) > 0. ( 6
)
If the functions g i can be decomposed as products of the form

g i (u i ; θ) = ξ i (θ)η i (u i ) for i = 1, . . . , n, then conditions (A.
3), (A.4), and ( 6) can be simplified. For example, an exponential random variable with mean θ can be generated bylog(U i )/θ where U i ∼ U(0, 1). When this decomposition holds, we can write

r i (u; θ) = d log ξ i (θ) dθ d log η i (u i ) du i , i = 1, . . . , n, and
d(u; θ) = n i=1 d log ξ i (θ) dθ η i (u i )η i (u i ) (η i (u i )) 2 -1 .
(A.3') Boundedness and integrability conditions on functions of input random numbers:

inf u i ∈(0,1) d log η i (u i ) du i > 0 and E η i (U i )η i (U i ) (η i (U i )) 2 < ∞, i = 1, . . . , n.
(A.4') Boundary condition on functions of input random numbers:

lim u i →1 - d log η i (u i ) du i = lim u i →0 + d log η i (u i ) du i = ∞, i = 1, . . . , n. Corollary 1 Suppose that g i (u i ; θ) = ξ i (θ)η i (u i ), i = 1, . . . , n, and ϕ(•) is bounded and max i=1,...,n sup θ∈Θ ∂ log ξ i (θ) ∂θ < ∞.
Then conditions (A.3) and (A.4) can be replaced by (A.3') and (A.4'), respectively, and condition (6) in Proposition 1 also simplifies to

inf θ∈Θ ξ i (θ)| > 0, inf u i ∈[ε,1-ε] η i (u i ) > 0, inf u i ∈(0,1)\[ε,1-ε] η i (u i ) > 0, i = 1, . . . , n.
Unbiasedness of the GLR estimators in many examples of this paper can be justified by verifying these simplified conditions.

The Independent Case

Let us return to the independent case of Example 1 and suppose that each X i is continuous with density f i (•; θ). Our goal is to estimate

∂E[ϕ(F -1 1 (U 1 ; θ), . . . , F -1 n (U n ; θ))] ∂θ
, for which the Jacobian is

J g (u; θ) =      1 f 1 (X 1 (u 1 ;θ);θ) 0 • • • 0 0 1 f 2 (X 2 (u 2 ;θ);θ) • • • 0 . . . . . . . . . . . . 0 0 • • • 1 fn(Xn(un;θ);θ)     
, and

∂ θ g(u; θ) = ∂X 1 (u 1 ; θ) ∂θ , . . . , ∂X n (u n ; θ) ∂θ T , where X i (u i ; θ) := F -1 i (u i ; θ) and ∂X i (u i ; θ) ∂θ := - ∂F i (x i ; θ) ∂θ f i (x i ; θ) x i =X i (u i ;θ)
.

Then the weight functions in the GLR estimator are:

r i (u; θ) = - ∂F i (x i ; θ) ∂θ x i =X i (u i ;θ) and d(u; θ) = n i=1 ∂ log f i (x i ; θ) ∂θ x i =X i (u i ;θ)
, so lim

u i →1 - r i (u; θ) = lim u i →0 + r i (u; θ) = 0. Therefore, ∂E[ϕ(F -1 1 (U 1 ; θ), . . . , F -1 n (U n ; θ))] ∂θ = E ϕ(F -1 1 (U 1 ; θ), . . . , F -1 n (U n ; θ)) d(U ; θ) = E ϕ(X) n i=1 ∂ log f i (X i ; θ) ∂θ .
The expression inside the last expectation coincides with the classic LR derivative estimator in the case where the LR method is applicable, i.e., when there are no structural parameters in the sample performance [START_REF] Glynn | Likelihood ratio gradient estimation for stochastic systems[END_REF]. From this perspective, the GLR method generalizes the LR method by allowing the appearance of structural parameters.

Archimedean Copulas

We consider estimating

∂ ∂θ E ϕ ψ θ - log U 1 Y θ , . . . , ψ θ - log U n Y θ ,
where the expectation is with respect to both Y θ and the independent U i , i = 1, . . . , n. By conditioning, we can use a mixture of LR and GLR:

∂ ∂θ E ϕ ψ θ - log U 1 Y θ , . . . , ψ θ - log U n Y θ = E ϕ ψ θ - log U 1 Y θ , . . . , ψ θ - log U n Y θ ∂ log f Y (y; θ) ∂θ y=Y θ + E     ∂E ϕ ψ θ -log U 1 y , . . . , ψ θ -log Un y ∂θ y=Y θ     , (7)
where f Y (•; θ) is the density function of Y θ . This equality follows from Theorem 1 of L' Ecuyer (1990) with ω in that theorem replaced by y and h(θ, ω) replaced by h(y, θ) := E[ϕ(ψ θ (-ln U 1 /y), . . . , ψ θ (-ln U n /y))], under the assumption that for all y, this last expectation is continuous in θ and differentiable except perhaps on a countable set. Specifically, (7) can be rewritten as

∂ ∂θ R h(y, θ)F Y (dy; θ) = ∂ ∂θ R h(y, θ)L(y, θ, θ 0 )F Y (dy; θ 0 ) = R h(y, θ) ∂L(y, θ, θ 0 ) ∂θ + ∂h(y, θ) ∂θ L(y, θ, θ 0 ) F Y (dy; θ 0 ),
where L(y, θ, θ 0 ) := f Y (y; θ)/f Y (y; θ 0 ). The first term on the right-hand side of ( 7) can be dealt with by the LR method straightforwardly if f Y (•; θ) admits an analytical form. [START_REF] Glasserman | Sensitivity estimates from characteristic functions[END_REF] show how to apply the LR method with only the Laplace transform ψ θ (•).

We now show how to use GLR to handle the second term on the right-hand side of ( 7) with Y θ fixed and generated from other uniform random numbers. The Archimedean copula model falls into the special case where g(u; θ) = (g 1 (u 1 ; θ), . . . , g n (u n ; θ)), discussed after Theorem 1. The Jacobian in this case is

J g (u; θ, y) =      -1 u 1 y ψ θ -log u 1 y 0 • • • 0 0 -1 u 2 y ψ θ -log u 2 Y θ • • • 0 . . . . . . . . . . . . 0 0 • • • -1 uny ψ θ -log un y     
, and

∂ θ g(u; θ, y) = ∂ψ θ (x 1 ) ∂θ x 1 =- log u 1 y , . . . , ∂ψ θ (x n ) ∂θ xn=-log un y T .
The weight functions in the GLR estimator are

r i (u; θ, y) = - u i y ψ θ (x i ) ∂ψ θ (x i ) ∂θ x i =- log u i y , d(u; θ, y) = n i=1 - 1 ψ θ (x i ) ∂ψ θ (x i ) ∂θ + ψ θ (x i ) (ψ θ (x i )) 2 ∂ψ θ (x i ) ∂θ + ∂ψ θ (x i ) ∂θ y ψ θ (x i ) x i =-log u i y .
Example 3 The Clayton Copula. The generator function for the Clayton copula is

ψ θ (x) = (1 + x) -1 θ , θ ∈ (0, ∞).
Then,

∂ψ θ (x) ∂θ = 1 θ 2 log(1 + x)(1 + x) -1 θ , ∂ψ θ (x) ∂θ = 1 θ 2 (1 + x) -1 θ -1 1 - 1 θ log(1 + x) , and
ψ θ (x) = - 1 θ (1 + x) -1 θ -1 , ψ θ (x) = 1 θ 1 θ + 1 (1 + x) -1 θ -2 .
By the inverse Laplace transformation, we find that Y θ ∼ Γ(1/θ, 1), the gamma distribution with

density f Y (y; θ) = y 1 θ -1 e -y
Γ(1/θ) , where Γ(s) := ∞ 0 t s-1 e -t dt, and the LR term is

∂ log f Y (y; θ) ∂θ = - d log Γ(1/θ) dθ - 1 θ 2 log y.
The weight functions in the GLR estimator are

r i (u; θ, y) = - 1 θ u i y 1 - log u i y log 1 - log u i y , d(u; θ, y) = n i=1 1 θ [1 + (1 -(x i + 1)y) log(1 + x i )] x i =- log u i y .
In addition, we have lim

u i →0 + r i (u; θ, y) = lim u i →1 -r i (u; θ, y) = 0.
GLR for Ali-Mikhail-Haq copulas can be found in the online Appendix A. For both the Clayton and Ali-Mikhail-Haq copulas, conditions (A.2) and (A.4) in Theorem 1 are satisfied. If ϕ(•) is bounded, condition (A.3) in Theorem 1 can also be verified straightforwardly for any y = Y θ .

Applications

We apply the GLR method to distribution sensitivity estimation, and estimate sensitivities for stopping time problems and credit risk derivatives, with specific forms for the function ϕ(•).

Distribution Sensitivities

For g(•; θ) : (0, 1) n → R, we estimate the following two first-order distribution sensitivities:

∂F (z; θ) ∂θ = ∂E[1{g(U ; θ) -z ≤ 0}] ∂θ = E ∂E[1{g(U i , U -i ; θ) -z ≤ 0}|U -i ] ∂θ , f (z; θ) = ∂E[1{g(U ; θ) -z ≤ 0}] ∂z = E ∂E[1{g(U i , U -i ; θ) -z ≤ 0}|U -i ] ∂z , where f (•; θ) is the density function of Z(θ) = g(U ; θ) and U -i := (U 1 , . . . , U i-1 , U i+1 , . . . , U n ), i = 1, . . . , n.
By applying GLR, we obtain

E ∂E[1{g(U i , U -i ; θ) -z ≤ 0}|U -i ] ∂θ = E[G 1,i (U ; z, θ)],
where

G 1,i (U ; z, θ) := 1{g(U i ; θ) -z ≤ 0}r i (U i ; θ) -1{g(U i ; θ) -z ≤ 0}r i (U i ; θ) + 1{g(U ; θ) -z ≤ 0}d(U ; θ), r i (u; θ) = ∂g(u; θ) ∂u i -1 ∂g(u; θ) ∂θ , and
d(u; θ) = ∂g(u; θ) ∂u i -1 ∂g(u; θ) ∂u i -1 ∂g(u; θ) ∂θ ∂ 2 g(u; θ) ∂u 2 i - ∂ 2 g(u; θ) ∂u i ∂θ .
We also obtain

E ∂E[1{g(U i , U -i ; θ) -z ≤ 0}|U -i ] ∂z = E[G 2,i (U ; z, θ)],
where

G 2,i (U ; z, θ) :=1{g(U i ; θ) -z ≤ 0}r i (U i ; θ) -1{g(U i ; θ) -z ≤ 0}r i (U i ; θ) + 1{g(U ; θ) -z ≤ 0} d(U ; θ), with ri (u; θ) = - ∂g(u; θ) ∂u i -1 and d(u; θ) = - ∂g(u; θ) ∂u i -2 ∂ 2 g(u; θ) ∂u 2 i .
To establish the unbiasedness of G 1,i (U ; z, θ) and G 2,i (U ; z, θ) conditional on U -i , condition (A.1)

in Theorem 1 can be justified by checking the conditions in Proposition 1. Given this conditional unbiasedness, the unconditional unbiasedness of G 1,i (U ; z, θ) and G 2,i (U ; z, θ) will follow from

E sup θ∈Θ |E[G 1,i (U ; z, θ)|U -i ] < ∞ and E sup z∈Z |E[G 2,i (U ; z, θ)|U -i ]| < ∞,
where Z is a neighborhood of z. Since the indicator function is bounded, the conditions of Proposition 1 follow from the integrability condition on the weight functions: for i = 1, . . . , n,

E sup θ∈Θ |r i (U i ; θ)] < ∞, E sup θ∈Θ |r i (U i ; θ)] < ∞, E sup θ∈Θ |d(U ; θ)] < ∞, and
E |r i (U i ; θ)] < ∞, E [|r i (U i ; θ)]] < ∞, E | d(U ; θ)] < ∞.
The GLR estimator for estimating the distribution sensitivities is not unique. We can consider the above GLR estimator for each i and construct the following linear combination of these n GLR estimators with real-values weights w i , as in Hammersley and Handscomb (1964, p.19):

n i=1 w i G r,i (U ; z, θ) subject to n i=1 w i = 1, r = 1, 2.
An optimal GLR estimator, which minimizes the variance, can be obtained by solving arg min (w 1 ,...,wn)

Var

n i=1 w i G r,i (U ; z, θ) subject to n i=1 w i = 1.
This leads to the optimal weights

w * i = e T i Σ -1 e e T Σ -1 e , i = 1, . . . , n, (GLR-Opt)
where e = (1, . . . , 1) T , e i is a d-dimensional unit vector in ith direction, and Σ = (Σ i i ) n×n is the covariance matrix of (G r,1 (U ; z, θ), . . . , G r,n (U ; z, θ)). In practice, w * i 's must be estimated, and such estimators wiil be correlated with the G r,i . This linear combination idea is equivalent to a control variate formulation.

Example 4 Distribution Sensitivities for Quantiles. For 0 ≤ α ≤ 1, the α-VaR (or α-quantile) of a random variable Z(θ) = 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 Z j := g(U (j) ; θ), j = 1, . . . , m. The empirical α-quantile F -1 m (α), which is the inverse of the empirical distribution evaluated at α, is simply

Z ( αm ) , where Z (1) < • • • < Z (m)
are the realizations of Z 1 , . . . , Z m sorted in increasing order (the order statistics), and x denotes the smallest integer greater than or equal to x. 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 α (θ); θ) .
Traditionally, batching and sectioning techniques are used to estimate the asymptotic variance to construct a confidence interval on the empirical quantile, and these methods lead to subcanonical convergence rates [START_REF] Nakayama | Confidence intervals for quantiles using sectioning when applying variance-reduction techniques[END_REF]. With the GLR density estimator, however, we can estimate the asymptotic variance by

mα(1 -α) m j=1 G 2,i (U (j) ; z, θ)| z= F -1 m (α)
, using the same realizations of the uniform random variables U (j) as in the quantile estimator F -1 m (α). It follows from [START_REF] Peng | On the asymptotic analysis of quantile sensitivity estimation by Monte Carlo simulation[END_REF] that this asymptotic variance estimator is consistent. We can also estimate the quantile sensitivity q α (θ) = ∂q α (θ)/∂θ by estimating distribution sensitivities. By the implicit function theorem,

q α (θ) = - ∂F (z; θ) ∂θ z=qα(θ)
f (q α (θ); θ), so the quantile sensitivity can be estimated by the ratio

D(θ) := - m j=1 G 1,i (U (j) ; z, θ) m j=1 G 2,i (U (j) ; z, θ) z= F -1 m (α)
.

This ratio estimator is biased, because the expectation of the ratio is not equal to the ratio of expectations and also because of the bias in the quantile estimator, but from [START_REF] Peng | On the asymptotic analysis of quantile sensitivity estimation by Monte Carlo simulation[END_REF], it is consistent and obeys a central limit theorem when m → ∞. Its asymptotic variance can also be estimated by the GLR estimator for higher-order distribution sensitivities (Peng et al., 2020a).

Stopping Time Problems

In this subsection, we consider estimating the derivative of the expectation of a sample performance that depends on a stopping time N :

∂E [ϕ N (X 1 , . . . , X N )] ∂θ , (8)
where {X i , i ≥ 1} is a Markov chain defined by the following stochastic recurrence:

X 1 = g 1 (U 1 ; θ), X i = κ(X i-1 , g i (U i ; θ)), (9)
and N is the first time the Markov chain hits a set Ω, i.e.,

N = min {n ∈ N : X n ∈ Ω} = min {n ∈ N : κ(X n-1 , g n (U n ; θ)) ∈ Ω} .
The expectation of this sample performance can be rewritten as

E [ϕ N (X 1 , . . . , X N )] = ∞ n=1 E ϕ n (X 1 , . . . , X n ) n-1 i=1 1{X i / ∈ Ω}1{X n ∈ Ω} . ( 10
)
For each expectation in the summation on the right-hand side of the equation, the stochastic model falls into the special case g(u; θ) = (g 1 (u 1 ; θ), . . . , g n (u 1 ; θ)) discussed after Theorem 1, and the derivative of can be estimated by the GLR method. We have

∂E [ϕ N (X 1 , . . . , X N )] ∂θ = ∞ n=1 n i=1 E ϕ n (X +i 1 , . . . , X +i n ) n-1 j=1 1{X +i j / ∈ Ω}1{X +i n ∈ Ω}r i (1 -; θ) - ∞ n=1 n i=1 E ϕ n (X -i 1 , . . . , X -i n ) n-1 j=1 1{X -i j / ∈ Ω}1{X -i n ∈ Ω}r i (0 + ; θ) + ∞ n=1 E ϕ n (X 1 , . . . , X n ) n-1 i=1 1{X i / ∈ Ω}1{X n ∈ Ω}d n (U 1 , . . . , U n ; θ) ,
where {X +i j , j ≥ 1} and {X -i j , j ≥ 1} are Markov chains generated by ( 9) with U i replaced by 1 and 0, respectively, and d n (•; θ) denotes the second part in the GLR estimator for the nth term on the right-hand side of (10). We use the randomization technique of [START_REF] Rhee | Unbiased estimation with square root convergence for SDE models[END_REF] to obtain a single-run unbiased estimator of ( 8), which is given by

1 p(N ) S N (U 1 , . . . , U N ) + D N (U 1 , . . . , U N ), ( 11
)
where N is a discrete random variable supported on Z + with probability mass function p(•),

S n (U 1 , . . . , U n ) := n i=1 ϕ n (X +i 1 , . . . , X +i n ) n-1 j=1 1{X +i j / ∈ Ω}1{X +i n ∈ Ω}r i (1 -; θ) - n i=1 ϕ n (X -i 1 , . . . , X -i n ) n-1 j=1 1{X -i j / ∈ Ω}1{X -i n ∈ Ω}r i (0 + ; θ), and
D n (U 1 , . . . , U n ) := ϕ n (X 1 , . . . , X n ) n-1 i=1 1{X i / ∈ Ω}1{X n ∈ Ω}d n (U 1 , . . . , U n ; θ).
This stopping time problem generalizes those in [START_REF] Glynn | Likelihood ratio gradient estimation: an overview[END_REF] and [START_REF] Heidergott | Gradient estimation for a class of systems with bulk services (a problem in public transportation)[END_REF] by allowing the distribution of the stopping time N to depend on parameter θ. The classical IPA and LR methods do not cover this case. In the next example, we discuss sensitivity analysis for control charts, widely used in statistical quality control.

Example 5 Controls Charts. Control charts aim to detect (statistically) when a manufacturing or business process goes out of control. The system is assumed to output samples having different statistical distributions when in control versus when out of control. The systems goes out of control at the (unobservable) random time χ with cdf F 0 (•). The output Z i of the ith sample has cdf F 1 (•) when in control and F 2 (•) when out of control. For an exponentially weighted moving average (EWMA) chart, the (observable) test statistic after the ith sample is

Y i = αZ i + (1 -α)Y i-1 .
The stopping time N is the time when the system is declared out of control:

N = min{i : Y i / ∈ [θ 1 , θ 2 ]} = min{i : X i / ∈ [0, 1]},
where θ 1 and θ 2 are fixed lower and upper control limits for the test statistic, and

X i = Y i -θ 1 θ 2 -θ 1 . Note that {X i , i ≥ 1} is a Markov chain that follows the recursion X i = (1 -α)X i-1 + g i (U i ; θ), where
g i (U i ; θ) = α θ 2 -θ 1 1{i < χ/∆}F -1 1 (U i ) + 1{i ≥ χ/∆}F -1 2 (U i ) -θ 1 ,
and ∆ is the sampling period (the time between any two successive monitoring epochs).

Here we consider sensitivity analysis with respect to θ = θ 2 . This model falls into the special case of Corollary 1 where

g i (u i ; θ) = ξ i (θ)η i (u i ), with ξ i (θ) = α/(θ -θ 1 ) and η i (u i ) = 1{i < χ/∆}F -1 1 (u i ) + 1{i ≥ χ/∆}F -1 2 (u i ) -θ 1 , for i = 1, . . . , n. Then we have r i (u; θ) = - 1 θ 2 -θ 1 1{i < χ/∆}F -1 1 (u i )f 1 (F -1 1 (u i )) + 1{i ≥ χ/∆}F -1 2 (u i )f 2 (F -1 2 (u i )) -θ 1 , and
d(u; θ) = 1 θ 2 -θ 1 n i=1 1{i < χ/∆} F -1 1 (u i ) f 1 (F -1 1 (u i )) + 1{i ≥ χ/∆} F -1 2 (u i ) f 2 (F -1 2 (u i )) + 1 .
Suppose we pay a cost of c per unit of time when the system is out of control and this is not yet detected, and a one-time cost C to fix the system when it is declared out of control. This model is a regenerative process, which regenerates each time we fix the system. The expected cost over one regenerative cycle is cE[(Nχ/∆ ) + ] + C, and therefore the average cost per unit of time over

an infinite horizon is cE[(N -χ/∆ ) + ] + C E[N ] . ( 12
)
The goal might be to select the control limits θ 1 and θ 2 to minimize this average cost. Widening the gap θ 2 -θ 1 would reduce the frequency of intervention, so we would pay the fixed cost C less often, but then the penalty c would be paid over longer periods of time on average. The optimal control limits achieve an optimal balance between these two types of costs. The sample performance of the expectation in the numerator of ( 12) is different from those treated by SPA in [START_REF] Fu | Efficient design and sensitivity analysis of control charts using Monte Carlo simulation[END_REF] and Fu et al. (2009b), the development of which depends on the specific structure of the problem.

The GLR method in our work provides unbiased derivative estimators for the expectations in both the numerator and denominator of (12).

Credit Risk Derivatives

We consider two important types of credit risk derivatives: basket default swaps (BDSs) and collateralized debt obligations (CDOs). In a BDS contract, the buyer pays fixed premia p 1 , ..., p k to the protection seller at dates 0 < T 1 < . . . < T k < T , and if the ith default time τ (i) occurs before T , i.e., τ (i) < T , these premium payments stop, and the seller undertakes the loss of the ith default and makes a payment to the buyer. Let L i be the loss of the ith default. The discounted value of the ith default swap is the difference between the discounted payments made by the seller and those made by the buyer:

V bds (τ ) = V value (τ ) -V prot (τ ),
where V prot (τ ) is the discounted premium payed by the buyer:

V prot (τ ) = j=1 p j exp(-rT j ) + p +1 exp(-rτ (i) ) τ (i) -T T +1 -T , if T ≤ τ (i) ≤ T +1 , k j=1 p i exp(-rT j ), if τ (i) > T ,
and V value (X) the discounted payment made by the seller:

V value (X) = L (i) exp(-rτ (i) )1{τ (i) < T }.
In a CDO, the losses caused by the defaults of the assets in the portfolio are packaged together and then tranched. The tranches are ordered so that losses are absorbed sequentially following the order of the tranches. For example, a tranche of a CDO absorbs the loss above a threshold L -and below a threshold L + , i.e.,

V cdo (τ ) = (L -L -) • 1{L > L -} -(L -L + ) • 1{L > L + }, where L = n i=1 L i • 1{τ i < T }.
Suppose that the vector of default times (τ 1 , . . . , τ n ) have a joint distribution with marginal cdf's

F i (•), i = 1, .
. . , n, and a dependence structure modeled by an Archimedean copula. Then we can generate τ i 's by generating V = (V 1 , . . . , V n ) from the copula, and putting

τ i = F -1 i (V i ), i = 1, . . . , n.
The sample performances V value (τ ), V prot (τ ), and V cdo (τ ) may be discontinuous with respect to the structural parameter θ in the copula model due to the presence of indicator functions and order indices. As a result, neither IPA nor LR can be applied directly for this model. On the other hand, V value (τ ), V prot (τ ), and V cdo (τ ) are all of the form ϕ(g(U ; θ)) that fits our framework, due to the generality of the measurable function ϕ(•). Unlike in [START_REF] Lei | Sensitivity analysis of portfolio credit derivatives by conditional Monte Carlo simulation[END_REF] where a separate CMC technique needs to be derived for each type of cash flow, the GLR method in this work can estimate the derivative of the expectation for all three types of cash flows.

Numerical Experiments

In this section, we present numerical examples to demonstrate broad applicability and flexibility of the proposed GLR method to estimate sensitivities in various situations. The examples include an indicator function applied to a linear combination of two exponential random variables, a stochastic activity network (SAN), control charts, and a CDO model.

Distribution Sensitivities for a Linear Model

We estimate distribution sensitivities where ϕ(•) is an indicator function for a linear combination of two independent exponential random variables with means 1/λ 1 and 1/λ 2 , i.e.,

ϕ(g(U ; θ)) = 1{g(U ; θ) ≤ z} where g(U ; θ) = - θ λ 1 log(U 1 ) - 1 λ 2 log(U 2 ).
This sample performance falls into the special case stated in Corollary 1, and the conditions in Proposition 1 and the integrability condition on the weight function discussed in Section 4.1 can be checked straightforwardly. We have Then the GLR estimators for distribution sensitivities in Section 4.1 are

∂g(u; θ) ∂u 1 = - θ λ 1 u 1 , ∂g(u; θ) ∂u 2 = - 1 λ 2 u 2 , ∂g(u; θ) ∂θ = - 1 λ 1 log u 1 , and
∂ 2 g(u; θ) ∂u 2 1 = θ λ 1 u 2 1 , ∂ 2 g(u; θ) ∂u 2 2 = 1 λ 2 u 2 2 , ∂ 2 g(u; θ) ∂θ∂u 1 = - 1 λ 1 u 1 , ∂ 2 g(u; θ) ∂θ∂u 2 = 0.
G 1,1 (U ; z, θ) = - 1 θ 1{g(U ; θ) ≤ z} (log U 1 + 1) , G 1,2 (U ; z, θ) = λ 2 λ 1 1 - θ λ 1 log U 1 ≤ z log U 1 -1{g(U ; θ) ≤ z} log U 1 , G 2,1 (U ; z, θ) = λ 1 θ 1 - 1 λ 2 log U 2 ≤ z -1{g(U ; θ) ≤ z} , G 2,2 (U ; z, θ) = λ 2 1 - θ λ 1 log U 1 ≤ z -1{g(U ; θ) ≤ z} .
For our numerical experiments, we take λ 1 = 1, λ 2 = 1, and θ = 1, and we estimate the density of g(U ; θ) at z, as a function of z from z = 0.01 to 10 with a step size of 0.01. We perform 10 6 independent simulation runs with GLR. The curves of the estimated density and estimated variances as a function of z are given in Figure 1. We also estimate the sensitivity of E[ϕ(g(U ; θ))]

with respect to θ at z, as a function of z from z = 0.01 to 10 with a step size of 0.01. The numerical observations are similar, and the details can be found in the online Appendix B. with weights given by (GLR-Opt). Due to page limit, the variance comparison between GLR and the finite difference method with common random numbers (FDC) is relegated to online Appendix B. The variance of FDC is much larger than those of the GLR estimators, and the variance of FDC(0.01) is about 10 times of that of FDC(0.1), which indicates that FDC suffers from a biasvariance tradeoff issue. (b) Coverage rates of confidence intervals.

Figure 2 90%-confidence intervals estimated by GLR for quantiles and coverage rates in example of Section 5.1.

In Figure 2(a), the solid line in the center is the estimate for the α-quantiles q α of g(U ; θ), α = 0.1, . . . , 0.9, using a batch of 10 4 independent simulation replications, and the upper and lower dotted lines are calculated respectively by adding and subtracting twice the estimated standard deviation to the estimated quantile values, which are estimated by GLR using the same batch of 10 4 independent simulation replications for estimating the quantiles. Figure 2(b) presents the coverage rates of the 90%-confidence intervals for quantiles q α , α = 0.1, . . . , 0.9, by 10 4 independent macro experiments. The true quantile values can be calculated by inverting a hypoexponential distribution (L 'Ecuyer et al., 2019, Section 4.6). We can see the coverage rates of the 90%-confidence intervals match the target value statistically. As the quantile gets closer to the tail of the distribution, the variances of the quantile estimates and confidence interval estimates are larger because there are many fewer samples in the tail.

Distribution Sensitivities for a Stochastic Activity Network

We estimate distribution sensitivities for the sample performance of a small SAN depicted in Figure 3. There are five nodes representing different stages of activity. The nodes are connected by the arcs representing the activities in each stage. The durations of activities follow independent exponential distributions, i.e., X i = -1 λ i log(U i ), i = 1, . . . , 6. Let θ = λ 6 . There are three different paths representing the tasks to reach the final stage of a project, i.e., π 1 = (1, 4, 6), π 2 = (2, 5, 6), π 3 = (1, 3, 5, 6), and the completion time for each path is additive, i.e, j∈π i X j , i = 1, 2, 3. The completion time for the entire project is max(X 1 + X 4 + X 6 , X 2 + X 5 + X 6 , X 1 + X 3 + X 5 + X 6 ), and the sample performance for the distribution function of completion time is where ϕ(v 1 , v 2 , v 3 ) = 3 i=1 1{v i ≤ 0}, and we have

ϕ(g(U ; θ)) =1 {max(X 1 + X 4 + X 6 , X 2 + X 5 + X 6 , X 1 + X 3 + X 5 + X 6 ) ≤ z} =1 {X 1 + X 4 + X 6 ≤ z} 1 {X 2 + X 5 + X 6 ≤ z} 1 {X 1 + X 3 + X 5 + X 6 ≤ z} ,
g i (u; θ) = - j∈π i 1 λ j log u j -z, i = 1, 2, 3, ∂ z g(u; θ) = -1, 1, 1 T and ∂ θ g(u; θ) = log u 6 θ 2 1, 1, 1 T .
This sample performance goes beyond the setting in Section 4.1, but it can be put under the more general stochastic model (1). In the theory of GLR, the dimension of the vector of the input uniform random numbers is assumed to be the same as that of the argument vector of function g. For the SAN model, there are six input uniform random numbers, while the dimension of the argument vector of g is three. Therefore, we can arbitrarily choose three uniform random numbers to condition on and treat the remaining three as the inputs to the stochastic models for deriving GLR. We condition on (U 4 , U 5 , U 6 ) and treat (U 1 , U 2 , U 3 ) as the input uniform random numbers for deriving GLR. The Jacobian matrix is

J g (u; θ, z) = -   1 λ 1 u 1 0 0 0 1 λ 2 u 2 0 1 λ 1 u 1 0 1 λ 3 u 3   .
The GLR estimator for E[ϕ(g(U ; θ))] with respect to θ can be obtained by multiplying the GLR estimator for the density bylog U 6 /θ 2 . Here we only report on the performance of the GLR estimators for the density of g(U ; θ). Additional results for quantiles can be found in online Appendix B. We have

r 1 (u; θ, z) = λ 1 u 1 , r 2 (u; θ, z) = λ 2 u 2 , r 3 (u; θ, z) = 0, and d(u; θ, z) = -λ 1 -λ 2 .
In this case, condition (A.1) in Theorem 1 can be justified by checking conditions in Proposition 1. Conditions (A.2) and (A.3) can also be checked straightforwardly.

In the experiments, we set λ i = 1, i = 1, . . . , 6, in the stochastic model, and estimate distribution sensitivity curves from z = 0.01 to 15 with a step size of 0.01. Figure 4 presents the density and variance curves estimated by the GLR method using 10 6 independent simulation replications.

Sensitivities of Control Charts

We estimate sensitivities of the expectations in the numerator and denominator of ( 12) for control charts with respect to upper control limit θ = θ 2 discussed in Section 4.2. When the system is in control, the output sample is assumed to follow a uniform distribution on [-1, 1] and we define

Z i = 2U i -1.
When the system is out of control, the output sample is assumed to follow a uniform distribution on [0, 2] and we define Z i = 2U i . Then the transition function of the Markov chain is

X i = (1 -α)X i-1 + α θ 2 -θ 1 [1{i < χ/∆}(2U i -1) + 1{i ≥ χ/∆}2U i -θ 1 ] ,
and the weight functions in the GLR estimator are

r i (u; θ) = - 1 θ 2 -θ 1 1{i < χ/∆} u i - 1 2 + 1{i ≥ χ/∆}u i - θ 1 2 , d(u; θ) = n θ 2 -θ 1 .
In the experiment, we set α = 1/2, θ 1 = -1, θ 2 = 1, ∆ = 1, X 0 = 0, and assume χ = 1 + 3 log U .

The randomized horizon N in (11) follows a geometric distribution with parameter 0.1. In the random horizon problem, it is not easy to synchronize the two sample paths for the finite difference (FD) method, because perturbing parameter θ affects the stopping time N , so it would require substantially more computational overhead to generate the sample paths using common random numbers to implement FDC. In general, decreasing the perturbation size δ would reduce the bias of FD(δ) but increase the variance. In Table 1, we can see that the sensitivity results estimated by GLR match those estimated by FD(0.01), while GLR has smaller variance than FD(0.01), whereas the sensitivity results estimated by FD(0.1) are significantly biased.

Sensitivities of Collateralized Debt Obligations

We estimate the sensitivity with respect to the parameter θ that governs the dependence in the copula model for the expectation of the loss absorbed by the tranche that covers the first 30% of the total losses for 10 assets if there are defaults, i.e., L -= 0 and L + = 0.3 × ( 10 i=1 L i ). We set r = 0.1 and T = 1. The marginal distributions of the defaults are assumed to be exponential, so

τ i = -1 λ i log(X i ), i = 1, .
. . , 10. The parameters λ i and loss L i , i = 1, . . . , n, are randomly generated from the uniform distribution over (0, 1) in the experiments.

We compare the GLR estimator with FDC(δ), where δ is the perturbation size. Due to the simplicity of the weight function of GLR, the computational time of GLR barely increases relative to that required to run the simulation model itself, so the sensitivity estimate by GLR is almost a free byproduct that can be obtained simultaneously during the simulation. Table 2 shows the sensitivity estimates with sample sizes m = 10 4 , m = 10 5 , and m = 10 6 under the Clayton copula, for θ = 0.5. We can see that the variances with GLR are comparable to those with FDC(0.1), but smaller than those with FDC(0.01). For sample size m = 10 6 , the estimate with FDC(0.1) lies outside of the 90% confidence interval of the GLR estimate, whereas the estimate with GLR lies in the 90% confidence interval of the FDC(0.01) estimate with the sample size m = 10 7 , which is -0.179 ± 0.002. This indicates that FDC suffers from the bias-variance tradeoff, while GLR is accurate under a relatively small sample size. Numerical results for sensitivities of CDOs under the Ali-Mikhail-Haq copula and sensitivities of BDS under both the Clayton and Ali-Mikhail-Haq copulas can be found in online Appendix B. The observations are similar to those in this example.

Conclusions

In this paper, a GLR method is proposed for a family of stochastic models with uniform random numbers as inputs. The framework studied in this work covers a large range of discontinuities, and it includes many applications such as density estimation, and sensitivity analysis for statistical quality control and credit risk financial derivatives. Since uniform random numbers are the basic building blocks for generating other random variables, our new method significantly relaxes the limitations on the input random variables in Peng et al. (2018) and Peng et al. (2020a). The technical conditions for justifying unbiasedness of GLR are relatively easy to satisfy in practice, and we show how to verify them on illustrative examples. Ongoing work includes combining GLR with randomized quasi-Monte Carlo methods and CMC.

The interchange of the derivative and integration can be justified under condition (A.3), which implies that ∂ θ g(u; θ) is bounded in Ω ε × Θ. By the Gauss-Green Theorem,

[ε,1-ε] n ∇ u ϕ (g(u; θ))J -1 g (u; θ)∂ θ g(u; θ) du 1 • • • du n = n i=1 [ε,1-ε] n-1 ϕ (g(u; θ)) J -1 g (u; θ)∂ θ g(u; θ) T e i j =i du j 1-ε u i =ε - [ε,1-ε] n ϕ (g(u; θ))div J -1 g (u; θ)∂ θ g(u; θ) du 1 • • • du n ,
where for n = 1, the integration in first term on the right-hand side of the equation is absent, and

for h(u) = (h 1 (u), . . . , h n (u)), div(h(u)) := n i=1 ∂h(u) ∂u i . Then div J -1 g (u; θ)∂ θ g(u; θ) = n i=1 ∂ ∂u i e T i J -1 g (u)∂ θ g(u; θ) = n i=1 e T i ∂ u i J -1 g (u; θ)∂ θ g(u; θ) + trace(J -1 g (u)∂ θ J g (u; θ)).
By differentiating equation J -1 g (u; θ)J g (u; θ) = I with respect to u i on both sides, we have 0 =∂ u i J -1 g (u; θ)J g (u; θ)

=∂ u i J -1 g (u; θ)J g (u; θ) + J -1 g (u; θ)∂ u i J g (u; θ),
which leads to

∂ u i J -1 g (u; θ) = -J -1 g (u; θ)∂ u i J g (u; θ)J -1 g (u; θ).
Therefore, we have

d(u; θ) = -div J -1 g (u; θ)∂ θ g(u; θ) .
With the discussion above and condition (A.2),

∂ ∂θ [ε,1-ε] n ϕ (g(u; θ)) du = n i=1 [ε,1-ε] n-1 ϕ (g(u; θ))r i (u; θ) j =i du j 1-ε u i =ε + [ε,1-ε] n ϕ (g(u; θ))d(u; θ) du. (1)
With condition (A.2), r i (u; θ) and d(u; θ) are bounded in Ω ε × Θ, and with condition (A.1), lim

→0 sup θ∈Θ [ε,1-ε] n (ϕ (g(u; θ)) -ϕ(g(u; θ)))d(u; θ) du ≤ sup θ∈Θ [ε,1-ε] n |ϕ (g(u; θ)) -ϕ(g(u; θ))| p du 1/p [ε,1-ε] n |d(u; θ)| q du 1/q ≤C lim →0 [ε,1-ε] n |ϕ (g(u; θ)) -ϕ(g(u; θ))| p du 1/p = 0,
where the first inequality holds by applying Hölder's inequality for 1/p + 1/q = 1, and C is a positive constant. Similarly, we can show the uniform convergence of the first term on the right-hand side of (1) for n > 1, and lim

→0 [ε,1-ε] n ϕ (g(u; θ)) du - [ε,1-ε] n ϕ(g(u; θ)) du = 0.
For n = 1, the first term on the right-hand side of (1) becomes

ϕ (g(1 -ε; θ))r 1 (1 -ε; θ) -ϕ (g(ε; θ))r 1 (ε; θ),
; θ)) -ϕ(g(1 -ε; θ)))r 1 (1 -ε; θ) -(ϕ (g(ε; θ)) -ϕ(g(ε; θ)))r 1 (ε; θ)| = 0. Therefore, ∂ ∂θ [ε,1-ε] n ϕ(g(u; θ)) du = lim →0 ∂ ∂θ [ε,1-ε] n ϕ (g(u; θ)) du = n i=1 [ε,1-ε] n-1 ϕ(g(u; θ)) r i (u; θ) j =i du j 1-ε u i =ε + [ε,1-ε] n ϕ(g(u; θ)) d(u; θ) du.
With condition (A.3), uniform convergence of the integrals on the right-hand side of equation above can be established over θ ∈ Θ as ε → 0, and we have

∂ ∂θ (0,1) n ϕ(g(u; θ)) du = lim ε→0 ∂ ∂θ [ε,1-ε] n ϕ(g(u; θ)) du = n i=1 (0,1) n-1 lim u i →1 - ϕ(g(u; θ)) r i (u; θ) -lim u i →0 + ϕ(g(u; θ)) r i (u; θ) j =i
[ε,1-ε] n-1 ϕ ,L (g(u; θ))r i (u; θ) j =i du j 1-ε u i =ε = 0.
With condition (A.2), g(u; θ) and det(J g (u; θ)) are bounded in Ω ε × Θ, By change of variables and Hölder's inequality, lim

→0 sup θ∈Θ [ε,1-ε] n (ϕ ,L (g(u; θ)) -ϕ L (g(u; θ)))d(u; θ) du = lim →0 sup θ∈Θ S ε (ϕ ,L (x) -ϕ L (x))| det(J g (u; θ))| d(u; θ)| u=g -1 (x;θ) dx ≤ lim →0 ϕ ,L -ϕ L p sup θ∈Θ S ε | det(J g (u; θ)) d(u; θ)| q u=g -1 (x;θ) dx 1/q = 0,
where

S ε := {x ∈ R n : x = g(u; θ), u ∈ [ε, 1 -ε] n }.
With condition (A.3), we have

lim ε→0 lim →0 sup θ∈Θ [ε,1-ε] n ϕ ,L (g(u; θ))d(u; θ) du - (0,1) n ϕ L (g(u; θ))d(u; θ) du = 0. Therefore, ∂ ∂θ (0,1) n ϕ L (g(u; θ)) du = lim ε→0 lim →0 ∂ ∂θ [ε,1-ε] n ϕ L (g(u; θ)) du = (0,1) n ϕ L (g(u; θ)) d(u; θ) du.
With condition (A.3), sup θ∈Θ (0,1) n (ϕ L (g(u; θ)) -ϕ(g(u; θ))) d(u; θ) du = 0.

Then we have

∂ ∂θ (0,1) n ϕ L (g(u; θ)) du = lim L→∞ ∂ ∂θ (0,1) n ϕ L (g(u; θ)) du = (0,1) n ϕ(g(u; θ)) d(u; θ) du,
which proves the theorem.

Proof of Proposition 1

Proof. Define

χ (x) :=                1 x < -, 1 - (x + ) 2 -≤ x ≤ , 0 x > , and
χ (x) := 1 R φ (x -y)φ(y/ )dz,
where φ(•) is the density of the standard normal distribution. By construction, χ (•) is smooth.

From the condition of the proposition, g i (u i ; θ) is strictly monotone with respect to u i on [ε, 1 -ε].

Without loss of generality, we assume g i (u i ; θ) is strictly increasing with respect to u i on [ε, 1 -ε].

Then we have

1 0 χ (g i (u i ; θ)) -1{g i (u i ; θ) ≤ 0} p du i ≤ 1 1 0 R 1{min(y -, 0) ≤ g i (u i ; θ) ≤ max(y + , 0)}φ(y/ )dy p du i ≤ 1 1 0 R 1{min(y -, 0) ≤ g i (u i ; θ) ≤ max(y + , 0)} p φ(y/ )dy du i ≤ 1 R\[- √ , √ ] φ(y/ )dy + 1 0 1{-≤ g i (u i ; θ) ≤ + √ } du i = 1 R\[- √ , √ ] φ(y/ )dy + g -1 i ( + √ ; θ) -g -1 i (-; θ),
where the second inequality holds by applying Jensen's inequality. We have

g -1 i ( + √ ; θ) -g -1 i (-; θ) = (2 + √ ) ∂g i (u i ; θ) ∂u i -1 u i =g -1 i (ξ ;θ)
,

where ξ ∈ (-, + √ ). Since inf θ∈Θ,u i ∈(0,1)\[ε,1-ε] |g i (u i ; θ)| > 0, g -1 (ξ ; θ) ∈ [ε, 1 -ε] for a sufficiently small . Therefore, lim →0 sup θ∈Θ 1 0 ϕ (g(u; θ) -z) -1{g(u; θ) ≤ z} p du i ≤ lim →0 sup θ∈Θ (2 + √ ) ∂g i (u i ; θ) ∂u i -1 u i =g -1 (ξ ;θ) ≤ lim →0 (2 + √ ) inf θ∈Θ,u i ∈[ε,1-ε] ∂g i (u i ; θ) ∂u i -1 = 0. From inf θ∈Θ,u i ∈(0,1)\[ε,1-ε] |g i (u i ; θ)| > 0, for u i ∈ (0, 1) \ [ε, 1 -ε], lim →0 sup θ∈Θ χ (g i (u i ; θ)) -1{g i (u i ; θ) ≤ 0} = lim →0 1{|g i (u i ; θ)| ≤ } = 0.
In addition,

1 0 • • • 1 0 n j=1 χ (g j (u j ; θ)) - n j=1 1{g j (u j ; θ) ≤ 0} p du 1 • • • du n = (0,1) n n i=1 i j=1 χ (g j (u j ; θ)) n j=i+1 1{g j (u j ; θ) ≤ 0} - i-1 j=1 χ (g j (u j ; θ)) n j=i 1{g j (u j ; θ) ≤ 0} p du ≤n p-1 n i=1 (0,1) n i j=1 χ (g j (u j ; θ)) n j=i+1 1{g j (u j ; θ) ≤ 0} - i-1 j=1 χ (g j (u j ; θ)) n j=i 1{g j (u j ; θ) ≤ 0} du ≤n p-1 n i=1 (0,1) χ (g i (u i ; θ)) -1{g i (u i ; θ) ≤ 0} du i ,
where the first inequality holds by applying Jensen's inequality. Then the rest of the proof is straightforward.

GLR for the Ali-Mikhail-Haq Copula.

The generator function for this copula is

ψ θ (x) = 1 -θ e x -θ , θ ∈ [0, 1).
Then, we have

∂ψ θ (x) ∂θ = - 1 e x -θ + 1 -θ (e x -θ) 2 , ∂ψ θ (x) ∂θ = e x (e x -θ) 2 - 2e x (1 -θ) (e x -θ) 3 , and ψ θ (x) = - e x (1 -θ) (e x -θ) 2 , ψ θ (x) = 2e 2x (1 -θ) (e x -θ) 3 - e x (1 -θ) (e x -θ) 2 .
Y θ has a geometric distribution with parameter θ, with probability mass function

(1 -θ)θ y-1 for y = 1, 2, . . . , so the LR term is

∂ log f Y (y; θ) ∂θ = - 1 1 -θ + y -1 θ .
The weight functions in the GLR estimator are

Appendix B: Supplements of Section 5

Distribution Sensitivities for a Linear Model For the numerical experiments in Section 5.1 of the main body of the paper, we estimate the sensitivity of E[ϕ(g(U ; θ))] with respect to θ at z, as a function of z from z = 0.01 to 10 with a step size of 0.01. We perform 10 6 independent simulation runs with GLR. Figure 1 This alternative estimator leads to a comparable performance to the one in Section 4.2.

The setting of the problem is the same as that in Section 5.2 of the main body of the paper. Similar to Figure 2 in the example in Section 5.1 of the main body of the paper, Figure 3 presents confidence intervals using a batch of 10 4 independent simulation replications and the coverage rates of the 0.9-confidence intervals by 10 4 independent macro experiments for quantiles q α , α = 0.1, . . . , 0.9.

Again the coverage ratios of the 0.9-confidence intervals match the target value statistically.

Figure 4 shows that the variances of FDC(0.01) and FDC(0.1) are much larger than those of GLR, and variances of FDC(0.01) are about 10 times of those of FDC(0.1) throughout the curve.

Sensitivities of Credit Risk Derivatives

For CDOs, Table 1 shows the sensitivity estimates for sample sizes m = 10 4 , m = 10 5 , and m = 10 6 under the Ali-Mikhail-Haq copula, with θ = 0.5. The numerical observations are similar to those in Section 5.4 of the main body of the paper. Again for sample size m = 10 6 , the estimate with FDC(0.1) falls outside of the 90% confidence interval of the GLR estimate, whereas the estimate with GLR lies in the 90% confidence interval of the FDC(0.01) estimate with sample size m = 10 7 , which is -0.193 ± 0.002.

For BDS, we estimate the sensitivities of the expectation of the discounted payment to the fifth default of 10 assets made by the seller V value (τ ) = L (5) exp(-rτ (5) )1{τ (5) < T }. The marginal

  in which case a smoothed function ϕ (•) can be constructed explicitly. The performance function in the distribution sensitivities discussed in Section 4.1 is an indicator function. The distribution sensitivities for the completion time in an SAN in Section 5.2 and the sensitivities of a control chart in Section 5.3 have performance functions which are products of n indicators. The proof of the following proposition can be found in the online Appendix A. Proposition 1 Consider the stochastic model

Figure 1

 1 Figure 1Density curve estimated by GLR in the example of Section 5.1.

Figure 1

 1 Figure 1(a) shows the estimated density curve, and Figure 1(b) presents the sample variance curves of three distribution sensitivity estimators: G 2,1 (•) (GLR-1), G 2,2 (•) (GLR-2), and a combined GLR estimator that minimizes the variance of the linear combination of GLR-1 and GLR-2

Figure 3 A

 3 Figure 3A SAN with six activities.

Figure 4

 4 Figure 4Density estimation by GLR in the SAN example of Section 5.2.

  1) n ϕ(g(u; θ)) d(u; θ) du, which completes the proof. Proof of Theorem 1 under conditions (A.2)-(A.4) Proof. From Peng et al. (2018), there exists a sequence of bounded functions |ϕ L (x)| ≤ ϕ(x) such that lim L→∞ ϕ L (x) = ϕ(x), and there exists a sequence of bounded and smooth functions ϕ ,L (•) such that lim L→∞ ϕ ,L -ϕ L p = 0, where p > 1, and h p := R n |h(x)| p dx 1/p . Except for replacing ϕ (•) with ϕ ,L (•), the procedures before (1) are the same as in the proof for Theorem 1 under conditions (A.1)-(A.3

  + r i (u; θ, y) = lim u i →1 - r i (u; θ, y) = 0.

  (a) shows the curves of the distribution sensitivities with respect to θ, and Figure1(b) presents the sample variance curves of three distribution sensitivity estimators: G 1,1 (•) (GLR-1), G 1,2 (•) (GLR-2), and a combined GLR estimator minimizing variance in a family of linear combinations of GLR-1 and GLR-2 given by (GLR-Opt).

Figure 1

 1 Figure 1Distribution sensitivities with respect to θ estimated by GLR in the example of Section 5.1 in the main body of the paper.

Figure 2

 2 Figure2presents the variance curves for two distribution sensitivities estimated by FDC using a batch of 10 6 independent simulation replications. FDC(δ) denotes the FDC with perturbation

Figure 2

 2 Figure 2Variance curves for distribution sensitivities estimated by FDC in the example of Section 5.1 of the main body the paper.

  Coverage rates of confidence intervals.

Figure 3

 3 Figure390%-confidence intervals estimated by GLR for quantiles and coverage rates in example of Section 5.2 of the main body of the paper.

Figure 4

 4 Figure 4Variance curves for density estimates by FDC in the example of Section 5.2 of the main body of the paper.

Table 1

 1 Derivatives of control charts with respect to upper control limit θ = 1, based on 10 6 independent replications (mean ± standard error).

			GLR	FD(0.01)	FD(0.1)
		∂E[N ]/∂θ	8.6 ± 0.2	8.6 ± 0.4	10.5 ± 0.04
	∂E[(N -χ/∆ ) + ]/∂θ 8.8 ± 0.07	8.6 ± 0.4	10.5 ± 0.04
		m = 10 4		m = 10 5	m = 10 6
	GLR	-0.187 ± 0.01	-0.185 ± 0.004	-0.181 ± 0.002
	FDC(0.1)	-0.171 ± 0.007	-0.177 ± 0.003	-0.176 ± 0.002
	FDC(0.01)	-0.155 ± 0.07	-0.186 ± 0.02	-0.189 ± 0.006

Table 2

 2 Sensitivity estimates of CDO with 10 assets governed by the Clayton copula with θ = 0.5 based on 10 2

experiments (mean ± standard error).

Table 3

 3 Sensitivity estimates of BDS with 10 assets governed by the Ali-Mikhail-Haq copula with θ = 0.5 based on 10 2 experiments (mean ± standard error).

		n = 10 4	n = 10 5	n = 10 6
	GLR	0.051 ± 0.002	0.050 ± 0.001	0.051 ± 0.0009
	FDC(0.1)	0.048 ± 0.003	0.052 ± 0.001	0.051 ± 0.001
	FDC(0.01)	0.074 ± 0.03	0.039 ± 0.01	0.049 ± 0.003

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grants 71901003, 71720107003, 71690232, 91846301, 71790615, by the National Science Foundation (NSF) under Grant CMMI-1434419, by Discovery Grant RGPIN-2018-05795 from NSERC-Canada, and by the Key Research and Development Program of Beijing Municipal Science and Technology Commission.

Appendix A: Theoretical Supplements Proof of Theorem 1 under conditions (A.1)-(A.3) Proof. We have distributions of the defaults and the parameters are set in the same way as those in Section 5.1 in the main body of the paper. Tables 2 and3 show the respective Clayton copula and Ali-Mikhail-Haq copula sensitivity estimates with θ = 0.5 for sample sizes m = 10 4 , m = 10 5 , and m = 10 6 .

We also implement FDC(0.01) with sample size m = 10 7 , which leads to 0.034 ± 0.001 for the Clayton copula and 0.051 ± 0.002 for the Ali-Mikhail-Haq copula. The results are similar to those for estimating CDO sensitivities.


Joaquim Casellas
email: joaquim.casellas@uab.cat

Cecilia Esquivelzeta

Andrés Legarra

Keywords: accuracy, genomic selection, mutation, relationship matrix 2

come

Short Communication

Current availability of massive genotyping technologies based on SNP has provided high-throughput data for the prediction of genomic breeding values by appropriate analytical approaches [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF]. These genomic breeding values try to capture the additive genetic variability from QTL, although recent studies suggested that part of the genetic variability could escape inclusion in genomic predictions [START_REF] Maher | Personal genomes: The case of the missing heritability[END_REF]. Indeed, [START_REF] Casellas | Effect of mutation age on genomic predictions[END_REF] showed that the statistical performance of genomic prediction models was highly influenced by both the selective history of the trait and the age of each mutation involved in a given QTL, recent mutations systematically escaping inclusion in analyses. Note that genomic analyses would typically not take into account the effect of new mutations. Genomic BLUP (GBLUP) models [START_REF] Nejati-Javaremi | Effect of total allelic relationship on accuracy of evaluation and response to selection[END_REF] account for additive genetic relationships by the genomic relationship matrix (G;[START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF]. Despite the expectation of G being the regular numerator relationship matrix (A), several authors highlighted the gain by using G within the context of animal breeding [START_REF] Legarra | Performance of genomic selection in mice[END_REF][START_REF] Vanraden | Reliability of genomic predictions for North American Holstein bulls[END_REF]. Nevertheless, G does not properly account for new mutations contributing novel additive genetic variability after the founder generation. This limitation was overcome by [START_REF] Wray | Accounting for mutation effects in the additive genetic variance-covariance matrix and its inverse[END_REF] for the standard matrix A, and the main objective of this short communication was to adapt and evaluate the approach of [START_REF] Wray | Accounting for mutation effects in the additive genetic variance-covariance matrix and its inverse[END_REF] and [START_REF] Casellas | Within-generation mutation variance for litter size in inbred mice[END_REF] into a genomics context.

Although G can be constructed by several different methods [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF], we focused as an example on the first approach outlined by [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF], given the benefits reported by [START_REF] Toro | A note on the rationale for estimating genealogical coancestry from molecular markers[END_REF]. Nevertheless, our approach can be straightforwardly adapted to other algorithms for G. In this research, procedures for calculating G were adapted to accommodate the occurrence of new mutations in the autosomal genome [i.e., the mutational genomic relationship matrix (G m )]. Note that G m can be defined as G m = G 1 + G 2 + . . . + G g , G q being a genomic relationship matrix where individuals from generation q were treated as founders and the rows and columns linked to ancestors from previous generations were fixed to 0. Assuming σ m 5399 Two different hierarchical mixed linear models were assumed and alternative parameterization of mixed model equations by [START_REF] Henderson | Applications of Linear Models in Animal Breeding[END_REF] was adopted to avoid the inversion of matrices G and G m . The standard GBLUP parameterization (model A) was y = μ + Z(Ga ) + e, where y was the vector of phenotypic records, μ was the vector of population means, e was the vector of residuals, and a = (Ga ) was the vector of additive genomic effects with appropriate incidence matrix Z. From a Bayesian point of view, the conditional distribution of y was assumed multivariate normal i.e., MVN e μ + ( )

⎡ ⎣ ⎢ ⎤ ⎦ ⎥ { } Z Ga I • , , σ 2
and the a priori dis-

tribution of a was p M V N a a a G 0 G • | ~, . σ σ 2 1 2 ( ) ( )
-Note that I was an incidence matrix with dimensions equal to the number of elements in y, and σ a 2 and σ e 2 were additive genomic and residual variances, respectively. Remaining a priori distributions were assumed flat with appropriate boundaries. Estimated allelic frequencies in the base generation were used to compute G, as used in North American dairy cattle genomic evaluations (P. M. VanRaden, Animal Improvement Programs Laboratory, ARS, USDA, Beltsville, MD, personal communication). On the other hand, an expanded version of the previous model (model AM) was proposed to account for new mutational genomic effects (m) as follows: σ As for model A, the remaining a priori distributions were assumed flat.

y = μ + Z(Ga ) + Z(G m m ) + e,
Simulations relied on a diploid genome composed by 10 chromosomes (100 cM each) with 2,000 SNP and 200 QTL per chromosome. Each replicate evolved during 1,000 non-overlapping generations with effective size (N e ) 100. Founder individuals were homozygous throughout the whole genome for the wild-type allele (i.e., allele 1), and mutation rates of 2.5 × 10 -3 (SNP) and 2.5 × 10 -5 (QTL; [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF] were applied during the first 1,000 generations, switching the allele state from 1 to 2, or vice versa. After that, generations expanded to 500 individuals, with 50 males and 450 females under random mating; mutation rate for SNP changed to 2.5 × 10 -8 [START_REF] Hickey | Simulated data for genomic selection and genome-wide association studies using a combination of coalescent gene drop methods[END_REF], whereas this parameter did not modify for QTL markers. Analyses were performed on 3 different scenarios with 3 (G3), 5 (G5), or 10 (G10) generations with both genomic (i.e., genotypes from polymorphic SNP) and phenotypic data (see below), and a last generation only contributing genomic data. Although the number of generations was a rough simplification of the real structure of current livestock populations, scenario G3 could be viewed as a reasonable approximation to a long-generation interval population with few generations contributing data (e.g., dairy cattle), whereas scenarios G5 and G10 tried to mimic medium-(e.g., intensive pig breeds) and short-generation interval populations (e.g., rabbits and broilers), respectively. On the other hand, 3 different heritabilities (h 2 ) were assumed for the phenotypic trait (h 2 = 0.1, 0.25 and 0.5), with σ e 2 = 1, and QTL effects were appropriately scaled. Note that the wild-type QTL had a null contribution to the phenotype, whereas the additive effect of the mutate allele was sampled from a gamma distribution with shape and scale parameters 0.4 and 1.6 [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF]. Moreover, the sign of every QTL effect was sampled to be positive or negative with probability 0.5. A total of 100 replicates were analyzed for each scenario and heritability. Model performance was evaluated in terms of goodness of fit by the deviance information criterion (DIC) statistic [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF] and the precision of predicted breeding values (u), being u = a (model A) or u = a + m (model AM). Precision was calculated for each replicate as the correlation coefficient (r u,û ) between simulated and predicted breeding values. Within a Bayesian framework, all replicates were analyzed under model A and model AM. For each analysis, a unique Monte Carlo Markov chain with 100,000 iterations was launched, after discarding the first 10,000 iterations as burn-in [START_REF] Raftery | How many iterations in the Gibbs sampler? Pages 763-773 in Bayesian Statistics IV[END_REF].

Given that this research focused on the analysis of simulated populations, a preliminary examination of the estimated σ m 2 became mandatory to corroborate the validity of further results about model performance. As shown in = σ σ , where σ T 2 is the total phenotypic variance) averaged between 0.04 and 0.25%, agreeing with the range of values previously reviewed by [START_REF] Lynch | The rate of polygenic mutation[END_REF] and [START_REF] Houle | Comparing mutational variabilities[END_REF]; indeed, even the most extreme estimates fit the range, although some of them were close to the lower boundary. This highlighted the conservative behavior of the simulation procedure where the impact of new mutations was not artificially inflated, preventing unrealistic advantages for the model from accounting for new mutational variance.

Both GBLUP models were compared in terms of goodness of fit by the DIC statistic, which suggested important influences from heritability and pedigree depth (Table 2). Note that the DIC evaluates both model fit and complexity, and average differences larger than 3 to 5 DIC units are generally assumed statistically relevant [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF]. Small heritabilities and shallow pedigrees favored model A, whereas moderate-to-high heritabilities or deeper pedigrees took benefit of the most complex structure of model AM accounting for new mutational variability. Similar evidences were revealed by r u,û , where higher accuracies (P < 0.05) were reported for model AM under simulation scenarios with moderate-to-high heritabilities or deeper pedigrees (Tables 3 and4), albeit the gains were marginal. These results suggested that model AM could be ineffective to properly capture mutational variability when lowly heritable traits were analyzed on shallow pedigrees, whereas heritability and pedigree depth showed synergistic behavior in less extreme scenarios [i.e., traits with moderate-to-high heritabilities (even under shallow pedigrees) or moderate-to-deeper pedigrees (even for lowly heritable traits) would marginally benefit from the inclusion of new mutational effects in GBLUP models].

Regardless of the GBLUP parameterization, r u,û showed very similar trends when comparing different heritabilities, pedigree depths, or individuals with and without phenotypic data. The accuracy of the genomic predictions significantly increased with heritability as anticipated by [START_REF] Calus | Accuracy of genomic selection using different methods to define haplotypes[END_REF], and demonstrated by previous studies without the inclusion of mutational effects [START_REF] Meuwissen | Accuracy of breeding values of 'unrelated' individuals predicted by dense SNP genotyping[END_REF][START_REF] Zhong | Factors affecting accuracy from genomic selection in populations derived from multiple inbred lines: A barley case study[END_REF]. Accuracy was higher for phenotyped individuals than for those without contributing phenotypic data to the analysis. Indeed, this decline in precision for the unphenotyped individuals from recent generations has been widely demonstrated in the scientific literature (Meuwissen = σ σ , where σ a 2 is the additive variance and σ T 2 is the total phenotypic variance.

3

h m m T 2 2 2
= σ σ , where σ m 2 is the mutational variance and σ T 2 is the total phenotypic variance. [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Legarra | Performance of genomic selection in mice[END_REF]. It is important to note that r u,û increased with the number of generations involved in G, this being probably due to the growth of the genomic/phenotypic data set [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Goddard | Genomic selection: Prediction of accuracy and maximisation of long term response[END_REF] and the rise of the degree of relationship among individuals from late generations in small populations [START_REF] Habier | The impact of genetic relationship information on genome-assisted breeding values[END_REF]. In summary, maximum DIC and r u,û differences were reported in scenario G10 and h 2 = 0.5. The accuracy of genomic breeding values increased 0.57% (phenotyped individuals) or 1.07% (unphenotyped individuals) under model AM. This suggests that model AM could be of interest for genomic selection in livestock populations, although each case must be evaluated on its merits, depending on population structure and the genetic model underlying the trait. Moreover, other genetics-related phenomena such as selection, population stratification, or even genetic drift must be carefully considered because they could vary the real contribution of mutational effects.

A delicate point is the consideration of generations when they overlap. [START_REF] Wray | Accounting for mutation effects in the additive genetic variance-covariance matrix and its inverse[END_REF] pointed out that the decisive notion in her model was time span, not generations; the model applies equally well to overlapping generations pedigrees. In this case, generations would be cohorts born in the same time span (e.g., year-byyear groups). The only requirement is that no individual within one cohort can be descendant from another from the same cohort. Even in our model using genomic relationships, a notion of pedigree or time flow is absolutely For this partitioning of G m , the diagonal element of individual i in G i is fixed to 0, and this is so because the covariance matrix G i must be centered around this element. The G i -specific mutational genetic value of i is exactly defined as 0, has no uncertainty, and this individual becomes the genetic basis of the genetic covariance matrix G i . However, the ith element in matrix G m will not be 0 because of the ith diagonal contributions of remaining cohorts.

To the best of our knowledge, this is the first approach contributing a reliable way to account for new additive mutational variability in genomic prediction models, with this GBLUP parameterization being able to provide small advantages in terms of accuracy of the genomic predictions, a way of handling mutational variance in selection schemes, and possibly more accurate estimates of mutational variance than those based on pedigrees. Although the relevance of new mutations could be small when individuals from few generations are involved (i.e., first steps in genomic selection programs), this model could be of special relevance as times goes on and genotyped individuals accumulate generation by generation. Also, it allows better handling of new genetic variability, such as the selection of "original" animals or loci, which are penalized by standard animal breeding techniques [START_REF] Goddard | Genomic selection: Prediction of accuracy and maximisation of long term response[END_REF].

Table 1

 1

			, estimated mutational heritabilities
	(h m 2	m 2	T 2

Table 1 .

 1 Estimated heritabilities under the different simulation scenarios by assuming both direct and mutational additive genetic sources of variation

				Standard heritability 2	Mutational heritability 3 (× 100)
	Item 1			Mean	Range	Mean	Range
	Scenario G3 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	0.12 0.24 0.48	0.04-0.21 0.13-0.35 0.37-0.66	0.04 0.09 0.14	0.01-1.57 0.02-2.30 0.04-2.75
	Scenario G5 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	0.12 0.23 0.52	0.05-0.16 0.17-0.31 0.43-0.66	0.04 0.11 0.19	0.02-1.50 0.04-1.99 0.07-2.32
	Scenario G10 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	0.11 0.25 0.51	0.06-0.15 0.19-0.29 0.42-0.58	0.06 0.14 0.25	0.02-0.91 0.07-1.09 0.13-0.89
	1 G3, G5, and G10 = 3, 5, and 10 generations, respectively.
	2 h a 2	a 2	T 2

Table 2 .

 2 Average deviance information criterion (DIC) estimates (means ± SE)

	Item 1	Model A 2	Model AM 3	Difference 4
	Scenario G3 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	4,807 ± 4 4,829 ± 4 4,877 ± 5	4,815 ± 4 4,832 ± 4 4,867 ± 4	-8 ± 2** -3 ± 3 NS 10 ± 4*
	Scenario G5 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	7,431 ± 5 7,472 ± 5 7,501 ± 6	7,429 ± 5 7,462 ± 5 7,479 ± 7	2 ± 3 NS 10 ± 3* 22 ± 4***
	Scenario G10 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	13,192 ± 8 13,236 ± 9 13,260 ± 9	13,183 ± 8 13,216 ± 8 13,225 ± 9	9 ± 4* 20 ± 4*** 35 ± 5***

Table 3 .

 3 Correlation coefficients (× 100; mean ± SE) between simulated and predicted breeding values for individuals with phenotypic data Consider thus cohorts i = 1, 2, . . . , q as the starting point. According to properties of G[START_REF] Gianola | Additive genetic variability and the Bayesian alphabet[END_REF] and under Hardy-Weinberg equilibrium, G i would contain an average of 1 in the diagonal and 0 off-diagonal for individuals from cohort i. Within this context, we could reduce the notion of cohort to the extreme where each individual defines a new cohort and, therefore,G m = G 1 + G 2 + . . . + G t ,where t is the total number of individuals in the pedigree file. Only descendants of the ith individual would have values different from 0 in G i , as previously suggested by García-Cortés et al. (2010) for matrix A.

	Item 1	Model A 2	Model AM 3	Difference 4
	Scenario G3 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	52.49 ± 0.23 64.25 ± 0.24 82.19 ± 0.24	52.51 ± 0.22 64.29 ± 0.25 82.27 ± 0.26	0.02 ± 0.05 NS 0.04 ± 0.04 NS 0.08 ± 0.05 NS
	Scenario G5 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	57.66 ± 0.24 69.21 ± 0.25 85.12 ± 0.23	57.69 ± 0.25 69.31 ± 0.23 85.38 ± 0.25	0.03 ± 0.04 NS 0.10 ± 0.04* 0.26 ± 0.03***
	Scenario G10 h 2 = 0.10 h 2 = 0.25 h 2 = 0.50	63.03 ± 0.22 73.81 ± 0.21 88.60 ± 0.19	63.10 ± 0.23 74.03 ± 0.25 89.17 ± 0.18	0.07 ± 0.05 NS 0.24 ± 0.04*** 0.57 ± 0.03***
	1 G3, G5, and G10 = 3, 5, and 10 generations, respectively.
	2 Genomic BLUB (GBLUP) model accounting for standard additive genetic effects.
	3 Genomic BLUP model accounting for both standard additive and new mutational genetic effects.
	4 Differences between model A and model AM were tested by a paired-samples t-test.
	*P < 0.05; ***P < 0.001.

Journal of Dairy Science Vol. 96No. 8, 2013

ACKNOWLEDGMENTS

The research contract of J. Casellas was partially financed by Spain's Ministerio de Ciencia e Innovación (program Ramón y Cajal, reference RYC-2009-04049). This research was funded by grant AGL2010-21176/ GAN (Ministerio de Ciencia e Innovación, Madrid, Spain).