**CHAPTER 9 **

**Output Data Analysis for a Single System **

9.1 Introduction ...2

9.2 Transient and Steady-State Behavior of a Stochastic Process ...10

9.3 Types of Simulations with Regard to Output Analysis ...12

9.4 Statistical Analysis for Terminating Simulations...14

9.4.1 Estimating Means ...16

9.4.2 Estimating Other Measures of Performance...22

9.4.3 Choosing Initial Conditions ...24

9.5 Statistical Analysis for Steady-State Parameters ...25

9.5.1 The Problem of the Initial Transient...26

9.5.2 Replication/Deletion Approaches for Means ...27

9.5.3 Other Approaches for Means ...28

9.5.4 Estimating Other Measures of Performance...40

9.6 Statistical Analysis for Steady-State Cycle Parameters...41

9.7 Multiple Measures of Performance ...43

9.8 Time Plots of Important Variables...45

**9.1 Introduction **

Basic, most serious disadvantage of simulation:

With a stochastic simulation, don’t get exact “answers” from run Two different runs of same model ⇒ different numerical results Random input:

Random numbers Random variates

### →

Simulationmodel/code

### →

Random output:Performance measures Thus, the output performance measures:

Really observations from their probability distribution Ask questions about this distribution:

*Mean (expected value): E(average WIP) *
Variance: Var(average WIP)
Probabilities: *P(average WIP > 250) *

Quantiles: *What value of x is such that P(avg. WIP > x) *≤ 0.02?

*Interpreting simulation output: statistical analysis of output data *

Failure to recognize, deal with randomness in simulation output can lead to serious errors, misinterpretation, bad decisions

Also, must take care to use appropriate statistical methods, since simulation output data are usually nonstationary, autocorrelated, and non-normal, contrary to assumptions behind classical IID statistical methods

Enhanced computer power and speed is making it much easier to carry out

**The Statistical Nature of Simulation Output **

*Let Y*_{1}*, Y*_{2}*, ... be an output process from a single simulation run *
*Y**i** = delay in queue of ith arriving customer *

*Y*_{i}* = production in ith hour in a factory *

*Y** _{i}*’s are random variables that are generally neither independent nor identically
distributed (nor normally distributed), so classical IID normal-theory statistical

*methods don’t apply directly to the Y*

*’s*

_{i}*Let y*11*, y*12*, ..., y**1m** be a realization of the random variables Y*1*, Y*2*, ..., Y**m* resulting
*from making a single simulation run of length m observations, using a particular *
stream of underlying U(0, 1) random numbers.

If we use a separate stream of random numbers for another simulation run of this
*same length, we get a realization y*21*, y*22*, ..., y**2m* that is independent of, but
*identically distributed to, y*_{11}*, y*_{12}*, ..., y*_{1m}

*Make n such independent runs, each using “fresh” random numbers, to get *
*y*11*, y*12*, ..., y**1m *

*y*_{21}*, y*_{22}*, ..., y*_{2m}

**. **
**. **
**. **

*y*_{n1}*, y*_{n2}*, ..., y*_{nm}

Within a row: not IID
*Across the ith column: IID *

*realizations of the r.v. Y** _{i}* (but still not
necessarly normally distributed)

Can compute a summary measure within a run, and then the summary measures across the runs are IID (but still not necessarily normally distributed)

Bank with 5 tellers, one FIFO queue, open 9am-5pm, flush out before stopping Interarrivals ~ expo (mean = 1 min.), service times ~ expo (mean = 4 min.) Summary measures from 10 runs (replications):

Clearly, there’s variation across runs;

need appropriate statistical analysis

**Types of Output Performance Measures **
What do you want to know about the system?

Average time in system

Worst (longest) time in system Average, worst time in queue(s)

Average, worst, best number of “good” pieces produced per day

Variability (standard deviation, range) of number of parts produced per day Average, maximum number of parts in the system (WIP)

Average, maximum length of queue(s)

Proportion of time a machine is down, up and busy, up and idle Ask the same questions of the model/code

Think ahead: Asking for additional output performance measures can change how the simulation code is written, or even how the system is modeled

Simple queueing model:

Server

Customers in queue

Customer in service Customer

arrivals Customer

departures

Want: Average number of customers in queue Proportion of time server is busy

Maybe: Average time customers spend in queue

Question: How does wanting the average time in queue affect how the model is coded, data structures, etc.?

*As simulation progresses through time, there are typically two kinds of processes *
that are observed:

*Discrete-time process: There is a natural “first” observation, “second” *

observation, etc.—but can only observe them when they “happen”

*S*_{i}* = time in system of ith part produced, i *∈ {1, 2, ...}

*Suppose there are M parts produced during the simulation *

*i*

1 2 3 *...* *M *

* *
*S**i*

Typical discrete-time output performance measures:

Average time in system:

*M*
*S*
*M*

*S*

*M*

*i*

### ∑

*i*

= =^{1}

) (

Maximum time in system: _{i}

*M*

*i* *S*

*M*
*S*

,..., 2 , 1

max )

(

* = =

Proportion of parts that were in system more than 60 minutes:

*M*
*S*
*I*

*M*
*P*

*M*

*i*

### ∑

*i*

= ∞

= ^{1} ^{(}^{60}^{,} ^{)}

60

) ( )

