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

CONSTRUCTING SIMULATION OUTPUT INTERVALS UNDER INPUT UNCERTAINTY VIADATA SECTIONINGPeter W. Glynn

Henry Lam

Department of MS&E
Stanford University
Huang Engineering Center 359A
Stanford, CA 94305, USA

Department of IEOR
Columbia University
500 W. 120th Street
New York, NY 10027, USA

ABSTRACT
We study the problem of constructing confidence intervals (CIs) for simulation outputs in the presenceof input uncertainty, where the constructed CIs capture both the statistical noises from the simulation
replications and the input data. We present a simple technique based on sectioning input data that provides
exact asymptotic confidence guarantees. Unlike some existing approaches, our technique bypasses the needto consistently estimate variances that could be computationally demanding. It can be flexibly applied todependent data and to both parametric and nonparametric input models.
1

INTRODUCTION

Stochastic simulation consists of repeated generation of random variates from input models, fed through thesystem logic, to estimate system outputs. When the input models are not fully known but estimated fromrelevant historical data (the “input data”), valid output analysis needs to account for the input statisticalnoise, in addition to the Monte Carlo noise from the random variate generation. This issue, sometimes
known as the input uncertainty, has been actively studied in recent years (see, e.g., the surveys Barton et al.
2002; Henderson 2003; Chick 2006; Barton 2012; Song et al. 2014; Lam 2016).
This paper focuses on a standard output-analytic task of constructing confidence intervals (CIs) fortarget performance measures under input uncertainty. When only simulation noise is considered, these CIs
can be constructed by conventional techniques based on, e.g., asymptotic normality (Asmussen and Glynn
2007 Chapter 3). When input noise is substantial, this conventional construction will under-cover the truemeasure, and the CIs need correspondingly inflated to retain enough coverage. In the literature, several
approaches have been studied to approximate this inflation. First is to utilize the delta method, which
requires estimating the sensitivities of the performance measure with respect to the input parameters andcombining them with the parameter sampling variance (e.g., Cheng and Holland 2004; Lin et al. 2015).
Second is to bootstrap the input data and use them to drive simulation runs to obtain bootstrapped quantiles(Barton and Schruben 1993; Barton and Schruben 2001) or variance (Cheng and Holland 1997). The
third approach follows a Bayesian perspective that samples from the posterior distributions of the input
parameters (e.g., Glynn 1986; Chick 2001; Zouaoui and Wilson 2004; Xie et al. 2014; Biller and Corlu
2011). These approaches are closely related and some of them can be integrated with metamodeling tospeed up computation (Barton et al. 2013).
While these methods can be effective in some situations, they face computational challenges in terms of
the required simulation efforts, unless more restrictive assumptions are imposed. The root of these challenges
lies in the need to consistently estimate key quantities in the CI, such as the parameter sensitivities in the delta
method, and the variance or quantile in the bootstrap. This estimation typically requires heavy simulation
effort to “outwash” the Monte Carlo noise from the input noise. Moreover, along this effort comes theGlynn and Lamneed to devise simulation allocation rules or procedural tuning, e.g., the distribution of nested simulation
efforts to the resampling of the input data and the simulation replications per resample in the bootstrap,
and the selection of the finite-difference perturbation parameter in the delta method. These procedural
configurations can determine the performance of the CIs in a delicate fashion and are not always easy to
optimize. Approaches that avoid these challenges, on the other hand, are based on additional, parametricassumptions. These include, e.g., Gaussian assumptions on the response surface in metamodel-assisted
bootstrap and on the simulation outputs in Bayesian approaches.
Motivated by these challenges, this paper proposes a simple alternative to construct CIs under minimal
simulation requirements and structural assumptions. Our method is based on sectioning the input data,
where each input section drives a batch of simulation runs which are then integrated via a t-statistic toconstruct the CI. Our distinction is that, since using the involved t-statistic does not hinge on a consistentestimation of variance, we no longer need an overwhelming simulation effort to deconvolute the Monte
Carlo noise from the input noise. As long as the output estimate in each section is approximately normal(with a variance including both sources of noise which could be unknown), we can construct a pivotal
t-statistic that leads to a valid CI. This applies to both parametric and nonparametric input regimes, and todependent data under suitable mixing conditions.
Sectioning and related techniques (on output samples) have been widely used in simulation output
analysis, such as in steady-state estimation (e.g., Glynn and Iglehart 1990) and estimating nonlinear
performance measures, e.g., quantiles, conditional value-at-risk (Asmussen and Glynn 2007; Muñoz and
Glynn 1997; Nakayama 2014). In the input uncertainty literature, these are related to the so-called direct
resampling in Barton and Schruben (2001) that divides input data and uses the quantiles of the resulting
simulation runs to construct CIs. However, this latter procedure is heuristic and serves to motivate thebootstrap schemes suggested by the authors.
In the remainder of this paper, we will give an overview of our approach (Section 2), present our main
assumptions to justify a joint central limit theorem that underpins our t-asymptotic (Section 3), and discuss
our main performance guarantees (Section 4). Due to space limit, we only give brief discussion on someof our results, and leave further details and other generalizations, including a variant of our approach using“sectioned jackknife”, to a full journal paper.
2

