6. Output Analysis

How to draw statistically valid conclusions from simulation output: transient versus steady-state behavior, confidence intervals for terminating simulations, and initial-data deletion, replication/deletion, batch means, regenerative and sequential methods for steady-state parameters.

flowchart LR
    ID(["input distributions"]) --> RV["random variates (rv's)"] --> P["simulation program"] --> Y["output random variables Y₁, Y₂, …"]

The inputs are random, so the outputs $Y_1, Y_2, \dots$ are random variables too. A performance measure is a statistic, i.e. a function of the output r.v.s:

$$h(Y_1, Y_2, \dots) \;\to\; \text{performance measure}, \qquad \text{e.g. } h = \frac{Y_1 + \cdots + Y_n}{n} .$$

Wrong procedure in practice: perform a single simulation run of arbitrary length and treat the simulation estimates as the true characteristics of the system. Why it is wrong:

• One run gives one realization of $h(Y_1, Y_2, \dots)$, which is a random variable: a different seed gives a different answer, and a single number says nothing about how far it may be from the truth.
• Within a run the $Y_i$’s are not i.i.d. (see the example below), so the formulas of §2.6, which need i.i.d. data, cannot be applied to them.
• With an arbitrary run length, the result can depend heavily on the initial conditions (§6.1).

Q: What is the correct way to estimate the true performance of the system based on simulation outputs?

Ex: queueing simulation. $Y_i$ = the delay of the $i$-th customer. Lindley’s equation (§2.5):

$$Y_1 = 0, \qquad Y_{i+1} = \max\{0,\ Y_i + S_i - X_{i+1}\},$$

where $S_i$ is the service time of customer $i$ and $X_{i+1}$ the interarrival time between customers $i$ and $i+1$.

The $Y_i$’s are not i.i.d.:

• Not independent: $Y_{i+1}$ is computed from $Y_i$. A customer who waits long is usually followed by another who waits long (positive correlation).
• Not identically distributed: $Y_1 = 0$ always (the system starts empty), $Y_2$ is often $0$, and the distribution of $Y_i$ keeps changing with $i$ (§6.1).

Solution: perform multiple replications, i.e. independent runs, each with the same initial conditions but its own independent random numbers (e.g. a different stream per run, §3.3).

output process	$Y_1$	$Y_2$	$\cdots$	$Y_m$
run 1	$Y_{11}$	$Y_{12}$	$\cdots$	$Y_{1m}$
run 2	$Y_{21}$	$Y_{22}$	$\cdots$	$Y_{2m}$
$\vdots$	$\vdots$	$\vdots$		$\vdots$
run $n$	$Y_{n1}$	$Y_{n2}$	$\cdots$	$Y_{nm}$

$Y_{ik}$ = the $k$-th observation in run $i$, a random observation obtained from a particular run. Along a row (within one run) the observations are dependent. Down a column (fixed $k$, different runs), $Y_{1k}, Y_{2k}, \dots, Y_{nk}$ are i.i.d., because the runs are independent copies of the same experiment. So

$$\frac{1}{n}\sum_{i=1}^n Y_{ik} \quad \text{is an unbiased estimator of } \mathbf{E}[Y_k],$$

and the usual i.i.d. statistics (CIs, §2.6) apply to each column.

Idea: use the observations in multiple runs to draw statistical inference about $Y_1, Y_2, \dots, Y_m$.

6.1. Transient and Steady-State Behavior of a Stochastic Process

Output process $Y_1, Y_2, \dots, Y_j, \dots$ with $Y_j \sim F_j(\cdot \mid I)$, where $F_j(y \mid I) = P(Y_j \le y \mid I)$ and $I$ is the initial condition (e.g. the number of customers present at time 0).

The fact that $F_j(\cdot \mid I)$ is different for different $Y_j$’s is called the transient behavior; $F_j(\cdot \mid I)$ is the transient distribution of $Y_j$. If there is a c.d.f. $F(y)$ such that

$$F_j(y \mid I) \to F(y) \quad \text{as } j \to \infty, \text{ for all } y \text{ and all initial conditions } I,$$

then $F(y)$ is called the steady-state distribution, and $Y \sim F(y)$ is the steady-state random variable.

This is convergence in distribution (weak convergence, §2.4): $Y_j \xrightarrow{d} Y$. What it means, and what it does not:

• The limit $F$ does not depend on $I$: every initial condition leads to the same steady-state distribution. But how fast $F_j(\cdot \mid I)$ approaches $F$ does depend on $I$.
• Strictly, $F$ is only reached as $j \to \infty$. In practice there is a $k$ such that $F_j \approx F$ for all $j > k$, and steady state is said to “start” around $k$.
• Steady state does not mean that the $Y_j$’s stop changing or become equal within a run. It means that they have (approximately) the same distribution; they are still dependent on each other.

E.g. M/M/1 queue with $\rho = 0.9$ (arrival rate 1, mean service time 0.9), averaged over 20,000 replications. Started empty and idle, the mean delay $\mathbf{E}[Y_j \mid I]$ rises from 0; started with a long first wait, it falls. Both approach the same steady-state mean $\mathbf{E}[Y] = 8.1$, but only after several hundred customers. (The average delay of the first 25 customers from an empty start is about 2.12: a transient quantity, far from 8.1. It reappears in §6.3.2.)

6.2. Types of Simulations

• terminating
• nonterminating
  • steady-state parameters: $Y_1, Y_2, \dots \to Y$; characteristics of $Y$, e.g. $\mathbf{E}[Y]$, $P(Y > s)$, $\operatorname{Var}(Y)$
  • steady-state cycle parameters

Terminating simulation

• There is a natural event $E$ that determines the length of each simulation run.
• Ex: simulate the operation of a retail store over an 8-hr shift; $E$ = {8 hours of simulated time have elapsed}.
• Ex: a manufacturer receives a contract to produce 100 planes; $E$ = {100 planes have been completed}.
• It never reaches steady state, and it does not need to: the run stops at $E$, and the behavior up to $E$, transient part included, is exactly what is being measured.
• Both the initial conditions and the stopping event should be representative of the real system (e.g. the store opens empty every morning, so each run starts empty).

Nonterminating simulation

• Any system that operates continuously, with no natural event $E$ to stop a run (e.g. a production line running around the clock).
• Interested in long-run / steady-state behavior.
• The initial condition does not matter: the steady-state distribution is the same for every $I$ (§6.1).
• Question: how long should the simulation be performed to ensure that the influence of the initial condition can be ignored? (In the figure of §6.1, several hundred customers.)
• When the system is cyclic, we are interested in steady-state cycle parameters.

Ex: $Y_1, Y_2, \dots, Y_8$, where $Y_i$ = server utilization in hour $i$ of an 8-hr shift at a bank (one that operates continuously, shift after shift). The hourly utilizations follow the same pattern in every shift (busy hours, quiet hours), so $Y_i$ never settles into one distribution: there is no steady state. Instead, look at whole cycles (shifts):

$$\bar Y_1^c = \frac{Y_1 + \cdots + Y_8}{8} \ \ \text{(average server utilization over shift 1)}, \qquad \bar Y_2^c = \frac{Y_9 + \cdots + Y_{16}}{8}, \ \dots$$

The cycle averages $\bar Y_1^c, \bar Y_2^c, \dots$ can have a steady-state distribution, and

$$\lim_{n \to \infty}\mathbf{E}\big[\bar Y_n^c\big] \quad \text{is a steady-state cycle parameter.}$$

6.3. Statistical Analysis for Terminating Simulations

