{'fraction_non_alphanumeric': 0.03921016957822851, 'fraction_numerical': 0.018479678747730006, 'mean_word_length': 4.19326559075513, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 61, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Subsampling to Enhance Efficiency in Input Uncertainty Quantification

In stochastic simulation, input uncertainty refers to the output variability arising from the statistical noise in specifying the input models. This uncertainty can be measured by a variance contribution in the output, which, in the nonparametric setting, is commonly estimated via the bootstrap. However, due to the convolution of the simulation noise and the input noise, the bootstrap consists of a two-layer sampling and typically requires substantial simulation effort. This paper investigates a subsampling framework to reduce the required effort, by leveraging the form of the variance and its estimation error in terms of the data size and the sampling requirement in each layer. We show how the total required effort can be reduced from an order bigger than the data size in the conventional approach to an order independent of the data size in subsampling. We explicitly identify the procedural specifications in our framework that guarantee relative consistency in the estimation, and the corresponding optimal simulation budget allocations. We substantiate our theoretical results with numerical examples.


Introduction
Stochastic simulation is one of the most widely used analytic tools in operations research. It provides a flexible means to approximate complex models and to inform decisions. See, for instance, Law et al. (1991) for applications in manufacturing, revenue management, service and operations systems etc. In practice, the simulation platform relies on input models that are typically observed or calibrated from data. These statistical noises can propagate to the output analysis, leading to significant errors and suboptimal decision-making. In the literature, this problem is commonly known as input uncertainty or extrinsic uncertainty.
In conventional simulation output analysis where the input model is completely pre-specified, the statistical errors come solely from the Monte Carlo noises, and it suffices to account only for such noises in analyzing the output variability. When input uncertainty is present, such an analysis will undermine the actual variability. One common approach to quantify the additional uncertainty is to estimate the variance in the output that is contributed from the input noises (e.g., Song et al. (2014)); for convenience, we call this the input variance. This quantity acts as an uncertainty measure which, when added together with the Monte Carlo variance, gives rise to the overall variance in the outputs. A refined decomposition of input variance across multiple input sources can be used to identify models that are overly ambiguous and flag the need of more data collection (e.g., Song et al. (2014)). Input variance also provides a building block to construct valid output confidence intervals (CIs) that account for combined input and simulation errors (e.g., Cheng and Holland (2004)). Motivated by its central role in quantifying input uncertainty, this paper aims to study the efficient estimation of input variance.
In the literature, bootstrap resampling is a common approach for the above purpose. This applies most prominently in the nonparametric regime, namely when no assumptions are placed on the input parametric family. It could also be used in the parametric case (where more alternatives are available). For example, Cheng and Holland (1997) proposes the variance bootstrap,  studies the consistency of this strategy on a random-effect model that describes the uncertainty propagation. A bottleneck with using bootstrap resampling in estimating input variances, however, is the need to "outwash" the simulation noise, which often places substantial burden on the required simulation effort. More precisely, to handle both the input and the simulation noises, the bootstrap procedure typically comprises a two-layer sampling that first resamples the input data (i.e., outer sampling), followed by running simulation replications using each resample (i.e., inner replications). Due to the reciprocal relation between the magnitude of the input variance and the input data, the input variance becomes increasingly small as the input data size increases. This deems the control of the relative estimation error increasingly expensive, and requires either a large outer bootstrap size or inner replication size to extinguish the effect of simulation noises.
The main goal of this paper is to investigate subsampling as a simulation saver for input variance estimation. This means that, instead of creating distributions by resampling a data set of the full size, we only resample (with or without replacement) a set of smaller size. We show that a judicious use of subsampling can reduce the total simulation effort from an order bigger than the data size in the conventional two-layer bootstrap to an order independent of the data size, while retaining the estimation accuracy. This approach leverages the interplay between the form of the input variance and its estimation error, in terms of the data size and the sampling effort in each layer of the bootstrap. On a high level, the subsample is used to estimate an input variance as if less data are available, followed by a correction of this discrepancy in the data size by properly rescaling the input variance. We call this approach proportionate subsampled variance bootstrap. We explicitly identify the procedural specifications in our approach that guarantee estimation consistency, including the minimally required simulation effort in each layer. We also study the theoretical behavior of our estimation error, in relation to the simulation effort allocation in these layers as well as the input data and subsample sizes, which in turn reveals the optimal configurations and provides implementation guidance.
In the statistics literature, subsampling has been used as a remedy for situations where the fullsize bootstrap does not apply, due to a lack (or undeterminability) of uniform convergence required for its statistical consistency, which relates to the functional smoothness or regularity of the estimators (e.g., Politis and Romano (1994)). Subsampling has been used in time series and dependent data (e.g., Politis et al. (1999), Hall et al. (1995), Datta and McCormick (1995)), extremal estimation (e.g., Bickel and Sakov (2008)), shape-constrained estimation (e.g., Sen et al. (2010)) and other econometric contexts (e.g., Abadie and Imbens (2008), Guggenberger (2009, 2010)). In contrary to these works, our subsampling approach is introduced to reduce the simulation effort faced by the two-layer sampling necessitated from the presence of both the input and simulation noises. In other words, we are not concerned about the issue of uniform convergence, but instead, we aim to distort the relation between the required simulation effort and data size in a way that allows more efficient deconvolution of the effects of the two noises. We also note that, as we will use resampling with replacement (instead of without replacement), our approach is closer to the so-called m out of n bootstrap (Bickel et al. (1997), Bickel and Sakov (2008)). For coherence, throughout the paper we use the term subsampling broadly to indicate a bootstrap with a smaller resample size than the original data size.
We close this introduction with a brief review of other related work in input uncertainty. In the nonparametric regime (the focus of this paper), besides Cheng and Holland (1997) and Song and Nelson (2015) that study bootstrap-based estimation of the input variance, Barton and Schruben (1993) and Barton and Schruben (2001) investigate the percentile bootstrap to construct CIs (i.e., the CI limits are determined from the quantiles of the bootstrap distributions). Like variance bootstrap, percentile bootstrap also encounters two-layer sampling that requires substantial simulation efforts. Yi and Xie (2017) investigates adaptive budget allocation policies based on ranking and selection to reduce simulation cost in the percentile bootstrap, and empirically shows the computational advantage of their approach. On the other hand, contrary to this paper, they do not investigate the required simulation efforts in relation to the input data size. Qian (2016, 2017) study the use of empirical likelihood as an optimization-based alternative to the percentile bootstrap, which requires simulation efforts to estimate the gradient information that remain substantial. Beyond the frequentist regime considered in this paper, Xie et al. (2019) studies nonparametric Bayesian methods based on Dirichlet process mixtures to estimate the variance contributed from input uncertainty and construct CIs. Glasserman and Xu (2014), Hu et al. (2012), Lam (2016b) and Ghosh and Lam (2019) study input uncertainty from a robust optimization viewpoint, where they compute worst-case bounds subject to constraints or so-called uncertainty sets that represent partial beliefs on unknown distributions. In the parametric regime, Barton et al. (2013) and Xie et al. (2016) investigate the basic bootstrap with a metamodel built in advance, a technique known as the metamodel-assisted bootstrap. Xie et al. (2016) and Biller and Corlu (2011) study multivariate input uncertainty assuming a parametric dependency structure in the form of product-moment correlations. Cheng and Holland (1997) studies the delta method, and Holland (1998, 2004) reduce its computation burden via the so-called two-point method. Lin et al. (2015) and Song and Nelson (2019) study regression approaches to estimate sensitivity coefficients which are used to apply the delta method, generalizing the gradient estimation method in Wieland and Schmeiser (2006). Zhu et al. (2020) studies risk criteria and computation to quantify parametric uncertainty. Finally, Chick (2001), Zouaoui and Wilson (2003), Zouaoui and Wilson (2004) and Xie et al. (2014) study variance estimation and interval construction from a Bayesian perspective. We comment that although the exposition in this paper focuses on the nonparametric setting, the same idea of subsampling can be adapted naturally to the parametric setting, with similar advantages in computational efficiency. For general surveys on input uncertainty, readers are referred to Barton et al. (2002), Henderson (2003), Chick (2006), Barton (2012), Song et al. (2014), Lam (2016a), and Nelson (2013) Chapter 7.
The remainder of the paper is as follows. Section 2 introduces the input uncertainty problem and explains the simulation complexity bottleneck in the existing bootstrap schemes. Section 3 presents our subsampling idea, procedures and the main statistical results. Section 4 discusses the key steps in our theoretical developments. Section 5 reports our numerical experiments. Section 6 concludes the paper. All proofs are relegated to the Appendix.

Problem Motivation
This section describes the problem and our motivation. Section 2.1 first describes the input uncertainty problem, Section 2.2 presents the existing bootstrap approach, and Section 2.3 discusses its computational barrier, thus motivating our subsampling investigation. We aim to provide intuitive explanations in this section, and defer mathematical details to later sections.

The Input Uncertainty Problem
Suppose there are m independent input processes driven by input distributions F 1 , F 2 , . . . , F m .
We consider a generic performance measure ψ(F 1 , . . . , F m ) that is simulable, i.e., given the input distributions, independent unbiased replications of ψ can be generated in a computer. As a primary example, think of F 1 and F 2 as the interarrival and service time distributions in a queue, and ψ is some output measure such as the mean queue length averaged over a time horizon. Our study also applies when the F i 's are multivariate distributions.
The input uncertainty problem arises in situations where the input distributions F 1 , . . . , F m are unknown but real-world data are available. One then has to use their estimates F 1 , . . . , F m to drive the simulation. Denote a point estimate of ψ(F 1 , . . . , F m ) asψ( F 1 , . . . , F m ), where typically we takē . . , F m ) being a conditionally unbiased simulation replication driven by F 1 , . . . , F m .
This point estimate is affected by both the input statistical noises and the simulation noises. By conditioning on the estimated input distributions (or viewing the point estimate as a random effect model with uncorrelated input and simulation noises), the variance ofψ( F 1 , . . . , F m ) can be expressed as is the input variance, and is the variance contributed from the simulation noises. Assuming that the estimates F i 's are consistent in estimating F i 's, then, as input data sizes grow, σ 2 S is approximately Var[ψ r (F 1 , . . . , F m )]/q and can be estimated by taking the sample variance of all simulation replications (see, e.g., Cheng and Holland (1997)). On the other hand, σ 2 I signifies the output variance contributed solely from the input data noises, assuming a fully accurate evaluation of the performance measure ψ. Estimating σ 2 I is the key and the challenge in quantifying input uncertainty, which is the focus of this paper.
Before going into details, we discuss two conceptual properties on σ 2 I that would be relevant in motivating and pinpointing our study. First, suppose further that for each input model i, we have . . , X i,n i } generated from the distribution F i . When n i 's are large, typically the overall input variance σ 2 I is decomposable into where σ 2 i /n i is the variance contributed from the data noise for model i, with σ 2 i being a constant. In the parametric case where F i comes from a parametric family containing the estimated parameters, this decomposition is well known from the delta method (Asmussen and Glynn (2007), Chapter 3).
where ∇ i ψ is the collection of sensitivity coefficients, i.e., the gradient, with respect to the parameters in model i, and Σ i is the asymptotic estimation variance of the point estimates of these parameters (scaled reciprocally with n i ). In the nonparametric case where the empirical distribution F i (x) := n i j=1 δ X i,j (x)/n i is used (where δ X i,j denotes the delta measure at X i,j ), (2) still holds under mild conditions (e.g., Propositions 1 and 2 in the sequel). In this setting the quantity σ 2 i is equal to Var F i [g i (X i )], where g i (·) is the influence function (Hampel (1974)) of ψ with respect to the distribution F i , whose domain is the value space of the input variate X i , and Var F i [·] denotes the variance under F i . The influence function can be viewed as a functional derivative taken with respect to the probability distributions F i 's (see Serfling (2009), Chapter 6), and dictates the first-order asymptotic behavior of the plug-in estimate of ψ. Although the mathematical form of σ 2 i 's is known, it relies on gradient information that needs to be estimated via simulation itself. Moreover, in the nonparametric case, the gradient dimension in a sense grows with the data size. Thus directly using the delta method in this case could be challenging. In our subsequent developments, we focus on the nonparametric case, both because this is more challenging, and also that this can be viewed as a generalization of the parametric case by viewing the "parameter" simply as a function of F i 's.
Second, under further regularity conditions, a Gaussian approximation holds forψ( F 1 , . . . , F m ) so thatψ is an asymptotically tight (1 − α)-level CI for ψ(F 1 , . . . , F m ), where z 1−α/2 is the standard normal 1 − α/2 quantile. This CI, which provides a bound-based alternative to quantify input uncertainty, again requires a statistically valid estimate of σ 2 I or m i=1 σ 2 i /n i (and σ 2 S ). In this paper we primarily focus on the estimation of σ 2 I and how our proposed approach substantially improves upon previous methods in this regard. Naturally, the improved estimate of σ 2 I also translates into a better CI when using (3). We caution, however, that an optimal procedural configuration to estimate σ 2 I does not necessarily correspond to an optimal configuration in constructing the CI, as the performance of the latter is measured by different criteria such as coverage or half-width (such a difference in optimally estimating variance versus CI has also been observed in other contexts such as time series (Sun et al. (2008))). Nonetheless, we will show that a direct plug-in of our new estimator of σ 2 I into (3) is already enough to significantly outperform conventional bootstrap-based CIs suggested in the literature, both theoretically and also supported by consistent empirical evidence.
Next we will discuss bootstrap resampling, the commonest estimation technique that forms the basis of our comparison.

