$\Lambda^{*}=\begin{bmatrix}\Lambda_{\mathcal{J}}&0\\\ 0&0\end{bmatrix}$ and
$\Lambda_{\mathcal{J}}$ denotes sub-matrix of $\Lambda$ with respect to subset
$\mathbf{X}_{\mathcal{J}}$.

Since the condition of Theorem 13.1 in [28] still holds, we have

$\sqrt{m}(m_{p}({\bf
X}^{*}_{\mathcal{J}})-m_{p}(\mathbf{x}_{c}))\stackrel{{\scriptstyle
D}}{{\rightarrow}}\mbox{N}(0,\sigma^{2}_{\mathcal{J}})$

where $\sigma_{\mathcal{J}}^{2}=\nabla_{\mathcal{J}}m_{p}({\bf
x}_{c})^{\top}\Lambda_{\mathcal{J}}\nabla_{\mathcal{J}}m_{p}({\bf x}_{c})>0$,
$\nabla_{\mathcal{J}}$ is gradient with respect to subset
$\mathbf{X}_{\mathcal{J}}$ This establishes the constant. Moreover, the
Condition (14) of Theorem 3.8 of [29] holds. Thus, the cost function estimator
$\widehat{c}(\mathcal{J})$ is strongly consistent:
$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{c}(\mathcal{J})=\sigma_{\mathcal{J}}^{2}$
almost surely.

Similarly, we can prove the cost function estimator
$\widehat{c}(\mathcal{J}\cup\\{\ell\\})$ is also strongly consistent:
$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{c}(\mathcal{J}\cup\\{\ell\\})=\sigma_{\mathcal{J}\cup\\{\ell\\}}^{2}$
almost surely, where
$\sigma_{\mathcal{J}\cup\\{\ell\\}}^{2}=\nabla_{\mathcal{J}\cup\\{\ell\\}}m_{p}({\bf
x}_{c})^{\top}\Lambda_{\mathcal{J}\cup\\{\ell\\}}\nabla_{\mathcal{J}\cup\\{\ell\\}}m_{p}({\bf
x}_{c})>0$, $\nabla_{\mathcal{J}\cup\\{\ell\\}}$ is gradient with respect to
$\mathbf{X}_{\mathcal{J}\cup\\{\ell\\}}$, and
$\Lambda_{\mathcal{J}\cup\\{\ell\\}}$ is sub-matrix of $\Lambda$ with respect
to $\mathbf{X}_{\mathcal{J}\cup\\{\ell\\}}$.

Consequently, we have Sharpley Value estimator $\widehat{s}_{\ell}$ is
strongly consistent:

$\lim_{m\rightarrow\infty}\lim_{B\rightarrow\infty}m\widehat{s}_{\ell}=\sigma_{s}^{2},~{}~{}\text{with}~{}~{}\sigma_{s}^{2}=\sum_{\mathcal{J}\subset\mathcal{L}/\\{\ell\\}}\dfrac{(L-|\mathcal{J}|-1)!|\mathcal{J}|!}{L!}\left[\sigma_{\mathcal{J}\cup\\{\ell\\}}^{2}-\sigma_{\mathcal{J}}^{2}\right]$

almost surely.

Next, we will show:

$\displaystyle\lim_{m\rightarrow\infty}ms_{\ell}$
$\displaystyle=\sigma_{s}^{2}$

First, we need to show
$\lim_{m\rightarrow\infty}mc(\mathcal{J})=\sigma_{\mathcal{J}}^{2}$. Let
$L({\bf X}^{*}_{\mathcal{J}},{\bf x}_{c})$ denote the line segment joining
${\bf X}^{*}_{\mathcal{J}}$ and ${\bf x}_{c}$. According to the Multivariate
Taylor Formula [27],

$\displaystyle m_{p}({\bf X}^{*}_{\mathcal{J}})=m_{p}({\bf x}_{c})+\nabla
m_{p}({\bf x}_{c})^{\top}({\bf X}^{*}_{\mathcal{J}}-{\bf x}_{c})$
$\displaystyle+\frac{1}{2}({\bf X}^{*}_{\mathcal{J}}-{\bf
x}_{c})^{\top}\nabla^{2}m_{p}({\bf x}_{c})({\bf X}^{*}_{\mathcal{J}}-{\bf
x}_{c})+\mathcal{R}({\bf X}^{*}_{\mathcal{J}},{\bf x}_{c}).$

