Output analysis is the examination of data generated by a simulation. Its purpose is either to predic the performance ofa system or to compare the performance oftwo or more alternative system designs. This chapter deals with estimating absolute performance, by which we mean estimating the value of one or more system performance measures as opposed to the comparison of two or more systems, in other words relative performance.

The need for statistical output analysis is based on the observation that the output data from a simulation exhibits random variability when random-number generators are used to produce the values of the input variables—that is, two different streams or sequences of random numbers wil produce two sets ofoutputs, which (probably) will differ. Ifthe performance ofthe system is measured by a parameter $\theta ,$ the result ofa set ofsimulation experiments will be an estimator θ ofθ. The precision of the estimator $\widehat {\theta}$ can be measured by the standard error of $\widehat {\theta}$-or by the width of a confidence interva - -for θ. The purpose of the statistical analysis is either to estimate this standard error or confidence interval or to figure out the number of observations required to achieve a standard error or confidence interval of a given size—or both.

Consider a typical output variable Y, the total cost per week of an inventory system; Y should be treated as a random variable with an unknown distribution. A simulation run of length 1 week provides a single sample observation from the population of all possible observations on Y. By increasing the run length, the sample size can be increased to n observations,$Y_{1} , Y_{2} , \ldots , Y_{n} .$, based on a run length of n weeks. However, these observations do not constitute a random sample, in the classic sense, because they are not statistically independent. In this case, the inventory on hand a the end of one week is the beginning inventory on hand for the next week, and so the value of $Y_{i}$ has some influence on the value of $Y_{i + 1}$. Thus, the sequence of random variables $Y_{1} , Y_{2} , \dots , Y_{n}$ could be autocorrelated (i.e., correlated with itself). This autocorrelation, which is a measure of a lack of statistical independence, means that classic methods of statistics, which assume independence, are not directly applicable to the analysis of these output data. The methods must be properly modified and the simulation experiments properly designed for valid inferences to be made.

In addition to the autocorrelation present in most simulation output data, the specification o the initial conditions of the system at time 0 can pose a problem for the simulation analyst and could influence the output data. By “time $0^{\circ}$ we mean whatever point in time the beginning of the simulation run represents. For example, the inventory on hand and the number of backorders at time 0 (Monday morning) would most likely influence the value of $Y_{1}$, the total cost for week 1. Because of the autocorrelation, these initial conditions would also influence the costs $Y_{2} , \dots , Y_{n}$ for subsequent weeks. The specified initial conditions, if not chosen well, can have an effect on estimation of the steady-state (long-run) performance of a simulation model. For purposes of statistical analysis, the effect of the initial conditions is that the output observations might not be identically distributed and that the initial observations might not be representative of the steady-state behavior of the system.

Section 1 distinguishes between two types of simulation—transient versus steady state—and defines commonly used measures of system performance for each type of simulation. Section 2 illustrates by example the inherent variability in a stochastic (i.e., probabilistic) discrete-event simulation and thereby demonstrates the need for a statistical analysis of the output. Section 3 covers the statistical estimation of performance measures. Section 4 discusses the analysis of terminating simulations, and Section 5 the analysis of steady-state simulations. To support the calculations in this chapter a spreadsheet called SimulationTools.xls is provided at www.bcnn.net. SimulationTools.xls is a menu-driven application that allows the user to paste in their own data—generated by any simulation application—and perform statistical analyses on it. A detailed user guide is also available on the book web site, and SimulationTools.xls has integrated help.

1 Types of Simulations with Respect to Output Analysis

In analyzing simulation output data, a distinction is made between terminating or transient simulations and steady-state simulations.

A terminating simulation is one that runs for some duration of time $T_{E}$, where E is a specified event (or set of events) that stops the simulation. Such a simulated system opens at time 0 under well-specified initial conditions and closes at the stopping time $T_{E}$. The next four examples are terminating simulations.

Example 1

The Shady Grove Bank opens at 8:30 A.M. (time 0) with no customers present and 8 of the 11 tellers working (initial conditions) and closes at 4:30 P.M. (time $T_{E} = 480$ minutes). Here, the event E is merely the fact that the bank has been open for 480 minutes. The simulation analyst is interested in modeling the interaction between customers and tellers over the entire day, including the effect of starting up and closing down at the end of the day.

Example 2

Consider the Shady Grove Bank of Example 1, but restricted to the period from 11:30 A.M. (time 0) to 1:30 P.M., when it is especially busy. The simulation run length is $T_{E} = 120$ minutes. The initial conditions at time 0 (11:30 A.M.) could be specified in essentially two ways: (a) the real system could be observed at 11:30 on a number of different days and a distribution of number of customers in system (at 11:30 A.M.) could be estimated, then these data could be used to load the simulation model with customers at time 0; or (b) the model could be simulated from 8:30 A.M. to 11:30 A.M. without collecting output statistics, and the ending conditions at 11:30 A.M. used as initial conditions for the 11:30 A.M. to 1:30 P.M. simulation.

Example 3

