<figure><img src="image_1.png" /></figure>

HAL
open science

# Efficiency of Estimating Smooth Functions of Parameters in Rare-Event Contexts

Marvin K Nakayama, Bruno Tuffin

► To cite this version:

Marvin K Nakayama, Bruno Tuffin. Efficiency of Estimating Smooth Functions of Parameters in Rare-Event Contexts. 2025. hal-05338278

HAL Id: hal-05338278

https://inria.hal.science/hal-05338278v1

Preprint submitted on 30 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficiency of Estimating Smooth Functions of Parameters in Rare-Event Contexts

Marvin K. Nakayama

Bruno Tuffin

New Jersey Institute of Technology
Computer Science Department
Newark, New Jersey, USA

Inria, Univ Rennes, CNRS, IRISA
Rennes, France

October 30, 2025

## Abstract

Consider an estimand that is a known smooth function of an unknown vector of parameters, such as means, probabilities, quantiles, etc. Some but not necessarily all of the parameters relate to rare events, which typically make them difficult to estimate with naive Monte Carlo, so variance-reduction techniques must be applied to obtain efficient estimators of those parameters. We develop conditions ensuring that the overall estimand can be efficiently estimated, providing both sufficient and necessary- and sufficient conditions. The paper illustrates the applicability of the theory through many examples and numerical results.

## 1 Introduction

Many rare-event performance/risk metrics can be expressed as a known function of unknown parameters, such as means, probabilities, or quantiles. Some but not necessarily all of the parameters relate to the rarity of events of interest, making them difficult to estimate via naive Monte Carlo (MC) simulation, i.e., by using a sample generated from the (direct) simulation of the underlying random variables or stochastic process. For example, for a stable GI/GI/1 queue, the expected time for the queue length to first reach a specified threshold (e.g., a buffer overflow) can be written as a ratio, where the numerator is the expectation of the minimum of the busy-cycle length and the time for the queue length to hit the threshold, and the denominator is the probability of the queue length hitting the threshold within a busy cycle. When the threshold is large, the denominator is a rare-event probability, and estimating it via naive MC is typically inefficient, motivating the application of variance-reduction techniques (Asmussen and Glynn, 2007, Chapters V and VI), such as importance sampling. In contrast, the numerator can typically be estimated well through naive MC.

Building on prior work in (Nakayama and Tuffin, 2023), which provided sufficient conditions to ensure the efficiency of estimating a known function of unknown *means* in a rare-event context, we now expand the previous theory by considering the estimation of a function of unknown *parameters*, which may include means, probabilities, quantiles, etc. We assume the parameters have MC estimators obeying a joint *central limit theorem* (CLT) as a budget parameter $n \to \infty$. We obtain a plug-in estimator by computing the known function at the vector of estimators. While (Nakayama and Tuffin, 2023) provides only sufficient conditions to ensure that the plug-in estimator is efficient, as measured by its relative error, we further identify various settings in which the conditions are also necessary, which give a deeper understanding by allowing one to prove when certain estimators are *not* efficient.

An outline of the rest of the paper follows. Section 2 reviews different notions of efficiency for simulation estimators for problems involving rare events. Section 3 describes the setting of the estimand as a function of parameters, and illustrates in Example 3 that efficiency of the parameters' and of the estimand's estimators are not equivalent: the estimand can be efficiently estimated despite inefficient estimators of the parameters,
and efficient estimators of the parameters do not guarantee an efficient estimand estimator. Section 3 then presents Theorem 1, which provides sufficient conditions to ensure the estimand can be estimated efficiently. Section 4 gives several examples that demonstrate the broad applicability of Theorem 1. But the sufficient conditions of Theorem 1 provide only a partial picture of efficiency since its conditions are not necessary in general, so Section 5 identifies situations in which the conditions are also necessary. Section 6 gives numerical results complementing the theory, and concluding remarks appear in Section 7.

## 2 Efficiency of Simulation Estimators

We first review various notions of statistical efficiency of estimators in rare-event simulations. For a given stochastic model, let $\alpha$ be an estimand, such as a mean, quantile, or other unknown parameter, and the goal is to estimate $\alpha$ via a Monte Carlo method, e.g., naive MC (also known as simple random sampling) or a variance-reduction technique. Let $\hat{\alpha}_n$ be an MC estimator of $\alpha$ based on a computation budget (e.g., sample size) $n$, and suppose that the following CLT holds:

$$ \sqrt{n}[\hat{\alpha}_n - \alpha] \Rightarrow \mathcal{N}(0, \tau^2) \quad \text{as } n \to \infty, \qquad (1) $$

where $\Rightarrow$ represents convergence in distribution (Billingsley, 1995, Section 25), $\mathcal{N}(a, b^2)$ denotes a normal random variable with mean $a$ and variance $b^2$, and $\tau^2 \in (0, \infty)$ is the CLT's asymptotic variance, which depends on the MC method applied and the form of the estimator $\hat{\alpha}_n$. If $\hat{\tau}_n^2$ consistently estimates $\tau^2$ in the sense that $\hat{\tau}_n^2 \Rightarrow \tau^2$ as $n \to \infty$, then we obtain $[\hat{\alpha}_n \pm 1.96\hat{\tau}_n/\sqrt{n}]$ as a 95% confidence interval (CI) for $\alpha$. Assuming that $\alpha \neq 0$, the CI's relative half-width for a fixed $n$ is approximately $1.96 \text{RE}[\hat{\alpha}_n]/\sqrt{n}$, where

$$ \text{RE}[\hat{\alpha}_n] = \frac{\tau}{|\alpha|} \qquad (2) $$

denotes the relative error (RE; e.g., (L'Ecuyer et al., 2010)) of the estimator $\hat{\alpha}_n$. Thus, obtaining a CI with relative half-width $\epsilon > 0$ requires a computation budget that is roughly

$$ n = (1.96 \text{RE}[\hat{\alpha}_n]/\epsilon)^2, \qquad (3) $$

so $\text{RE}[\hat{\alpha}_n]$ provides a measure of the efficiency of the estimator $\hat{\alpha}_n$.

**Example 1** (Estimating an expectation via a sample mean). Suppose that $\alpha = \mathbb{E}[Z]$ is the expectation of a random variable $Z$, and we estimate $\alpha$ via the sample average $\hat{\alpha}_n = (1/n) \sum_{k=1}^n Z_k$ of $n$ i.i.d. copies $Z_1, Z_2, \dots, Z_n$ of $Z$, so the computation budget $n$ is the sample size. If $\text{Var}[Z] \in (0, \infty)$, where $\text{Var}[\cdot]$ denotes the variance operator, then the estimator $\hat{\alpha}_n$ obeys the CLT (1) with $\tau^2 = \text{Var}[Z]$ (Billingsley, 1995, Theorem 27.1). This framework accommodates importance sampling (IS, Section V.1 of (Asmussen and Glynn, 2007)) as follows. For $Q$ a random variable with distribution $F_0$, suppose that $\alpha = \mathbb{E}_0[Q]$, where $\mathbb{E}_0[\cdot]$ represents the expectation operator under $F_0$. Let $F_1$ be another distribution such that the measure corresponding to $F_0$ is absolutely continuous with respect to that of $F_1$, and let $\mathbb{E}_1[\cdot]$ denote the expectation operator under $F_1$. Let $dF_0/dF_1$ be the likelihood ratio (LR), also known as the Radon-Nykodim derivative, which we assume can be computed. Then we may write $\alpha = \int Q dF_0 = \int Q \frac{dF_0}{dF_1} dF_1 = \mathbb{E}_1[QL]$, so $\alpha = \mathbb{E}[Z]$ with $Z = QL$ and $\mathbb{E}$ corresponding to $\mathbb{E}_1$. $\square$

(L'Ecuyer et al., 2010) examines the RE of only estimators fitting Example 1 (estimating a mean through a sample average of i.i.d. observations of $Z$), so (L'Ecuyer et al., 2010) defines RE slightly differently by replacing $\tau$ in (2) with the exact standard deviation $\sqrt{\text{Var}[Z]}$. But for estimation settings outside of Example 1, the exact variance may be quite complicated or intractable, such as when estimating smooth functions of parameters (Section 3) or quantiles (Example 5). This is the reason why (2) uses instead the square root of the CLT's asymptotic variance in defining RE, as is also done in (Shahabuddin, 1994), (Bassamboo et al., 2008), and (Li et al., 2024).

We focus primarily on situations when $\alpha$ is related to a rare event, in which case naive MC can result in an inefficient estimator of $\alpha$. A typical situation arises in Example 1 when $Z = I(W \in A)$, where $W$ is a
real-valued random variable, $A \subset \mathbb{R}$ represents a set of outcomes of interest, and $I(\cdot)$ denotes the indicator function, which equals 1 (resp., 0) when its argument is true (resp., false). Then $\alpha = \mathbb{E}[Z] = \mathbb{E}[I(W \in A)] = \mathbb{P}(W \in A)$ is the probability of $W$ lying in set $A$, and we assume that $\{W \in A\}$ is a rare event in the sense that $\alpha$ is close to 0. We study this mathematically by considering a family of stochastic models indexed by a *rarity parameter* $r > 0$, so $\alpha = \alpha_{\langle r \rangle}$, $Z = Z_{\langle r \rangle}$, $W = W_{\langle r \rangle}$, and $A = A_{\langle r \rangle}$, where we assume that $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$. For example, for a stable GI/GI/1 queue, let $W_{\langle r \rangle}$ equal the maximum queue length before the queue first empties, and let $A_{\langle r \rangle} = \{r, r+1, \dots\}$ correspond to a set of “large” queue lengths. Then $\alpha_{\langle r \rangle}$ is the probability of the queue-length process reaching level $r$ within the first busy period, and $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$ (Sadowsky, 1991, Theorem 1). In this case, (1) holds for each fixed $r$, where $\tau^2 \equiv \tau^2_{\langle r \rangle} = \alpha_{\langle r \rangle}(1 - \alpha_{\langle r \rangle})$, so $\tau^2_{\langle r \rangle} \to 0$ as $r \to \infty$. While this indicates that the error in $\alpha_{\langle r \rangle}$ shrinks to 0 as $r \to \infty$, its RE in (2) blows up:

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\tau_{\langle r \rangle}}{\alpha_{\langle r \rangle}} = \frac{\sqrt{\alpha_{\langle r \rangle}(1 - \alpha_{\langle r \rangle})}}{\alpha_{\langle r \rangle}} = \frac{\sqrt{1 - \alpha_{\langle r \rangle}}}{\sqrt{\alpha_{\langle r \rangle}}} \to \infty \quad \text{as } r \to \infty. $$

Therefore, the sample size in (3) required to obtain a CI with a specified half-width $\epsilon$ has to grow with $r$. This exemplifies the primary difficulty of rare-event simulations.

Applying a variance-reduction technique may lead to an estimator $\hat{\alpha}_{\langle r \rangle, n}$ whose relative error vanishes, remains bounded, or grows slowly as $r \to \infty$. Such methods have been developed for a variety of estimation problems involving rare events (e.g., (Rubino and Tuffin, 2009b), (Asmussen and Glynn, 2007, Chapter VI), (Blanchet and Lam, 2011), and (Budhiraja and Dupuis, 2019)). To explain this, we introduce the notation $b_1(r) = O(b_2(r))$ as $r \to \infty$ to mean that $|b_1(r)/b_2(r)|$ is bounded for all sufficiently large $r$. Also, $b_1(r) = \Theta(b_2(r))$ denotes that both $b_1(r) = O(b_2(r))$ and $b_2(r) = O(b_1(r))$ as $r \to \infty$. Now if $\tau^2_{\langle r \rangle} = O(\alpha_{\langle r \rangle})$ as $r \to \infty$, then

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\sqrt{O(\alpha_{\langle r \rangle}^2)}}{\alpha_{\langle r \rangle}} = \frac{O(\alpha_{\langle r \rangle})}{\alpha_{\langle r \rangle}} = O(1) \quad \text{as } r \to \infty, \quad (4) $$

and we say the estimator has *bounded relative error* (BRE). In this case, the sample size $n$ needed to obtain a CI for $\alpha_{\langle r \rangle}$ with a fixed relative width remains bounded as the event $\{W_{\langle r \rangle} \in A_{\langle r \rangle}\}$ of interest becomes rarer with $r \to \infty$. Instead suppose now that $\tau^2_{\langle r \rangle} = o(\alpha_{\langle r \rangle})$ as $r \to \infty$, where $b_1(r) = o(b_2(r))$ means that $b_1(r)/b_2(r) \to 0$ as $r \to \infty$. Then

$$ \mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\sqrt{o(\alpha_{\langle r \rangle}^2)}}{\alpha_{\langle r \rangle}} = \frac{o(\alpha_{\langle r \rangle})}{\alpha_{\langle r \rangle}} = o(1) \quad \text{as } r \to \infty, \quad (5) $$

and the estimator $\hat{\alpha}_{\langle r \rangle, n}$ exhibits a stronger property than BRE known as *vanishing RE* (VRE).

In some rare-event settings, the best performance that a user may practically be able to hope for is that the estimator $\hat{\alpha}_{\langle r \rangle, n}$ possesses a slightly weaker property than BRE known as *logarithmic efficiency*, also called *asymptotic optimality*. To describe this, for each constant $\zeta \in \mathbb{R}$, let

$$ \mathrm{RE}_{\zeta}[\hat{\alpha}_{\langle r \rangle, n}] = \frac{\tau_{\langle r \rangle}}{|\alpha_{\langle r \rangle}|^{\zeta}}, \quad (6) $$

which for $\zeta = 1$ is the same as RE in (2). We say that the estimator $\hat{\alpha}_{\langle r \rangle, n}$ of $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$ is *logarithmically efficient* (LE) if

$$ \text{for all } \delta > 0, \mathrm{RE}_{1-\delta}[\hat{\alpha}_{\langle r \rangle, n}] = o(1) \quad \text{as } r \to \infty. \quad (7) $$

This property may be satisfied when applying importance sampling to estimate a probability $\alpha_{\langle r \rangle}$ that shrinks to 0 at an exponential rate established by a large-deviations principle (Billingsley, 1995, p. 148). For example, this occurs for the queue-length process of a stable GI/GI/1 queue with light-tailed interarrival
and service times when $\alpha_{\langle r \rangle}$ represents the probability that the queue length hits a large threshold $r$ before emptying out and the so-called exponential tilting is used as importance sampling (Sadowsky, 1991). When $\alpha_{\langle r \rangle} > 0$ with $\alpha_{\langle r \rangle} \to 0$ and $\tau_{\langle r \rangle} \to 0$ as $r \to \infty$, (7) can be equivalently stated as $\liminf_{r \to \infty} \frac{|\ln \tau_{\langle r \rangle}|}{|\ln(\alpha_{\langle r \rangle})|} \ge 1$, whence the terminology LE. The previous discussion focuses on estimands satisfying $\alpha_{\langle r \rangle} \to 0$ as $r \to \infty$, but there are other situations in which $|\alpha_{\langle r \rangle}|$ instead becomes large as $r$ grows. For example, this happens when $\alpha_{\langle r \rangle}$ is the expected hitting time to a large threshold $r$ for the queue-length process a stable GI/GI/1 queue. To handle such settings with $\alpha_{\langle r \rangle} \to \pm\infty$ as $r \to \infty$, we alternatively define LE to denote that

$$ \text{for all } \delta > 0, \quad \text{RE}_{1+\delta}[\hat{\alpha}_{\langle r \rangle, n}] = o(1) \quad \text{as } r \to \infty, \qquad (8) $$

which is equivalent to $\limsup_{r \to \infty} \frac{|\ln \tau_{\langle r \rangle}|}{|\ln(\alpha_{\langle r \rangle})|} \le 1$ when both $\alpha_{\langle r \rangle} \to \infty$ and $\tau_{\langle r \rangle} \to \infty$. LE was originally defined to study the estimation of exponentially decreasing probabilities, and our analysis of LE does not include situations in which $\alpha_{\langle r \rangle} \to \chi_0$ as $r \to \infty$ for a non-zero constant $\chi_0 \in (-\infty, \infty)$ since it is not at all obvious how LE should then be defined.

We note that there are other forms of efficiency discussed in the literature but not considered here. For example, the “work-normalized RE” (Glynn et al., 2009, Section 4.2) further takes into account the expected CPU time to compute $\hat{\alpha}_{\langle r \rangle, n}$. Also, we do not examine efficiency measures based on moments higher than order 2, which can be useful in analyzing the quality of the resulting CI (Glynn et al., 2009, Section 4.4).

In what follows, we sometimes omit the subscript $\langle r \rangle$ on variables to simplify the notation.

## 3 SMOOTH FUNCTIONS OF PARAMETERS

We consider an estimand $\alpha$ having the form

$$ \alpha = g(\boldsymbol{\theta}) \qquad (9) $$

for some known “smooth” (defined in Theorem 1) function $g: \mathbb{R}^d \to \mathbb{R}$ and vector $\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_d)^\top \in \mathbb{R}^d$ of unknown parameters, for some fixed $d \ge 1$, where a superscript $\top$ denotes transpose. Let $\mathcal{I} = \{1, 2, \dots, d\}$ be the set of coordinate indices of $\boldsymbol{\theta}$. For each $i \in \mathcal{I}$, $\theta_i$ is an estimatable parameter, such as a mean (e.g., a probability) or a non-mean quantity (e.g., a quantile). We assume that the value of $\boldsymbol{\theta} = \theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2}, \dots, \theta_{\langle r \rangle, d}^\top$ may depend on the rarity parameter $r$, but $g$ and $d$ do not. Then the estimand $\alpha = \alpha_{\langle r \rangle}$ in (9) also depends on $r$, and we study when $\alpha_{\langle r \rangle}$ can be estimated with BRE as in (4), with VRE as in (5), or with LE as in (7) or (8).

To estimate $\alpha_{\langle r \rangle} = g(\boldsymbol{\theta}_{\langle r \rangle})$ in (9), assume that we have a computing budget (e.g., sample size) $n$ used to construct an estimator

$$ \hat{\boldsymbol{\theta}}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2}, \dots, \hat{\theta}_{\langle r \rangle, n, d})^\top \qquad (10) $$

of $\boldsymbol{\theta}_{\langle r \rangle}$, leading to

$$ \hat{\alpha}_{\langle r \rangle, n} = g(\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}) \qquad (11) $$

as a plug-in estimator of $\alpha_{\langle r \rangle}$. (The computing budget $n$ may represent something other than the sample size in i.i.d. sampling, such as the simulated time of a stochastic process, as we will see later in Example 10 in Section 4.) We impose the following:

**Assumption A1.** For each $r$, the estimator $\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}$ of $\boldsymbol{\theta}_{\langle r \rangle} \in \mathbb{R}^d$ obeys a joint CLT

$$ \sqrt{n} [\hat{\boldsymbol{\theta}}_{\langle r \rangle, n} - \boldsymbol{\theta}_{\langle r \rangle}] \Rightarrow \mathcal{N}_d(\mathbf{0}, \Sigma_{\langle r \rangle}) \quad \text{as } n \to \infty, \qquad (12) $$