And the remainder term

$\mathcal{R}({\bf X}^{*}_{\mathcal{J}},{\bf
x}_{c})=\frac{1}{3!}\sum_{i_{1}=1}^{d_{\mathcal{J}}}\sum_{i_{2}=1}^{d_{\mathcal{J}}}\sum_{i_{3}=1}^{d_{\mathcal{J}}}\left.\frac{\partial^{3}m_{p}(x_{1},\ldots,x_{d})}{\partial
x_{i_{1}}\partial x_{i_{2}}\partial x_{i_{3}}}\right|_{{\bf x}={\bf
z}}\prod_{j=1}^{3}(X_{\mathcal{J},i_{j}}^{*}-x_{c,i_{j}})$

where $d_{\mathcal{J}}=\sum_{\ell\in\mathcal{J}}h_{\ell}$, ${\bf z}$ denotes a
value in the interior of $L({\bf X}^{*}_{\mathcal{J}},{\bf x}_{c})$, and
$X_{\mathcal{J},i}$ and $x_{c,i}$ denote the $i$th components of the vectors
${\bf X}^{*}_{\mathcal{J}}$ and ${\bf x}_{c}$.

Following the same procedure of Theorem 4.2.’s proof , we can obtain

$[\mbox{E}(m_{p}({\bf X}^{*}_{\mathcal{J}})]^{2}=m_{p}^{2}({\bf
x}_{c})+m_{p}({\bf x}_{c})\mbox{E}[({\bf X}^{*}_{\mathcal{J}}-{\bf
x}_{c})^{\top}\nabla^{2}m_{p}({\bf x}_{c})({\bf X}^{*}_{\mathcal{J}}-{\bf
x}_{c})]+O\left(\frac{1}{m^{2}}\right).$

and

$\displaystyle\mbox{E}[m_{p}^{2}({\bf X}^{*}_{\mathcal{J}})]$
$\displaystyle=m_{p}^{2}({\bf x}_{c})+\mbox{E}\Big{[}({\bf
X}^{*}_{\mathcal{J}}-{\bf x}_{c})^{\top}\nabla m_{p}({\bf x}_{c})\nabla
m_{p}({\bf x}_{c})^{\top}({\bf X}^{*}_{\mathcal{J}}-{\bf x}_{c})$
$\displaystyle~{}~{}+({\bf X}^{*}_{\mathcal{J}}-{\bf x}_{c})^{\top}m_{p}({\bf
x}_{c})\nabla^{2}m_{p}({\bf x}_{c})({\bf X}^{*}_{\mathcal{J}}-{\bf
x}_{c})\Big{]}+O\left(\frac{1}{m^{2}}\right).$

Then,

$\displaystyle\mbox{Var}[m_{p}({\bf
X}^{*}_{\mathcal{J}})]=\mbox{E}[m_{p}^{2}({\bf
X}^{*}_{\mathcal{J}})]-\Big{(}\mbox{E}[m_{p}({\bf
X}^{*}_{\mathcal{J}})]\Big{)}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{m}\nabla m_{p}({\bf x}_{c})^{\top}\Lambda^{*}\nabla
m_{p}({\bf x}_{c})+O\left(\frac{1}{m^{2}}\right).$ $\displaystyle=$
$\displaystyle\frac{1}{m}\nabla_{\mathcal{J}}m_{p}({\bf
x}_{c})^{\top}\Lambda_{\mathcal{J}}\nabla_{\mathcal{J}}m_{p}({\bf
x}_{c})+O\left(\frac{1}{m^{2}}\right).$

Therefore, we have
$\lim_{m\rightarrow\infty}mc(\mathcal{J})=\sigma_{\mathcal{J}}^{2}$.
Similarly, we obtain
$\lim_{m\rightarrow\infty}mc(\mathcal{J}\cup\\{\ell\\})=\sigma_{\mathcal{J}\cup\\{\ell\\}}^{2}$.
Finally, by applying continuous mapping theorem, we have
$\lim_{m\rightarrow\infty}ms_{\ell}=\sigma_{s}^{2}$

Remark: The independent variables in our stochastic kriging metamodel consist
of central moments and standardized central moments, rather than raw moments.
However, Theorem 3 can easily be extended to central and standardized central
moments as follows.

Since standardized moments are continuous functions of raw moments, denoted
generically as $g(\cdot)$, we can consider the composite function $(m_{p}\circ
g)(\cdot)$ and follow steps analogous to those in the proof of Theorem 3. Up
to the third derivatives we have

$\displaystyle(m_{p}\circ g)^{\prime}(t)$ $\displaystyle=$ $\displaystyle
m_{p}^{\prime}(g(t))g^{\prime}(t)$ $\displaystyle(m_{p}\circ
g)^{\prime\prime}(t)$ $\displaystyle=$ $\displaystyle
m_{p}^{\prime\prime}(g(t))[g^{\prime}(t)]^{2}+m_{p}^{\prime}(g(t))g^{\prime\prime}(t)$
$\displaystyle(m_{p}\circ g)^{(3)}(t)$ $\displaystyle=$ $\displaystyle
m_{p}^{(3)}(g(t))[g^{\prime}(t)]^{3}+2m_{p}^{\prime\prime}(g(t))g^{\prime}(t)g^{\prime\prime}(t)$
$\displaystyle+m_{p}^{\prime\prime}(g(t))g^{\prime}(t)g^{\prime\prime}(t)+m_{p}^{\prime}(g(t))g^{(3)}(t).$

Let $u$ denote the mean, $u_{i}^{\prime}$ denote the $i$th order raw moment
and $u_{i}$ denote the $i$th order central moment. Then the first three
central moments can be expressed as functions of raw moments as follows:

$\displaystyle u_{1}$ $\displaystyle=$ $\displaystyle u,$ $\displaystyle
u_{2}$ $\displaystyle=$ $\displaystyle u_{2}^{\prime}-u^{2}$ $\displaystyle
u_{3}$ $\displaystyle=$ $\displaystyle
u_{3}^{\prime}-3uu_{2}^{\prime}+2u^{3}.$

The first three standardized central moments are $u_{1},\sqrt{u_{2}}$ and
$u_{3}/u_{2}^{3/2}$. For a non-degenerate distribution, the second central
moment is positive and bounded away from $0$. Thus, the first three
derivatives $g^{\prime},g^{\prime\prime},g^{(3)}$ exist and are finite.

##  Appendix D Experiment Design

To fit SK metamodels we recommend the experiment design developed in [3] which
demonstrated robust performance over a number of test examples. In this
section, we briefly review the basic methodology; for detailed information
please refer to [3].

The experiment design is not specified a priori; instead the design space,
denoted by $\mathcal{D}$, depends on the real-world data
$\mathbf{z}_{\mathbf{m}}^{(0)}$ that will eventually be resampled. In this way
the design is adaptive.

At a high level, this is the approach: Generate a large number of bootstrap
samples from the real-world data $\mathbf{z}_{\mathbf{m}}^{(0)}$ and compute
the corresponding sample moments. Find a regular region that encompasses a
large fraction of this sample; this will be the design space. Generate
additional bootstrap samples to test that the regular region does indeed cover
the desired fraction of the feasible space of sample moments, and refine if
necessary. Once satisfied, embed a space-filling design in the regular region.
These design points correspond to input distribution moments at which to run
simulation experiments to fit the SK metamodel. We provide some more details
below.

Suppose we are interested in a $(1-\alpha)100\%$ CI; we set $\alpha=0.05$ in
our empirical study. We want the experiment design to lead to a metamodel that
is accurate for moments ${\bf x}$ that are the most likely bootstrap moment
vectors generated from $\mathbf{z}_{\mathbf{m}}^{(0)}$; by “likely” we mean,
for instance, covering $q=99\%>(1-\alpha)100\%=95\%$ of the feasible bootstrap
moments.

To this end we find an ellipsoid that will contain an independent bootstrap
moment vector obtained by random sampling from $\mathbf{z}_{\mathbf{m}}^{(0)}$
with probability at least $q$. We then generate a space-filling experiment
design inside this ellipsoid. The procedure for constructing the design is as
follows:

  1. 1.

Generate $B_{0}$ bootstrap resamples from $\mathbf{z}_{\mathbf{m}}^{(0)}$ and
compute the corresponding sample moments to generate a set of sample moments
$D_{T}=\\{\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)},b=1,2,\ldots,B_{0}\\}$.

  2. 2.