Make $n$ independent replications of a terminating simulation: each replication starts with the same initial conditions, stops at the event $E$, and uses its own independent random numbers. Let $X_j$ be the performance measure from replication $j$, i.e. a random variable computed from the outputs of run $j$, for example

$$X_j = h(Y_{j1}, \dots, Y_{jm}) = \text{average delay in queue of the customers in replication } j$$

(one $X_j$ per row of the replication table at the start of the chapter). Since the replications are independent copies of the same experiment,

$X_1, X_2, \dots, X_n$ are i.i.d.

Estimating means

Goal: to estimate $\mu = \mathbf{E}[X_j]$.

$\bar X(n) = \dfrac{1}{n}\sum_{j=1}^n X_j$ is an unbiased estimator for $\mu$.

A $100(1-\alpha)\%$ CI for $\mu$ is

$$\left[\bar X(n) - t_{\alpha/2,\,n-1}\sqrt{\frac{S^2(n)}{n}},\;\; \bar X(n) + t_{\alpha/2,\,n-1}\sqrt{\frac{S^2(n)}{n}}\,\right]$$

(Case 2 of §2.6.1). It is exact if the $X_j$’s are normal, and approximate otherwise (see §6.3.2). Here

$$S^2(n) = \frac{\sum_{j=1}^n \big(X_j - \bar X(n)\big)^2}{n-1}$$

is the sample variance of the $X_j$’s.

Ex: ten independent replications of a single-server queue model; $X_j$ = the average delay in queue in replication $j$:

$$X_j = 2.02,\ 0.73,\ 3.20,\ 6.23,\ 1.76,\ 0.47,\ 3.89,\ 5.45,\ 1.44,\ 1.23 .$$

For $\alpha = 0.1$: $\bar X(10) = 2.64$, $S^2(10) = 3.96$, $t_{0.05,\,9} = 1.833$, half-width $1.833\sqrt{3.96/10} = 1.15$, so a 90% CI is approximately $[1.49,\ 3.79]$.

6.3.1. Interpretation

“I am 95% sure that $[9, 11]$ covers $\mu$”? No.

$[9, 11]$ is an estimate: a fixed interval computed from one set of $n$ replications. $\mu$ is a fixed (unknown) number, so $[9, 11]$ either covers $\mu$ or it does not; nothing random is left. The formula $\big[\bar X(n) \pm t_{\alpha/2,\,n-1}\sqrt{S^2(n)/n}\big]$ is an estimator: a random interval, because $\bar X(n)$ and $S^2(n)$ are random quantities. If the whole procedure is repeated many times (each time with $n$ new replications), you get $[9, 11]$, $[8.5, 10.5]$, …, and 95% of these intervals cover $\mu$. The 95% describes the procedure, not one particular interval.

6.3.2. Quality of CIs

The CI assumes normal $X_j$’s; when they are not normal, the actual coverage can be below $1 - \alpha$. To check, use test models whose $\mu$ is known.

Ex: an M/M/1 queue ($\rho = 0.9$, starting empty and idle; $X_j$ = average delay in queue of the first 25 customers in replication $j$):

$$\mu = \mathbf{E}[\text{average delay in queue}] = 2.12 .$$

Construct 500 90% CIs for estimating $\mu$, each CI constructed based on $n$ simulation replications, and count the percentage of these 500 intervals that actually cover $\mu$.

$n$	actual coverage (500 CIs)	simulation (20,000 CIs)
5	88%	83.5%
10	86%	85.8%
20	89%	88.1%
40	92%	89.1%

Ex: reliability example: a system of three components that fails when component 1 fails or when both components 2 and 3 have failed; the component lifetimes are i.i.d. Weibull (shape 0.5, scale 1). $X_j$ = time to failure of the system in replication $j$, and $\mu = 0.78$ (exactly $7/9$). Perform the same experiment:

$n$	actual coverage (500 CIs)	simulation (20,000 CIs)
5	71%	71.4%
10	75%	75.9%
20	80%	80.1%
40	84%	83.7%

A coverage value based on only 500 intervals is accurate to about $\pm 2.6\%$ (the standard error of a proportion near 0.9 from 500 trials is $\sqrt{0.9 \times 0.1/500} = 1.3\%$). So the small ups and downs in the first table (88% → 86% → 89% → 92%) are within this noise; with 20,000 intervals the coverage rises steadily toward 90% as $n$ grows.

The reliability coverage is lower. Why is there such a big difference?

Fact: the quality of the confidence interval depends on both $n$ and the nature of $X_j$. If $X_j$ itself is the average of a large number of quantities, then the resulting CI in general has good coverage.

Why: the $t$-interval is exact only for normal $X_j$; otherwise it relies on the CLT to make $\bar X(n)$ approximately normal, and the further the distribution of $X_j$ is from normal (especially the more skewed it is), the larger $n$ must be. In the M/M/1 example, $X_j$ is an average of 25 delays and is already fairly close to normal (skewness 1.7). In the reliability example, $X_j$ is a single time to failure, not an average: most systems fail almost immediately but a few last very long (skewness 5.3), far from normal.

(The 25 delays within a run are positively correlated, so averaging them helps less than averaging 25 independent values would, but it still helps.)

6.3.3. Obtaining a specified precision

In the $100(1-\alpha)\%$ CI $\ \bar X(n) \pm t_{\alpha/2,\,n-1}\sqrt{S^2(n)/n}$:

• $\bar X(n)$ is the center;
• $1 - \alpha$ is the coverage probability;
• $n$ is the sample size (number of replications);
• $\delta = t_{\alpha/2,\,n-1}\sqrt{S^2(n)/n}$ is the half-width of the CI, because the CI is $[\text{center} - \text{half-width},\ \text{center} + \text{half-width}]$. The half-width measures the precision: the smaller, the more precise.

If I increase $1 - \alpha$ ($\alpha$ goes down), the tail area $\alpha/2$ gets smaller, so the critical value $t_{\alpha/2,\,n-1}$ is pushed to the right:

However, if $t$ gets larger, the interval becomes wider and the precision is lower. If the range is too large, it is not informative: e.g. a 100% CI for a salary is $[0, \infty)$, which is useless.

So for a fixed $n$, more confidence costs precision. The way to get both is to increase $n$: the half-width shrinks like $1/\sqrt n$ (made precise in the procedures below).

Ex: two systems

1. M/M/1, 2. M/D/1 (deterministic: every service time equals the mean service time of system 1). I want to do a comparison. Let $D_i$ = average delay in queue in the $i$-th system, $i = 1$ or $2$. Estimating $D_1 - D_2$ is the comparison: which system has the smaller delay? Do this by simulation.

$$\mu = \mathbf{E}[D_1 - D_2] = \mathbf{E}[D_1] - \mathbf{E}[D_2] .$$

Make $n$ replications of each system. Replication $j$ gives $D_{1j}$ and $D_{2j}$; let

$$Z_j = D_{1j} - D_{2j}, \qquad j = 1, \dots, n .$$

The replications are independent, so $Z_1, \dots, Z_n$ are i.i.d. with mean $\mu$: this is exactly the setting of §6.3 (estimating means) with $X_j$ replaced by $Z_j$. So $\bar Z(n)$ is an unbiased estimator of $\mu$, and with $S_Z^2(n)$ the sample variance of the $Z_j$’s, a

95% CI for $\mu$ is

$$\left[\bar Z(n) - t_{0.025,\,n-1}\sqrt{\frac{S_Z^2(n)}{n}},\;\; \bar Z(n) + t_{0.025,\,n-1}\sqrt{\frac{S_Z^2(n)}{n}}\,\right],$$