where $\mathcal{N}_d(\mathbf{0}, \Sigma_{\langle r \rangle})$ is a $d$-dimensional normal random vector with mean vector $\mathbf{0} = (0, 0, \dots, 0)^\top \in \mathbb{R}^d$ and a finite and positive-definite covariance matrix $\Sigma_{\langle r \rangle} = (\Sigma_{\langle r \rangle, i, j} : i, j \in \mathcal{I})$.
As will be explained later in Example 2, Assumption A1 holds, e.g., when $\theta_{\langle r \rangle}$ is a vector of means that is estimated by a sample average of $n$ i.i.d. observations, but Assumption A1 also allows for more general situations, which will be discussed in Section 4. We call $\Sigma_{\langle r \rangle}$ in (12) the asymptotic covariance matrix in the CLT for the vector estimator $\hat{\theta}_{\langle r \rangle, n}$ in (10). The $d$ components of $\hat{\theta}_{\langle r \rangle, n}$ may be dependent, as is often the case when they are generated from the same simulation, or they may be independent, generated from different simulations applying different MC methods, such as naive MC, IS, conditional Monte Carlo, splitting, etc. Moreover, Assumption A1 allows for the estimators to be constructed from non-i.i.d. data, as in steady-state simulations (Example 10) or randomized quasi-Monte Carlo (RQMC) (Nakayama and Tuffin, 2024; Nakayama et al., 2025), when the estimators obey a joint CLT (12).

To analyze the asymptotic efficiency of $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n})$ in (11), define the gradient (when it exists) of $g(\theta)$ by $\nabla g(\theta) = (g_1(\theta), g_2(\theta), \dots, g_d(\theta))^{\top}$, where $g_i(\theta) = \frac{\partial}{\partial \theta_i} g(\theta)$, $i \in \mathcal{I}$. Assume that $g$ has a nonzero differential at $\theta_{\langle r \rangle}$, which holds, e.g., when each $g_i(\cdot)$, $i \in \mathcal{I}$, exists in a neighborhood of $\theta_{\langle r \rangle}$ and is continuous at $\theta_{\langle r \rangle}$ with $\nabla g(\theta_{\langle r \rangle}) \neq 0$ (Apostol, 1974, Theorem 12.11). If Assumption A1 holds, then the delta method (Serfling, 1980, p. 124) secures the CLT

$$ \sqrt{n} [\hat{\alpha}_{\langle r \rangle, n} - \alpha_{\langle r \rangle}] \Rightarrow \mathcal{N}(0, \tau_{\langle r \rangle}^2) \quad \text{as } n \to \infty, \quad (13) $$

$$ \text{with } \tau_{\langle r \rangle}^2 = [\nabla g(\theta_{\langle r \rangle})]^{\top} \Sigma_{\langle r \rangle} [\nabla g(\theta_{\langle r \rangle})] = \sum_{i=1}^{d} \sum_{j=1}^{d} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} \quad (14) $$

as the asymptotic variance. Thus, we get $\mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}] = \tau_{\langle r \rangle} / |\alpha_{\langle r \rangle}|$ by (2).

We next provide a general setting satisfying Assumption A1; other examples will be covered later. The example considers a vector of means estimated through the average of $n$ i.i.d. observations, which is the sole setting considered in (Nakayama and Tuffin, 2023), and is specialized later to particular estimands $g(\theta_{\langle r \rangle})$ for specific functions $g$ of practical interest in Examples 8, 9, 11, 12, 13, 14, and 15.

**Example 2** ($\theta_{\langle r \rangle}$ is a vector of means). Suppose that

$$ \theta_{\langle r \rangle} = \mathbb{E}[\mathbf{X}_{\langle r \rangle}] \quad \text{for some } \mathbb{R}^d\text{-valued random vector } \mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2}, \dots, X_{\langle r \rangle,d})^{\top}, \quad (15) $$

so $\theta_{\langle r \rangle}$ is a vector of means. Then we may define an estimator of $\theta_{\langle r \rangle}$ as

$$ \hat{\theta}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2}, \dots, \hat{\theta}_{\langle r \rangle, n, d})^{\top} = \frac{1}{n} \sum_{k=1}^{n} \mathbf{X}_{\langle r \rangle, k}, \quad (16) $$

where $\mathbf{X}_{\langle r \rangle, k} = (X_{\langle r \rangle, k, 1}, \dots, X_{\langle r \rangle, k, d})^{\top}$, $k = 1, 2, \dots, n$, are i.i.d. observations of $\mathbf{X}_{\langle r \rangle}$, whose joint distribution may depend on $r$. Let $\Sigma_{\langle r \rangle} = (\Sigma_{\langle r \rangle, i, j} : i, j = 1, 2, \dots, d)$ be the covariance matrix of the random vector $\mathbf{X}_{\langle r \rangle}$, with $\Sigma_{\langle r \rangle, i, j} = \mathrm{Cov}[X_{\langle r \rangle, i}, X_{\langle r \rangle, j}] = \mathbb{E}[(X_{\langle r \rangle, i} - \mathbb{E}[X_{\langle r \rangle, i}))(X_{\langle r \rangle, j} - \mathbb{E}[X_{\langle r \rangle, j}])]$. If $\Sigma_{\langle r \rangle}$ is finite and positive definite, then the CLT (12) holds (Billingsley, 1995, Theorem 29.5).

The different components of $\mathbf{X}_{\langle r \rangle}$ may be generated using different MC methods. For example, $X_{\langle r \rangle, 1}$ may be generated using naive MC, and $X_{\langle r \rangle, 2} = Q_{\langle r \rangle} L_{\langle r \rangle}$ is generated (independently of $X_{\langle r \rangle, 1}$) using IS, as in Example 1. The expectation operator $\mathbb{E}[\cdot]$ in (15) then applies $\mathbb{E}_0[\cdot]$ (resp., $\mathbb{E}_1[\cdot]$) to the first (resp., second) coordinate of $\mathbf{X}$. (Goyal et al., 1992) call this approach measure-specific IS (MSIS). $\square$

The next simple example shows that the estimators $\hat{\theta}_{\langle r \rangle, n}$ and $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n})$ can have quite different asymptotic efficiencies as $r \to \infty$, highlighting the interest of obtaining conditions on the asymptotic efficiency of $\hat{\theta}_{\langle r \rangle, n}$ for getting the same property for $\hat{\alpha}_{\langle r \rangle, n}$.

**Example 3** ($\hat{\theta}_{\langle r \rangle, n}$ and $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n})$ with differing asymptotic behaviors). As in Example 2 with dimension $d = 1$, let $\mathbf{X}_{\langle r \rangle} \sim \mathcal{N}(\theta_{\langle r \rangle}, \sigma_{\langle r \rangle}^2)$, with $\sigma_{\langle r \rangle}^2 \in (0, \infty)$, and we estimate $\theta_{\langle r \rangle}$ via $\hat{\theta}_{\langle r \rangle, n}$ in (16) averaging $n$ i.i.d. copies of $\mathbf{X}_{\langle r \rangle}$. Then $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = \sigma_{\langle r \rangle} / |\theta_{\langle r \rangle}|$ by (2). Now consider the function $g$ with
$\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle}) = \exp(\theta_{\langle r \rangle})$, so the derivative $g'(\theta_{\langle r \rangle}) = g(\theta_{\langle r \rangle})$, and (13) holds with $\tau_{\langle r \rangle}^2 = [g'(\theta_{\langle r \rangle})]^2 \sigma_{\langle r \rangle}^2$ in (14). By (2) we then get

$$ \mathrm{RE}\left[g(\hat{\theta}_{\langle r \rangle, n})\right] = \frac{|g'(\theta_{\langle r \rangle})|\sigma_{\langle r \rangle}}{|g(\theta_{\langle r \rangle})|} = \sigma_{\langle r \rangle}. \qquad (17) $$

For different choices of $\theta_{\langle r \rangle}$ and $\sigma_{\langle r \rangle}^2$, we next compare the asymptotic behaviors of the estimators $\hat{\theta}_{\langle r \rangle, n}$ and $g(\hat{\theta}_{\langle r \rangle, n})$, finding that their efficiencies can vastly differ.

* **Case 1:** Suppose that $\theta_{\langle r \rangle} = r$ and $\sigma_{\langle r \rangle} = \sqrt{\theta_{\langle r \rangle}}$. Then $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = \sigma_{\langle r \rangle}/|\theta_{\langle r \rangle}| = 1/\sqrt{r} \to 0$ as $r \to \infty$, so $\hat{\theta}_{\langle r \rangle, n}$ has VRE by (5). But (17) implies that $\mathrm{RE}[g(\hat{\theta}_{\langle r \rangle, n})] = \sigma_{\langle r \rangle} = \sqrt{r} \to \infty$ as $r \to \infty$, so $g(\hat{\theta}_{\langle r \rangle, n})$ has unbounded (and polynomially increasing) RE. So even if $\hat{\theta}_{\langle r \rangle, n}$ is asymptotically efficient, the resulting $\hat{\alpha}_{\langle r \rangle, n}$ may be inefficient.

* **Case 2:** Now suppose that $\theta_{\langle r \rangle} = 1/r$ and $\sigma_{\langle r \rangle} = \sqrt{\theta_{\langle r \rangle}}$. Then $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = \sigma_{\langle r \rangle}/|\theta_{\langle r \rangle}| = \sqrt{r} \to \infty$ as $r \to \infty$, so $\hat{\theta}_{\langle r \rangle, n}$ has unbounded (polynomially increasing) RE. But $\mathrm{RE}[g(\hat{\theta}_{\langle r \rangle, n})] = \sigma_{\langle r \rangle} = 1/\sqrt{r} \to 0$ as $r \to \infty$, so $g(\hat{\theta}_{\langle r \rangle, n})$ has VRE. Compared to Case 1, we are now in the opposite situation: $\hat{\alpha}_{\langle r \rangle, n}$ is asymptotically efficient even though $\hat{\theta}_{\langle r \rangle, n}$ is not.

* **Case 3:** Suppose that $\theta_{\langle r \rangle} = c_0 e^{-c_1 r}$ and $\sigma_{\langle r \rangle} = c_3 r^{c_2} e^{-c_1 r}$ for some positive constants $c_0, c_1, c_2$, and $c_3$. Then $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n}] = \sigma_{\langle r \rangle}/|\theta_{\langle r \rangle}| = [c_3 r^{c_2} e^{-c_1 r}]/c_0 e^{-c_1 r} = [c_3/c_0]r^{c_2} \to \infty$ as $r \to \infty$, so on a linear scaling, $\hat{\theta}_{\langle r \rangle, n}$ has unbounded RE. However, for each $\delta > 0$, $RE_{1-\delta}[\hat{\theta}_{\langle r \rangle, n}] = [c_3 r^{c_2} e^{-c_1 r}]/[c_0 e^{-c_1 r}]^{1-\delta} = [c_3/c_0^{1-\delta}]r^{c_2} e^{-\delta c_1 r} \to 0$ as $r \to \infty$, so $\hat{\theta}_{\langle r \rangle, n}$ is LE by (7). But $\mathrm{RE}[g(\hat{\theta}_{\langle r \rangle, n})] = \sigma_{\langle r \rangle} = c_3 r^{c_2} e^{-c_1 r} \to 0$ as $r \to \infty$, so $g(\hat{\theta}_{\langle r \rangle, n})$ has VRE. In other words, asymptotic efficiency properties of $\hat{\theta}_{\langle r \rangle, n}$ and $\hat{\alpha}_{\langle r \rangle, n}$ can be different even when considering logarithmic behavior.

□

Example 3 shows that the asymptotic efficiency $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n})$ is not determined by the efficiency of each $\hat{\theta}_{\langle r \rangle, n, i}$, so we now develop conditions that ensure the efficiency of $\hat{\alpha}_{\langle r \rangle, n}$. Specifically, we want to study when $\alpha_{\langle r \rangle}$ can be estimated with BRE as in (4), with VRE as in (5), or with LE as in (7) or (8). To do this, for each $\psi, \zeta > 0$, define

$$ h_{i, \psi, \zeta}(\theta_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle, i}^{\psi} g_i(\theta_{\langle r \rangle})}{g^{\zeta}(\theta_{\langle r \rangle})} \right|. \qquad (18) $$

For analyzing LE in (7) and (8) (but unneeded for BRE or VRE), we further impose the following.

**Assumption A2.** As $r \to \infty$, $\alpha_{\langle r \rangle} \to \chi_0 \in \{0, \infty, -\infty\}$ and $\theta_{\langle r \rangle, i} \to \chi_i \in [-\infty, \infty] \forall i \in \mathcal{I}$.

Under Assumption A2 for establishing LE, the estimand $\alpha_{\langle r \rangle}$ is not permitted to converge to a nonzero finite constant, as noted before in the discussion after (8). But for $i \in \mathcal{I}$, we allow some of the $\theta_{\langle r \rangle, i}$ to converge to nonzero finite constants $\chi_i$ so we can handle situations such as estimating the expected time for the queue-length process of a stable GI/GI/1 queue to first hit a high level $r$; this will be explained later in Examples 8 and 15. For each $\delta > 0$ and $i = 0, 1, \dots, d$, define $\delta_i$ based on Assumption A2 as follows:

$$ \delta_i = \begin{cases} -\delta & \text{if } \chi_i = 0, \\ 0 & \text{if } \chi_i \in (-\infty, \infty) \text{ with } \chi_i \neq 0, \\ \delta & \text{if } |\chi_i| = \infty. \end{cases} \qquad (19) $$

**Theorem 1.** Assume that for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{I}$, and the function $g$ in (9) has a nonzero differential at $\theta_{\langle r \rangle}$. Then under Assumption A1, the following hold as $r \to \infty$:
(i) $\hat{\alpha}_{\langle r \rangle, n}$ has BRE if

$$h_{i,1,1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] = O(1) \text{ as } r \to \infty \text{ for all } i \in \mathcal{I}. \quad (20)$$

(ii) $\hat{\alpha}_{\langle r \rangle, n}$ has VRE if

$$h_{i,1,1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] = o(1) \text{ as } r \to \infty \text{ for all } i \in \mathcal{I}. \quad (21)$$

(iii) When additionally Assumption A2 holds, $\hat{\alpha}_{\langle r \rangle, n}$ is LE if

$$\text{for all } i \in \mathcal{I} \text{ and all } \delta > 0, h_{i,1+\delta_i,1+\delta_0}(\theta_{\langle r \rangle}) \mathrm{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle, n, i}] = o(1), \quad (22)$$

with $\delta_i$ in (19).

(Nakayama and Tuffin, 2023) considers only the special case of Theorem 1 for the setting of Example 2, where $\theta_{\langle r \rangle}$ is a vector of means estimated through a sample average of $n$ i.i.d. random vectors with mean $\theta_{\langle r \rangle}$. Our more general setting allows for $\theta_{\langle r \rangle}$ to consist of non-mean parameters whose estimators satisfy Assumption A1, and conditions (20), (21), and (22) clearly show that even if $\hat{\theta}_{\langle r \rangle, n, i}$ has BRE, VRE, or LE for each $i \in \mathcal{I}$, this does not guarantee that $g(\hat{\theta}_{\langle r \rangle, n})$ is similarly efficient, as seen in Example 3. The factors $h_{i,\psi,\zeta}(\theta_{\langle r \rangle})$ from (18), which capture certain joint behavior of all indices in $\mathcal{I}$, could overwhelm the RE[$\hat{\theta}_{\langle r \rangle, n, i}$] or RE$_{1+\delta_i}$[$\hat{\theta}_{\langle r \rangle, n, i}$].

**Proof of Theorem 1.** Under our condition that $g$ in (9) has a nonzero differential at $\theta_{\langle r \rangle}$ and Assumption A1, the delta method (Serfling, 1980, p. 124) yields the CLT (13) with asymptotic variance $\tau_{\langle r \rangle}^2$ in (14). By the Cramer-Wold device (e.g., (Serfling, 1980, Theorem 1.5.2)), the multivariate CLT (12) implies that each individual component $\hat{\theta}_{\langle r \rangle, n, i}$ of the vector $\hat{\theta}_{\langle r \rangle, n}$ in (10) satisfies a univariate CLT $\sqrt{n}[\hat{\theta}_{\langle r \rangle, n, i} - \theta_{\langle r \rangle, i}] \Rightarrow N_1(0, \Sigma_{\langle r \rangle, i, i})$ as $n \to \infty$, where $\Sigma_{\langle r \rangle, i, i}$ is the asymptotic variance in the univariate CLT. By Assumption A1, $\Sigma_{\langle r \rangle}$ is a covariance matrix of a multivariate normal $Z = (Z_1, Z_2, \dots, Z_d)^\top$, so each $\Sigma_{\langle r \rangle, i, i}$ is a variance (of the normal $Z_i$), and $\Sigma_{\langle r \rangle, i, j}$ is a covariance (of the two normals $Z_i$ and $Z_j$). Thus, the Cauchy-Schwarz inequality implies that $|\Sigma_{\langle r \rangle, i, j}| \le \sqrt{\Sigma_{\langle r \rangle, i, i} \Sigma_{\langle r \rangle, j, j}}$, and also we have that RE[$\hat{\theta}_{\langle r \rangle, n, i}$] = $\sqrt{\Sigma_{\langle r \rangle, i, i}}/|\theta_{\langle r \rangle, i}|$, as in (2), because $\theta_{\langle r \rangle, i} \ne 0$. Thus, the RE of $\hat{\alpha}_{\langle r \rangle, n}$ then satisfies

$$\begin{align} (\mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}])^2 &= \frac{\tau_{\langle r \rangle}^2}{\alpha_{\langle r \rangle}^2} = \sum_{i=1}^d \sum_{j=1}^d \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \Sigma_{\langle r \rangle, i, j} \le \sum_{i=1}^d \sum_{j=1}^d \left| \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \Sigma_{\langle r \rangle, i, j} \right| \nonumber \\ &\le \sum_{i=1}^d \sum_{j=1}^d \left| \frac{g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \right| \sqrt{\sum_{\langle r \rangle, i, i} \sum_{\langle r \rangle, j, j}} \tag{23} \\ &= \sum_{i=1}^d \sum_{j=1}^d \left| \frac{\theta_{\langle r \rangle, i} g_i(\theta_{\langle r \rangle}) \theta_{\langle r \rangle, j} g_j(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \right| \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, j}] \nonumber \\ &= \left( \sum_{i=1}^d h_{i,1,1}(\theta_{\langle r \rangle}) \mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, i}] \right)^2 \end{align}$$

by (18) since $\theta_{\langle r \rangle, i} \ne 0$ and $\theta_{\langle r \rangle, j} \ne 0$, so parts (i) and (ii) hold by (20) and (21), respectively.

For part (iii) under Assumption A2, we simplify the discussion by considering the square of (6) for
$\zeta = 1 + \delta_0$. As in (23), we get