Find the smallest ellipsoid $E$ such that it contains the fraction $q$ of the
data in $D_{T}$ when the ellipsoid’s center and shape are the sample mean and
covariance matrix, respectively, of the elements of $D_{T}$.

  3. 3.

Perform a hypothesis test where the null hypothesis is that a bootstrap moment
will be contained in this ellipsoid with probability at least $q$. This
requires computing the number of bootstrap moment resamples, denoted by
$B_{1}$, and the constant $c$ that defines the rejection region to attain the
desired Type I error and power for the test.

  4. 4.

Generate $B_{1}$ additional independent bootstrap resamples from
$\mathbf{z}_{\mathbf{m}}^{(0)}$ and compute the moments
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)},b=B_{0}+1,B_{0}+2,\ldots,B_{0}+B_{1}$.
If more than $c$ of these $B_{1}$ resamples are contained in the ellipsoid,
then accept the current $E$ as the design space. Otherwise, add these
bootstrap resamples to $D_{T}$, let $B_{0}\leftarrow B_{0}+B_{1}$ and go to
Step 2 to update the ellipsoid.

  5. 5.

Generate $k$ space-filling design points in the ellipsoid $E$. To place design
points into this space, we employ an algorithm due to [32], §3.2.1, for
generating points uniformly distributed in an ellipsoid. The algorithm first
generates the polar coordinates of a point uniformly distributed in a
hypersphere, then transforms it to Cartesian coordinates, and finally
transforms it again to a point uniformly distributed in an ellipsoid. The
advantage of this approach is that each element of the initial polar
coordinates are independently distributed, allowing them to be generated
coordinate by coordinate via their inverse cumulative distribution function.
Rather than use randomly chosen points, however, we begin with a Latin
hypercube sample on $(0,1)^{d}$.

  6. 6.

