Marginal Data Density

Definition

The Marginal Data Density (MDD) — also called the marginal likelihood or model evidence — is the probability of observing the data under a given model, with parameters integrated out. It is the canonical Bayesian criterion for comparing models of any type (nested or non-nested).

Science as data reduction (Sims 1996). A model that compresses data without sacrificing predictive content is a good model; the MDD is the formal measure of this compression. Just as Newton's mechanics compressed Kepler's laws by subsuming them under fewer principles, a model with high MDD achieves the same explanatory content with less parameter freedom. This framing motivates MDD as the epistemologically correct criterion for Bayesian model evaluation in macroeconomics, where a single historically given time series — not repeated experiments — is the object of inference. See Sims (1996).

Key Ideas

• The MDD penalizes overfitting automatically: a model that achieves a high likelihood only in a narrow region of parameter space is penalised by the prior volume.
• For unrestricted reduced-form vector autoregressions (VARs), the MDD has a closed-form expression.
• For restricted (identified or over-identified) structural VARs, it requires numerical approximation via the Chib (1995) Gibbs-sampler algorithm.
• MDD comparisons between VARs and Dynamic Stochastic General Equilibrium (DSGE) models are the standard for evaluating whether theoretical restrictions compromise fit.

How It Works

Definition

$$p(Y_T \mid \mathcal{M}) = \int_\Theta L(\phi \mid Y_T)\, p(\phi \mid \mathcal{M})\, d\phi$$

where $\phi$ collects all model parameters, $\Theta$ is their domain, $L(\phi \mid Y_T)$ is the likelihood, and $p(\phi \mid \mathcal{M})$ is the prior under model $\mathcal{M}$.

Bayes factor between two models $\mathcal{M}_1$ and $\mathcal{M}_2$:

$$\text{BF}_{12} = \frac{p(Y_T \mid \mathcal{M}_1)}{p(Y_T \mid \mathcal{M}_2)}$$

$\text{BF}_{12} > 1$ favors $\mathcal{M}_1$. Because the integral over $\Theta$ is typically evaluated in log space, the comparison is often reported as $\log p(Y_T \mid \mathcal{M}_1) - \log p(Y_T \mid \mathcal{M}_2)$.

Closed Form for Unrestricted Reduced-Form VARs

With the Sims-Zha Prior, the reduced-form VAR has a Normal-Wishart posterior structure. The MDD has a closed-form expression involving the ratio of normalizing constants:

$$\log p(Y_T) = \text{const} + \frac{T \cdot n}{2} \log(2\pi) + \log \frac{\Gamma_n(\nu_T/2)}{\Gamma_n(\nu_0/2)} + \frac{\nu_0}{2}\log|\Psi_0| - \frac{\nu_T}{2}\log|\Psi_T| + \frac{1}{2}\log\frac{|V_0|}{|V_T|}$$

where $(\nu_0, \Psi_0, V_0)$ are prior parameters and $(\nu_T, \Psi_T, V_T)$ are posterior parameters. No Markov Chain Monte Carlo (MCMC) is needed.

Chib (1995) Algorithm for Restricted VARs

For over- or exactly-identified structural VARs with linear restrictions $Q_j a_j = 0$, $R_j f_j = 0$, no closed form exists. The Chib algorithm exploits the identity:

$$\log p(Y_T) = \log L(\phi^* \mid Y_T) + \log p(\phi^*) - \log p(\phi^* \mid Y_T)$$

evaluated at any point $\phi^*$ (typically the posterior mode). The three terms on the right are:

1. Log-likelihood at $\phi^*$: direct evaluation.
2. Log prior at $\phi^*$: closed form given the Gaussian Sims-Zha Prior.
3. Log posterior at $\phi^*$: requires evaluation of $p(\phi^* \mid Y_T)$, which is itself approximated via the Gibbs sampler output using conditional posteriors.

Gibbs Sampler for Structural VARs (Waggoner-Zha 2003a)

Partition the free parameters into blocks $\{b_j, g_j\}_{j=1}^n$ where $a_j = U_j b_j$ and $f_j = V_j g_j$ (see Structural Identification). The sampler alternates:

1. Draw $g_j \mid b_j, \{b_i, g_i\}_{i \neq j}, Y_T$ — posterior of $g_j$ given $b_j$ is Gaussian (closed form).
2. Draw $b_j \mid \{b_i\}_{i \neq j}, \{g_i\}_{i=1}^n, Y_T$ — posterior of $b_j$ given all other blocks has a closed-form pdf and can be simulated exactly.

Because each conditional pdf is available in closed form, Chib's algorithm can evaluate $p(\phi^* \mid Y_T)$ accurately without long chains. In practice, a 10-variable, 13-lag VAR completes in under one minute.

Carlin-Chib (1995) Augmented Sampler

An alternative to computing $p(y|M_j)$ separately for each model is to run a single MCMC that visits all $K$ models simultaneously and reads off posterior model probabilities directly from visit frequencies.

The absorbing state problem. In a naive augmented sampler with model indicator $M \in \{1,\ldots,K\}$, parameters of inactive models are "dropped" at each iteration. When a parameter is forced out of the model, the chain enters an absorbing state (the dropped parameter can never be regenerated by a standard Gibbs step), violating the irreducibility condition for convergence.

The fix: pseudo-priors. Introduce pseudo-priors (linking densities) $p(\theta_j | M \neq j)$ — proper distributions for the parameters of model $j$ when $j$ is not the active model. Since $y \perp \theta_j | (M \neq j)$, these densities are irrelevant to $p(y|M=j)$ and can be chosen freely. With pseudo-priors, the complete joint model is:

$$p(y, \theta_1,\ldots,\theta_K, M=j) = f(y|\theta_j)\,p(\theta_j|M=j)\,\prod_{i \neq j} p(\theta_i|M \neq i)\,\pi_j$$

All full conditionals are now well-defined:

• $\theta_j | M=j, y$: standard model-$j$ posterior; when $M \neq j$, draw from pseudo-prior $p(\theta_j|M\neq j)$.
• $M | \theta, y$: discrete distribution proportional to $f(y|\theta_j)\,p(\theta_j|M=j)\,\pi_j$.

No absorbing states exist; standard Gibbs convergence conditions apply.

Bayes factor estimation. At convergence, visit frequencies estimate posterior model probabilities:

$$\hat{P}(M=j|y) = G_j / G, \qquad B_{12} = \frac{\hat{P}(M=1|y)}{\hat{P}(M=2|y)} \cdot \frac{\pi_2}{\pi_1}$$

Standard errors follow from batching the autocorrelated $M$-chain.

Pseudo-prior tuning. Mixing depends on how closely pseudo-priors approximate the model-specific posteriors $p(\theta_j|M=j, y)$. Recommended workflow: run $K$ individual model-specific chains first; use first-order normal approximations to these posteriors as pseudo-priors. Diffuse pseudo-priors produce competitive proposals only rarely → high autocorrelation in $M$-chain → slow convergence. Prior weights $\pi_j$ can be tuned to rebalance a badly skewed chain; the estimated Bayes factor adjusts accordingly.

Contrast with Chib (1995). Chib computes $\log p(y|M_j)$ via an identity evaluated at a single point from one model's MCMC output; Carlin-Chib runs a combined sampler spanning all models. Both target the same quantity. Carlin-Chib is more direct for small $K$ (2–5 models); Chib scales better to large $K$ and avoids pseudo-prior design. For the cointegrated vector error correction model (VECM) application of Warne (2006), the Chib identity is preferred because $K$ is the full $(r, k)$ grid.

Sensitivity note. When the likelihood does not identify a model's parameters well (e.g., a 4-component mixture where components 2 and 3 nearly overlap), Bayes factors are sensitive to the prior. In the galaxy example (Carlin-Chib 1995), replacing Dirichlet(1,1,1,1.5) with Dirichlet(3,8,3,3) shifts $\hat B_{12}$ from 0.57 to 284 — a reversal of the model ranking driven entirely by prior specification.

Lewis-Raftery (1997): Laplace-Metropolis Estimator

Lewis and Raftery (1997) apply the Laplace approximation to the marginal likelihood integral using quantities estimated from existing MCMC output — no additional simulation runs required.

Basic Laplace-Metropolis (LM) estimator (fixed-effects models, eq. 4): $$\log\hat{f}(Y) = \log f(Y|\theta^*) + \log f(\theta^*) + \frac{P}{2}\log(2\pi) + \frac{1}{2}\log|H^*|$$