A communications system consists of several components plus several backup components. It is represented schematically in Figure 1. Consider the system over a period of time $T_{E} ,$, until the system fails. The stopping event E is defined by $E = \{A$ fails, or D fails, or (B and C both fail) . Initia conditions are that all components are new at time 0.

Figure 1 Example of a communications system.

Notice that, in the bank model of Example 1, the stopping time $T_{E} = 480$ minutes is known, but in Example 3, the stopping time $T_{E}$ is generally unpredictable in advance; in fact,$T_{E}$ is probably the output variable of interest, as it represents the total time until the system breaks down. One goal o the simulation might be to estimate $E ( T_{E} )$, the mean time to system failure.

Example 4

A manufacturing process runs continuously from Monday mornings until Saturday mornings. The first shift of each work week is used to load inventory buffers and chemical tanks with the components

and catalysts needed to make the final product. These components and catalysts are made continually throughout the week, except for the last shift Friday night, which is used for cleanup and maintenance. Thus, most inventory buffers are near empty at the end of the week. During the first shift on Monday, a buffer stock is built up to cover the eventuality of breakdown in some part of the process. It is desired to simulate this system during the first shift (time 0 to time $T_{E} = 8$ hours) to study various scheduling policies for loading inventory buffers.

In the simulating of a terminating system, the initial conditions of the system at time 0 must be specified, and the stopping time $T_{E} - \mathrm {{o r}}$, alternatively, the stopping event E—must be well defined. Although it is certainly true that the Shady Grove Bank in Example 1 will open again the next day, the simulation analyst has chosen to consider it a terminating system because the object of interest is one day’s operation, including start up and close down. On the other hand, if the simulation analys were interested in some other aspect of the bank’s operations, such as the flow of money or operation of automated teller machines, then the system might be considered as a nonterminating one. Similar comments apply to the communications system of Example 3. If the failed component were replaced and the system continued to operate, and, most important, if the simulation analyst were interested in studying its long-run behavior, it might be considered as a nonterminating system. In Example 3, however, interest is in its short-run behavior, from time 0 until the first system failure at time $T_{E}$ Therefore whether a simulation is considered to be terminating depends on both the obiectives of the simulation study and the nature of the system.

Example 4 is a terminating system, too. It is also an example of a transient (or nonstationary) simulation: the variables of interest are the in-process inventory levels, which are increasing from zero or near zero (at time 0) to full or near full (at time 8 hours).

A nonterminating system is a system that runs continuously, or at least over a very long period of time. Examples include assembly lines that shut down infrequently, continuous production systems of many different types, telephone and other communications systems such as the Internet, hospita emergency rooms, police dispatching and patrolling operations, fire departments, and continuously operating computer networks.

A simulation of a nonterminating system starts at simulation time 0 under initial conditions defined by the analyst and runs for some analyst-specified period of time $T_{E}$. (Significant problems arise concerning the specification of these initial and stopping conditions, problems that we discuss later.) Usually, the analyst wants to study steady-state, or long-run, properties of the system—tha is, properties that are not influenced by the initial conditions of the model at time 0. A steadystate simulation is a simulation whose objective is to study long-run, or steady-state, behavior of a nonterminating system.

The next two examples are steady-state simulations.

Example 5

Consider the manufacturing process ofExample 4, beginning with the second shift, when the complete production process is under way. It is desired to estimate long-run production levels and production efficiencies. For the relatively long period of 13 shifts, this may be considered as a steady-state simulation. To obtain sufficiently precise estimates of production efficiency and other response variables, the analyst could decide to simulate for any length of time,$T_{E}$(even longer than 13 shifts). That is,$T_{E}$ is not determined by the nature of the problem (as it was in terminating simulations); rather, it is set by the analyst as one parameter in the design of the simulation experiment.

Example 6

HAL Inc., a large web-based order-processing company, has many customers worldwide. Thus, its large computer system with many servers, workstations, and peripherals runs continuously, 24 hours per day. To handle an increased work load, HAL is considering additional CPUs, memory, and storage devices in various configurations. Although the load on HAL’s computers varies throughou the day, management wants the system to be able to accommodate sustained periods of peak load. Furthermore, the time frame in which HAL’s business will change in any substantial way is unknown, so there is no fixed planning horizon. Thus, a steady-state simulation at peak-load conditions is appropriate. HAL systems staff develops a simulation model of the existing system with the current peak work load and then explores several possibilities for expanding capacity. HAL is interested in long-run average throughput and utilization of each computer. The stopping time $T_{E}$ is determined not by the nature of the problem, but rather by the simulation analyst, either arbitrarily or with a certain statistical precision in mind.

2 Stochastic Nature of Output Data

Consider one run of a simulation model over a period of time [0,$T_{E} ]$. Since the model is an inputoutput transformation, and since some of the model input variables are random variables, it follows that the model output variables are random variables. Two examples are now given to illustrate the nature of the output data from stochastic simulations and to give a preliminary discussion of severa important properties of these data. Do not be concerned if some of these properties and the associated terminology are not entirely clear on a first reading. They will be explained carefully later in the chapter.

Example 7: Software Made Personal

Software Made Personal (SMP) customizes software products for clients in two areas: financia tracking and contact management. They have a customer support call center that handles questions for owners of their software from 8 A.M. to 4 P.M. Eastern Time. When a customer calls they use an automated system to select among the two product lines. Each product line has its own operators, and if an appropriate operator is available then the call is immediately routed to the operator; otherwise, the caller is placed in a hold queue. SMP is hoping to reduce the total number of operators they need by cross-training them so that they can answer calls for any product line. This is expected to increase the time to process a call by about 10%. Before considering reducing operators, however, SMP wants to know what their quality of service would be if they simply cross-train the operators they currently have. Using historical data, incoming calls have been modeled as a nonstationary Poisson arriva process and distributions have been fit for operator time to answer questions. Their quality of service measures are caller response time (time from initiating a call until speaking to an operator) and the number of callers on hold. Response-time differences of more than 30 seconds matter.

This is clearly a terminating simulation experiment with time 0 corresponding to 8 A.M. and $T_{E}$ corresponding to 4 P.M. (another way to terminate the simulation is when the last call that arrives before 4 P.M. is terminated, in which case $T_{E}$ is a random value). Table 1 shows the results from 4 replications of the SMP simulation; each row presents the average number of minutes callers spen on hold, and the average number of callers on hold, during the course of a simulated 8-hour day.

Table 1 Results of Four Independent Replications of the SMP Call Center

The table illustrates that there is significant variation in these averages from day to day, implying that a substantial number of replications may be required to precisely estimate, say, the mean values $w_{Q} = E ( \widehat {w} _ {Q} )$ and $L _ {\mathcal {Q}} = E ( \widehat {L} _ {\mathcal {Q}} )$. Even with a large number of replications there will still be some - -estimation error, which is why we want to report not just a single value (the point estimate), but also a measure of its error in the form of a confidence interval or standard error.

These questions are addressed in Section 4 for terminating simulations such as Example 7. Classic methods of statistics may be used because $\widehat {w} _ {Q 1} , \widehat {w} _ {Q 2} , \widehat {w} _ {Q 3}$, and $\widehat {w} _ {Q 4}$ constitute a random sample—that is, they are independent and identically distributed. In addition,$w_{Q} = E ( \widehat {w} _ {Q r} )$ is the parameter being estimated, so each $\widehat {w} _ {Q r}$ is an unbiased estimate of the true mean waiting time $w_{Q}$ The analysis of Example 7 is considered in Example 10 of Section 4. A survey of statistical methods applicable to terminating simulations is given by Law [1980]. Additional guidance may be found in Alexopoulos and Seila [1998], Kleijnen [1987], Law [2007], and Nelson [2001].

The next example illustrates the effects of correlation and initial conditions on the estimation o long-run mean measures of performance of a system.

Example 8

Semiconductor wafer fabrication is an amazingly expensive process involving hundreds of millions of dollars in processing equipment, material handling, and human resources. FastChip, Inc. wants to use simulation to design a new “fab” by evaluating the steady-state mean cycle time (time from wafer release to completion) that they can expect under a particular loading of two distinct products, called C-chip and D-chip. The process of wafer fabrication consists of two basic steps, diffusion and lithography, each of which contains many substeps, and the fabrication process requires multiple passes through these steps. Product will be released in cassettes at the rate of 1 cassette/hour, 7 days a week, 24-hours a day, to achieve a desired throughput of 1 cassette/hour. The hope is that the curren design will achieve long-run average cycle times of less than 45 hours for the C-chip and 30 hours for the D-chip.

Figure 2 shows the first 30 Chip-C cycle times from one run of the simulation model of the fab. The goal is to estimate the long-run average cycle time as, conceptually, the number of cassettes produced goes to infinity. Notice that the first few cycle times seem lower than the ones that follow; this occurs because the simulation was started with the fab “empty and idle,” which is not representative of long-run or steady-state conditions. Also notice that when there is a low cycle time (say, around 40 hours), it tends to be followed by several additional low cycle times. Thus, the cycle time of successive cassettes may be correlated. These features suggest that the cycle-time data are not a random sample, so classic statistical methods may not apply. The analysis of Example 8 is considered in Section 5.

Figure 2 First 30 Chip-C cycle times from one simulation run.

3 Absolute Measures of Performance and Their Estimation

Consider the estimation of a performance parameter θ (or φ) of a simulated system. It is desired to have a point estimate and an interval estimate of θ (or φ). The length of the interval estimate is a measure of the error in the point estimate. The simulation output data are of the form $\{Y_{1} , Y_{2} , \ldots , Y_{n} \}$ for estimating $\theta ;$ we refer to such output data as discrete-time data, because the index n is discrete valued. The simulation output data are of the form $\{Y ( t ) , 0 \leq t \leq T_{E} \}$ for estimating φ; we refer to such output data as continuous-time data, because the index t is continuous valued. For example,$Y_{i}$ might be the delay of customer i, or the total cost in week $i ; Y ( t )$ might be the queue length at time t, or the number of backlogged orders at time t. The parameter $\theta$ is an ordinary mean;$\phi$ will be referred to as a time-weighted mean. Whether we call the performance parameter $\theta$ or $\phi$ does not really matter; we use two different symbols here simply to provide a distinction between ordinary means and time-weighted means. To make the discussions more concrete, we consider output measures of queueing systems.

3.1 Point Estimation

The point estimator of θ, based on the data $\{Y_{1} , \ldots , Y_{n} \}$, is defined by

where $\widehat {\theta}$ is a sample mean based on a sample of size n. Computer simulation languages may refer to -this as a “discrete-time,” “collect,” “tally,” or “observational” statistic.

The point estimator θ is said to be unbiased for θ if its expected value is $\theta ;$ that ${\mathrm{i s}} ,$ if

In general, however,

and $E ( \widehat {\theta} ) - \theta$ is called the bias in the point estimator ${\widehat {\theta}} .$ It is desirable to have estimators that are -unbiased,$^ {\mathrm{o r ,}}$-if this is not possible, to have a small bias relative to the magnitude of θ. Examples of estimators of the form of Equation (1) include $\widehat {w}$ and $\widehat {w} _ {Q}$ of equations $\widehat {w} = \textstyle \frac {1} {N} \sum_{i = 1} ^ {N} W_{i}$ and $\widehat {w} _ {Q} =$$\frac {1} {N} \sum_{i = 1} ^ {N} \ : W_{i} ^ {Q} \longrightarrow w_{Q}$ as $N \longrightarrow \infty$, in which case $Y_{i}$ is the time spent in the (sub)system by customer i. The point estimator of $\phi ,$ based on the data $\{Y ( t ) , 0 \leq t \leq T_{E} \}$, where $T_{E}$ is the simulation run length, is defined by

and is called a time average of $Y ( t )$ over $[ 0 , T_{E} ]$. Simulation languages may refer to this as a “continuous-time,” “discrete-change,” or “time-persistent” statistic. In general,

and $\widehat {\phi}$ is said to be biased for φ. Again, we would like to obtain unbiased or low-bias estimators. -Examples of time averages include $\widehat {L}$ and $\widehat {L} _ {Q}$ of Equations $\begin{array} {r} {\widehat {L} = \frac {1} {T} \int_{0} ^ {T} L ( t ) d t \longrightarrow L} \end{array}$ as $T \longrightarrow \infty$ and $\begin{array} {r} {\widehat {L} _ {Q} = \frac {1} {T} \displaystyle \sum_{i = 0} ^ {\infty} i T_{i} ^ {Q} = \frac {1} {T} \int_{0} ^ {T} L_{Q} ( t ) d t \longrightarrow L_{Q}} \end{array}$ as $T \longrightarrow \infty$ and $Y_{j}$ of Equation (11.4)

Generally, θ and φ are regarded as mean measures of performance of the system being simulated. Other measures usually can be put into this common framework. For example, consider estimation o the proportion of days on which sales are lost through an out-of-stock situation. In the simulation, le

With n equal to the total number of days ${\widehat {\theta}} ,$, as defined by Equation (1), is a point estimator of $\theta_{;}$, the -proportion of out-of-stock days. For a second example, consider estimation of the proportion of time queue length is greater than 10 customers. If $L_{Q} ( t )$ represents simulated queue length at time $t ,$ then (in the simulation) define

Then ${\widehat {\phi}} ,$ as defined by Equation (4), is a point estimator of $\phi ,$, the proportion of time that the queue -length is greater than 10 customers. Thus, estimation of proportions or probabilities is a special case of the estimation of means.

A performance measure that does not fit this common framework is a quantile or percentile. Quantiles describe the level of performance that can be delivered with a given probability $p .$ For instance, suppose that Y represents the delay in queue that a customer experiences in a service system, measured in minutes. Then the $p = 0 . 85$ quantile of Y is the value θ such that

As a percentage,$\theta$ is the 100pth or 85th percentile of customer delay. Therefore, 85% of all customers will experience a delay of θ minutes or less. Stated differently, a customer has only a 0.15 probabilit of experiencing a delay of longer than θ minutes. A widely used performance measure is the median, which is the 0.5 quantile or 50th percentile.

The problem of estimating a quantile is the inverse of the problem of estimating a proportion or probability. In estimating a proportion,$\theta$ is given and the probability $p$ is to be estimated; however, in estimating a quantile,$p$ is given and θ is to be estimated.

The most intuitive method for estimating a quantile is to form a histogram of the observed values of Y, then find a value $\widehat {\theta}$ such that 100p% of the histogram is to the left of (smaller than) ${\widehat {\theta}} .$ For instance, if we observe $n = 250$ customer delays $\{Y_{1} , \ldots , Y_{250} \}$-, then an estimate of the 85th percentile of delay is a value $\widehat {\theta}$ such that $( 0 . 85 ) ( 250 ) = 212 . 5 \approx 213$ of the observed values are less than or equal to $\theta .$-An obvious estimate is, therefore, to set $\widehat {\theta}$ equal to the 213th smallest value -in the sample (this requires sorting the data). When the output is a continuous-time process, such as the queue-length process $\{L_{Q} ( t ) , 0 \leq t \leq T_{E} \}$, then a histogram gives the fraction of time tha the process spent at each possible level (queue length, in this example). However, the method for quantile estimation remains the same: Find a value $\widehat {\theta}$ such that 100p% of the histogram is to the left of ${\widehat {\theta}} .$.

3.2 Confidence-Interval Estimation

To understand confidence intervals fully, it is important to understand the difference between a measure of error and a measure of risk. One way to make the difference clear is to contrast a confidence interval with a prediction interval (which is another useful output-analysis tool).

Both confidence intervals and prediction intervals are based on the premise that the data being produced by the simulation is represented well by a probability model. Suppose that model is the normal distribution with mean θ and variance $\sigma^{2}$, both unknown. To make the example concrete, let $Y_{i}$ be the average cycle time for parts produced on the ith replication (representing a day of production) of the simulation. Therefore,$\theta$ is the mathematical expectation of $Y_{i}$, and σ is represents the day-to-day variation of the average cycle time.

Suppose our goal is to estimate θ. If we are planning to be in business for a long time, producing parts day after day, then $\theta$ is a relevant parameter, because it is the long-run mean daily cycle time. Our average cycle time will vary from day to day, but over the long run the average of the averages will be close to θ.

The natural estimator for $\theta$ is the overall sample mean of n independent replications,$\bar {Y} \ =$$\scriptstyle \sum_{i = 1} ^ {n} Y_{i} / n$. But $\bar {Y}$ is not θ, it is an estimate, based on a sample, and it has error. A confidence interva is a measure of that error. Let

be the sample variance across the R replications. The usual confidence interval, which assumes the $Y_{i}$ are normally distributed, is

where $t_{\alpha / 2 , R - 1}$ is the quantile of the t distribution with $R - 1$ degrees of freedom that cuts off $\alpha / 2$ the area of each tail. (See Table A.5 in the Appendix). We cannot know for certain exactly how far $\bar {Y}$ is from $\theta ,$ but the confidence interval attempts to bound that error. Unfortunately, the confidence interval itself may be wrong. A confidence level, such as 95%, tells us how much we can trust the interval to actually bound the error between Y¯ and θ. The more replications we make, the less error there is in ${\bar {Y}} ,$ and our confidence interval reflects that because $t_{\alpha / 2 , n - 1} S / {\sqrt {n}}$ will tend to get smaller as n increases, converging to 0 as n goes to infinity.

Now suppose we need to make a promise about what the average cycle time will be on a particular day. A good guess is our estimator ${\bar {Y}} ,$, but it is unlikely to be exactly right. Even θ itself, which is the center of the distribution, is not likely to be the actual average cycle time on any particular day, because the daily average cycle time varies. A prediction interval, on the other hand, is designed to be wide enough to contain the actual average cycle time on any particular day, with high probability.

The normal-theory prediction interval is

The length of this interval will not go to 0 as n increases. In fact, in the limit it becomes

to reflect the fact that, no matter how much we simulate, our daily average cycle time still varies.

In summary, a prediction interval is a measure of risk, and a confidence interval is a measure of error. We can simulate away error by making more and more replications, but we can never simulate away risk, which is an inherent part of the system. We can, however, do a better job of evaluating risk by making more replications.

Example 9

Suppose that the overall average of the average cycle time on 120 replications of a manufacturing simulation is 5.80 hours, with a sample standard deviation of 1.60 hours. Since $t_{0 . 025 , 119} = 1 . 98 ,$, a 95% confidence interval for the long-run expected daily average cycle time is $5 . 80 {\pm} 1 . 98 ( 1 . 60 / \sqrt {120} )$ or $5 . 80 \pm 0 . 29$ hours. Thus, our best guess of the long-run average of the daily average cycle times is 5.80 hours, but there could be as much as 0.29 hours error in this estimate.

On any particular day, we are 95% confident that the average cycle time for all parts produced on that day will be

or $5 . 80 \pm 3 . 18$ hours. The 3.18 hours reflects the inherent variability in the daily average cycle times and the fact that we want to be 95% confident of covering the actual average cycle time on a particular day, rather than simply covering the long-run average.

A caution: In neither case (prediction interval or confidence interval) can we make statements about the cycle time for individual products because we observed the daily average. To make valid statements about the cycle times of individual products, the analysis must be based on individua cycle times.

4 Output Analysis for Terminating Simulations

Consider a terminating simulation that runs over a simulated time interval $[ 0 , T_{E} ]$ and results in observations $Y_{1} , \dots , Y_{n}$. The sample size n may be a fixed number, or it may be a random variable (say, the number of observations that occur during time $T_{E} )$. A common goal in simulation is to estimate

When the output data are of the form $\{Y ( t ) , 0 \leq t \leq T_{E} \}$, the goal is to estimate

The method used in each case is the method of independent replications. The simulation is repeated a total of R times, each run using a different random number stream and independently chosen initia conditions (which includes the case that all runs have identical initial conditions). We now address this problem.

4.1 Statistical Background

Perhaps the most confusing aspect of simulation output analysis is distinguishing within-replication data from across-replication data, and understanding the properties and uses of each. The issue can be further confused by the fact that simulation languages often provide only summary measures, like sample means, sample variances, and confidence intervals, rather than all of the raw data. Sometimes these summary measures are all that the simulation language provides, without a lot of extra work.

To illustrate the key ideas, think in terms of the simulation of a manufacturing system and two performance measures of that system, the cycle time for parts (time from release into the factory until completion) and the work in process (WIP, the total number of parts in the factory at any time). In computer applications, these two measures could correspond to the response time and the length of the task queue at the CPU; in a service application, they could be the time to fulfill a customer’s request and the number of requests on the “to do” list; in a supply-chain application, they could be the order fill time and the inventory level. Similar measures appear in many systems.

Here is the usual set up for something like cycle time: Let $Y_{i j}$ be the cycle time for the jth part produced in the ith replication. If each replication represents two shifts of production, then the number of parts produced in each replication might differ. Table 2 shows, symbolically, the results of R replications.

Table 2 Within and Across-Replication Cycle-Time Data

The across-replication data are formed by summarizing within-replication data:^1${\bar {Y}} _ {i \cdot}$ is the sample mean of the $n_{i}$ cycle times from the ith replication,$S_{i} ^ {2}$ is the sample variance of the same data, and

is a confidence-interval half-width based on this dataset.

From the across-replication data, we compute overall statistics, the average of the daily cycle time averages

the sample variance of the daily cycle time averages

and finally, the confidence-interval half-width

The quantity $S / \sqrt {R}$ is the standard error, which is sometimes interpreted as the average error in ${\bar {Y}} .$ as an estimator of θ. Notice that $S^{2}$ is not the average of the within-replication sample variances $S_{i} ^ {2}$; rather, it is the sample variance of the within-replication averages $\bar {Y} _ {1 \cdot} , \bar {Y} _ {2 \cdot} , \dots , \bar {Y} _ {R \cdot}$ Table 3 Within and Across-Replication WIP Data

Within a replication, work in process is a continuous-time output, denoted $Y_{i} ( t )$. The stopping time for the ith replication $T _ {E_{i}}$ could be a random variable, in general; in this example, it is the end of the second shift. Table 3 is an abstract representation of the data produced

The within-replication sample mean and variance are defined appropriately for continuous-time data:

and

A definition for $H_{i}$ is more problematic, but, to be concrete, take it to be

It is difficult to conceive of a situation in which $H_{i}$ is relevant, a topic we discuss later. Although the definitions of the within-replication data change for continuous-time data, the across-replication statistics are unchanged, and this is a critical observation.

Here are the key points that must be understood:

The overall sample average ${\bar {Y}} _ {\cdot \cdot}$ and the individual replication sample averages ${\bar {Y}} _ {i}$ are always unbiased estimators of the expected daily average cycle time or daily average WIP.

Across-replication data are independent, since they are based on different random numbers; are identically distributed, since we are running the same model on each replication; and tend to be normally distributed, if they are averages of within-replication data, as they are here. This implies that the confidence interval ${\bar {Y}} . . \pm H$ is often pretty good.

Within-replication data, on the other hand, might have none of these properties. The individua cycle times may not be identically distributed, if the first few parts of the day find the system empty; they are almost certainly not independent, because one part follows another; and whether they are normally distributed is difficult to know in advance. For this reason,$S_{i} ^ {2}$ and $H_{i \cdot}$, which are computed under the assumption of independent and identically distributed $\mathrm{( i . i . d . )}$ data, tend not to be useful, although there are exceptions.

There are situations in which Y¯ and $Y_{i \cdot}$ are valid estimators of the expected cycle time for an individual part or the expected WIP at any point in time, rather than the daily average. (See Section 5 on steady-state simulations.) Even when this is the case, the confidence interva ${\bar {Y}} . . \pm H$ is valid, and $\bar {Y} _ {i \cdot} \pm H_{i}$ is not. The difficulty occurs because $S_{i} ^ {2}$ is a reasonable estimator of the variance of the cycle time, but $S_{i} ^ {2} / n_{i}$ and $S_{i} ^ {2} / T _ {E_{i}}$ are not good estimators of the $\mathrm{V a r} [ \bar {Y} _ {i \cdot} ] ;$ more on this in Section 5.2.

Example 10: Software Made Personal (cont.)

Consider Example 7, the SMP call center problem, with the data for $R = 4$ replications given in Table 1. The four number-on-hold estimates $\widehat {L} _ {Q r}$ are time averages of the form of Equation (11), where $Y_{r} ( t ) = L _ {\mathcal {O}} ( t )$-, the number of customers on hold at time t. Similarly, the four average system times $\widehat {w} _ {Q 1} , \dots , \widehat {w} _ {Q 4}$ are analogous to ${\bar {Y}} _ {r} .$ of Table 2, where $Y_{r i}$ is the actual time spent on hold by - -customer i on replication $r .$ Suppose that the analyst desires a 95% confidence interval for the mean time spent on hold,$w_{Q}$ Using Equation (8), compute an overall point estimator

Using Equation (9), compute its estimated variance:

Thus, the standard error of $\widehat {w} _ {Q} = 2 . 64$ is estimated by $\mathrm{s . e .} ( \widehat {w} _ {O} ) = S / \sqrt {4} = 1 . 14$. Obtain $t_{0 . 025 , 3} =$- -3.18 from Table A.5, and compute the 95% confidence interval half-width by Equation (10) as

giving $2 . 64 \pm 3 . 62$ minutes with 95% confidence. Clearly 4 replications are far too few to make any sort of useful estimate of $w_{Q}$, since a negative mean waiting time is impossible. Later in the chapter we return to this problem and determine how many replications are required to get a meaningfu estimate.

4.2 Confidence Intervals with Specified Precision

By Expression (10), the half-length H of a $100 ( 1 - \alpha ) \%$ confidence interval for a mean θ, based on the t distribution, is given by

where $S^{2}$ is the sample variance and R is the number of replications. In simulation we seldom need to select R and just accept the confidence interval half-length that results. Instead, we can drive the number of replications to be large enough so that H is small enough to facilitate the decision that the simulation is supposed to support. The acceptable size of H depends on the problem at hand: i might be $5000 in a financial simulation, 6 minutes for product cycle time, or 1 customer on hold in a call center.

Suppose that an error criterion of  is specified; in other words, it is desired to estimate θ by ${\bar {Y}} .$ to within  with high probability—say, at least 1 α. Thus, it is desired that a sufficiently large sample size R be taken to satisfy

When the sample size R is fixed, no guarantee can be given for the resulting error. But if the sample size can be increased, an error criterion (such as $\epsilon = \$5000$, 6 minutes, or 1 customer) can be specified.

Assume that an initial sample of size $R_{0}$ replications has been observed—that is, the simulation analyst initially makes $R_{0}$ independent replications. We must have $R_{0} \ge 2$, with 10 or more being desirable. The $R_{0}$ replications will be used to obtain an initial estimate $S_{0} ^ {2}$ of the population variance $\sigma^{2}$. To meet the half-length criterion, a sample size R must be chosen such that $R \geq R_{0}$ and

Solving for R in Inequality (14) shows that R is the smallest integer satisfying $R \geq R_{0}$ and

An initial estimate for R is given by

where $z_{\alpha / 2}$ is the $100 ( 1 - \alpha / 2 )$ percentage point of the standard normal distribution from Table A.3 in the Appendix. And since $t_{\alpha / 2 , R - 1} \approx z_{\alpha / 2}$ for large R, say,$R \geq 50$, the second inequality for R is adequate when R is large. After determining the final sample size R, collect $R - R_{0}$ additional observations (i.e., make $R - R_{0}$ additional replications, or start over and make R total replications) and form the 100(1 α)% confidence interval for θ by

where ${\bar {Y}} .$ and $S^{2}$ are computed on the basis of all R replications,${\bar {Y}} .$ by Equation (8), and $S^{2}$ by Equation (9). The half-length of the confidence interval given by Inequality (17) should be approximately  or smaller; however, with the additional $R - R_{0}$ observations, the variance estimator $S^{2}$ could differ somewhat from the initial estimate $S_{0} ^ {2} ,$ possibly causing the half-length to be greater than desired. If the confidence interval from Inequality (17) is too large, the procedure may be repeated, using Inequality (15), to determine an even larger sample size.

Example 11

Suppose that it is desired to estimate the mean waiting time in SMP’s hold queue in Example 7 to within 0.5 minutes with 95% confidence. We will illustrate this calculation using the “Get Sample $\mathrm {S i z e^{\gamma}}$ utility of SimulationTools.xls. The input dialog box is shown in Figure 3 and the results are displayed in Figure 4. Notice that “Absolute Precision” means that we directly specify the half width of the confidence interval ( 0.5 minutes in this case); “Relative Precision” specifies the length of the half width as a fraction of the mean.

An initial sample of size $R_{0} ~ = ~ 20$ is taken, and the program determines that $R \ = \ 175$ is the smallest integer satisfying Inequality (16). Therefore,$R - R_{0} = 175 - 20 = 155$ additional replications are needed. After obtaining the additional outputs, we would again need to compute the half-width H to ensure that it is as small as is desired. In this case the result—which can also be found at www.bcnn.net—was $2 . 66 \pm 0 . 34$ minutes

Figure 3 Input Dialog for the ‘‘Get Sample Size’’ utility of SimulationTools.xls

Figure 4 Sample Size Calculation using ‘‘Get Sample Size’’ in SimulationTools.xl

4.3 Quantiles

To present the interval estimator for quantiles, it is helpful to review the interval estimator for a mean in the special case when the mean represents a proportion or probability $p .$ In this book, we have chosen to treat a proportion or probability as just a special case of a mean. However, in many statistics texts, probabilities are treated separately.

When the number of independent replications $Y_{1} , \dots , Y_{R}$ is large enough that $t_{\alpha / 2 , n - 1} \doteq z_{\alpha / 2}$ the confidence interval for a probability $p$ is often written as

where $\widehat {p}$ is the sample proportion (algebra shows that this formula for the half-width is precisely -equivalent to Equation (10) when used in estimating a proportion).

As mentioned in Section 3, the quantile-estimation problem is the inverse of the probabilityestimation problem: Find θ such that $\operatorname* {P r} \{Y \leq \theta \} = p .$. Thus, to estimate the $p$ quantile, we find tha value $\widehat {\theta}$ such that 100p% of the data in a histogram of Y is to the left of ${\widehat {\theta}} ,$, or, stated differently, the -npth smallest value of $Y_{1} , \dots , Y_{R}$ Extending this idea, an approximate $( 1 - \alpha ) 100 \%$ confidence interval for $\theta$ can be obtained by finding two values:$\theta_{\ell}$ that cuts off $100 p_{\ell}$% of the histogram and $\theta_{u}$ that cuts off $100 p_{u} \%$ of the

histogram, where

(Recall that we know $p . )$ In terms of sorted values,$\widehat {\theta_{\ell}}$ is the $R p_{\ell}$ smallest value (rounded down) and $\widehat {\theta_{u}}$ is the $R p_{u}$ smallest value (rounded up), of $Y_{1} , \dots , Y_{R}$

Example 12

Suppose that we want to estimate the 0.75 quantile of average daily waiting time in the hold queue for SMP and form a 95% confidence interval for it, based on $R = 300$ replications. The data may be found in the spreadsheet Quantile.xls at www.bcnn.net. The point estimator is $\widehat \theta = 3 . 3 \dot {9}$ minutes, because that is the $300 \times 0 . 75 = 225 \mathrm{t h}$-smallest value of the sorted data. The interpretation is that, over the long run, on 75% of the days the average time customers spend on hold will be less than or equal to 3.39 minutes.

To obtain the confidence interval, we first compute

The lower bound ofthe confidence interval is $\widehat {\theta_{\ell}} = 3 . 13$ minutes (the $300 \times p_{\ell} = 2$10th smallest value, -rounding down); the upper bound of the confidence interval is $\widehat {\theta_{u}} = 3$.98 minutes (the $300 \times p_{u} =$-240th smallest value, rounding up). Thus, with 95% confidence, the point estimate of 3.39 minutes is no more than max $\{3 . 98 - 3 . 39 , 3 . 39 - 3 . 13 \} = 0 . 59$ minutes off the true 0.75 quantile.

It is important to understand the 3.39 minutes is an estimate of the 0.75 quantile (75th percentile) of daily average customer delay, not an individual customer’s delay. Because the call center starts the day empty and idle, and the customer load on the system varies throughout the day, the 0.75 quantile of individual customer waiting time varies throughout the day and cannot be given by a single number.

4.4 Estimating Probabilities and Quantiles from Summary Data

Knowing the equation for the confidence interval half-width is important if all that the simulation software provides is Y¯ and H and you need to work out the number of replications required to get a prespecified precision, or if you need to estimate a probability or quantile. You know the number of replications, so the sample standard deviation can be extracted from H by using the formula

With this information, the method in Section 4.2 can be employed.

The more difficult problem is estimating a probability or quantile from summary data. When al we have available is the sample mean and confidence-interval halfwidth (which gives us the sample standard deviation), then one approach is to use a normal-theory approximation for the probabilities or quantiles we desire, specifically

and

The following example illustrates how this is done.

Example 13

From 25 replications of a manufacturing simulation, a 90% confidence interval for the daily average WIP is 218 32. What is the probability that the daily average WIP is less than 350? What is the 85th percentile of daily average WIP?

First, we extract the standard deviation:

Then, we use the normal approximations and Table A.3 to get

and

There are shortcomings to obtaining our probabilities and quantiles this way. The approximation depends heavily on whether the output variable of interest is normally distributed. If the outpu variable itself is not an average, then this approximation is suspect. Therefore, we expect the approximation to work well for statements about the average daily cycle time, for instance, but very poorly for the cycle time of an individual part.

5 Output Analysis for Steady-State Simulations

Consider a single run of a simulation model whose purpose is to estimate a steady-state, or long-run, characteristic of the system. Suppose that the single run produces observations $Y_{1} , Y_{2} , \ldots$, which, generally, are samples of an autocorrelated time series. The steady-state (or long-run) measure o performance, θ, is defined by

with probability 1, where the value of θ is independent of the initial conditions. (The phrase “with probability $1^{\circ}$ means that essentially all simulations of the model, using different random numbers, will produce a series $Y_{i} , i = 1 , 2 , . .$. whose sample average converges to θ.) For example, if $Y_{i}$ was the time customer i spent talking to an operator, then θ would be the long-run average time a customer spends talking to an operator; and, because θ is defined as a limit, it is independent of the call center’s conditions at time 0. Similarly, the steady-state performance for a continuous-time output measure $\{Y ( t ) , t \geq 0 \}$, such as the number of customers in the call center’s hold queue, is defined as

with probability 1.

Of course, the simulation analyst could decide to stop the simulation after some number of observations, say, n, have been collected; or the simulation analyst could decide to simulate for some length of time $T_{E}$ that determines n (although n may vary from run to run). The sample size n (or $T_{E} )$ is a design choice; it is not inherently determined by the nature of the problem. The simulation analyst will choose simulation run length (n or $T_{E} )$ with several considerations in mind:

1. Any bias in the point estimator that is due to artificial or arbitrary initial conditions. (The bias can be severe if run length is too short, but generally it decreases as run length increases.)

2. The desired precision of the point estimator, as measured by the standard error or confidence interval half-width.

3. Budget constraints on the time available to execute the simulation.

The next subsection discusses initialization bias and the following subsections outline two methods of estimating point-estimator variability. For clarity of presentation, we discuss only estimation of θ from a discrete-time output process. Thus, when discussing one replication (or run), the notation

will be used; if several replications have been made, the output data for replication r will be denoted by

5.1 Initialization Bias in Steady-State Simulations

There are several methods of reducing the point-estimator bias caused by using artificial and unrealistic initial conditions in a steady-state simulation. The first method is to initialize the simulation in a state that is more representative of long-run conditions. This method is sometimes called intelligent initialization. Examples include

1. setting the inventory levels, number of backorders, and number of items on order and their arrival dates in an inventory simulation;

2. placing customers in queue and in service in a queueing simulation;

3. having some components fail or degrade in a reliability simulation.

There are at least two ways to specify the initial conditions intelligently. If the system exists, collect data on it and use these data to specify more nearly typical initial conditions. This method sometimes requires a large data-collection effort. In addition, if the system being modeled does no exist—for example, if it is a variant of an existing system—this method is impossible to implement. Nevertheless, it is recommended that simulation analysts use any available data on existing systems to help initialize the simulation, as this will usually be better than assuming the system to be “completely stocked,” “empty and idle,” or “brand new” at time 0.

A related idea is to obtain initial conditions from a second model of the system that has been simplified enough to make it mathematically solvable. The queueing models are very useful for this purpose. The simplified model can be solved to find long-run expected or most likely conditions— such as the expected number of customers in the queue—and these conditions can be used to initialize the simulation.

A second method to reduce the impact of initial conditions, possibly used in conjunction with the first, is to divide each simulation run into two phases: first, an initialization phase, from time 0 to time $T_{0} ,$ followed by a data-collection phase from time $T_{0}$ to the stopping time $T_{0} + T_{E} ;$ that is, the simulation begins at time 0 under specified initial conditions $I_{0}$ and runs for a specified period of time $T_{0}$. Data collection on the response variables of interest does not begin until time $T_{0}$ and continues until time $T_{0} + T_{E}$. The choice of $T_{0}$ is quite important, because the system state at time $T_{0} ,$ denoted by $I ,$ should be more nearly representative of steady-state behavior than are the origina initial conditions $I_{0}$ at time 0. In addition, the length $T_{E}$ of the data-collection phase should be long enough to guarantee sufficiently precise estimates of steady-state behavior. Notice that the system state I at time $T_{0}$ is a random variable and to say that the system has reached an approximate steady state is to say that the probability distribution of the system state at time $T_{0}$ is sufficiently close to the steady-state probability distribution as to make the bias in point estimates of response variables negligible. Figure 5 illustrates the two phases of a steady-state simulation. The effect of starting a simulation run of a queueing system in the empty and idle state, as well as a useful plot to aid the simulation analyst in choosing an appropriate value for $T_{0} ,$, are given in the following example.

Figure 5 Initialization and data collection phases of a steady-state simulation run.

Consider the FastChip wafer fabrication problem discussed in Example 8. Suppose that a total of R 10 independent replications were made, each run long enough so that 250 cycle times were collected for each of Chip-C and Chip-D. We will focus on Chip-C.

Normally we average all data within each replication to obtain a replication average. However, our goal at this stage is to identify the trend in the data due to initialization bias and find out when i dissipates. To do this, we will average corresponding cycle times across replications and plot them (this idea is usually attributed to Welch [1983]). Such averages are known as ensemble averages. Specifically, let $Y_{r 1} , Y_{r 2} , \dots , Y_{r 250}$ be the 250 Chip-C cycle times, in order of completion, from replication r. For the jth cycle time, define the ensemble average across all R replications to be

(R 10 here). The ensemble averages are plotted as the solid line in Figure 6 (the dashed lines will be discussed below). This plot was produced by the SimulationTools.xls spreadsheet, and we will refer to it as a mean plot.

Notice the gradual upward tread at the beginning. The simulation analyst may suspect that this is due to the downward bias in these estimators, which in turn is due to the fab being empty and idle at time 0. As time becomes larger (so the number of cycle times recorded increases), the effect of the initial conditions on later observations lessens and the observations appear to vary around a common mean. When the simulation analyst feels that this point has been reached, then the data-collection phase begins. Here approximately 120 Chip-C cycle times might be adequate.

Although we have determined a deletion amount of 120 cycle times, the initialization phase is almost always specified in terms of simulation time $T_{0} ,$ because we need to consider the initialization phase for each performance measure, and they may differ. In the FastChip example we would also

Figure 6 Mean Plot Across 10 Replications of 250 Product C Cycle Times

need to examine the cycle times for Chip-D. If we specified the initialization phase as a count for each output, then we would have to track these counts. When the initialization phase is specified by time, on the other hand, then a single event can be scheduled at time $T_{0} ,$, and data only retained beyond that point.

To specify $T_{0}$ by time we have to convert the appropriate counts into approximate times and use the largest of these. For instance, for the FastChip simulation 120 Chip-C cycle times corresponds to approximately 200 hours of simulation time, since cassettes are released at 1 per hour and 60% o the releases are Chip-C: 120 cycle times/0.6 releases/hour 200 hours (we could also increase this somewhat to account for the time it takes the 120th Chip-C release to traverse the fab). This turns out to be larger than the initialization phase for Chip-D, so we set $T_{0} = 200$ hours.

Plotting ensemble averages is extra work (although SimulationTools.xls makes it easier), so it might be tempting to take a shortcut. For instance, one might plot the output from only replication 1,$Y_{11} , Y_{12} , \dots , Y_{1 n}$, or the cumulative average from a replication 1

for $\ell = 1 , 2 , \ldots , n .$. Figure 7 shows both. Notice that the raw data from a single replication can be too highly variable to detect the trend, while the cumulative plot is clearly biased low because it retains all of the data from the beginning of the run. Thus, these shortcuts should be avoided and the mean plot used instead.

Figure 7 Plot of 250 Product C Cycle Times from 1 Replication, Raw and Cumulative Average

Some additional factors to keep in mind when using the mean plot to determine the initialization phase:

1. When first starting to examine the initialization phase, a run length and number of replications will have to be guessed. Ensemble averages, such as Figure 6, will reveal a smoother and more precise trend as the number of replications R is increased, and it may have to be increased several times until the trend is clear. It is also possible that the run length guess is not long enough to get past the initialization phase (this would be indicated if the trend in the mean plot does not disappear), so the run length may have to be increased as well.

2. Ensemble averages can also be smoothed by plotting a moving average rather than the original ensemble averages. In a moving average, each plotted point is actually the average of severa adjacent ensemble averages. Specifically, the jth plot point would be

for some $m \geq 1$, rather than the original ensemble average ${\bar {Y}} _ {\cdot j}$. The value of m is typically chosen by trial and error until a smooth plot is obtained. The mean plot in Figure 6 was obtained with $m = 2$3. Since each ensemble average (or moving average) is the sample mean of i.i.d. observations across R replications, a confidence interval based on the t distribution can be placed around each point, as shown by the dashed lines in Figure $^ {6 ,}$ and these intervals can be used to judge whether or not the plot is precise enough to decide that bias has diminished. This is the preferred method to determine a deletion point.

4. Cumulative averages, such as in Figure 7, become less variable as more data are averaged. Therefore, it is expected that the left side of the curve will always be less smooth than the righ side. Remember that cumulative averages tend to converge more slowly to long-run performance than do ensemble averages, because cumulative averages contain all observations, including the most biased ones from the beginning of the run. For this reason, cumulative averages should be used only if it is not feasible to compute ensemble averages, such as when only a single replication is possible.

There has been no shortage of solutions to the initialization-bias problem. Unfortunately, for every “solution” that works well in some situations, there are other situations in which either it is not applicable or it performs poorly. Important ideas include testing for bias (e.g., Kelton and Law [1983], Schruben [1980], Goldsman, Schruben, and Swain [1994]); modeling the bias (e.g., Snel and Schruben [1985]); and randomly sampling the initial conditions on multiple replications (e.g., Kelton [1989]).

5.2 Error Estimation for Steady-State Simulation

If $\{Y_{1} , \ldots , Y_{n} \}$ are not statistically independent, then $S^{2} / n ,$, given by Equation $( 9 ) ,$, is a biased estimator of the true variance $\nabla ( \widehat {\theta} )$. This is almost always the case when $\{Y_{1} , \ldots , Y_{n} \}$ is a sequence of output observations from within a single replication. In this situation,$Y_{1} , Y_{2} , \ldots$. is an autocorrelated sequence, sometimes called a time series.

Suppose that our point estimator for θ is the sample mean $\textstyle {\bar {Y}} = \sum_{i = 1} ^ {n} Y_{i} / n . \operatorname{A}$ general result^2 from mathematical statistics is that the variance of $\bar {Y}$ is

where cov $( Y_{i} , Y_{i} ) = \operatorname{V} ( Y_{i} )$. To construct a confidence interval for θ, an estimate of $\mathrm{V} ( {\bar {Y}} )$ is required. But obtaining an estimate of Equation (22) is futile, because each term cov $( Y_{i} , Y_{j} )$ could be different, in general. Fortunately, systems that have a steady state will, if simulated long enough to pass the transient phase, produce an output process that is approximately covariance stationary. Intuitively, stationarity implies that $Y_{i + k}$ depends on $Y_{i + 1}$ in the same manner as $Y_{k}$ depends on $Y_{1}$. In particular, the covariance between two random variables in the time series depends only on the number o observations between them, called the $l a g$ For a covariance-stationary time series,$\{Y_{1} , Y_{2} , \ldots \}$, define the lag-k autocovariance by

which, by definition of covariance stationarity, is not a function of i. For $k = 0 ,$, γ becomes the population variance $\sigma^{2} ;$ that is,

The lag-k autocorrelation is the correlation between any two observations k apart. It is defined by

and has the property that

$$

• 1 \leq \rho_{k} \leq 1, \quad k = 1, 2, \dots $$

If a time series is covariance stationary, then Equation (22) can be simplified substantially. Algebra shows that

where $\rho_{k}$ is the lag-k autocorrelation given by Equation (25).

When $\rho_{k} > 0$ for all k (or most k), the time series is said to be positively autocorrelated. In this case, large observations tend to be followed by large observations, small observations by smal ones. Such a series will tend to drift slowly above and then below its mean. Figure 8(a) is an example of a stationary time series exhibiting positive autocorrelation. The output data from most queueing simulations are positively autocorrelated.

(a)

(b)

(c) Figure 8 (a) Stationary time series $\pmb {Y_{j}}$ exhibiting positive autocorrelation; (b) stationary time series $\pmb {Y_{j}}$ exhibiting negative autocorrelation; (c) nonstationary time series with an upward trend.

On the other hand, if some of the $\rho_{k} < 0 ,$ the series $Y_{1} , Y_{2} , \ldots$. will display the characteristics of negative autocorrelation. In this case, large observations tend to be followed by small observations, and vice versa. Figure 8(b) is an example of a stationary time series exhibiting negative autocorrelation. The output of certain inventory simulations might be negatively autocorrelated.

Figure 8(c) also shows an example of a time series with an upward trend. Such a time series is not stationary; the probability distribution of $Y_{i}$ is changing with the index i.

Why does autocorrelation make it difficult to estimate ${\mathrm{V}} ( {\bar {Y}} ) {\mathrm{?}}$ Recall that the standard estimator for the variance of a sample mean is $S^{2} / n$. By using Equation (26), it can be shown [Law, 1977] that the expected value of the variance estimator $S^{2} / n$ is

where

and c is the quantity in brackets in Equation (26). The effect of the autocorrelation on the estimator $S^{2} / n$ is derived by an examination of Equations (26) and (28). There are essentially three possibilities:

Case 1

If the $Y_{i}$ are independent, then $\rho_{k} \ = \ 0$ for $k = 1 , 2 , 3 , \ldots$. Therefore,$\begin{array} {r} {c = 1 + 2 \sum_{k = 1} ^ {n - 1} ( 1 -} \end{array}$$k / n ) \rho_{k} = 1$ and Equation (26) reduces to the familiar $\sigma^{2} / n$. Notice also that $B = 1$, so $S^{2} / n$ is an unbiased estimator of V(Y¯ ). The $Y_{i}$ will always be independent when they are obtained from differen replications; that independence is the primary reason that we prefer experiment designs calling for multiple replications.

Case 2

If the autocorrelations $\rho_{k}$ are primarily positive, then $\begin{array} {r} {c = 1 + 2 \sum_{k = 1} ^ {n - 1} ( 1 - k / n ) \rho_{k} > 1} \end{array}$, so that $n / c < n ,$, and hence $B < 1$. Therefore,$S^{2} / n$ is biased low as an estimator of $\mathrm{V} ( {\bar {Y}} )$. If this correlation were ignored, the nominal $100 ( 1 - \alpha ) \%$ confidence interval given by Expression (10) would be too short, and its true confidence coefficient would be less than $1 - \alpha .$. The practical effect would be that the simulation analyst would have unjustified confidence in the apparent precision of the poin estimator due to the shortness of the confidence interval. If the correlations $\rho_{k}$ are large, B could be quite small, implying a significant underestimation.

Case 3

If the autocorrelations $\rho_{k}$ are substantially negative, then $0 \leq c < 1$, and it follows that $B > 1$ and $S^{2} / n$ is biased high for $\mathrm{V} ( {\bar {Y}} )$. In other words, the true precision of the point estimator $\bar {Y}$ would be greater than what is indicated by its variance estimator $S^{2} / n$, because

As a result, the nominal $100 ( 1 - \alpha ) \%$ confidence interval of Expression (10) would have true confidence coefficient greater than $1 - \alpha$. This error is less serious than Case 2, because we are unlikely to make incorrect decisions if our estimate is actually more precise than we think it is.

A simple example demonstrates why we are especially concerned about positive correlation: Suppose you want to know how students on a university campus will vote in an upcoming election. To estimate their preferences, you plan to solicit 100 responses. The standard experiment is to randomly select 100 students to poll; call this experiment A. An alternative is to randomly select 20 students and ask each of them to state their preference 5 times in the same day; call this experiment B. Both experiments obtain 100 responses, but clearly an estimate based on experiment B will be less precise (will have larger variance) than an estimate based on experiment A. Experiment A obtains 100 independent responses, whereas experiment B obtains only 20 independent responses and 80 dependent ones. The five opinions from any one student are perfectly positively correlated (assuming a student names the same candidate all five times). Although this is an extreme example, i illustrates that estimates based on positively correlated data are more variable than estimates based on independent data. Therefore, a confidence interval or other measure of error should account correctly for dependent data, but $S^{2} / n$ does not.

Two methods for eliminating or reducing the deleterious effects of autocorrelation upon estimation of a mean are given in the following sections. Unfortunately, some simulation languages either use or facilitate the use of $S^{2} / n$ as an estimator of $\mathrm{V} ( {\bar {Y}} )$), the variance of the sample mean, in all situations. If used uncritically in a simulation with positively autocorrelated output data, the downward bias in $S^{2} / n$ and the resulting shortness of a confidence interval for θ will convey the impression of much greater precision than actually exists. When such positive autocorrelation is present in the output data, the true variance of the point estimator $\bar {Y}$ can be many times greater than is indicated by $S^{2} / n$