$$
\begin{align*}
(\mathrm{RE}_{1+\delta_0}[\hat{\alpha}_{\langle r \rangle, n}])^2 &= \frac{\tau_{\langle r \rangle}^2}{|\alpha_{\langle r \rangle}|^{2+2\delta_0}} \le \sum_{i=1}^d \sum_{j=1}^d \left| \frac{g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})}{g^{2+2\delta_0}(\theta_{\langle r \rangle})} \right| \sqrt{\sum_{\langle r \rangle, i, n} \sum_{\langle r \rangle, j, n} \tau_{\langle r \rangle}^2} \\
&= \sum_{i=1}^d \sum_{j=1}^d \left| \frac{\theta_{\langle r \rangle, i}^{1+\delta_i} g_i(\theta_{\langle r \rangle}) \theta_{\langle r \rangle, j}^{1+\delta_j} g_j(\theta_{\langle r \rangle})}{g^{1+\delta_0}(\theta_{\langle r \rangle}) g^{1+\delta_0}(\theta_{\langle r \rangle})} \right| \mathrm{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle, n, i}] \mathrm{RE}_{1+\delta_j}[\hat{\theta}_{\langle r \rangle, n, j}] \\
&= \left( \sum_{i=1}^d h_{i,1+\delta_i,1+\delta_0}(\theta_{\langle r \rangle}) \mathrm{RE}_{1+\delta_i}[\hat{\theta}_{\langle r \rangle, n, i}] \right)^2
\end{align*}
$$

by (18) since $\theta_{\langle r \rangle, i} \neq 0$ and $\theta_{\langle r \rangle, j} \neq 0$, so part (iii) holds by (7) and (8) under (22). $\square$

<figure><img src="image_3.png" /></figure>

4 Examples Related to Assumption A1 and Theorem 1

The following example generalizes MSIS from Example 2 to applying other MC methods to estimate different
components of $\theta_{\langle r \rangle}$ and with different sample sizes in Assumption A1.

**Example 4 (measure-specific sampling (MSS)).** Let $\theta = \theta_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2}, \dots, \theta_{\langle r \rangle,d})^\top$ be a vector of unknown parameters, not necessarily means, for fixed $d \ge 1$. For each $i \in \mathcal{I}$, we estimate $\theta_{\langle r \rangle, i}$ via some MC method $\mathfrak{M}_i$. Define positive sampling weights $\varrho = (\varrho_i : i \in \mathcal{I})$ with $\sum_{i \in \mathcal{I}} \varrho_i = 1$, and for a total sampling budget $n$, the estimator $\hat{\theta}_{\langle r \rangle, n, i}$ of $\theta_{\langle r \rangle, i}$ applies method $\mathfrak{M}_i$ with a sample size $n_i \equiv [\varrho_i n]$, with $\cdot$ the floor function. This generalizes MSIS by allowing MC methods other than naive MC and IS, and also permits different sample sizes for the various methods. Example 2 and (Nakayama and Tuffin, 2023) restrict to the setting of using the same sample size $n$ for each estimator $\hat{\theta}_{\langle r \rangle, n, i}$, $i \in \mathcal{I}$. $\square$

<figure><img src="image_4.png" /></figure>

Assumption A1 also allows for components of $\theta_{\langle r \rangle}$ to be non-mean parameters, such as a quantile, and
we now further consider some specific functions $g$ to get estimands as in (9).

**Example 5** ($\theta_{\langle r \rangle}$ is a vector of quantiles). Suppose that $F_{\langle r \rangle}$ is a cumulative distribution function (CDF) of some random variable $Y_{\langle r \rangle}$, which we denote by $Y_{\langle r \rangle} \sim F_{\langle r \rangle}$. For $0 < p_1 < p_2 < \dots < p_d < 1$, where each $p_j = p_{j,\langle r \rangle}$ may also depend on $r$, let $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2}, \dots, \theta_{\langle r \rangle,d})^\top = (F_{\langle r \rangle}^{-}(p_1), F_{\langle r \rangle}^{-}(p_2), \dots, F_{\langle r \rangle}^{-}(p_d))^\top$, where $F^{-}(p) = \inf\{y : F(y) \ge p\}$ is the p-quantile of a CDF $F$ for $p \in (0, 1)$. Given a sample $Y_{\langle r \rangle,1}, Y_{\langle r \rangle,2}, \dots, Y_{\langle r \rangle,n}$ of $n$ i.i.d. copies of $Y_{\langle r \rangle}$, we estimate each $\theta_{\langle r \rangle, i} = F_{\langle r \rangle}^{-}(p_i)$ as follows. Based on the sample, let $\hat{F}_{\langle r \rangle, n}$ be the empirical distribution, defined by $\hat{F}_{\langle r \rangle, n}(y) = (1/n) \sum_{k=1}^n I(Y_{\langle r \rangle, k} \le y)$ for each $y \in \mathbb{R}$. Then we may estimate $\theta_{\langle r \rangle, i}$ via $\hat{\theta}_{\langle r \rangle, n, i} = \hat{F}_{\langle r \rangle, n}^{-}(p_i)$. Assume that $F_{\langle r \rangle}$ has a density $f_{\langle r \rangle}$ in neighborhoods of $F_{\langle r \rangle}^{-}(p_i)$, $i \in \mathcal{I}$, and that $f_{\langle r \rangle}$ is positive and continuous at $F_{\langle r \rangle}^{-}(p_i)$, $i \in \mathcal{I}$. The vector estimator $\hat{\theta}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n,1}, \hat{\theta}_{\langle r \rangle, n,2}, \dots, \hat{\theta}_{\langle r \rangle, n,d})^\top$ then obeys the joint CLT (12), with $\sum_{\langle r \rangle, i, j} = p_i(1-p_j)/[f_{\langle r \rangle}(F_{\langle r \rangle}^{-}(p_i))f_{\langle r \rangle}(F_{\langle r \rangle}^{-}(p_j))]$ for $i \le j$, and $\sum_{\langle r \rangle, i, j} = \sum_{\langle r \rangle, j, i}$ for $i > j$, where we assume that $\sum_{\langle r \rangle}$ is positive definite so that Assumption A1 holds; see (Serfling, 1980, Theorem 2.3.3B, p. 80). One example of $\alpha_{\langle r \rangle}$ from (9) with $d = 2$ in this setting is the interquartile range (Serfling, 1980, Section 2.3.6), given by $\alpha_{\langle r \rangle} = g(\theta) = \theta_{\langle r \rangle,1} - \theta_{\langle r \rangle,2}$ with $\theta_{\langle r \rangle,1} = F_{\langle r \rangle}^{-}(0.75)$ and $\theta_{\langle r \rangle,2} = F_{\langle r \rangle}^{-}(0.25)$; this is sometimes used as an alternative to the standard deviation as a measure of a distribution's dispersion. Another example that arises in studies of economic inequality is the decile ratio (De Maio, 2007), where $g(\theta) = \theta_{\langle r \rangle,1}/\theta_{\langle r \rangle,2}$ with $\theta_{\langle r \rangle,1} = F_{\langle r \rangle}^{-}(0.9)$ and $\theta_{\langle r \rangle,2} = F_{\langle r \rangle}^{-}(0.1)$. $\square$

<figure><img src="image_5.png" /></figure>

Assumption A1 also permits the vector $\theta_{\langle r \rangle}$ to contain different types of parameters, such as a combination
of a quantile and mean.

**Example 6** ($\theta_{\langle r \rangle}$ is a vector of a quantile and a mean). For $d = 2$, suppose that $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2})^\top = (F_{\langle r \rangle}^{-}(p), \mu_{\langle r \rangle})^\top$, where $\mu_{\langle r \rangle} = \mathbb{E}[Y_{\langle r \rangle}] = \int y dF_{\langle r \rangle}(y)$ with $Y_{\langle r \rangle} \sim F_{\langle r \rangle}$ and $p \in (0,1)$. Given a sample
$Y_{\langle r \rangle, 1}, Y_{\langle r \rangle, 2}, \dots, Y_{\langle r \rangle, n}$ of n i.i.d. copies of $Y_{\langle r \rangle}$, we can estimate $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2})^\top$ via

$$ \hat{\theta}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2})^\top = (\hat{F}_{\langle r \rangle, n}^{-}(p), (1/n) \sum_{k=1}^{n} Y_{\langle r \rangle, k})^\top, $$

where $\hat{F}_{\langle r \rangle, n}$ is the empirical CDF as in Example 5. Assuming that $Y_{\langle r \rangle}$ has variance $\text{Var}[Y_{\langle r \rangle}] \in (0, \infty)$ and $F_{\langle r \rangle}$ has a density $f_{\langle r \rangle}$ that is continuous and positive at $F_{\langle r \rangle}^{-}(p)$, (Ferguson, 1999) establishes the joint CLT (12), with $\Sigma_{\langle r \rangle, 1, 1} = p(1-p)/f_{\langle r \rangle}^2(F_{\langle r \rangle}^{-}(p))$, $\Sigma_{\langle r \rangle, 2, 2} = \text{Var}[Y_{\langle r \rangle}]$, and $\Sigma_{\langle r \rangle, 1, 2} = \Sigma_{\langle r \rangle, 2, 1} = (p\mu_{\langle r \rangle} - \mathbb{E}[Y_{\langle r \rangle}I(Y_{\langle r \rangle} < F_{\langle r \rangle}^{-}(p))]/f_{\langle r \rangle}(F_{\langle r \rangle}^{-}(p));$ also, see (Lin et al., 1980). For this example, defining the function $g$ in (9) as $g(\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2}) = \theta_{\langle r \rangle, 1} - \theta_{\langle r \rangle, 2}$ leads to $\alpha_{\langle r \rangle}$ (typically with $p \approx 1$) being the economic capital (EC), which is a risk measure used in finance; e.g., see (Klaassen and van Eeghen, 2009, p. 5), (Lütkebohmert, 2009, Section 2.4), (Scandizzo, 2016, p. 194), and (Li et al., 2024). For $Y_{\langle r \rangle}$ as an i.i.d. sum of m light-tailed random variables, Example 7 below describes asymptotic properties ($a.s. \ r \to \infty$) of EC estimators established in (Li et al., 2024). $\square$

**Example 7 (Economic capital).** We now specialize Example 6 to $Y_{\langle r \rangle} = \sum_{j=1}^m Z_{\langle r \rangle, j}$ as the sum of $r = m$ i.i.d. light-tailed random variables, where each i.i.d. summand $Z_{\langle r \rangle, j}$ has a marginal distribution $G_0$, not depending on $r = m$. As in (Glynn, 1996), let $p = p_{\langle r \rangle} = 1 - e^{-\beta m}$ be the quantile level, for a constant $\beta > 0$. (Li et al., 2024) analyzes the RE ($a.s. \ m \to \infty$) of different combinations of naive MC and IS to estimate the EC $\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2}) = \theta_{\langle r \rangle, 1} - \theta_{\langle r \rangle, 2}$ for $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2})^\top = (F_{\langle r \rangle}^{-}(p), \mu_{\langle r \rangle})^\top$, where $\mu_{\langle r \rangle} = \mathbb{E}[Y_{\langle r \rangle}] = \int y dF_{\langle r \rangle}(y)$ with $Y_{\langle r \rangle} \sim F_{\langle r \rangle}$. For $\vartheta \in \mathbb{R}$, let $M_0(\vartheta) = \int e^{\vartheta z} dG_0(z)$ be the moment generating function of $G_0$, and let $Q_0(\vartheta) = \ln M_0(\vartheta)$ be the corresponding cumulant generating function. Define $\Delta = \{\vartheta : M_0(\vartheta) < \infty\}$ as the domain of $M_0(\cdot)$, and assume that its interior $\Delta^\circ$ contains the origin (so $G_0$ is light-tailed). For IS, we use an exponentially twisted distribution $\tilde{G}_{0, \vartheta}$ with $d\tilde{G}_{0, \vartheta}(z) = e^{\vartheta z - Q_0(\vartheta)} dG_0(z)$, $z \in \mathbb{R}$, for $\vartheta \in \Delta^\circ$, so IS samples $Z_{\langle r \rangle, 1}, Z_{\langle r \rangle, 2}, \dots, Z_{\langle r \rangle, m}$ as i.i.d. from $\tilde{G}_{0, \vartheta}$. Then when estimating just the mean $\mu_{\langle r \rangle}$, naive MC leads to an estimator of $\mu_{\langle r \rangle}$ with VRE as $m \to \infty$, whereas the IS estimator with any $\vartheta \neq 0$ has exponentially increasing RE as $m \to \infty$ (Li et al., 2024). When estimating just the $p_m$-quantile $F_{\langle r \rangle}^{-}(p_m)$, the naive MC estimator of $F_{\langle r \rangle}^{-}(p_m)$ has exponentially increasing RE as $m \to \infty$, and when applying IS, assume that there exists $\vartheta = \vartheta_* \in \Delta^\circ$ with $\vartheta_* > 0$ such that $-\vartheta_* Q'_0(\vartheta_*) + Q_0(\vartheta_*) = -\beta$ (Glynn, 1996). With twisting parameter $\vartheta_*$, the resulting IS estimator of $F_{\langle r \rangle}^{-}(p_m)$ has VRE as $m \to \infty$ when the characteristic function $C_0(\vartheta) = \int e^{\vartheta \sqrt{-1} z} dG_0(z)$, $\vartheta \in \mathbb{R}$, satisfies $\int_{\mathbb{R}} |C_0(\vartheta)|^{q_0} d\vartheta < \infty$ for some $q_0 \ge 1$. In (18) and Theorem 1, $h_{i, 1, 1}(\theta_{\langle r \rangle}) = 1$ for $i = 1, 2$. Now consider applying MSIS, with IS using an exponential twist with twisting parameter $\vartheta = \vartheta_*$ to estimate $F_{\langle r \rangle}^{-}(p_m)$, and independently employing naive MC to estimate $\mu_{\langle r \rangle}$. The resulting MSIS estimator of the EC $\alpha_{\langle r \rangle}$ has VRE as $m \to \infty$ by Theorem 1(ii), which is also established in (Li et al., 2024, Theorem 5(iii)). While our Theorem 1 establishes only sufficient (but not necessary) conditions to ensure an efficient estimator of $\alpha_{\langle r \rangle}$, (Li et al., 2024, Theorem 5(i) and (ii)) further prove that using naive MC only or IS only to estimate $\alpha_{\langle r \rangle}$ leads to unbounded RE growing exponentially in $m \to \infty$, so LE does not hold. That LE does not hold in these cases can also be established using Theorem 3(iii), which appears later in Section 5.2. $\square$

Other possibilities that can fit into the framework of (9) and Assumption A1 are distortion risk measures and their sensitivities (Glynn et al., 2021); estimated functions at estimated parameters (Glynn et al., 2020), which includes as a special case the CVaR = $\mathbb{E}[Y | Y > F^{-}(p)]$ for $Y \sim F$; etc.

Examples 8, 9, 11, 12, 13, 14, and 15 below, which are also considered in (Nakayama and Tuffin, 2023), provide other instances of estimands fitting the framework of (9) and Assumption A1.

**Example 8 (Expected hitting time).** Let $Y = (Y(t) : t \ge 0)$ be a nondelayed regenerative process (Kalashnikov, 1994) living on a state space $S \subset \mathbb{R}^v$ for some $v \ge 1$. Thus, there is an infinite sequence of times $0 = T_0 < T_1 < \dots$ at which $Y$ “probabilistically restarts”, and we split the sample path of $Y$ into i.i.d. cycles defined by successive regeneration times. The goal is to estimate $\alpha = \mathbb{E}[T_A]$, where $T_A = \inf\{t > 0 : Y(t) \in A\}$ is the hitting time to some set $A \subset S$ of states. For example, $Y$ may be the queue-length (including any
customer in service) process of a stable GI/GI/1 queue, where customer arrivals to an empty queue are regenerations, and the first customer arrives at time 0 to an empty system. Also, let $A = A_{\langle r \rangle} = \{r, r+1, \dots\}$ denote a set of large queue lengths (e.g., buffer overflows), so $\alpha = \alpha_{\langle r \rangle}$ is the expected time the queue length first reaches a large level $r$. Alternatively, $Y$ may be an irreducible, positive-recurrent continuous-time Markov chain (CTMC), so successive returns to a fixed state $s_0 \in S$ are regenerations, where $Y(0) = s_0$. As will be discussed later in Example 11, when $Y$ represents a Markovian reliability system, $A$ may be the set of “failed states”, where we assume that $s_0 \notin A$, so $\alpha$ denotes the mean time to failure (MTTF). Applying well-known regenerative theory (e.g., (Goyal et al., 1992) and (Glynn et al., 2017)) enables us to express $\alpha$ as a ratio $\alpha = \eta/\gamma$, with $\eta = \mathbb{E}[\min(T_A, T_1)]$ and $\gamma = \mathbb{E}[I(T_A < T_1)] = \mathbb{P}(T_A < T_1)$, where both $\min(T_A, T_1)$ and $I(T_A < T_1)$ are measurable with respect to the sigma-field of the process $Y$ up to time $T_1$. This problem setting fits into the framework of Example 2 with $d=2$, where the random vector $\mathbf{X} = (X_1, X_2)^\top$ in (15) has $X_1 = \min(T_A, T_1)$ and $X_2 = I(T_A < T_1)$. Thus, we obtain $\alpha = g(\boldsymbol{\theta})$ in (9) with $\boldsymbol{\theta} = (\eta, \gamma)^\top$, and $g(x_1, x_2) = x_1/x_2$. Also, the estimator $\hat{\boldsymbol{\theta}}_{\langle r \rangle, n}$ in (16) averages i.i.d. copies of $\mathbf{X}$ from simulating $Y$ for $n$ cycles. When $Y$ is a CTMC, we can obtain reduced-variance estimators of $\eta, \gamma$, and $\alpha$ via conditional MC (CMC; e.g., Section V.4 of (Asmussen and Glynn, 2007)) by conditioning on the CTMC’s embedded discrete-time Markov chain (DTMC). A simulation applying CMC then generates only the embedded DTMC, with the exponential holding times in visited states replaced by their conditional means. $\square$

**Example 9 (Conditional expectation).** For a random vector $(Z_1, Z_2)^\top \in \mathbb{R}^2$ and a set $A \subset \mathbb{R}$, consider the conditional expectation $\alpha = \mathbb{E}[Z_1 | Z_2 \in A] = \mathbb{E}[Z_1 I(Z_2 \in A)]/\mathbb{E}[I(Z_2 \in A)]$. Then $\alpha$ is in the form of (9) with $g(\theta_1, \theta_2) = \theta_1/\theta_2$, $\theta = (\mathbb{E}[X_1], \mathbb{E}[X_2])^\top$ for $\mathbf{X} = (X_1, X_2)^\top$ in (15), $X_1 = Z_1 I(Z_2 \in A)$, and $X_2 = I(Z_2 \in A)$. $\square$

For $d=2$ in Example 2, using $g(x_1, x_2) = x_1/x_2$ covers the case of ratio estimators, such as those arising when applying the regenerative method (Glynn and Iglehart, 1993) with $n$ the deterministic number of i.i.d. cycles simulated. This includes estimating the mean hitting time in Example 8. But the next example shows that our framework also can handle situations when the computing budget $n$ represents, e.g., the simulated time or the number of transitions. In these cases, the number of completed cycles is then typically random, with dependence in the sample paths of the regenerative process.