Assign $n=N/k$ replications to each design point, where $N$ denotes total
computational budget. Together the transformed Latin hypercube design points
and the number of replications $n$ define the experiment design $\mathcal{D}$.

In our experiments we set Type I error of the hypothesis test to $0.005$ and
its power to $0.95$ when the true probability is $q=0.97$.

##  Appendix E Queueing Network Example

In this section we use an queueing network example in Figure 2 to evaluate the
performance of our uncertainty analysis framework. Consider estimating the
steady state expected number of customers in this network. The interarrival
times follow a gamma distribution, $A\sim\mbox{gamma}(\alpha_{A},\beta_{A})$,
and the service times at the $i$th station also follow a gamma distribution,
$S_{i}\sim\mbox{gamma}(\alpha_{S_{i}},\beta_{S_{i}})$. Customers finishing
service at stations $1,2,3$ must make decisions about their next station.
These routing decisions follow Bernoulli distributions
$P_{i}\sim\mbox{Ber}(p_{i}),i=1,2,3$. The parameters of the input
distributions, $\alpha_{A},\beta_{A}$,
$\\{(\alpha_{S_{i}},\beta_{S_{i}}),i=1,2,3,4\\}$ and $\\{p_{i},i=1,2,3\\}$ are
all unknown and estimated from real-world data. Our goal is to estimate the
steady-state expected number of customers in the system when the input
parameters assume their true but unknown values.

Figure 2: Queueing network example.

Both interarrival and service times follow gamma distributions and the routing
decisions follow Bernoulli distributions. Thus, it is a 13-dimensional problem
with $L=8$ input processes that include both continuous and discrete
distributions. The true model parameters are $\alpha_{A}=1,\beta_{A}=0.25$,
$\alpha_{S_{i}}=1$, $\beta_{S_{i}}=0.2$ for $i=1,2,3,4$ and
$p_{1}=p_{2}=0.5,p_{3}=0.75$. These parameter values imply a tractable Jackson
network with steady-state number of customers in system $\mu({\bf
x}_{c})=12.67$.

In the experiments we assume that all input model parameters are unknown and
are estimated from a finite sample of real-world data. Notice that
$\alpha_{A},\beta_{A}$, $\alpha_{S_{i}}$, $\beta_{S_{i}}$ for $i=1,2,3,4$ are
estimated from continuous measurements, while the routing probabilities
$p_{1},p_{2},p_{3}$ are estimated from 0 or 1 observations that would
correspond to customer routing decisions. The model with estimated input
parameters is almost surely not a Jackson network and it could be unstable.
Our measure of uncertainty is a $95\%$ CI for $\mu({\bf x}_{c})$ as defined by
(3) because this is the objective desired in practice.