Bootstrap Resampling
Let F * i represent the empirical distribution constructed using a bootstrap resample from the orig- is given by To explain, the first term in (4) is an unbiased estimate of the variance ofψ b , which is denotes the expectation on F * i 's conditional on F i 's), sinceψ b incurs both the bootstrap noise and the simulation noise.
More generally, to estimate the variance contribution from the data noise of model i only, namely . . , F m )] by the bootstrap principle, in turn asymptotically equal to σ 2 i /n i introduced in (2). This observation appeared in, e.g., Song et al. (2014); in Section 4 we give further justifications.
In subsequent discussions, we use the following notations. For any sequences a and b, both depending on some parameter, say, n, we say that a = O(b) if |a/b| ≤ C for some constant C > 0 for all sufficiently large n, and a = o(b) if a/b → 0 as n → ∞. Alternately, we say a = Ω(b) if |a/b| ≥ C for some constant C > 0 for all sufficiently large n, and a = ω(b) if |a/b| → ∞ as n → ∞. We say that a = Θ(b) if C ≤ |a/b| ≤ C as n → ∞ for some constants C, C > 0. We use A = O p (b) to represent a random variable A that has stochastic order at least b, i.e., for any > 0, there exists M, N > 0 such that P (|A/b| ≤ M ) > 1 − for n > N . We use A = o p (b) to represent a random variable A that has stochastic order less than b, i.e., A/b p → 0. We use A = Θ p (b) to represent a random variable A that has stochastic order exactly at b, i.e., A satisfies A = O p (b) but not A = o p (b).