and where it lies relative to 0 gives the answer:

1. The CI lies entirely to the right of 0: $\mu > 0$, i.e. $\mathbf{E}[D_1] > \mathbf{E}[D_2]$, so system 2 has the smaller expected delay. The difference is statistically significant.
2. The CI lies entirely to the left of 0: $\mu < 0$, so system 1 has the smaller expected delay.
3. The CI contains 0: statistically insignificant. The data are consistent with $\mu = 0$, so we cannot tell which system is better. This does not show that $\mu = 0$; the interval may simply be too wide. The remedy is a smaller half-width, i.e. more replications. How many? That is what the procedures below answer.

E.g. both queues with $\rho = 0.9$ (arrival rate 1, mean service time 0.9), started empty; $D_i$ = average delay of the first 25 customers. With $n = 10$ replications of each system: $\bar Z(10) = 1.48$, $S_Z^2(10) = 6.90$, half-width $2.262\sqrt{6.90/10} = 1.88$, so the 95% CI is $[-0.40,\ 3.36]$. It contains 0: statistically insignificant (case 3), even though the true difference is $\mu = 0.71$ (system 2 is better).

Procedures for determining the number of replication runs

First case: $\alpha$ and $\beta$ (the desired half-width) are fixed; determine $n$. Set the half-width equal to $\beta$,

$$\beta = t_{\alpha/2,\,n-1}\sqrt{\frac{S^2(n)}{n}},$$

and solve for $n$:

$$n = \frac{S^2(n)\, t^2_{\alpha/2,\,n-1}}{\beta^2} .$$

Both sides contain $n$ ($S^2(n)$ and $t_{\alpha/2,\,n-1}$ depend on it), so in practice: make $n_0$ pilot replications, use $S^2(n_0)$ in place of $S^2(n)$, take the smallest $n$ with $t_{\alpha/2,\,n-1}\sqrt{S^2(n_0)/n} \le \beta$, and make $n - n_0$ more replications. (For large $n$, $t_{\alpha/2,\,n-1} \approx z_{\alpha/2}$, which gives $n \approx S^2(n_0)\, z^2_{\alpha/2}/\beta^2$ directly.)

Second case: $\alpha$ is fixed; reduce the half-width from $\beta_1$ (obtained with $n_1$ replications) to $\beta_2$. What is the additional number of replications needed?

$$\beta_1 = t_{\alpha/2,\,n_1-1}\sqrt{\frac{S^2(n_1)}{n_1}}, \qquad \beta_2 = t_{\alpha/2,\,n_2-1}\sqrt{\frac{S^2(n_2)}{n_2}} .$$

Assume $S^2(n_2) \approx S^2(n_1)$ (both estimate the same variance) and $t_{\alpha/2,\,n_2-1} \approx t_{\alpha/2,\,n_1-1}$. Dividing,

$$\frac{\beta_1}{\beta_2} \approx \sqrt{\frac{n_2}{n_1}} \quad\Rightarrow\quad n_2 \approx n_1\left(\frac{\beta_1}{\beta_2}\right)^2, \qquad \text{additional replications} = n_2 - n_1 \approx n_1\left[\left(\frac{\beta_1}{\beta_2}\right)^2 - 1\right].$$

E.g. halving the half-width ($\beta_2 = \beta_1/2$) needs $n_2 = 4n_1$, i.e. $3n_1$ more replications. In fact $t_{\alpha/2,\,n-1}$ decreases as $n$ grows, so this rule errs on the safe side.

E.g. (continued: $\alpha = 0.05$, target half-width $0.5$):

• First case, with the pilot $S_Z^2(10) = 6.90$: the smallest $n$ with $t_{0.025,\,n-1}\sqrt{6.90/n} \le 0.5$ is $n = 109$ (there $t_{0.025,\,108} = 1.98$, and $6.90 \times 1.98^2/0.5^2 = 108.4$).
• Second case, from $\beta_1 = 1.88$ with $n_1 = 10$ to $\beta_2 = 0.5$: $n_2 \approx 10\,(1.88/0.5)^2 = 141.2$, i.e. about 132 more replications. More than 109, because the rule keeps $t$ at $2.262$ instead of letting it drop to about $1.98$.

Running $n = 109$ replications gives $\bar Z(109) = 0.95$, $S_Z^2(109) = 4.44$, half-width $0.40$, and the 95% CI $[0.55,\ 1.35]$. It lies entirely to the right of 0 (case 1): the M/D/1 queue has the smaller expected delay. (The half-width came out below the target 0.5 because the pilot’s $S^2 = 6.90$ happened to be high: the procedures are only as good as the pilot’s variance estimate.)

Third scenario: precision is given in terms of the relative error. Fix $\alpha$ and the relative error $\gamma$; determine $n$. Set the half-width relative to the size of the estimate equal to $\gamma$:

$$\gamma = \frac{t_{\alpha/2,\,n-1}\sqrt{S^2(n)/n}}{\left|\bar X(n)\right|} \quad\Rightarrow\quad n = \frac{S^2(n)\, t^2_{\alpha/2,\,n-1}}{\gamma^2\, \bar X(n)^2} .$$

What it means: the half-width is required to be at most a fraction $\gamma$ of the estimate itself, so the CI is $\bar X(n)\,(1 \pm \gamma)$; with $\gamma = 0.1$, “$\mu$ is estimated to within about $\pm 10\%$”. This is useful when you do not know the scale of $\mu$ in advance: an absolute half-width $\beta = 0.5$ is very precise if $\mu \approx 100$ but useless if $\mu \approx 0.3$, while a relative error scales automatically. It is the first case with $\beta$ replaced by $\gamma\,|\bar X(n)|$, and it is solved the same way, using $\bar X(n_0)$ and $S^2(n_0)$ from a pilot run.

E.g. M/M/1 alone, $\mu = \mathbf{E}[D_1] = 2.12$: a pilot of $n_0 = 10$ replications gave $\bar X(10) = 2.17$ and $S^2(10) = 3.44$. For $\alpha = 0.05$ and $\gamma = 0.10$ the target half-width is $0.10 \times 2.17 = 0.217$, and the smallest $n$ with $t_{0.025,\,n-1}\sqrt{3.44/n} \le 0.217$ is $n = 282$ (formula: $3.44 \times 1.968^2/(0.1^2 \times 2.17^2) = 281.5$). In a simulation check, with 282 replications $\bar X$ was within 10% of $\mu$ in 96% of repetitions.

Two cautions. (1) The relative error is measured against $|\bar X(n)|$, not the unknown $|\mu|$; to guarantee $|\bar X(n) - \mu|/|\mu| \le \gamma$, the textbook uses the slightly smaller adjusted value $\gamma/(1+\gamma)$ in the formula. (2) Do not use a relative error for a quantity that may be close to 0, such as $\mathbf{E}[D_1 - D_2]$ above: as $\bar X(n) \to 0$, the target half-width $\gamma|\bar X(n)| \to 0$ and the required $n$ explodes.

Common random numbers. Above, the two systems used independent random numbers. Feeding both systems the same interarrival times in each replication (common random numbers, §4.2.1 advantage 4) makes $D_{1j}$ and $D_{2j}$ positively correlated (here the correlation is about 0.5), and

$$\operatorname{Var}(Z_j) = \operatorname{Var}(D_1) + \operatorname{Var}(D_2) - 2\operatorname{cov}(D_1, D_2)$$

