Simulation of Stochastic Systems

The Central Idea

Suppose a stochastic model produces a random output \(Y\), and we want to know a quantity such as

If exact analysis is difficult, generate independent simulated outputs

from the model and use the sample average

to estimate the true expectation.

Simulation Is an Experiment on a Model

A physical experiment observes a real system. A stochastic simulation performs repeated experiments on a mathematical model of that system.

Simulation therefore answers questions about the model as specified.

If the model assumptions are inaccurate, highly precise simulation can still produce an inaccurate description of reality.

Three Layers of Simulation

A useful simulation workflow separates three layers:

1. Model: the probability assumptions and state-transition rules.
2. Random generation: the mechanism used to draw random inputs.
3. Estimation: the statistics computed from simulated outputs.

Errors at any layer can invalidate the result.

Pseudorandom Numbers

Computer simulations usually begin with a pseudorandom number generator.

The generator creates a deterministic sequence designed to behave like random draws for computational purposes.

A basic source is modeled as

Uniform Random Numbers as Building Blocks

Uniform random variables on \((0,1)\) can be transformed into many other distributions.

This makes

a fundamental computational building block.

Simulating a Bernoulli Variable

To simulate

generate

and define

Simulating a Finite Discrete Distribution

Suppose

Partition the unit interval into subintervals having lengths \(p_1,\ldots,p_m\). Generate \(U\sim\operatorname{Uniform}(0,1)\) and return the value corresponding to the interval containing \(U\).

Cumulative-Probability Method

Define cumulative probabilities

Then return \(x_j\) when

Inverse Transform Sampling

Let \(F\) be a continuous cumulative distribution function.

If \(F^{-1}\) is available and

then

has distribution \(F\).

Example: Simulating an Exponential Variable

For

the cumulative distribution is

Setting \(U=F(T)\) and solving gives

Since \(1-U\) is also uniform on \((0,1)\), one often writes

Random Seeds

A random seed determines the starting state of a pseudorandom number generator.

Reusing the same seed with the same implementation usually reproduces the same sequence of pseudorandom values.

This is useful for debugging, auditing, and reproducibility.

Reproducibility Does Not Mean True Randomness

Pseudorandom sequences are reproducible precisely because they are generated algorithmically.

Statistical simulation requires that the generator behave suitably for the intended application, not that each number arise from a physically random mechanism.

Sample Paths

A stochastic process simulation produces a sample path:

in discrete time, or an event/state trajectory through continuous time.

Repeating the simulation produces many possible paths through the same model.

One Path Is Not the Distribution

A single sample path can be informative visually, but it does not reveal the complete probability distribution.

Rare events may not appear, and unusually favorable or unfavorable paths can be misleading.

Statistical conclusions usually require many replications or a sufficiently long path under justified ergodic assumptions.

Monte Carlo Estimation

Monte Carlo simulation estimates mathematical quantities by random sampling.

If

and \(X_1,\ldots,X_N\) are simulated independently from the model, then

is the standard Monte Carlo estimator.

Estimating a Probability

To estimate

define an indicator variable

Then

Example: Estimating an Extinction Probability

Suppose \(10{,}000\) independent branching-process simulations are run and \(3{,}180\) become extinct within a sufficiently justified simulation horizon.

The empirical extinction frequency is

A finite time horizon can still introduce additional truncation bias when eventual extinction is the target.

The Law of Large Numbers

Under standard independence and integrability assumptions,

converges toward

as the number of replications grows.

This is the theoretical foundation of ordinary Monte Carlo estimation.

Monte Carlo Error

Even when the estimator is unbiased, a finite simulation has random error.

If

then

Therefore,

The Square-Root Rule

Monte Carlo standard error decreases proportionally to

To cut the standard error in half, approximately four times as many independent replications are needed.

To reduce it by a factor of ten, roughly one hundred times as many replications are needed.

Estimating the Standard Error

The true standard deviation \(\sigma\) is often unknown.

Use the sample standard deviation

Then estimate

Approximate Confidence Intervals

For a large independent Monte Carlo sample, an approximate \(95\%\) confidence interval for \(E[Y]\) is

This interval quantifies simulation uncertainty, not model uncertainty.

Simulation Error Versus Model Error

Two distinct uncertainties must be separated:

