Proceedings of the 2018 Winter Simulation Conference
M. Rabe, A. A. Juan, N. Mustafee, A. Skoogh, S. Jain, and B. Johansson, eds.

SUBSAMPLING VARIANCE FOR INPUT UNCERTAINTY QUANTIFICATIONHenry Lam
Huajie Qian
Department of Industrial Engineering and Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027, USA

ABSTRACTIn stochastic simulation, input uncertainty refers to the output variability arising from the statistical noisein specifying the input models. This uncertainty can be measured by a variance contribution in the output,
which is estimated commonly via the bootstrap. However, due to the convolution of the simulation noise andthe input noise, the computation effort required in the existing bootstrap schemes are typically substantial.
This paper investigates a subsampling framework as a computation saver. We demonstrate the strengthsof our subsampled bootstrap in terms of theoretical computation requirements, and substantiate them withnumerical illustrations.
1

INTRODUCTION

Stochastic simulation is one of the most widely used analytic tools in operations research. It provides aflexible mean to approximate complex models and subsequently 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 andsuboptimal decision-making. In the literature, this problem is commonly known as input uncertainty or
extrinsic uncertainty.
This paper concerns the quantification of input uncertainty. In conventional simulation output analysis
where the input model is completely pre-specified, the statistical errors come solely from the Monte Carlonoises, and it suffices to account 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 theadditional uncertainty is to estimate the variance in the output that is contributed from the input noises(Song et al. 2014); for convenience, we call this the input variance. Moreover, the input variance can beused to identify models that are overly ambiguous and inform more data collection. They also collectively
provide a building block to construct valid output confidence intervals (CIs) that account for combined
input and simulation errors (Cheng and Holland 2004).
Bootstrap resampling is a common approach to estimate input variances. This applies most prominently
in the nonparametric regime, namely when no assumptions are placed on the input parametric family. But
it could also be used in the parametric case (where more alternatives are available). For example, Cheng
and Holland (1997) proposes the variance bootstrap, and Song and Nelson (2015) studies the consistencyof this strategy on a random-effect model that describes the uncertainty propagation. A bottleneck withusing bootstrap resampling in estimating input variances, however, is the need to “outwash” the simulationnoise, which often places substantial burden on the required simulation effort. More precisely, with boththe input and the simulation noises, the bootstrap procedure typically comprises a two-layer sampling that

Lam and Qian
first resamples the input data (i.e., outer sampling), followed by running simulation replications using eachresample (i.e., inner replications). On the other hand, the magnitude of the input variance typically scaleswith the amount of input data. It turns out that, relative to the data size, the procedure will require either
a large outer bootstrap size or a large inner replication size in order to extinguish the effect of simulationnoises. This requirement is often ignored in the existing literature by assuming the power to generate
plenty of bootstrap samples (there are exceptions, e.g., the recent work of Song and Nelson (2018)). Whensimulating high-fidelity models, where each simulation takes hours or days, one may not afford to havesuch computation resources.
The main goal of this paper is to investigate subsampling as a simulation saver for input varianceestimation. In particular, we show how a judicious use of subsampling can lead to estimation errors that
depend much more favorably on the outer and inner sample sizes. This approach leverages the relation
between the structure of the input variance in terms of the data size and the estimation error incurred in therequired two-layer sampling. In the statistics literature, subsampling has been proposed to tackle situations
where the conventional bootstrap fails, due to a lack of uniform convergence required for consistency, or intime series where data are non-i.i.d. (Politis et al. 1999). It typically involves sampling without replacement
so that each resample consists of truthful observations. In contrary, our subsampling approach is introduced
to reduce the simulation effort faced by the two-layer sampling. It serves to more efficiently deconvolutethe effect of the simulation noise from the input noise. We will illustrate how to incorporate our scheme
in the variance bootstrap that was studied in Cheng and Holland (1997) and Song and Nelson (2015), andwill analyze the involved simulation complexity (i.e., minimally required total simulation replication size)
and allocation rules.
We close this introduction with a brief review of other related work in input uncertainty. In the
nonparametric case, Barton and Schruben (1993) and Barton and Schruben (2001) use the basic bootstrap
to construct CIs, where the CI limits are determined from the quantiles of the bootstrap distributions.
Lam and Qian (2016), Lam and Qian (2017) study the use of empirical likelihood, and Xie et al. (2018)studies nonparametric Bayesian methods. In the parametric case, Barton et al. (2013) studies the percentile
bootstrap with a metamodel built in advance, a technique known as the metamodel-assisted bootstrap.
Cheng and Holland (1997) studies the delta method, and Cheng and Holland (1998), Cheng and Holland(2004) reduce its computation burden via the so-called two-point method. Lin et al. (2015) and Song and
Nelson (2018) study regression approaches to estimate sensitivity coefficients which are used to apply thedelta method. Finally, Chick (2001), Zouaoui and Wilson (2004) and Xie et al. (2014) study varianceestimation and interval construction from a Bayesian perspective. For general surveys on input uncertainty,
readers are referred to Barton et al. (2002), Henderson (2003), Barton (2012), Song et al. (2014) and Lam(2016).
Due to space limit, this paper only gives a concise discussion on some of our results, and leaves further
details and extensions, including a related approach that uses subsampling on a bootstrap aggregating
scheme, to a full journal paper.
The remainder of the paper is as follows. Section 2 introduces the input uncertainty problem andexplains the computation bottleneck in the existing bootstrap scheme. Section 3 presents our subsamplingidea, procedures and the main statistical results. Section 4 reports some numerical experiments.
2

