1 Modelling context

The analysis in Lessons 1 and 2 is exact: given the transition matrix, the state distributions at all times can be computed by matrix multiplication, and the stationary distribution can be solved for exactly. This raises a fair question: why simulate?

Simulation serves two purposes. First, it verifies the implementation: an exact calculation and a simulated estimate should agree, and a discrepancy indicates an error in one or the other. Second, it generalises: many real systems are too complex for exact analysis (inhomogeneous chains, chains with very large state spaces, chains coupled to other processes), and simulation is the only available method. The techniques developed here carry over to those more complex settings.

2 Learning objectives

3 Prerequisite review

Inverse transform

The method for generating samples from a discrete distribution: draw U ~ Uniform(0, 1), then find the smallest index i such that the cumulative probability exceeds U.

Pseudorandom generation

The generation of uniform values, covered in Module 4. Each simulation run should record its seed.

Sampling error

The random discrepancy between a simulated estimate and the true value. Decreases as 1/√N.

Ergodic theorem

For an ergodic Markov chain, the empirical frequency of visits to a state converges to the stationary probability as the number of steps grows.

4 Mathematical formulation

Simulation algorithm

To simulate a Markov chain with transition matrix P from an initial distribution π^(0) for N steps:

1. Choose the initial state. Draw U ~ Uniform(0, 1) and use inverse transform with π^(0) to select the starting state.
2. For each step t = 1, …, N:
   • Draw U ~ Uniform(0, 1).
   • Use inverse transform with the current state's row of P to determine the next state: find the smallest j such that Σ_{k≤j} P_{i,k} ≥ U, where i is the current state.
   • Record the next state.

The result is a path (X_0, X_1, X_2, …, X_N).

Estimating the stationary distribution

By the ergodic theorem, the empirical frequency of each state converges to its stationary probability:

π*_j ≈ (1/N) Σ_{t=1}^N 1{ X_t = j }

For finite N, the estimate differs from π*_j by sampling error of order 1/√N. Replication with independent seeds quantifies the error.

Estimating the n-step distribution

To estimate the distribution at step n, run many independent replications, each starting from the same initial state, and record the state at step n:

π̂_j^(n) = (1/R) Σ_{r=1}^R 1{ X_n^(r) = j }

where R is the number of replications. By the law of large numbers, this converges to π_j^(n) as R grows. Compare against the exact P^n as a verification.

Choosing between simulation and exact analysis

Exact analysis is preferred when it is available: it is exact, deterministic, and requires no seed management. Simulation is preferred when:

• The state space is very large (e.g. thousands of states) and the matrix power P^n is expensive to compute.
• The transition probabilities depend on time or on external covariates, so no single matrix P describes the chain.
• The chain is coupled to other processes (e.g. queueing, agent-based models) that complicate the analysis.
• The quantity of interest is a path functional (e.g. total time spent above a threshold) that is easier to estimate by simulation than analytically.

5 Worked example

       S    C    R
S   [ 0.7  0.2  0.1 ]
C   [ 0.3  0.4  0.3 ]
R   [ 0.2  0.3  0.5 ]

6 Verification check

7 Interpretation

Simulation of a Markov chain is a form of Monte Carlo estimation (Module 4) applied to a stochastic process. The path is the sample; the state frequencies are the estimator; the stationary distribution is the target.

Three interpretive points:

• Comparison with the theoretical distribution is a verification. When the exact P^n is available, agreement between simulation and theory verifies both the theory and its implementation in code.
• The empirical frequency is a Monte Carlo estimate. It carries the same sampling error as any other Monte Carlo estimator: order 1/√N.
• Simulation scales to problems that theory cannot handle. Once the algorithm is verified against a simple case, it can be applied to problems with inhomogeneous transitions, large state spaces, or coupling to other processes.

8 Uncertainty and limitations

The uncertainty in a simulated estimate of the stationary distribution comes from two sources.

9 Common implementation errors

10 Independent modelling task

11 Reflection prompt

12 Lesson summary

§ Self-quiz, practice, and homework