• simulation error: error caused by using finitely many random runs;
• model error: error caused by inaccurate assumptions about the real system.

Increasing \(N\) reduces simulation error but does not fix a poor model.

Bias

An estimator \(\widehat{\theta}\) has bias

Simulation bias can arise from truncation, initialization, incorrect random generation, discarded events, or implementation errors.

Variance

Two unbiased simulation procedures can have very different variances.

Lower variance means less random uncertainty for the same computational effort.

More advanced Monte Carlo methods deliberately reduce variance without changing the target quantity.

Mean Squared Error

A common overall measure is

A low-variance estimator can still be poor if it is substantially biased.

Simulating a Random Walk

For a random walk with

simulate each step from a Bernoulli draw and update

The resulting sequence is one random-walk sample path.

Random-Walk Simulation Questions

Repeated walks can estimate:

• \(P(X_n=k)\);
• probability of hitting a boundary;
• first-passage-time distributions;
• maximum displacement;
• return frequencies;
• effects of changing the drift.

Simulating a Finite Markov Chain

Suppose the transition matrix is

To simulate one step from current state \(i\):

1. read row \(i\) of the transition matrix;
2. treat that row as a categorical probability distribution;
3. sample the next state from those probabilities;
4. repeat using the newly sampled state.

Markov-Chain Path Simulation

If the initial state is sampled from

and every later transition is sampled from the appropriate row of \(P\), the resulting path has the probability law specified by the Markov model.

Estimating a Stationary Distribution

For an irreducible finite Markov chain, a long simulated path can estimate stationary probabilities using empirical state frequencies:

The successive observations are dependent, so ordinary independent-sample standard-error formulas require modification.

Burn-In

If a Markov-chain simulation begins far from equilibrium, early observations can reflect the initial state strongly.

Analysts may discard an initial burn-in segment before estimating stationary behavior.

Burn-in is a practical device, not a substitute for checking convergence.

Autocorrelation

Successive observations along one stochastic-process path are often dependent.

If nearby observations are strongly correlated, \(N\) recorded values contain less information than \(N\) independent observations would.

This affects uncertainty estimates.

Effective Sample Size

The effective sample size is a way to describe how much independent-information equivalent is contained in a dependent simulation.

Strong positive autocorrelation generally lowers effective sample size.

Exact formulas depend on the correlation structure and are beyond this introductory treatment.

Simulating a Poisson Process

One direct method uses exponential interarrival times.

Generate

and form arrival times

Stop when the next arrival exceeds the simulation horizon.

Alternative Poisson Simulation Over a Fixed Interval

To simulate all arrivals on \([0,t]\), one may:

1. sample \(N\sim\operatorname{Poisson}(\lambda t)\);
2. conditional on \(N=n\), sample \(n\) independent uniform points on \([0,t]\);
3. sort those points to obtain arrival times.

This uses the conditional uniform-arrival property of the Poisson process.

Simulating a Branching Process

Begin with \(Z_0\). For each current individual, independently generate an offspring count from the specified distribution.

Sum the offspring counts:

Stop at extinction or at a chosen generation horizon.

Branching Simulation and Truncation

If the target is eventual extinction, stopping a simulation after a finite number of generations can classify some processes as surviving even though they would become extinct later.

This creates horizon-dependent bias unless the remaining probability is accounted for.

Simulating Reliability Lifetimes

If component lifetime has distribution \(F\), generate a lifetime using an appropriate random-generation method.

For exponential lifetime,

Compare \(T\) with the mission duration to determine whether the component survives the mission.

Simulating Repairable Reliability

A two-state repairable system can alternate between:

and

Generate exponential up-times and repair times and track the fraction of the simulated horizon spent in the operational state.

Estimating Availability

If a simulated repairable system is operational for total time \(T_U\) over a long simulation horizon \(T\), estimate availability by

Simulating an \(M/M/1\) Queue

For an \(M/M/1\) queue:

• interarrival times are exponential with rate \(\lambda\);
• service times are exponential with rate \(\mu\).

The simulation tracks arrival times, service-start times, completion times, queue lengths, and waiting times.

Recursive Waiting-Time Update

Let \(A_n\) be the arrival time of customer \(n\), \(S_n\) the service time, and \(D_n\) the departure time.

In a single-server FCFS queue,