To evaluate the robustness of the metamodel-assisted bootstrapping approach,
we systematically examine the effect of the quantity of real-world data and
the number of design points and replications per design point used to fit the
metamodel; We consider a wide range for the quantity of real-world data
$m=50,500,5000$, letting $m_{\ell}=m$ for $\ell=1,2,\ldots L$. The levels for
the number of design points are $k=20,40,80,130$. For a 13-dimensional problem
$k=20$ is a very small design. The studies by [14] and [19] recommend that the
number of design points should be $10$ times the dimension of the problem for
kriging; we take this as the maximum number of design points. The same number
of replications are assigned to all design points and we try $n=10,50,100$.

[3] demonstrated that $\mbox{CI}_{0}$ has good performance when the impact of
metamodel or simulation uncertainty is negligible. In this empirical study we
focus on situations where metamodel uncertainty may be significant. However,
rather than creating a problem that actually takes hours or days to run, we
instead construct a problem with high metamodel uncertainty by using short run
lengths for each replication: 20 time units after the warm up, which is
roughly equivalent to 80 finished customers. To avoid the influence from
initial bias, all simulations start loaded with the number of customers at
each station being their steady-state expected values (rounded) under ${\bf
x}_{c}$. Furthermore, a long warmup period of 200 time units is used. The net
effect is that the point estimators of the steady-state number in the network
have low bias, but may be quite variable.

###  E.1 Performance of CIs

A fundamental assumption of simulation is that the expectation $\mu({\bf
x}_{c})$ exists. This assumption does not imply, however, that it exists for
all possible values of ${\bf x}$, ${\bf X}_{\bf m}$ or $\widehat{{\bf X}}_{\bf
m}^{(b)}$ that might be realized. The prototype example is a congestion-
related performance measure of a queueing system as time goes to infinity when
congestion increases without bound for some values of its interarrival-time
and service-time parameters. We refer to systems for which $\mu({\bf x})$ is
$\pm\infty$ for some values of ${\bf x}$ as potentially unstable. The
conditional probability that a bootstrap resampled moment
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}$ is located in the unstable region,
denoted by $U$, given the real-world data is

$P_{U}\equiv\mbox{Pr}\left\\{\left.\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}\in
U\right|\mathbf{z}_{\mathbf{m}}^{(0)}\right\\}.$ (34)

Since $P_{U}$ only depends on $m$ and ${\bf x}_{c}$, we ran a side experiment
to estimate it using

$\widehat{P}_{U}=\frac{1}{B}\sum_{b=1}^{B}\mbox{I}\left(\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}\in
U\right),$ (35)

where $\mbox{I}(\cdot)$ is the indicator function. The means and standard
deviations (SD) of $\widehat{P}_{U}$ for $m=50,500,5000$ were estimated based
on 1000 macro-replications and are displayed in Table 6. In each macro-
replication we independently generated a sample of size $m$ of “real-world
data.” Then, conditional on these data, we drew $B=2000$ bootstrap resampled
moments.

Table 6: Percentage of unstable bootstrap resampled moments. | $m=50$ | $m=500$ | $m=5000$
---|---|---|---
mean of $\widehat{P}_{U}$ | 44.4% | 2.3% | 0
SD of $\widehat{P}_{U}$ | 31.7% | 7.9% | 0

As $m$ increases the bootstrap resampled moments become more closely centered
around ${\bf x}_{c}$. Thus, both the mean and SD of $\widehat{P}_{U}$ decrease
with increasing $m$ as shown in Table 6. When $m=50$, ${P}_{U}$ appears to be
much larger than $\alpha/2$ so the bootstrap moments
$\widehat{\mathbf{X}}_{\mathbf{m}}^{(b)}$ that correspond to the upper
confidence bound are located in the unstable region $U$ with high probability.
When $m=500$, ${P}_{U}$ appears to be close to $\alpha/2=2.5\%$, while when
$m=5000$ there is little chance of getting unstable bootstrap moments.

Tables 7–8 show the results for $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ when
$m=50,500,5000$, including the probability of covering $\mu({\bf x}_{c})$, and
the mean and SD of the interval widths. All results are based on $1000$ macro-
replications. When $m=50$, $P_{U}$ is much greater than $\alpha/2$ according
to Table 6. This explains the very large CI widths in Table 7. Nevertheless,
both $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ have reasonable coverage overall.
Notice that $\mbox{CI}_{0}$ does exhibit undercoverage when we use a very
small experiment design of $k=20$ points, while the coverage of
$\mbox{CI}_{+}$ is much closer to the nominal value of $95\%$ in this case. If
we fix the number of replications $n$ and increase the number of design points
$k$, the coverage of $\mbox{CI}_{0}$ improves. For a fixed $k$ the effect of
increasing $n$ is not as obvious.

