Estimation of Expectations and Variance Components in Two-Level Nested Simulation Experiments

1. Introduction and Notation

Simulation is widely recognized as an effective technique to produce forecasts, evaluate risk (see, e.g., [1]), animate and illustrate the performance of a system over time (see, e.g., [2]). When there is uncertainty in a component of a simulation model, it is said to be a random component, and it is modeled using a probability distribution and/or a stochastic process that is generated during the simulation run, to produce a stochastic simulation. Random component typically depends on the value of certain parameters, and we will denote by $\theta$ a particular value for the vector of parameters of the random components of a stochastic simulation, and $\mathsf{\Theta }$ will denote the random vector that corresponds to the parameter values when there is uncertainty on the value of these parameters.

In general, the output of a stochastic (dynamic) simulation can be regarded as a stochastic process $\left\{Y\right(s):s\ge 0;\mathsf{\Theta }\}$, where $Y\left(s\right)$ is a random vector (of arbitrary dimension d) representing the state of the simulation at time $s\ge 0$. The term transient simulation applies to a dynamic simulation that has a well-defined termination time, so that the output of a transient simulation can be viewed as a stochastic process $\left\{Y\right(s):0\le s\le T;\mathsf{\Theta }\}$, where T is a stopping time (may be deterministic), see, e.g., [3] for a definition of stopping time. Note that this notation includes the case of a discrete-time output $Z_0,Z_1,\cdots$, if we assume that $Y\left(s\right)=Z_{\lfloor s\rfloor }$ , where $\lfloor s\rfloor$ denotes the integer part of s.

A performance variable W in transient simulation is a real-valued random variable (r.v.) that depends on the simulation output up to time T, i.e., $W=f\left(Y\right(s),0\le s\ge T;\mathsf{\Theta })$, and the expectation of a performance variable W is a performance measure that we usually estimate through experimentation with the simulation model. When there is no uncertainty in the parameters of the random components, the standard methodology that is used to estimate a performance measure in transient simulation is the method of independent replications, that consists on running the simulation model to produce n replications $W_1,W_2,\cdots ,W_n$ that can be regarded as independent and identically distributed (i.i.d.) random variables (see Figure 1) .

Under the method of independent replications, a point estimator for the expectation $\alpha =E\left[W_1\right]$ is the average $\widehat{\alpha }\left(n\right)=\frac{{\sum }_{i=1}^n{W}_i}{n}$. If $E\left[\right|W_1\left|\right]<\infty$, it follows from the classical Law of Large Numbers (LLN), that $\widehat{\alpha }\left(n\right)$ is consistent, i.e., it satisfies $\widehat{\alpha }\left(n\right)\Rightarrow \alpha$ , as $n\to \infty$ (where ⇒ denotes weak convergence of random variables), see, e.g., [3] for a proof. Consistency guarantees that the estimator approaches the parameter as the number of replications n increases, and the accuracy of the simulation-based estimator $\widehat{\alpha }\left(n\right)$ is typically assessed by an asymptotic confidence interval (ACI) for the parameter. The expression for an ACI for a parameter of a stochastic simulation is usually obtained through a Central Limit Theorem (CLT) for the estimator (see, for example, chapter 3 of [4]). For the case of the expectation $\alpha$ in the algorithm of Figure 1, if $E\left[W_1^2\right]<\infty$, the classical CLT implies that

$$\frac{\sqrt{n}(\widehat{\alpha }\left(n\right)-\alpha )}{\sigma }\Rightarrow N(0,1),$$

(1)

as $n\to \infty$, where ${\sigma }^2=E\left[{(W_1-\alpha )}^2\right]$ and $N(0,1)$ denotes a r.v. distributed as normal with mean 0 and variance 1. Then, if $E\left[W_1^2\right]<\infty$, it follows from (1) and Slutsky’s Theorem (see the Appendix) that

$$\frac{\sqrt{n}(\widehat{\alpha }\left(n\right)-\alpha )}{\widehat{\sigma }\left(n\right)}\Rightarrow N(0,1),$$

as $n\to \infty$, where $\widehat{\sigma }\left(n\right)$ denotes the sample standard deviation, i.e., ${\widehat{\sigma }}^2\left(n\right)=\frac{{\sum }_{i=1}^n{(W_i-\widehat{\alpha }\left(n\right))}^2}{n-1}$. This CLT implies that