**Example 10 (Regenerative process on the simulated time scale).** As in Example 8, let $Y = [Y(t) : t \ge 0]$ be a (nondelayed) continuous-time regenerative process on a state space $S$, with regeneration epochs $0 = T_0 < T_1 < T_2 < \dots$. We allow for $Y = Y_{\langle r \rangle}$. For a reward function $f: S \to \mathbb{R}$, we want to estimate the long-run time-average reward $\alpha_{\langle r \rangle} = \lim_{t \to \infty} (1/t) \int_0^t f(Y(s)) \, ds$. For each $k \ge 1$, let $V_k = \int_{T_{k-1}}^{T_k} f(Y(s)) \, ds$ be the cumulated reward over the $k$th cycle. Let $N(t) = \sup\{k \ge 1 : T_k \le t\}$ be the number of cycles completed by time $t$. Letting $n=t$, define $\hat{\alpha}_{\langle r \rangle, n} = \sum_{k=1}^{N(t)} V_k/T_{N(t)}$ as the ratio estimator based on the sample path up to time $t$ (since $T_{N(t)} = \sum_{k=1}^{N(t)} (T_k - T_{k-1})$). Assume that $\mathbb{E}[T_1] < \infty$ and $\mathbb{E}[(V_1 - \alpha_{\langle r \rangle} T_1)^2] \in (0, \infty)$. For $\theta_{\langle r \rangle, 1} = \mathbb{E}[V_1]$ and $\theta_{\langle r \rangle, 2} = \mathbb{E}[T_1]$, we have that (Glynn and Iglehart, 1993) $\alpha_{\langle r \rangle} = \theta_{\langle r \rangle, 1}/\theta_{\langle r \rangle, 2}$, and Assumption A1 holds with $\Sigma_{\langle r \rangle, 1, 1} = \mathbb{E}[T_1] \text{Var}[V_1]$, $\Sigma_{\langle r \rangle, 1, 2} = \mathbb{E}[T_1] \text{Var}[T_1]$, and $\Sigma_{\langle r \rangle, 1, 2} = \mathbb{E}[T_1] \text{Cov}[V_1, T_1]$. In contrast, when $n$ represents the (deterministic) number of cycles simulated, the setting in Example 2 holds, where Assumption A1 is valid but now with $\Sigma_{\langle r \rangle, 1, 1} = \text{Var}[V_1]$, $\Sigma_{\langle r \rangle, 1, 2} = \text{Var}[T_1]$, and $\Sigma_{\langle r \rangle, 1, 2} = \text{Cov}[V_1, T_1]$. $\square$

**Example 11 (Mean time to failure (MTTF)).** A highly reliable Markovian system (HRMS) is a well-known model involving rare events; e.g., see (Goyal et al., 1992) and (Shahabuddin, 1994). An HRMS comprises a fixed number of dependable components subject to random failures and repairs, where component failure and repair times are exponentially distributed. As in Example 8, let $Y = (Y(t) : t \ge 0)$ denote a CTMC of the evolution of the HRMS over time on a finite state space $S$, where each state $s \in S$ keeps track of which components are up and down, as well as any other information needed (e.g., queueing of failed components awaiting repair) to ensure that $Y$ has the Markov property. The system is failed when certain combinations of components are failed, and let $A \subset S$ be the set of failed states. At time 0, the system starts in a state $s_0 \in S \setminus A$ in which all components are operational, and let $\mathbb{E}_0$ denote the expectation operator under
the original system dynamics, which corresponds to naive MC. The HRMS literature (e.g., (Rubino and Tuffin, 2009a)) considers the system as made up of highly reliable components, modeled by parameterizing the failure rates of components as positive powers of $\epsilon = 1/r$, keeping repair rates as fixed constants, and analyzes the behavior of estimators of dependability measures of the HRMS as $r \to \infty$. Denote the CTMC as $Y \equiv Y_{\langle r \rangle}$ under this parameterization, and let $T_{A, \langle r \rangle} = \inf\{t \ge 0 : Y_{\langle r \rangle}(t) \in A\}$ and $T_{1, \langle r \rangle}$ be the first return time to $s_0$. We can then express the MTTF in Example 8 as $\alpha_{\langle r \rangle} = \mathbb{E}_0[T_{A, \langle r \rangle}] = \eta_{\langle r \rangle}/\gamma_{\langle r \rangle}$, with $\eta_{\langle r \rangle} = \mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_{1, \langle r \rangle})]$ and $\gamma_{\langle r \rangle} = \mathbb{E}_0[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})] = \mathbb{E}_1[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}]$, where $\mathbb{E}_1[\cdot]$ represents expectation under IS, with $L_{\langle r \rangle}$ as the resulting likelihood ratio. To put this in the framework of Example 2 with $d=2$, define the random vector $\mathbf{X} = (X_1, X_2)^\top = (\min(T_{A, \langle r \rangle}, T_{1, \langle r \rangle}), I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle})^\top$, where the expectation operator $\mathbb{E}[\cdot]$ applies $\mathbb{E}_0[\cdot]$ (resp., $\mathbb{E}_1[\cdot]$) to the first (resp., second) coordinate of $\mathbf{X}$. Moreover, we write $\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle})$ in (9) with $\theta_{\langle r \rangle} = (\eta_{\langle r \rangle}, \gamma_{\langle r \rangle})^\top$ and $g(x_1, x_2) = x_1/x_2$. When combining CMC (by conditioning on the embedded DTMC; see Example 8) with an IS technique known as “balanced failure biasing” (see Example 20), (Shahabuddin, 1994) establishes the BRE of the MTTF MSIS ratio estimator, applying the IS for estimating only the denominator $\gamma_{\langle r \rangle}$. We can similarly show the BRE of the MSIS ratio estimator through Theorem 1(i), using that both the numerator and denominator in the ratio estimator have BRE, and $h_{1,1,1}(\theta_{\langle r \rangle}) = h_{2,1,1}(\theta_{\langle r \rangle}) = 1$ in (18). $\square$

**Example 12 (Derivative of expected hitting time).** For the expected hitting time $\beta \equiv \mathbb{E}[T_A]$ in Example 8, suppose that we want to compute its sensitivity with respect to a model parameter $\lambda$. Examples of such a $\lambda$ include a component's failure rate in a reliability system, and a parameter in the service distribution in a GI/GI/1 queue. Sensitivity information can be useful for identifying parts of a system to upgrade to obtain the greatest improvement in the overall system performance. As $\beta$ now depends on the value of $\lambda$, we write $\beta(\lambda) = \eta(\lambda)/\gamma(\lambda)$, and our goal is to estimate $\beta'(\lambda) \equiv \frac{d}{d\lambda}\beta(\lambda) = [\eta'(\lambda)\gamma(\lambda)-\eta(\lambda)\gamma'(\lambda)]/\gamma^2(\lambda)$, where $\eta'(\lambda) = \frac{d}{d\lambda}\eta(\lambda)$ and $\gamma'(\lambda) = \frac{d}{d\lambda}\gamma(\lambda)$. Employing the likelihood ratio derivative method (Glynn, 1990; Glynn et al., 1991), we can express the derivatives as expectations $\eta'(\lambda) = \mathbb{E}[\min(T_A, T_1)L']$ and $\gamma'(\lambda) = \mathbb{E}[I(T_A < T_1)L']$, where $L'$ denotes the derivative of the likelihood ratio $L$ with respect to $\lambda$ up to time $\min(T_A, T_1)$; (Nakayama, 1998) further applies CMC (see Example 8) in the estimation for an HRMS. We then obtain an estimator of $\beta'(\lambda)$ fitting into the framework of Example 2 with $d=4$ by taking $\theta = \mathbb{E}[\mathbf{X}] = (\eta(\lambda), \gamma(\lambda), \eta'(\lambda), \gamma'(\lambda))^\top$ with $\mathbf{X} = (X_1, X_2, X_3, X_4)^\top$ for $X_1 = \min(T_A, T_1)$, $X_2 = I(T_A < T_1)$, $X_3 = \min(T_A, T_1)L'$, and $X_4 = I(T_A < T_1)L'$, and $g(x_1, x_2, x_3, x_4) = [x_3x_2 - x_4x_1]/x_2^2$ in (9). $\square$

**Example 13.** For an HRMS as in Example 11, (Nakayama, 1998) specializes Example 12 to estimate the derivative of the MTTF with respect to a component failure rate, by applying a combination of variance-reduction techniques. CMC conditions on the embedded DTMC. Also, MSIS estimates independently the four components in $\theta = (\eta(\lambda), \gamma(\lambda), \eta'(\lambda), \gamma'(\lambda))^\top$ in Example 12, using balanced failure biasing (resp., naive MC) to estimate independently $\eta'(\lambda)$, $\gamma(\lambda)$, and $\gamma'(\lambda)$ (resp., $\eta(\lambda)$). (Nakayama, 1998) establishes BRE for the resulting estimator of the derivative of the MTTF. Example 20 provides numerical results from using either naive MC or only IS to estimate independently all components of $\theta$. $\square$

**Example 14 (Alternative MTTF estimator).** Instead of Example 11's ratio estimator for the HRMS MTTF $\alpha_{\langle r \rangle} = \mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_{1, \langle r \rangle})]/\mathbb{E}_0[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})]$, (Nakayama and Tuffin, 2019) develops an alternative approach that is more efficient in certain contexts and also applies to Example 8's broader setting. Write

$$ \alpha_{\langle r \rangle} = \frac{\mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_{1, \langle r \rangle})I(T_{A, \langle r \rangle} > T_{1, \langle r \rangle})] + \mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_{1, \langle r \rangle})I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})]}{\mathbb{E}_0[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})]} \\ = \frac{\mathbb{E}_0[T_{1, \langle r \rangle}I(T_{A, \langle r \rangle} > T_{1, \langle r \rangle})] + \mathbb{E}_1[T_{A, \langle r \rangle}I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}]}{\mathbb{E}_1[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}]}, $$

which fits (9) for function $g(\theta_{\langle r \rangle}) = (\theta_{\langle r \rangle,1}+\theta_{\langle r \rangle,2})/\theta_{\langle r \rangle,3}$ and $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2}, \theta_{\langle r \rangle,3})^\top$ of dimension $d=3$ with $\theta_{\langle r \rangle,1} = \mathbb{E}_0[T_{1, \langle r \rangle}I(T_{A, \langle r \rangle} > T_{1, \langle r \rangle})]$, $\theta_{\langle r \rangle,2} = \mathbb{E}_1[T_{A, \langle r \rangle}I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}]$, and $\theta_{\langle r \rangle,3} = \mathbb{E}_1[I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}]$. We estimate $\theta_{\langle r \rangle}$ using MSIS as in Example 2 by generating i.i.d. copies of the random vector $\mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2}, X_{\langle r \rangle,3})^\top$ with $X_{\langle r \rangle,1} = T_{1, \langle r \rangle}I(T_{A, \langle r \rangle} > T_{1, \langle r \rangle})$, $X_{\langle r \rangle,2} = T_{A, \langle r \rangle}I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle})L_{\langle r \rangle}$,
and $X_{\langle r \rangle, 3} = I(T_{A, \langle r \rangle} < T_{1, \langle r \rangle}) L_{\langle r \rangle}$, where $X_{\langle r \rangle, 1}$ (resp., $(X_{\langle r \rangle, 2}, X_{\langle r \rangle, 3})$) is generated with naive MC (resp., IS). Hence, $\theta_{\langle r \rangle} = \mathbb{E}[\mathbf{X}_{\langle r \rangle}]$ as in (15), where the expectation operator $\mathbb{E}[\cdot]$ applies $\mathbb{E}_0[\cdot]$ (resp., $\mathbb{E}_1[\cdot]$) to the first (resp., second and third) coordinate(s) of $\mathbf{X}_{\langle r \rangle}$. Note that (18) has $h_{1,1,1}(\theta_{\langle r \rangle}) = \frac{\theta_{\langle r \rangle, 1}}{\theta_{\langle r \rangle, 1} + \theta_{\langle r \rangle, 2}} \le 1$, $h_{2,1,1}(\theta_{\langle r \rangle}) = \frac{\theta_{\langle r \rangle, 2}}{\theta_{\langle r \rangle, 1} + \theta_{\langle r \rangle, 2}} \le 1$ and $h_{3,1,1}(\theta_{\langle r \rangle}) = 1$, which are all bounded, so by Theorem 1(i), the estimator of $\alpha_{\langle r \rangle}$ has BRE if the estimators of each of the three components of $\theta_{\langle r \rangle}$ have the same property. We can instead show that the estimator of $\alpha_{\langle r \rangle}$ is LE (which is more appropriate in, e.g., the GI/GI/1 setting) by applying an analysis developed in Example 15. $\square$

**Example 15** (Expected time for GI/GI/1 queue length to reach a large threshold). Consider again the queue-length process of a stable GI/GI/1 queue from Example 8. For $\mathbb{E}_0$ and $\mathbb{P}_0$ denoting the expectation operator and probability measure under the original system dynamics, we want to estimate $\alpha_{\langle r \rangle} = \mathbb{E}_0[T_{A, \langle r \rangle}]$, which is the expected time that the queue length first enters the set $A_{\langle r \rangle} = \{r, r+1, \dots\}$. Denote a generic interarrival (resp., service) time by $U$ (resp., $V$), with $U$ and $V$ independent and having distributions not depending on $r$ and light-tailed (in the sense that their moment-generating functions exist in a neighborhood of 0, ensuring $U$ and $V$ have finite moments of all orders). Assume that $0 < \mathbb{E}_0[V] < \mathbb{E}_0[U]$ and $\mathbb{P}_0(V > U) > 0$. We then can express $\alpha_{\langle r \rangle}$ as a ratio $\alpha_{\langle r \rangle} = \eta_{\langle r \rangle}/\gamma_{\langle r \rangle} \in (0, \infty)$, with $\eta_{\langle r \rangle} = \mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_1)]$ and $\gamma_{\langle r \rangle} = \mathbb{E}_0[I(T_{A, \langle r \rangle} < T_1)] = \mathbb{E}_1[I(T_{A, \langle r \rangle} < T_1)L]$, where $T_1$ is the first regeneration epoch after time 0, and $\mathbb{E}_1[\cdot]$ denotes expectation under IS with $L$ as the resulting likelihood ratio. Under $\mathbb{P}_0$, the distribution of $T_{A, \langle r \rangle}$ depends on $r$, but the distribution of $T_1$ does not. In the framework of Example 2, we take $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2})^\top = (\eta_{\langle r \rangle}, \gamma_{\langle r \rangle})^\top$ having dimension $d=2$, and $g(x_1, x_1) = x_1/x_2$ in (9). We use MSIS to estimate $\theta_{\langle r \rangle} = \mathbb{E}[\mathbf{X}]$ for $\mathbf{X} = (\min(T_{A, \langle r \rangle}, T_1), I(T_{A, \langle r \rangle} < T_1)L)^\top$, where the expectation $\mathbb{E}[\cdot]$ corresponds to $\mathbb{E}_0[\cdot]$ (resp., $\mathbb{E}_1$) for the first (resp., second) coordinate.

We now apply Theorem 1(iii) to prove that the MSIS ratio estimator of $\alpha_{\langle r \rangle}$ has LE by verifying (22) for $i=1,2$. The first component of $\mathbf{X}$ satisfies $0 \le \min(T_{A, \langle r \rangle}, T_1) \le T_1$ for all $r$, so $\mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_1)] \le \mathbb{E}_0[T_1] < \infty$, where the finiteness follows from Corollary 1(c) of (Thorisson, 1985). Also, $\mathbb{P}_0(\lim_{r \to \infty} T_{A, \langle r \rangle} = \infty) = 1$ and $\mathbb{P}_0(T_1 < \infty) = 1$ ensure that $\mathbb{P}_0(\lim_{r \to \infty} \min(T_{A, \langle r \rangle}, T_1) = T_1) = 1$, so $\eta_{\langle r \rangle} = \mathbb{E}_0[\min(T_{A, \langle r \rangle}, T_1)] \to \mathbb{E}_0[T_1]$ as $r \to \infty$ by the dominated convergence theorem (Billingsley, 1995, Theorem 16.4). Moreover, because $0 < \mathbb{E}_0[U] \le \mathbb{E}_0[T_1] < \infty$, the numerator of $\alpha_{\langle r \rangle}$ satisfies $\lim_{r \to \infty} \eta_{\langle r \rangle} = \mathbb{E}_0[T_1] = \chi_1 \in (0, \infty)$, i.e., $\delta_1 = 0$ in (19). The second component of $\theta_{\langle r \rangle}$ satisfies $\gamma_{\langle r \rangle} \to 0$ as $r \to \infty$ by Theorem 1 of (Sadowsky, 1991), leading to $\delta_2 = -\delta$ in (19) for each $\delta > 0$. We then get that $\delta_0 = \delta$ in (19) because the estimand $\alpha_{\langle r \rangle} = \eta_{\langle r \rangle}/\gamma_{\langle r \rangle} \to \infty$ as $r \to \infty$. To prove (22) for $i=1,2$, observe that $g_1(\theta_{\langle r \rangle}) = 1/\theta_{\langle r \rangle, 2} = 1/\gamma_{\langle r \rangle}$ and $g_2(\theta_{\langle r \rangle}) = -\theta_{\langle r \rangle, 1}/\theta_{\langle r \rangle, 2}^2 = -\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}^2$. For $i=1$, using $\delta_1 = 0$ and $\delta_0 = \delta$ in (18) leads to

$$h_{1,1+\delta_1,1+\delta_0}(\theta_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle,1}^{1+\delta_1} g_1(\theta_{\langle r \rangle})}{g^{1+\delta_0}(\theta_{\langle r \rangle})} \right| = \left| \frac{\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}}{(\eta_{\langle r \rangle}/\gamma_{\langle r \rangle})^{1+\delta}} \right| = \frac{1}{\alpha_{\langle r \rangle}^{\delta}} \to 0 \quad \text{as } r \to \infty \quad (24)$$

for each $\delta > 0$. We next show that the naive MC estimator of $\theta_{\langle r \rangle, 1} = \eta_{\langle r \rangle}$ has BRE. Observe that $\mathbb{E}_0[(\min(T_{A, \langle r \rangle}, T_1))^2] \le \mathbb{E}_0[T_1^2]$ for all $r$, and $\mathbb{E}_0[T_1^2] < \infty$ does not depend on $r$, where Corollary 1(c) of (Thorisson, 1985) secures the finiteness because $\mathbb{E}_0[U^2] < \infty$ and $\mathbb{E}_0[V^2] < \infty$. Thus, the variance under $\mathbb{P}_0$ of $\min(T_{A, \langle r \rangle}, T_1)$ is bounded in $r$, and $\theta_{\langle r \rangle, 1} > 0$ is bounded away from 0 since $\eta_{\langle r \rangle} \ge \mathbb{E}_0[U] > 0$, securing our claim that we can estimate $\theta_{\langle r \rangle, 1}$ with BRE via naive MC. Consequently, (22) holds for $i=1$ by combining (24) with $RE_{1+\delta_1}[\hat{\theta}_{\langle r \rangle, n, 1}] = RE[\hat{\theta}_{\langle r \rangle, n, 1}] = O(1)$ as $r \to \infty$ by (6). To next verify (22) for $i=2$, where $\delta_2 = -\delta$ and $\delta_0 = \delta$ in (19), note that (18) becomes