A Complexity Barrier
We explain intuitively the total number of simulation runs needed to ensure that the variance bootstrap depicted above can meaningfully estimate the input variance. For convenience, we call this number the simulation complexity. This quantity turns out to be of order bigger than the data size. On a high level, it is because the input variance scales reciprocally with the data size (recall (2)). Thus, when the data size increases, the input variance becomes smaller and increasingly difficult to estimate with controlled relative error. This in turn necessitates the use of more simulation runs.
To explain more concretely, denote n as a scaling of the data size, i.e., we assume n i all grow linearly with n, which in particular implies that σ 2 I is of order 1/n. We analyze the error ofσ 2 BV from Algorithm 1 in estimating σ 2 I . Sinceσ 2 BV is unbiased for Var * [ψ( F * 1 , . . . , F * m )] which is in turn close to σ 2 I , roughly speaking it suffices to focus on the variance ofσ 2 BV . To analyze this later quantity, we denote a generic simulation run in our procedure,ψ r ( F * 1 , . . . , F * m ), aŝ are the errors arising from the bootstrap of the input distributions and the simulation respectively.
If ψ is sufficiently smooth, δ elicits a central limit theorem and is of order Θ p (1/ √ n). On the other hand, the simulation noise ξ is of order Θ p (1).
Via an ANOVA-type analysis as in Sun et al. (2011), we have Now, putting δ = Θ p (1/ √ n) and ξ = Θ p (1) formally into (5), and ignoring constant factors, results in The two terms in (6) correspond to the variances coming from the bootstrap resampling and the simulation runs respectively.
Since σ 2 I is of order 1/n, meaningful estimation of σ 2 I needs measured by the relative error. In other words, we want to achieveσ 2 BV /σ 2 I p → 1 as the simulation budget grows. This property, which we call relative consistency, requiresσ 2 BV to have a variance of order o(1/n 2 ) in order to compensate for the decreasing order of σ 2 I .
We argue that this implies unfortunately that the total number of simulation runs, BR, must be ω(n), i.e., of order higher than the data size. To explain, note that the first term in (6) forces one to use B = ω(1), i.e., the bootstrap size needs to grow with n, an implication that is quite natural. The second term in (6), on the other hand, dictates also that BR 2 = ω(n 2 ). Suppose, for the sake of contradiction, that B and R are chosen such that BR = O(n). Then, because we need BR × R = BR 2 = ω(n 2 ), R must be ω(n) which, combining with B = ω(1), implies that BR = ω(n) and leads to a contradiction.
We summarize the above with the following result. Let N be the total simulation effort, and recall n as the scaling of the data size. We have: Theorem 1 (Simulation complexity of the variance bootstrap). Under Assumptions 1-7 to be stated in Section 4.1, the required simulation budget to achieve relative consistency in estimating σ 2 I by Algorithm 1, i.e.,σ 2 BV /σ 2 Though out of the scope of this paper, there are indications that such a computational barrier occurs in other types of bootstrap. For instance, the percentile bootstrap studied in Barton andSchruben (1993, 2001) appears to also require an inner replication size large enough compared to the data size in order to obtain valid quantile estimates (the authors actually used one inner replication, but Barton (2012) commented that more is needed). Yi and Xie (2017) provides an interesting approach based on ranking and selection to reduce the simulation effort, though they do not investigate the order of the needed effort relative to the data size. The empirical likelihood framework studied in Lam and Qian (2017) requires a similarly higher order of simulation runs to estimate the influence function. Nonetheless, in this paper we focus only on how to reduce computation load in variance estimation.

Procedures and Guarantees in the Subsampling Framework
This section presents our methodologies and results on subsampling. Section 3.1 first explains the rationale and the subsampling procedure. Section 3.2 then presents our main theoretical guarantees, deferring some elaborate developments to Section 4.

Proportionate Subsampled Variance Bootstrap
As explained before, the reason why theσ 2 BV in Algorithm 1 requires a huge simulation effort, as implied by its variance (6), lies in the small scale of the input variance. In general, in order to estimate a quantity that is of order 1/n, one must use a sample size more than n so that the estimation error relatively vanishes. This requirement manifests in the inner replication size in constructingσ 2 BV .
To reduce the inner replication size, we leverage the relation between the form of the input variance and the estimation variance depicted in (6) as follows. The approximate input variance contributed from model i, with data size n i , has the form σ 2 i /n i . If we use the variance bootstrap directly as in Algorithm 1, then we need an order more than n total simulation runs due to (6). Now, pretend that we have fewer data, say s i , then the input variance will be σ 2 i /s i , and the required simulation runs is now only of order higher than s i . An estimate of σ 2 i /s i , however, already gives us enough information in estimating σ 2 i /n i , because we can rescale our estimate of σ 2 i /s i by s i /n i to get an estimate of σ 2 i /n i . Estimating σ 2 i /s i can be done by subsampling the input distribution with size s i . With this, we can both use fewer simulation runs and also retain correct estimation via multiplying by a s i /n i factor.
To make the above argument more transparent, the bootstrap principle and the asymptotic approximation of the input variance imply that The subsampling approach builds on the observation that a similar relation holds for where F * s i ,i denotes a bootstrapped input distribution of size s i (i.e., an empirical distribution of size s i that is uniformly sampled with replacement from {X i,1 , . . . , X i,n i }). If we let s i = θn i for some θ > 0 so that s i → ∞ (where · is the floor function, i.e. the largest integer less than or equal to ·), then we have Multiplying both sides with θ, we get Note that the right hand side above is the original input variance of interest. This leads to our proportionate subsampled variance bootstrap: We repeatedly subsample collections of input distributions from the data, with size θn i for model i, and use them to drive simulation replications.
We then apply the ANOVA-based estimator in (4) on these replications, and multiply it by a factor of θ to obtain our final estimate. We summarize this procedure in Algorithm 2. The term "proportionate" refers to the fact that we scale the subsample size for all models with a single factor θ. For convenience, we call θ the subsample ratio.
say F * s i ,i whose uncertainty is of interest, while fixing all the other distributions as the original empirical distributions, i.e., F j , j = i. All the remaining steps in Algorithm 2 remain the same (thus the "proportionate" part can be dropped). This procedure is depicted in Algorithm 3.

Statistical Guarantees
Algorithm 2 provides the following guarantees. Recall that N = BR is the total simulation effort, and n is the scaling of the data size. We have the following result: Theorem 2 (Procedural configurations to achieve relative consistency). Under Assumptions 1-7 to be stated in Section 4.1, if the parameters B, R, θ of Algorithm 2 are chosen such that then the variance estimateσ 2 SV B is relatively consistent, i.e.σ 2 SV B /σ 2 I p → 1.

Algorithm 3 Subsampled Variance Bootstrap for Variance Contribution from the i-th Input
Model Theorem 2 tells us what orders of the bootstrap size B, inner replication size R and subsample ratio θ would guarantee a meaningful estimation of σ 2 I . Note that θ ≈ s i /n i for each i, so that θ = ω(1/n) is equivalent to setting the subsample size s i = ω(1). In other words, we need the natural requirement that the subsample size grows with the data size, albeit can have an arbitrary rate.
Given a subsample ratio θ specified according to (7), the configurations of B and R under (7) that achieve the minimum overall simulation budget is B = ω(1) and R = Ω(θn). This is because to minimize N = BR while satisfying the second requirement in (7), it is more economical to allocate as much budget to R instead of B. This is stated precisely as: Corollary 1 (Minimum configurations to achieve relative consistency). Under the same conditions of Theorem 2, given θ = ω(n −1 ), the values of B and R to achieve (7) and hence relative consistency that requires the least order of effort are B = ω(1) and R = Ω(θn), leading to a total simulation budget N = ω(θn).
Note that θn is the order of the subsample size. Thus Corollary 1 implies that the required simulation budget must be of higher order than the subsample size. However, since the subsample size can be chosen to grow at an arbitrarily small rate, this implies that the total budget can also grow arbitrarily slowly. Therefore, we have: Corollary 2 (Simulation complexity of proportionate subsampled variance bootstrap).
Under the same conditions of Theorem 2, the minimum required simulation budget to achieve relative consistency in estimating σ 2 I by Algorithm 2, i.e.,σ 2 SV B /σ 2 Compared to Theorem 1, Corollary 2 stipulates that our subsampling approach reduces the required simulation effort from a higher order than n to an arbitrary order, i.e., independent of the data size. This is achieved by using a subsample size that grows with n at an arbitrary order, or equivalently a subsample ratio θ that grows faster than 1/n.
The following result describes the configurations of our scheme when a certain total simulation effort is given. In particular, it shows, for a given total simulation effort, the range of subsample ratio for which Algorithm 2 can possibly generate valid variance estimates by appropriately choosing B and R: Theorem 3 (Valid subsample ratio given total budget). Assume the same conditions of Theorem 2. Given a total simulation budget N = ω(1), if the subsample ratio satisfies ω(1/n) ≤ θ ≤ o(N/n) ∧ 1, then the bootstrap size B and the inner replication size R can be appropriately chosen according to criterion (7) to achieve relative consistency, i.e.,σ 2 SV B /σ 2 The next result is on the optimal configurations of our scheme in minimizing the Monte Carlo error. To proceed, define as the perfect form of our proportionate subsampled variance bootstrap introduced in Section 3.1, namely without any Monte Carlo noises, and 0 < θ ≤ 1 is the subsample ratio. We have: Theorem 4. Assume the same conditions of Theorem 2. Given a simulation budget N and a subsample ratio θ such that N = ω(θn) and θ = ω(n −1 ), the optimal outer and inner sizes that minimize the order of the conditional mean squared error Note that the mean squared error, i.e.
random because the underlying resampling is conditioned on the input data, therefore the bound at the end of Theorem 4 contains a stochastically vanishing term o p (1).
We next present the optimal tuning of the subsample ratio. This requires a balance of the tradeoff between the input statistical error and the Monte Carlo simulation error. To explain, the overall error ofσ 2 SV B by Algorithm 2 can be decomposed aŝ The first term is the Monte Carlo error for which the optimal outer size B, inner size R and the resulting mean squared error are governed by Theorem 4. In particular, the mean squared error there shows that under a fixed simulation budget N and the optimal allocation R = Θ(θn), the Monte Carlo error gets larger as θ increases. The second term is the statistical errors due to the finiteness of input data and θ. Since θ measures the amount of data contained in the resamples, we expect this second error to become smaller as θ increases. The optimal tuning of θ relies on balancing such a trade-off between the two sources of errors.
We have the following optimal configurations of B, R and θ altogether given a budget N : Theorem 5 (Optimal subsample size and budget allocation). Suppose Assumptions 1, 3-7 in Section 4.1 and Assumptions 10-12 in Section 4.3 hold. For a given simulation budget N = ω(1), if the subsample ratio θ and outer and inner sizes B, R for Algorithm 2 are set to then the gross errorσ 2 SV B − σ 2 I = E + o p (N −1/3 n −1 + n −3/2 ), where the leading term has a mean squared error Moreover, if R = Θ((ns) −1 ) and at least one of the Σ i 's are positive definite, where R and Σ i are as defined in Lemma 2, then (12) holds with an exact order (i.e., O(·) becomes Θ(·)) and the configuration (10), (11) is optimal in the sense that no configuration gives rise to a gross error Note from (12) that, if the budget N = ω(1), our optimal configurations guarantee the estimation mean squared error decays faster than 1/n 2 . Recall that the input variance is of order 1/n, and thus an estimation error of order higher than 1/n 2 ensures that the estimator is relatively consistent in the senseσ 2 SV B /σ 2 I p → 1. This recovers the result in Corollary 2. We also comment that the algorithmic configuration given in Theorem 5 is chosen to optimize the mean squared error of the input variance estimate, but does not necessarily generates the most accurate CI. There exists evidence (e.g., Sun et al. (2008)) that the optimal choice to minimize the mean squared error of the variance estimate can be different from the one that is optimal for statistical inference, although in our experiments they seem to match closely with each other.
We comment that all the results in this section hold if one estimates the individual variance contribution from each input model i, namely by using Algorithm 3. In this case we are interested in estimating the variance σ 2 i /n i , and relative consistency meansσ 2 SV B,i /(σ 2 i /n i ) p → 1. The data size scaling parameter n can be replaced by n i in all our results.
Finally, we also comment that the complexity barrier described in Section 2.3 and our framework presented in this section applies in principle to the parametric regime, i.e., when the input distributions are known to lie in parametric families with unknown parameters. The assumptions and mathematical details would need to be catered to that situation, which could be done naturally by viewing the "parameter" as a function of F i 's.

