M.7.5 Markov Chain Monte Carlo

Basic Prerequisites: 03 · Posterior Inference and Numerical Methods · Monte Carlo.

---

1. Intuition & Practical Objective

Conjugate priors give closed-form posteriors - but the interesting models (logistic regression, stochastic volatility, hierarchical models, anything with a nonlinear link) have no closed form, and the normalizing constant $m(x)=\int f(x\mid \theta )\pi (\theta )d\theta$ is an intractable integral in more than a few dimensions. Markov Chain Monte Carlo (MCMC) sidesteps the integral entirely: it constructs a Markov chain whose stationary distribution is exactly the posterior, and just runs it. Every sample is a draw from $p(\theta \mid x)$ without ever computing $m(x)$.

The practical objective is to build and diagnose the two workhorse samplers: Metropolis–Hastings (propose, then accept or reject) and Gibbs (cycle through the full conditionals). Knowing when each is available - and how to tell that the chain has not converged - is the difference between a Bayesian analysis and a random-number generator.

The one-sentence essence. “Design a Markov chain with the posterior as its stationary distribution - MH does it with a proposal and an accept/reject step, Gibbs does it by sampling each full conditional - then prove convergence before you believe a single number.”

---

2. Mathematical Ground Truth & Derivations

The target and the idea (Tsay §12.1). We cannot sample $p(\theta \mid x)\propto g(\theta )=f(x\mid \theta )\pi (\theta )$ directly, but if we build a chain ${\theta }_0,{\theta }_1,\dots$ with transition kernel $T(\theta \to {\theta }^{\prime })$ satisfying the detailed-balance (reversibility) condition with target $p$,

$$p(\theta )T(\theta \to {\theta }^{\prime })=p({\theta }^{\prime })T({\theta }^{\prime }\to \theta ),$$

then $p$ is stationary for the chain: if ${\theta }_t\sim p$ then ${\theta }_{t+1}\sim p$. Running the chain long enough (after a burn-in $m$) yields approximately iid draws from $p$, and the ergodic theorem replaces the intractable integral with a time average:

$$\hat{\mathbb{E}}[h(\theta )]=\frac{1}{N-m}\sum _{t=m+1}^{N}h({\theta }_t)\stackrel{a.s.}{\to }\int h(\theta )p(\theta \mid x)d\theta .$$

Metropolis–Hastings (Tsay §12.4.1–12.4.2). Propose ${\theta }^{\ast }$ from a proposal kernel $q({\theta }^{\ast }\mid {\theta }_{t-1})$, then accept with probability

$$\alpha =\min \left(1,\frac{p({\theta }^{\ast }\mid x)}{p({\theta }_{t-1}\mid x)}\cdot \frac{q({\theta }_{t-1}\mid {\theta }^{\ast })}{q({\theta }^{\ast }\mid {\theta }_{t-1})}\right).$$

• Symmetric proposal ($q$ symmetric, e.g. random walk ${\theta }^{\ast }={\theta }_{t-1}+ϵ$, $ϵ\sim N(0,s^2)$): the ratio collapses to $\alpha =\min (1,p({\theta }^{\ast })/p({\theta }_{t-1}))$ - this is Metropolis.
• Asymmetric proposal: the Hastings correction $q({\theta }_{t-1}\mid {\theta }^{\ast })/q({\theta }^{\ast }\mid {\theta }_{t-1})$ restores detailed balance; e.g. independence MH or proposals on a constrained space.
• The normalizing constant cancels, because $p$ appears as a ratio - this is the whole reason MCMC exists.

Gibbs sampling (Tsay §12.2). Partition $\theta =({\theta }_1,\dots ,{\theta }_k)$. Iterate: for each block $i$, draw

$${\theta }_i^{(t)}\sim p\left({\theta }_i\right|{\theta }_1^{(t)},\dots ,{\theta }_{i-1}^{(t)},{\theta }_{i+1}^{(t-1)},\dots ,{\theta }_k^{(t-1)},x).$$

Each full conditional is itself a posterior at the block level and is often a conjugate family (Normal, Inverse-Gamma, Beta, Dirichlet, …). Gibbs is MH with the full conditional as the proposal, and the acceptance probability is exactly 1 - no rejections, no tuning. After burn-in $m$, the mean is Tsay eq. (12.3):

$${\bar{\theta }}_i=\frac{1}{N-m}\sum _{t=m+1}^N{\theta }_i^{(t)}.$$

Gibbs needs the full conditionals to be samplable; when one is not (e.g. an MA or GARCH coefficient), one replaces that step with a Metropolis or Griddy Gibbs update (evaluate the univariate conditional on a grid, invert the CDF numerically - Tsay §12.4.3). MCMC for latent variables is often written as data augmentation (Tanner–Wong), with EM as its deterministic special case.

Convergence diagnostics (BDA3 Ch 11). MCMC gives correlated draws, so the naive standard error $\sigma /\sqrt{N}$ is wrong. Two standard checks:

• Gelman–Rubin: run $S$ chains from overdispersed starts; compare between-chain variance $B$ and within-chain variance $W$:

$$\hat{R}=\sqrt{\frac{\frac{N-1}{N}W+\frac{1}{N}B}{W}},\text{converged when }\hat{R}\approx 1.0(\text{e.g.}<1.01).$$