OVERVIEW OF OUR APPROACH

Consider estimating a simulatable performance measure driven by s input models. Notationally, we denote
this performance measure as ψ = ψ(νν ) where ν = (ν 1 , . . . , ν s ), with ν j each representing the “parameter”
of the j-th input model. Given any ν̃ν , ψ(ν̃ν ) can be approximated through simulation runs. As a genericexample, ν 1 and ν 2 can denote the rates of exponential interarrival and service time distributions respectively,and ψ(ν 1 , ν 2 ) the expected average queue length over some time horizon of a queue that can be simulated.
In this example, ν denote real-valued parameters. However, our framework is equally applicable in the
nonparametric regime, in which case ν can denote the collection of the whole probability distributions forthe input models, e.g., the interarrival and service time distributions, without restricting to any parametricfamily.
We consider the situations where ν needs to be estimated from input data, say the estimate is νb .
Assuming ψ is sufficiently smooth, the delta method implies that


σ2
b
ν
ν
(1)
ψ( ) − ψ(ν ) ≈ N 0,nwhere n is some scaling of the overall sample size (the reason for using a “scaling” is that there can bedifferent sample sizes for different input models). The quantity σ 2 /n is the variance contributed from all
the input noises.
On top of these input noises, the output bears the Monte Carlo noise incurred in simulating ψ. We denoteψ̂(νb ) as a simulation-based estimate of ψ(νb ) using input models calibrated with νb . We focus on the typicalGlynn and Lam
setting where the simulation replications are conditionally unbiased given ν , with ψ̂(νb ) = (1/r) ∑rl=1 ψ̂l (νb )and ψ̂l (νb ), l = 1, . . . , r being r i.i.d. copies of simulation replications. Under a central limit theorem (CLT),
we haveτ2ψ̂(νb ) − ψ(νb ) ≈ N 0,(2)rwhere τ 2 /r is the variance from the simulation noise, using the implicit property that νb well approximates
ν in some probabilistic sense.
We then combine (1) and (2) to obtainσ 2 τ2+(3)ψ̂(νb ) − ψ(νν ) ≈ N 0,nrwhere the variance in the normal approximation is the sum of variances contributed from the input noise
and the simulation noise.
Existing approaches in constructing CIs rely on estimating σ 2 /n, either by bootstrapping or by direct
estimation (They also estimate τ 2 /r, which is more straightforward). When using the bootstrap, consistentestimation of σ 2 /n typically requires simulation load that is of higher order than n (this can be seen by
tracing the mean square error of the bootstrapped variance using, e.g., the formula in Sun et al. 2011).
Direct estimation of σ 2 via the delta method, on the other hand, would require estimating sensitivity
coefficients whose involved effort grows with the number of parameters (e.g., Cheng and Holland 1997),
unless one uses conservative approximations (e.g., Cheng and Holland 2004). The challenge becomes more
substantial in the nonparametric case.
The key of our idea is to bypass the consistent estimation of σ 2 /n by sectioning the input data. Namely,
we divide the input data into k sections (of equal length). Using each section, we estimate νb and run
independent simulation replications to obtain an estimate ψ̂(νb ), say they are Y1 , . . . ,Yk . We then constructa pivotal statisticȲ − ψ(νν )
√
(4)
S/ k
where Ȳ and S2 are the sample mean and variance of Y1 , . . . ,Yk . If each section contains enough data anddrives enough simulation replications, the Yi ’s are approximately normal according to (3) and moreover
are i.i.d. (exactly if the data are i.i.d., and approximately if the data are serially dependent). Thus, byusing a standard relation of normal variables with t-statistics, (4) is distributed approximately as tk−1 , thet-distribution with degree of freedom k − 1. This implies thatSSȲ − tk−1,1−α/2 √ , Ȳ + tk−1,1−α/2 √kkis a (1 − α)-level CI, where tk−1,1−α/2 is the (1 − α/2)-quantile of tk−1 . Constructing this CI does not
require any knowledge about variances σ 2 /n and τ 2 /r, because the term σ 2 /n + τ 2 /r is “canceled out”
when forming the pivotal statistic (4).
Regarding the assumptions needed in the above derivation, we essentially only need the normal
approximation (3), a property satisfied for a wide range of settings when n and r are sufficiently large,
including cases where the input data are dependent under suitable mixing conditions. We note thatapproximation (3) has been used in some previous work (e.g., Cheng and Holland 1997). Since it is thebasis for our sectioning method, we first study in detail the conditions that guarantee such an approximation.
3

CENTRAL LIMIT THEOREM FOR JOINT INPUT AND SIMULATION NOISESThis section presents our assumptions and justifies in more details the central limit theorem (3). We divide
our discussion into parametric and nonparametric regimes, which are conceptually similar but involve
different sets of assumptions/notations.

