MCMC Maximum Likelihood

Definition

MCMC (Markov Chain Monte Carlo) Maximum Likelihood is a simulation-based method for computing the frequentist maximum likelihood estimator (MLE) in latent state models where $L(\theta) = \int p(Y|X,\theta)p(X|\theta)\,dX$ has no closed form. $J$ independent copies of the latent state path are augmented into the Gibbs target; raising the joint density to the $J$ -th power concentrates the marginal over $\theta$ at the MLE as $J\to\infty$ , without requiring gradient methods or direct likelihood evaluation.

Key Ideas

The augmented target is $p_J^m(\theta, \tilde{X}^J) \propto \bigl[\prod_{j=1}^J p(Y|\theta,X^j)p(X^j|\theta)\bigr]\cdot m(\theta)$ , where $m(\theta)$ is a dominating measure. The $\theta$ -marginal satisfies $p_J(\theta) \propto L(\theta)^J\cdot m(\theta)/m_J$ — a power-raised likelihood, as in simulated annealing (SA).
The SA connection: as $J\to\infty$ with $m(\theta)$ smooth and positive at $\hat\theta$ , $p_J(\theta)$ concentrates at the MLE. Unlike SA, $J$ is a fixed integer (not a temperature schedule), and the likelihood need not be directly computable — only samples from $p(X|\theta,Y)$ are needed.
Special cases:
- $J=1$ , $m(\theta)=p(\theta)$ : standard Bayesian MCMC (posterior distribution).
- $J=1$ , $m(\theta)\propto 1$ : may fail to integrate (improper posterior).
- $J$ large: $\theta^{(g)} \to$ MLE.
CPU time is linear in $J$ and in the number of Gibbs draws $G$ — identical order to Bayesian MCMC at the same $J$ .

How It Works

Algorithm

At each Gibbs step $g$ :

Draw $J$ latent state copies: $X^{j,(g+1)} \mid \theta^{(g)}, Y \sim p(X^j \mid \theta^{(g)}, Y)$ , independently for $j=1,\ldots,J$ . These are smoothed draws — conditioned on the full data $Y$ . Metropolis-Hastings (MH) steps used when direct sampling is infeasible.
Draw parameter: $\theta^{(g+1)} \mid \tilde{X}^{J,(g+1)}, Y \sim p_J^m(\theta \mid \tilde{X}^J, Y) \propto \bigl[\prod_j p(Y|\theta,X^{j,(g+1)})p(X^{j,(g+1)}|\theta)\bigr]\cdot m(\theta)$ . For conjugate models, stacking the $J$ copies yields a regression with $JT$ observations: exact Normal/Inverse-Gamma (IG) draws.

Convergence Theorem (Jacquier-Johannes-Polson [JJP] 2007)

Under mild regularity conditions on $L(\theta)$ and $m(\theta)$ near $\hat\theta$ :

$\sqrt{J}\bigl(\theta^{(g)} - \hat\theta\bigr) \;\Rightarrow\; \mathcal{N}\!\bigl(0,\; [-\ell''(\hat\theta)]^{-1}\bigr)$

where $-\ell''(\hat\theta)$ is the observed Fisher information. Two practical consequences:

Standard errors: Var of draws $\approx [-\ell''(\hat\theta)]^{-1}/J$ ; multiply draws by $\sqrt{J}$ and compute empirical variance.
Convergence diagnostic: normality of $\sqrt{J}(\theta^{(g)} - \hat\theta)$ — checked via Jarque-Bera or normal probability plot — signals whether $J$ is large enough, without knowing the true $\hat\theta$ .

Dominating Measure $m(\theta)$

Required when $L(\theta)^J$ fails to integrate for any $J$ . Examples:

Some Cauchy-type models: likelihood integrates only at $J=0$ .
Jump-diffusion models: $\sigma_v\to 0$ drives $L\to\infty$ ; $m(\theta)$ bounds $\sigma_v$ away from zero.
Two-factor stochastic volatility (SV) model with latent variance: $L(t,\sigma) \sim t^{-1}$ in the right tail for any fixed $J=1$ ; raising to $J=2$ integrates.

Often $m(\theta)\propto 1$ works for $J\geq 2$ , making the method prior-free.

Empirical Performance (JJP 2007)

Model	$J$	$G$	CPU	Variance ratio ( $J$ : $1\to 20$ )
SV: $\log V_t = a+d\cdot\log V_{t-1}+\sigma_v v_t$	20	25,000	$\approx 20$ min (Sun workstation)	$1/20$ for $d$
multivariate normal (MVN) jump-diffusion ( $K=3$ , Merton 1976)	20	5,000	$\approx 30\text{s} \times J$ per draw	$4.4\times$ reduction for $\lambda$

Smoothed volatility estimates at the MLE are nearly identical to Bayesian smoothed estimates: the MCMC smoother is preserved regardless of $J$ .

Why It Matters

Provides gradient-free MLE in models where gradient methods face local maxima, boundary issues, or expensive numerical Hessians.
Delivers the optimal smoothing distribution of latent states at the MLE — unlike approximations that substitute point estimates into a Kalman filter.
Unifies Bayesian and frequentist inference in a single $J$ -parameterized family, enabling sensitivity analysis of how posterior shape changes as prior information is removed.

Open Questions

Choice of $J$ is problem-specific; the normality diagnostic provides a practical rule but no closed-form formula.
Extension to boundary MLEs or non-stationary latent processes requires additional regularity conditions.
Cost scales linearly in $J$ ; very large state spaces may require impractically large $J$ for tight concentration.