Table 8 shows the results for $m=500,5000$. Compared with the results for
$m=50$, the mean and SD of the interval widths drop dramatically. The effects
of $k$ and $n$ are easier to discern especially when $m=5000$, which has no
unstable bootstrap moments. Specifically, for a fixed quantity of real-world
data $m$, if either the number of design points $k$ or replications per design
point $n$ is small then $\mbox{CI}_{0}$ tends to have undercoverage because it
fails to account for substantial simulation uncertainty. The most troubling
observation about $\mbox{CI}_{0}$ is that, for fixed $(n,k)$, as the amount of
input data $m$ increases its undercoverage becomes more serious. The
diminished coverage occurs because as $m\rightarrow\infty$ the width of
$\mbox{CI}_{0}$ shrinks to zero, which is not appropriate when there is still
simulation uncertainty. As the interval, $\mbox{CI}_{+}$, is able to account
for the effect of the remaining simulation estimation error, it can work under
more general situations where the simulated systems are complex and the
simulation budget is tight. As $n$ and $k$ increase, the coverages of
$\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ become closer to each other.

Table 7: Results for $\mbox{CI}_{0}$, $\mbox{CI}_{+}$ and $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=50$.  $m=50$ | $k=20$ | $k=40$
---|---|---
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 91.9% | 92.3% | 91.5% | 93.8% | 94.4% | 93.4%
Coverage of $\mbox{CI}_{+}$ | 93.9% | 94.9% | 93.7% | 94.9% | 95.6% | 95.9%
$\mbox{CI}_{0}$ Width (mean) | 326.4 | 332.4 | 339.5 | 319.1 | 328.6 | 326.5
$\mbox{CI}_{+}$ Width (mean) | 344.1 | 348.8 | 357.1 | 332.3 | 342.3 | 341.2
$\mbox{CI}_{0}$ Width (SD) | 183.1 | 173.6 | 180.7 | 176.4 | 167.6 | 175
$\mbox{CI}_{+}$ Width (SD) | 188 | 175.7 | 183.8 | 178.2 | 169.2 | 176.1
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.963 | 0.965 | 0.964 | 0.973 | 0.973 | 0.971
$m=50$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 94.6% | 96.3% | 95.4% | 94.2% | 95.1% | 95.4%
Coverage of $\mbox{CI}_{+}$ | 95.9% | 96.7% | 96.1% | 94.5% | 96% | 96.1%
$\mbox{CI}_{0}$ Width (mean) | 312.1 | 314.8 | 322.7 | 322 | 321.86 | 320
$\mbox{CI}_{+}$ Width (mean) | 322 | 325.7 | 334 | 330.2 | 331 | 329.4
$\mbox{CI}_{0}$ Width (SD) | 169.7 | 159.1 | 164.7 | 171.5 | 169.3 | 172.3
$\mbox{CI}_{+}$ Width (SD) | 171.2 | 159.4 | 165 | 172.7 | 169.5 | 172.7
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.982 | 0.98 | 0.978 | 0.985 | 0.985 | 0.983
Table 8: Results of the queueing network example for $\mbox{CI}_{0}$, $\mbox{CI}_{+}$ and $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=500$ and $m=5000$.  $m=500$ | $k=20$ | $k=40$
---|---|---
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 90.5% | 94.6% | 95.1% | 94.9% | 96.7% | 96.4%
Coverage of $\mbox{CI}_{+}$ | 95.7% | 97.7% | 97.8% | 96.6% | 98.3% | 97.8%
$\mbox{CI}_{0}$ Width (mean) | 24.8 | 28.1 | 29.4 | 27.1 | 28.5 | 28.7
$\mbox{CI}_{+}$ Width (mean) | 28.9 | 30.8 | 32.2 | 29.6 | 30.3 | 30.5
$\mbox{CI}_{0}$ Width (SD) | 19.9 | 19.4 | 20.6 | 19.1 | 19.2 | 19.9
$\mbox{CI}_{+}$ Width (SD) | 20.6 | 20.4 | 21.7 | 19.7 | 19.9 | 20.6
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.88 | 0.932 | 0.933 | 0.932 | 0.957 | 0.958
$m=500$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 96.5% | 97.5% | 95.8% | 95.4% | 96.5% | 95.9%
Coverage of $\mbox{CI}_{+}$ | 98% | 98.3% | 97.3% | 97.5% | 97.1% | 96.9%
$\mbox{CI}_{0}$ Width (mean) | 26.3 | 28 | 28.7 | 26.4 | 27.9 | 27.6
$\mbox{CI}_{+}$ Width (mean) | 28 | 29 | 29.7 | 27.9 | 28.6 | 28.2
$\mbox{CI}_{0}$ Width (SD) | 17.4 | 18 | 19.3 | 18.8 | 19.6 | 19.3
$\mbox{CI}_{+}$ Width (SD) | 17.7 | 18.4 | 19.6 | 18.9 | 19.9 | 19.5
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.952 | 0.977 | 0.978 | 0.957 | 0.984 | 0.987
$m=5000$ | $k=20$ | $k=40$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 70.7% | 89.2% | 93.1% | 81.5% | 94.3% | 94.8%
Coverage of $\mbox{CI}_{+}$ | 91.3% | 96.3% | 95.6% | 96.5% | 96.1% | 96.3%
$\mbox{CI}_{0}$ Width (mean) | 3.29 | 3.97 | 4.14 | 3.93 | 4.23 | 4.3
$\mbox{CI}_{+}$ Width (mean) | 5.85 | 4.8 | 4.56 | 6.08 | 4.64 | 4.52
$\mbox{CI}_{0}$ Width (SD) | 1.89 | 1.2 | 1 | 1.64 | 0.87 | 0.83
$\mbox{CI}_{+}$ Width (SD) | 2.12 | 1.13 | 1 | 1.52 | 0.89 | 0.85
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.588 | 0.85 | 0.924 | 0.664 | 0.924 | 0.959
$m=5000$ | $k=80$ | $k=130$
| $n=10$ | $n=50$ | $n=100$ | $n=10$ | $n=50$ | $n=100$
Coverage of $\mbox{CI}_{0}$ | 88.9% | 93.6% | 94.9% | 89.5% | 93.7% | 94.8%
Coverage of $\mbox{CI}_{+}$ | 98.1% | 95% | 96% | 98% | 95.6% | 95.5%
$\mbox{CI}_{0}$ Width (mean) | 4.54 | 4.29 | 4.29 | 4.52 | 4.35 | 4.32
$\mbox{CI}_{+}$ Width (mean) | 6.1 | 4.56 | 4.42 | 5.98 | 4.64 | 4.45
$\mbox{CI}_{0}$ Width (SD) | 1.37 | 0.85 | 0.77 | 1.28 | 0.9 | 0.79
$\mbox{CI}_{+}$ Width (SD) | 1.27 | 0.85 | 0.78 | 1.13 | 0.87 | 0.77
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ | 0.757 | 0.946 | 0.974 | 0.766 | 0.945 | 0.974