Service begins when both the customer and the server are available.

Queue Waiting Time from the Recursion

The waiting time before service is

The total time in the system is

Discrete-Event Simulation

Queueing and reliability simulations are often discrete-event simulations.

The simulated clock jumps directly to the next event rather than advancing through every tiny time increment.

This can be far more efficient than fixed-step simulation when events are sparse.

Event Calendar

A discrete-event simulator often maintains an event calendar containing future scheduled events.

The algorithm repeatedly:

1. selects the earliest event;
2. advances the clock to that event time;
3. updates the system state;
4. schedules any newly created future events.

Simulating a Martingale

For a symmetric random walk martingale, generate many paths and compute the empirical mean at each time:

The values should fluctuate around

This illustrates martingale behavior but does not prove it.

Independent Replications

One common simulation design runs the entire system repeatedly from its initial condition.

If the random-number streams are independent across replications, then the resulting outputs can often be analyzed with ordinary independent-sample methods.

One Long Run

Another design uses one long stochastic-process trajectory.

This is natural for estimating steady-state behavior, but observations along the path are generally dependent.

Long-run averages require structural assumptions such as ergodicity.

Transient Simulation

A transient simulation studies system behavior over a specified finite horizon from a particular initial condition.

Examples include:

• queue length during the first two hours after opening;
• probability of system failure during a 10-hour mission;
• population size after 20 generations;
• random-walk position after 100 steps.

Steady-State Simulation

A steady-state simulation seeks long-run quantities after transient effects have become negligible.

Examples include:

• long-run server utilization;
• stationary state probabilities;
• long-run availability;
• average waiting time in an operating system.

Warm-Up Bias

A steady-state simulation often begins from an artificial initial state such as an empty queue.

Early observations can therefore be systematically unrepresentative of equilibrium behavior.

Removing a warm-up segment is one way to reduce this initialization bias.

Choosing a Warm-Up Length

There is no universal warm-up length.

It should be chosen using model knowledge, convergence diagnostics, repeated runs, or other evidence that initialization effects have sufficiently diminished.

Choosing a Run Length

A run should be long enough to estimate the target quantity with acceptable precision.

Longer runs reduce random error, but highly persistent systems may require especially long simulations because observations remain correlated.

Rare Events

Ordinary simulation can be inefficient for extremely rare events.

If

then millions of independent runs may be needed just to observe a modest number of occurrences.

Specialized rare-event methods are used in advanced simulation.

Variance Reduction

Advanced Monte Carlo methods can reduce estimator variance without increasing the number of runs proportionally.

Examples include:

• antithetic variables;
• control variates;
• common random numbers;
• importance sampling;
• stratified sampling.

These ideas are introduced here only as extensions.

Common Random Numbers

When comparing two system designs, using matched random inputs can reduce the variance of the estimated difference.

The goal is not to make the two systems independent, but to expose them to comparable simulated conditions.

Verification

Verification asks:

Did we implement the mathematical model correctly?

Verification checks code logic, equations, units, indexing, state updates, and event scheduling.

Validation

Validation asks:

Is the mathematical model an adequate representation of the real system for the intended purpose?

A simulation can be perfectly verified and still fail validation.

Use Known Cases as Tests

One of the strongest implementation checks is to compare simulation output with cases having known exact answers.

Examples include:

• binomial probabilities for a simple random walk;
• stationary probabilities of a small Markov chain;
• Poisson count means;
• \(M/M/1\) queue formulas;
• exponential reliability formulas.

Sanity Checks

Useful basic checks include:

• probabilities lie between \(0\) and \(1\);
• probability vectors sum to \(1\);
• counts are nonnegative integers;
• event times are ordered;
• queue lengths never become negative;
• absorbing states remain absorbing;
• units remain consistent.

Convergence Diagnostics

A simulation estimate should be monitored as computational effort increases.

Useful diagnostics include:

• running averages;
• estimates from several independent replications;
• estimated standard errors;
• comparison across run lengths;
• autocorrelation diagnostics for long-run simulations.

Running Average

The running average after \(n\) observations is

Plotting or recording \(\overline{Y}_n\) can reveal whether an estimate is still changing substantially as more simulation output is added.

Batch Means

For one long dependent simulation, output can be divided into large batches.