Glynn and Lam
We recall that there are s input models. For each input model j, we have an available data set. Our
approach involves dividing each of these data sets into k sections. In other words, letting {Xlj }l=1,...,kn j(where Xlj ∈ X j for some space X j ) be the full data set for input model j, we batch the data into
{Xlj }l=(i−1)n j +1,...,in j , i = 1, . . . , k each of size n j . We call the collection of data {Xlj }l=(i−1)n j +1,...,in j , j =1, . . . , s the i-th section. Throughout our exposition, we ignore the issue of divisibility and assume forconvenience that the data size for each input model is divisible by k.
3.1 Conditions in the Parametric Regime
We first discuss the parametric case. Assume a parametric family for each input model, with a parameterjvector θ j ∈ Rd , j = 1, . . . , s. Denote d = ∑sj=1 d j as the dimension of the concatenated parameter vector.
Consider a simulatable performance measure ψ = ψ(θθ ) where θ = (θ 1 , . . . , θ s ). We make the following
assumption regarding the parameter estimators:
Assumption 1 (Normality of parameter estimators) We can construct an estimator θ̂θ i for θ using section
i. There is a scaling parameter m for the section sample size such that as m goes to ∞, we have√
( m(θ̂θ i − θ ))i=1,...,k ⇒ (Wi )i=1,...,ki.i.d.

where Wi ∼ N(0, Σ) for some covariance matrix Σ ∈ Rd .
To get a sense of the meaning of m and the quantity Σ, consider the simple example where eachindependent input model j is estimated from i.i.d. data of size kn j = kmp j , so that each section consists
j
j
j
j
dj
of
√ n =j mp observations. Suppose each parameter vector θ ∈ R is estimated with θ̂i such thatn j (θ̂i − θ j ) ⇒ N(0, Σ j ). This holds if we use, e.g., maximum likelihood estimator (MLE) under standard
conditions, in which case Σ j is the inverse of the Fisher information matrix (see, e.g., Serfling 2009 Section
4.2.2). In this case, if we fix p j and let m → ∞, Assumption 1 holds with Σ = diag((1/p j )Σ j , j = 1, . . . , s)approx.
where diag(·) denotes the block diagonal matrix with entries (·). In other words, we have θ̂θ i ∼N(θθ , diag((1/n j )Σ j , j = 1, . . . , k)) as n j grow large for all j.
Assumption 1 holds more generally for data exhibiting cross-sectional dependence (where one cansimply think of the collection of input models as a single giant input model), serially dependent data under
stationarity and suitable mixing conditions, and combinations of these aforementioned situations.
Next we make an assumption on the function ψ(·):
Assumption 2 (Parametric smoothness of ψ) The function ψ(·) has a non-zero differential at θ .
The differential defined for ψ in Assumption 2 means the function∇ψ(θθ )0 t, t ∈ Rd
where ∇ψ is the gradient of ψ, given that the partial derivatives of ψ all exist. As a direct consequenceof Serfling (2009) Section 1.12.2, the following is a handy sufficient condition for Assumption 2:
Lemma 3 Suppose all the partial derivatives of ψ exist in a neighborhood of θ and are continuous at θ .
Moreover, ∇ψ(θθ ) is not a zero vector. Then Assumption 2 holds.
Roughly speaking, the existence and continuity of partial derivatives in Lemma 3 hold when the input
densities are differentiable with respect to θ , and that an exchange of derivative and expectation (appearing
typically in defining ψ) is valid; general conditions of the latter can be found in, e.g., L’Ecuyer (1995).
Assumptions 1 and 2 together imply asymptotic normality via the delta method (see, e.g., Serfling2009 Section 3.3 Theorem A):
Lemma 4 Under Assumptions 1 and 2, we have√
( m(ψ(θ̂θ i ) − ψ(θθ )))i=1,...,k ⇒ (Zi )i=1,...,kGlynn and Lam
i.i.d.