$$h_{2,1+\delta_2,1+\delta_0}(\theta_{\langle r \rangle}) = \left| \frac{\theta_{\langle r \rangle,2}^{1+\delta_2} g_2(\theta_{\langle r \rangle})}{g^{1+\delta_0}(\theta_{\langle r \rangle})} \right| = \left| \frac{\gamma_{\langle r \rangle}^{1-\delta} (-\eta_{\langle r \rangle}/\gamma_{\langle r \rangle}^2)}{(\eta_{\langle r \rangle}/\gamma_{\langle r \rangle})^{1+\delta}} \right| = \frac{1}{\eta_{\langle r \rangle}^{\delta}}, \quad (25)$$

which, for each $\delta > 0$, is bounded in $r$ due to $\eta_{\langle r \rangle} \ge \mathbb{E}_0[U] > 0$. Under conditions from (Sadowsky, 1991), an IS estimator (via exponential twisting) for $\theta_{\langle r \rangle, 2} = \gamma_{\langle r \rangle}$ is LE, so for each $\delta > 0$, (6) yields $RE_{1+\delta_2}[\hat{\theta}_{\langle r \rangle, n, 2}] = RE_{1-\delta}[\hat{\theta}_{\langle r \rangle, n, 2}] = o(1)$ as $r \to \infty$ as in (7). Combining this with (25) proves (22) for $i=2$, so the MSIS estimator of $\alpha_{\langle r \rangle}$ is LE by Theorem 1(iii). $\square$
5 Making the Sufficient Conditions of Theorem 1 Also Necessary

We now identify situations that allow us to make the sufficient conditions of Theorem 1 to also be necessary, providing a more detailed understanding of efficiency and enabling us to prove the inefficiency of some estimators. One key reason why the conditions in Theorem 1 are only sufficient but not necessary is that, because each of $g_i(\theta_{\langle r\rangle})$, $g_j(\theta_{\langle r\rangle})$, and $\Sigma_{\langle r\rangle,i,j}$ in (14) could be negative or positive, the summands in (14) may have different signs, which could possibly lead in (14) to the asymptotic variance $\tau_{\langle r\rangle}^2$ having complete cancellation of its highest-order terms as $r \to \infty$, as we will see below in Example 16. Thus, although BRE, VRE, or LE may not hold if there is not complete cancellation of the highest-order terms in (14), when they do completely cancel out, the next-highest-order terms in (14) may result in BRE, VRE, or LE. The proof of Theorem 1 avoids this issue through two key ideas for BRE and VRE (similar comments also apply for LE): first, the triangle inequality is used in the line before (23) to take absolute values of each summand in (14); and second, the Cauchy-Schwarz inequality in (23) bounds the absolute value of each (possibly negative) covariance in (14) by the square root of the product of (positive) variances. The result then is that (23) sums only nonnegative terms, precluding the possibility of any cancellations, but this is at the expense of Theorem 1 providing only sufficient but not necessary conditions.

We next give a simple example in which the asymptotic variance $\tau_{\langle r\rangle}^2$ in (14) has complete cancellation of its highest-order terms.

**Example 16** (Complete cancellation of highest-order terms in asymptotic variance). As in Example 2 with $d=2$ so $\mathcal{J} = \{1,2\}$, let $\theta_{\langle r\rangle} = (\theta_{\langle r\rangle,1}, \theta_{\langle r\rangle,2})^\top$ with $\theta_{\langle r\rangle} = \mathbb{E}[\mathbf{X}_{\langle r\rangle}]$ and $\mathbf{X}_{\langle r\rangle} = (X_{\langle r\rangle,1}, X_{\langle r\rangle,2})^\top$. The goal is to estimate $\alpha_{\langle r\rangle} = g(\theta_{\langle r\rangle}) = \theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2}$, and $\hat{\theta}_{\langle r\rangle,n}$ averages n i.i.d. copies of $\mathbf{X}_{\langle r\rangle}$, as in (16). The partial derivatives of $g(\theta_{\langle r\rangle}) = \theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2}$ are $g_1(\theta_{\langle r\rangle}) = 1$ and $g_2(\theta_{\langle r\rangle}) = -1$, so (18) leads to $h_{1,1,1}(\theta_{\langle r\rangle}) = |\theta_{\langle r\rangle,1}/(\theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2})|$ and $h_{2,1,1}(\theta_{\langle r\rangle}) = |-\theta_{\langle r\rangle,2}/(\theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2})|$. Now suppose that $X_{\langle r\rangle,1} = A_{\langle r\rangle} + B_{\langle r\rangle}$ and $X_{\langle r\rangle,2} = A_{\langle r\rangle}$, where $A_{\langle r\rangle} \sim \mathcal{N}(m_{\langle r\rangle,A}, v_{\langle r\rangle,A})$ and $B_{\langle r\rangle} \sim \mathcal{N}(m_{\langle r\rangle,B}, v_{\langle r\rangle,B})$ are independent; we later specify values for $(m_{\langle r\rangle,i}, v_{\langle r\rangle,i})$, $i = A, B$. (Other distributions for $A_{\langle r\rangle}$ and $B_{\langle r\rangle}$ are also possible, as long as they allow for the means $m_{\langle r\rangle,i}$ and variances $v_{\langle r\rangle,i}$ given below.) Then $\theta_{\langle r\rangle,1} = m_{\langle r\rangle,A} + m_{\langle r\rangle,B}$ and $\theta_{\langle r\rangle,2} = m_{\langle r\rangle,A}$, so $\alpha_{\langle r\rangle} = m_{\langle r\rangle,B}$, which does not depend on the value of $m_{\langle r\rangle,A}$. Also, $\Sigma_{\langle r\rangle,1,1} = v_{\langle r\rangle,A} + v_{\langle r\rangle,B}$, $\Sigma_{\langle r\rangle,2,2} = v_{\langle r\rangle,A}$, and $\Sigma_{\langle r\rangle,1,2} = \mathbb{E}[X_{\langle r\rangle,1}X_{\langle r\rangle,2}] - \mathbb{E}[X_{\langle r\rangle,1}]\mathbb{E}[X_{\langle r\rangle,2}] = \mathbb{E}[A_{\langle r\rangle}^2] + \mathbb{E}[A_{\langle r\rangle}]\mathbb{E}[B_{\langle r\rangle}] - (m_{\langle r\rangle,A} + m_{\langle r\rangle,B})m_{\langle r\rangle,A} = \mathbb{E}[A_{\langle r\rangle}^2] - m_{\langle r\rangle,A}^2 = v_{\langle r\rangle,A}$. Thus, we get $g_1^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,1} = v_{\langle r\rangle,A} + v_{\langle r\rangle,B}$, $g_2^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,2,2} = v_{\langle r\rangle,A}$, and $g_1(\theta_{\langle r\rangle})g_2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,2} = -v_{\langle r\rangle,A}$, so (14) yields

$$ \tau_{\langle r\rangle}^2 = g_1^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,1} + g_2^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,2,2} + 2g_1(\theta_{\langle r\rangle})g_2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,2} = v_{\langle r\rangle,B} \quad (26) $$

as the asymptotic variance of $g(\hat{\theta}_{\langle r\rangle,n})$.

Now suppose that $v_{\langle r\rangle,A} = r^6$ and $v_{\langle r\rangle,B} = r^2$, so $g_1^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,1} = r^6 + r^2$, $g_2^2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,2,2} = r^6$, and $g_1(\theta_{\langle r\rangle})g_2(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,1,2} = -r^6$ have highest-order term $r^6$ of differing signs, but $\tau_{\langle r\rangle}^2 = r^2$ by (26). The highest-order terms in $\tau_{\langle r\rangle}^2$ then completely cancel out, and $\mathrm{RE}[g(\hat{\theta}_{\langle r\rangle,n})] = r/|m_{\langle r\rangle,B}|$. For any $m_{\langle r\rangle,A}$, if $m_{\langle r\rangle,B} = r^2$, then $g(\hat{\theta}_{\langle r\rangle,n})$ has VRE, but we now show that Theorem 1(ii) cannot secure this in certain situations. Suppose that $m_{\langle r\rangle,A} = r^2$, so $\mathrm{RE}[\hat{\theta}_{\langle r\rangle,n,1}] = \sqrt{\Sigma_{\langle r\rangle,1,1}/|\theta_{\langle r\rangle,1}|} = \sqrt{r^6 + r^2}/(2r^2) = \Theta(r) \to \infty$ as $r \to \infty$. Also, $h_{1,1,1} = |(m_{\langle r\rangle,A} + m_{\langle r\rangle,B})/m_{\langle r\rangle,B}| = 2$, so $h_{1,1,1} \mathrm{RE}[\hat{\theta}_{\langle r\rangle,n,1}] = \Theta(r)$; i.e., the sufficient condition (21) in Theorem 1(ii) does not hold for $i=1$; thus, $g(\hat{\theta}_{\langle r\rangle,n})$ has VRE, but Theorem 1(ii) cannot establish this. Example 21 in Section 6 gives numerical results. $\square$

Example 16 is a simple (linear) toy example of complete cancellation of highest-order terms of the asymptotic variance $\tau_{\langle r\rangle}^2$. But we can also observe similar behavior with more complicated models, such as $g(\theta_{\langle r\rangle}) = \theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2} + e^{-\theta_{\langle r\rangle,2}}$ or $g(\theta_{\langle r\rangle}) = \theta_{\langle r\rangle,1} - \theta_{\langle r\rangle,2} + e^{-(\theta_{\langle r\rangle,2}+\theta_{\langle r\rangle,2})^2}$ with the same $\theta_{\langle r\rangle} = (\theta_{\langle r\rangle,1}, \theta_{\langle r\rangle,2})^\top$ and $\mathbf{X}_{\langle r\rangle}$ as in Example 16, through more tedious and opaque calculations.

We can also have partial (but not complete) cancellation of highest-order terms of $\tau_{\langle r\rangle}^2$.
**Example 17** (Partial but not complete cancellation of highest-order terms in asymptotic variance). Consider the setting of the first paragraph of Example 16, but now $v_{\langle r \rangle,A} = v_{\langle r \rangle,B} = r^6$. Then $g_1^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle,1,1} = 2r^6$, $g_2^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle,2,2} = r^6$, and $g_1(\theta_{\langle r \rangle})g_2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle,1,2} = -r^6$ have highest-order terms of $\Theta(r^6)$ with differing signs. In contrast to Example 16, now $\tau_{\langle r \rangle}^2 = r^6$ by (26), so the different signs lead to partial (but not complete) cancellation of highest-order terms in $\tau_{\langle r \rangle}^2$. Now determining if $\hat{\theta}_{\langle r \rangle,n}$ has VRE, BRE or LE depends on the values of $m_{\langle r \rangle,i}$, $i = A, B$. For example, if $m_{\langle r \rangle,A} = r^2$ and $m_{\langle r \rangle,B} = r^3$, then Theorem 1(i) implies that $\hat{\theta}_{\langle r \rangle,n}$ has BRE because condition (20) holds since $h_{1,1,1} \text{RE}[\hat{\theta}_{\langle r \rangle,n,1}] = \sqrt{2}$ and $h_{2,1,1} \text{RE}[\hat{\theta}_{\langle r \rangle,n,2}] = 1$. In contrast, changing $m_{\langle r \rangle,B}$ to $r^2$ and keeping $m_{\langle r \rangle,A} = r^2$ lead to $\hat{\theta}_{\langle r \rangle,n}$ not having BRE because now $\text{RE}[\hat{\theta}_{\langle r \rangle,n}] = r$. But even though $h_{1,1,1} \text{RE}[\hat{\theta}_{\langle r \rangle,n,1}] = \sqrt{2}r^3$ and $h_{2,1,1} \text{RE}[\hat{\theta}_{\langle r \rangle,n,2}] = r^3$ show that (20) is violated, the fact that (20) is only a sufficient but not necessary condition in Theorem 1(i) prevents us from establishing the inefficiency of $\hat{\theta}_{\langle r \rangle,n}$. Example 22 in Section 6 provides numerical results. $\square$

To further illustrate the difficulties that can arise when the highest-order terms in $\tau_{\langle r \rangle}^2$ completely cancel, consider again the estimand $\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle}) = \theta_{\langle r \rangle,1} - \theta_{\langle r \rangle,2}$ in Example 16 for $\theta_{\langle r \rangle} = (\theta_{\langle r \rangle,1}, \theta_{\langle r \rangle,2})^\top$ with $\theta_{\langle r \rangle} = \mathbb{E}[\mathbf{X}_{\langle r \rangle}]$ and $\mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2})^\top$. But now suppose that $X_{\langle r \rangle,1} = A_{\langle r \rangle} + B_{\langle r \rangle} + C_{\langle r \rangle} + D_{\langle r \rangle}$ and $X_{\langle r \rangle,2} = A_{\langle r \rangle} + B_{\langle r \rangle} + C_{\langle r \rangle}$, where $A_{\langle r \rangle}, B_{\langle r \rangle}, C_{\langle r \rangle}$ and $D_{\langle r \rangle}$ are independent, with corresponding variances $v_{\langle r \rangle,A} \gg v_{\langle r \rangle,B} \gg v_{\langle r \rangle,C} \gg v_{\langle r \rangle,D}$. In this case, the estimator $g(\hat{\theta}_{\langle r \rangle,n})$ has asymptotic variance $\tau_{\langle r \rangle}^2 = \text{Var}[D_{\langle r \rangle}] = v_{\langle r \rangle,D}$, which is of the fourth-highest order. In general, the complete cancellation in $\tau_{\langle r \rangle}^2$ could continue for the first $k$ highest-order terms but not for $k+1$, and handling such a situation from the asymptotic efficiency properties of the estimators of the $\theta_{\langle r \rangle,i}$ ($i = 1, \dots, d$) becomes quite complicated, so we will not consider it.

To develop necessary conditions for BRE/VRE/LE, which will allow us to establish when BRE, VRE, or LE does not hold, we will restrict ourselves to settings that ensure the highest-order terms of $\tau_{\langle r \rangle}^2$ in (14) do not completely cancel out. To do this, it is sometimes helpful to rewrite (14) as

$$ \tau_{\langle r \rangle}^2 = \sum_{i=1}^{d} g_i^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, i} + 2 \sum_{i=1}^{d} \sum_{j=i+1}^{d} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} \quad (27) $$

to obtain below some specific situations that avoid such total cancellation of highest-order terms. Specifically, Sections 5.1 and 5.2 will identify settings in which the highest-order terms are all nonnegative (so there is no possible cancellation of highest-order terms), and Section 5.3 will allow for partial but not complete cancellation of the highest-order terms. The estimand $g(\theta_{\langle r \rangle})$ can similarly have cancellations of highest-order terms, but this is accounted for by $g^c(\theta_{\langle r \rangle})$ in (18).

## 5.1 Condition (28) Makes Theorem 1's Sufficient Conditions Also Necessary By Ensuring No Cancellation of the Highest-Order Terms in the Asymptotic Variance

A simple condition that ensures no cancellations of highest-order terms in (27) is that

$$ \Sigma_{\langle r \rangle, i, j} = 0 \text{ for all } i, j \in \mathcal{I} \text{ with } i \neq j. \quad (28) $$

One situation in which (28) is satisfied occurs when $\theta_{\langle r \rangle,i}$, $i \in \mathcal{I}$, are estimated independently, as with MSS in Example 4 using a different independent MC method for each $i$. But (28) can also hold when $\hat{\theta}_{\langle r \rangle,n,i}$, $i \in \mathcal{I}$, are dependent but the correlations vanish as $n \to \infty$. Specifically, let $\mathbf{Z} = (Z_1, Z_2, \dots, Z_d)^\top$ be a random vector having the limiting multivariate normal distribution in (12) with the covariance matrix $\Sigma_{\langle r \rangle}$. Adding Condition (28) then implies that $Z_i$ and $Z_j$ are uncorrelated for $i \neq j$, which means that $Z_i$ and $Z_j$ are independent because of the normality. For $i \neq j$, (28) allows for the estimators $\hat{\theta}_{\langle r \rangle,n,i}$ and $\hat{\theta}_{\langle r \rangle,n,j}$ to be dependent for finite $n$, but they have to be asymptotically uncorrelated as $n \to \infty$. Under (28), the sufficient conditions for BRE, VRE, or LE in Theorem 1 become necessary and sufficient in the following result.
**Theorem 2.** Assume that for each $r$, $\theta_{\langle r\rangle,i} \neq 0$ for each $i \in \mathcal{I}$, and $g$ in (9) has a nonzero differential at $\theta_{\langle r\rangle}$. Also, assume that Assumption A1 and (28) hold. Then the following hold as $r \to \infty$:

(i) $\hat{\alpha}_{\langle r\rangle,n}$ has BRE if and only if (20) holds.

(ii) $\hat{\alpha}_{\langle r\rangle,n}$ has VRE if and only if (21) holds.

(iii) When additionally Assumption A2 holds, $\hat{\alpha}_{\langle r\rangle,n}$ is LE if and only if (22) holds.

*Proof of Theorem 2.* Under (28), the double sum vanishes in (27), so $\tau_{\langle r\rangle}^2 = \sum_{i=1}^d g_i^2(\theta_{\langle r\rangle}) \Sigma_{\langle r\rangle,i,i}$, which sums only nonnegative terms, precluding the possibility of cancellations. Thus,

$$ (\mathrm{RE}[\hat{\alpha}_{\langle r\rangle,n}])^2 = \sum_{i=1}^{d} \frac{g_i^2(\theta_{\langle r\rangle})}{g^2(\theta_{\langle r\rangle})} \Sigma_{\langle r\rangle,i,i} = \sum_{i=1}^{d} \frac{\theta_{\langle r\rangle,i}^2 g_i^2(\theta_{\langle r\rangle})}{g^2(\theta_{\langle r\rangle})} \frac{\Sigma_{\langle r\rangle,i,i}}{\theta_{\langle r\rangle,i}^2} = \sum_{i=1}^{d} h_{i,1,1}^2(\theta_{\langle r\rangle}) (\mathrm{RE}[\hat{\theta}_{\langle r\rangle,n,i}])^2 $$

by (18) since each $\theta_{\langle r\rangle,i} \neq 0$. Parts (i) and (ii) then follow from (20) and (21), respectively.

For part (iii), (6) and (19) lead to $(\mathrm{RE}_{1+\delta_0}[\hat{\alpha}_{\langle r\rangle,n}])^2$ equaling

$$ \sum_{i=1}^{d} \frac{g_i^2(\theta_{\langle r\rangle})}{g^{2+2\delta_0}(\theta_{\langle r\rangle})} \Sigma_{\langle r\rangle,i,i} = \sum_{i=1}^{d} \frac{\theta_{\langle r\rangle,i}^{2+2\delta_i} g_i^2(\theta_{\langle r\rangle})}{g^{2+2\delta_0}(\theta_{\langle r\rangle})} \frac{\Sigma_{\langle r\rangle,i,i}}{\theta_{\langle r\rangle,i}^{2+2\delta_i}} = \sum_{i=1}^{d} h_{i,1+\delta_i,1+\delta_0}^2(\theta_{\langle r\rangle}) (\mathrm{RE}_{1+\delta_i}[\hat{\theta}_{\langle r\rangle,n,i}])^2 $$

by (18) since each $\theta_{\langle r\rangle,i} \neq 0$. Thus, part (iii) follows from (22). $\square$

<figure><img src="image_6.png" /></figure>