5.3 Replication Method for Steady-State Simulations

If initialization bias in the point estimator has been reduced to a negligible level (through some combination of intelligent initialization and deletion), then the method of independent replications can be used to estimate point-estimator variability and to construct a confidence interval. The basic idea is simple: Make R replications, initializing and deleting from each one the same way.

If, however, significant bias remains in the point estimator and a large number of replications are used to reduce point-estimator variability, the resulting confidence interval can be misleading. This happens because bias is not affected by the number of replications R; it is affected only by deleting more data (i.e., increasing $T_{0} )$ or extending the length of each run (i.e., increasing $T_{E} )$ Thus, increasing the number of replications R could produce shorter confidence intervals around the “wrong point.” Therefore, it is important to do a thorough job of investigating the initial-condition bias.

If the simulation analyst decides to delete d observations of the total of n observations in a replication, then the point estimator of θ is

That is, the point estimator is the average of the remaining data. The basic raw output data $\{Y_{r j} , r =$$1 , \dots , R ; j = 1 , \dots , n \}$ are exhibited in Table 4. For instance,$Y_{r j}$ could be the delay of customer j in queue, or the response time of job j in a job shop, on replication r. The number d of deleted observations and the total number of observations n might vary from one replication to the next, in which case replace d by $d_{r}$ and n by $n_{r} .$. For simplicity, assume that d and n are constant over replications.