where Zi ∼ N(0, σ 2 ) and σ 2 = ∇ψ(θθ )0 Σ∇ψ(θθ ), with ∇ψ being the gradient of ψ that exists under
Assumption 2.
β ) be an unbiased simulation
Next we state some assumptions about the simulation noise. Let ψ̂l (ββ )|ββ ] = ψ(ββ ). Denoterun (l is used as an index when running multiple replications) so that E[ψ̂l (ββ ) = Var(ψ̂l (ββ )|ββ ) as the simulation variance for each replication conditional on the input parameter.
τ 2 (ββ ) = E[|ψ̂l (ββ ) − ψ(ββ )|3 |ββ ] as the conditional third-order absolute central moment.
Similarly, denote κ3 (βdβ ∈ R : kββ − θ k2 < η} be an η-neighborhood of θ .
Let Nη (θθ ) = {βAssumption 5 (Parametric simulation noise) We have
1. τ 2 (·) is continuous at θ with 0 < τ 2 (θθ ) < ∞.
β ) is uniformly bounded for any β ∈ Nη (θθ ) for some η > 0.
2. κ3 (βSufficient conditions guaranteeing Assumption 5 are similar to the ones in Lemma 3, but needing ahigher order moment condition on the performance function that guarantees the existence of τ 2 and κ3 .
3.2 Conditions in the Nonparametric Regime
We now turn to the nonparametric regime. Instead of using/assuming parametric families for the input
distributions, we utilize the empirical distributions of the data as the input estimates. In this case, consider
a performance measure ψ = ψ(F) where F = (F 1 , . . . , F s ) denotes the vector of distribution functions, eachfor one of the s independent input models.
With collected data divided into k sections each of size n j for model j (like before), we constructF̂i = (F̂i1 , . . . , F̂is ), the vector of empirical distributions obtained from section i, wherejin1F̂i (·) = j
δ j (·)∑
n l=(i−1)n j +1 Xljand δX j (·) denotes the delta measure at Xlj .
lFor any linear map L = (L1 , . . . , Ls ) where L j : X j → R, and any G = (G1 , . . . , Gs ) where eachj
G : X j → R is a probability distribution, we defineshL , Gi = ∑Z

L j dG j .

j=1

Hence we havejin1
∑ j L j (Xlj ).
jnj=1
l=(i−1)n +1s

hL , F̂i i = ∑In contrast to the parametric case, here in the nonparametric setting we depict local assumptions (i.e.,
differentiability) in a probabilistic way. In other words, instead of separating the input estimation errors(similar to Assumption 1) and the smoothness property of ψ (similar to Assumption 2), we impose anassumption that enforces ψ applied on F̂i is linearizable with an error that is statistically small. This is
described as follows:
Assumption 6 (Nonparametric input noise) Assume ψ is differentiable in the sense that
ψ(F̂i ) = ψ(F) + hDψ, F̂i − Fi + ξi
where Dψ is a linear map that could depend on F̂i and F, such that there exists a scaling parameter m withthe property that as m → ∞,

Glynn and Lam
1.

√( mhDψ, F̂i − Fi)i=1,...,k ⇒ (Zi )i=1,...,ki.i.d.

where Zi ∼ N(0, σ 2 ) for some σ 2 < ∞.
2.

√pmξi → 0 for all i = 1, . . . , k.

The linear map Dψ in Assumption 6 can be viewed as a derivative with respect to the collection of
distribution functions (hence the notation). The reason why we use such a probabilistic assumption is to
avoid the introduction of a metric in defining the error in using F̂i to approximate F, which can lead to theneed of more regularity condition on ψ. Note that, in using the latter approach, the error term ξi depends
on both the chosen metric and the continuity of ψ with respect to this metric. When a weak norm is usedthen F̂i → F (e.g., the supremum norm, in which case F̂i → F as a consequence of the Glivenko-Cantellilemma), but more conditions are needed for the function ψ to be continuous (e.g., if ψ = E[h(X)], thenh needs to have bounded variation). The imposition of a probabilistic statement like Assumption 6 avoids
this dilemma, and in many situations a direct verification of this assumption can be done.
For example, consider a finite-horizon performance measure ψ(F) = EF [h(X j (t), j = 1, . . . , s, t =1, . . . , T j )] where h is some performance function, T j is the time horizon for each input sequence, X j (t)are i.i.d. replications generated from independent input model j and EF [·] denotes the expectation under
the collection of input distributions F. Suppose that
E[h(X j (ltj ), j = 1, . . . , s, t = 1, . . . , T j )2 ] < ∞(5)where ltj ≥ 1 can be any integers. Suppose also that the data size of each section for model j satisfies
n j = mp j for a scaling parameter m and p j > 0 are fixed constants. Then one can verify directly that
Assumption 6 holds.
Note that the parametric conditions in Section 3.1 can be translated into the nonparametric conditions inAssumption 6 under appropriate modifications. For instance, if ψ(θθ ) satisfies Assumption 2 and θ = θ (F)is differentiable with respect to F in a similar sense as Assumption 6, then ψ also satisfies Assumption 6
when viewed as a function of F.
Next we state some assumptions about the simulation noise. Let ψ̂l (·) be a conditionally unbiased
simulation replication for ψ(·), such that E[ψ̂l (G)|G] = ψ(G). Denote τ 2 (G) = Var(ψ̂l (G)|G) as theconditional variance, and for convenience we denote τ 2 = τ 2 (F). Similarly, denote κ3 (G) = E[|ψ̂l (G) −ψ(G)|3 |G] as the conditional third-order absolute central moment. We make the following assumption on
the simulation noise:
Assumption 7 (Nonparametric simulation noise) We havep

1. τ 2 (F̂i ) → τ 2 with 0 < τ 2 < ∞.
2. P(κ3 (F̂i ) ≤ M) → 1 for some fixed M > 0.
Assumption 7 holds for, e.g., the finite-horizon performance measure discussed above under a higher
order moment condition than in (5).
3.3 A Joint Central Limit Theorem
The following proposition describes a CLT held jointly between the input noise and the simulation noise:
Proposition 8 Consider a collection of estimators ν̂ν i , i = 1, . . . , k for ν that satisfies√( m(γ(ν̂ν i ) − γ(νν )))i=1,...,k ⇒ (Zi )i=1,...,kGlynn and Lam
i.i.d.

where Zi ∼ N(0, σ 2 ) for some function γ and σ 2 > 0. Consider another (possibly dependent) collectionof estimators ν̂ν 0u , u = 1, . . . , v. Let
1 ru
ψ̂(ν̂ν 0u ) = ∑ ψ̂l (ν̂ν 0u )ru l=1
where ru = mqu for some fixed qu > 0, ψ̂l (ν̂ν 0u ) are i.i.d. unbiased simulation runs for ψ(ν̂ν 0u ) given ν̂ν 0u ,β ) = Var(ψ̂l (ββ )|ββ ) as the conditional
and are conditionally independent across u = 1, . . . , v. Denote τ 2 (β22β) =simulation variance for each replication, and for convenience denote τ = τ (νν ). Similarly, denote κ3 (βpβ ) − ψ(ββ )|3 |ββ ] as the conditional third-order absolute central moment. Assume that τ 2 (ν̂ν 0u ) → τ 2
E[|ψ̂l (β0and P(κ3 (ν̂ν u ) ≤ M) → 1 for u = 1, . . . , v, as m → ∞.
Then,

√
√
( m(γ(ν̂ν i ) − γ(νν ))i=1,...,k , ( m(ψ̂(ν̂ν 0u ) − ψ(ν̂ν 0u )))u=1,...,v ⇒ ((Zi )i=1,...,k , (Wu )u=1,...,v )(6)i.i.d.

i.i.d.

where Zi ∼ N(0, σ 2 ) and Wu ∼ N(0, τ 2 /qu ) are independent.
The proposition applies whether ν is a real-valued parameter vector θ (parametric case) or the collectionof probability distributions F (nonparametric case), or a mix of both. To apply Proposition 8 to justify oursectioning scheme, we will see that the function γ(·) can be taken either as the performance function ψ(·) or its
plinearization. The assumption on the simulation replications, namely τ 2 (ν̂ν 0u ) → τ 2 and P(κ3 (ν̂ν 0u ) ≤ M) → 1,
can be satisfied by using a law of large numbers under the presented assumptions in the parametric caseand Assumption 7 directly in the nonparametric case.
Note that similar results hold if ru is not of the same order as m. Suppose ru = ω(m), then, under(otherwise) the same assumptions as Proposition 8, the limit in (6) can be replaced by ((Zi )i=1,...,k , (0)u=1,...,v ).
On the other hand, if ru = mpqu where p = o(m), then (6) is replaced by√√( mp(γ(ν̂ν i ) − γ(νν ))i=1,...,k , ( mp(ψ̂(ν̂ν 0u ) − ψ(ν̂ν 0u )))u=1,...,v ⇒ ((0)i=1,...,k , (Wu )u=1,...,v ).
The sectioning scheme we depict next works for these imbalanced cases with little modification.
4

SECTONING INPUT DATA FOR CONFIDENCE INTERVAL ESTIMATIONThis section presents our sectioning approach to construct CIs. Section 4.1 states our procedure and shows
the asymptotic confidence guarantee. Section 4.2 discusses other performance aspects of our proposed
procedure.
4.1 Procedure and Main Guarantee
We revisit our procedure discussed in Section 2 using our introduced notations. To estimate ψ = ψ(νν )(where ν can be θ or F), we divide the data set of each input model into k sections of equal length. For
section i that consists of {Xlj }l=(i−1)n j +1,...,in j , j = 1, . . . , s, we estimate ν̂ν i (can be θ̂θ i or F̂i ), and thenYi = ψ̂(ν̂ν i ) where
1 rψ̂(ν̂ν i ) = ∑ ψ̂l (ν̂ν i ).
r l=1Here ψ̂l (ν̂ν i ) are conditionally unbiased i.i.d. simulation replications for ψ(ν̂ν i ), and independent across i.
Let Ȳ = (1/k) ∑ki=1 Yi and S2 = (1/(k − 1)) ∑si=1 (Yi − Ȳ )2 be the sample mean and variance of the k sectioned
estimates. We propose to use the CISSIsec = Ȳ − tk−1,1−α/2 √ , Ȳ + tk−1,1−α/2 √(7)kkGlynn and Lam
where tk−1,1−α/2 is the (1 − α/2)-quantile of the tk−1 distribution. This construction holds as long as k ≥ 2.
We have the following guarantees for (7). For this, we introduce the scaling r = mq, where q is fixed
while m goes to ∞. We have:
Theorem 9 Under either:
1. (parametric case) Assumptions 1, 2 and 5; or
2. (nonparametric case) Assumptions 6 and 7,
we haveȲ − ψ
√ ⇒ tk−1
S/ k(8)

for any fixed k ≥ 2. Consequently, Isec is an asymptotically exact (1 − α)-level CI for ψ, i.e.,P(ψ ∈ Isec ) → 1 − α as m → ∞.
√Proof.
We first develop a limit theorem for m(Yi − ψ) for both the parametric and the nonparametriccase. Then we proceed to justify the t-asymptotic (8) and conclude the theorem.
Parametric case: We can write
Yi = ψ +Ui +Vi
where Ui = ψ(θ̂θ i ) − ψ and Vi = ψ̂(θ̂θ i ) − ψ(θ̂θ i ). Under Assumptions 1 and 2, we use Lemma 4 to deduce√
i.i.d.
that m(Ui )i=1,...,k ⇒ (Qi )i=1,...,k where Qi ∼ N(0, σ 2 ) with σ 2 = ∇ψ(θθ )0 Σ∇ψ(θθ ). Using the notationsin Assumption 5, a law of large numbers for θ̂θ i implied from Assumption 1 and the local properties inpAssumption 5 conclude that τ 2 (θ̂θ i ) → τ 2 (θθ ) and P(κ3 (θ̂θ i ) ≤ M) → 1 for some M > 0. Hence Proposition√
i.i.d.
8 implies that m((Ui )i=1,...,k , (Vi )i=1,...,k ) ⇒ ((Qi )i=1,...,k , (Wi )i=1,...,k ) where Wi ∼ N(0, τ 2 /q) and Qi and
Wi are all independent. Therefore,√m(Yi − ψ)i=1,...,k ⇒ (Zi )i=1,...,ki.i.d.

where Zi ∼ N(0, σ 2 + τ 2 /q).
Nonparametric case: Under Assumption 6, we write
Yi = ψ +Ui + ξi +Vi
where Ui = hDψ, F̂i − Fi, ξi defined as in Assumption 6, and Vi = ψ̂(F̂i ) − ψ(F̂i ). Assumption 6 states that√√pi.i.d.
( mUi )i=1,...,k ⇒ (Qi )i=1,...,k where Qi ∼ N(0, σ 2 ) for some σ 2 < ∞ and mξi → 0 for all i = 1, . . . , k. Notepthat Assumption 7 enforces τ 2 (F̂i ) → τ 2 (F) and P(κ3 (F̂i ) ≤ M) → 1 for some M > 0. Hence, Proposition√
i.i.d.
8 implies that m((Ui )i=1,...,k , (Vi )i=1,...,k ) ⇒ ((Qi )i=1,...,k , (Wi )i=1,...,k ) where Wi ∼ N(0, τ 2 /q) and Qi and
Wi are all independent. Therefore, by Slutsky’s theorem,√m(Yi − ψ)i=1,...,k ⇒ (Zi )i=1,...,ki.i.d.

where Zi ∼ N(0, σ 2 + τ 2 /q).
Finally, in both cases
Ȳ − ψZ̄√ ⇒q
S/ k(1/(k − 1)) ∑ki=1 (Zi − Z̄)2 /k(9)Glynn and Lam
where Z̄ = (1/k) ∑ki=1 Zi , by the continuous mapping theorem. Note that, since Zi are Gaussian, Z̄ is2
2
is a χ 2 random variable with degree of
where χk−1
independent of Zi − Z̄, and ∑ki=1 (Zi − Z̄)2 ∼ χk−1
freedom k − 1. This shows that the limit in (9) is a tk−1 random variable. The last assertion in the theorem
follows by a standard pivotal argument.
4.2 Conservativeness and Variability of the Confidence Intervals
We discuss the performance of our sectioning-based CIs. Ideally, one would use as many sections as
possible, as this would improve the length of the interval (in terms of reducing the average length and also
its variability); in the limiting case, this would give rise to a CI constructed from normality asymptotic,
which can be viewed as the best possible. However, with fixed data and simulation budgets, using more
sections means a smaller data size and simulation replications per section, so that the quality of the Gaussian
approximation within each section, and consequently the performance of the pivotal t-statistic, deteriorates.
The number of sections and the size per section thus result in a tradeoff between desirable half-width
properties and the statistical validity of the CIs.
We use the mean half-width to measure the conservativeness of the CI, and its standard deviation tomeasure its variability. In the following, we do a rough-cut analysis on the effect of section numbers on these
two measurements, which reveals the minimal section numbers needed to achieve a certain performancecompared to the benchmark normality CI. A full investigation of the optimal section number, however,would require analyzing the approximation error of the involved asymptotic t-limit.
Suppose that we use the full data set and all the simulation budget to get an estimate ψ̂(ν̂ν ) where ν̂ν isthe input estimate and ψ̂ denotes the average of all simulation runs. Then, under the same assumption as
in Theorem 9 (with k set to 1), we have, asymptotically, ψ̂(ν̂ν ) ∼ N(ψ, σ02 + τ02 ) for some σ02 and τ02 , whereσ02 denotes the variance from the input noise and τ02 the simulation noise. For convenience, we denote
λ 2 = σ02 + τ02 , so that ψ̂(ν̂ν ) ∼ N(ψ, λ 2 ) approximately. If we instead divide the data and the simulationbudget into k sections, then for each section i our estimate becomes ψ̂(ν̂ν i ), and they are asymptoticallyi.i.d.

∼ N(ψ, kλ 2 ).
√The half-width of Isec in (7) is tk−1,1−α/2 S/ k. Here S2 , the sample variance of the ψ̂(ν̂ν i )’s, follows
2 /(k − 1) approximately. Therefore, the half-width H of I
kλ 2 · χk−1
sec isss22
tk−1,1−α/2χχk−1√· kλ 2 k−1 = tk−1,1−α/2 λ
H≈
.
k−1
k−1
k
q
√
2 = 2Γ(k/2)/Γ((k − 1)/2), the expected half-width satisfies
Using the fact that E χk−1
Γ(k/2)
EH ≈ tk−1,1−α/2 λΓ((k − 1)/2)

r

2
.
k−1

2 = k − 1, the second moment of the half-width satisfies
On the other hand, using Eχk−1
2
EH 2 ≈ tk−1,1−α/2λ2

so that the variance is2
Var(H) ≈ tk−1,1−α/2
λ2

Γ(k/2)
1−
Γ((k − 1)/2)

2

2
k−1

!
.

p
It can be shown that Γ(k/2)/Γ((k − 1)/2) 2/(k − 1) → 1 as k → ∞. Combined with the fact thattk−1 ⇒ N(0, 1), the approximate EH → z1−α/2 λ , where z1−α/2 is the 1 − α/2 standard normal quantile,

Glynn and Lam
and the approximate Var(H) → 0. For example, when α = 5% (so that we are looking at the 95% CI),
EH/λ ≈ 1.98 at k = 100, compared with z1−α/2 = 1.96. Note that the quantity z1−α/2 λ is the half-widthof a normality CI. Therefore, as discussed above, the normality CI can be viewed as the best possible (Ofcourse, to get this, one would need information about λ , which as discussed in the introduction can bechallenging to obtain in the input uncertainty setting).
The approximate EH, EH 2 and Var(H) all monotonically decrease as k increases. Fixing α = 5%,
the decreasesof all the considered quantities occur rapidly from k = 2 to 6 (EH/λ ≈ 10.1, 3.8, 2.9, 2.6, 2.4
pand Var(H)/λ ≈ 7.7, 2.0,
p1.2, 0.9, 0.8 at k = 2 to 6), continues to decrease at a slower rate from k = 6 to
10 (e.g., EH/λ ≈ 2.2 andp Var(H)/λ ≈ 0.5 at k = 10), and afterwards the gain becomes quite negligible(e.g., EH/λ ≈ 2.0 and Var(H)/λ ≈ 0.2 at k = 50). These numbers suggest that increasing k from 2 to 3
leads to a significant gain, whereas increasing k further to 4 to 6 can improve sensibly further. Dependingon the data/simulation size, a reasonable k could be chosen in the range 3 to 6.
5

A NUMERICAL EXAMPLE

We test our approach with a numerical example of M/M/1 queue. In this queue, the true arrival rate is 1
and service rate is 1.2. We initiate from an empty system. We consider two performance measures: theexpected waiting time averaged over the first 20 customers, denoted as Ewaittime , and the expected timeaveraged queue length over time period [0, 5], denoted as Equeuelength . The true values of Ewaittime = 1.53
and Equeuelength = 0.63.
We use a data size 100 on both the interarrival and service times (i.e., these data are sampled fromthe exponential distributions with rates 1 and 1.2 respectively). We execute our sectioning scheme toconstruct 95% CIs for measures Ewaittime and Equeuelength , using various numbers of sections, but under atotal simulation budget of 100 runs. In each setting, we repeat the calculation 1000 times, where each time
we sample a new set of observations and construct a CI. From these repetitions, we estimate the coverage
probability and the average CI length. We consider the nonparametric setting (i.e., using the empirical
distributions for the interarrival and service times).
Table 1 shows the estimated coverage probabilities (both the point estimate and the 95% confidence
interval) and the average CI lengths in estimating Ewaittime and Equeuelength as the number of sections k varies(e.g., when k = 2, each section constitutes 50 observations for both the interarrival and service times and
50 replications). The coverages of both performance measures seem to be valid and quite stable until aboutk = 14, after which they start to deteriorate from the nominal 95%. This shows, at least in this example,
that the coverage probability is relatively robust to the data size in each section. On the other hand, theaverage CI length drops significantly from k = 2 to 3, continues to drop from k = 3 to 6, and afterwards
levels off. This shows a huge gain in reducing conservativeness from 2 to 3 sections, and further gains till
about 6 sections. These observations reconcile with our discussion in Section 4.2.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020
and CAREER CMMI-1653339/1834710.
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Stochastic Simulation: Algorithms and Analysis, Volume 57. Springer
Science & Business Media, New York.
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., 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”.

Glynn and Lam
Table 1: Nonparametric CIs using 100 observations on the interarrival times and 100 observations on theservice times, with total simulation budget 100.
Section #Coverage (wait time)

CI length (wait time)

Coverage (queue length)

CI length (queue length)

2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

0.940 (0.925, 0.955)
0.943 (0.929, 0.957)
0.928 (0.912, 0.944)
0.949 (0.935, 0.963)
0.930 (0.914, 0.946)
0.941 (0.926, 0.956)
0.942 (0.927, 0.957)
0.936 (0.921, 0.951)
0.953 (0.940, 0.966)
0.947 (0.933, 0.961)
0.971 (0.961, 0.981)
0.947 (0.933, 0.961)
0.947 (0.933, 0.961)
0.945 (0.931, 0.959)
0.924 (0.908, 0.940)
0.936 (0.921, 0.951)
0.928 (0.912, 0.944)
0.915 (0.898, 0.932)
0.922 (0.905, 0.939)

10.024
3.995
3.034
2.745
2.558
2.553
2.451
2.423
2.430
2.362
2.465
2.513
2.397
2.551
2.429
2.628
2.512
2.454
2.398

0.936 (0.921, 0.951)
0.947 (0.933, 0.961)
0.946 (0.932, 0.960)
0.952 (0.939, 0.965)
0.945 (0.931, 0.959)
0.941 (0.926, 0.956)
0.960 (0.948, 0.972)
0.939 (0.924, 0.954)
0.949 (0.935, 0.963)
0.950 (0.936, 0.964)
0.963 (0.951, 0.975)
0.943 (0.929, 0.957)
0.949 (0.935, 0.963)
0.936 (0.921, 0.951)
0.928 (0.912, 0.944)
0.918 (0.901, 0.935)
0.908 (0.890, 0.926)
0.900 (0.881, 0.919)
0.900 (0.881, 0.919)

3.704
1.485
1.136
1.046
1.011
0.973
0.995
0.980
0.973
0.989
1.019
1.098
1.006
1.122
1.080
1.164
1.153
1.100
1.080

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.
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 ofthe 2001 Winter Simulation Conference, edited by B. A. Peters et al., Volume 1, 372–378. Piscataway,
New Jersey: IEEE.
Biller, B., and C. G. Corlu. 2011. “Accounting for Parameter Uncertainty in Large-Scale Stochastic
Simulations with Correlated Inputs”. Operations Research 59(3):661–673.
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. 2004. “Calculation of Confidence Intervals for Simulation Output”. ACM
Transactions on Modeling and Computer Simulation 14(4):344–362.
Chick, S. E. 2001. “Input Distribution Selection for Simulation Experiments: Accounting for Input Uncertainty”. Operations Research 49(5):744–758.
Chick, S. E. 2006. “Bayesian Ideas and Discrete Event Simulation: Why, What and How”. In Proceedingsof the 2006 Winter Simulation Conference, edited by L. F. Perrone et al., 96–106. Piscataway, New
Jersey: IEEE.
Glynn, P. W. 1986. “Problems in Bayesian Analysis of Stochastic Simulation”. In Proceedings of the 1986
Winter Simulation Conference, edited by J. Wilson et al., 376–379. New York, New York: ACM.
Glynn, P. W., and D. L. Iglehart. 1990. “Simulation Output Analysis Using Standardized Time Series”.
Mathematics of Operations Research 15(1):1–16.

Glynn and Lam
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 Jersey: IEEE.
L’Ecuyer, P. 1995. “Note: On the Interchange of Derivative and Expectation for Likelihood Ratio Derivative
Estimators”. Management Science 41(4):738–747.
Lin, Y., E. Song, and B. Nelson. 2015. “Single-Experiment Input Uncertainty”. Journal of Simulation 9(3):249–259.
Muñoz, D. F., and P. W. Glynn. 1997. “A Batch Means Methodology for Estimation of a Nonlinear Functionof a Steady-State Mean”. Management Science 43(8):1121–1135.
Nakayama, M. K. 2014. “Confidence Intervals for Quantiles Using Sectioning When Applying VarianceReduction Techniques”. ACM Transactions on Modeling and Computer Simulation (TOMACS) 24(4):19.
Serfling, R. J. 2009. Approximation Theorems of Mathematical Statistics, Volume 162. John Wiley & Sons.
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.
Zouaoui, F., and J. R. Wilson. 2004. “Accounting for Input-Model and Input-Parameter Uncertainties inSimulation”. IIE Transactions 36(11):1135–1151.
AUTHOR BIOGRAPHIESPETER W. GLYNN is the Thomas Ford Professor in the Department of Management Science and Engineering (MS&E) at Stanford University. He is a Fellow of INFORMS and of the Institute of Mathematical
Statistics, has been co-winner of Best Publication Awards from the INFORMS Simulation Society in
1993, 2008, and 2016, and was the co-winner of the John von Neumann Theory Prize from INFORMS
in 2010. In 2012, he was elected to the National Academy of Engineering. His research interests lie instochastic simulation, queueing theory, and statistical inference for stochastic processes. His email address
is glynn@stanford.edu.
HENRY 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.