The following example describes a situation where (28) holds, so if $\hat{\alpha}_{\langle r\rangle,n}$ does not have BRE, then Theorem 2 can be applied to show this. In contrast, Theorem 1 provides only sufficient conditions for BRE, so Theorem 1 cannot establish this, showing the additional usefulness of Theorem 2.

**Example 18 (MSS ratio estimator).** Consider applying MSS from Example 4 to estimate a ratio $\alpha_{\langle r\rangle} = \theta_{\langle r\rangle,1}/\theta_{\langle r\rangle,2} = g(\theta_{\langle r\rangle,1}, \theta_{\langle r\rangle,2})$ for $g(x_1, x_1) = x_1/x_2$, with one MC method applied to estimate the numerator $\theta_{\langle r\rangle,1} \neq 0$ and another MC method used to estimate the denominator $\theta_{\langle r\rangle,2} \neq 0$, where the two MC methods are applied independently. The two MC methods do not necessarily have to be naive MC and IS, as with MSIS. Example 11 considers the specific case of applying MSIS for estimating MTTF of an HRMS, but here we allow more generally for MSS and ratios other than the MTTF, as in Examples 8 and 9, and not necessarily for an HRMS or a GI/GI/1 queue, as in Example 15. Now $h_{i,1,1}(\theta_{\langle r\rangle}) = 1$ for $i = 1,2$, in (18). Thus, by Theorem 2, the MSS ratio estimator $\hat{\alpha}_{\langle r\rangle,n} = g(\hat{\theta}_{\langle r\rangle,n,1}, \hat{\theta}_{\langle r\rangle,n,2})$ has BRE (resp., VRE) **if** and **only if** the estimators $\hat{\theta}_{\langle r\rangle,n,1}$ and $\hat{\theta}_{\langle r\rangle,n,2}$ both have BRE (resp., VRE). While Theorem 1 also covers the sufficiency, it does not handle the necessity. Thus, we can use Theorem 2 to show that when the MSS ratio estimator does **not** have BRE or VRE. For example, for the MTTF ratio in an HRMS, applying an ineffective IS (or any other) method to estimate the denominator may lead to the denominator's estimator not having BRE, so then the MSS estimator of the ratio does **not** have BRE. $\square$

<figure><img src="image_7.png" /></figure>

## 5.2 Less Restrictive Conditions (34) and (38) Also Ensure No Cancellation of the Highest-Order Terms in the Asymptotic Variance

We now weaken condition (28) that each asymptotic (non-variance) covariance is zero. The basic idea here is that in the double sum of (27) for the asymptotic variance $\tau_{\langle r\rangle}^2$, we will instead require that $g_i(\theta_{\langle r\rangle})g_j(\theta_{\langle r\rangle})\Sigma_{\langle r\rangle,i,j}$ for $i \neq j$ is “small” for $i$ and $j$ corresponding to the highest-order terms in the first sum of (27). To make this idea precise when studying BRE and VRE (we will later handle LE), define $\mathcal{A} \subseteq \mathcal{I}$ and its complement $\mathcal{B} = \mathcal{I} \setminus \mathcal{A}$ as sets of coordinates such that as $r \to \infty$,

$$ g_i(\theta_{\langle r\rangle}) \sqrt{\Sigma_{\langle r\rangle,i,i}} = O \left( g_k(\theta_{\langle r\rangle}) \sqrt{\Sigma_{\langle r\rangle,k,k}} \right) \quad \text{for all } i \in \mathcal{I} \text{ and } k \in \mathcal{A}, \qquad (29) $$

$$ g_l(\theta_{\langle r\rangle}) \sqrt{\Sigma_{\langle r\rangle,l,l}} = o \left( g_k(\theta_{\langle r\rangle}) \sqrt{\Sigma_{\langle r\rangle,k,k}} \right) \quad \text{for all } l \in \mathcal{B} \text{ and } k \in \mathcal{A}. \qquad (30) $$
The interpretation of (29) is that $\mathcal{A}$ is the subset of coordinates $k$ such that $g_k(\theta_{\langle r \rangle})\sqrt{\Sigma_{\langle r \rangle,k,k}}$ has a rate of the highest order as $r \to \infty$, and (30) specifies that indices $l \in \mathcal{B}$ have that $g_l(\theta_{\langle r \rangle})\sqrt{\Sigma_{\langle r \rangle,l,l}}$ is asymptotically strictly smaller than that for $k \in \mathcal{A}$. In the first sum of (27), each term $g_i^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,i,i} \ge 0$, and we claim that

$$g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k} \text{ is a term in (27) of the highest order as } r \to \infty \text{ for each } k \in \mathcal{A}. \quad (31)$$

To see why, note that (30) implies that as $r \to \infty$,

$$g_l^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,l,l} = o(g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}) \quad \text{for all } l \in \mathcal{B} \text{ and } k \in \mathcal{A}, \quad (32)$$

which verifies (31) for the terms in the single sum in (27). In the double sum in (27), each $\Sigma_{\langle r \rangle,i,j}$ is a non-variance covariance because $i \neq j$, and by the Cauchy-Schwarz inequality,

$$|g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,i,j}| \le |g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})|\sqrt{\Sigma_{\langle r \rangle,i,i}\Sigma_{\langle r \rangle,j,j}} \\ = \begin{cases} O(g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}) & \text{for each } k \in \mathcal{A} \text{ and any } i, j \in \mathcal{I}, \\ o(g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}) & \text{for each } k \in \mathcal{A} \text{ if } i \in \mathcal{B} \text{ or } j \in \mathcal{B} \text{ (or both)}, \end{cases} \quad (33)$$

as $r \to \infty$ by (29) and (30). Thus, in the double sum in (27), each summand (which has $i \neq j$) asymptotically (as $r \to \infty$) does not dominate $g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}$ for $k \in \mathcal{A}$, and is asymptotically negligible when $i \in \mathcal{B}$ or $j \in \mathcal{B}$. Thus, (31) holds, where the highest-order terms in (27) may also include some terms from the double sum with $i \neq j$ when $i \in \mathcal{A}$ and $j \in \mathcal{A}$ because of the first line of (33), but the second line of (33) and (32) ensure that $i \in \mathcal{B}$ or $j \in \mathcal{B}$ results in a term in (27) that does not have the highest order.

Even though (33) holds under (29) and (30), we may still have some $i, j \in \mathcal{A}$ with $i \neq j$ for which $g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,i,j}$ is negative and of the same asymptotic order as $g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}$ for $k \in \mathcal{A}$, and this could possibly lead to complete cancellation of the highest-order terms in (27). To preclude this possibility, we now relax (28) to instead assume that

$$g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,i,j} = o(g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle,k,k}) \quad \text{for all } i, j, k \in \mathcal{A} \text{ with } i \neq j, \quad (34)$$

which strengthens the first line of (33). (Note that if (28) holds, then the left side of (34) is zero, so (34) trivially holds, making it clear that (34) is weaker than (33).) In contrast to (28), the condition (34) allows for some types of non-zero covariance among the components of $\mathbf{Z} = (Z_1, Z_2, \dots, Z_d)^\top$ having the multivariate normal limit distribution in (12). We will see in Theorem 3 below that under (34), the sufficient conditions for BRE and VRE in Theorem 1(i)-(ii) become also necessary.

For a similar extension of Theorem 1(iii) for LE, we introduce additional notation. Assume that

$$\lim_{r \to \infty} \frac{1}{r} \ln |g_i(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, i, i}}| \quad \text{exists (and is finite) for all } i \in \mathcal{I}, \quad (35)$$

and define the set $\mathcal{A}_{\text{LE}} \subseteq \mathcal{I}$ and its complement $\mathcal{B}_{\text{LE}} = \mathcal{I} \setminus \mathcal{A}_{\text{LE}}$ such that

$$\lim_{r \to \infty} \frac{1}{r} \ln |g_i(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, i, i}}| \le \lim_{r \to \infty} \frac{1}{r} \ln |g_k(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, k, k}}| \quad \text{for all } i \in \mathcal{I} \text{ and } k \in \mathcal{A}_{\text{LE}}, \quad (36)$$

$$\lim_{r \to \infty} \frac{1}{r} \ln |g_l(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, l, l}}| < \lim_{r \to \infty} \frac{1}{r} \ln |g_k(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, k, k}}| \quad \text{for all } l \in \mathcal{B}_{\text{LE}} \text{ and } k \in \mathcal{A}_{\text{LE}}. \quad (37)$$

Thus, on a logarithmic scale, the contributions from coordinates $k \in \mathcal{A}_{\text{LE}}$ to the first sum in the asymptotic variance (27) in the limit (as $r \to \infty$) dominate the contributions from coordinates $l \in \mathcal{B}_{\text{LE}}$, and the proof of Theorem 3(iii) establishes that the same holds on a linear scale, i.e., (36) and (37) imply that (30) holds when the sets $\mathcal{A}$ and $\mathcal{B}$ are replaced with $\mathcal{A}_{\text{LE}}$ and $\mathcal{B}_{\text{LE}}$; see (41).

To get some insight into how $\mathcal{A}_{\text{LE}}$ and $\mathcal{B}_{\text{LE}}$ in (36) and (37) compare to $\mathcal{A}$ and $\mathcal{B}$ in (29) and (30), it is helpful to examine some specific examples, where we define $f_i(r) = g_i(\theta_{\langle r \rangle})\sqrt{\Sigma_{\langle r \rangle,i,i}}$ to simplify the notation in (29), (30), (36), and (37). Each example has $\theta_{\langle r \rangle} = \mathbb{E}[\mathbf{X}_{\langle r \rangle}]$ for $\mathbf{X}_{\langle r \rangle} = (X_{\langle r \rangle,1}, X_{\langle r \rangle,2}, \dots, X_{\langle r \rangle,d})^\top$, where
$X_{\langle r \rangle, i} \sim \mathcal{N}(\theta_{\langle r \rangle}, \Sigma_{\langle r \rangle, i, i})$ for each $i \in \mathcal{I}$, with $\theta_{\langle r \rangle}$ estimated by averaging $n$ i.i.d. copies of $\mathbf{X}_{\langle r \rangle}$, as in Example 2. Also, let $g(\theta_{\langle r \rangle}) = \sum_{i=1}^d \theta_{\langle r \rangle, i}$, so $g_i(\theta_{\langle r \rangle}) = 1$ for each $i \in \mathcal{I}$, in which case $f_i(r) = \sqrt{\sum_{\langle r \rangle, i, i}}$. We first take $d=5$ functions $f_i(r)$, $i \in \mathcal{I} = \{1, 2, \dots, 5\}$, shrinking to 0 at different rates (or remaining constant) as $r \to \infty$: $f_1(r) = e^{-r}$, $f_2(r) = r^p e^{-r}$ for some fixed $p \ge 0$, $f_3(r) = e^{-r+\sqrt{r}}$, $f_4(r) = e^{-2r}$, and $f_5(r) = r^{-q}$ for some fixed $q \ge 0$. Inspired by (Asmussen and Glynn, 2007, p. 159), the first three of the $f_i(r)$ all share the same asymptotic logarithmic behavior: $\lim_{r \to \infty} \frac{1}{r} \ln |f_i(r)| = -1$ for $i=1, 2, 3$. The other two functions have $\lim_{r \to \infty} \frac{1}{r} \ln |f_4(r)| = -2$ and $\lim_{r \to \infty} \frac{1}{r} \ln |f_5(r)| = 0$, so $\mathcal{A}_{\text{LE}} = \{5\}$ and $\mathcal{B}_{\text{LE}} = \{1, 2, 3, 4\}$ in (36) and (37). If we instead consider the linear scaling in (29) and (30), then $\mathcal{A} = \{5\} = \mathcal{A}_{\text{LE}}$ and $\mathcal{B} = \{1, 2, 3, 4\} = \mathcal{B}_{\text{LE}}$.

Now instead consider just the first $d=4$ functions, so $\mathcal{I} = \{1, 2, 3, 4\}$. Then $\mathcal{A}_{\text{LE}} = \{1, 2, 3\}$ and $\mathcal{B}_{\text{LE}} = \{4\}$ in (36) and (37), but $\mathcal{A} = \{3\}$ and $\mathcal{B} = \{1, 2, 4\}$ in (29) and (30). Even though the logarithmic and linear scalings lead to different dominant and nondominant sets (i.e., $\mathcal{A}_{\text{LE}} \neq \mathcal{A}$ and $\mathcal{B}_{\text{LE}} \neq \mathcal{B}$), we still have that on a linear scale, all functions with indices in $\mathcal{A}_{\text{LE}}$ dominate all functions with indices in $\mathcal{B}_{\text{LE}}$.

Next consider another example but with some functions growing to $\infty$ but others shrinking to 0 (or remaining constant) as $r \to \infty$. For $d=5$ functions, let $f_1(r) = e^r$, $f_2(r) = r^{-p}e^r$ for fixed $p \ge 0$, $f_3(r) = e^{r-\sqrt{r}}$, $f_4(r) = e^{-2r}$, and $f_5(r) = r^q$ for fixed $q \ge 0$. Then the first three all share the same asymptotic logarithmic behavior: $\lim_{r \to \infty} \frac{1}{r} \ln |f_i(r)| = 1$ for $i=1, 2, 3$. The other two functions have $\lim_{r \to \infty} \frac{1}{r} \ln |f_4(r)| = -2$ and $\lim_{r \to \infty} \frac{1}{r} \ln |f_5(r)| = 0$, so $\mathcal{A}_{\text{LE}} = \{1, 2, 3\}$ and $\mathcal{B}_{\text{LE}} = \{4, 5\}$ in (36) and (37). If we instead use the linear scaling of (29) and (30), then $\mathcal{A} = \{1\} \neq \mathcal{A}_{\text{LE}}$ and $\mathcal{B} = \{2, 3, 4, 5\} \neq \mathcal{B}_{\text{LE}}$, but we still have that all functions with indices in $\mathcal{A}_{\text{LE}}$ dominate all functions with indices in $\mathcal{B}_{\text{LE}}$ on a linear scale.

Now rather than (34) used before to handle BRE and VRE, assume instead for LE that

$$g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle, i, j} = o\left(\sum_{k \in \mathcal{A}_{\text{LE}}} g_k^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle, k, k}\right) \quad \text{for all } i, j \in \mathcal{A}_{\text{LE}} \text{ with } i \neq j. \quad (38)$$

**Theorem 3.** Assume that the function $g$ in (9) has a nonzero differential at $\theta_{\langle r \rangle}$, and that Assumption A1 holds. Then the following hold as $r \to \infty$:

(i) If (34) holds and if for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}$, where $\mathcal{A}$ is defined by (29) and (30), then $\hat{\alpha}_{\langle r \rangle, n}$ has BRE **if and only if** (20) holds with "for all $i \in \mathcal{A}$" replacing "for all $i \in \mathcal{I}$".

(ii) If (34) holds and if for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}$, then $\hat{\alpha}_{\langle r \rangle, n}$ has VRE **if and only if** (21) holds with "for all $i \in \mathcal{A}$" replacing "for all $i \in \mathcal{I}$".

(iii) Suppose now that Assumption A2 also holds, (35) and (38) hold instead of (34), and for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}_{\text{LE}}$, where $\mathcal{A}_{\text{LE}}$ is defined by (36) and (37). Then $\hat{\alpha}_{\langle r \rangle, n}$ is LE **if and only if** (22) holds with "for all $i \in \mathcal{A}_{\text{LE}}$" replacing "for all $i \in \mathcal{I}$".

*Proof of Theorem 3.* The single sum in (27) satisfies as $r \to \infty$

$$\sum_{i=1}^{d} g_i^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, i} = \sum_{k \in \mathcal{A}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k} + \sum_{l \in \mathcal{B}} g_l^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, l, l} = \left[ \sum_{k \in \mathcal{A}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k} \right] [1 + o(1)] \quad (39)$$

by (32) since $|\mathcal{I}| = d$ is fixed. Using the second line of (33) and (34) on the double sum in (27) gives

$$2 \sum_{i=1}^{d} \sum_{j=i+1}^{d} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} = o(g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}) \quad \text{for each } k \in \mathcal{A}. \quad (40)$$

Putting (39) and (40) into (27) yields $\tau_{\langle r \rangle}^2 = [\sum_{k \in \mathcal{A}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}] [1 + o(1)]$ as $r \to \infty$, where the sum
includes only nonnegative summands. Thus, as $r \to \infty$,

$$ (\mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}])^2 = \frac{1 + o(1)}{g^2(\boldsymbol{\theta}_{\langle r \rangle})} \sum_{k \in \mathcal{A}} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k} = [1 + o(1)] \sum_{k \in \mathcal{A}} \frac{\theta_{\langle r \rangle, k}^2 g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}}{g^2(\boldsymbol{\theta}_{\langle r \rangle}) \theta_{\langle r \rangle, k}^2} \\ = [1 + o(1)] \sum_{k \in \mathcal{A}} h_{k, 1, 1}^2(\boldsymbol{\theta}_{\langle r \rangle}) (\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, k}])^2 $$

by (18) since $\theta_{\langle r \rangle, k} \neq 0$ for each $k \in \mathcal{A}$. Parts (i) and (ii) then follow from (20) and (21), respectively, with “for all $i \in \mathcal{I}$” replaced by “for all $i \in \mathcal{A}$.”

For part (iii), (35) is equivalent to having for each $i \in \mathcal{I}$ that $\frac{1}{r} \ln |g_i(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, i, i}}| = c_i + o(1)$ as $r \to \infty$ for a constant $c_i \in (-\infty, \infty)$, so $|g_i(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, i, i}}| = e^{[c_i + o(1)]r}$. Also, for each $k \in \mathcal{A}_{\mathrm{LE}}$ and $l \in \mathcal{B}_{\mathrm{LE}}$, (37) ensures that $c_k > c_l$, so

$$ \frac{|g_l(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, l, l}}|}{|g_k(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, k, k}}|} = \frac{e^{[c_l + o(1)]r}}{e^{[c_k + o(1)]r}} = e^{-[c_k - c_l + o(1)]r} \to 0 $$

as $r \to \infty$ because $c_k - c_l > 0$ is fixed and $o(1) \to 0$. Hence, as $r \to \infty$,

$$ g_l(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, l, l}} = o(g_k(\boldsymbol{\theta}_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, k, k}}) \quad \text{for all } l \in \mathcal{B}_{\mathrm{LE}} \text{ and } k \in \mathcal{A}_{\mathrm{LE}}, \quad (41) $$

so the single sum in (27) satisfies

$$ \sum_{i=1}^{d} g_i^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, i} = \left[ \sum_{k \in \mathcal{A}_{\mathrm{LE}}} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k} \right] [1 + o(1)]. \quad (42) $$

In the double sum in (27), each term is a non-variance covariance because $i \neq j$, so

$$ |g_i(\boldsymbol{\theta}_{\langle r \rangle}) g_j(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}| = o\left(\sum_{k \in \mathcal{A}_{\mathrm{LE}}} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}\right) \quad \text{for all } i, j \in \mathcal{I} \text{ with } i \neq j, \quad (43) $$