Table 4 Raw Output Data from a Steady-State Simulation

When using the replication method, each replication is regarded as a single sample for the purpose of estimating θ. For replication $r ,$ define

as the sample mean of all (nondeleted) observations in replication r. Because all replications use different random-number streams and all are initialized at time 0 by the same set of initial conditions $I_{0} ,$, the replication averages

are independent and identically distributed random variables; that is, they constitute a random sample from some underlying population having unknown mean

The overall point estimator, given in Equation (29), is also given by

as can be seen from Table 4 or from using Equation (21). Thus, it follows that

also. If d and n are chosen sufficiently large, then $\theta_{n , d} \approx \theta$, and ${\bar {Y}} _ {\cdot \cdot} ( n , d )$ is an approximately unbiased estimator of θ. The bias in ${\bar {Y}} _ {\cdot \cdot} ( n , d )$ is $\theta_{n , d} - \theta$ For convenience, when the value of n and d are understood, abbreviate $\bar {Y} _ {r} . ( n , d )$(the mean of the undeleted observations from the rth replication) and ${\bar {Y}} _ {\cdot \cdot} ( n , d )$[the mean of $\bar {Y} _ {1 \cdot} ( n , d ) , \dots , \bar {Y} _ {R \cdot} ( n , d ) ]$ by ${\bar {Y}} _ {r} .$ and ${\bar {Y}} .$, respectively. To estimate the standard error of ${\bar {Y}} . .$, first compute the sample variance,

The standard error of ${\bar {Y}} .$ is given by

A $100 ( 1 - \alpha ) \%$ confidence interval for $\theta ,$, based on the t distribution, is given by

where $t_{\alpha / 2 , R - 1}$ is the $100 ( 1 - \alpha / 2 )$ percentage point of a t distribution with $R - 1$ degrees of freedom. This confidence interval is valid only if the bias of ${\bar {Y}} .$ is approximately zero.

As a rough rule, the length of each replication, beyond the deletion point, should be at least ten times the amount of data deleted. In other words,$( n - d )$ should at least 10d (or more generally,$T_{E}$ should be at least $10 T_{0} )$. Given this run length, the number of replications should be as many as time permits, up to about 25 replications. Kelton [1986] established that there is little value in dividing the available time into more than 25 replications, so, if time permits making more than 25 replications of length $T_{0} + 10 T_{0}$, then make 25 replications of longer than $T_{0} + 10 T_{0} \mathrm{.}$, instead. Again, these are rough rules that need not be followed slavishly.