drops from 3.79 to 2.23 in this example. Since the required $n$ is proportional to the variance (first case), that is about 40% fewer replications for the same half-width.

6.3.4. Choosing initial conditions in terminating simulations

The initial conditions affect the results of a terminating simulation, so they must be representative of the real system (§6.2). In general there are two approaches.

Ex: estimate the expected average delay in queue of bank customers who arrive and complete their delays between noon and 1 pm.

Noon to 1 pm is the lunch rush, so at noon the bank is usually already busy. Starting the simulation empty and idle at noon would miss the customers already waiting then, and would underestimate the delays (as in the figure of §6.1, a queue started empty has smaller delays for a while).

Approach 1: start the simulation empty at 9 am (when the bank really opens, so “empty” is realistic), run the model for 3 hrs, reset the statistical counters at 12:00 noon and collect data for the next hr. The 9 am to noon part is a warm-up that builds up a realistic state at noon; its statistics are discarded.

Approach 2: collect data on the number of customers present at 12 noon (over many days), use the data to fit a distribution, and draw from that distribution to initialize the simulation at noon in each replication; then collect data from noon to 1 pm.

• Approach 1 needs no extra data, but every replication spends 3 simulated hours only to reach noon, and the arrivals during those hours must be modeled correctly (e.g. with a time-varying arrival rate, §4.6).
• Approach 2 saves those 3 hours, but needs real data on the state at noon, and the number of customers may not describe the whole state (e.g. how much service the customers already in service still need; for exponential service times this does not matter, by the memoryless property).

In both approaches the replications stay i.i.d.: in approach 2, each replication draws its own initial number of customers independently.

6.3.5. Summary for terminating simulations

1. We must perform multiple independent simulation runs (replications); a single run is not enough (start of this chapter).
2. It is important to control the precision of the estimator, i.e. determine $n$ (§6.3.3).
3. Initial conditions must be carefully chosen (§6.3.4).

6.4. Statistical Analysis for Steady-State Parameters (Non-Terminating)

Let $Y_1, Y_2, \dots$ be the output r.v.s of a nonterminating simulation. $Y_i$ converges to $Y$ in distribution; more formally,

$$F_{Y_i}(y) \to F_Y(y) \quad \text{as } i \to \infty, \text{ for all } y,$$

where $F_{Y_i}(y) = P(Y_i \le y)$ is the transient distribution of $Y_i$ (written $F_i(y \mid I)$ in §6.1) and $F_Y$ is the steady-state distribution.

$\mathbf{E}[Y]$, $P(Y > s)$, $\operatorname{Var}(Y)$ are called steady-state parameters (§6.2): they are characteristics of the steady-state distribution $F_Y$. Goal: estimate a steady-state parameter, e.g. the steady-state mean

$$\mu = \mathbf{E}[Y] = \lim_{i \to \infty}\mathbf{E}[Y_i] .$$

Why is it difficult?

• Problem of initial transience: $\mathbf{E}[Y_i] \ne \mathbf{E}[Y]$, because the early $Y_i$’s come from transient distributions $F_{Y_i} \ne F_Y$. So for the plain average of a run, $\bar Y(m) = \frac{1}{m}\sum_{i=1}^m Y_i$, we get $\mathbf{E}[\bar Y(m)] \ne \mathbf{E}[Y]$: $\bar Y(m)$ is a biased estimator of $\mu$ (initialization bias). E.g. for the M/M/1 queue of §6.1 started empty ($\mu = 8.1$): $\mathbf{E}[\bar Y(100)] = 4.13$, $\mathbf{E}[\bar Y(1000)] = 7.30$, $\mathbf{E}[\bar Y(10000)] = 8.02$. The bias shrinks only slowly as the run gets longer.
• Within one run the $Y_i$’s are correlated, so the i.i.d. CI of §2.6 cannot be applied to them directly. (§6.4.2–§6.4.4 are three ways around this.)

6.4.1. Initial-data deletion

Given observations $Y_1, Y_2, \dots, Y_m$ of a run, delete the first $l$ observations and use

$$\bar Y(m, l) = \frac{1}{m - l}\sum_{i=l+1}^{m} Y_i, \qquad 1 \le l \le m - 1 .$$

The deleted observations $Y_1, \dots, Y_l$ are the warm-up period (the same idea as the warm-up in approach 1 of §6.3.4); $l = 0$ would give back the plain average $\bar Y(m)$.

Is this a better estimator? Why? The deleted observations are the ones most affected by the initial condition ($Y_1 = 0$ always, $Y_2$ is often $0$, …), i.e. the ones whose means $\mathbf{E}[Y_i]$ are farthest from $\mu$. If the process is (approximately) in steady state after observation $l$, every remaining $Y_i$ has $\mathbf{E}[Y_i] \approx \mu$, so $\mathbf{E}[\bar Y(m, l)] \approx \mu$: the bias is (almost) removed. The price: $\bar Y(m, l)$ averages only $m - l$ observations instead of $m$, so its variance is larger. So it is better in bias, worse in variance.

Want $\mathbf{E}[\bar Y(m, l)] \approx \mathbf{E}[Y]$. How to choose $l$? Trade-off:

• if $l$ is small ⇒ high bias (transient observations are still included);
• if $l$ is large ⇒ high variance (only a few observations are left).

E.g. the M/M/1 queue above, runs of length $m = 2000$ (bias computed numerically, standard deviations from 100,000 simulated runs):

$l$	bias $\mathbf{E}[\bar Y(m, l)] - \mu$	standard deviation of $\bar Y(m, l)$
0	$-0.405$	3.3
100	$-0.217$	3.4
250	$-0.110$	3.6
500	$-0.043$	3.9
1000	$-0.009$	4.7
1500	$-0.003$	5.9

The bias matters more than these numbers suggest, because the variance can be reduced by making more runs and the bias cannot. With $n$ independent runs, the estimates $\bar Y(m, l)$ of the runs are i.i.d., so they can be averaged and turned into a CI as in §6.3 (the replication/deletion approach, §6.4.2). Averaging divides the variance by $n$, but the expectation stays $\mathbf{E}[\bar Y(m, l)]$ for every $n$, so more runs only shrink the CI around the wrong value. E.g. with $l = 0$, the 95% CI covered $\mu = 8.1$ in 83% of experiments with $n = 10$ runs, 68% with $n = 100$ and 26% with $n = 400$; with $l = 500$ it covered $\mu$ in about 93% of experiments for both $n = 100$ and $n = 400$. So choose $l$ just large enough that the bias is negligible. That is what Welch’s procedure is for.

It is very difficult to determine $l$ from a single replication, because the observations are so variable: a value of $l$ that works well for one replication may be very bad for another. So Welch’s procedure uses several replications.

Welch’s procedure (for determining $l$). Welch’s procedure chooses $l$ so that $\mathbf{E}[Y_i] \approx \mu$ for $i > l$. The idea is to perform multiple runs, average across runs and then within runs (over neighboring observations) to construct a moving average process, plot the result and “eyeball” the value of $l$.

1. Make $n$ runs of the simulation ($n \ge 5$), each of length $m$ ($m$ much larger than the anticipated $l$). As in the table at the start of the chapter, $Y_{ji}$ = the $i$-th observation of run $j$:

$$\begin{array}{c|cccc} & Y_1 & Y_2 & \cdots & Y_m \\ \hline \text{run } 1 & Y_{11} & Y_{12} & \cdots & Y_{1m} \\ \text{run } 2 & Y_{21} & Y_{22} & \cdots & Y_{2m} \\ \vdots & \vdots & \vdots & & \vdots \\ \text{run } n & Y_{n1} & Y_{n2} & \cdots & Y_{nm} \end{array}$$