Developments of Theoretical Results
We present our main developments leading to the algorithms and results in Section 3. Section 4.1 first states in detail our assumptions on the performance measure. Section 4.2 presents the theories leading to estimation accuracy, simulation complexity and optimal budget allocation in the proportionate subsampled variance bootstrap. Section 4.3 investigates optimal subsample sizes that lead to overall best configurations.

Regularity Assumptions
We first assume that the data sets for all input models are of comparable size.
Recall in Sections 2 and 3 that we have denoted n as a scaling of the data size. More concretely, we take n = (1/m) m i=1 n i as the average input data size under Assumption 1.
We next state a series of general assumptions on the performance measure ψ. These assumptions hold for common finite-horizon measures, as we will present. For each i let Ξ i be the support of the i-th true input model F i , and the collection of distributions P i be the convex hull spanned by F i and all Dirac measures on Ξ i , i.e.
We assume the following differentiability of the performance measure.

Assumption 2 (First order differentiability). For any distributions
. . , m and as all ν i 's approach zero The differentiability described above is defined with respect to a particular direction, namely Q i −P i , in the space of probability measures, and is known as Gateaux differentiability or directional differentiability (e.g., Serfling (2009), Van der Vaart (2000)). Assumption 2 therefore requires the performance measure ψ to be Gateaux differentiable when restricted to the convex set P 1 ×· · ·×P m .
The functions g i 's are also called the influence functions (e.g., Hampel (1974)) that play analogous roles as standard gradients in the Euclidean space. The condition of g i 's having vanishing means is without loss of generality since such a condition can always be achieved by centering, i.e., subtracting the mean. Note that doing this does not make any difference to the first term of expansion (13) because both Q i and P i are probability measures. Taking each ν i = 1 in (13), one informally obtains the Taylor expansion of ψ around P i 's When each P i is set to be the true input model F i and Q i to be the empirical input model F i , the above linear expansion is expected to be a reasonably good approximation as the data size grows.
The next assumption imposes a moment bound on the error of this approximation: Assumption 3 (Smoothness at true input models). Denote by g i (·) := g i (F 1 , . . . , F m ; ·) the influence functions at the true input distributions F i , i = 1, . . . , m. Assume that the remainder in the Taylor expansion of the performance measure satisfies E[ 2 ] = o(n −1 ), and the influence functions g i 's are non-degenerate, i.e. σ 2 i := Assumption 3 entails that the error of the linear approximation formed by influence functions is negligible in the asymptotic sense. Indeed, the linear term in (14) is asymptotically of order Θ p (n −1/2 ) by the central limit theorem, whereas the error is implied by Assumption 3 to be o p (n −1/2 ). Hence the variance of the linear term contributes dominantly to the overall input variance as n i 's are large. Then, thanks to the independence among the input models, the input variance can be expressed in the additive form described in (2) together with a negligible error.
Proposition 1. Under Assumptions 1-3, the input variance σ 2 I defined in (1) takes the form is the variance of the i-th influence function.
As the higher order o 1/n error suggests, the additive decomposition m i=1 σ 2 i n i is guaranteed to be accurate only in the large-sample regime. Note that this decomposition is used solely as a theoretical vehicle for asymptotic analysis rather than the actual input variance estimator in our procedure, the latter using bootstrapping schemes that could exhibit better finite-sample performances.
As mentioned before, consistent estimation of input variance σ 2 I relies on the bootstrap principle, for which we make the following additional assumptions. The assumption states that the error of the linear approximation (14) As the data sizes n i 's grow, the empirical input distributions F i converge to the true ones F i . Hence the empirical influence functionsĝ i 's are expected to approach the influence functions g i 's associated with the true input distributions, which explains the convergence condition in Assumption 4. The fourth moment condition on the remainder * is needed for controlling the variance of our variance estimator. Since the fourth moment is with respect to the resampling measure and thus depends on the underlying input data, the condition is described in terms of stochastic order. Note that we require (15) to hold not just when F i = F * s i ,i for all i but also when some F i = F i . This allows us to estimate the variance contributed from an arbitrary group of input models and in particular an individual input model.
Assumptions 2-4 are on the performance measure ψ itself. Next we impose assumptions on the simulation noise, i.e. the stochastic errorψ r − ψ whereψ r is an unbiased simulation replication for ψ. We denote by τ 2 (P 1 , . . . , P m ) the variance ofψ r when simulation is driven by arbitrary input models P 1 , . . . , P m , i.e.
Assumption 6 (Convergence of bootstrapped variance). For every Assumption 7 (Boundedness of the fourth moment). For every Assumptions 5 and 6 stipulate that the variance of the simulation replicationψ r as a functional of the underlying input models is smooth enough in the inputs. Conceptually Assumption 5 is in line with Assumption 3 in the sense that both concern smoothness of a functional around the true input models, whereas Assumption 6 is similar to Assumption 4 since both are about smoothness property around the empirical input models. Assumption 7 is a fourth moment condition like in Assumption 4 used to control the variance of the variance estimator. Similar to Assumption 4, we impose Assumptions 6 and 7 for each F i = F i or F * s i ,i so that the same guarantees remain valid when estimating input variances from individual input models, i.e., Algorithm 3.
Although the above assumptions may look complicated, they can be verified, under minimal conditions, for generic finite-horizon performance measures in the form distributed under F i , each T i being a deterministic time, and h is a performance function. An unbiased simulation replicationψ r of the performance measure is h(X 1 , . . . , X m ).
Suppose we have the following conditions for the performance function h: Assumption 9 (Parameter k). For each i let I i = (I i (1), . . . , I i (T i )) be a sequence of indices such The conditional expectation in Assumption 8 is in fact the influence function of the performance measure (16) under the true input models. So Assumption 8 is precisely the non-degenerate variance condition in Assumption 3. All other parts of Assumptions 2-7 are consequences of the moment condition in Assumption 9: Theorem 6. Under Assumptions 1, 8 and Assumption 9 with k = 4, we have Assumptions 2-7 hold for the finite-horizon performance measure ψ given by (16).

