Checking Markov Chain Simulations Against Exact Probabilities

Exact Probabilities as Simulation Benchmarks

A simulated Markov chain produces one random path, not a list of probabilities. Its observed transitions will generally differ from the probabilities used to generate them. To decide whether a simulation is behaving as intended, compare its frequencies with exact quantities computed from the transition matrix, and ask whether the observed discrepancy is plausible as sampling variation. This is a diagnostic, not a requirement that simulated counts match their expected values exactly.

Let \((X_n)_{n\ge0}\) be a time-homogeneous Markov chain on a finite state space \(S\), defined on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Its transition matrix is \(P=(p_{ij})\), where \(p_{ij}=\mathbb{P}(X_{n+1}=j\mid X_n=i)\). The Markov property says that this conditional probability is unchanged when the earlier history is also given. Each entry is nonnegative and every row sums to 1. Write \(\mu\) for the initial distribution as a row vector, so the distribution at time \(n\) is \(\mu P^n\).

For a simulated path \(X_0,\ldots,X_N\), define the number of departures from state \(i\), and the number of observed transitions from \(i\) to \(j\), by

When \(C_i>0\), the empirical one-step frequency from \(i\) to \(j\) is \(C_{ij}/C_i\). It can be compared with \(p_{ij}\). If \(C_i=0\), that empirical frequency is undefined: the path did not visit \(i\) at a time when a transition could be observed. The denominator is therefore part of the diagnosis, not a detail to omit.

Proof. Use the filtration \(\mathcal{F}_n=\sigma(X_0,\ldots,X_n)\), which records the path through time \(n\). For \(1\le r\le N\), put \(Y_r=\mathbf{1}_{\{X_{r-1}=i\}}\bigl(\mathbf{1}_{\{X_r=j\}}-p_{ij}\bigr)\). Then \(D_{ij}=\sum_{r=1}^N Y_r\). By the Markov property, conditional on \(\mathcal{F}_{r-1}\), the probability that \(X_r=j\) is \(p_{ij}\) if \(X_{r-1}=i\). If \(X_{r-1}\ne i\), the indicator multiplying the parenthesis is zero. Thus \(\mathbb{E}[Y_r\mid\mathcal{F}_{r-1}]=0\), and taking expectations gives \(\mathbb{E}[D_{ij}]=0\).

For \(r<s\), \(Y_r\) is determined by \(\mathcal{F}_{s-1}\). Therefore \(\mathbb{E}[Y_rY_s]=\mathbb{E}\bigl[Y_r\mathbb{E}[Y_s\mid\mathcal{F}_{s-1}]\bigr]=0\). The cross terms in the variance of the sum vanish. Conditional on \(\mathcal{F}_{r-1}\), the squared term has expectation \(\mathbf{1}_{\{X_{r-1}=i\}}p_{ij}(1-p_{ij})\), since an indicator with success probability \(p_{ij}\) has variance \(p_{ij}(1-p_{ij})\). Summing these expectations gives \(\operatorname{Var}(D_{ij})=p_{ij}(1-p_{ij})\mathbb{E}[C_i]\). Finally, \(0\le C_i\le N\), so \(\mathbb{E}[C_i]\le N\). This proves the stated bound. \(\square\)

The result explains one useful feature of a one-step diagnostic: after accounting for the number of times the chain departed from \(i\), the observed count has the correct transition probability as its reference value. But the theorem concerns the residual \(C_{ij}-p_{ij}C_i\), not the ratio \(C_{ij}/C_i\). Also, the variance depends on the expected number of departures, which can be small when the chain rarely visits \(i\). A large discrepancy from a rarely visited state deserves a different interpretation from the same discrepancy after many departures.

Worked One-Step Comparisons

One-step frequencies are a direct diagnostic of how transitions are generated. Multi-step frequencies check a different feature: whether the simulated chain moves through time according to repeated application of the same transition rule. The exact benchmark follows by matrix multiplication, not by treating successive transitions as independent.

Proof. At \(m=0\), the distribution is \(\mu=\mu P^0\), since \(P^0\) is the identity matrix. Suppose the distribution at time \(m\) is \(\mu P^m\). If its \(i\)th entry is the probability of state \(i\) at time \(m\), then the probability of state \(j\) at time \(m+1\) is \(\sum_i(\mu P^m)_i p_{ij}\), by conditioning on \(X_m=i\). This sum is the \(j\)th entry of \((\mu P^m)P=\mu P^{m+1}\). Induction proves the result for every \(m\). For a chain started at \(i\), take \(\mu\) to put probability 1 on \(i\); the \(j\)th entry of \(\mu P^m\) is then \((P^m)_{ij}\). \(\square\)

Worked Multi-Step Comparisons

For independent chains run from a specified initial distribution, the distribution of each endpoint at a fixed lag is exactly the matrix-based benchmark. For one long simulated path, however, endpoints at different times—or overlapping windows of length \(m\)—are generally dependent. A collection of lagged observations from one path must not be treated as if it were automatically a collection of independent simulations. The exact probability remains a valid reference, but dependence affects how much variation to expect in the collected frequencies.

What Discrepancies Can Tell You

A comparison is most useful when the simulation and the benchmark describe the same quantity. For a one-step check, count transitions out of a specified state and compare the row frequencies with that state’s row of \(P\). For a lag-\(m\) check, use \(P^m\) and specify the starting distribution. If the simulation begins from a fixed state \(i\), use row \(i\); if its initial states follow \(\mu\), use \(\mu P^m\). Ignoring the starting distribution can make a correct simulation appear incorrect.

Several patterns are informative. A transition frequency that differs slightly from its benchmark may reflect ordinary randomness, particularly when there were few departures from that state. A row of empirical frequencies that does not sum to 1 suggests a counting or denominator error. A persistent discrepancy in a frequently visited row, repeated across independent runs, is more concerning than one isolated difference. A mismatch in lagged frequencies, even when one-step frequencies look reasonable, can point to a problem in how the code advances the chain or records its state.

Common mistake: squaring entries instead of multiplying matrices. For the first row of \(P\), squaring entries individually gives \((\tfrac{9}{16},\tfrac1{16})\), which sums to \(\tfrac{10}{16}\), not 1. These are not the two-step probabilities. The correct first entry is \(\tfrac34\cdot\tfrac34+\tfrac14\cdot\tfrac12=\tfrac{11}{16}\), because both possible intermediate states must be included. The other entry is \(\tfrac34\cdot\tfrac14+\tfrac14\cdot\tfrac12=\tfrac5{16}\); together the correct probabilities sum to 1.

Before interpreting any discrepancy, verify that the transition matrix is the intended one and that every row sums to 1. Then verify the arithmetic for the benchmark, including the multiplication order and initial distribution. Finally, check how many relevant transitions or endpoints were observed and whether the observations come from independent runs or from one dependent path. These checks separate basic implementation mistakes from the ordinary variation that simulation is designed to produce.