$$\underset{n\to \infty }{lim}P\left[\right|\widehat{\alpha }\left(n\right)-\alpha |\le z_{\beta }\widehat{\sigma }\left(n\right)/\sqrt{n}]=1-\beta ,$$

for $0<\beta <1$, where $z_{\beta }$ denotes the ($1-\beta /2$)-quantile of a N(0,1), which is sufficient to establish a $(1-\beta )100\%$ ACI for $\alpha$ with halfwidth

$$H{W}_{\alpha }=z_{\beta }\widehat{\sigma }\left(n\right)/\sqrt{n}.$$

(2)

A halfwidth in the form of (2) is the typical measure used in simulation software (e.g., Simio, see [2]) to assess the accuracy of $\widehat{\alpha }\left(n\right)$ for the estimation of expectation $\alpha$.

In contrast to the estimation of (output) performance measures, parameters of (input) random components of a simulation model are usually estimated from real-data observations (x) and, while most applications covered in the relevant literature assume that no uncertainty exists in the value of these parameters, the uncertainty can be significant when little data is available. In these cases, Bayesian statistics can be used to incorporate this uncertainty in the output analysis of simulation experiments via the use of a posterior distribution $p\left(\theta \right|x)$. A methodology currently proposed for the analysis of simulation experiments under parameter uncertainty, is a two-level nested simulation algorithm (see, e.g., [6,7,8]. In the outer level, we simulate (n) observations for the parameters from a posterior distribution $p\left(\theta \right|x)$, while in the inner level we simulate (m) observations for the performance variable with the parameters fixed at the value ($\theta$) generated in the outer level (see Figure 2). In this paper, we focus on the output analysis of two-level simulation experiments, for the case where the observations at the inner level are independent, showing how the variance of a simulated observation can be decomposed into parametric and stochastic variance components. Afterwards, we derive a CLT for both the estimator of the point forecast and the estimators of the variance components. Our CLTs allow us to compute an ACI for each estimator. Our results are validated through experiments with a forecast model for sporadic demand reported in [10]. This paper is an extended version of results initially reported in [11] and the missing proofs in [11] are provided.

Following this introduction, we present the proposed methodology for the construction of an ACI for the point forecast and the variance components in a two-level simulation experiment. Afterwards, we present an illustrative example that has an analytical solution for the parameters of interest in this paper. This example is used in the next section to illustrate the application and validity of our proposed methodologies for the construction of an ACI. Finally, in the last section, we present conclusions and directions for future research.

2. Theoretical Results

To identify the variance components in each observation $W_{ij}$ of the algorithm illustrated in Figure 2, let $\mu \left(\mathsf{\Theta }\right)=E\left[W_{11}\right|\mathsf{\Theta }]$, and ${\sigma }^2\left(\mathsf{\Theta }\right)=E\left[W_{11}^2\right|\mathsf{\Theta }]-{\mu }^2\left(\mathsf{\Theta }\right)$. Under this notation, the point forecast is $\alpha =E\left[\mu \right(\mathsf{\Theta }\left)\right]$, and the variance of each $W_{ij}$ is:

$$V\left[W_{ij}\right]\stackrel{def}{=}E\left[W_{ij}^2\right]-E{\left[W_{ij}\right]}^2=E[E\left[W_{ij}^2\right|\mathsf{\Theta }]-\mu {\left(\mathsf{\Theta }\right)}^2]+E\left[\mu {\left(\mathsf{\Theta }\right)}^2\right]-E{\left[\mu \left(\mathsf{\Theta }\right)\right]}^2={\sigma }_S^2+{\sigma }_P^2,$$

(3)

for $i=1,...,n$; $j=1,...,m$, where ${\sigma }_P^2=V\left[\mu \left(\mathsf{\Theta }\right)\right]\stackrel{def}{=}E\left[\mu {\left(\mathsf{\Theta }\right)}^2\right]-E{\left[\mu \left(\mathsf{\Theta }\right)\right]}^2$, and ${\sigma }_S^2=E\left[{\sigma }^2\left(\mathsf{\Theta }\right)\right]$. It is worth mentioning that, in the relevant literature, ${\sigma }_S^2$ is commonly referred to as stochastic variance and ${\sigma }_P^2$ is commonly referred to as parametric variance.

2.1. Point Estimators

In this paper, we are interested in both the estimation of the point forecast $\alpha =E\left[\mu \right(\mathsf{\Theta }\left)\right]$ and the estimators of the variance components of every observations generated in the algorithm of Figure 2 and defined in (3), thus we first consider the natural point estimators

$$\widehat{\alpha }\left(n\right)=\frac{1}{n}\sum _{i=1}^n{\widehat{\alpha }}_i,{\widehat{\sigma }}_T^2\left(n\right)=\frac{1}{n-1}\sum _{i=1}^n{({\widehat{\alpha }}_i-\widehat{\alpha }\left(n\right))}^2,{\widehat{\sigma }}_S^2\left(n\right)=\frac{1}{n}\sum _{i=1}^n{S}_i^2,$$

(4)

where ${\widehat{\alpha }}_i=m^{-1}{\sum }_{j=1}^m{W}_{ij}$, and $S_i^2={(m-1)}^{-1}{\sum }_{j=1}^m{(W_{ij}-{\widehat{\alpha }}_i)}^2$, $i=1,...,m$. Note that the ${\widehat{\alpha }}_i$’s are i.i.d. with expectation $E\left[{\widehat{\alpha }}_1\right]=\alpha$ and variance

$$\begin{array}{ccc}\hfill {\sigma }_T^2& \stackrel{def}{=}& E\left[{({\widehat{\alpha }}_1-\alpha )}^2\right]=m^{-2}(mE\left[{(W_{11}-\alpha )}^2\right]+m(m-1)E\left[(W_{11}-\alpha )(W_{12}-\alpha )\right])\hfill \\ & =& m^{-1}({\sigma }_S^2+{\sigma }_P^2)+m^{-1}(m-1){\sigma }_S^2={\sigma }_S^2+m^{-1}{\sigma }_P^2.\hfill \end{array}$$

(5)

On the other hand, the $S_i^2$ are i.i.d. with expectation $E\left[S_1^2\right]={\sigma }_S^2$. Thus, the next proposition follows from the classical LLN .

Proposition 1.

Given $m\ge 1$, if $E\left[W_{11}^2\right]<\infty$ then $\widehat{\alpha }\left(n\right)$ and ${\widehat{\sigma }}_T^2\left(n\right)$ are unbiased and consistent (as $n\to \infty$) estimators for α and ${\sigma }_T^2$ (as defined in (5)), respectively. Furthermore, if $m\ge 2$ and $E\left[W_{11}^2\right]<\infty$, then ${\widehat{\sigma }}_S^2\left(n\right)$ is an unbiased and consistent (as $n\to \infty$) estimator for ${\sigma }_S^2$ (as defined in (3)).

2.2. Accuracy of the Point Estimators

As we established in Proposition 1, under mild assumptions the point estimators proposed in (4) are consistent, and thus converge to the corresponding parameter value (as $n\to \infty$). Nonetheless, to establish the level of accuracy of these estimators, we must establish a CLT for each estimator to derive a valid expression for the corresponding ACI. Note that both $\widehat{\alpha }\left(n\right)$ and ${\widehat{\sigma }}_S^2\left(n\right)$ are averages of i.i.d observations, thus the next proposition follows from the classical CLT for i.i.d. observations.

Proposition 2.

Given $m\ge 1$, if $E\left[W_{11}^2\right]<\infty$ then

$$\frac{\sqrt{n}(\widehat{\alpha }\left(n\right)-\alpha )}{{\sigma }_T}\Rightarrow N(0,1),$$

as $n\to \infty$. Furthermore, if $m\ge 2$ and $E\left[W_{11}^4\right]<\infty$, then

$$\frac{\sqrt{n}({\widehat{\sigma }}_S^2\left(n\right)-{\sigma }_S^2)}{\sqrt{V_S}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\sigma }_S^2$ is defined in (3), ${\sigma }_T^2$ is defined in (5), $\widehat{\alpha }\left(n\right)$, ${\widehat{\sigma }}_S^2\left(n\right)$, $S_1^2$ are defined in (4), and $V_S=E\left[{(S_1^2-{\sigma }_S^2)}^2\right]$.

Since we have consistent estimators for ${\sigma }_T^2$ and $V_S$ (under mild assumptions), the next corollary follows from Proposition 1 and Slutsky’s Theorem, details of a proof are given in the Appendix.

Corollary 1.

Under the same notation and assumptions as in Proposition 2, for $m\ge 1$ we have

$$\frac{\sqrt{n}(\widehat{\alpha }\left(n\right)-\alpha )}{\sqrt{{\widehat{\sigma }}_T^2\left(n\right)}}\Rightarrow N(0,1),$$

as $n\to \infty$, and for $m\ge 2$ we have

$$\frac{\sqrt{n}({\widehat{\sigma }}_S^2\left(n\right)-{\sigma }_S^2)}{\sqrt{{\widehat{V}}_S\left(n\right)}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\widehat{\sigma }}_S^2\left(n\right)$ and ${\widehat{\sigma }}_T^2\left(n\right)$ are defined in (4), and

$$\widehat{V_s}\left(n\right)=\frac{1}{n-1}\sum _{i=1}^n{(S_i^2-S^2)}^2,S^2=\frac{1}{n}\sum _{i=1}^n{S}_i^2$$

In order to obtain a CLT for ${\widehat{\sigma }}_T^2\left(n\right)$, note that this estimator is the sample variance of a set of i.i.d. observations, thus we can use the following Lemma. A proof using the Delta Method (see, e.g., Proposition 2 of [9] for a proof) is provided in the Appendix.

Lemma 1.

If $X_1,X_2,...$ is a sequence of i.i.d. random variables with $E\left[X_1^4\right]<\infty$, then

$$\frac{\sqrt{n}(S^2\left(n\right)-{\sigma }_1^2)}{\sqrt{{\sigma }_2^2}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\sigma }_1^2={\mu }_2-{\mu }_1^2$, ${\sigma }_2^2={\mu }_1^2{\mu }_2-4{\mu }_1^4-4{\mu }_1{\mu }_3+{\mu }_4-{\mu }_2^2$, ${\mu }_k=E\left[X_1^k\right]$, $k=1,2,3,4$; $S^2\left(n\right)={(n-1)}^{-1}{\sum }_{i=1}^n{(X_i-{\widehat{\mu }}_1)}^2$, ${\widehat{\mu }}_1=n^{-1}{\sum }_{i=1}^n{X}_i$.

Corollary 2.

Under the same assumptions as in Lemma 1 we have

$$\frac{\sqrt{n}(S^2\left(n\right)-{\sigma }_S^2)}{\sqrt{{\widehat{\sigma }}_2^2\left(n\right)}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\widehat{\sigma }}_2^2=8{\widehat{\mu }}_1^2{\widehat{\mu }}_2-4{\widehat{\mu }}_1^4-4{\widehat{\mu }}_1{\widehat{\mu }}_3+{\widehat{\mu }}_4-{\widehat{\mu }}_2^2$, ${\widehat{\mu }}_k=n^{-1}{\sum }_{i=1}^n{X}_i^k$.

Corollary 2 follows from the fact that ${\widehat{\mu }}_k$ is an unbiased and consistent estimator of ${\mu }_k$, and the next corollary follows from the fact that ${\widehat{\sigma }}_T^2\left(n\right)$ is the sample variance of the ${\widehat{\alpha }}_i$.

Corollary 3.

Given $m\ge 1$, if $E\left[W_{11}^4\right]<\infty$ then

$$\frac{\sqrt{n}({\widehat{\sigma }}_T^2\left(n\right)-{\sigma }_T^2)}{\sqrt{{\widehat{V}}_T\left(n\right)}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\widehat{V}}_T\left(n\right)=8{\overline{\alpha }}_1^2{\overline{\alpha }}_2-4{\overline{\alpha }}_1^4-4{\overline{\alpha }}_1{\overline{\alpha }}_3+{\overline{\alpha }}_4-{\overline{\alpha }}_2^2$, ${\overline{\alpha }}_k=n^{-1}{\sum }_{i=1}^n{\widehat{\alpha }}_i^k$.

Let $0<\beta <1$, and using corollaries 1 and 3 we can establish a $100(1-\beta )\%$ ACI for the point forecast $\alpha$, and variance components ${\sigma }_S^2$ and ${\sigma }_T^2={\sigma }_S^2+m^{-1}{\sigma }_P^2$; each ACI is centered in the corresponding point estimator ($\widehat{\alpha }\left(n\right)$, ${\widehat{\sigma }}_S^2\left(n\right)$ or ${\widehat{\sigma }}_T^2\left(n\right)$) and the corresponding halfwidth is given by:

$$H{W}_{\alpha }=z_{\beta }\frac{\sqrt{{\widehat{\sigma }}_T^2\left(n\right)}}{\sqrt{n}},H{W}_{{\sigma }_S^2}=z_{\beta }\frac{\sqrt{{\widehat{V}}_S\left(n\right)}}{\sqrt{n}},\mathit{and}H{W}_{{\sigma }_T^2}=z_{\beta }\frac{\sqrt{{\widehat{V}}_T\left(n\right)}}{\sqrt{n}},$$

(6)

for $\alpha$, ${\sigma }_2^2$ and ${\sigma }_T^2$, respectively, where ${\widehat{\sigma }}_T^2\left(n\right)$ is defined in (4), ${\widehat{V}}_S\left(n\right)$ and ${\widehat{V}}_T\left(n\right)$ are defined in Corollary 1, and in Corollary 3, respectively.

Note that the ACIs proposed in (6) assume that the value of m in the algorithm of Figure 2 is fixed and the accuracy of the estimator improves as n (the number of observations in the outer level) increases (in turn, the halfwidth of the ACI gets smaller). Given that we can build a valid ACI for any value of m, a relevant question is how to find an adequate value of m to get an acceptable level of accuracy in a reasonable amount of running time. In order to answer this question for the case of the point estimator of $\alpha$, let us fix the total number of iterations in the algorithm of Figure 2 to k = $nm$, and note from (5) and Proposition 2 that the asymptotic variance of $(\widehat{\alpha }\left(n\right)-\alpha )$ is

$$n^{-1}{\sigma }_T^2=k^{-1}(m{\sigma }_S^2+{\sigma }_P^2),$$

(7)

and takes its minimal value for $m=1$, suggesting that the point estimator $\widehat{\alpha }\left(n\right)$ defined in (4) is more accurate as m approaches the value of 1. Note that for $m=1$, a fixed number of iterations $k=nm$ is convenient (from the point of view of running time), when the computation of $W_{ij}$ requires the same or more computation time as ${\mathsf{\Theta }}_i$, as suggested in the relevant literature (see, for example, [6]). Furthermore, if we allow m to increase with n, we can obtain the following proposition (a proof using Lindeberg-Feller Theorem is provided in the Appendix).

Proposition 3.

Given $0<p\le 1$, if $m=*n^{-1+1/p}$ and $E\left[W_{11}^2\right]<\infty$ then

$$\frac{\sqrt{n}(\widehat{\alpha }\left(n\right)-\alpha )}{\sqrt{{\sigma }_T^2}}\Rightarrow N(0,1),$$

as $n\to \infty$, where ${\sigma }_T^2$ is defined in (5).

Note that the last proposition implies that the ACI defined in equation (6) for the point forecast $\alpha$ is also valid under the assumptions of Proposition 3. If, once again, we set the total number of iterations in the algorithm of Figure 2 to $k=nm$, we let $n\approx k^p$, $m\approx k^{1-p}$, and nm = k, it follows from Proposition 3 that the asymptotic variance of $\widehat{\alpha }\left(n\right)$ is $n^{-1}{\sigma }_T^2\approx k^{-p}(k^{1-p}{\sigma }_S^2+{\sigma }_P^2)$ for every $0\le p\le 1$. Note that, for fixed k, $n^{-1}{\sigma }_T^2$ reaches its minimum value when $p=1$, that is, when $n=k$ and $m=1$. However, note that we need $m\ge 2$ in order to estimate ${\sigma }_S^2$. In the following section we report some empirical results that confirm our theoretical results. It is worth mentioning that the case $n=k$ and $m=1$ has been reported in the literature as the posterior sampling algorithm (see, e.g., [12,13])

3. An Example with Analytical Solution

The following model (reported in [10]) has been proposed to forecast sporadic demand by incorporating data on times between arrivals and customer demand; where uncertainty on the model parameters is incorporated using a Bayesian approach. For this model, we will show analytical expressions for the performance measures defined in Section 2. These expressions are used in the following section to illustrate the validity of the ACIs proposed in the previous section.

Customer arrivals for a particular item in a shop follow a Poisson process, yet there is uncertainty in the arrival rate ${\mathsf{\Theta }}_0$, so that given $[{\mathsf{\Theta }}_0={\theta }_0]$, interarrival times between customers are i.i.d. with exponential density:

$$f\left(y\right|{\theta }_0)=\left\{\begin{array}{cc}\hfill {\theta }_0{e}^{-{\theta }_{0}y},& \hfill y>0,\\ \hfill 0,& \hfill \mathrm{otherwise},\end{array}\right.$$

(8)

where ${\theta }_0\in S_{00}=(0,\infty )$. Every client can order j units of this item with probability ${\mathsf{\Theta }}_{1j}$, $j=1,...,q$, $q\ge 2$. Let ${\mathsf{\Theta }}_1=({\mathsf{\Theta }}_{11},...{\mathsf{\Theta }}_{1(q-1)})$ and ${\mathsf{\Theta }}_{1q}=1-{\sum }_{j=1}^{q-1}{\mathsf{\Theta }}_{1j}$, then $\mathsf{\Theta }=({\mathsf{\Theta }}_0,{\mathsf{\Theta }}_1)$ is the parameter vector, and $S_0=S_{00}⨂S_{01}$ is the parameter space, where $S_{01}=\{({\theta }_{11},...,{\theta }_{1(q-1)}):{\sum }_{j=1}^{q-1}{\theta }_{1j}\le 1;{\theta }_{1j}\ge 0,j=1,...,q-1\}$.

Total demand during a period of length T is

$$D=\left\{\begin{array}{cc}\hfill \sum _{i=1}^{N\left(T\right)}{U}_i,& \hfill N\left(T\right)>0\\ \hfill 0,& \hfill \mathrm{otherwise},\end{array}\right.$$

(9)

where $N\left(s\right)$ is the number of customer arrivals during the interval $[0,s]$, $s\ge 0$, and $U_1,U_2,...$ are the individual demands (conditionally independent relative to $\mathsf{\Theta }$). The information about $\mathsf{\Theta }$ consists of i.i.d. observations $v=(v_1,...,v_r)$, $u=(u_1,...,u_r)$ of past customers, where $v_i$ is the interarrival time between customer i and customer ($i-1$), and $u_i$ is the number of units ordered by client i. By taking Jeffrey’s non-informative prior as the prior density for $\mathsf{\Theta }$, we obtain the posterior density (see [10] for details) $p\left(\theta \right|x)=p\left({\theta }_0\right|v)p\left({\theta }_1\right|u)$, where $x_i=(v_i,u_i)$, $i=1,...,r$ , $x=(x_1,...,x_r)$, $\theta =({\theta }_0,{\theta }_1)$ as

$$p\left({\theta }_0\right|v)=\frac{{\theta }_0^{r-1}{\left({\sum }_{i=1}^r{v}_i\right)}^r{e}^{-{\theta }_0{\sum }_{i=1}^r{v}_i}}{(r-1)!},p\left({\theta }_1\right|u)=\frac{{(1-{\sum }_{j=1}^{q-1}{\theta }_{1j})}^{c_q-1/2}{\mathsf{\Pi }}_{j=1}^{q-1}{\theta }_{1j}^{c_j-1/2}}{B(c_1+1/2,...,c_q+1/2)},$$

(10)

where $c_j={\sum }_{i=1}^{r}I[u_i=j]$, and $B(a_1,...a_q)={\mathsf{\Pi }}_{j=1}^q\Gamma \left(a_j\right)/\Gamma \left({\sum }_{j=1}^q{a}_j\right)$, for $a_1,...,a_q>0$. Using this notation, we can show that (see [1] for details)

$$\alpha =E\left[T{\mathsf{\Theta }}_0\right]\sum _{j=1}^{q}j{p}_j,$$

$${\sigma }_P^2=\frac{E\left[T^2{\mathsf{\Theta }}_0^2\right]}{(q_0+1)}\sum _{j=1}^q{j}^2{p}_j+\frac{E{\left[T{\mathsf{\Theta }}_0\right]}^2[(q_0/n)-1]}{(q_0+1)}{\left(\sum _{j=1}^{q}j{p}_j\right)}^2,$$

$${\sigma }_S^2=E\left[T{\mathsf{\Theta }}_0\right]\sum _{j=1}^q{j}^2{p}_j,$$

where $E\left[T{\mathsf{\Theta }}_0\right]=Tr{\left({\sum }_{i=1}^r{v}_i\right)}^{-1}$, $E\left[T^2{\mathsf{\Theta }}_o^2\right]=T^{2}r(1+r){\left({\sum }_{i=1}^r{v}_i\right)}^{-2}$, $p_j=q_j/q_0$, $q_j=c_j+1/2$, $j=1,...,q$, $q_0={\sum }_{j=1}^q{q}_j$, and $c_j$ are defined in (10).

4. Empirical Results

To validate the ACIs proposed in (4), we conducted some experiments with the Bayesian model of the previous section to illustrate the estimation of $\alpha$, ${\sigma }_S^2$ and ${\sigma }_T^2$. We considered the values $T=15$, $r=20$, ${\sum }_{i=1}^r{x}_i=10$, $q=5$, $c_1=5$, $c_2=3$, $c_3=2$, $c_4=3$, $c_5=7$. With this data, the point forecast is $\alpha \approx 95.333$, and the variance components are ${\sigma }_S^2\approx 380.667$, ${\sigma }_P^2\approx 568.598$. The empirical results that we report below illustrate a typical behavior that we should experiment for any other feasible data set.

In all the experiments reported in this Section we considered 1000 independent replications of the algorithm of Figure 2 for different number of observations in the outer level (n) and in the inner level (m); in each replication we computed the point estimators for $\alpha$, ${\sigma }_S^2$, and ${\sigma }_T^2$, and the corresponding halfwidths of 90% ACI’s according to equations (6). Since we know the value of the parameters we are estimating, we were able to report (for n and m given), the empirical coverage (i.e., the fraction of independent replications that the corresponding ACI covered the true parameter value), the average and standard deviation of halfwidths, and the squared root of the empirical mean squared error defined by

$$RMSE=\sqrt{\frac{1}{n_0}{\displaystyle \sum _{i=1}^{n_0}}{({\widehat{\theta }}_i-\theta )}^2},$$

where ${\widehat{\theta }}_i$ denotes the value obtained in the i-th replication for the estimation of a parameter $\theta ,i=1,2,\cdots ,n_0$ ($n_0=1000$ in our experiments).

In a first set of experiments we considered $nm=200,2000,20000$, and $m=2,4,8$ for each value of $nm$, to compare the effect of increasing the number of observations in the inner level for a given value value of $nm$. The results of this set of experiments are summarized in Figure 3, Figure 4 and Figure 5. Note that we are not considering $m=1$ in this set of experiments to be able to construct an ACI for the stochastic variance ${\sigma }_S^2$.

In Figure 3 we illustrate the performance measures for the quality of the estimation procedure that we obtained for the estimation of the point forecast $\alpha$. As we observe from Figure 3, the coverages are acceptable (very close to the nominal value of 0.9, even for $n=100$). These results validate the ACI defined in (6) for the point forecast $\alpha$. We also observe from Figure 3 that the RMSE, average halfwidth and standard deviation of halwidths improve (decrease) as the number of observations in the outer level (n) increases, as suggested by Corollary 1. Note also from Figure 3 that a smaller value of m provides smaller RMSE, average halfwidths and standard deviations of halwidths, validating our theoretical results.

In Figure 4 we illustrate the performance measures for the quality of the estimation procedure that we obtained for the estimation of the stochastic variance ${\sigma }_S^2$. As we observe from Figure 4, the coverages are acceptable (very close to the nominal value of 0.9, even for $n=100$). These results validate the ACI defined in (6) for the stochastic variance ${\sigma }_S^2$. We also observe from Figure 3 that the RMSE, average halfwidth and standard deviation of halwidths improve (decrease) as the number of observations in the outer level (n) increases, as suggested by Corollary 2. However, contrary to what we observe for the estimation of $\alpha$, a larger value of m provides smaller RMSE, average halfwidths and standard deviations of halwidths, suggesting that, for a fixed value of $nm$, the quality of the estimation for the stochastic variance ${\sigma }_S^2$ improves as the number of the observations in the inner loop (m) increases.

For the estimation of the total variance ${\sigma }_T^2$ (illustrated n Figure 5) we obtained similar results for the quality of the estimation as for the estimation of the point forecast $\alpha$, except that a larger values of n is required to obtain reliable coverages. As we observe from Figure 3, the coverages are acceptable (very close to the nominal value of 0.9, for $n=1000$ and 10000). These results validate the ACI defined in (6) for the total variance ${\sigma }_T^2$. We also observe from Figure 3 that the RMSE, average halfwidth and standard deviation of halwidths improve (decrease) as the number of observations in the outer level (n) increases, as suggested by Corollary 3. Note also from Figure 5 that a smaller value of m provides smaller RMSE, average halfwidths and standard deviations of halwidths, validating our theoretical results.

Figure 4. Performance of the estimation of stochastic variance ${\sigma }_S^2$ for $nm$ fixed comparing different values of m.

Figure 4. Performance of the estimation of stochastic variance ${\sigma }_S^2$ for $nm$ fixed comparing different values of m.

Figure 5. Performance of the estimation of total variance ${\sigma }_T^2$ for $nm$ fixed comparing different values of m.

Figure 5. Performance of the estimation of total variance ${\sigma }_T^2$ for $nm$ fixed comparing different values of m.

In a second set of experiments we considered $nm=100,1000,10000$, with $m=1$ and $m\approx {\left(nm\right)}^{1/3}$ for each value of $nm$, to compare the quality of the estimation procedures using the value of m that we suggest as optimal for the estimation of point forecast $\alpha$ with the value of m suggested in [6] as an adequate choice for m in the case of biased estimators in the inner level of the algorithm of Figure 2. The results of this set of experiments are summarized in Figure 6 and Figure 7. Note that we are not considering the estimation of the stochastic variance ${\sigma }_S^2$ in this set of experiments because $m\ge 2$ is required to construct an ACI for the stochastic variance ${\sigma }_S^2$. Note also that we considered ${100}^{1/3}\approx 5$, ${1000}^{1/3}\approx 10$, and ${10000}^{1/3}\approx 20$, and we are using the same color for $m=5,10,20$ in Figure 6 and Figure 7.

In Figure 6 we illustrate the performance measures for the quality of the estimation procedure that we obtained for the estimation of the point forecast $\alpha$ in our second set of experiments . As we observe from Figure 6, the coverages are acceptable (very close to the nominal value of 0.9, even for $n=100$). These results validate the ACI defined in (6) for the point forecast $\alpha$, and the ACI suggested by Proposition 3. We also observe from Figure 6 that the RMSE, average halfwidth and standard deviation of halwidths are worse for $m\approx {\left(nm\right)}^{1/3}$, confirming our finding that, for the same number of replications $nm$, $m=1$ produces better point estimators for $\alpha$ than $m\approx {\left(nm\right)}^{1/3}$ confirming the result of Proposition 3.

Finally, in Figure 7 we show the results of our second set of experiments for the estimation of the total variance ${\sigma }_T^2$. We found similar resulta as for the case of the estimation of the point forecast $\alpha$, coverages are very good (even for n = 100), and all performance measure for the ACI (RMSE, average and standard deviation of halfwidths) are worse for $m\approx {\left(nm\right)}^{1/3}$, suggesting that, for the same number of replications $nm$, $m=1$ produces better point estimators for ${\sigma }_T^2$ than $m\approx {\left(nm\right)}^{1/3}$.

5. Conclusions

In this paper, we propose methodologies to calculate point estimators (and their corresponding halfwidths), for both the point forecast and the variance components in two-level nested stochastic simulation experiments, for the case where the observations at both levels of the algorithm are independent. These methods can be applied to the construction of Bayesian forecasts based on experiments using a simulation model under parameter uncertainty.

Both our theoretical and our experimental results confirm that the proposed point estimators and their corresponding halfwidths are asymptotically valid, i.e., the point estimators converge to the corresponding parameter values and the halfwidths converge to the nominal coverage as the number of replications (n) of the outer level increases.

Furthermore, given a fixed number of total observations ($nm$), we show that the choice of only one replication in the inner level ($m=1$) provides more accurate estimators for both the point forecast ($\alpha$), and the variance of the point forecast (${\sigma }_T^2$). However, $m\ge 2$ is required for the estimation of ${\sigma }_S^2$.

Directions for future research on this topic includes experimentation with other point estimators, such as, quasi Monte Carlo or Simpson integration, with the objective of finding more accurate point estimators for the parameters considered in this paper.