PROBLEM MOTIVATION

We describe the problem and our motivation. Section 2.1 first describes the input uncertainty problem,
Section 2.2 discusses the existing bootstrap approach, and Section 2.3 discusses its computational barrier,thus motivating our subsampling investigation.

Lam and Qian
2.1 The Input Uncertainty Problem
Suppose there are m independent input processes driven by input distributions F1 , F2 , . . . , Fm . We consider ageneric performance measure ψ(F1 , . . . , Fm ) 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 F1 and F2 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.
The input uncertainty problem arises in situations where the input distributions F1 , . . . , Fm are unknown
but real-world data are available. One then has to use their estimates Fb1 , . . . , Fbm to drive the simulation.
Denote a point estimate of ψ(F1 , . . . , Fm ) as ψ̄(Fb1 , . . . , Fbm ), where typically we takeR1
ψ̄(Fb1 , . . . , Fbm ) = ∑ ψ̂r (Fb1 , . . . , Fbm )R r=1
with ψ̂r (Fb1 , . . . , Fbm ) being a conditionally unbiased simulation replication driven by Fb1 , . . . , Fbm . This point
estimate is affected by both the input statistical noises and the simulation noises. By conditioning on
the estimated inputs (or viewing the point estimate as a random effect model with uncorrelated input andsimulation noises), the variance of ψ̄(Fb1 , . . . , Fbm ) can be expressed asVar[ψ̄(Fb1 , . . . , Fbm )] = σI2 + σS2
whereσI2 = Var[ψ(Fb1 , . . . , Fbm )](1)

is interpreted as the overall input variance, andσS2 =E[Var[ψ̂r (Fb1 , . . . , Fbm )|Fb1 , . . . , Fbm ]]R

as the variance contributed from the simulation noises. Assuming that the estimates Fbi ’s are consistentin estimating Fi ’s, then, as ni grows, σS2 is approximately Var[ψ̂r (F1 , . . . , Fm )]/R and can be estimated bytaking the sample variance of all simulation replications (see, e.g., Cheng and Holland 1997). The key andthe challenge in quantifying input uncertainty is to estimate σI2 .
To this end, suppose further that for each input model i, we have ni i.i.d. data {Xi,1 , . . . , Xi,ni } generated
from the distribution Fi . Typically, the overall input variance is decomposable, when ni ’s are large, intom

σi2
i=1 niσI2 ≈ ∑

(2)

where σi2 /ni is the variance contributed from the data noise for model i. In the parametric case where Fbi
comes from a parametric family containing the estimated parameters, this decomposition is well known
from the delta method (Asmussen and Glynn 2007, Chapter 3). Here, σi2 /ni is typically ∇i ψ 0 Σi ∇i ψ, 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 reciprocallyiwith ni ). In the nonparametric case where the empirical distribution Fbi (x) := ∑nj=1
1{Xi, j ≤ x}/ni is used,(2) still holds under mild conditions (e.g., Assumption 1 in the sequel). The quantity σi2 will be equal to
VarFi [gi (Xi )], where gi (·) is the influence function (Hampel 1974) of ψ with respect to the distribution Fi ,defined on the value space of the input variate Xi , and VarFi [·] denotes the variance under Fi . The influencefunction can be viewed as a Gateaux derivative taken with respect to the probability distributions Fi ’s (seeSerfling 2009, Chapter 6), and dictates the first-order asymptotic variance of the plug-in estimate of ψ.
Related derivations can be found in, e.g., Lam and Qian (2017).

Lam and Qian
Under further regularity conditions, a Gaussian approximation holds for ψ̄(Fb1 , . . . , Fbm ) so thatqψ̄(Fb1 , . . . , Fbm ) ± z1−α/2 σI2 + σS2(3)