( , where

≤

= >

∞ 0 if 60

60 if

) 1

)(

, 60 (

*i*
*i*

*i* *S*

*S* *S*
*I*

Other examples of discrete-time processes:

*D*_{i}* = delay of ith customer in queue *

*Y*_{i}* = throughput (production) during ith hour *
*B*_{i}* = 1 if caller i gets a busy signal, 0 otherwise *

*Continuous-time process: Can jump into system at any point in time (real, *
continuous time) and take a “snapshot” of something — there is no natural

“first” or “second” observation

*Q(t) = number of parts in a particular queue at time t *∈ [0, ∞)
*Run simulation for T units of simulated time *

* Q(t )*

0 1 2 3

*t* *T*

Typical continuous-time output performance measures:

Time-average length of queue:

*T*
*dt*
*t*
*Q*
*T*

*Q*

### ∫

*T*

= ^{0} ( )
)

(

Maximum length of queue: *( ) max ( )

0 *Q* *t*

*T*
*Q*

*T*
*t*≤

= ≤

Proportion of time that there were more than two in the queue:

*T*

*dt*
*t*
*Q*
*I*

*T*
*P*

### ∫

*T*

^{∞}

= ^{0} ^{(}^{2}^{,} ^{)}

2

)) ( ( )

(

*Another important kind of continuous-time statistic: utilizations *

Let =

*t*
*t* *t*

*B* 0 if server isidle at time
at time
busy

is server if

) 1 (

0 1

*t* *T*

* B(t )*

Server utilization (proportion of time busy):

*T*
*dt*
*t*
*B*
*T*

*U*

### ∫

*T*

= ^{0} ( )
)

(

Other examples of continuous-time processes:

*N(t) = number of parts in shop at time t (WIP) *

*D(t) = number of machines down at time t *

Typically, we want to observe several (maybe lots of) different performance measures from the same system/model

Usually low additional cost/hassle to do so, can always ignore later
*But not getting a particular output measure could imply rerunning *

Difficulty in statistical analysis of output with several performance measures:

May want to make several simultaneous estimates, statements Be careful how this is done, what is said

*Multiple-comparisons problem in statistics literature (more in Sec. 9.7) *

**9.2 Transient and Steady-State Behavior of a ** **Stochastic Process **

*Output process (discrete-time) Y*1*, Y*2, ...

*Let F*_{i}*(y | I) = P(Y** _{i}*≤

*y | I) be the transient (cumulative) distribution of the process*

*at (discrete) time i*

*In general, F*_{i}* depends on both the time i and the initial condition I *
Corresponding transient density functions:

*If there is a CDF F(y) such that F*_{i}*(y | I) *→* F(y) as i *→ ∞* for all y and for all initial *
*conditions I, then F(y) is the steady-state distribution of the output process *
*F(y) may or may not exist *

*F(y) must be independent of the initial conditions — same for all I *

*Roughly speaking, if there is a time index k such that for i > k F*_{i}*(y | I) *≈* F(y) in *
*some sense, then we say that the process is “in steady state” after time k *

*Even though the distribution of the Y**i**’s after time k is not appreciably changing, *
*observations on the Y** _{i}*’s could still have large variance and thus “bounce
around” a lot — they’re just not systematically trending any more

Steady-state distribution does not depend on initial conditions, but the nature and rate of convergence of the transient distributions can depend heavily on the initial conditions

*M/M/1 queue, E(delay in queue), different number of customers s present initially: *

*Inventory system, E(cost in month i): *

**9.3 Types of Simulations with Regard to ** **Output Analysis **

*Terminating: Parameters to be estimated are defined relative to specific initial and *
stopping conditions that are part of the model

There is a “natural” and realistic way to model both the initial and stopping conditions

Output performance measures generally depend on both the initial and stopping conditions

*Nonterminating: There is no natural and realistic event that terminates the model *

Interested in “long-run” behavior characteristic of “normal” operation If the performance measure of interest is a characteristic of a steady-state

*distribution of the process, it is a steady-state parameter of the model *
Theoretically, does not depend on initial conditions

Practically, must ensure that run is long enough so that initial-condition effects have dissipated

Not all nonterminating systems are steady-state: there could be a periodic

*“cycle” in the long run, giving rise to steady-state cycle parameters *

Focus on terminating vs. steady-state simulations:

Which is appropriate?

Depends on goals of the study

*Statistical analysis for terminating simulations is a lot easier *
Is steady-state relevant at all? Maybe:

24 hours/day, “lights-out” manufacturing Global telecommunications

Design conservatively for peak load of infinite duration Some examples:

Physical model

Terminating estimand Steady-state estimand Single-server

queue

Expected average delay in queue of first 25 customers, given empty-and-idle initial conditions

Long-run expected delay in queue of a customer

Manufacturing system

Expected daily production, given some number of

workpieces in process initially

Expected long-run daily production

Reliability system

Expected lifetime, or

probability that it lasts at least a year, given all components initially new, working

Probably not sensible

Battlefield model

Probability that attacking force loses half its strength before defending force loses half its strength

Probably not sensible

**9.4 Statistical Analysis for Terminating ** **Simulations **

*Make n IID replications of a terminating simulation *
Same initial conditions for each replication
Same terminating event for each replication
Separate random numbers for each replication

*Let X*_{j}* be a summary measure of interest from the jth replication *

*e.g., X**j** = the average delay in queue of all customers in the jth replication *

*Then X*_{1}*, X*_{2}*, ..., X** _{n}* are IID random variables, can apply classical statistical analysis
to them

Rely on central limit theorem to justify normality assumption even though it’s generally not true

*So basic strategy is replication of the whole simulation some number n of times *
*One simulation run is a sample of size one (not worth much statistically) *

**What About Classical Statistics? **

Classical statistical methods don’t work directly within a simulation run, due to
*autocorrelation usually present *

*Example: Delays in a queue of individual jobs: D*_{1}*, D*_{2}*, D*_{3}*, ..., D** _{m}*
Want to estimate

*µ = E(average delay of the m jobs)*

Sample mean *D* *m* *D* *m*

*m*

*i*

### ∑

*i*

=

=

1

)

( is an unbiased estimator for *µ *

Need to estimate Var(*D(m*)) for confidence intervals on *µ, test hypotheses like H*0:
*µ = µ*_{0}

But “sample variance”

### (

^{(}

^{)}

### ) [

^{(}

^{1}

^{)}

### ]

1

2 −

### ∑

−=

*m*
*m*
*m*

*D*
*D*

*m*

*i*

*i* may be severely biased for

Var(*D(m*))
Reason:

*Corr(D*_{i}*, D** _{i + l}*) ≠ 0, in general

Unbiasedness of variance estimators follows from independence of data, which
*is not true within a simulation *

Usual situation:

*Positive autocorrelation: Corr(D*_{i}*, D** _{i + l}*) > 0

Causes

### ( ) [ ]

### ∑

− −=

) 1 ( )

(

1

2 *m* *m*

*m*
*D*
*D*
*E*

*m*

*i*

*i* < Var(*D(m*)*) — maybe far too small *
Intuition:

*D*_{i + 1}* is pretty much the same as D*_{i}

*D** _{i}*’s are more stable than if they were independent
Their sample variance is understated

Thus, must take care to estimate variances properly: understating variances Have too much faith in our point estimates

Believe our simulation results too much

**9.4.1 Estimating Means **

Want: Estimate of some parameter *µ of the process *
Often (not always): *µ = E(something) *

*µ = E(average delay in queue of m customers) *
*µ = E(time-average WIP) *

*µ = E(machine utilization) *

*µ = P(average delay > 5 minutes) = E[I*_{(5,∞)}(average delay)]

Point estimate: *µˆ* = 12.3

How close is *µˆ* to the true unknown value of *µ? *

*Customary, useful method for assessing precision of estimator: confidence interval *
for *µ *

*Pick confidence level 1 – α (typically 0.90, 0.95, etc.) *

*Use simulation output to construct an interval [A, B] that covers µ with *
probability 1 – *α *

Interpretation: 100(1 – *α)% of the time, the interval formed in this way will *
cover *µ *

*µ (unknown)*

*Wrong interpretation: “I’m 95% sure that µ is between 9.4 and 11.1” *

Common approach to statistical analysis of simulation output:

*Can’t do “classical” (IID, unbiased) statistics within a simulation run *
Try to modify setup, design, to get back to classical statistics

In terminating simulations, this is conceptually easy:

*Make n independent replications of the whole simulation *
*Let X*_{j}* be the performance measure from the jth replication *

*X** _{j}* = average of the delays in queue

*X*

*= time-average WIP*

_{j}*X** _{j}* = utilization of a bottleneck machine

*Then X*_{1}*, X*_{2}*, ..., X** _{n}* are IID and unbiased for

*µ = E(X*

*j*)

*Apply classical statistics to X*_{j}*’s, not to observations within a run *
Approximate 100(1 – *α)% confidence interval for µ: *

*n*
*X*
*n*

*X*

*n*

*j*

### ∑

*j*

= =^{1}

)

( is an unbiased estimator of *µ *

1 )) ( (

) (

2 2

−

−

=

### ∑

*n*

*n*
*X*
*X*
*n*

*S*

*j*

* is an unbiased estimator of Var(X** _{j}*)

*n*
*n*
*t* *S*

*n*

*X* * _{n}* ( )

)

( ± _{−}_{1}_{,}_{1}_{−}_{α}_{/}_{2} covers *µ with approximate probability 1 – α *

*(t*_{n–1,1–}* _{α/2}* = point below which is area (probability) 1–

*α/2 in Student’s t*

*distribution with n – 1 d.f.)*

Most important point:

The “basic ingredients” to the statistical analysis are the performance measures from the different, independent replications

*One whole simulation run = a “sample” of size one (not worth much) *

*Example: n = 10 replications of single-server queue *
*X*_{j}* = average delay in queue from jth replication *

*X** _{j}*’s: 2.02, 0.73, 3.20, 6.23, 1.76, 0.47, 3.89, 5.45, 1.44, 1.23
Want 90% confidence interval, i.e.,

*α = 0.10*

) 10 (

*X* * = 2.64, S*^{2}*(10) = 3.96, t*_{9, 0.95} = 1.833

Approximate 90% confidence interval is 2.64 ± 1.15, or [1.49, 3.79]

Other examples in text:

Inventory model

Estimation of proportions

**Why “approximate” 90% confidence interval? **

*Assumes X** _{j}*’s are normally distributed — never true, but does it matter?

*Central-limit theorem: *

*As n (number of replications) grows, coverage probability *→ 1 – *α *
*Robustness study with this model: *

Know *µ = 2.12 from queueing theory *
*Pick an n (like n = 10) *

Make 500 “90%” confidence intervals (total no. runs = 10 × 500)
Count % of these 500 intervals covering *µ = 2.12: *

*n * Estimated coverage (want 90%)
5

10 20 40

88%

86%

89%

92%

Bad news: actual coverages can depend (a lot) on the model
Reliability model: ^{1}

2

3

Components fail independently

*Times T** _{i}* to failure in components ~ Weibull(

*α = 0.5, β = 1.0)*

*Time to system failure = T = min(T*

_{1}

*, max(T*

_{2}

*, T*

_{3}))

*µ = E(time to system failure) = 0.78 *

*n * Estimated coverage (want 90%)
5

10 20 40

71%

75%

80%

84%

**Obtaining a Specified Precision **

*If the number n of replication is chosen too small, the confidence intervals might *
too wide to be useful

*M/M/1 example: *

*90% confidence interval from n = 10 replications: 2.64 ± 1.15, or [1.49, 3.79] *

Half width (1.15) is 44% of point estimate (2.64)
*Equivalently: 2.64 ± 44%, not very precise *

Half-width expression: *δ(α, n) = *

*n*
*n*
*t*_{n}*S*( )

2
/
1
,
1 −*α*

−

To decrease: *α * ↑: undesirable, since *α = probability of missing *
*S(n) *↓: estimates Var(*X* * _{j}*), which is fixed (maybe ...)

*n*↑

*: more replications*

*Sequential sampling: Increase n until δ(α, n) is “small enough” *

Two versions of what “small enough” might mean (more details in text):

*Absolute precision: *

Specify *β > 0, want n big enough so that δ(α, n) < β *

Requires at least some knowledge of context to set meaningful *β *
*Relative precision: *

Specify *γ (0 < γ < 1), want n big enough so that δ*(*α*,*n*) *X*(*n*)<*γ*

Need not know much about context to set meaningful *γ *
Notes:

Above description leaves out a few technical details; see text

“Fixes” robustness issue: As *β or γ →* 0, coverage probability → 1 – *α *

Can be dangerous for small *β or γ : Required n increases quadratically as β or *
*γ decrease *

May be difficult to automate with a simulation language, depending on modeling constructs available, what automated statistical capabilities present, and what access the user has to internal software variables

**9.4.2 Estimating Other Measures of Performance **

*Sometimes can miss important system/model characteristics if we look only at *
averages

Other measures: Proportions, variances, quantiles

**Proportions **

Compare two operating policies for queueing system with five servers

vs.

Estimates

Performance measure Five queues Single queue Average delay in queue

Average number in queue(s) Number of delays ≥ 20 minutes

5.57 minutes 5.52

33

5.57 minutes 5.52

6

**Variances (or Standard Deviations) **

*Interested in process variability *
*X** _{j}* = daily production of good items
Want to estimate Var(

*X*

*)*

_{j}*Make n replications, compute S(n) as before *

Confidence interval on Var(*X* * _{j}*): use chi-square distribution

**Quantiles **

*Inventory system, X** _{j}* = maximum inventory level during the horizon
Want to determine storage capacity that is sufficient with probability 0.98

*Want to find x such that P(X*

*≤*

_{j}*x) = 0.98*

One approach (more details in text):

*Make n replications, observe X*_{1}*, X*_{2}*, ..., X*_{n}*Sort the X** _{j}*’s into increasing order

*Estimate x to be a value below which are 98% of the X** _{j}*’s

**9.4.3 Choosing Initial Conditions **

For terminating simulations, the initial conditions can affect the output performance measure, so the simulations should be initialized appropriately

Example: Want to estimate expected average delay in queue of bank customers who arrive and complete their delay between noon and 1:00pm

Bank is likely to be crowded already at noon, so starting empty and idle at noon will probably bias the results low

Two possible remedies:

If bank actually opens at 9:00am, start the simulation empty and idle, let it run for 3 simulated hours, clear the statistical accumulators, and observe statistics for the next simulated hour

Take data in the field on number of customers present at noon, fit a (discrete) distribution to it, and draw from this distribution to initialize the simulation at time 0 = noon. Draw independently from this distribution to initialize multiple replications.

Note: This could be difficult in simulation software, depending on the modeling constructs available

**9.5 Statistical Analysis for Steady-State ** **Parameters **

Much more difficult problem than analysis for terminating simulations Want to estimate (e.g. for discrete- and continuous-time processes)

=

= =

∞

→

∞

→

*t*
*t*

*Q*
*t*
*Q*
*E*

*i*
*D*

*D*
*E*

*t*

*i*
*i* *i*

at time queue

in number )

( )) ( ( lim

customer th

of queue in

delay )

(
*ν* lim

Basic question for designing runs:

Many short

runs One long

run
*X*_{1}

*X*_{2}

*X*3

*X*_{4}

*X*_{5}

*X*_{1}

Many short runs One long run

Good Simple (same as terminating) Get IID data

*Less point-estimator bias *
No restarts

**9.5.1 The Problem of the Initial Transient **

If steady-state is the goal, initial conditions will generally bias the results of the simulation for some initial period of time

*Most common technique is to warm up the model, also called initial-data deletion *
*Identify index l (for discrete-time processes) or time t*_{0} (for continuous-time

processes) beyond which the output appears not to be drifting any more Clear statistical accumulators at that time

Start over with data collection, and “count” only data observed past that point After warmup period, observations will still have variance, so will bounce around

— they are just not “trending” any more

Facility for doing this in most simulation software (but the user must specify the warmup period)

Challenge — identifying a warmup period that is long enough, yet no so long as to be excessively wasteful of data — see text for details and examples

Some statistical-analysis tests have been devised

Most practical (and widely-used) method is to make plots, perhaps averaged across and within replications to dampen the noise, and “eyeball” a cutoff If there are multiple output processes, and if they disagree about what the

appropriate warmup period is, a decision must be made whether to use

different warmups for each process, or to use the same one for each process

— which would have to be the maximum of the individual warmups, to be save, and so would be conservative for most of the output processes

A different approach: Try to find “smarter” initialization states or distributions that

**9.5.2 Replication/Deletion Approaches for Means **

Assume that an appropriate warmup period has been determined

*X*_{j}* = output measure on jth replication, collected only past warmup point *
Proceed with statistical analysis exactly as in terminating case

Make independent replications, each warmed up

Compute mean, variance estimates across replications, confidence intervals Advantages (compared to methods to be discussed below):

Simple — aside from warmup, the same as for terminating simulations

Get truly IID observations — important not only for variance estimation and confidence-interval construction, but also for more sophisticated statistical goals and techniques to be discussed in Chap. 10–12

Disadvantages:

No completely reliable method to identify an appropriate warmup Too long ⇒ wasteful of data

Too short ⇒ point-estimator bias, which can have serious consequences, especially if used in concert with a sequential procedure:

)
*(n*

*X* is no longer an unbiased estimator of *ν *

*Confidence interval is centered in the –wrong” place— at E[X(n*)] ≠ *ν *
*As n *↑, confidence interval shrinks down around the wrong point, causing

*coverage to drop *

* n*=10
* n*=20
* n*=30
* n*=40

Work harder,

do worse (in coverage sense)

**9.5.3 Other Approaches for Means **

Make just one “replication,” presumably to ameliorate initial bias

Point estimator of *ν: average Y(m*)* of all the data Y*_{1}*, Y*_{2}*, ..., Y*_{m}* in the run *
Problem: How to estimate Var(*Y(m*)), needed to get c.i.’s, etc.?

*Know one way not to do this: *

### (

^{(}

^{)}

### ) [

^{(}

^{1}

^{)}

### ]

1

2 −

### ∑

−=

*m*
*m*
*m*

*Y*
*Y*

*m*

*i*
*i*

Several methods to estimate Var(*Y(m*)):

Batch means

Time-series models Spectral analysis

Standardized time series

A different one-long-run approach (different point estimator):

Regenerative method

Two alternative modes of operation:

*Fixed-sample size procedures — select run length m in advance, precision not *
controlled

*Sequential procedures — prespecify precision, increase run length m as needed *

**Batch Means **

*Divide the run (of length m) into n contiguous “batches” of length k each (m = nk) *
Let *Y*_{j}*(k*)* be the (batch) mean of the jth batch of size k *

*i*

* *
*Y*_{i}

*k* *k* *k* *k* *k*

*Y *_{1} *Y *_{ } * *_{2} *Y *_{ } * *_{3} *Y *_{ } * *_{4} *Y *_{ } * *_{5} *m* * * =*nk*
Pretend: *Y*_{j}*(k*)’s are IID, unbiased for *ν *

Note: ( )

) (

1 1

*m*
*m* *Y*

*Y*
*n*

*k*
*Y*

*m*

*i*
*i*
*n*

*j*
*j*

=

=

### ∑

### ∑

= =, i.e., mean of batch means is the “grand” mean

Form “sample variance” among the batch means,

### ( )

1 ) ( ) ( )

( ^{1}

2

2

−

−

=

### ∑

=

*n*

*m*
*Y*
*k*
*Y*
*n*

*S*

*n*

*j*
*j*
*Y*

Approximate 100(1 – *α)% confidence interval for ν: *

*n*
*n*
*t* *S*

*m*

*Y* _{n}* ^{Y}* ( )

)

( ± _{−}_{1}_{,}_{1}_{−}_{α}_{/}_{2}

*(k) * *(k) * *(k) * *(k) * *(k) *

Appeals of batch means:

Simple (relatively)

Often works fairly well (in terms of coverage)

Automatically implemented in some simulation software

Issues with batch means (see text for more details and references):

Choose batches big enough so that*Y*_{j}*(k*)’s are approximately uncorrelated
Otherwise, *S*_{Y}^{2}(*n*) can be biased (usually low) for Var(*Y*_{j}*(k*)), causing

undercoverage

*How to choose batch size k? Equivalently, how many batches n? *

*Evidence: It may never pay to have more than n = 20 or 30 batches *

Reason: Due to autocorrelation, splitting run into a larger number of smaller batches, while increasing degrees of freedom, degrades the quality

(variability) of each individual batch Generalizations of batch means

Overlapping batch means Separated batches

Weighting points within the batches

Sequential batch-means methods to control confidence-interval width
*Fix the number n of batches *

*Increase overall run length m *

*Increase batch size k proportionally to increase in m *

**Time-Series Models **

*Assume a relatively simple statistical model for the output process *
*Discrete-time output process Y*_{1}*, Y*_{2}, ... with steady-state mean *ν *

As with batch means, use overall average *Y(m*) as point estimator of *ν *
Examples of statistical models for output process:

*Autoregressive of order 2 (AR(2)): *

*Y** _{i}* =

*ν + φ*1

*(Y*

*–*

_{i–1}*ν) + φ*2

*(Y*

*–*

_{i–2}*ν) + ε*

*i*,

*ε*

*i*’s IID N(0,

*σ*

^{2})

*Autoregressive moving average of order (2, 1) (ARMA(2, 1)):*

*Y** _{i}* =

*ν + φ*1

*(Y*

*–*

_{i–1}*ν) + φ*2

*(Y*

*–*

_{i–2}*ν) + θ*1

*ε*

*i–1*+

*ε*

*i*,

*ε*

*i*’s IID N(0,

*σ*

^{2})

*In general, for ARMA(p, q) model:*

*Use simulation output to identify (estimate) p and q — get * *pˆ* and *qˆ* (several
methods exist)

Estimate parameters via regression (least-squares) — (“fit the model”):

2

ˆ 2 1

ˆ 2 1

ˆ

,..., ˆ , ˆ ˆ

,..., ˆ , ˆ ˆ

*σ*

*θ*
*θ*

*θ*

*φ*
*φ*

*φ*

*q*
*p*

Under this model, Var(*Y(m*)*) = g(φ*1, *φ*2, ..., *φ**p*, *θ*1, *θ*2, ..., *θ**q*, *σ*^{2}) (the
*function g is known but messy) *

Estimate Var(*Y(m*)*) by g(φ*ˆ_{1},*φ*ˆ_{2},...,*φ*ˆ_{p}_{ˆ},*θ*ˆ_{1},*θ*ˆ_{2},...,*θ*ˆ_{q}_{ˆ},*σ*ˆ^{2})
Appeals of time-series models:

Physical intuition

Some theoretical support — most processes can be approximated by an
*ARMA(p, q) if we’re willing to admit large p and q *

Issues with time-series models:

**Spectral Analysis **

Assume that output process is covariance-stationary:

*Cov(Y*_{i}*, Y*_{i + j}*) = C*_{j}*, i.e., it depends only on j, and not on i *
Fairly mild assumption, especially after some warmup

Then

*m*

*C*
*m*
*j*
*C*

*m*
*Y*

*m*

*j*

### ∑

^{−}

*j*

=

− +

=

1

1

0 (1 / )

)) ( ( Var

*Estimate C*_{j}*, j = 1, 2, ..., m – 1 from output data: *

### ( ) ( )

*j*
*m*

*m*
*Y*
*Y*

*m*
*Y*
*Y*
*C*

*j*
*m*

*i*

*j*
*i*
*i*

*j* −

−

−

=

### ∑

^{−}

= +

1

) ( )

( ˆ

Plug *Cˆ ’s into formula for Var(*_{j}*Y(m*)*) in place of corresponding C** _{j}*’s
Appeals of spectral analysis:

Using an exact variance formula

Relationship to other methods (overlapping batch means) Issues with spectral analysis:

*For large m (as we’d expect), computationally burdensome to get all the* *Cˆ ’s — *_{j}*number of operations is about m*^{2}/2

*Remedy: Computational device — Fast Fourier Transform *

*For j near the end (m – 1), Cˆ will be based on only a few terms, so will itself be *_{j}*highly variable—in fact, for j = m – 1, C*^{ˆ}_{m}_{−}1 =

### (

*Y*1 −

*Y*

^{(}

*m*

^{)}

### )(

*Y*

*−*

_{m}*Y*

^{(}

*m*

^{)}

### )

, which is based on only one term in the sum*Remedy: Use different weights — spectral window — in sum for estimate *
of Var(*Y(m*))

**Regenerative Method **

*Different approach to gathering data — do not use Y(m*) as point estimator for *ν *
*Assume output process is regenerative: *

Can identify a sequence of random indices 1 ≤* B*_{1}* < B*_{2} < ... such that:

*Starting from each index B** _{j}*, the process from then on follows the same

*probability distribution as does the process from any other B*

_{j}*The process from index B*_{j}* on is independent of the process before B*_{j}*This divides the process into IID cycles (or tours): *

M

1 1

1 1

3 2

2

2 1

1

,..., ,

,..., ,

− +

− +

*B*
*B*

*B*

*B*
*B*

*B*

*Y*
*Y*

*Y*

*Y*
*Y*

*Y*

Most familiar example:

Single- or multiple-server queue

Any interarrival, service-time distributions

*Y*_{i}* = D*_{i}* = delay in queue of ith arriving customer *

*B*_{j}* = index of the jth customer who arrives to find the whole system empty of *
other customers and all servers idle

For single-server case, a customer begins a regeneration cycle if and only if delay in queue is zero (not true for multiple servers)

*i*

* *
*D*_{i}

* B*1=1 * B*2=6 * B*3=8 * B*4=13 * B*5=14 6=20

Comparable random variables defined on successive cycles are IID
*N*_{j}* = B*_{j + 1}* – B*_{j}* = length of jth cycle *

### ∑

^{−}

=

= + 1 1

*j*

*j*

*B*

*B*
*i*

*i*

*j* *Y*

*Z* * = sum of the observations in the jth cycle *

Can show (renewal reward theory): *ν = E(Z**j**)/E(N** _{j}*)

*Make a simulation run of n*′* cycles (don’t stop in the middle of a cycle) *
*For j = 1, 2, ..., n*′**, let U***j** = (Z**j**, N**j*)^{T}*; these are IID random vectors *

Point estimator for *ν: *

) (

) ˆ (

*n*
*N*

*n*
*Z*

′

= ′

*ν* *, which is biased, but strongly consistent *

Variance estimation, confidence interval:

Let **ΣΣ = **

22 12

12 11

*σ*
*σ*

*σ*

*σ* ** be the covariance matrix of U**_{j}

*Define V*_{j}* = Z** _{j}* –

*ν N*

*j,*which is unobservable, but has mean zero, variance

22 2 12 11

2 *σ* 2*νσ* *ν* *σ*

*σ** _{V}* = − +

*Central-limit theorem applied to V** _{j}*’s:

→

′

− ′

′ _{D}

*V* *n*

*n*
*N*
*n*

*Z*

2

) ( )

(
*σ*

*ν* * N(0, 1) as n*′→ ∞

Need to estimate *σ : *_{V}^{2}

Estimate **ΣΣ by **

### [ ][ ]

1

) ( )

( )

( ˆ ) ( ˆ

) ˆ ( )

ˆ ( _{1}

22 12

12 11

′−

− ′

− ′

=

′

′

′

′

### ∑

^{′}

=

*n*

*n*
*n*

*n*
*n*

*n*
*n*

*n*

*j*

*T*
*j*

*j* **U** **U** **U**

**U**
*σ*

*σ*

*σ*
*σ*

Let ^{(} ^{)} ^{ˆ} ^{(} ^{)} ^{2}^{ˆ}^{(} ^{)}^{ˆ} ^{(} ^{)}

### [

^{ˆ}

^{(}

^{)}

### ]

^{ˆ}22

^{(}

^{)}2

12 11

2 *n* *n* *n* *n* *n* *n*

*V* ′ =*σ* ′ − *ν* ′*σ* ′ + *ν* ′ *σ* ′

*σ* , which is strongly

consistent for *σ*_{V}^{2}

Can replace *σ in above CLT by *_{V}^{2} *σ*ˆ_{V}^{2}(*n*′)* (continuous mapping thm.): *

→

′

− ′

′ _{D}

*V* *n*

*n*
*N*
*n*

*Z*
ˆ2

) ( )

(
*σ*

*ν* * N(0, 1) as n*′→ ∞

*Manipulate to get asymptotically valid 100(1 – α)% confidence interval for ν: *

)
(
ˆ ^{2}

2 /

1 *n* *n*

*z* * _{V}* ′ ′

′ ± ^{−} *σ*

*ν* ^{α}

Appeals of regenerative method:

Firm mathematical foundation—asymptotic validity guarantees proper coverage
probability (1 – *α) as n′*→ ∞

Issues with regenerative method:

Underlying process must be regenerative — not universally true

Must identify regeneration points (proof), code their recognition into program, gather data differently

Asymptotic validity depends on having a lot of cycles — if cycles tend to be long (as they often do in complicated models) we may be able to observe only a few cycles, and asymptotic validity doesn’t kick in

*There are other ways to get a confidence interval (notably jackknifing; see text) *
Sequential-sampling versions have been developed—keep simulating more

*cycles until the confidence interval is small enough *

**Standardized Time Series **

Classical univariate statistics:

*Take IID sample X*_{1}*, X*_{2}*, ..., X** _{n}*
Want to estimate

*µ = E(X*

*i*)

“Standardize” univariate data:

) 1 , 0 ( ) N

( of deviation standard

of Estimate

unknown )

( − →*D*

*n*
*X*
*n*

*X* *µ*

Basis for inference (confidence intervals, hypothesis tests, ...)

Observing a process (via simulation):

*Observe the process Y*_{1}*, Y*_{2}*, ..., Y** _{m}*
Want to estimate lim (

*)*

_{i}*i* *E* *Y*

∞

= →

*ν*

“Standardize” process data:

→

− *D*

*i*

*Y*
*Y*

of i

deviation standard

of Estimate

unknown *ν*

*Brownian bridge process *
(Brownian bridge process is fully understood, like N(0, 1)

Basis for inference (confidence intervals, hypothesis tests, ...)

Specifics:

*Observe process Y*_{1}*, Y*_{2}*, ..., Y*_{m}

Point estimator: *Y(m*), like batch means, time-series models, spectral analysis
*Form n batches of size k each (m = nk), as for batch means *

*For large m, Y(m*) is approximately normal with mean *ν and variance τ*^{2}*/m *
where *τ*^{2} = lim Var(*Y*(*m*))

*m*→∞

Let

_{∑ ∑∑} ( )

= = = + −

−

= − ^{n}

*j*
*k*

*s*
*s*

*i*

*k*
*j*
*i*

*j* *k* *Y*

*k* *Y*
*A* *k*

1

2

1 1

) 1

312 ( ) (

*As k *→ ∞*, A/τ*^{2} → *χ*^{2}* distribution with n d.f., and is independent of Y(m*)
*So for large k, *

### ( )

) (

) ) (

(

2 2

*mn*
*A*

*m*
*Y*

*n*
*A*

*m*
*m*

*Y* *ν*

*τ*
*τ*

*ν* = −

−

*has an approximate t distribution with n d.f., so an asymptotically valid *
confidence interval for *ν is Y*(*m*)±*t*_{n}_{,}_{1}_{−}_{α}_{/}_{2} *A* (*mn*)

Appeals of standardized time series:

Firm mathematical foundation—asymptotically valid intervals
*Assumptions are much weaker than for regenerative method *
Relatively simple

Issues with standardized time series:

As with all asymptotic, methods, how long is long enough?

*Above is called “area” approach — A represents area under the standardized *
Brownian bridge

*Other approaches look instead at the maximum attained by the standardized *

**Summary of Steady-State Estimation Procedures **

Nothing will work well if computational budget is unduly limited

Batch means, spectral analysis, standardized time series display generally good performance (with respect to coverage probability, efficiency of data usage) Simplicity might argue against spectral analysis

Sequential-sampling versions exist for most of the main methods, to control confidence-interval width

Performance measures other than means (variances, proportions, quantiles) have been investigated

It may be difficult to implement these methods in the context of some existing simulation software, though some software does allow for, and even builds in and automates, some of these methods

**9.5.4 Estimating Other Measures of Performance **

Means do not always provide the one and only appropriate measure of performance

*Probabilities: For some set B, estimate p = P(Y *∈* B), where Y is the steady-state *
*random variable of the process Y*_{1}*, Y*_{2}, ...

*e.g., Y = delay of a message, B = (0, 5 minutes), so p is the probability that a *
message is delayed by less than 5 minutes

This is a special case of estimating means, if we define the indicator random variable

∈

= 0 otherwise if

1 *Y* *B*

*Z* *, then p = P(Y *∈* B) = P(Z = 1) = 1*×*P(Z = 1) + *
0×*P(Z = 0) = E(Z) *

Thus, can use all the methods described above for means to estimate a proportion

*Quantiles: The q-quantile y**q** is the value such that P(Y *≤* y**q**) = q *

*e.g., Y = delay of a message, y*_{0.75} is the value below which are 75% of the
message delays

See text for details on methods for estimating quantiles based on order statistics, batch means, and regenerative methods

**9.6 Statistical Analysis for Steady-State ** **Cycle Parameters **

To dichotomize simulations into terminating vs. steady-state is a bit of an oversimplification:

Manufacturing model, has 8-hour shifts

Simulating in detail what goes on within a shift Performance might fluctuate widely within a shift

*But what matters is the production (say) over the whole shift *

*Aggregate the data so that a “basic” observation Y*_{i}* is the production on the ith *
*shift *

*Steady-state cycle (e.g., cycle = shift) parameter: steady-state mean of process *
defined over a cycle as a “unit” of time

Use above steady-state methods, except applied to random variables defined over a cycle, rather than individually at their finest level of detail

Example:

Production process, shift = 8 hours, lunch break 4 hours into each shift
*N*_{i}* = production in hour i; would expect it to drop in 4th hour; plot averaged *

over 10 replications:

*C*

*N = production in shift i; plot averaged over 10 replications: **i*

No steady- state

Seems to have a steady- state; use earlier methods

**9.7 Multiple Measures of Performance **

Usually: Want to observe several performance measures from a large simulation Average length of queue(s)

Maximum length of queue(s) Utilization(s)

Throughput(s) Difficulty:

Estimate each (expected) performance measure with a confidence interval Any of the intervals could “miss” its expected performance measure

*Must be careful about overall statements of coverage (i.e., that all intervals *
*contain their expected performance measures simultaneously) *

*Sometimes called the problem of multiple comparisons *

*Have k output performance measures, want overall (familywise) probability of at *
least 1 – *α that the confidence intervals for them all contain their target expected *
performance measures

*For s = 1, 2, ..., k, suppose the confidence interval for performance measure s is at *
confidence level 1 – *α*s

*Then P(all intervals contain their respective performance measures) *≥ 1 –

### ∑

=
*k*

*s*
*s*
1

*α *
*(Bonferroni inequality) *

Thus, pick *α**s*’s so that

### ∑

*α*=

*α*

=
*k*

*s*
*s*
1

Could pick *α**s* = *α/k for all s, or pick α**s*’s differently with smaller *α**s*’s for the more
important performance measures

**Alternative Approach **

*Use multivariate statistical methods to get a joint confidence region (not necessarily *
*a rectangle) that covers expected-performance-measure vector with probability 1 *
– *α *

Example:

*k = 2 performance measures µ*1 and *µ*2

Bonferroni approach: Separate intervals for *µ*1 and *µ*2

Multivariate approach: Ellipse containing (*µ*1, *µ*2) with probability 1 – *α *

*µ*1

*µ*_{2}

100(1−*α*1)% confidence
interval for *µ*_{1}

100(1−*α*_{2})% confidence
interval for *µ*2

Bonferroni " box "

≥100(1−*α)% confidence*

region ^{ }

100(1−*α)% *

confidence ellipse

Specific multivariate methods:

Multivariate batch means Multivariate spectral analysis Multivariate time-series methods Appeal of multivariate methods:

Smaller area (volume) with multivariate methods

**9.8 Time Plots of Important Variables **

So far, have concentrated on performance measures over the course of the whole simulation run

Averages Variances Extrema Proportions Quantiles

But these may mask important dynamic (within-run) behavior patterns Periodicities

Explosions Learning

Ways to pick up such dynamic behavior

Plot output processes (discrete or continuous) Animation