as $r \to \infty$, which we next verify by examining all possible cases for $i, j \in \mathcal{I}$. First, when $i, j \in \mathcal{A}_{\mathrm{LE}}$ with $i \neq j$, (43) holds by (38). On the other hand, when $i \in \mathcal{B}_{\mathrm{LE}}$ or $j \in \mathcal{B}_{\mathrm{LE}}$, the Cauchy-Schwarz inequality bounds the left side of (43) by $|g_i(\boldsymbol{\theta}_{\langle r \rangle}) g_j(\boldsymbol{\theta}_{\langle r \rangle})| \sqrt{\Sigma_{\langle r \rangle, i, i} \Sigma_{\langle r \rangle, j, j}}$, so (43) follows from (41). Hence, (42) and (43) imply that $\tau_{\langle r \rangle}^2 = [\sum_{k \in \mathcal{A}_{\mathrm{LE}}} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}] [1 + o(1)]$ as $r \to \infty$. Consequently, as $r \to \infty$, by (19),

$$ (\mathrm{RE}_{1+\delta_0}[\hat{\alpha}_{\langle r \rangle, n}])^2 = \frac{1 + o(1)}{g^{2+2\delta_0}(\boldsymbol{\theta}_{\langle r \rangle})} \sum_{k \in \mathcal{A}_{\mathrm{LE}}} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k} = [1 + o(1)] \sum_{k \in \mathcal{A}_{\mathrm{LE}}} \frac{\theta_{\langle r \rangle, k}^{2+2\delta_k} g_k^2(\boldsymbol{\theta}_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}}{g^{2+2\delta_0}(\boldsymbol{\theta}_{\langle r \rangle}) \theta_{\langle r \rangle, k}^{2+2\delta_k}} \\ = [1 + o(1)] \sum_{k \in \mathcal{A}_{\mathrm{LE}}} h_{k, 1+\delta_k, 1+\delta_0}^2(\boldsymbol{\theta}_{\langle r \rangle}) (\mathrm{RE}_{1+\delta_k}[\hat{\theta}_{\langle r \rangle, n, k}])^2 $$

by (18) since $\theta_{\langle r \rangle, k} \neq 0$ for each $k \in \mathcal{A}_{\mathrm{LE}}$, Part (iii) then follows from (22) with “for each $i \in \mathcal{I}$” replaced by “for each $i \in \mathcal{A}$.” $\square$

We next present an example where showing that $\hat{\alpha}_{\langle r \rangle, n}$ does not have BRE or VRE can be handled by Theorem 3 but not by Theorems 1 or 2.

**Example 19 (Ratio estimator).** As in Example 18, consider estimating a ratio $\alpha_{\langle r \rangle} = \theta_{\langle r \rangle, 1}/\theta_{\langle r \rangle, 2} = g(\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2})$ for $g(x_1, x_1) = x_1/x_2$, but now we do not require independently estimating the numerator and denominator using MSS, so the estimators $\hat{\theta}_{\langle r \rangle, n, 1}$ and $\hat{\theta}_{\langle r \rangle, n, 2}$ may be dependent. Again, we have that $h_{i, 1, 1}(\boldsymbol{\theta}_{\langle r \rangle}) = 1$ in (18) for $i = 1, 2$. Thus, by Theorem 3(i) and (ii), if (34) holds, then the ratio estimator
$\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2})$ has BRE (resp., VRE) **if and only if** the estimators $\hat{\theta}_{\langle r \rangle, n, 1}$ and $\hat{\theta}_{\langle r \rangle, n, 2}$ both have BRE (resp., VRE). While Theorem 1 also covers the sufficiency, it does not handle the necessity. Also, assume that $\theta_{\langle r \rangle, 1}$ and $\theta_{\langle r \rangle, 2}$ are both estimated using naive MC with the same simulation runs, so the estimators $\hat{\theta}_{\langle r \rangle, n, 1}$ and $\hat{\theta}_{\langle r \rangle, n, 2}$ are dependent, which will typically violate (28), so Theorem 2 does not apply.

We now specialize $\alpha_{\langle r \rangle} = g(\theta_{\langle r \rangle, 1}, \theta_{\langle r \rangle, 2})$ as the MTTF for an HRMS, as in Example 11. In this case as $r \to \infty$, (Shahabuddin, 1994, Theorem 3 with $\epsilon = 1/r$) shows that $\theta_{\langle r \rangle, 1} = \Theta(r^q)$ and $\theta_{\langle r \rangle, 2} = \Theta(1/r^s)$ for some fixed $q \ge 1$ and $s > 0$, so $\alpha_{\langle r \rangle} = \Theta(r^{q+s})$. Also, when using CMC (by conditioning on the embedded DTMC, as described at the end of Example 8) but no other variance-reduction method (such as IS), for each fixed $r$, the estimator $\hat{\theta}_{\langle r \rangle, n} = (\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2})^\top$ satisfies the joint CLT (12), with as $r \to \infty$, $\Sigma_{\langle r \rangle, 1, 1} = O(1/r^t)$ for some fixed $t \ge 0$, $\Sigma_{\langle r \rangle, 2, 2} = \Theta(1/r^s)$ by (Shahabuddin, 1994, Theorem 1), and $\Sigma_{\langle r \rangle, 1, 2} = \Sigma_{\langle r \rangle, 2, 1} = O(1/r^{(s+t)/2})$ by the Cauchy-Schwarz inequality. For the partial derivatives $g_i(\theta_{\langle r \rangle})$ in (27), we have $g_1(\theta_{\langle r \rangle}) = 1/\theta_{\langle r \rangle, 2} = \Theta(r^s)$ and $g_2(\theta_{\langle r \rangle}) = -\theta_{\langle r \rangle, 1}/(\theta_{\langle r \rangle, 2})^2 = \Theta(r^{q+2s})$ as $r \to \infty$, so $g_1(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, 1, 1}} = \Theta(r^s)O(1/r^{t/2}) = O(r^{s-(t/2)})$ and $g_2(\theta_{\langle r \rangle}) \sqrt{\Sigma_{\langle r \rangle, 2, 2}} = \Theta(r^{q+2s})\Theta(1/r^{s/2}) = \Theta(r^{q+(3s/2)})$. Hence, (29) and (30) yield $\mathcal{A} = \{2\}$ and $\mathcal{B} = \{1\}$ since $q + (3s/2) > s - (t/2)$. Because $\mathcal{A}$ has only a single element, (34) vacuously holds because there does not exist $i, j \in \mathcal{A}$ with $i \neq j$. Now $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, 1}] = \sqrt{\Sigma_{\langle r \rangle, 1, 1}}/|\theta_{\langle r \rangle, 1}| = O(1/r^{t/2})/\Theta(r^q) = O(1/r^{(t/2)+q})$, so $\hat{\theta}_{\langle r \rangle, n, 1}$ has BRE (actually, VRE) because $t \ge 0$ and $q \ge 1$. But $\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, 2}] = \sqrt{\Sigma_{\langle r \rangle, 2, 2}}/|\theta_{\langle r \rangle, 2}| = \Theta(1/r^{s/2})/\Theta(1/r^s) = \Theta(r^{s/2}) \to \infty$ because $s > 0$, so $\hat{\theta}_{\langle r \rangle, n, 2}$ does not have BRE. Thus, Theorem 3 implies that when applying conditional MC but no other variance-reduction technique, the ratio estimator $\hat{\alpha}_{\langle r \rangle, n} = g(\hat{\theta}_{\langle r \rangle, n, 1}, \hat{\theta}_{\langle r \rangle, n, 2})$ does not have BRE, which (Shahabuddin, 1994, Theorem 4) also shows. $\square$

## 5.3 An Even Less Restrictive Condition (44) Allowing for Partial but Not Complete Cancellation of Highest-Order Terms in the Asymptotic Variance

Recall that the condition (34) strengthened the first line of (33), which led to Theorem 3 providing a necessary and sufficient condition for $\hat{\alpha}_{\langle r \rangle, n}$ to have BRE or VRE when (34) holds. The goal here is to come up with another condition even weaker than (34) that will allow us to determine situations when $\hat{\alpha}_{\langle r \rangle, n}$ does not have BRE or VRE (LE will also be handled separately) in the case when there is partial but not complete cancellation of the highest-order terms of $\tau_{\langle r \rangle}^2$, as occurs in Example 17. From the definitions in (29) and (30) of the index sets $\mathcal{A}$ and $\mathcal{B}$, we have that (32) and (33) hold, which show that $i \in \mathcal{B}$ or $j \in \mathcal{B}$ leads to non-dominant terms in (14). But in (14), we could still have for $i, j \in \mathcal{A}$ that $g_i(\theta_{\langle r \rangle})g_j(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle, i, j}$ is negative and of the same (highest) order as $g_i^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle, i, i}$ and $g_j^2(\theta_{\langle r \rangle})\Sigma_{\langle r \rangle, j, j}$, which was precluded by (34) in Section 5.2. Thus, these mixed signs in (14) could lead to all of the highest-order terms completely cancelling out, and we want to replace (34) with another (weaker) condition that prevents this. We thus assume that

$$R_{\langle r \rangle} \equiv \frac{\left| \sum_{i \in \mathcal{A}} \sum_{j \in \mathcal{A}} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} \right|}{\sum_{i \in \mathcal{A}} \sum_{j \in \mathcal{A}} |g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}|} \to c_0 \in (0, 1] \quad \text{as } r \to \infty, \quad (44)$$

where the double sums consider only $i, j \in \mathcal{A}$. Also, note that (44) is more related to the double sum (14) than to the double sum in (27) because (44) has double sums that allow $j = i$, as in (14), but the double sum in (27) requires $j \neq i$. Thus, (44) includes both variances and non-variance covariances. The triangle inequality implies that $R_{\langle r \rangle} \le 1$, and if the highest-order terms do completely cancel out, then $R_{\langle r \rangle} \to 0$ as $r \to \infty$, so (44) precludes this possibility (but there could still be a partial cancellation of highest-order terms, as in, e.g., Example 17). Under (44), we then have that the highest-order term of $\tau_{\langle r \rangle}^2$ is completely determined by indices $i, j \in \mathcal{A}$ in (14), so

$$\frac{\tau_{\langle r \rangle}^2}{\sum_{i \in \mathcal{A}} \sum_{j \in \mathcal{A}} |g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}|} \to c_0 \in (0, 1] \quad \text{as } r \to \infty, \quad (45)$$
where $c_0$ in (44) and (45) are the same. Also, under (44), we have that

$$ \tau_{\langle r \rangle}^2 = \left[ \sum_{i \in \mathcal{A}} \sum_{j \in \mathcal{A}} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} \right] [1 + o(1)] \quad \text{as } r \to \infty, \quad (46) $$

where we note that (46) does not have absolute values, which are present in (45).

While the condition (44) is appropriate for analyzing BRE and VRE, which consider asymptotics on a linear scale, we need to make adjustments to handle LE, which instead examines logarithmic asymptotics. Assuming that (35) holds as well as

$$ \lim_{r \to \infty} \frac{1}{r} \ln |g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}| \text{ exists (and is finite) for all } i, j \in \mathcal{A}_{\text{LE}}, \quad (47) $$

where $\mathcal{A}_{\text{LE}}$ is defined by (36) and (37), we impose the condition that

$$ \lim_{r \to \infty} \frac{\frac{1}{r} \ln \left( \left| \sum_{i \in \mathcal{A}_{\text{LE}}} \sum_{j \in \mathcal{A}_{\text{LE}}} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} \right| \right)}{\frac{1}{r} \ln \left( \sum_{i \in \mathcal{A}_{\text{LE}}} \sum_{j \in \mathcal{A}_{\text{LE}}} |g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}| \right)} = 1, \quad (48) $$

which ensures that the highest-order terms in $\sum_{i \in \mathcal{A}_{\text{LE}}} \sum_{j \in \mathcal{A}_{\text{LE}}} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j}$ do not completely cancel out on a logarithmic scale. Thus, the exponential rates of the terms inside the logs of the numerator and denominator of (48) are the same.

**Theorem 4.** Assume that the function $g$ in (9) has a nonzero differential at $\theta_{\langle r \rangle}$, and that Assumption A1 holds. Then the following hold as $r \to \infty$:

(i) If (44) holds and if for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}$, where $\mathcal{A}$ is defined by (29) and (30), then $\hat{\alpha}_{\langle r \rangle, n}$ has BRE **if and only if** (20) holds with “for all $i \in \mathcal{A}$” replacing “for all $i \in \mathcal{I}$”.

(ii) If (44) holds and if for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}$, then $\hat{\alpha}_{\langle r \rangle, n}$ has VRE **if and only if** (21) holds with “for all $i \in \mathcal{A}$” replacing “for all $i \in \mathcal{I}$”.

(iii) Suppose now that Assumption A2 also holds, along with (35), (47), and (48). In addition, assume that for each $r$, $\theta_{\langle r \rangle, i} \neq 0$ for each $i \in \mathcal{A}_{\text{LE}}$, where $\mathcal{A}_{\text{LE}}$ is defined by (36) and (37). Then $\hat{\alpha}_{\langle r \rangle, n}$ is LE **if and only if** (22) holds with “for all $i \in \mathcal{A}_{\text{LE}}$” replacing “for all $i \in \mathcal{I}$”.

**Proof of Theorem 4.** Because (44) ensures that there is not complete cancellation of the highest-order terms in $\tau_{\langle r \rangle}^2$, the first line of (33) yields as $r \to \infty$, $\sum_{i \in \mathcal{A}} \sum_{j \in \mathcal{A}} g_i(\theta_{\langle r \rangle}) g_j(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, i, j} = \Theta(\sum_{k \in \mathcal{A}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k})$. Thus, $g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}$ for $k \in \mathcal{A}$ determine the exact order of magnitude of $\tau_{\langle r \rangle}^2$ in (46), so $\tau_{\langle r \rangle}^2 = \Theta(\sum_{k \in \mathcal{A}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}) = \sum_{k \in \mathcal{A}} \Theta(g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k})$ because all summands are nonnegative and $|\mathcal{A}| \le d < \infty$. We then get

$$ \begin{aligned} (\mathrm{RE}[\hat{\alpha}_{\langle r \rangle, n}])^2 &= \frac{1}{g^2(\theta_{\langle r \rangle})} \sum_{k \in \mathcal{A}} \Theta(g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}) = \sum_{k \in \mathcal{A}} \Theta \left( \frac{\theta_{\langle r \rangle, k}^2 g_k^2(\theta_{\langle r \rangle})}{g^2(\theta_{\langle r \rangle})} \frac{\sum_{\langle r \rangle, k, k}}{\theta_{\langle r \rangle, k}^2} \right) \\ &= \sum_{k \in \mathcal{A}} \Theta(h_{k, 1, 1}^2(\theta_{\langle r \rangle}) (\mathrm{RE}[\hat{\theta}_{\langle r \rangle, n, k}])^2) \end{aligned} $$

by (18) since we assumed that $\theta_{\langle r \rangle, k} \neq 0$ for each $k \in \mathcal{A}$. Hence, parts (i) and (ii) follow from (20) and (21), respectively, with “for all $i \in \mathcal{I}$” replaced by “for all $i \in \mathcal{A}$.”

For part (iii), as shown before in the proof of Theorem 3(iii), (35) implies the single sum in (27) satisfies (42). Also, (47) and (48) ensure that there is not complete cancellation of the highest-order terms in the
numerator of (48), which implies that $\tau_{\langle r \rangle}^2 = \Theta(\sum_{k \in \mathcal{A}_{\text{LE}}} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}) = \sum_{k \in \mathcal{A}_{\text{LE}}} \Theta(g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k})$. Consequently, as $r \to \infty$, by (19),

$$ (\text{RE}_{1+\delta_0}[\hat{\alpha}_{\langle r \rangle, n}])^2 = \frac{1}{g^{2+2\delta_0}(\theta_{\langle r \rangle})} \sum_{k \in \mathcal{A}_{\text{LE}}} \Theta(g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}) = \sum_{k \in \mathcal{A}_{\text{LE}}} \Theta\left(\frac{\theta_{\langle r \rangle, k}^{2+2\delta_k} g_k^2(\theta_{\langle r \rangle}) \Sigma_{\langle r \rangle, k, k}}{g^{2+2\delta_0}(\theta_{\langle r \rangle}) \theta_{\langle r \rangle, k}^{2+2\delta_k}}\right) \\ = \sum_{k \in \mathcal{A}_{\text{LE}}} \Theta\left(h_{k,1+\delta_k,1+\delta_0}^2(\theta_{\langle r \rangle})(\text{RE}_{1+\delta_k}[\hat{\theta}_{\langle r \rangle, n, k}])^2\right) $$

by (18) since $\theta_{\langle r \rangle, k} \neq 0$ for each $k \in \mathcal{A}_{\text{LE}}$. Thus, part (iii) follows from (22) with “for each $i \in \mathcal{I}$” replaced by “for each $i \in \mathcal{A}_{\text{LE}}.” $\square$

# 6 NUMERICAL RESULTS

Our previous paper (Nakayama and Tuffin, 2023), which contains only Theorem 1 (restricted to the setting of Example 2) but not the other ones in Section 5, gives numerical results for two examples fitting that framework: the mean hitting time to a failed state for a semi-Markov process modeling a reliability system with two components, and also the MTTF for an HRMS, as described in Example 11. We rather focus here on Examples 13, 16, and 17 related to the new results in Section 5 not covered in (Nakayama and Tuffin, 2023).

**Example 20.** As in Example 11, we consider an HRMS made of $c = 3$ types of component and 5 components of each type. Each component has an exponentially distributed lifetime with rate $\epsilon = 1/r$, with all nonfailed components subject to failure. The time to repair a failed component is exponential with rate 1, where repair starts as soon as the component has failed. Failure and repair times are all independent. The system is down whenever fewer than two components of any one type are operational. Our goal is to compute, as described in Example 13, the derivative of the MTTF with respect to the failure rate $\lambda$ of the first component type. We simulate the embedded DTMC (i.e., CMC; see Example 8) to estimate the four components in $\theta = (\eta(\lambda), \gamma(\lambda), \eta'(\lambda), \gamma'(\lambda))^{\top}$, with the components of $\theta$ estimated independently to illustrate Theorem 2. We considered two different approaches, naive MC and IS, for estimating the vector $\theta = \theta_{\langle r \rangle}$ by averaging $n = 10^6$ i.i.d. copies a random vector $\mathbf{X}_{\langle r \rangle}$ with independent components and mean $\theta$, putting us in the framework of Example 2, with computation budget $n$. For naive MC, all elements of $\theta$ are estimated independently without IS. For IS, all components of $\theta$ are estimated independently using the IS technique balanced failure biasing (BFB) (Shahabuddin, 1994). Under BFB, in a state from which both failures and repairs are possible, the cumulative probability assigned to all failure transitions is $\rho = 0.8$, equally allocating $\rho$ across the individual failure transitions. Repair transitions get a total probability of $1-\rho$, distributed to individual repair transitions in proportion to their original probabilities. The precise description of each estimator and of the variance computation can be found in (Nakayama et al., 1994; Nakayama, 1995). Table 1 gives the estimated relative error $\hat{\text{RE}}$, obtained from the estimated variance and mean, for each estimator $\hat{\theta}_{\langle r \rangle, n, i}$, $i \in \mathcal{I} = \{1, 2, 3, 4\}$, and $\hat{\alpha}_{\langle r \rangle, n}$. Also, the table presents $h\text{RE}_i = h_{i,1,1}(\hat{\theta}_{\langle r \rangle, n})\hat{\text{RE}}[\hat{\theta}_{\langle r \rangle, n, i}]$, $i \in \mathcal{I}$, as estimates of (20) used in Theorem 2(i). As $r \to \infty$ (equivalently $\epsilon \to 0$), the naive MC estimators of $\hat{\theta}_{\langle r \rangle, n, i}$, $i = 2, 4$, have unbounded RE, which combined with the behavior of the corresponding $h_{i,1,1}(\hat{\theta}_{\langle r \rangle, n})$ lead to the same for $\hat{\alpha}_{\langle r \rangle, n}$, as in Theorem 2(i). In contrast, all of the estimators have BRE when applying IS. It is also interesting to note that IS is not as efficient as naive MC for $\epsilon = 10^{-1}$, due to the fact that a system failure within a cycle is not a very rare event, so applying BFB is not appropriate. For $\epsilon = 10^{-2}$, the naive MC and IS confidence intervals do not intersect, which occurs from an underestimation of the variance with naive MC, a well-known rare-event issue; e.g., see (Smith, 2001), (Glynn et al., 2009, Section 4.4.1), (Asmussen and Glynn, 2007, pp. 134–135), and (Arief et al., 2021). $\square$