2. Average across runs (down each column):

$$\bar Y_i = \frac{1}{n}\sum_{j=1}^{n} Y_{ji}, \qquad i = 1, \dots, m, \qquad \mathbf{E}[\bar Y_i] = \mathbf{E}[Y_i] .$$

Also $\operatorname{Var}(\bar Y_i) = \operatorname{Var}(Y_i)/n$, since the runs are independent. So the averaged process $\bar Y_1, \bar Y_2, \dots, \bar Y_m$ has the same transient mean curve $\mathbf{E}[Y_i]$ as a single run, which is what we want to see, but only $1/n$ of the variance.
3. Construct a moving average process. $\bar Y_i$ still fluctuates a lot from one $i$ to the next; averaging neighboring values smooths out these fast fluctuations but keeps the slow trend $\mathbf{E}[Y_i]$. Specify the window $w > 0$ (an integer, $w \le \lfloor m/4 \rfloor$) and let

$$\bar Y_i(w) = \begin{cases} \dfrac{1}{2w+1}\displaystyle\sum_{s=-w}^{w} \bar Y_{i+s}, & i = w+1, \dots, m-w,\\[2.5ex] \dfrac{1}{2i-1}\displaystyle\sum_{s=-(i-1)}^{i-1} \bar Y_{i+s}, & i = 1, \dots, w. \end{cases}$$

For $i = w+1, \dots, m-w$, $\bar Y_i(w)$ is the simple average of the $2w+1$ values $\bar Y_{i-w}, \dots, \bar Y_{i+w}$ centered at $i$. For $i \le w$ there are only $i - 1$ values to the left of $i$, so the window shrinks to $i - 1$ values on each side: $2i - 1$ values, still centered at $i$, so that the smoothed curve is not shifted. E.g. $\bar Y_1(w) = \bar Y_1$, $\bar Y_2(w) = \frac{1}{3}(\bar Y_1 + \bar Y_2 + \bar Y_3)$, $\bar Y_3(w) = \frac{1}{5}(\bar Y_1 + \cdots + \bar Y_5)$. For $i > m - w$ there are fewer than $w$ values to the right, so $\bar Y_i(w)$ is not defined there.

4. Plot $\bar Y_i(w)$ for $i = 1, 2, \dots, m - w$ and eyeball a cutoff value for $l$: e.g. the value of $i$ such that beyond $i$, $\bar Y_i(w)$ does not change much (it appears to have converged).

E.g. the M/M/1 queue of §6.1 ($\rho = 0.9$, started empty), $m = 2000$, $w = 100$. This queue is very variable: with $n = 10$ runs the moving average still swings between roughly 4 and 12, too much to see where it levels off, so $n = 100$ runs were made. $\bar Y_i(w)$ rises until about $i = 600$ and after that only fluctuates around a constant level, so take $l \approx 600$. Since this is a test model, the choice can be checked: the bias of $\bar Y(2000, 600)$ is only $-0.03$, against $-0.405$ without deletion (table above).

Choosing $w$ and $n$: plot $\bar Y_i(w)$ for several values of $w$ and use the smallest $w$ for which the plot is reasonably smooth. If $w$ is too small, the plot is ragged and it is hard to see where it levels off; if $w$ is too large, the rise is smeared out and the shape of the transient is lost. If no $w$ gives a smooth plot, make more runs and repeat step 2 with all of them ($\operatorname{Var}(\bar Y_i) = \operatorname{Var}(Y_i)/n$ shrinks). Eyeballing is subjective; when in doubt, take $l$ a little larger, since extra variance can be fixed with more runs but bias cannot.

6.4.2. Replication/deletion approach for means

Once a good cutoff value (warm-up period) $l$ has been determined, the analysis proceeds like the terminating-simulation analysis of §6.3:

• make $n$ independent replications, each of length $m'$ and each with warm-up period $l$;
• compute the mean and variance estimates across replications;
• construct confidence intervals;
• often, use $n$ pilot runs to determine the warm-up period (Welch’s procedure, §6.4.1) and then throw away the first $l$ observations from the new production runs.

In formulas: from replication $j$ keep only the observations after the warm-up and average them,