Simulation Complexity and Allocation
This section presents theoretical developments on our proportionate subsampled variance bootstrap. We first establish relative consistency assuming infinite computation resources. Recall (8) as the proportionate subsampled variance bootstrap estimator without any Monte Carlo errors.
The following theorem gives a formal statement on the performance of this estimator discussed in Section 3.1.
Theorem 7. Under Assumptions 1-4, if the subsample ratio θ = ω(n −1 ), then the proportionate subsampled variance bootstrap without Monte Carlo error, namely (8), is relatively consistent as The requirement θ = ω(n −1 ) implies that s i → ∞, which is natural as one needs minimally an increasing subsample size to ensure the consistency of our estimator. It turns out that this minimal requirement is enough to ensure consistency even relative to the magnitude of σ 2 I .
and R is arbitrary, then the conditional mean squared error In addition to the condition θ = ω(n −1 ) which has appeared in Theorem 7, we also require B = ω(1) in Lemma 1. As the proof reveals, with such a choice of B, we can extract the leading term of the conditional mean squared error shown in (18), which takes a neat form and is easy to analyze.
Note that σ 2 I here is of order n −1 by Proposition 1. Hence the Monte Carlo noise of the variance estimate output by our algorithm has to vanish faster than n −1 in order to achieve relative consistency. Combining Theorem 7 and Lemma 1, we obtain the simulation complexity ofσ 2 SV B in Theorem 2. To establish the theoretical optimal allocation on the outer and inner sizes B, R, for given data sizes n i , subsample ratio θ, and total simulation budget N , we minimize the conditional mean square error (18) subject to the budget constraint BR = N . This gives rise to the following result that gives a more precise (theoretical) statement than Theorem 4.
Theorem 8. Suppose Assumptions 1-7 hold. Given a simulation budget N and a subsample ratio θ such that N = ω(θn) and θ = ω(n −1 ), the optimal outer and inner sizes that minimize the condi- Theorem 8 gives the exact choices of B and R that minimize the Monte Carlo error. However, this is more of theoretical interest because the optimal R * involves the desired input variance m i=1 σ 2 i /n i . Having said that, we can conclude from the theorem that the optimal inner size R is of order Θ(θn), the same as the subsample size, because the input variance is of order Θ(1/n) by Proposition 1 and τ 2 is a constant. This results in Theorem 4 in Section 3.2.