is an asymptotically tight (1 − α)-level CI for ψ(F1 , . . . , Fm ), where z1−α/2 is the standard normal 1 − α/2quantile. This CI, which provides a bound-based alternative to quantify input uncertainty, again requiresa statistically valid estimate of σI2 (and σS2 ).
2
We will investigate how to estimate ∑mi=1 σi /ni . We focus primarily on the nonparametric case. The
most common estimation technique is bootstrap resampling, which we discuss next.
2.2 Bootstrap Resampling
Let Fbi∗ represent the empirical distribution constructed using a bootstrap resample from the original
data {Xi,1 , . . . , Xi,ni } for input Fi , i.e., ni points drawn by uniformly sampling with replacement from{Xi,1 , . . . , Xi,ni }. The bootstrap variance estimator is Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )], where Var∗ [·] denotes the varianceover the bootstrap resamples from the data, under Fb1 , . . . , Fbm .
The principle of bootstrap entails that Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )] ≈ Var[ψ(Fb1 , . . . , Fbm )] and consistently approximates the asymptotic variance σI2 under mild conditions. Here Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )] is obtained froma (hypothetical) infinite number of bootstrap resamples and simulation runs per resample. In practice,
however, one would need to use a finite bootstrap size and a finite simulation size. This comprises B
conditionally independent bootstrap resamples of {Fb1∗ , . . . , Fbm∗ }, and R simulation replications driven byeach realization of resampled input distributions. This generally incurs two layers of Monte Carlo errors.
Denote ψ̂r (Fb1b , . . . , Fbmb ) as the r-th simulation run driven by the b-th bootstrap resample {Fb1b , . . . , Fbmb }.
Denote ψ̄ b as the average of the R simulation runs driven by the b-th resample, and ψ̄¯ as the grand sample
average from all the BR runs. An unbiased estimator for Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )] is given byV
1 B∑ (ψ̄ b − ψ̄¯ )2 − R
B − 1 b=1
where
V=

(4)

B R
1
∑ ∑ (ψ̂r (Fb1b , . . . , Fbmb ) − ψ̄ b )2 .
B(R − 1) b=1
r=1

To explain, the first term in (4) is an unbiased estimate of the variance of ψ̄ b , which is Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )] +(1/R)E∗ [Var[ψ̂r (Fb1∗ , . . . , Fbm∗ )|Fb1∗ , . . . , Fbm∗ ]] (where E∗ [·] denotes the expectation on Fbi∗ ’s under Fbi ’s) since ψ̄ b
incurs both the bootstrap noise and the simulation noise. In other words, the variance of ψ̄ b is upward biased
for Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )]. The second term in (4), namely V /R, removes this bias. This bias adjustment canbe viewed as arising from Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )] being the variance of a conditional expectation. Alternately,ψ(Fb1∗ , . . . , Fbm∗ ) can be viewed as a random effect model where each “group” corresponds to each realizationof Fb1∗ , . . . , Fbm∗ , and (4) estimates the “between-group” variance in an analysis-of-variance (ANOVA). Formula(4) appeared in the input uncertainty literature, e.g., Cheng and Holland (1997), Song and Nelson (2015),
Lin et al. (2015), and also in Zouaoui and Wilson (2004) in the Bayesian context. Algorithm 1 summarizes
the procedure.
2.3 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 thesimulation complexity. This turns out to be a lot bigger than the data size.