$$X_j = \bar Y_j(m', l) = \frac{1}{m' - l}\sum_{i=l+1}^{m'} Y_{ji}, \qquad j = 1, \dots, n,$$

then compute $\bar X(n)$ and $S^2(n)$ from $X_1, \dots, X_n$, and the $100(1-\alpha)\%$ CI for $\mu$ is

$$\bar X(n) \pm t_{\alpha/2,\,n-1}\sqrt{\frac{S^2(n)}{n}} .$$

Why this works: the $X_j$’s are i.i.d., because the runs are independent copies of the same experiment (same initial conditions, separate random-number streams), and $\mathbf{E}[X_j] \approx \mu$, because the warm-up period was deleted. So the $X_j$’s play exactly the role of the $X_j$’s of §6.3, and everything there applies: the CI, its interpretation (§6.3.1), its quality (§6.3.2), and the procedures for choosing $n$ for a desired precision (§6.3.3). The correlation within a run no longer causes trouble: it is hidden inside each $X_j$, and only the $X_j$’s, one per run, enter the CI.

The production runs should be much longer than the warm-up ($m' \gg l$); they can be longer than the pilot runs.

E.g. the M/M/1 queue ($\rho = 0.9$, started empty, $\mu = 8.1$) with $l = 600$ from Welch’s procedure. $n = 10$ production runs of length $m' = 5000$ gave

$$X_j = 6.68,\ 9.13,\ 6.34,\ 10.06,\ 6.36,\ 6.90,\ 6.71,\ 9.70,\ 9.50,\ 5.61,$$

$\bar X(10) = 7.699$, $S^2(10) = 2.836$, half-width $2.262\sqrt{2.836/10} = 1.20$, so the 95% CI is $[6.49,\ 8.90]$, which covers $\mu$. Repeating the whole experiment 2000 times, the CI covered $\mu$ in 91.6% of experiments with $n = 10$ and 93.6% with $n = 20$ (without deletion, $l = 0$: 89.0% and 89.2%). The remaining shortfall from 95% is the skewness of the $X_j$’s, as in §6.3.2.

Pros and cons:

• + Simple: it is the familiar i.i.d. analysis, and the $X_j$’s are truly independent.
• − Every run has to pass through the transient: $n \times l$ observations are thrown away (here $10 \times 600 = 6000$ out of $50{,}000$).
• − If $l$ is chosen too small, the remaining bias is in every $X_j$ and does not average out over the runs (§6.4.1).

The second point motivates the batch means method: use one run, so the transient is passed only once.

6.4.3. The batch means method

The batch means method is the most popular alternative to replication/deletion. The idea is to perform one very long simulation run and estimate the steady-state parameters from that single run by dividing the output into batches.

Make one simulation run of length $m$, with $m$ very large. Divide the observations $Y_1, Y_2, \dots, Y_m$ into $n$ batches of size $k$, so that $m = nk$:

$$\underbrace{Y_1, \dots, Y_k}_{\text{batch } 1},\ \underbrace{Y_{k+1}, \dots, Y_{2k}}_{\text{batch } 2},\ \dots,\ \underbrace{Y_{(n-1)k+1}, \dots, Y_{nk}}_{\text{batch } n} .$$

(Here $n$ is the number of batches, not the number of runs: there is only one run.) Let

$$\bar Y_j(k) = \frac{1}{k}\sum_{i=(j-1)k+1}^{jk} Y_i$$

be the sample mean of the $k$ observations in the $j$-th batch. Then $\bar Y_1(k), \bar Y_2(k), \dots, \bar Y_n(k)$ are the batch means. Let

$$\bar Y(n, k) = \frac{1}{n}\sum_{j=1}^n \bar Y_j(k)$$

be the grand sample mean. This is used as a point estimator of $\mathbf{E}[Y]$. Note that

$$\bar Y(n, k) = \frac{1}{n}\sum_{j=1}^n \bar Y_j(k) = \frac{1}{n}\sum_{j=1}^n \frac{1}{k}\sum_{i \in \text{batch } j} Y_i = \frac{1}{nk}\sum_{i=1}^m Y_i = \frac{1}{m}\sum_{i=1}^m Y_i :$$

since all batches have the same size, the average of the batch means is just the average of all $m$ observations. Batching does not change the point estimate; it is needed only to estimate its variance.

When $m$ is large, the influence of the initial condition dies out: the biased early observations are only a small fraction of the $m$ observations. (One can also delete a warm-up period first and batch the rest, as in the figure and the example below.)

If $k$ is large enough, then the batch means $\bar Y_j(k)$ will be approximately uncorrelated and approximately normally distributed. Why:

• Uncorrelated: $Y_i$ and $Y_{i+h}$ are correlated, but the correlation dies out as the lag $h$ grows. Two neighboring batches are correlated only through the observations near their common boundary (the end of one batch and the beginning of the next). When $k$ is much longer than the “memory” of the process, these are a small part of each batch, so the correlation between batch means is small.
• Normal: each batch mean is an average of $k$ observations; by a central limit theorem for correlated observations, it is approximately normal for large $k$ (the same reason as in §6.3.2, where an $X_j$ that is an average of many observations gave good coverage).

For normal random variables, uncorrelated plus jointly normal implies independence. Thus, for sufficiently large $k$, it is reasonable to treat the batch means as i.i.d. normal random variables and construct confidence intervals for $\mathbf{E}[Y]$ exactly as in §6.3, with the batch means in the role of the $X_j$’s. An estimate of the variance (of one batch mean) is

$$S_{\bar Y}^2 = \frac{\sum_{j=1}^n \big(\bar Y_j(k) - \bar Y(n, k)\big)^2}{n - 1},$$

and the $100(1-\alpha)\%$ CI for $\mathbf{E}[Y]$ is

$$\bar Y(n, k) \pm t_{\alpha/2,\,n-1}\sqrt{\frac{S_{\bar Y}^2}{n}} .$$

Why not simply treat the $m$ observations $Y_i$ themselves as i.i.d., and use their sample variance divided by $m$? Because they are positively correlated, $\operatorname{Var}(\bar Y(m))$ is much larger than $\operatorname{Var}(Y)/m$ (positive covariances add to the variance of a sum: $\operatorname{Var}(Y_1 + \cdots + Y_m) = \sum_i \operatorname{Var}(Y_i) + 2\sum_{i < i'}\operatorname{cov}(Y_i, Y_{i'})$, §2.3), so that CI would be far too narrow. In the M/M/1 example below ($m = 120{,}000$), the naive i.i.d. 95% CI has an average half-width of only 0.05, while the actual standard deviation of $\bar Y(m)$ is 0.49: the naive standard error is about 19 times too small, and the naive CI covered $\mu$ in only 9% of experiments. Batching fixes this because each batch mean already contains the effect of the correlation within its batch.

Choosing the batch size. The advantage of batch means is that the initial transient is encountered only once. The difficulty is choosing the batch size $k$ (for a fixed $m$, choosing $k$ is the same as choosing the number of batches $n = m/k$):

• If $k$ is too small, the batch means are highly correlated and $S_{\bar Y}^2$ can be a poor estimator. (Positively correlated batch means look less spread out than independent ones would, so $S_{\bar Y}^2$ underestimates the variance and the CI is too narrow.)
• If $k$ is too large, there are fewer batches and the variance of the variance estimator is large: $S_{\bar Y}^2$ is based on only a few batch means and $t_{\alpha/2,\,n-1}$ is large, so the CI is wide and unstable.

A rule of thumb is that there is little benefit from more than 30 batches, but there should not be fewer than 10:

$$10 \le n \le 30 .$$

Why: M/M/1 ($\rho = 0.9$), warm-up 600 deleted, then one run of $m = 120{,}000$ cut into $n$ batches; 2000 repetitions, 95% CIs:

batches $n$	batch size $k$	correlation of neighboring batch means	$t_{0.025,\,n-1}$	average half-width	coverage
5	24,000	0.01	2.776	1.26	94.2%
10	12,000	0.02	2.262	1.04	93.3%
20	6,000	0.01	2.093	0.97	93.4%
30	4,000	0.03	2.045	0.95	93.2%
50	2,400	0.05	2.010	0.92	93.0%
100	1,200	0.11	1.984	0.86	91.1%
300	400	0.33	1.968	0.71	85.3%
600	200	0.51	1.964	0.58	76.9%

From 5 to about 20–30 batches, the CI gets narrower (the critical value drops from 2.78 to about 2.05, and $S_{\bar Y}^2$ is based on more batch means) at no cost in coverage. Beyond about 30, $t_{\alpha/2,\,n-1}$ hardly changes any more (it is already close to $z_{0.025} = 1.96$), so there is little left to gain; the further narrowing comes from correlated batch means making $S_{\bar Y}^2$ too small, and the coverage drops.

E.g. one run of the M/M/1 queue: warm-up 600 deleted, then $m = 120{,}000$ in $n = 20$ batches of $k = 6000$: $\bar Y(20, 6000) = 7.518$, $S_{\bar Y}^2 = 2.268$, half-width $2.093\sqrt{2.268/20} = 0.70$, so the 95% CI is $[6.81,\ 8.22]$, which covers $\mu = 8.1$.

6.4.4. Regenerative method

The regenerative method is similar in spirit to batch means. It tries to construct independent replications from a single run, using regeneration points. The difference: batch means cut the run at fixed places (every $k$ observations), so the pieces are only approximately independent; the regenerative method cuts the run at random places chosen so that the pieces are exactly i.i.d.

A regeneration point is a random time at which the process probabilistically restarts itself: from that time on, the process behaves (in distribution) exactly as if it had just been started, independently of what happened before. The process between two successive regeneration points is called a regeneration cycle.

Ex: single-server queue. A customer who arrives to an empty system finds the same situation every time: nobody waiting, the server idle, and the customer starts service at once. Everything after that (the following interarrival times and service times) is new and independent of the past. So the arrival times of customers who find the system empty are regeneration points. In terms of the delays, these are exactly the customers with $Y_i = 0$. A run started empty has $Y_1 = 0$, so it starts at a regeneration point.

Suppose the output process $Y_1, Y_2, \dots$ has regeneration points

$$1 \le B_1 < B_2 < \cdots .$$

At each regeneration point, the process after $B_j$ is independent of the process prior to $B_j$. Let

$$N_j = B_{j+1} - B_j$$

be the length of the $j$-th cycle (the number of observations in it), and let

$$Z_j = \sum_{i=B_j}^{B_{j+1}-1} Y_i$$

be the sum of the observations in the $j$-th cycle.

In the queue, cycle $j$ consists of the customers served in one busy period: $N_j$ = number of customers served in the busy period, $Z_j$ = their total delay. A cycle with $N_j = 1$ is a customer who arrives to an empty system and leaves before the next arrival, so $Z_j = 0$.

Then

$$U_j = (Z_j, N_j)^T$$

are i.i.d. bivariate random variables: different cycles are independent copies (each cycle starts from a fresh restart), but within a cycle $Z_j$ and $N_j$ are dependent (a long busy period, large $N_j$, also has large delays, large $Z_j$).

By the renewal-reward theorem,

$$\mu = \mathbf{E}[Y] = \lim_{m \to \infty}\frac{\sum_{i=1}^m Y_i}{m} = \frac{\mathbf{E}[Z_j]}{\mathbf{E}[N_j]}, \qquad \text{for all } j.$$

This is the long-run average cost or reward per unit time: each cycle earns a “reward” $Z_j$ and takes “time” $N_j$ (counted in observations). Why it holds: over the first $n'$ cycles there are $N_1 + \cdots + N_{n'}$ observations with total $Z_1 + \cdots + Z_{n'}$, so their average is

$$\frac{Z_1 + \cdots + Z_{n'}}{N_1 + \cdots + N_{n'}} = \frac{\frac{1}{n'}\sum_j Z_j}{\frac{1}{n'}\sum_j N_j} \;\to\; \frac{\mathbf{E}[Z_j]}{\mathbf{E}[N_j]} \quad \text{as } n' \to \infty,$$

by the SLLN (§2.4) applied separately to the i.i.d. $Z_j$’s and the i.i.d. $N_j$’s. E.g. in the M/M/1 queue, a busy period serves on average $\mathbf{E}[N_j] = 1/(1-\rho) = 10$ customers, so $\mathbf{E}[Z_j] = \mu\,\mathbf{E}[N_j] = 81$.

Assume we simulate exactly $n'$ regeneration cycles. Let

$$\bar Z(n') = \frac{1}{n'}\sum_{j=1}^{n'} Z_j, \qquad \bar N(n') = \frac{1}{n'}\sum_{j=1}^{n'} N_j .$$

A point estimator of $\mu$ is

$$\hat\mu(n') = \frac{\bar Z(n')}{\bar N(n')} .$$

This estimator is biased, because the expected value of a ratio is not the ratio of the expected values:

$$\mathbf{E}\left[\frac{\bar Z(n')}{\bar N(n')}\right] \ne \frac{\mathbf{E}[\bar Z(n')]}{\mathbf{E}[\bar N(n')]} = \mu,$$

but it is strongly consistent:

$$\hat\mu(n') \to \mu \quad \text{with probability 1 as } n' \to \infty$$

(this is the calculation above). In the M/M/1 example, the average of $\hat\mu(n')$ over many runs was 6.97 for $n' = 100$ cycles, 7.93 for $n' = 1000$ and 8.08 for $n' = 5000$: a large bias for few cycles, disappearing as $n'$ grows.

Confidence interval. The problem: $\hat\mu$ is a ratio of two averages, so the CLT cannot be applied to it directly. The trick is to turn the ratio into an average. Let $\Sigma$ be the covariance matrix of $U_j = (Z_j, N_j)^T$:

$$\Sigma = \begin{pmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{pmatrix}, \qquad \sigma_{11} = \operatorname{Var}(Z_j),\quad \sigma_{22} = \operatorname{Var}(N_j),\quad \sigma_{12} = \sigma_{21} = \mathbf{E}\big[(Z_j - \mathbf{E}[Z_j])(N_j - \mathbf{E}[N_j])\big].$$

Let

$$V_j = Z_j - \mu N_j .$$

The $V_j$’s are i.i.d. (functions of the i.i.d. $U_j$’s), and

$$\mathbf{E}[V_j] = \mathbf{E}[Z_j] - \mu\,\mathbf{E}[N_j] = 0$$

because $\mu = \mathbf{E}[Z_j]/\mathbf{E}[N_j]$, and

$$\sigma_V^2 = \operatorname{Var}(V_j) = \operatorname{Var}(Z_j - \mu N_j) = \sigma_{11} - 2\mu\sigma_{12} + \mu^2\sigma_{22}$$

(the variance of a linear combination: $\operatorname{Var}(aZ + bN) = a^2\operatorname{Var}(Z) + 2ab\operatorname{cov}(Z, N) + b^2\operatorname{Var}(N)$ with $a = 1$, $b = -\mu$). By the central limit theorem,

$$\frac{\bar V(n')}{\sqrt{\sigma_V^2/n'}} \approx N(0, 1), \qquad \bar V(n') = \frac{1}{n'}\sum_{j=1}^{n'} V_j .$$

Why this gives a CI for $\mu$: the average of the $V_j$’s is

$$\bar V(n') = \bar Z(n') - \mu\,\bar N(n') = \bar N(n')\big(\hat\mu(n') - \mu\big),$$

so the statement above says

$$\frac{\hat\mu(n') - \mu}{\sqrt{\sigma_V^2/n'}\,\big/\,\bar N(n')} \approx N(0, 1),$$

i.e. $\hat\mu(n')$ is approximately normal around $\mu$ with standard deviation $\sqrt{\sigma_V^2/n'}/\bar N(n')$.

Since $\sigma_V^2$ is unknown, estimate it using the sample covariance matrix

$$\hat\Sigma(n') = \frac{\sum_{j=1}^{n'}\big(U_j - \bar U(n')\big)\big(U_j - \bar U(n')\big)^T}{n' - 1} = \begin{pmatrix} \hat\sigma_{11}(n') & \hat\sigma_{12}(n') \\ \hat\sigma_{21}(n') & \hat\sigma_{22}(n') \end{pmatrix},$$

where $\bar U(n') = (\bar Z(n'), \bar N(n'))^T$; so $\hat\sigma_{11}$ is the sample variance of the $Z_j$’s, $\hat\sigma_{22}$ the sample variance of the $N_j$’s, and $\hat\sigma_{12}$ their sample covariance. Then, with $\hat\mu(n')$ in place of the unknown $\mu$,

$$\hat\sigma_V^2(n') = \hat\sigma_{11}(n') - 2\hat\mu(n')\hat\sigma_{12}(n') + \big[\hat\mu(n')\big]^2\hat\sigma_{22}(n') .$$

An approximate $100(1-\alpha)\%$ confidence interval for $\mu$ is

$$\left[\hat\mu(n') - z_{\alpha/2}\frac{\sqrt{\hat\sigma_V^2(n')/n'}}{\bar N(n')},\;\; \hat\mu(n') + z_{\alpha/2}\frac{\sqrt{\hat\sigma_V^2(n')/n'}}{\bar N(n')}\right].$$

($z$, not $t$: the whole result is a large-$n'$ approximation.)

E.g. one run of the M/M/1 queue, started empty, $n' = 1000$ cycles (8438 customers). The first cycles were $(Z_j, N_j) = (77.35, 18), (0, 1), (0, 1), (0.35, 2), (30.99, 11), (3.07, 7), \dots$ Over all 1000 cycles: $\bar Z = 66.695$, $\bar N = 8.438$, so $\hat\mu = 66.695/8.438 = 7.904$; $\hat\sigma_{11} = 396{,}988$, $\hat\sigma_{12} = 22{,}976$, $\hat\sigma_{22} = 1431.8$, so

$$\hat\sigma_V^2 = 396{,}988 - 2(7.904)(22{,}976) + 7.904^2(1431.8) \approx 123{,}230,$$

$$\text{half-width} = 1.96\,\frac{\sqrt{123{,}230/1000}}{8.438} = 2.58,$$

and the 95% CI is $[5.33,\ 10.48]$. It covers $\mu = 8.1$ but is wide. Over many runs, this CI covered $\mu$ in 55% of runs with $n' = 100$ cycles, 81% with $n' = 1000$ and 91% with $n' = 5000$. The cycles of this queue are extremely variable (most busy periods are short, a few are very long with huge $Z_j$), so the ratio estimator is biased and the normal approximation is poor unless there are many cycles.

Difficulty of the regenerative method: many processes may not have regeneration points. Even if they do, the expected cycle length may be so large that only a few cycles can be simulated (e.g. in a large network of queues, the moment when all queues are empty at the same time is rare).

Ex: single-server queue regeneration points. Possible regeneration points include the arrival times of customers who find the system empty (above). Two choices that do not work:

• Times immediately after busy cycles end (the moment the system becomes empty) are not good regeneration points if the time until the next arrival is not identically distributed. At that moment the next arrival is already “on its way”: the remaining time until it depends on how long ago the previous arrival happened, which differs from one busy period to the next. (Only for exponential interarrival times does this not matter, by the memoryless property.) An arrival to an empty system, by contrast, starts a brand-new interarrival time.
• Arrival times of customers who find one customer in the system are also not regeneration points, because the remaining service time of the customer in service is differently distributed: it depends on how long that service has already been going on. The process after such an arrival is not a fresh copy.

What the regenerative method offers in exchange: if the run starts at a regeneration point (e.g. an empty queue), there is no initial transient problem at all, there is no batch size to choose, and the cycles are exactly i.i.d.

6.4.5. Sequential procedure

A sequential procedure sequentially determines the length of a single run needed to construct an acceptable confidence interval.

Explanation. In §6.3.3 we chose the number of replications $n$ to reach a desired precision. For a single long run (e.g. batch means), the corresponding question is: how long must the run be? This cannot be computed in advance, because the variance of $\bar Y(m)$ depends on the unknown correlation structure of the process. A sequential procedure decides while the run is going. One simple version:

1. Run the simulation for an initial length $m_1$ (after the warm-up).
2. Construct a CI from the data so far, e.g. by batch means with a fixed number of batches.
3. If the CI is acceptable (e.g. its relative half-width is at most $\gamma$, as in §6.3.3), stop. Otherwise continue the same run to a longer length (e.g. double it) and go back to step 2.

(The procedures in the textbook are more refined, e.g. they also check that the batch means are approximately uncorrelated before accepting the CI.)

E.g. M/M/1 queue, warm-up 600, 20 batches, start at $m = 10{,}000$ and double $m$ until the 95% CI has relative half-width $\le \gamma = 0.1$. Over 2000 repetitions the procedure stopped at $m = 80{,}000$ (9% of runs), $160{,}000$ (41%), $320{,}000$ (44%) or $640{,}000$ (5%), and rarely earlier. The estimate was within 10% of $\mu$ in 95.8% of runs, but the final CI covered $\mu$ in only 90.2% (with a fixed $m = 320{,}000$: 94.8%). The stopping rule is slightly optimistic: it tends to stop at moments when $S_{\bar Y}^2$ happens to be small, i.e. when the CI happens to look narrow.

6.4.6. Summary for steady-state simulations

The problem (§6.4): within one run the observations are correlated and the early ones are biased. All methods construct (approximately) i.i.d. pieces and then use the i.i.d. analysis of §6.3:

method	runs	warm-up	i.i.d. pieces	CI	main difficulty
replication/deletion (§6.4.2)	$n$ independent runs	deleted in every run ($l$ from Welch, §6.4.1)	$X_j = \bar Y_j(m', l)$, exactly i.i.d.	$t$-CI from the $X_j$’s	choosing $l$; the transient is paid $n$ times
batch means (§6.4.3)	one long run	passed once	batch means $\bar Y_j(k)$, approximately i.i.d.	$t$-CI from the batch means	choosing $k$ ($10 \le n \le 30$ batches)
regenerative (§6.4.4)	one run	none, if started at a regeneration point	cycles $(Z_j, N_j)$, exactly i.i.d.	ratio CI with $z$	needs regeneration points; many cycles needed
sequential (§6.4.5)	one run, extended as needed	passed once	e.g. batch means	stops when the CI is precise enough	stopping rule is slightly optimistic

6.5. Multiple Performance Measures

In practice, we are often interested in several performance measures simultaneously (e.g. the average delay in queue, the average number in queue and the server utilization of the same system). Thus we may construct several confidence intervals:

$$P(V_1 \in I_1) = 1 - \alpha_1, \qquad P(V_2 \in I_2) = 1 - \alpha_2, \qquad \dots, \qquad P(V_k \in I_k) = 1 - \alpha_k,$$

where $V_s$ is the true value of the $s$-th performance measure and $I_s$ is its CI. The question is how these individual confidence intervals affect the overall confidence level:

$$P(V_s \in I_s \text{ for all } s = 1, \dots, k) = \ ?$$

This is the probability that all $k$ confidence intervals simultaneously contain their corresponding performance measures. It matters whenever conclusions are drawn from all the intervals together (“the delay is in $I_1$ and the utilization is in $I_2$”).

If the performance measures (more precisely, their CIs) are independent, then

$$P(V_s \in I_s,\ \forall s) = \prod_{s=1}^k P(V_s \in I_s) = (1 - \alpha_1)(1 - \alpha_2)\cdots(1 - \alpha_k) .$$

E.g. three independent 95% CIs: $0.95^3 = 0.857$, so all three are correct together only with probability 85.7%, not 95%.

In practice, however, performance measures are unlikely to be independent (they are computed from the same simulation runs). A conservative approach uses Bonferroni’s inequality, which holds for any dependence:

$$\begin{aligned} P(V_s \in I_s,\ \forall s) &= P(V_1 \in I_1 \cap V_2 \in I_2 \cap \cdots \cap V_k \in I_k) \\ &= 1 - P(V_1 \notin I_1 \cup V_2 \notin I_2 \cup \cdots \cup V_k \notin I_k) \\ &\ge 1 - \sum_{s=1}^k P(V_s \notin I_s) \\ &\ge 1 - \sum_{s=1}^k \alpha_s . \end{aligned}$$

Step by step:

1. “All intervals are correct” is the intersection of the events $\{V_s \in I_s\}$.
2. Complement: its opposite, “at least one interval is wrong”, is the union of the events $\{V_s \notin I_s\}$ (De Morgan’s law), so $P(\text{all correct}) = 1 - P(\text{at least one wrong})$.
3. Union bound (Boole’s inequality): $P(A_1 \cup \cdots \cup A_k) \le \sum_s P(A_s)$, because the sum counts the overlaps of the events more than once. Subtracting a larger number gives “$\ge$”.
4. Each interval is wrong with probability $P(V_s \notin I_s) = \alpha_s$ (at most $\alpha_s$ when the CI is approximate).

Therefore, the overall confidence level is at least

$$100\left(1 - \sum_{s=1}^k \alpha_s\right)\% .$$

E.g. three 95% CIs: overall at least $100(1 - 3 \times 0.05)\% = 85\%$ (compare 85.7% in the independent case). In the other direction, to get an overall level of at least $100(1-\alpha)\%$, choose the $\alpha_s$’s so that $\sum_s \alpha_s = \alpha$, e.g. $\alpha_s = \alpha/k$: for $k = 2$ measures and overall 90%, make two 95% CIs; for $k = 5$ and overall 90%, five 98% CIs. The price is that each individual CI is wider (a larger critical value $t_{\alpha_s/2,\,n-1}$), so with many measures the individual CIs can become too wide to be useful.