Optimal Subsample Ratio
In this section we further establish the optimal subsample ratio θ or equivalently subsample sizes s i that balance the two sources of errors in (9). For this, we need more regularity conditions on the performance measure. The first assumption we need is third order Gateaux differentiability in the convex set P 1 × · · · × P m : Assumption 10 (Third order differentiability). Using the same notations P i , Q i , P ν i i as in Assumption 2, assume that there exist second order influence functions g i 1 i 2 (P 1 , . . . , P m ; ·) : Ξ i 1 × Ξ i 2 → R and third order influence functions g i 1 i 2 i 3 (P 1 , . . . , P m ; ·) : Ξ i 1 × Ξ i 2 × Ξ i 3 → R for i 1 , i 2 , i 3 = 1, . . . , m which are symmetric under permutations, namely g i 1 i 2 (P 1 , . . . , P m ; x 1 , x 2 ) = g i 2 i 1 (P 1 , . . . , P m ; x 2 , x 1 ) and for all x, y satisfy Moreover, as all ν i 's approach zero the following Taylor expansion holds Assumption 10 complements and strengthens Assumption 2 in that it imposes stronger differentiability property. Similarly, the following two assumptions strengthen Assumptions 3 and 4 respectively by considering cubic expansions.
Assumption 11 (Third order smoothness at true input models). Denote by g i 1 i 2 (·) := g i 1 i 2 (F 1 , . . . , F m ; ·) and g i 1 i 2 i 3 (·) := g i 1 i 2 i 3 (F 1 , . . . , F m ; ·) the second and third order influence functions under the true input models. Assume the remainder in the Taylor expansion of the plug-in , and the high order influence functions satisfy the moment conditions for all i 1 , i 2 , i 3 and j 2 ≤ 2, j 3 ≤ 3, where X i,j is the j-th data point from the i-th input model.
Similar to the remainder in Assumption 3, the moment condition on 3 here is used to control the error of the cubic approximation of ψ formed by up to third order influence functions. With these additional assumptions, the error term in Proposition 1 can be refined as follows: Proposition 2. Under Assumptions 1, 3 and 10-11, the overall input variance, as defined in (1), can be expressed as We also need third order differentiability around the empirical input models: Assumption 12 (Third order smoothness at empirical input models). Denote bŷ . . , F m ; ·) the second and third order influence functions under the empirical input models. Assume that the remainder in the Taylor expansion of the bootstrapped performance measure ψ( F * s 1 ,1 , . . . , F * sm,m ) . In addition, assume the high order empirical influence functionsĝ i 1 i 2 andĝ i 1 i 2 i 3 converge in mean square error, i.e.
for all i 1 , i 2 , i 3 and j 2 ≤ 2, j 3 ≤ 3, where X i,j is the j-th data point from the i-th input model. For the first order influence functionĝ i , assume the remainder in the Taylor expansion As for Assumptions 3 and 4, finite-horizon performance measures under mild conditions satisfy the above two assumptions: Theorem 9. Under Assumptions 1, 8 and Assumption 9 with k = 4, we have Assumptions 10-12 hold for the finite-horizon performance measure ψ given by (16).
With Assumptions 11 and 12, we can identify the statistical error of our variance estimator assuming infinite computation resources, which we summarize in the following lemma.
Lemma 2. Under Assumptions 1, 3-4 and 10-12, the statistical error of the proportionate subsampled bootstrap variance is characterized by where Z is a random variable such that R is defined as where frac(x) := x − x denotes the fraction part of x ∈ R, and for each i, X i , X i are independent copies of the random variable distributed under F i .
Combining the statistical error (20), and the minimal Monte Carlo error (19) under the optimal budget allocation into the trade-off (9), we obtain the overall error of the outputσ 2 SV B of Algorithm 2: Theorem 10 (Overall error of the variance estimate). Suppose Assumptions 1, 3-7 and 10-12 hold. Given a simulation budget N and a subsample ratio θ such that N = ω(θn) and θ = ω(n −1 ), if outer and inner sizes B, R for Algorithm 2 are chosen to be R = Θ(θn), B = N/R, then the gross error of our Monte Carlo estimateσ 2 SV B − σ 2 I = E + o p (θ 1/2 (N n) −1/2 + θ −1 n −2 + n −3/2 ), where the leading term has a mean squared error where R, λ i 's and Σ i 's are defined in Lemma 2.
It is clear from their definitions in Lemma 2 that R = O(θ −1 n −2 ) and each (λ hence the mean squared error (21) is in general of order O(θ(N n) −1 + θ −2 n −4 + n −3 ). When R and at least one of the λ T i Σ i λ i 's satisfy the non-degeneracy condition in Theorem 5, this bound becomes tight in order, and the optimal subsample ratio can be established by minimizing the order of the leading overall error E.

Numerical Experiments
This section reports our experimental findings. We consider two examples with different scales and complexities: M/M/1 queue: The first example we consider is an M/M/1 queue that has true arrival rate 0.5 and service rate 1. Suppose the system is empty at time zero. The performance measure of interest is the probability that the waiting time of the 20-th arrival exceeds 2 units of time, whose true value is approximately 0.182. Specifically, the system has two input distributions, i.e., the inter-arrival time distribution F 1 = Exp(0.5) and the service time distribution F 2 = Exp(1), for which we have n 1 and n 2 i.i.d. data available respectively. If A t is the inter-arrival time between the t-th and (t + 1)-th arrivals, and S t is the service time for the t-th arrival, then the system output where the waiting time W 20 is calculated by the Lindley recursion W t+1 = max{W t + S t − A t , 0} for t = 1, . . . , 19 and W 1 = 0. To test the proposed approach under different levels of utilization, we also consider true arrival rate 0.9 and service rate 1, for which case the target performance measure is taken to be the probability that the waiting time of the 20-th arrival exceeds 6 units of time (true value 0.190). The data sizes n 1 , n 2 are chosen so that n 1 = 2n 2 in the experiments, so only the minimum min i n i is reported for convenience.
Computer network: We also consider a computer communication network borrowed from Cheng and Holland (1997) and Lin et al. (2015). The structure of the system is characterized by the undirected graph in Figure 1: Four message-processing units, which correspond to the nodes, are connected by four transport channels that are represented by the edges. For every pair i, j  Table 1 True arrival rates λi,j of messages to be transmitted from node i to node j.
0.001 seconds to process a message, and has unlimited storage capacity. The messages have lengths that are independent and follow an exponential distribution with mean 300 bits, and each channel has a capacity of 275000 bits, therefore there are queuing and transmission delays. The messages travel through the channels with a velocity of 150000 miles per second, and the i-th channel has a length of 100 · i miles for i = 1, 2, 3, 4, leading to a propagation delay of 100·i 150000 seconds along the i-th channel. The total time that a message of length l bits occupies the i-th channel is therefore l 275000 + 100·i 150000 seconds. Suppose the system is empty at time zero. The performance measure of interest is the average delay of the first 30 messages that arrive to the system, or mathematically, where D k is the time for the k-th message to be transmitted from its entering node to destination node. The true value of the performance measure is approximately 6.91 × 10 −3 seconds. In the experiment, we assume that the arrival rates of the different types of messages, as well as the distribution of the message length, are unknown, therefore there are 13 input models in total. Like in the example of M/M/1 queue, the data sizes across different input models are kept proportional to each other and only the minimum size is reported.
In the experiments we investigate the simulation efforts needed for our subsampling procedure to generate accurate estimates of the input variance, the impacts of the procedural parameters θ, B, R on the estimation accuracy, and practical guidelines on optimal choices of these parameters.
Regarding performance metrics of the method, we primarily focus on the mean squared error of the obtained input variance estimate. In addition, note that our estimated input variance can also be used to construct CIs by plugging into formula (3). We also examine the quality of these CIs, measured by coverage accuracy and width, as impacted by the estimation accuracy of the input variance.
We compare our subsampling approach with the variance bootstrap depicted in Algorithm 1 and the percentile bootstrap suggested by Barton andSchruben (1993, 2001). In converting our subsampled input variance estimate to CI, we also investigate the use of a "splitting" versus a "non-splitting" approach. In most part of this section, we use the splitting approach that divides the budget into two portions with one used to estimate the input variance and the other to compute the point estimator. To describe it in detail, suppose we have a total budget of N simulation runs. We allocate R v simulation runs to estimate σ 2 I using either Algorithm 1 or 2, and the remaining R e = N − R v simulation runs driven by the empirical input distributions to compute the point estimatorψ( F 1 , . . . , F m ). When constructing the CI in (3) accurately approximate the plug-in estimator ψ( F 1 , . . . , F m ) with an error that is negligible relative to the input variability. Using the former as a surrogate for the latter avoids splitting the budget; however, we will see later that this may introduce too much bias to maintain the desired coverage level when the subsample size is relatively small.
The rest of this section is organized as follows. Section 5.1 investigates practical guidelines for choosing the algorithmic parameters in our procedure. Using these guidelines, in Section 5.2 we compare the proposed procedure with the variance bootstrap and the percentile bootstrap. Section 5.3 studies further the conversion of input variance estimate into CI, and compares the associated splitting and non-splitting approaches.

Guidelines for Algorithmic Configuration
We examine the performances using a wide range of parameter choices for θ, B, R.  We first study and establish guidelines for the outer size B and inner size R for a given subsample size. Figure 2 shows how the estimation error changes as the inner replication size R grows from 5 to 40 (correspondingly the outer size B drops from 200 to 25) and the subsample size θ min i n i is fixed at a certain value. Each curve represents the results for one of the considered examples under a particular input data size. Although the precise optimal choice for B, R varies from one example to another even when the subsample size is chosen the same, the estimation error appears robust to the parameter choices, with a range of values that only slightly underperform the optimal. In particular, compared to the unknown optimal choice, an R between 1 6 θ min i n i and 1 3 θ min i n i seems to achieve a comparable accuracy level in the variance estimation, hence is recommended as a general choice. Now we turn to optimal choices for the subsample size. Provided that B, R is properly chosen as above, we examine the behavior of the variance estimation error as the subsample size varies. As we have discussed in Section 3.1, subsampling is preferred when the input data size is relatively large, and thus we consider input data sizes ≥ 500 for our M/M/1 queue and computer network, and for each considered data size we plot the variance estimation error versus the subsample size in Figure 3. We see that a too large size such as 120 always leads to a larger estimation error than moderate sizes like 30, whereas a too small size around 5 can lift the error by even more in some cases, which is consistent with the theoretical insight from the bound (21). Therefore, in general we recommend the use of a subsample size θ min i n i between 20 and 40 to optimize the estimation accuracy. Figure 3 shows that, under the suggested subsample size, the relative rmse is as low as 0.2-0.5 across all the cases.

Comparisons with the Variance Bootstrap and the Percentile Bootstrap
We compare our subsampling method with the standard variance bootstrap and the percentile bootstrap, under the same total budget of 1500 simulation runs. In addition to the relative rmse of the input variance estimate, we also report the actual coverage probability and width of the CI constructed by plugging in the input variance estimate. To estimate all these performance metrics, we construct 1000 95%-level CIs for the target performance measures, each from an independently generated input data set. The "splitting" approach that splits the total budget into R v = 1000, R e = 500 is adopted for the subsampling approach and the variance bootstrap, whereas for the percentile bootstrap all the 1500 simulation runs are used for the resamples. As suggested in Section 5.1, we use the parameter values θ = 30 min i n i , B = 100, R = 10 in our method in all the cases, whereas for the other two methods we vary the parameter configurations over a reasonable range constrained by the simulation budget and then report the best results generated by these considered configurations.  Table 4 shows those for the computer network. The shorthand "PSVB" stands for proportionate subsampled variance bootstrap, i.e., our subsampling approach.
For each method, the "coverage estimate" column displays estimates of the actual coverage probability based on 1000 independent CIs, and the "CI width" column shows their average width.
The second column of each table shows the ratio between the input standard error σ I and the simulation standard error σ S for different input data sizes in our "splitting" approach. A ratio close to or greater than 1 means that the input noise is a major source of uncertainty relative to the simulation noise, thus indicating the need to be taken into account in output analysis.
We compare the approaches based on Tables 2-4 Table 3 Results for the M/M/1 queue with arrival rate 0.9 and service rate 1.
estimates generated by our approach have a smaller relative error than those by the variance bootstrap in all considered cases, and the gap becomes more significant as the data size grows larger. In particular, as the data size grows from 30 to thousands, the estimation error keeps decreasing from 0.7 to 0.25 in our approach, whereas in variance bootstrap it keeps increasing from 0.7 to larger than 1, a level that makes the estimate too crude to be useful. These demonstrate the computational advantage and dictate the use of subsampling especially when the input data size is relatively large. Note that the same budget of 1000 simulation runs are used in input variance  estimation for all considered data sizes and that the estimation accuracy seems much better for large data sizes than for small sizes, and one may wonder whether more simulation runs should be used for small data sizes to further improve the estimation accuracy. It turns out that the estimation errors are mostly due to the inadequacy of the input data rather than the simulation budget, hence a budget of 1000 is already large enough and further increasing the budget does not bring much benefit. For instance, in the case of data size 30 in Table 2, the relative error of the input variance estimate remains as large as 0.69 even if the simulation budget is increased by 10 times.
Secondly, thanks to the high accuracy in the input variance estimates, our subsampling approach generates accurate CIs whose coverage probabilities quickly approach the nominal level 95% as the input data size grows. In contrast, the CIs using the variance bootstrap exhibit under-coverage, and the percentile bootstrap CIs significantly over-cover the truth. We see that the coverage of the variance bootstrap is below 90% in most considered cases, and in the very few cases where the CIs happen to have relatively good coverages, the intervals are much wider than those by our subsampling approach. For example, in the case of data size 2000 in Table 2, the variance bootstrap gives a fairly accurate coverage 92.8%, but on average the interval is 1.23 times as wide as that by our method. This shows that the better estimates of the input variance using subsampling translate to better CIs significantly compared to using the variance bootstrap, in terms of both coverage accuracy and width. The percentile bootstrap CIs show an overly high coverage probability close to 100% and are 2-3 times wider than those by subsampling for all considered input data sizes except 30. The over-coverage issue of the percentile CIs arises because the order statistics capture only the input noise but not the simulation noise in the resampled performance measures, a phenomenon that has been discussed in Barton et al. (2007Barton et al. ( , 2018. When one can afford a sufficiently large budget of simulation relative to the input data size, the simulation noise can be made negligible so that the CIs have the correct coverage. However, when simulation resources are relatively limited (e.g., when data size ≥ 100 in Tables 2-4), the CIs are unnecessarily widened by the extra simulation noise that leads to over-coverage. We also notice that the percentile bootstrap CIs do show more accurate coverage than the other two methods when the input data size is 30, which may suggest that the percentile bootstrap is the preferred approach to constructing CIs in small data cases.
However, this outperformance is a result of optimally choosing the parameters B, R in hindsight. In our experiments, this best parameter set varies from one case to another, and the actual coverage under different configurations varies in a range of 8%.
Thirdly, results across different input data sizes show that, the advantages of subsampling in both input variance estimation and CI construction are most significant in situations with relatively large input data size. Note that one may argue in such situations input uncertainty is negligible.
However, whether this is indeed the case relates to the error tolerance of the decision-maker and the magnitude of the target performance measure itself. For the large data sizes we consider, the input noise appears still relatively substantial. For instance, when the input data size is 2000 in Table   3, the average width of the CIs as a measure of the input uncertainty and simulation uncertainty combined amounts to as much as 57% of the target tail probability, and that the input uncertainty serves as a major component of the total uncertainty (a ratio of 1.3 relative to the simulation uncertainty).
Lastly, in situations with small input data size like 30 the CI coverage clearly falls below 95% in Tables 2 and 3. This under-coverage phenomenon may appear to stem from the nonlinear effect of the performance measure that is inadequately captured by the Gaussian-approximation-based CI given in (3). The real reason, as our experiments suggest, turns out to be the insufficient accuracy of the input variance estimates. In fact, if the true input variance σ 2 I (which can be accurately estimated by repeatedly generating independent input data sets) is plugged into (3) to construct CIs, the coverage probability under the data size 30 rises to 94%-95% for both the M/M/1 queue and the computer network. This indicates a positive impact of an accurate input variance estimate on the CI quality, a point that we will discuss further momentarily.

Approaches
We study in more depth the relation between the input variance estimation accuracy and CI quality, and compare the splitting approach for CI construction that has been used in previous subsections, with the alternate non-splitting approach described at the beginning of this section. Finally, we provide practical budget allocation strategies for the splitting approach.
First, to see how the estimation accuracy of the input variance affects the coverage accuracy of the CIs, we use the splitting approach to compute 95%-level CIs, with 1000 simulation runs assigned to input variance estimation and another 500 runs to point estimator evaluation. Figure   4a plots the coverage probability versus the relative rmse when the subsample size θ min i n i is chosen 30 in the M/M/1 queue example, where each point corresponds to a particular combination of the data size min i n i , the outer replication size B, and the inner replication size R. Figure 4b plots the same for the computer network example with subsample size 30. Both figures clearly show that, the more accurately the input variance is estimated, the closer to the nominal level 95% the coverage probability will be. Accurate estimation of the input variance thus appears to play a crucial role in the construction of accurate CIs.
Next we compare the splitting and non-splitting approaches under the same total budget of 1500 simulation runs. Like in the splitting approach, we use a subsample size θ min i n i = 30 for our non-splitting approach, but use B = 75, B = 20 to consume all the 1500 simulation runs. We  Monotonicity between coverage accuracy and input variance estimation accuracy.
find that the CIs generated from the two approaches have similar lengths, but the non-splitting approach underperforms in terms of coverage accuracy. Each plot in Figure 5 shows the coverage probabilities of the non-splitting CIs versus the splitting ones for each of the considered example systems, as the input data size grows from 30 to thousands. We see that when the data size is relatively small (e.g., below 500), the two approaches generate CIs with similar coverage accuracy.
When the data size grows larger, however, the coverage probability of the non-splitting CIs keeps dropping in all the three examples, especially in the M/M/1 queue with arrival rate 0.9 where a drop towards 86% is observed, whereas the splitting CIs exhibit almost exact 95% coverage.
A possible cause of the undercoverage is the overly small subsample size compared to the input data size, which leads to a high bias in the point estimator. With a subsample size s, the bias of the non-splitting point estimator E * [ψ] with respect to the truth ψ(F 1 , . . . , F m ) can be as large as O(1/s). Given that the input standard error is Θ(1/ √ n), E * [ψ] has a negligible bias only when the subsample size is large enough, namely when s = ω( √ n), indicating that a small subsample size relative to the data size can corrupt the CI. In our experiment, we find that the (supposedly unobservable) bias can be as large as 25% of the CI width when the input data size is 2000 in the M/M/1 queue with arrival rate 0.9, and that artificially removing the bias from the point estimator can improve the coverage to a similar level achieved by the splitting approach. Because of the bias and the consequent under-coverage issue, we caution the use of the non-splitting approach, that it should only be used when a relatively large subsample size is adopted.  Since the splitting approach is recommended, next we explore strategies of splitting a given budget. Our goal is to generate shortest possible CIs that have a sufficiently accurate coverage probability. As in the beginning of the section, denote by R v the number of simulation runs used to estimate the input variance, and by R e to construct the point estimator. Under a fixed total budget R v + R e = 1500, we try four different splits R v = 100, 250, 500, 1000 (accordingly R e = 1500 − R v ), and for each split the subsample size is fixed at θ min i n i = 30 and several choices of B, R are tested among which the one with the best coverage probability is reported. Figure 6 plots the coverage probability versus the CI width for the four considered splits, where the M/M/1 queue with arrival rate 0.9 is considered and input data size is 2000. We notice that the split controls a tradeoff between the coverage accuracy and the CI width. The more simulation runs one allocates to input variance estimation, the more accurate but wider CIs one would obtain, because the input variance is more accurately estimated while the point estimator becomes more noisy. The plot suggests that allocating 500-1000 replications to variance estimation achieves a good balance of accuracy and width, in the sense that the intervals from the split "500+1000" or "1000 + 500" are only slightly wider than those by other splits and that allocating less (say 250) to variance estimation results in a considerable drop in coverage probability from the nominal level 95%. The results from Tables 2-4, where the split "1000 + 500" is used, also validates the effectiveness of such a strategy. Therefore, for a given simulation budget, we recommend that the user allocate 500-1000 replications to input variance estimation with our subsampling approach and all the remaining budget to the construction of the point estimator. input data size min i n i are shown in Figure 7, where each line corresponds to either PSVB or the variance bootstrap applied to one of the four differently configured computer networks. The phenomena that we have observed in Tables 2-4 still persist for all the four computer networks. As the input data size grows, the estimation accuracy of the input variance improves to a level of 0.3 in relative rmse for PSVB, but deteriorates significantly for the variance bootstrap. Accordingly the CI coverage of PSVB stays within a 1% margin around the nominal level 95%, whereas the variance bootstrap shows a significantly lower coverage than the nominal level because of inaccurate input variance estimates. All these demonstrate that the proposed guidelines for using PSVB deliver superior and robust performance across different systems, and therefore can be used as a default algorithmic configuration in practice. However, the relatively low accuracy of the variance estimates and CIs appears again in our approach when the input data size is limited (around 30). The same limitation arises for the variance bootstrap. This reconciles with our observation from Section 5.2 that subsampling is most beneficial for cases with moderately large input data sizes where other approaches like the variance bootstrap start to become computationally demanding.

Conclusion
We have explained how estimating input variances in stochastic simulation can require large computation effort when using conventional bootstrapping. This arises as the bootstrap involves a two-layer sampling, which adds up to a total effort of larger order than the data size in order to achieve relative consistency. To alleviate this issue, we have proposed a subsampling method that leverages the relation between the structure of input variance and the estimation error from the two-layer sampling, so that the resulting total effort can be reduced to being independent of the data size. We have presented the theoretical results in this effort reduction, and the optimal choices of the subsample ratio and simulation budget allocation in terms of the data size and the budget. We have also demonstrated numerical results to support our theoretical findings, and provided guidelines in using our proposed methods to estimate input variances and also construct output CIs. Future work comprises a more comprehensive investigation of our subsampling scheme, including its generalization to input processes with serial dependence and potentially non-smooth performance measures such as quantiles and other risk measures.

Proofs of Statements
We first verify the proposed assumptions for the special case of finite-horizon performance measures in Section EC.1. Section EC.2 then proves results on the validity of the input variance decomposition (2). Section EC.3 proves the consistency of the proposed input variance estimate and analyzes its Monte Carlo error in relation to the parameters B, R. Lastly, Section EC.4 further analyzes the statistical error to obtain the overall error of the input variance estimate, and derives the optimal choices for θ, B, R that minimizes the overall error. In all the proofs, we write a ≈ b to mean a/b p → 1.

EC.1. Finite-Horizon Performance Measures
In this section, we show that Assumptions 2-7 and 10-12 hold for the finite-horizon performance measure (16), thereby proving Theorems 6 and 9. We first prove Assumptions 2 and 10, then present the useful Lemma EC.1 which will later be used to prove all other assumptions.
Proof of Assumptions 2 and 10. The finite horizon structure allows the following expansion of the performance measure ψ(P ν 1 1 , . . . , P νm m ) around the input models P 1 , . . . , P m ψ(P ν 1 1 , . . . , P νm m ) Here each x i,T i := (x i,t ) t∈T i . Expressing terms with d = 1, 2, 3 in a more explicit form gives Since each signed measure Q i − P i in the product measure in (EC.2) has zero total measure, adding to the integrand a function that is independent of at least one of the integration variables does not change the integral value. Hence one can replaceg's by the following centered versions for where X i , X i , X i denote independent variables distributed under F i , and replace the function where the order (i 1 , t 1 ) < (i 2 , t 2 ) is defined as either i 1 < i 2 , or i 1 = i 2 but t 1 < t 2 . This leads to the new Taylor expansion Note that now all the integrands above have zero marginal means due to centering, e.g.
h c T 1 ,...,Tm (x 1,T 1 , . . . , x m,Tm )dP i (x i,t ) = 0 for all i and t ∈ T i . (EC.5) However, the functionsg c i ,g c i 1 i 2 ,g c i 1 i 2 i 3 are not necessarily symmetric under permutations as required in Assumption 10, so we perform the following symmetrization to find the influence functions for all π let g i π(1) i π(2) i π(3) (x π(1) , x π(2) , x π(3) ) :=g c i 1 i 2 i 3 (x 1 , x 2 , x 3 ) for i 1 < i 2 < i 3 ec4 e-companion to Lam and Qian: Subsampling for Input Uncertainty Quantification where the dependence on P 1 , . . . , P m is suppressed and π denotes any permutation of (1, 2, 3). Then one can check that g i 1 i 2 and g i 1 i 2 i 3 not only retain the property of zero marginal means, but also satisfy the symmetry condition in Assumption 10. Permutation symmetry implies that We continue to verify other assumptions, for which we use the following lemma.
Lemma EC.1. Suppose Assumption 9 holds with positive and even k. For each i let F i ∈ {F i , F i } be either the i-th true or empirical input model. Then the following bounds hold uniformly for every where the influence functions g i 's are now under the true input models F 1 , . . . , F m . Each empirical influence functionĝ i satisfies Therefore E F 1 ,..., Fm [h k (X 1,I 1 , . . . , X m,Im )] = O p (1) for each I 1 , . . . , I m . Since there are finitely many of them, the maximum is also bounded in probability. This proves the first bound.
To explain the other bounds, we put ψ( F 1 , . . . , F m ) in the form of the expansion (EC.4) with where h c T 1 ,...,Tm (x 1,T 1 , . . . , x m,Tm )dF i (x i,t ) = 0 for all i and t ∈ T i , according to the property of zero marginal means (EC.5). To obtain a moment bound for h c T 1 ,...,Tm , observe that by Assumption 9 and Jensen's inequality any conditional expectation of the performance function h has a k-th By (EC.5) the expectation on the right hand side is zero if some data point X i,j appears only once.
Therefore the number of non-zero expectations is bounded above by C(k, m, Moreover, from (EC.12) each expectation satisfies the following by generalized Hölder's inequality This gives the second bound.
The third bound can be established by the same argument, but considering only the remainders for which d ≥ 2. We then prove the bounds on influence functions. According to the representation of g i (P 1 , . . . , P m ; ·) in the proof of Assumptions 2 and 10, the empirical influence functionĝ i iŝ First we derive the following Taylor expansion for each conditional expectation (i 1 ,t 1 )<(i 2 ,t 2 ),t 1 ∈T i 1 ,t 2 ∈T i 2 h (i,t),T 1 ,...,Tm dF i 1 (x i 1 ,t 1 )dF i 2 (x i 2 ,t 2 ) then we have the following parallel property of (EC.5) h c (i,t),T 1 ,...,Tm (x 1,T 1 , . . . , x i,T i ∪{t} , . . . , x m,Tm )dF i (x i ,t ) = 0 for all i and t ∈ T i and by comparing the first order remainders (EC.13) and (EC.14) ofĝ i with the second order influence functions g i 1 i 2 it is easy to establish that Proof of Assumption 4. The convergence ofĝ i to g i in fourth order moment is a direct con-  Proof of Assumption 6. We writeτ 2 = τ 2 (F 1 , . . . , F m ) for short. First rewrite Applying Lemma EC.1 to h 2 (k = 2) with the true distributions being F 1 , . . . , F m we get

EC.2. Proofs of Propositions 1 and 2
This section proves results concerning the validity of the additive decomposition 2 of the input variance. We first prove Proposition 1, and then provide the key Lemma EC.2 that will be used in the proof of Proposition 2 as well as many results in Section EC.4.
Proof of Proposition 1. Following the expansion (14)  This completes the proof.
The following important lemma on variance decomposition plays a crucial role in our analysis. Lemma EC.2 (ANOVA decomposition, adapted from Efron and Stein (1981)). Let Y i , i = 1, . . . , n be independent but not necessarily identically distributed random variables, and φ(y 1 , . . . , y n ) be a function such that E[φ 2 (Y 1 , . . . , Y n )] < ∞, then there exist functions φ i 1 ,...,i k for 1 ≤ i 1 < · · · < i k ≤ n and k ≤ n such that Moreover, the 2 n − 1 random variables in the decomposition have mean zero and are mutually uncorrelated.
With this lemma, we can prove Proposition 2: Proof of Proposition 2. The proof of Proposition 1 derives the following expression for input    Here a ≈ b means a/b p → 1 as aforementioned. Therefore the variance can be expressed as where the second equality holds because of the convergence ofσ 2 i ,τ 2 to σ 2 i , τ 2 . Proof of Theorem 8. One can easily verify that such B * and R * minimize the mean squared error (18) under the constraint that BR = N and B = ω(1). The mean square error (19) then follows from evaluating (18) at B * , R * .
Proof of Corollary 1. It is obvious that when B = ω(1) and R = Ω(θn) the configuration (7) is satisfied hence the estimateσ 2 SV B is relatively consistent under such allocation. To show that a simulation budget N = ω(θn) is necessary for (7) to hold, note that multiplying the first two requirements in (7) gives that B 2 R 2 = ω((θn) 2 ), hence BR = ω(θn) must hold true.

EC.4. Proofs for Results in Section 4.3 and Theorem 5
In this section we analyze the statistical error of σ 2 SV B relative to the true input variance σ 2 I , therefore, combined with the Monte Carlo errorσ 2 SV B − σ 2 SV B given in Lemma 1, provide the overall error of the estimateσ 2 SV B , and then minimize the overall error to obtain the optimal choices for the parameters θ, B, R. We first prove Lemma 2 using Lemma EC.1 and Proposition 2 which have been presented in Section EC.2, then use Lemma 2 to conclude Theorem 10. Lastly, Theorem 5 is derived from Theorem 10. As is the case in the proof of Proposition 2, the covariances can be simplified to This leads to From the above expression of the variance estimator one can verify that it suffices to show the following three results To see this, if the three equations hold then      True arrival rates λi,j for system #4.