The mean within each batch is computed, and the batch means can sometimes be treated approximately like less-correlated observations for uncertainty estimation.

Batch size must be chosen carefully.

Replication Versus Precision

More replications improve statistical precision, but computational cost also increases.

Simulation design is therefore an allocation problem: how much computational effort is required to answer the decision question accurately enough?

Sensitivity Analysis

A stochastic model should often be simulated under several plausible parameter values.

This reveals whether conclusions are robust or depend strongly on uncertain assumptions such as:

• arrival rates;
• service rates;
• failure rates;
• transition probabilities;
• offspring distributions.

Scenario Analysis

Different system designs or operating policies can be represented as scenarios.

Simulation can compare:

• staffing levels;
• redundancy configurations;
• maintenance policies;
• queue disciplines;
• capacity limits;
• transition rules.

Simulation Is Not Optimization by Itself

A simulation evaluates the performance of specified choices.

Optimization requires an additional procedure for searching among choices and defining an objective.

Simulation can support optimization, but it is not automatically an optimizer.

Simulation Is Not Prediction Without Assumptions

A simulation generates consequences of a probability model.

It does not establish that those probabilities will govern future reality.

Forecasting requires both a stochastic model and evidence that the model is appropriate for the forecasting context.

Common Error: Trusting One Simulation Run

One run is one random realization.

It may be unusually high, low, short, long, stable, or unstable.

Use repeated runs or justified long-run averaging when statistical inference is the goal.

Common Error: Reporting an Estimate Without Uncertainty

A simulation estimate such as

is incomplete when the precision matters.

Whenever practical, report an estimated standard error, confidence interval, or other measure of simulation uncertainty.

Common Error: Increasing \(N\) to Fix Model Bias

More simulation runs reduce Monte Carlo noise.

They do not correct:

• wrong transition probabilities;
• wrong arrival assumptions;
• wrong service distributions;
• incorrect independence assumptions;
• coding errors.

Common Error: Confusing Seed Choice with Statistical Replication

Reusing one seed reproduces one pseudorandom sequence.

Independent statistical replications require appropriately distinct random streams or generator states.

Common Error: Ignoring Initialization

Starting every queue empty or every Markov chain in one state can create short-run behavior that differs from equilibrium.

Decide explicitly whether the target is transient or steady-state behavior.

Common Error: Ignoring Dependence

Observations from one stochastic-process path are generally not independent.

Treating them as independent can understate uncertainty when positive autocorrelation is present.

Common Error: Stopping When the Result Looks Right

Ending a simulation because the current estimate happens to match an expected answer can introduce subjective bias.

Simulation length and stopping rules should be defined by precision, computational budget, or other principled criteria.

Common Error: Forgetting Rare Outcomes

If a rare event does not appear in a modest simulation, its probability is not necessarily zero.

Absence of observed occurrences provides limited information when the event is sufficiently rare.

A Reliable General Simulation Procedure

1. Define the stochastic model precisely.
2. Identify the target probability, expectation, distribution, or path quantity.
3. Choose a valid method for generating each random input.
4. Specify the initial state and simulation horizon.
5. Implement the state-transition or event-update rules.
6. Record the outputs needed for the target statistic.
7. Repeat independently or run long enough for the intended long-run analysis.
8. Compute the estimate and an uncertainty measure.
9. Verify the implementation using known cases and sanity checks.
10. Validate the model assumptions against the intended application.

A Reliable Procedure for Monte Carlo Probabilities

1. Define the target event \(A\).
2. Simulate one complete realization of the model.
3. Record \(I=1\) if \(A\) occurs and \(I=0\) otherwise.
4. Repeat for \(N\) independent replications.
5. Estimate \(P(A)\) by the fraction of successes.
6. Estimate uncertainty using the Bernoulli sample variance or an appropriate interval.
7. Increase \(N\) if the precision is insufficient.

A Reliable Procedure for Steady-State Simulation

1. Check that a meaningful steady state is expected to exist.
2. Choose an initial state.
3. Run the process long enough for initialization effects to diminish.
4. Apply a justified warm-up rule if needed.
5. Collect long-run observations.
6. Account for serial dependence when estimating uncertainty.
7. Repeat with longer runs or independent replications to check stability.

Example: Monte Carlo Estimate of a Random-Walk Event