• Effective sample size from the autocorrelation ${\rho }_k$: $\mathrm{ESS}=N/(1+2{\sum }_{k\ge 1}{\rho }_k)$ - the number of independent draws the chain is worth. Report $N/\mathrm{ESS}$, not $N$.

Forward-filtering backward-sampling (FFBS; Tsay §12.8). For state-space / stochastic-volatility models the state path $z_{1:n}$ is drawn jointly (a Gibbs step over the whole latent path):

1. Forward: run the Kalman filter to get filtered moments $p(z_t\mid F_t)$ for all $t$.
2. Backward: sample $z_n\sim p(z_n\mid F_n)$, then for $t=n-1,\dots ,1$ sample from the Markov property $p(z_t\mid z_{t+1},F_n)=p(z_t\mid z_{t+1},F_t)$.

Carter–Kohn (1994) / Frühwirth-Schnatter (1994). This is the Bayesian counterpart of the Kalman smoother (see Kalman Filtering) and the engine of Bayesian stochastic-volatility estimation.

---

3. Computational Implementation - MH and Gibbs, from scratch

Three experiments: (a) random-walk Metropolis on a bimodal target, swept over three proposal scales, showing the mixing/acceptance trade-off and a biased estimate when mixing fails; (b) Gibbs on a bivariate normal, checking the recovered correlation; (c) Gibbs on the normal model with unknown mean and variance (semi-conjugate prior), sampling an Inverse-Gamma by transform. Stdlib only.

What the output proves.

• PART A - the mixing trade-off is real. With a tiny step ($s=0.5$) the chain accepts 81% of proposals but cannot cross the gap between the two modes, so it reports $\bar{\theta }=-0.72$ (target $-0.50$) and understates the variance ($5.52$ vs $6.03$). With $s=1.0$ and $s=4.0$ it mixes across both modes and nails the truth ($-0.5057$, $-0.5005$). High acceptance is not convergence - a chain that never moves is perfectly “accepted.”
• PART B - Gibbs recovers the joint exactly: correlation $0.8007$ vs target $0.8$, variances $\approx 1$, means $\approx 0$. No tuning, no rejections.
• PART C - the semi-conjugate normal model yields posterior $\mathbb{E}[\mu ]=1.1373$ (shrunk below $\bar{y}=1.2333$ toward the prior mean 0) and $\mathbb{E}[{\sigma }^2]=0.5792$; the sampler handles the unknown variance that would break a fixed-$\sigma$ formula.

---

4. Failure Modes & First-Principles Breakdowns

1. Pseudo-convergence / poor mixing (the $s=0.5$ failure). A chain can sit in one mode for the whole run and report a confident, wrong answer. The diagnostic is not acceptance rate - it is multi-chain $\hat{R}$ (start chains far apart), trace plots, and comparing to a longer run. The bimodal example above is the textbook case.
2. Forgetting burn-in (or using too little). Starting values are not draws from the posterior. Discard the burn-in and, better, check that estimates are stable across burn-in lengths. “Started at the MLE” is not an excuse - the chain needs to explore.
3. Autocorrelation invalidating $\sigma /\sqrt{N}$. Gibbs draws are serially correlated, so the effective sample size is far below $N$. Report ESS and inflate intervals accordingly; naive intervals are too narrow.
4. Reducible chains. If the proposal never reaches part of the support (or the chain is periodic), the stationary distribution is not the target, and the ergodic theorem does not apply. Verify detailed balance is satisfied, and use enough chains to detect it.
5. Bad proposal scaling in high dimensions. A random-walk MH needs its proposal scale to decay as $O(d^{-1/2})$ (optimal acceptance $\approx 0.234$, Roberts–Gelman–Gilks 1997); a badly chosen $s$ gives near-zero or near-one acceptance in high dimensions. Gibbs or HMC replaces blind proposals with dimension-aware moves.
6. Gibbs fails when full conditionals are not available or are nearly deterministic. Strongly correlated parameters (e.g. $\mu$ and ${\sigma }^2$ in the normal model, or hierarchical centering) make Gibbs crawl; the funnel geometry of hierarchical models is the classic offender (§06). Reparameterize (non-centered) or use HMC.
7. Confusing the posterior with the MLE. MCMC targets the posterior, not the maximum-likelihood point. If you want MLE, use an optimizer (see Numerical Optimization); MCMC gives the whole distribution and the MLE is a by-product.

---

5. References

• Tsay, Ruey S., Analysis of Financial Time Series (3rd ed.), Ch 12
• Robert, C. P. & Casella, G., Monte Carlo Statistical Methods (2nd ed.)
• Gelman et al., Bayesian Data Analysis (3rd ed.)
• Hastie, Tibshirani & Friedman, ESL
• McElreath, Statistical Rethinking (2nd ed.)
• Glasserman, Monte Carlo Methods in Financial Engineering

---

6. Connected Graph Bridges

• Back: 03 · Posterior Inference · 04 · Bayesian & Regularization · Index Hub
• Continue: 06 · Advanced Extensions (hierarchical models where MCMC is mandatory)
• Base: Numerical Methods · Monte Carlo (sampling, MC error) · Probability & Measure Theory (Markov chains, stationarity) · Econometrics & Time Series (state-space models, Kalman filter - the forward half of FFBS)
• Forward: GMM) (Baum–Welch is EM, Bayesian HMM is Gibbs/FFBS) · Kalman Filtering · Interest-Rate Models (Bayesian estimation of latent-factor term-structure models)