• Mode $\theta^*$: the multivariate median ($L_1$ center) of the MCMC draws — robust to the occasional distant excursions that MCMC chains exhibit, where the posterior mean and direct maximization can fail.
• Variance $H^*$: the minimum volume ellipsoid weighted variance matrix (Rousseeuw-van Zomeren 1990) — also robust to outliers.

Compound LM estimator (hierarchical/random-effects models): For models with nuisance random effects $a_i$ per subject, marginalizing directly over all $a_i$ is infeasible. The compound LM uses a nested Laplace approximation (inner Laplace over each $a_i$ at its Newton-step conditional mode $a_i^*$, outer LM over the remaining parameters $(\eta, \xi)$). For the logistic hierarchical model this is fully closed-form.

Uncertainty via batch means: Divide MCMC output into $B$ batches; compute a separate LM estimate per batch; treat the $B$ estimates as approximately independent to form a 95% interval.

Validation — Iran Fertility Survey (29 women, mixed logistic): compound LM gives $\log\hat{f}\approx-220.5$ against gold-standard $-221.8$ (50 million prior draws); 95% credible interval (CI) $(−222.0,−219.0)$ covers the true value.

Practical position: requires no additional MCMC runs (unlike Chib 1995 and Chib-Jeliazkov [CJ] 2001); more stable than the harmonic mean estimator; extendable to hierarchical models (unlike basic Chib). Main limitation: assumes approximate posterior unimodality — fails for multimodal or heavy-tailed posteriors. See Lewis-Raftery (1997).

Alternative Marginal Likelihood (ML) Estimators (Amisano-Federico 2004)

When $\Omega$ is high-dimensional, Chib's method requires $(q-1)$ separate Gibbs runs (one per parameter block) and can return numerically unstable results. Two alternatives are:

Gelfand-Dey (1994) modified importance sampling: $$m(y_T \mid \mathcal{M})^{-1} \approx \frac{1}{N}\sum_{i=1}^N \frac{f(\theta^{(i)})}{p(\theta^{(i)} \mid \mathcal{M})\, \ell(\theta^{(i)}; y_T)}$$

where $f(\theta)$ is a truncated multivariate normal integrating to one over the parameter space (chosen to have thinner tails than the posterior). The estimate is a single pass over the Gibbs output. More numerically stable than Chib but slightly less precise.

Laplace method: $$\log m(y_T \mid \mathcal{M}) \approx \frac{k}{2}\log(2\pi) + \log|\tilde{\Sigma}|^{1/2} + \log L(\tilde\theta; y_T) + \log p(\tilde\theta \mid \mathcal{M})$$

where $\tilde\theta$ is the posterior mode and $\tilde\Sigma$ the inverse negative Hessian (approximated by the posterior mean and covariance from the Gibbs output). Requires no additional Gibbs runs and agrees closely with Chib for low-dimensional models; precision degrades as the parameter space grows.

Savage-Dickey density ratio: For a nested comparison $\mathcal{M}_0: \theta_1 = \theta_{1,0}$ vs. unrestricted $\mathcal{M}$: $$\text{BF}(\mathcal{M}_0, \mathcal{M}) = \frac{p(\theta_{1,0} \mid y_T, \mathcal{M})}{p(\theta_{1,0} \mid \mathcal{M})}$$

The prior density in the denominator is analytic; the posterior density in the numerator is estimated by Rao-Blackwellization from the Gibbs output. Exact for the nested case and requires no additional simulation.

Practical ranking (Amisano-Federico 2004): Chib $\approx$ Laplace $>$ Gelfand-Dey for low-dimensional models; all three break down (diverge) when $p(\theta^* \mid y_T, \mathcal{M})$ concentrates in a region the Gibbs sampler visits rarely — which happens for the Kronecker model M3.

Chib-Greenberg (1998): Kernel Density When Normalising Constants Are Unknown

When a full conditional density $\pi(\sigma|Z,\beta)$ has an intractable normalising constant — as in the multivariate probit where $\sigma$ lies in the convex set $C$ of valid correlation matrices — the Rao-Blackwell estimator of Chib (1995) is unavailable for that block. Chib and Greenberg (1998) estimate $\pi(\sigma^*|y)$ via kernel density smoothing (eq. 10):