Lam and Qian
Algorithm 1 ANOVA-based Variance Bootstrap
Parameters: B ≥ 2, R ≥ 2
for b = 1 to B dob , . . . , X b } for all i to obtain bootstrapped empirical distributions F
bib
Draw a resample {Xi,1i,nifor r = 1 to R do
Draw a ψ̂r (Fb1b , . . . , Fbmb )end forCompute ψ̄ b = R1 ∑Rr=1 ψ̂r (Fb1b , . . . , Fbmb )end for1
Compute V = B(R−1)∑Bb=1 ∑Rr=1 (ψ̂r (Fb1b , . . . , Fbmb ) − ψ̄ b )2 and ψ̄¯ = B1 ∑Bb=1 ψ̄ b
Output σ̂V2B = 1 ∑Bb=1 (ψ̄ b − ψ̄¯ )2 − V
B−1R

To facilitate discussion, denote n as a scaling of the data size (i.e., assume ni all grow linearly with n).
Note that σ̂V2B is unbiased for Var∗ [ψ(Fb1∗ , . . . , Fbm∗ )], which is approximately σI2 . Thus, to analyze the errorof σ̂V2B , it suffices to focus on its variance. An ANOVA-type analysis as in Sun et al. (2011) gives
1
1
Var∗ [σ̂V2B ] = Θ p+(5)
Bn2 BR2
where Θ p denotes stochastic boundedness (similarly, o p will denote stochastic asymptotic negligibility).
Stochastic boundedness is used here, as well as in later sections, since Var∗ [σ̂V2B ] depends on the input
data and hence is random. The two terms in (5) correspond to the variance coming from the bootstrap
resampling and from the simulation noise respectively. From the Θ(1/n) magnitude of σI2 and in turnVar∗ [ψ(Fb1∗ , . . . , Fbm∗ )], consistency requires σ̂V2B to have a variance o(1/n2 ). This implies unfortunately thatthe total number of simulation runs, BR, must be ω(n), i.e., of order higher than the data size, because
only then the second term in the bound (5) can be o(1/n2 ). We summarize with the following result.
Theorem 1 Under the assumptions in Section 3.2.1, the required simulation budget N = BR to achievep(relative) consistency in estimating σI2 by Algorithm 1, i.e., σ̂V2B /σI2 → 1, is N = ω(n).
Finally, we note that such a computational barrier does not only occur in variance bootstrap. For
instance, the basic bootstrap studied in Barton and Schruben (1993), Barton and Schruben (2001) appears
to also require a dominating inner replication size than the data size in order to obtain valid quantile
estimates (the authors actually used one inner replication, but Barton (2012) commented that more isneeded). Similarly, the regression approach in Song and Nelson (2018) requires simulation effort of ordernγ where 1 < γ < 2. Likewise, the empirical likelihood framework studied in Lam and Qian (2017) also
requires a 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.
3

PROCEDURES AND GUARANTEES IN THE SUBSAMPLING FRAMEWORKThis section presents our methodologies and results on subsampling. Section 3.1 first discusses how toincorporate subsampling in the variance bootstrap estimator depicted in Algorithm 1. Section 3.2 thenprovides theoretical guarantees of the subsampled procedure.
3.1 Proportionate Subsampled Variance Bootstrap
As explained before, the reason why the σ̂V2B in Algorithm 1 requires a huge simulation effort, as impliedby its mean square error (5), 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 standard
deviation relatively vanishes. This requirement manifests in the inner replication size in constructing σ̂V2B .

Lam and Qian
To reduce the inner replication size, we leverage the relation between the form of the input varianceand the estimation variance depicted in (5) as follows. The approximate input variance contributed frommodel i, with data size ni , has the form σi2 /ni . If we use the variance bootstrap directly as in Algorithm
1, then we need an order more than ni total simulation runs due to (5). Now, pretend that we have fewer
data, say si , then the input variance will be σi2 /si , and the required simulation runs is now only of orderhigher than si . An estimate of σi2 /si , however, already gives us enough information in estimating σi2 /ni ,because we can merely rescale our estimate of σi2 /si by si /ni to get an estimate of σi2 /ni . This estimation
can be done by subsampling the input distribution with size si . Thus, by subsampling to estimate σi2 /si ,we can both use fewer simulation runs and also retain correct estimation via multiplying by a si /ni factor.
To make the above argument more transparent, the bootstrap principle and the asymptotic approximationof the input variance imply thatm2σVar∗ [ψ(Fb1∗ , . . . , Fbm∗ )] = ∑ i (1 + o p (1)).
i=1 niThe subsampling approach builds on the observation that a similar relation holds form

σi2(1 + o p (1))
i=1 si∗∗, . . . , Fbm,sVar∗ [ψ(Fb1,s)] = ∑
m1∗ denotes a bootstrapped distribution of size s (i.e., the empirical distribution that is uniformly
where Fbi,siisampled with replacement from {Xi,1 , . . . , Xi,ni }). If we let si = bθ ni c for some θ > 0 so that si → ∞, thenwe havemσi2∗b∗(1 + o p (1)).
Var∗ [ψ(Fb1,bθ
,
.
.
.
,
F
)]
=
∑
n1 c
m,bθ nm c
i=1 θ ni

Multiplying both sides with θ , we getm

σi2(1 + o p (1)).
i=1 ni∗b∗θ Var∗ [ψ(Fb1,bθ
n1 c , . . . , Fm,bθ nm c )] = ∑Note that the right hand side above is the original input variance of interest. This leads to our proportionate
subsampled variance bootstrap summarized in Algorithm 2.
Algorithm 2 Proportionate Subsampled Variance Bootstrap
Parameters: B ≥ 2, R ≥ 2, 0 < θ ≤ 1Compute si = bθ ni c for all ifor b = 1 to B dob , . . . , X b } with replacement from data for all i, which induces theDraw a uniform subsample {Xi,1i,sibempirical distribution Fbi,si
for r = 1 to R do
b ,...,Fb )bm,sDraw a ψ̂r (Fb1,sm
1
end forb ,...,Fb )bm,sCompute ψ̄ b = R1 ∑Rr=1 ψ̂r (Fb1,sm
1
end for1b ,...,Fb ) − ψ̄ b )2 and ψ̄bm,s¯ = 1 ∑Bb=1 ψ̄ bCompute V = B(R−1)∑Bb=1 ∑Rr=1 (ψ̂r (Fb1,s
B
m
1
B
1
2
b − ψ̄
2−V)¯
Output σ̂SV
=
θ
(
(
ψ̄
)
∑
b=1
B
B−1

R

Lam and Qian
3.2 Statistical Guarantees
We present our key theoretical guarantees for Algorithm 2. Section 3.2.1 first states our assumptions on
the performance measure. Sections 3.2.2 then presents the theories on estimation consistency, simulation
complexity and budget allocation for the proportionate subsampled variance bootstrap.
3.2.1 Regularity Assumptions
We assume the data sizes of different input models are of the same order. To be specific, maxi ni / mini ni ≤ cfor some finite constant c > 0 as ni → ∞. We first make an assumption on the smoothness of the performancemeasure:Assumption 1 (Smoothness at the true input distributions) The performance measure ψ, as a functionalof the input distributions, has first order influence functions gi : R → R, i = 1, . . . , m at F1 , . . . , Fm such that
EFi [gi (Xi )] = 0, 0 < VarFi [gi (Xi )] < ∞ for all i and the plug-in estimator ψ(Fb1 , . . . , Fbm ) admits the Taylor
expansion
m
1 ni
(6)
ψ(Fb1 , . . . , Fbm ) = ψ(F1 , . . . , Fm ) + ∑ ∑ gi (Xi, j ) + ε
i=1 ni j=1
where the remainder ε satisfies E[ε 2 ] = o(n−1 ).
Assumption 1 entails that the performance measure is differentiable and the error of the linear approximation is negligible in an asymptotic sense. Specifically, the linear term in (6) is Θ p (n−1/2 ) by thecentral limit theorem, whereas the error ε is o p (n−1/2 ). Hence the variance of the linear term contributes
dominantly to the overall input variance as ni ’s are large, which gives rise to the additive form of the input
variance thanks to the independence among the input models. Thus we have:2
Proposition
1 Under Assumption 1, the input variance σI2 defined in (1) takes the form σI2 = ∑mi=1 σi /ni +−12o n , where each σi = VarFi [gi (Xi )] is the variance of the i-th influence function.
The next assumption concerns the bootstrap principle under subsampling:
Assumption 2 (Smoothness at the empirical input distributions) The performance measure ψ has firstorder influence functions ĝi : R → R, i = 1, . . . , m at Fb1 , . . . , Fbm such that
1 ni
1 ni
1 ni 4pĝi (Xi, j ) = 0, σ̂i2 := ∑ ĝ2i (Xi, j ) → σi2 ,∑∑ ĝi (Xi, j ) = O p (1) for all ini j=1
ni j=1
ni j=1(7)

and the bootstrapped performance measure with subsample size si admits the Taylor expansion
m

si

m

1 3 
.

1∗∗∗) + ε∗
ψ(Fb1,s, . . . , Fbm,s) = ψ(Fb1 , . . . , Fbm ) + ∑ ∑ ĝi (Xi,km1si
i=1 k=1
where the remainder ε ∗ satisfies
E∗ [(ε ∗ − E∗ [ε ∗ ])4 ] = o p

∑ si

i=1

The ĝi ’s in Assumption 2 are the empirical influence functions and are expected to approach the true
influence functions gi ’s in Assumption 1. This leads to the first and second conditions in (7). The fourth
moment conditions on the empirical influence functions, i.e., the last condition in (7), and on the remainderε ∗ are needed for controlling the variance of our variance estimator.
To introduce our assumptions on the simulation noise, denote by τ 2 the simulation variance of ψ̂ under
the true input models, i.e.,τ 2 = EF1 ,...,Fm [(ψ̂(F1 , . . . , Fm ) − ψ(F1 , . . . , Fm ))2 ].

Lam and Qian
Similarly, we use τ̂ 2 and τ̂∗2 to represent the variance of ψ̂ conditional on the empirical input models andthe bootstrapped input models∗∗∗∗
τ̂ 2 = EFb1 ,...,Fbm [(ψ̂(Fb1 , . . . , Fbm )−ψ(Fb1 , . . . , Fbm ))2 ], τ̂∗2 = EFb∗ ,...,Fb∗ [(ψ̂(Fb1,s, . . . , Fbm,s)−ψ(Fb1,s, . . . , Fbm,s
))2 ]
m
m
1
1
1,s1

m,sm

and use µ̂4∗ to represent the central fourth moment of ψ̂ under the bootstrapped input models∗∗∗∗µ̂4∗ = EFb∗ ,...,Fb∗ [(ψ̂(Fb1,s, . . . , Fbm,s) − ψ(Fb1,s, . . . , Fbm,s))4 ].
m
m
1
1
1,s1

m,sm

The assumptions on the simulation noise are:p

Assumption 3 (Convergence of variance) τ̂ 2 → τ 2 and E∗ [(τ̂∗2 − τ̂ 2 )2 ] = o p (1).
Assumption 4 (Boundedness of the fourth moment of ψ̂) E∗ [µ̂4∗ ] = O p (1).
Assumption 3 stipulates that the variance of the Monte Carlo estimator ψ̂ as a functional of theunderlying input models is smooth enough in the inputs. Assumption 4 is again a fourth moment condition
used to control the variance of our variance estimator.
Although the above assumptions may look complicated, one can verify that all of them hold for a wideclass of performance measures, namely those over finite time horizons in the formψ(F1 , . . . , Fm ) = EF1 ,...,Fm [h(X1 , . . . , Xm )](8)

where each Ti is a deterministic time, the i-th input process Xi = (Xi (1), . . . , Xi (Ti )) consists of Ti i.i.d. variables
distributed under Fi , and h is some performance function. An unbiased simulation estimate ψ̂r of theperformance measure simply outputs a copy of h(X1 , . . . , Xm ). Assume the performance function h satisfies
the following:
Ti
Assumption 5 For each i, 0 < VarFi [∑t=1
EF1 ,...,Fm [h(X1 , . . . , Xm )|Xi (t) = Xi ]] < ∞.
Assumption 6 For each i let Ii = (Ii (1), . . . , Ii (Ti )) be a sequence of indices such that 1 ≤ Ii (t) ≤ Ti , andXi,Ii = (Xi (Ii (1)), . . . , Xi (Ii (Ti ))). Assume maxI1 ,...,Im EF1 ,...,Fm [|h(X1,I1 , . . . , Xm,Im )|4 ] < ∞.
Then we have:
Proposition 2 Under Assumptions 5 and 6, we have Assumptions 1-4 hold for the finite-horizon performancemeasure ψ given by (8).
3.2.2 Simulation Complexity and Optimal Allocation
Recall from Section 3.1 that the proportionate subsampled variance takes the form
2
b∗
b∗
σSV
B = θ Var∗ [ψ(Fbθ n1 c,1 , . . . , Fbθ nm c,m )](9)where 0 < θ ≤ 1 is the subsample ratio, and b·c outputs the integer part of a number. In Section 3.1 wehave presented a heuristic argument on consistency of (9). Here we give a formal statement:
Theorem 2 Under Assumptions 1 and 2, if the subsample ratio θ = ω(n−1 ), then the subsampled variancep2 /σ 2 →(9) is consistent as ni → ∞, i.e. σSV
1.
I
B
−1
The requirement θ = ω(n ) implies that si → ∞, 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 isenough to ensure consistency even relative to the magnitude of σI2 .
Now we turn to the discussion on the Monte Carlo error, from both the outer bootstrapping and the2
inner simulation, for σ̂SV
B in Algorithm 2. The following lemma characterizes the level of Monte Carlonoise in terms of the conditional mean square error.

Lam and Qian
Lemma 1 Under Assumptions 1-4, if B = ω(1), θ = ω n−1 , then the mean square error conditional on
the input data is
2  m σi2 τ 2 θ 222
2E∗ [(σ̂SV
−
σ
)
]
=
(10)
∑ ni + R (1 + o p (1)).
B
SV B
B i=1
The target quantity σI2 here is of order n−1 by Proposition 1, and hence the Monte Carlo noise of ourvariance estimate has to vanish faster than n−1 in order to achieve consistency. This leads to the following
simulation complexity of Algorithm 2:p2
22
Theorem 3 Under Assumptions 1-4, the variance estimate σ̂SV
B is consistent, i.e. σ̂SV B /σI → 1 if andonly if the parameters B, R, θ of Algorithm 2 are chosen such that1.
B = ω(1), BR2 = ω (θ n)2 , θ = ω
n

(11)


In particular, the minimum requirement of total simulation budget is N = BR = ω(1) by choosing θ = ω n−1 .
Theorem 3 tells us what values of B, R would work in general. For a given subsample ratio θ , thefirst two conditions in (11) implies a simulation budget of ω(θ n) for consistent estimation, which recovers
the result in Theorem 1 by setting θ = Θ(1). Next, we establish the optimal configuration of B, R, θ tominimize the Monte Carlo error:
Theorem 4 Under Assumptions 1-4, given a simulation budget N and a subsample ratio θ such thatN = ω(θ n) and θ = ω(n−1 ), the optimal outer and inner sizes that minimize the conditional mean squareerror (10) areθ τ2NB∗ = ∗ (1 + o p (1)), R∗ = 2 (1 + o p (1))
RσI
2 − σ 2 )2 ] = (1 + o (1))8θ τ 2 σ 2 /N.
giving a conditional mean square error E∗ [(σ̂SV
p
I
B
SV B
Since σI2 = Θ(n−1 ), the theorem indicates that the optimal choice of inner size R is of order Θ(θ n),
the same as the subsample size, and the outer size B is then chosen accordingly.
22Lastly, we show a result on the overall errors in using σ̂SV
B to estimate σI . This comprises of both
2 , and the input error between σ 2
2
the Monte Carlo error in approximating σSV
B
SV B and σI . The latter is notaffected by the choices of B, R, but by the subsample ratio θ . The optimal choice of θ in minimizing thegross error balances a trade-off between the above two sources of errors. Due to space limit, we skip thefull details here. We show a result specialized to the case of finite-horizon performance measures:
Theorem 5 Suppose ψ is a finite-horizon performance measure (8), and Assumptions 5 and 6 hold. The
2 − σ 2 are respectively
optimal θ , R, B given simulation budget N = ω(1) in minimizing the gross error σ̂SVIB(θ ∗ = Θ N 1/3 n−1
if N ≤ n3/2
N
, R∗ = Θ(θ ∗ n), B∗ = ∗
−1/2
∗
−2
3/2
R
Θ(n
) ≤ θ ≤ Θ Nn ∧ 1
if N > ngiving the minimum mean square error E[(σ̂SV B − σI2 )2 ] = Θ 1/((N 2/3 ∧ n)n2 ) .

4

NUMERICAL EXPERIMENTS

We consider an M/M/1 queuing system with 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 F1 = Exp(0.5) and the service time distribution
F2 = Exp(1), for which we have n1 and n2 i.i.d. data available respectively. If At is the inter-arrival time
between the t-th and (t + 1)-th arrivals, and St is the service time for the t-th arrival, then the system output
ψ(F1 , F2 ) = EF1 ,F2 [1{W20 > 2}]

Lam and Qian
where the waiting time W20 is calculated by the Lindley recursion W1 = 0,Wt+1 = max{Wt + St − At , 0}.
We test performance of Algorithms 1 and 2 under different input data sizes. For each input data size,
1000 95%-level CIs in the form of (3) are constructed, each from an independently generated input data
set. We allocate a total number of 1000 simulation runs to estimate σI2 using either Algorithm 1 or 2,
and another 500 simulation runs driven by the empirical input distributions to compute the point estimatorψ̄(Fb1 , . . . , Fbm ). The simulation variance σS2 is calculated as (V /2 + τ̃ 2 /2)/500, where V is the within-groupvariance from Algorithm 1 or 2 and τ̃ 2 is the sample variance of the 500 simulation replications.
Tables 1, 2 and 3 show, for different algorithms and settings, the coverage, mean and standard deviationof the CI lengths, the estimated ratio between input and simulation standard deviations, and the number oftimes that the algorithm outputs a negative variance estimate (because of the debiasing, the estimate is notguaranteed to be always positive). The ratio between the input and simulation variance under each setting
demonstrates that the input uncertainty is not negligible relative to the stochastic noise.
We make a few observations. First, the use of subsampling reduces the variability of the varianceestimate when compared with standard variance bootstrap, and in turn gives rise to more accurate CIs.
This is supported by the large counts of negative variance estimates in Table 1 versus negligible counts
in Table 2 and 3, and the larger standard deviations of the CI lengths in Table 1 (7 times larger than inTable 2 and 14 times larger than in Table 3 when n1 = 4000, n2 = 2000). There are two cases in Table 1
that standard bootstrap has 90.7% and 94.0% coverage. Although they look close to 95%, this could be
because we have reset the variance estimate to 0 whenever it is negative (which happens in almost halfof the cases) and such an artificial correction makes the variance estimate upward biased, leading to alarger but possibly incorrect coverage probability. Second, with the simulation budget fixed, the coverage
probabilities in both Tables 2 and 3 increase from around 85% to the nominal level 95% as the data size
grows, validating Theorem 3. Thirdly, as shown in Table 3, under the subsample ratio (N 1/3 ∧ n1/2 )n−1
as suggested by Theorem 5 our subsampled bootstrap exhibits reasonably good performance. However,
the hidden constant there can be difficult to estimate in general. Table 2 indicates that simply using a
subsample size around 30 seems a good enough compromise for this particular example.
Table 1: Algorithm 1 with B = 50, R = 20.
input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2000, n2 = 1000
n1 = 4000, n2 = 2000

coverageestimate
84.3%
87.5%
84.5%
90.7%
94.0%

mean CI
length
0.419
0.245
0.150
0.122
0.113

std. CI
length
0.180
0.093
0.070
0.061
0.055

σI
σS

# neg var

6.96
3.77
2.26
1.76
1.56

5
43
236
409
458

Table 2: Algorithm 2 with B = 50, R = 20, θ = 30
n2 .
input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2000, n2 = 1000
n1 = 4000, n2 = 2000

coverageestimate
85.2%
89.9%
94.7%
93.8%
94.9%

mean CI
length
0.424
0.249
0.158
0.103
0.087

std. CI
length
0.183
0.069
0.028
0.012
0.008

σI
σS

# neg var

7.02
3.78
2.18
1.19
0.84

12
1
0
0
0

Lam and Qian
Table 3: Algorithm 2 with B = 100, R = 10, θ = (10 ∧input data sizes
n1 = 60, n2 = 30
n1 = 200, n2 = 100
n1 = 600, n2 = 300
n1 = 2000, n2 = 1000
n1 = 4000, n2 = 2000
5

coverageestimate
87.3%
92.0%
94.4%
94.5%
94.0%

mean CI
length
0.420
0.240
0.150
0.098
0.083

p
(n1 + n2 )/2)/((n1 + n2 )/2).

std. CI
length
0.123
0.043
0.017
0.007
0.004

σI
σS

# neg var

6.85
3.73
2.16
1.18
0.84

0
0
0
0
0

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 depends much less on the data size. We have presented some theoretical and numerical results to
support our claims. Future work comprises a more comprehensive investigation of our subsampling scheme,
including its generalization to a bootstrap aggregating approach, and more extensive numerical studies.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020,
CMMI-1523453 and CAREER CMMI-1653339/1834710.
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Stochastic Simulation: Algorithms and Analysis, Volume 57. Springer
Science & Business Media.
Barton, R. R. 2012. “Tutorial: Input Uncertainty in Output Analysis”. In Proceedings of the 2012 Winter
Simulation Conference, edited by C. Laroque et al., 1–12. Piscataway, New Jersey: IEEE.
Barton, R. R., and L. W. Schruben. 1993. “Uniform and Bootstrap Resampling of Empirical Distributions”.
In Proceedings of the 1993 Winter Simulation Conference, edited by G. W. Evans et al., 503–508.
Piscataway, New Jersey: IEEE.
Barton, R. R., and L. W. Schruben. 2001. “Resampling Methods for Input Modeling”. In Proceedings of the2001 Winter Simulation Conference, edited by B. A. Peters et al., 372–378. Piscataway, New Jersey:
IEEE.
Barton, R. R., S. E. Chick, R. C. Cheng, S. G. Henderson, A. M. Law, B. W. Schmeiser, L. M. Leemis,
L. W. Schruben, and J. R. Wilson. 2002. “Panel Discussion on Current Issues in Input Modeling”.
In Proceedings of the 2002 Winter Simulation Conference, edited by E. Yücesan et al., 353–369.
Piscataway, New Jersey: IEEE.
Barton, R. R., B. L. Nelson, and W. Xie. 2013. “Quantifying Input Uncertainty via Simulation ConfidenceIntervals”. INFORMS Journal on Computing 26(1):74–87.
Cheng, R. C., and W. Holland. 1997. “Sensitivity of Computer Simulation Experiments to Errors in Input
Data”. Journal of Statistical Computation and Simulation 57(1-4):219–241.
Cheng, R. C., and W. Holland. 1998. “Two-Point Methods for Assessing Variability in Simulation Output”.
Journal of Statistical Computation Simulation 60(3):183–205.
Cheng, R. C., and W. Holland. 2004. “Calculation of Confidence Intervals for Simulation Output”. ACM
Transactions on Modeling and Computer Simulation 14(4):344–362.

Lam and Qian
Chick, S. E. 2001. “Input Distribution Selection for Simulation Experiments: Accounting for Input Uncertainty”. Operations Research 49(5):744–758.
Hampel, F. R. 1974. “The Influence Curve and Its Role in Robust Estimation”. Journal of the American
Statistical Association 69(346):383–393.
Henderson, S. G. 2003. “Input Modeling: Input Model Uncertainty: Why Do We Care and What Should
We Do About It?”. In Proceedings of the 2003 Winter Simulation Conference, edited by S. Chick etal., 90–100. Piscataway, New Jersey: IEEE.
Lam, H. 2016. “Advanced Tutorial: Input Uncertainty and Robust Analysis in Stochastic Simulation”. In
Proceedings of the 2016 Winter Simulation Conference, edited by T. M. K. Roeder et al., 178–192.
Piscataway, New Jersy: IEEE.
Lam, H., and H. Qian. 2016. “The Empirical Likelihood Approach to Simulation Input Uncertainty”. In
Proceedings of the 2016 Winter Simulation Conference, edited by T. M. K. Roeder et al., 791–802.
Piscataway, New Jersy: IEEE.
Lam, H., and H. Qian. 2017. “Optimization-Based Quantification of Simulation Input Uncertainty via
Empirical Likelihood”. arXiv preprint arXiv:1707.05917.
Law, A. M., W. D. Kelton, and W. D. Kelton. 1991. Simulation Modeling and Analysis, Volume 2.
McGraw-Hill New York.
Lin, Y., E. Song, and B. Nelson. 2015. “Single-Experiment Input Uncertainty”. Journal of Simulation 9(3):249–259.
Politis, D. N., J. P. Romano, and M. Wolf. 1999. Subsampling. Springer.
Serfling, R. J. 2009. Approximation Theorems of Mathematical Statistics, Volume 162. John Wiley & Sons.
Song, E., and B. L. Nelson. 2015. “Quickly Assessing Contributions to Input Uncertainty”. IIE Transactions 47(9):893–909.
Song, E., and B. L. Nelson. 2018. “Input-Output Uncertainty Comparisons for Optimization via Simulation”.
Accepted in Operations Research.
Song, E., B. L. Nelson, and C. D. Pegden. 2014. “Advanced Tutorial: Input Uncertainty Quantification”. In
Proceedings of the 2014 Winter Simulation Conference, edited by A. Tolk et al., 162–176. Piscataway,
New Jersey: IEEE.
Sun, Y., D. W. Apley, and J. Staum. 2011. “Efficient Nested Simulation for Estimating the Variance of A
Conditional Expectation”. Operations Research 59(4):998–1007.
Xie, W., B. L. Nelson, and R. R. Barton. 2014. “A Bayesian Framework for Quantifying Uncertainty in
Stochastic Simulation”. Operations Research 62(6):1439–1452.
Xie, W., C. Li, and P. Zhang. 2018. “A Bayesian Nonparametric Hierarchical Framework for Uncertainty
Quantification in Simulation”. Working paper, available at homepages.rpi.edu.
Zouaoui, F., and J. R. Wilson. 2004. “Accounting for Input-Model and Input-Parameter Uncertainties inSimulation”. IIE Transactions 36(11):1135–1151.
AUTHOR BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations
Research at Columbia University. His research focuses on Monte Carlo simulation, risk and uncertainty
quantification, and stochastic optimization. His email address is khl2114@columbia.edu.
HUAJIE QIAN is a Ph.D. student in the Department of Industrial Engineering and Operations Researchat Columbia University. His research interest lies in simulation uncertainty quantification, data-driven
simulation analysis, and stochastic optimization. His email address is h.qian@columbia.edu.