###  E.2 Performance of $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$

Tables 7–8 also demonstrate that $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$
provides a good measure of the relative contribution of model uncertainty to
overall uncertainty. For a fixed amount of real-world data $m$, increasing the
number of design points and replications $(n,k)$ drives
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ toward 1, indicating a decrease in
simulation uncertainty. For fixed simulation effort $(n,k)$, increasing the
amount of real-world data $m$ decreases
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$, indicating that there is
relatively less model uncertainty. Notice, however, that the relationship is
not simple because as $m$ increases the design space over which we fit the
metamodel becomes smaller, so that even with the same simulation effort the
absolute level of simulation uncertainty will decrease somewhat. When
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ is near $1$, the behaviors
(coverage and width) of $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ are similar and
both have coverage close to the nominal level; this is illustrated in Figure
3. Recall that $\mbox{CI}_{0}$ does not account for simulation uncertainty,
and that $\widehat{\sigma}_{I}/\widehat{\sigma}_{T}\approx 1$ indicates that
model uncertainty is large relative to simulation uncertainty, which is when
$\mbox{CI}_{0}$ will do best. Figure 3 also illustrates the general robustness
of $\mbox{CI}_{+}$.

Figure 3: The coverage errors for $\mbox{CI}_{0}$ and $\mbox{CI}_{+}$ vs.
$\widehat{\sigma}_{I}/\widehat{\sigma}_{T}$ when $m=5000$ across all values of
$n$ and $k$.