$$\hat\pi(\sigma^*|y) = G^{-1}\sum_{g=1}^G \prod_{j=1}^p b_j^{-1}K\bigl\{b_j^{-1}(\sigma_j^* - \sigma_j^{(g)})\bigr\}$$

where $b_j = s_j G^{-1/(p+4)}$ (Scott 1992) and $K = \phi$ (standard Normal kernel). For $p = J(J-1)/2$ large (e.g., $p=21$ for $J=7$), the ordinate is factored into blocks and each estimated sequentially from reduced runs: $\pi(\sigma^*|y) = \pi_1(\sigma_1^*|y)\pi_2(\sigma_2^*|y,\sigma_1^*)\cdots$. The normalising constant of the prior $\pi(\sigma^*)$ is found by MC proportion of draws from $\phi_p(\sigma|\sigma_0,G_0^{-1})$ satisfying positive definiteness.

Limitation. Kernel density estimation is biased in high dimensions and yields numerical standard error (NSE) an order of magnitude larger than Rao-Blackwellisation. This motivated Chib-Jeliazkov (2001), which replaces kernel density with a two-expectation estimator derived from the detailed balance identity of the Metropolis-Hastings (M-H) kernel — applicable whenever M-H is used and requiring no tractable normalising constant.

Chib-Jeliazkov (2001): Extension to M-H Output

Chib and Jeliazkov (2001) resolve the central limitation of Chib (1995) — which fails when any parameter block uses a M-H step because the Rao-Blackwell estimator requires tractable full conditionals. The fix exploits the local reversibility of every M-H subkernel.

Local reversibility identity (eq. 7): for any M-H kernel with proposal $q$ and acceptance probability $\alpha$,

$$p(\theta, \theta^*|y)\,\pi(\theta|y) = \pi(\theta^*|y)\,p(\theta^*, \theta|y)$$

Integrating over $\theta$ and rearranging gives the CJ estimator for the posterior ordinate at the fixed point $\theta^*$:

$$\hat\pi(\theta^*|y) = \frac{E_1\bigl\{\alpha(\theta, \theta^*|y)\,q(\theta, \theta^*|y)\bigr\}}{E_2\bigl\{\alpha(\theta^*, \theta|y)\bigr\}}$$

• $E_1$: expectation over the posterior $\pi(\theta|y)$ — estimated by averaging $\alpha(\theta^{(g)}, \theta^*|y)\,q(\theta^{(g)}, \theta^*|y)$ over the $G$ draws from the main MCMC run (no extra simulation).
• $E_2$: expectation over the proposal $q(\theta^*, \cdot|y)$ at fixed $\theta^*$ — estimated by drawing $J$ proposal values $\tilde\theta^{(j)} \sim q(\theta^*, \cdot|y)$ in a short reduced run and averaging $\alpha(\theta^*, \tilde\theta^{(j)}|y)$. Cost: $J$ likelihood evaluations.

Log marginal likelihood (eq. 10):

$$\log\hat{m}(y) = \log f(y|\theta^*) + \log\pi(\theta^*) - \log\hat\pi(\theta^*|y)$$

Multi-block with latent variables (eq. 11–15): For a two-block sampler with latent data $z$ where block 1 uses M-H and block 2 uses Gibbs, the reduced conditional ordinate for block 1 is:

$$\hat\pi(\theta_1^*|y) = \frac{E_1\bigl\{\alpha(\theta_1, \theta_1^*|y,\theta_2,z)\,q(\theta_1, \theta_1^*|y,\theta_2,z)\bigr\}}{E_2\bigl\{\alpha(\theta_1^*, \theta_1|y,\theta_2,z)\bigr\}}$$

where $E_1$ averages over $\{(\theta_1^{(g)}, \theta_2^{(g)}, z^{(g)})\}$ from the main run and $E_2$ averages $J$ proposal draws at $\theta_1^*$ holding $(\theta_2^{(g)}, z^{(g)})$ fixed. Placing M-H blocks first avoids M-H steps during the reduced run.

General $B$-block algorithm (eq. 16–18):

$$\log\hat{m}(y) = \log f(y|\theta^*) + \log\pi(\theta^*) - \sum_{i=1}^B \log\hat\pi\!\left(\theta_i^*\,\big|\,y,\theta_1^*,\ldots,\theta_{i-1}^*\right)$$