Suppose a symmetric random walk starts at \(0\), and we want to estimate

One replication:

1. generate \(20\) independent \(\pm1\) steps;
2. sum them to obtain \(X_{20}\);
3. record \(I=1\) if \(X_{20}\ge4\), otherwise \(I=0\).

After \(N\) replications,

This can then be compared with the exact binomial calculation as a verification check.

Practice 1: Bernoulli Simulation

Describe how to use

to simulate a Bernoulli random variable with success probability

Practice 2: Exponential Simulation

Give an inverse-transform formula for simulating

from a uniform random variable \(U\).

Practice 3: Monte Carlo Mean

Simulated outputs are

Find the Monte Carlo estimate of the expected output.

Practice 4: Probability Estimate

A target event occurs in \(236\) of \(1000\) independent simulations.

Estimate its probability.

Practice 5: Square-Root Rule

A simulation uses \(2{,}500\) replications. Approximately how many replications are required to cut the Monte Carlo standard error in half?

Practice 6: Markov-Chain Simulation

A finite Markov chain is currently in state \(i\). Which part of the transition matrix is used to sample the next state under the row-stochastic convention?

Practice 7: Poisson-Process Simulation

What distribution should be used for consecutive interarrival times in a homogeneous Poisson process of rate \(\lambda\)?

Practice 8: Queue Simulation

In a single-server FCFS simulation, a customer arrives at time \(A_n\), the previous customer departs at \(D_{n-1}\), and the new service time is \(S_n\).

Write the recurrence for \(D_n\).

Practice 9: Verification Versus Validation

What is the difference between verifying a simulation and validating its model?

Practice 10: Identify the Error

A student runs one stochastic simulation, obtains a waiting time of \(12.4\) minutes, and reports:

“The expected waiting time is exactly \(12.4\) minutes.”

Explain the error.

Practice Answers

Practice 1

Generate \(U\sim\operatorname{Uniform}(0,1)\) and define

Then

Practice 2

For rate \(5\),

has an exponential distribution with rate \(5\).

Practice 3

Practice 4

Practice 5

Halving standard error requires about four times as many independent replications:

Practice 6

Use row \(i\) of the transition matrix. Its entries form the categorical distribution for the next state.

Practice 7

Consecutive interarrival times should be independent

random variables.

Practice 8

Practice 9

Verification asks whether the intended mathematical model was implemented correctly.

Validation asks whether that mathematical model is an adequate representation of the real system for the intended use.

Practice 10

One run produces one random realization, not the expected value.

The expected waiting time must be estimated from many replications or a justified long-run average, and the estimate should generally be accompanied by a measure of simulation uncertainty.

What to Remember

Simulation turns a stochastic model into artificial random realizations that can be summarized statistically.

• Uniform pseudorandom numbers are basic inputs for generating other random variables.
• A random seed supports reproducibility.
• Inverse-transform sampling converts uniform draws into many continuous distributions.
• Monte Carlo estimation uses simulated sample averages and event frequencies.
• The law of large numbers explains why independent Monte Carlo averages stabilize.
• Ordinary Monte Carlo standard error decreases like \(1/\sqrt N\).
• More simulation reduces sampling noise but does not correct a wrong model.
• Random walks, Markov chains, branching processes, Poisson processes, queues, reliability processes, and martingales can all be studied through simulation.
• Discrete-event simulation advances from event to event rather than through every instant.
• Transient and steady-state simulations answer different questions.
• Warm-up bias and serial dependence matter in long-run process simulation.
• Verification checks implementation; validation checks model adequacy.
• Known exact cases provide valuable tests of simulation code.
• Estimates should be accompanied by uncertainty when precision matters.
• Rare events can require specialized methods beyond ordinary Monte Carlo.

The central computational principle is: generate according to the model, repeat enough to quantify random uncertainty, and validate both the implementation and the assumptions before trusting the result.

Where We Go Next

This lesson completes the stochastic-process sequence by turning the major models of the section into computational experiments. The same simulation principles will reappear throughout modeling, statistics, numerical methods, optimization, and applied mathematical analysis.

The next stage is to use these stochastic tools as components inside broader mathematical models: defining assumptions, generating outcomes, measuring uncertainty, comparing scenarios, and testing how conclusions change when model structure or parameters change.