**Example 21.** We numerically illustrate Example 16 with full cancellation of highest-order terms of the asymptotic variance by simulating with $n = 10^6$ independent copies of A and B in the model for various values of $r$. Table 2 gives the results for $m_{\langle r \rangle, A} = m_{\langle r \rangle, B} = r^2$, $v_{\langle r \rangle, A} = r^6$, and $v_{\langle r \rangle, B} = r^2$.
<table><caption>Table 1: Results from estimating the sensitivity of the MTTF with respect to the individual failure rate, for the HRMS of Example 20, from $n = 10^6$ independent cycles for the estimation of each parameter, where “nan” means “not a number”. To simplify the table, we define $hRE_i = h_{i,1,1}(\hat{\theta}_{\langle r \rangle,n})\hat{RE}[\hat{\theta}_{\langle r \rangle,n,i}]$.</caption><thead><tr><th>Method</th><th>1/r</th><th>RE[θ̂<sub>⟨r⟩,n,1</sub>]<br>(hRE<sub>1</sub>)</th><th>RE[θ̂<sub>⟨r⟩,n,2</sub>]<br>(hRE<sub>2</sub>)</th><th>RE[θ̂<sub>⟨r⟩,n,3</sub>]<br>(hRE<sub>3</sub>)</th><th>RE[θ̂<sub>⟨r⟩,n,4</sub>]<br>(hRE<sub>4</sub>)</th><th>α̂<sub>⟨r⟩,n</sub></th><th>RE[α̂<sub>⟨r⟩,n</sub>]</th><th>95% CI</th></tr></thead><tbody><tr><td>Naive</td><td>10<sup>-1</sup></td><td>1.32e+00<br>(2.22e+00)</td><td>1.02e+01<br>(2.73e+01)</td><td>3.52e+00<br>(2.37e+00)</td><td>1.12e+01<br>(1.88e+01)</td><td>-9.17e+03</td><td>6.65e+01</td><td>(-9.77e+03, -8.57e+03)</td></tr><tr><td>Naive</td><td>10<sup>-2</sup></td><td>1.00e+00<br>(9.83e-01)</td><td>7.07e+02<br>(1.40e+03)</td><td>1.20e+00<br>(2.30e-02)</td><td>4.33e+02<br>(4.25e+02)</td><td>-5.09e+09</td><td>2.93e+03</td><td>(-1.97e+10, 9.52e+09)</td></tr><tr><td>Naive</td><td>10<sup>-3</sup></td><td>1.00e+00<br>(nan)</td><td>nan<br>(nan)</td><td>1.00e+00<br>(nan)</td><td>nan<br>(nan)</td><td>nan</td><td>nan</td><td>(nan, nan)</td></tr><tr><td>Naive</td><td>10<sup>-4</sup></td><td>1.00e+00<br>(nan)</td><td>nan<br>(nan)</td><td>1.00e+00<br>(nan)</td><td>nan<br>(nan)</td><td>nan</td><td>nan</td><td>(nan, nan)</td></tr><tr><td>IS</td><td>10<sup>-1</sup></td><td>6.15e+01<br>(8.40e+01)</td><td>2.00e+01<br>(4.74e+01)</td><td>2.20e+02<br>(8.05e+01)</td><td>2.36e+02<br>(3.22e+02)</td><td>-1.44e+04</td><td>6.92e+02</td><td>(-2.42e+04, -4.64e+03)</td></tr><tr><td>IS</td><td>10<sup>-2</sup></td><td>2.04e+00<br>(1.93e+00)</td><td>5.58e+00<br>(1.09e+01)</td><td>1.31e+00<br>(6.88e-02)</td><td>5.46e+00<br>(5.17e+00)</td><td>-8.47e+08</td><td>2.44e+01</td><td>(-8.67e+08, -8.27e+08)</td></tr><tr><td>IS</td><td>10<sup>-3</sup></td><td>2.21e+00<br>(2.06e+00)</td><td>7.01e+00<br>(1.36e+01)</td><td>1.00e+00<br>(6.68e-02)</td><td>7.09e+00<br>(6.62e+00)</td><td>-8.20e+13</td><td>3.04e+01</td><td>(-8.45e+13, -7.96e+13)</td></tr><tr><td>IS</td><td>10<sup>-4</sup></td><td>2.23e+00<br>(2.08e+00)</td><td>7.22e+00<br>(1.40e+01)</td><td>1.00e+00<br>(6.76e-02)</td><td>7.34e+00<br>(6.85e+00)</td><td>-8.17e+18</td><td>3.14e+01</td><td>(-8.42e+18, -7.92e+18)</td></tr><tr><td>IS</td><td>10<sup>-5</sup></td><td>2.24e+00<br>(2.08e+00)</td><td>7.24e+00<br>(1.40e+01)</td><td>1.00e+00<br>(6.77e-02)</td><td>7.37e+00<br>(6.87e+00)</td><td>-8.17e+23</td><td>3.15e+01</td><td>(-8.42e+23, -7.92e+23)</td></tr></tbody></table>

<table><caption>Table 2: Results of Example 16, for various values of $r$ from $n = 10^6$ independent copies of A and B.</caption><thead><tr><th>r</th><th>RE[θ̂<sub>⟨r⟩,n,1</sub>]</th><th>h<sub>1,1,1</sub>(θ<sub>⟨r⟩</sub>)RE[θ̂<sub>⟨r⟩,n,1</sub>]</th><th>RE[θ̂<sub>⟨r⟩,n,2</sub>]</th><th>h<sub>2,1,1</sub>(θ<sub>⟨r⟩</sub>)RE[θ̂<sub>⟨r⟩,n,2</sub>]</th><th>RE[α̂<sub>⟨r⟩,n</sub>]</th></tr></thead><tbody><tr><td>1</td><td>7.07e-01</td><td>1.41e+00</td><td>1.00e+00</td><td>1.00e+00</td><td>1.00e+00</td></tr><tr><td>5</td><td>2.50e+00</td><td>5.01e+00</td><td>4.99e+00</td><td>5.00e+00</td><td>2.00e-01</td></tr><tr><td>10</td><td>4.99e+00</td><td>1.00e+01</td><td>9.96e+00</td><td>1.00e+01</td><td>9.99e-02</td></tr><tr><td>100</td><td>4.88e+01</td><td>1.00e+02</td><td>9.54e+01</td><td>1.00e+02</td><td>9.99e-03</td></tr></tbody></table>

It can be checked that $\hat{\alpha}_{\langle r \rangle,n}$ has VRE even though $\hat{\theta}_{\langle r \rangle,n,1}$ and $\hat{\theta}_{\langle r \rangle,n,2}$ have unbounded RE and the $h_{i,1,1}(\theta_{\langle r \rangle})RE[\hat{\theta}_{\langle r \rangle,n,i}]$ (for $i \in \{1, 2\}$) are unbounded. This is due to the total cancellation of highest-order terms of the asymptotic variance. Establishing sufficient (or necessary-and-sufficient) conditions in this case requires analyzing the non-highest-order terms in a Taylor expansion of the variances, which is complicated and impractical in general.

**Example 22.** We similarly illustrate Example 17 with in this case partial (but not complete) cancellation of highest-order terms of the asymptotic variance, still by simulating $n = 10^6$ independent copies of A and B in the model with $m_{\langle r \rangle,A} = m_{\langle r \rangle,B} = r^2$ and $v_{\langle r \rangle,A} = v_{\langle r \rangle,B} = r^6$. Table 3 presents the results for various values of $r$. Condition (20) of Theorem 1(i) is violated, rendering the theorem unable to conclude anything. But from Theorem 4(i) the estimator of $\alpha$ has unbounded RE.

<table><caption>Table 3: Results of Example 17, for various values of $r$ from $n = 10^6$ independent copies of A and B.</caption><thead><tr><th>r</th><th>RE[θ̂<sub>⟨r⟩,n,1</sub>]</th><th>h<sub>1,1,1</sub>(θ<sub>⟨r⟩</sub>)RE[θ̂<sub>⟨r⟩,n,1</sub>]</th><th>RE[θ̂<sub>⟨r⟩,n,2</sub>]</th><th>h<sub>2,1,1</sub>(θ<sub>⟨r⟩</sub>)RE[θ̂<sub>⟨r⟩,n,2</sub>]</th><th>RE[α̂<sub>⟨r⟩,n</sub>]</th></tr></thead><tbody><tr><td>1</td><td>7.07e-01</td><td>1.41e+00</td><td>1.00e+00</td><td>1.00e+00</td><td>1.00e+00</td></tr><tr><td>5</td><td>3.54e+00</td><td>7.09e+00</td><td>4.99e+00</td><td>5.02e+00</td><td>5.01e+00</td></tr><tr><td>10</td><td>7.07e+00</td><td>1.42e+01</td><td>9.96e+00</td><td>1.01e+01</td><td>1.01e+01</td></tr><tr><td>100</td><td>7.12e+01</td><td>1.51e+02</td><td>9.54e+01</td><td>1.07e+02</td><td>1.07e+02</td></tr></tbody></table>
Although we presented only a few numerical examples illustrating our theory, the literature provides many other relevant results. For example, see (Goyal et al., 1992; Nakayama and Tuffin, 2019; Shahabuddin, 1994; Bassamboo et al., 2008; Li et al., 2024), among many others.

# 7 CONCLUDING REMARKS

This paper examined the efficiency of estimating a smooth function of unknown parameters, when some but not necessarily all of the parameters relate to rare events. The parameters could be means, quantiles, or other quantities, with the assumption that their estimators obey a joint central limit theorem. We considered different notions of efficiency: bounded relative error, vanishing relative error, and logarithmic efficiency. We illustrated that, at first sight counter-intuitively, efficiency of the parameters' estimators does not ensure efficiency of the function of parameters' estimators, and efficiency of these parameters' estimators is not always required for efficient function estimation. Based on this, Theorem 1 provided sufficient conditions to ensure BRE, VRE, and LE, and we provided many practical examples (Section 4) demonstrating the broad applicability of the theory. We also identified situations (Section 5) that allow us to strengthen Theorem 1's sufficient conditions to also be necessary (Theorems 2, 3, and 4), providing a deeper understanding of efficiency and enabling us to establish when certain estimators are inefficient.

# Acknowledgments

This work has been supported in part by the National Science Foundation under Grant No. CMMI-2345330. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

# References

* Apostol, T. M. (1974). *Mathematical Analysis*. Addison-Wesley, Reading, MA, 2nd edition.

* Arief, M., Huang, Z., Kumar, G. K. S., Bai, Y., He, S., Ding, W., Lam, H., and Zhao, D. (2021). Deep probabilistic accelerated evaluation: A robust certifiable rare-event simulation methodology for black-box safety-critical systems. In Banerjee, A. and Fukumizu, K., editors, *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics*, volume 130, pages 595–603. PMLR.

* Asmussen, S. and Glynn, P. (2007). *Stochastic Simulation: Algorithms and Analysis*. Springer, New York.

* Bassamboo, A., Juneja, S., and Zeevi, A. (2008). Portfolio credit risk with extremal dependence: Asymptotic analysis and efficient simulation. *Operations Research*, 56(3):593–606.

* Billingsley, P. (1995). *Probability and Measure*. John Wiley and Sons, New York, 3rd edition.

* Blanchet, J. and Lam, H. (2011). Rare event simulation techniques. In Jain, S., Creasey, R., Himmelspach, J., White, K., and Fu, M., editors, *Proceedings of the 2011 Winter Simulation Conference*, page 146–160, Piscataway, New Jersey. Institute of Electrical and Electronics Engineers.

* Budhiraja, A. and Dupuis, P. (2019). *Analysis and Approximation of Rare Events: Representations and Weak Convergence Methods*. Springer US, New York, NY.

* De Maio, F. G. (2007). Income inequality measures. *Journal of Epidemiology & Community Health*, 61(10):849–852.

* Ferguson, T. (1999). Asymptotic joint distribution of sample mean and a sample quantile. https://www.math.ucla.edu/~tom/papers/unpublished/meanmed.pdf.
Glynn, P. W. (1990). Likelihood ratio derivative estimators for stochastic systems. *Communications of the ACM*, 33:75-84.

Glynn, P. W. (1996). Importance sampling for Monte Carlo estimation of quantiles. In Ermakov, S. M. and Melas, V. B., editors, *Mathematical Methods in Stochastic Simulation and Experimental Design: Proceedings of the 2nd St. Petersburg Workshop on Simulation*, pages 180–185, St. Petersburg, Russia. Publishing House of St. Petersburg Univ.

Glynn, P. W., Fan, L., Fu, M. C., Hu, J.-Q., and Peng, Y. (2020). Technical note—central limit theorems for estimated functions at estimated points. *Operations Research*, 68(5):1557–1563.

Glynn, P. W. and Iglehart, D. L. (1993). Conditions for the applicability of the regenerative method. *Management Science*, 39:1108–1111.

Glynn, P. W., L'Ecuyer, P., and Ades, M. (1991). Gradient estimation for ratios. In *1991 Winter Simulation Conference Proceedings*, pages 986–993. IEEE.

Glynn, P. W., Nakayama, M. K., and Tuffin, B. (2017). On the estimation of the mean time to failure by simulation. In Chan, W. K. V., A.D'Ambrogio, Zacharewicz, G., Mustafee, N., Wainer, G., and Page, E., editors, *Proceedings of the 2017 Winter Simulation Conference*, pages 1844–1855, Piscataway, NJ. Institute of Electrical and Electronics Engineers.

Glynn, P. W., Peng, Y., Fu, M. C., and Hu, J.-Q. (2021). Computing sensitivities for distortion risk measures. *INFORMS Journal on Computing*, 33(4):1520–1532.

Glynn, P. W., Rubino, G., and Tuffin, B. (2009). *Robustness Properties and Confidence Interval Reliability Issues*, chapter Chapter 4, pages 63–84. Wiley.

Goyal, A., Shahabuddin, P., Heidelberger, P., Nicola, V., and Glynn, P. W. (1992). A unified framework for simulating Markovian models of highly dependable systems. *IEEE Transactions on Computers*, C-41(1):36–51.

Kalashnikov, V. (1994). *Topics on Regenerative Processes*. CRC Press, Boca Raton.

Klaassen, P. and van Eeghen, I. (2009). *Economic Capital: How It Works, and What Every Manager Needs to Know*. Elsevier, Burlington, MA.

L'Ecuyer, P., Blanchet, J. H., Tuffin, B., and Glynn, P. W. (2010). Asymptotic robustness of estimators in rare-event simulation. *ACM Transactions on Modeling and Computer Simulation*, 20(1):Article 6.

Li, Y., Kaplan, Z. T., and Nakayama, M. K. (2024). Monte carlo methods for economic capital. *INFORMS Journal on Computing*, 36(1):266–284.

Lin, P.-E., Wu, K.-T., and Ahmad, I. A. (1980). Asymptotic joint distribution of sample quantiles and sample mean with applications. *Communications in Statistics - Theory and Methods*, A9(1):51–60.

Lütkebohmert, E. (2009). *Concentration Risk in Credit Portfolios*. Springer, Berlin.

Nakayama, M. K. (1995). Asymptotics of likelihood ratio derivative estimators in simulations of highly reliable Markovian systems. *Management Science*, 41:524–554.

Nakayama, M. K. (1998). On derivative estimation of the mean time to failure in simulations of highly reliable Markovian systems. *Operations Research*, 46:285–290.

Nakayama, M. K., Goyal, A., and Glynn, P. W. (1994). Likelihood ratio sensitivity analysis for Markovian models of highly dependable systems. *Operations Research*, 42:137–157.
Nakayama, M. K. and Tuffin, B. (2019). Efficient estimation of the mean hitting time to a set of a regenerative system. In Mustafee, N., Bae, K.-H., Lazarova-Molnar, S., Rabe, M., Szabo, C., Haas, P., and Son, Y.-J., editors, *Proceedings of the 2019 Winter Simulation Conference*, pages 416–427, Piscataway, New Jersey. Institute of Electrical and Electronics Engineers.

Nakayama, M. K. and Tuffin, B. (2023). Efficiency of estimating functions of means in rare-event contexts. In *2023 Winter Simulation Conference (WSC)*, pages 351–362. IEEE.

Nakayama, M. K. and Tuffin, B. (2024). Sufficient Conditions for Central Limit Theorems and Confidence Intervals for Randomized Quasi-Monte Carlo Methods. *ACM Transactions on Modeling and Computer Simulation*, 34(3):1–38.

Nakayama, M. K., Tuffin, B., and L'Ecuyer, P. (2025). Central limit theorem for a randomized quasi-monte carlo estimator of a smooth function of means. In Azar, E., Djanatliev, A., Harper, A., Kogler, C., Ramamohan, V., Anagnostou, A., and Taylor, S. J. E., editors, *Proceedings of the 2025 Winter Simulation Conference*. IEEE. to appear.

Rubino, G. and Tuffin, B. (2009a). Markovian models for dependability analysis. In Rubino, G. and Tuffin, B., editors, *Rare Event Simulation using Monte Carlo Methods*, pages 125–144. John Wiley & Sons, Chichester, UK.

Rubino, G. and Tuffin, B., editors (2009b). *Rare Event Simulation using Monte Carlo Methods*. John Wiley, Chichester, UK.

Sadowsky, J. S. (1991). Large deviations theory and efficient simulation of excessive backlogs in a GI/GI/m queue. *IEEE Transactions on Automatic Control*, 36:1383–1394.

Scandizzo, S. (2016). *The Validation of Risk Models: A Handbook for Practitioners*. Palgrave Macmillan, New York.

Serfling, R. J. (1980). *Approximation Theorems of Mathematical Statistics*. John Wiley and Sons, New York.

Shahabuddin, P. (1994). Importance sampling for highly reliable Markovian systems. *Management Science*, 40:333–352.

Smith, P. J. (2001). Underestimation of rare event probabilities in importance sampling simulations. *SIMULATION*, 76(3):140–150.

Thorisson, H. (1985). The queue GI/G/1: Finite moments of the cycle variables and uniform rates of convergence. *Stochastic Processes and their Applications*, 19(1):85–99.