Each ordinate $\hat\pi(\theta_i^*|\cdot)$ is computed by Rao-Blackwell if block $i$ is Gibbs, or by the CJ two-expectation formula if block $i$ is M-H. The reduced run for block $i$ conditions on all previously estimated blocks $(\theta_1^*, \ldots, \theta_{i-1}^*)$.

NSE and efficiency (§2.4): Collect into the vector process $h^{(g)}$ all simulation quantities entering $\log\hat{m}(y)$. Estimate $\text{var}(\hat{h})$ by Newey-West with $m=40$ lags. NSE follows from the Delta method: $\sqrt{a'\,\widehat{\text{var}}(\hat{h})\,a}$. Validated by 50 independent replications. Lower MCMC inefficiency factors translate directly into smaller NSE.

Four empirical examples:

Example	Sampler	$\log\hat{m}$	NSE
Binary logit (Mroz 1987, $n=200$, 8 params)	Tailored M-H	$-144.756$	~0.03
Binary logit	Random-walk M-H	$-144.756$	~0.30 (10× larger)
AIDS CD4 longitudinal (467 subjects)	Multi-block Gibbs	$-3577.57$	small
AIDS CD4 longitudinal	Single-block M-H	$-3577.57$	larger
Poisson longitudinal (epilepsy, $n=58$)	CJ	$-915.23$	$\approx 0.1$
Poisson longitudinal	Kernel density	$-915.49$	$\approx 0.1$
Multivariate probit (PSID 520×7, 21 $\sigma$ params)	One $\sigma$-block	$-1550.533$	0.842
Multivariate probit	Four $\sigma$-blocks	$-1550.45$	0.278

Single likelihood evaluation for the multivariate probit ≈ 2 min at 550 MHz CPU, making Gelfand-Dey / Meng-Wong importance sampling prohibitively expensive.

Comparison to alternatives: Kernel density — biased in high dimensions; NSE order-of-magnitude larger than CJ. Gelfand-Dey / Meng-Wong — too expensive when likelihood evaluation is costly. Reversible-jump — essential when model dimension varies; unnecessary here. CJ is the preferred method for comparing fixed-dimensional models fit by any MCMC scheme.

Han-Carlin (2001): Comparative Benchmark of Six Methods

Han and Carlin (2001) evaluate six MCMC approaches to computing $B_{21} = p(y|M_2)/p(y|M_1)$ on three empirical datasets of increasing scale.

The six methods:

Method	Abbrev	Key idea
Carlin-Chib (1995) product space	CC	Single Gibbs chain over all $K$ model spaces with pseudo-priors; $M$ updated via full conditional (eq. 2)
Metropolized CC	MCC	Replace Gibbs step for $M$ with Metropolis step; draw only from proposed model's pseudo-prior (eq. 4)
Reversible Jump (Green 1995)	RJ	Dimension-matching bijection $\mathbf{g}_{j,j'}(\theta_j,u)$ with Jacobian; MH acceptance (eq. 5)
Partial Analytic Structure (Godsill 2001)	PAS	Integrate out analytically tractable subvectors within RJ acceptance ratio (eq. 6)
Chib (1995) Marginal Likelihood Identity (MLI)	Chib	$\log p(y)=\log f(y|\tilde\theta)+\log p(\tilde\theta)-\log p(\tilde\theta|y)$; posterior density by Rao-Blackwellization (eq. 7)
Chib-Jeliazkov (2001)	CJ	$p(\tilde\theta|y)=E_1\{\alpha q\}/E_2\{\alpha\}$ via two MH expectations; extends Chib to non-Gibbs samplers (eq. 11)

Empirical results:

• Radiata pine ($n=42$, linear regression, Table 2): all six methods give accurate estimates; true $B_{21}=4862$, estimates 4782–4943. No performance differences at this scale.
• AIDS CD4 counts (467 subjects, random-effects model): CC and MCC fail completely — after 1.5 million iterations the chain never mixes across models; RJ: $\hat{B}_{21}=67.83$; Chib and CJ succeed. Failure of CC/MCC is not a tuning problem; it is structural: the high-dimensional random-effect vectors make cross-model proposals essentially impossible under pseudo-priors.
• Prostate cancer binary ($n=53$, logistic regression with latent data): RJ with Albert-Chib latent variable augmentation gives $\hat{B}_{21}=6.67$ [6.61, 6.73].

Practical guidance: Marginal likelihood methods (Chib, CJ) are preferred for standard hierarchical models with fixed parameter dimension. Reversible jump is essential when parameter dimension differs across models. CC/MCC should not be applied to models with many random effects.

Chib MLI for Cointegrated VAR (Warne 2006)

The Chib (1995) identity takes a specific form when applied to the VECM at a fixed rank $r$. Evaluated at posterior estimates $(\tilde\alpha, \tilde\Psi)$ (typically the posterior mean from the Gibbs sampler):

$$\log p(Y|r) = \underbrace{\log p(Y|\tilde\alpha,\tilde\Psi,r)}_{\text{(i) analytic}} + \underbrace{\log p(\tilde\alpha,\tilde\Psi|r)}_{\text{(ii) analytic}} - \underbrace{\log p(\tilde\Psi|\tilde\alpha,Y,r)}_{\text{(iii) matrix-}t} - \underbrace{\log p(\tilde\alpha|Y,r)}_{\text{(iv) Rao-Blackwell}}$$

Each term:

(i) Log-likelihood. $p(Y|\tilde\alpha,\tilde\Psi,r) = \prod_t \mathcal{N}(\Delta x_t;\, \hat\mu_t,\, \tilde\Omega)$ evaluated at the Gibbs posterior mean of $\Omega$. Analytic.

(ii) Log prior. $\log p(\tilde\alpha,\tilde\Psi|r) = \log p(\tilde\Psi) + \log p(\tilde\alpha|\tilde\Psi,r)$. Both are analytic: $p(\tilde\Psi)$ is the matrix-$t$ Grassmann prior; $p(\tilde\alpha|\tilde\Psi,r)$ is Gaussian under the Villani/Warne prior.

(iii) Posterior of $\Psi$. $p(\tilde\Psi|\tilde\alpha,Y,r)$ is the full conditional from the Gibbs sampler, which is a matrix-$t$ distribution with known parameters. Analytic.

(iv) Posterior of $\alpha$. $p(\tilde\alpha|Y,r) = \int p(\tilde\alpha|\Psi,Y,r)\, p(\Psi|Y,r)\, d\Psi$ has no closed form. Estimated by Rao-Blackwellization over Gibbs draws: $\hat{p}(\tilde\alpha|Y,r) = G^{-1}\sum_{i=1}^G p(\tilde\alpha|\Psi^{(i)},Y,r)$, where each $p(\tilde\alpha|\Psi^{(i)},Y,r)$ is an analytic Gaussian evaluated at the draw $\Psi^{(i)}$.

The full posterior $p(r,k|Y)$ is then obtained by running the sampler for each $(r,k)$ combination and combining via Bayes' theorem with a (typically uniform) prior over $(r,k)$. The closed-form at full rank (Corollary 1 of Warne 2006) eliminates the need for MCMC when $r=n$, so MCMC cost is incurred only for $r < n$ candidate ranks.

Why It Matters

• MDD is the only model-comparison criterion that is simultaneously consistent, automatic in its complexity penalty, and applicable to non-nested models.
• Enables systematic evaluation of whether DSGE-implied restrictions degrade empirical fit (Smets-Wouters 2003; Del Negro-Schorfheide 2004).
• The Chib algorithm makes it feasible to compare large numbers of identified VARs with different restriction patterns.

Open Questions

• MDD is sensitive to the prior, especially its tails; poorly calibrated hyperparameters can reverse model rankings.
• For very large VAR systems, posterior sampling remains computationally demanding even with the Chib algorithm.

Related

• Vector Autoregression
• Time-Varying Parameter VAR
• Sims-Zha Prior
• Structural Identification
• Gibbs Sampler
• Gianni Amisano
• Tao Zha
• Daniel Waggoner
• Cointegration
• Anders Warne
• Siddhartha Chib
• Bradley P. Carlin
• Cong Han
• Ivan Jeliazkov
• Edward Greenberg
• Chib-Greenberg (1998)
• Lewis-Raftery (1997)
• Adrian E. Raftery