Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Identification–aware Markov Chain Monte Carlo
abstractLeaving posterior sensitivity concerns aside, non-identifiability of the parameters does not raise a difficulty for Bayesian inference as far as the posterior is proper, but multi-modality or flat regions of the posterior induced by the lack of identification leaves a challenge for modern Bayesian computation. Sampling methods often struggle with slow or non-convergence when dealing with multiple modes or flat regions of the target distributions. This paper develops a novel Markov chain Monte Carlo (MCMC) approach for non-identified models, leveraging the knowledge of observationally equivalent sets of parameters, and highlights an important role that identification plays in modern Bayesian analysis.We show that our identification-aware proposal eliminates mode entrapment, achieving a convergence rate uniformly bounded away from zero, in sharp contrast to the exponentially decaying rates characterizing standard Random Walk Metropolis and Hamiltonian Monte Carlo. Simulation studies show its superior performance compared to other popular computational methods including Hamiltonian Monte Carlo and sequential Monte Carlo. We also demonstrate that our method uncovers non-trivial modes in the target distribution in a structural vector moving-average (SVMA) application.
spacing{1.5}
\section{Introduction}
Many Bayesian economists and statisticians are inclined to adopt lindley1972bayesian’s assertion that “unidentifiability causes no real difficulty in the Bayesian approach,” a view further supported by gelman2014identifiability, who argues in his blog that “the concept of identification is less important in the Bayesian world than elsewhere.” In theory, if a well‐defined prior and likelihood yield a proper posterior, Bayesian inference can proceed without fundamental obstruction. However, in practice, modern applications often involve structural models with large parameter spaces, incomplete data, or weak identifying assumptions, leading to multi‐modal or flat posterior regions.
Applied economists often proceed under the assumption that a chosen sampler, such as Metropolis–Hastings (MH), Hamiltonian Monte Carlo (HMC) or sequential Monte Carlo (SMC), will eventually reveal the relevant features of the posterior, as is guaranteed asymptotically under certain strong regularity conditions.
However, real data sets are finite, computational budgets are limited, and posteriors are often analytically intractable, so the conditions required for fast convergence are difficult to verify and frequently fail. As a result, even a “proper” posterior can pose significant practical challenges for effective sampling.
Despite their widespread use in econometrics, popular algorithms such as MH, HMC, and SMC\footnote{Although SMC is not technically an MCMC method, the implementation by herbst2014sequential applies a local MCMC step to each particle, integrating MCMC techniques into the SMC framework.} often struggle to explore complex posterior surfaces thoroughly.
These methods can fail to traverse low-density regions separating multiple modes, can suffer from slow mixing in high-dimensional spaces, provide little guidance in flat or weakly identified regions. Several strategies have been proposed to alleviate these issues, including tempering methods that flatten the target distribution to facilitate movement between modes, and mode-jumping techniques that rely on either prior knowledge of mode locations or extensive presampling. Both approaches come with substantial costs: tempering often requires careful tuning of intermediate distributions and can remain computationally expensive in large-scale problems, while mode-jumping typically depends on substantial a priori information. Moreover, the challenges of multi-modality, high dimensionality, and flat or ill-conditioned likelihood regions compound one another, making it even harder to design robust, general-purpose samplers. As a result, determining how best to navigate unknown distributions that combine these features remains an open and active area of research.
This paper argues that identification analysis is crucial even for Bayesians, and shows that explicit knowledge of observationally equivalent parameter values can be used to overcome the challenges of sampling from complex multi-modal posteriors. We propose a new “identification–aware” MCMC approach that exploits the knowledge or computability of the set of observationally equivalent parameters and can attain faster convergence than widely used MCMC algorithms. The key innovation is a teleportation step: the chain moves directly to an observationally equivalent point in the state space, bypassing low-probability regions or distant valleys that incremental updates struggle to cross. By construction, this teleportation step preserves the target posterior while ensuring efficient exploration of flat regions, multiple modes, or disconnected neighborhoods. Our approach treats teleportation as a modular component that can be integrated into any Markov chain–based method, including MH, HMC, or the mutation steps in SMC. Our method combines teleportation with a standard MCMC-type transition (e.g., MH or HMC updates) to refine the position locally. This scheme leverages both directed global jumps and fine-tuned local adjustments, allowing the algorithm to balance global exploration with local optimization and to sample effectively from high-dimensional, multi-modal posteriors without requiring an exponentially increasing number of evaluations as dimension grows. The specific implementation depends on the structure of the observationally equivalent set, which may be finite, of varying cardinality, or a low-dimensional manifold. We provide theoretical guarantees that, in such settings, our identification–aware sampler improves convergence rates relative to standard MCMC approaches.
Simulation studies validate the superior performance of our proposal against popular computational methods. In multimodal mixture Gaussian settings, our algorithm successfully teleports across isolated modes where standard Random Walk Metropolis (RWM) consistently becomes trapped. Similarly, in set-identified scenarios, it achieves a more diffuse and representative exploration of the parameter space compared to SMC, which can exhibit clustering around initialization points. We further demonstrate the method's practical utility in Bayesian estimation for SVMA models as considered in plagborg2019bayesian. As shown in lippi1994var, SVMA models are not globally identified, with no closed-form expression for the impulse response identified-set. It makes the standard MCMC difficult to capture all the modes of the posterior distribution. In the complex environment with the number of parameters as many as 153, we demonstrate implementation of identification aware HMC algorithms and show that they can recover the dominant marginal modes found by the standard HMC, but also uncovers non-trivial high-density regions that remain elusive to the standard HMC. This suggests that incorporating teleportation step exploiting analytical characterization of identification of the model can help characterize the posterior landscape and provide a more robust and reliable numerical approximation for it.
Related Literature: There is a well-established literature in economics applying MCMC methods for Bayesian estimation and inference, including chib1996markov, an2007bayesian, flury2011bayesian, herbst2014sequential, kline2016bayesian, plagborg2019bayesian, meager2022aggregating, and muller2024locally. Beyond the Bayesian framework, these methods are also employed in frequentist settings. For example, chernozhukov2003mcmc uses standard MCMC to obtain minimizers of criterion functions for point-identified models, and chen2018monte applies SMC to construct confidence sets for identified sets defined by either likelihood-based conditions or moment equalities and inequalities, following approaches similar to herbst2014sequential and durham2014adaptive.
Most of the literature mentioned above treats these computational methods primarily as tools for estimation and inference rather than focusing on their theoretical properties such as convergence rates, central limit theorems (CLT), or laws of large numbers (LLN). One notable exception is herbst2014sequential, which provides some formal results, demonstrating that LLN and CLT can hold under certain conditions in their proposed SMC.
Outside of economics, a substantial literature has established theoretical guarantees for MCMC algorithms, but these results often rest on restrictive conditions that rarely hold in the complex settings typical of economic models. For example, classical mixing and convergence analyses for RWM (e.g. mengersen1996rates, roberts1996geometric, roberts2004general) typically assume log-concave or strongly unimodal targets, compact support, and uniform minorization, conditions that do not capture high-dimensional, multi-modal, or weakly identified likelihoods. Although a few exceptions can tackle more complicated cases, including multimodality (e.g., guan2007small, mangoubi2018does), the general applicability of these results to empirical economic problems remains unknown.
Similar issues arise in the theoretical treatment of more sophisticated methods, such as SMC and HMC. While several works chopin2004central, del2012adaptive, huggins2015convergence, durmus2017convergence, durmus2020irreducibility, mangoubi2021mixing offer insights into convergence rates, central limit theorems, and stability, they too rely on assumptions such as smoothness (e.g. Lipschitz-continuous gradients or bounded higher-order derivatives) and well-behaved high-dimensional scaling (e.g. controlled deterioration of error rates with dimension), assumptions that often fail in the heavily parameterized, data-driven models typical in econometrics.
For advanced algorithms explicitly designed to tackle multimodality, such as tempering-based samplers woodard2009sufficient and mode-jumping methods zhou2011multi, pompe2020framework, the theoretical underpinnings are even sparser. woodard2009sufficient demonstrates that tempering can still mix slowly in certain multi-modal regimes. Meanwhile, mode-jumping approaches often rely on heuristics or prior knowledge, or on computationally expensive approximations of mode locations, providing little theoretical guarantee beyond basic ergodicity. As a result, despite significant progress in algorithm development and analysis, sampling efficiently and reliably from complex, unknown distributions, especially those with multiple modes, high dimensionality, or flat identified sets, remains a major open challenge.
To the best of our knowledge, only a handful of studies exploit identification information directly during sampling.
For Bayesian estimation of finite mixture models, fruhwirth2006finite and geweke2007interpretation handle multimodality of the posterior by permuting or augmenting the observationally equivalent parameter values in a standard MCMC run.
Our identification–aware MCMC proposal includes their augmentation trick as a special case, while ours can cover more general cases in which cardinality of the observationally equivalent parameter values is infinite or varies over the parameter space.
Bayesian estimation for non-identified models has been studied extensively in the literature of structural vector autoregressions (SVAR). uhlig2005effects and rubio2010structural propose a Gibbs sampling algorithm with uniform draws of a non-identified orthonormal matrix. The step of drawing an orthonormal matrix can be viewed as a certain teleportation step in our framework.
In the SVAR context, our proposal of identification–aware MCMC can accommodate a non-uniform prior for the orthonormal matrix implied by a prior on the structural parameters along the proposal of baumeister2015sign.
Our analytical results on fast convergence of the identification–aware MCMC to the target posterior highlight a benefit of having a step of drawing the orthonormal matrix instead of performing standard MCMC procedures directly on the structural parameter space.
For locally-identified SVAR, bacchiocchi2025locally considers a transportation step in their Bayesian approach, while they do not investigate the convergence rate for approximating the posterior distribution.
The identification–aware proposal is designed to improve sampling efficiency in Bayesian inference, and should be viewed as complementary to the Bayesian sensitivity analysis to set-identified models.
Relying on computability of the mapping between reduced-form and structural parameters, kline2016bayesian and giacomini2021robust study how to draw Bayesian inference for the identified set.
chen2018monte develops criterion-based quasi-Bayesian procedures for frequentist-valid inference in partially identified models.
Their implementation builds on MCMC draws from the quasi-posterior that can have multiple modes or flat regions due to set-identification.
Our identification–aware MCMC can contribute to their inference approach by better approximating the flat quasi-posterior.
The method discussed in this paper also shares key structural similarities with Algorithm 1 in kuang2025robust, while its focus is on Bayesian robustness and inference on the identified set.
\paragraph{Roadmap} The rest of the paper is organized as follows: In Section (ref), we review standard MCMC methods (including MH and HMC) and SMC. We then introduce our algorithm and compare it against RWM and HMC in both multimodal and set-identified settings. Section (ref) presents two main simulation exercises comparing the sampling efficiency of our method with standard RWM, HMC and SMC. Section (ref) presents Bayesian analyses of an MA(1) process, and a 3-variable, 17-lag structural vector moving average model (SVMA), showing that our method delivers superior performance in large models.
\section{Markov Chain Monte Carlo and Mixing Time}
\subsection{MCMC Methods}\
Markov chain Monte Carlo methods are a class of algorithms used to sample from a target probability distribution \(\pi\) over a state space \(\Omega\), which can be either finite or general (e.g., \(\mathbb{R}^d\)). The goal of MCMC is to construct a Markov chain with a stationary distribution equal to \(\pi\), ensuring that, as the chain evolves, its samples approximate the desired target distribution.
When \(\Omega\) is finite, the Markov chain can be described by a transition matrix \(P(x, y)\), where \(P(x, y)\) represents the probability of transitioning from state \(x\) to state \(y\) conditional on that the current state is at \(x\). In the case of a general state space, \(\Omega\) is often a subset of \(\mathbb{R}^d\), and the Markov chain is defined through a transition kernel \(P(x, \cdot)\), which specifies a conditional probability measure given the current state \(x\). Popular MCMC algorithms, such as MH and HMC, operate as Markov chains with transition kernels constructed to preserve the target distribution and satisfy the necessary conditions for convergence.
\paragraph{Metropolis–Hastings.}
In MH, one proposes candidate states from a distribution \(q(x'|x)\) and either accepts or rejects these proposals based on an acceptance probability that ensures detailed balance. A simple version of the MH procedure is shown in Algorithm (ref). MH is appealing for its generality: any proposal distribution \(q(\theta'|\theta)\) that adequately explores the parameter space produces a Markov chain with \(\pi(\theta)\) as its stationary distribution. Common proposals include $\delta$-ball random walks, Gaussian random walks, and adaptive schemes that adjust proposal scales to target appropriate acceptance rates. The efficiency of MH depends on balancing exploration against excessive rejections, with optimal asymptotic acceptance rates around 0.234 gelman1997weak.
\begin{algorithm}[MH]
\begin{enumerate}
• Given target distribution $\pi(\theta)$, proposal distribution $q(\theta'|\theta)$, set initial state $\theta_0$ and number of iterations $N$.
• For $t = 1$ to $N$,
\begin{enumerate}
• Sample $\theta' \sim q(\cdot|\theta_{t-1})$.
• Compute the acceptance probability:
\[
\alpha := \min\left\{1, \frac{\pi(\theta') q(\theta_{t-1}|\theta')}{\pi(\theta_{t-1}) q(\theta'|\theta_{t-1})}\right\}.
\]
• Accept $\theta_t=\theta'$ with probability $\alpha$, otherwise, $\theta_t=\theta_{t-1}$.
\end{enumerate}
• Return $\{\theta_t\}_{t=1}^N$.
\end{enumerate}
\end{algorithm}
An important special case of the MH framework is the Gibbs sampler, which samples directly from conditional distributions when available, thereby avoiding acceptance steps. Gibbs can be efficient in hierarchical or conjugate settings but may be infeasible when conditionals lack convenient closed forms.
\paragraph{Hamiltonian Monte Carlo.}
Despite its generality, MH often suffers from slow random-walk behavior in high-dimensional spaces. HMC addresses this limitation by incorporating gradient information to make longer, directed moves through the parameter space. By introducing auxiliary momentum variables and simulating Hamiltonian dynamics, HMC proposals reduce random walks and achieve better mixing, particularly for correlated or high-dimensional targets neal2011mcmc. A standard version of HMC is given in Algorithm (ref).
\begin{algorithm}[HMC]
\begin{enumerate}
• Input: target density $\pi(\theta) \propto e^{-U(\theta)}$, step size $\epsilon$, number of leapfrog steps $L$, initial state $\theta_0$, number of iterations $N$.
• For $t = 1,\dots,N$:
\begin{enumerate}
• Sample momentum $p^{(0)} \sim \mathcal{N}(0,I)$.
• Initialize $(\theta^{(0)}, p^{(0)}) = (\theta_{t-1}, p^{(0)})$.
• Perform $L$ leapfrog steps:
\begin{enumerate}
• Half-step momentum update:
\[
p^{(0)} \leftarrow p^{(0)} - \tfrac{\epsilon}{2} \nabla_\theta U(\theta^{(0)}).
\]
• For $i = 1,\dots,L$:
\begin{enumerate}
• Position update:
\[
\theta^{(i)} = \theta^{(i-1)} + \epsilon p^{(i-1)}.
\]
• If $i < L$, full momentum update:
\[
p^{(i)} = p^{(i-1)} - \epsilon \nabla_\theta U(\theta^{(i)}).
\]
\end{enumerate}
• Final half-step momentum update:
\[
p^{(L)} = p^{(L-1)} - \tfrac{\epsilon}{2} \nabla_\theta U(\theta^{(L)}).
\]
\end{enumerate}
• Denote the proposal as $(\theta^\ast, p^\ast) = (\theta^{(L)}, p^{(L)})$, compute the acceptance probability:
\[
\alpha = \min\Bigl\{1,
\exp\bigl(-U(\theta^\ast) - \tfrac{1}{2}\|p^\ast\|^2
+ U(\theta_{t-1}) + \tfrac{1}{2}\|p^{(0)}\|^2\bigr)\Bigr\}.
\] Set $\theta_t = \theta^\ast$ with probability $\alpha$, and $\theta_{t-1}$ otherwise.
\end{enumerate}
• Return $\{\theta_t\}_{t=1}^N$.
\end{enumerate}
\end{algorithm}
In practice, performance depends on the choice of step size \(\epsilon\) and the number of leapfrog steps \(L\), often tuned adaptively or with algorithms such as NUTS hoffman2014no. Well-tuned HMC typically achieves acceptance rates of 60–80% betancourt2017conceptual, offering substantial efficiency gains relative to random-walk MH.
Although both MH and HMC are powerful, they remain local in nature, relying on proposals that stay reasonably close to the current state. This can be problematic in multi-modal settings where local moves risk becoming trapped in a single mode. SMC methods, discussed in Appendix (ref), take a more global approach using populations of particles to explore the state space in parallel, thereby providing a complementary alternative to standard MCMC.
\subsection{ Identification–aware MCMC}
Denote the target distribution $\pi(\theta)$, such as a posterior distribution of the form $p(\theta | Y)$ in Bayesian analysis. We define an associated partition $K: \Omega \rightarrow 2^\Omega$ which maps each parameter $\theta\in \Omega$ to the set of parameters $\theta'$ that are observationally equivalent to $\theta$. Specifically, this equivalence relationship is characterized by
$
L(\cdot | \theta)=L(\cdot | \theta')
$
for all $\theta' \in K(\theta)$ and possible data outcomes \(y\), where \(L(\cdot | \theta)\) denotes the likelihood function of the observed data given \(\theta\).
In the case of a mixture of two Gaussian distributions, with parameters $\theta=\left(p, (1-p), \mu_1, \mu_2, \sigma^2_1, \sigma^2_2\right)$, each equivalence class $K(\theta)$ captures the inherent label-switching problem: swapping the component labels does not affect the likelihood. Thus, for this mixture Gaussian case, each $K(\theta)$ comprises exactly two elements: $\left(p, (1-p), \mu_1, \mu_2, \sigma^2_1, \sigma^2_2\right)$ and $\left((1-p), p, \mu_2, \mu_1, \sigma^2_2, \sigma^2_1\right)$. However, unlike the setting in geweke2007interpretation, our $K(\theta)$ is not necessarily a finite set, nor does it need to have a fixed number of elements for all $\theta$. This flexibility allows for more general equivalence relationships.
When the parameters $\theta$ are not identified, the shape of the posterior for $\theta$ over $K(\theta)$ is fully determined by the prior for $\theta$ along $K(\theta)$ since the likelihood is flat on $K(\theta)$. A uniform prior on $K(\theta)$ can be a benchmark analysis if the researcher wishes to summarize the shape of the likelihood or put an unbiased belief for $\theta$. An alternative scenario is to specify a non-uniform prior on $K(\theta)$ to reflect an available prior knowledge for $\theta$. Our identification–aware MCMC algorithms can accommodate either cases so that it is useful no matter whether the purpose of posterior analysis is to perform subjective Bayesian inference or summarize the shape of likelihood.
The teleport kernel is the posterior conditional law restricted to the observationally equivalent set $K(\theta)$:
\begin{equation}
T(\theta,A) = \frac{\int_{A \cap K(\theta)} \pi(u) \nu(du)}{\int_{K(\theta)} \pi(u) \nu(du)}, \qquad 0<\int_{K(\theta)} \pi(u) \nu(du)<\infty, \pi - a.e.
\end{equation}
The reference measure $\nu$ is the natural one for the geometry of $K(\theta)$: counting measure if $K(\theta)$ is finite, Lebesgue measure if it has full dimension, and $r$-dimensional Hausdorff measure if it is an $r$-dimensional manifold.
We can interpret the posterior $\pi$ restricted to $K(\theta)$ as the conditional prior over observationally equivalent parameters given they belongs to identified set $K(\theta)$. Since the likelihood $L(y |\cdot)$ is constant on $K(\theta)$, it is uniform only if this conditional prior density is constant along $K(\theta)$.
When $\pi$ is constant on each $K(\theta)$, $T(\theta,\cdot)$ reduces to: \begin{equation}
T(\theta, A) = \frac{\nu(A \cap K(\theta))}{\nu(K(\theta))}, \qquad 0<\nu(K(\theta))<\infty, \pi - a.e.
\end{equation}
Note that for the teleport kernel to preserve the target distribution, $K(\theta)$ in ((ref)) does not strictly need to be the identified set for $\theta$ and suffices to be any subset of the identified set of $\theta$. That is, the teleportation kernel $T$ need not be supported on the full observationally equivalent set $K(\theta)$ and one can restrict the teleport move to a known or computationally convenient subset of the identified set.
All invariance and reversibility statements given below hold for any target $\pi$. For quantitative mixing bounds we will later impose a mild regularity condition on the restriction of $\pi$ to each $K(\theta)$. See Assumption (ref)(ref) in Section (ref).
To motivate and illustrate our proposal, consider the following toy example with 4 states.
\begin{example}
Consider two parameters \(\theta_1\) and \(\theta_2\), each taking binary values 0 and 1. The parameters have the following joint posterior distribution (we obscure its dependence on data $Y$ for simplicity, and write it as $p(\theta_1, \theta_2)$):
\[
p(\theta_1, \theta_2) =
\begin{cases}
p_{00}, & (\theta_1, \theta_2) = (0, 0), \\
p_{01}, & (\theta_1, \theta_2) = (0, 1), \\
p_{10}, & (\theta_1, \theta_2) = (1, 0), \\
p_{11}, & (\theta_1, \theta_2) = (1, 1).
\end{cases}
\]
In this setup, assume the target distribution depends only on $\theta_1$ and $\theta_2$ through their difference $\theta_1-\theta_2$. Therefore, \((\theta_1, \theta_2) = (1, 1)\) is observationally equivalent to \((\theta_1, \theta_2) = (0, 0)\), that is, under uniform prior,
$$p_{00} = p_{11} \quad \text{for any realization of data.}$$
We construct a Markov transition matrix \(P\) based on Gibbs sampler over the four possible states \(\{(0,0), (0,1), (1,0), (1,1)\}\), listed in that order. Starting from \((0,0)\), we run 100,000 Gibbs iterations under the setting \(p_{11} = p_{00} = 0.49999\) and \(p_{10} = p_{01} = 0.00001\) to simulate an extremely bimodal scenario.\footnote{The equality
$p_{11} = p_{00}$
follows from recognized observational equivalences and is assumed to be known to the researcher. In contrast, the choice
$p_{10} = p_{01}$
is introduced purely for computational simplicity and does not stem from any observational equivalence. Breaking this equality will not affect our results.} As shown in Figure (ref), the sampler is stuck at \((0,0)\) for the entire run and never visits \((1,1)\).\footnote{Extending the chain length to 1,000,000 increases the ratio of visits to $(1,1)$ to 0.53, but it still fails to accurately reflect the true ratio implied by the target distribution. Additionally, this resulting ratio varies significantly across different simulations, indicating instability in Gibbs's performance.}
\end{example}
\begin{figure}[htbp]
\begin{subfigure}[t]{0.48\textwidth}
\caption{Gibbs-sampled distribution.}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\caption{ identification–aware Gibbs.}
\end{subfigure}
\caption{Comparison of Gibbs vs. identification–aware Gibbs, 100,000 samples.}
\end{figure}
A natural remedy is to make the Gibbs sampler identification–aware.
Before each round of update, we first teleport to an observationally equivalent point,
chosen uniformly from $K(\theta)$.
This ensures that, in Example (ref), the chain moves freely between $(0,0)$ and $(1,1)$
according to their posterior probabilities, rather than remaining stuck in whichever state it started.
The same idea extends beyond Gibbs to any MCMC method.
Let $P(\theta,\cdot)$ denote the transition kernel of a standard algorithm such as MH, HMC, or Gibbs.
We compose the teleportation kernel \(T(\theta,\cdot)\) in (ref) with \(P\) to obtain an identification–aware kernel \(\widetilde P = PT\).\footnote{We write kernel composition by
\[
(PT)(x,A):=\int P(y,A)\,T(x,dy),
\]
so \(PT=P \circ T \) and \(PTf=P(Tf)\). Order can also matter: in general \(PT\) and \(TP\) need not have the same mixing rate unless the kernels commute.}
Concretely, here is one approach we recommend for implementing identification–aware MH:
\begin{algorithm}[Identification–aware MH]
\begin{enumerate}
• Given target distribution \(\pi(\theta)\), identification structure \(K(\cdot)\), symmetric proposal distribution \(q(\theta'|\theta)\), set initial state \(\theta_0\) and number of iterations \(N\).
• For \(t = 1\) to \(N\):
\begin{enumerate}
• Draw $\theta' \sim T(\theta_{t-1},\cdot)$.
• Sample \(\theta'' \sim q(\cdot|\theta')\).
• The acceptance probability should be defined as:
$$
\alpha\left(\theta',\theta^{\prime \prime}\right)=\min \Bigg\{1, \frac{\pi(\theta'') }{\pi(\theta') }\Bigg\}
$$
where:
• Accept \(\theta_t = \theta''\) with probability \(\alpha\), otherwise \(\theta_t = \theta'\).
\end{enumerate}
• Return \(\{\theta_t\}_{t=1}^N\).
\end{enumerate}
\end{algorithm}
Our algorithm blends two elements: a global move that draws a point from the observationally-equivalent class $K(\theta)$, step 2 (a), and a standard local transition, steps 2 (b) - (d).
The teleport step in Algorithm (ref) can be introduced also to the HMC algorithms shown in Algorithm (ref).
Specifically, before step 2 in Algorithm (ref), we teleport $\theta_{t-1}$ within $K(\theta_{t-1})$ to obtain $\theta^{\prime}$ as done in step 2 of Algorithm (ref), and run Algorithm (ref) step 2 (a) onward with $\theta^{\prime}$ in place of $\theta_{t-1}$.
By construction, $ \widetilde{P}$ inherits the advantages of the original MCMC method, and it preserves the stationary distribution, formalized below.
\begin{proposition}[Stationary Distribution of the Composite Kernel]
Let $\pi$ be the target probability distribution. If the transition kernel $P(\theta, A)$ is $\pi$--invariant, then the teleportation kernel $T(\theta, A)$ is $\pi$-preserving (i.e., $\int T(\theta, A) \pi(d\theta) = \pi(A)$), and the composite kernel $\widetilde{P} = PT $ is also $\pi$--invariant.
\end{proposition}
However, the $\pi$-invariance in Proposition (ref) does not, by itself, imply detailed balance\footnote{A Markov kernel $P$ is reversible (or satisfies detailed balance) with respect to a measure $\pi$ if it fulfills the condition $\pi(d\theta) P(\theta, d\theta') = \pi(d\theta') P(\theta', d\theta)$.} for the composed kernel $\widetilde P = PT $. In this paper, reversibility is needed only as a technical device for deriving mixing-rate bounds. When reversibility is required, we replace $\widetilde P$ by the order-randomized envelope $\bar P=\tfrac12(PT )+\tfrac12(TP)$, which is $\pi$-reversible. See Proposition (ref) in Appendix (ref).
The same principle generates a family of identification-aware samplers. One may replace $T$ by a finite-group label-switch operator, embed $T$ in a single Metropolis step via an exact or pseudo-marginal mixture proposal, or post-process a stored chain with conditional draws on each $K(\theta)$ (as in geweke2007interpretation), which is particularly helpful in high dimensions when likelihood evaluation is costly. Another option is the convex mixture $\widetilde P_H=(1-\varepsilon)P+\varepsilon T$. Implementation details, advantages, and limitations of these variants are collected in Appendix (ref). Empirically, no single variant dominates across settings. In our simulations and applications we select among $\widetilde P$, the order–randomized $\bar P$, batch augmentation and the convex mixture $\widetilde P_H$ according to computational cost and problem structure. For theoretical mixing bounds we analyze $\bar P$, which is $\pi$–reversible by construction. The convergence arguments in Section (ref) extend to all reversible variants with at most constant–factor changes in the conductance and hence the spectral–gap bounds. See Proposition (ref) in Appendix (ref) for details.
From a computational perspective, the identification–aware samplers add a fixed cost. Each sample from $T(\theta,\cdot)$ must (i) identify the observationally-equivalent set \(K(\theta)\) and (ii) draw from the conditional distribution \(\pi(\cdot | K(\theta))\). When \(K(\theta)\) is obtained by simple algebra, such as label permutations, sign flips, or orthogonal rotations, the overhead is trivial relative to one likelihood call, yet it eliminates the spectral bottlenecks that stall purely local chains, a gain most visible for highly multi-modal or high-dimensional posteriors.
If $K(\theta)$ can be computed only through expensive numerics (e.g., solving a high–dimensional nonlinear system), that gain may be partly offset by run time. Two compromises are then useful. One is to run several inexpensive local updates per teleport: perform one $T$ move then apply $P$ for $m$ consecutive iterations. This periodic schedule is time–inhomogeneous but remains $\pi$–stationary because both $P$ and $T$ leave $\pi$ invariant.\footnote{When $T$ is much cheaper than $P$, one may instead apply $m$ successive teleports followed by one local update, i.e.\ $PT ^{m}$, optionally retaining the $m$ intermediate teleport states. This also preserves $\pi$. See the batch augmentation in Appendix (ref).}
Alternatively, the convex combination $\widetilde P_H = (1-\varepsilon)P + \varepsilon T$ keeps the kernel time–homogeneous.
If we further let $\varepsilon \approx \frac{1}{m+1}$, this matches the average teleport frequency of the periodic schedule while keeping a fixed kernel each iteration. In either case, they preserve improvements on mixing rate while holding down the cost of computing \(K(\theta)\).
The practical rule is therefore to engage more identification–aware moves whenever \(K(\theta)\) is algebraic and the likelihood dominates computation. Otherwise, adopt a sparse or hybrid schedule that balances the mixing benefit against the price of constructing \(K(\theta)\).
\addtocounter{example}{-1}
\begin{example}[Continued]
Define an identification–aware transition kernel \(\widetilde{P} = T \cdot P\), where
$$T = \begin{pmatrix}
1/2 & 0 & 0 & 1/2\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
1/2 & 0 & 0 & 1/2
\end{pmatrix},$$
and \(P\) is the original Gibbs transition. As shown in Figure (ref), the identification–aware Markov chain mixes much better.
\end{example}
\subsection{Convergence and Mixing rate}
In general, we are interested in multiple aspects of the performance of these samplers.
First, we would like to have convergence between the sample distribution and target distribution, either in total variation distance, Kullback-Leibler divergence, or other discrepancy measures, which may or may not depend on the initial point/distribution. Then, we would like to understand how fast they converge.
Although other asymptotic properties such as the law of large numbers (LLN) and central limit theorems (CLT) can be of primary interest in their own right, particularly because most applications focus on the performance of a mean estimator of the form \(\frac{1}{N}\sum_{i=1}^N f(\theta_i)\) for a given function \(f\), these results are consequences of the convergence of the sample distribution.
In this paper, we measure convergence with the total–variation norm\footnote{If the chain has a unique stationary law $\pi$ and
$P^n(\theta,\cdot)\to\pi$ in total variation for every start $\theta$, then the ergodic theorem gives: for any $f\in L^1(\pi)$, the sample mean
$\bar f_N := \frac{1}{N}\sum_{t=1}^N f(\theta_t)$
satisfies $\bar f_N \to \mathbb{E}_\pi[f]$ almost surely. If, in addition, the chain is geometrically ergodic (defined below) and $f\in L^{2+\delta}(\pi)$ for some $\delta>0$, then a CLT holds:
$\sqrt{N} (\bar f_N-\mathbb{E}_\pi[f])\Rightarrow \mathcal N(0,\sigma_f^2)$ with
$\sigma_f^2=\mathrm{Var}_\pi \big(f(\theta_0)\big)+2\sum_{k\ge1}\mathrm{Cov}_\pi \big(f(\theta_0),f(\theta_k)\big)$.
Here $\theta_0 \sim \pi$, and all expectations/covariances are taken under $\pi$. See, e.g., meyn2009markov,jones2004markov,roberts2004general.}
\[
\|\mu_1 - \mu_2\|_{\mathrm{TV}} = \sup_{A \in \mathcal{F}}\bigl|\mu_1(A)-\mu_2(A)\bigr|,
\]
where \(\mu_{1}\) and \(\mu_{2}\) are probability measures on the same measurable space \((\Omega,\mathcal F)\). For a Markov kernel \(P\) with stationary distribution \(\pi\) and any initial law \(\mu_{0}\), we monitor
\(\|\mu_{0}P^{n}-\pi\|_{\mathrm{TV}}\) to quantify convergence.
The chain is called \textit{geometrically ergodic} if there exists \(0<\gamma<1\) such that
\[
\|\mu_0 P^n - \pi \|_{\mathrm{TV}}
\le
(1-\gamma)^n \|\mu_0 - \pi\|_{\mathrm{TV}},
\]
implying an exponential decay of rate \(1-\gamma\) in TV distance.
In a finite state space, assume \(P\) is irreducible, aperiodic, and \(\pi\)–reversible. Then \(1\) is a simple eigenvalue of \(P\) and all other eigenvalues lie strictly inside the unit circle. Let \(\lambda_2\) be the second–largest eigenvalue in absolute value. The \textit{spectral gap} is
\(\gamma(P)=1-|\lambda_2|\), and the spectral theorem yields
\[
\|\mu_0 P^n - \pi\|_{\mathrm{TV}} \le C (1-\gamma(P))^n
\]
for a constant \(C\) depending only on \(\mu_0\) and \(\pi\).
This finite–state characterization extends to general state spaces when $P$ is reversible. Let
\[
L_2(\pi)=\Bigl\{f:\Omega\to\mathbb{R}:\int f^2 d\pi<\infty\Bigr\},\qquad
\langle f,g\rangle=\int f g d\pi,\qquad
L_2^0(\pi)=\{f\in L_2(\pi):\langle f,1\rangle=0\}.
\]
Reversibility makes the Markov operator $(Pf)(x)=\int f(y) P(x,dy)$ self–adjoint on $L_2(\pi)$ and a contraction on $L_2^0(\pi)$. Its \textit{(absolute) spectral gap} is
\[
\gamma(P) = 1-\sup_{\substack{f\in L_2^0(\pi)\\ \|f\|=1}}\bigl|\langle Pf,f\rangle\bigr|
= 1-\sup_{f\in L_2^0(\pi)}\frac{\|Pf\|}{\|f\|}.
\]
In finite state spaces this coincides with $1-|\lambda_2|$. For the rest of the paper we use this spectral–gap viewpoint to compare convergence rates. Our theoretical results cover three settings: (i) a finite state space with $m^2$ states and $m$ modes; (ii) a general state space where $\pi$ is a mixture of $m$ log–concave components; and (iii) a general state space with local non-identification. The finite state space case is attractive by its simplicity and clear intuition, while it is less common in practice, so we place the results for that case in the Appendix (ref).
\subsubsection{Multi-modality}
Multimodal posteriors arise frequently in applied econometrics, even beyond the mixture models that have been extensively discussed in statistics fruhwirth2006finite. For example, regime-switching models deliver separate likelihood peaks that correspond to distinct combinations of structural parameters and latent states (diebold2001long). In SVAR models identification of structural parameters could hold locally but not globally depending on identifying restrictions (bacchiocchi2025locally). Taken together, these examples show that multimodality is a pervasive feature, making it a natural and important focus of our discussion.
In this section, we analyze a canonical example inspired by guan2007small to illustrate the challenges of sampling from multimodal distributions and the effectiveness of an identification–aware approach. Employing a circular topology, which is boundary-free, circumvents the endpoint treatments of a bounded interval.
\begin{example}
The state space $\Omega$ is a one-dimensional circle of circumference $4L$, represented as the interval $[-2L, 2L]$ with its endpoints identified. The target distribution $\pi(\theta)$ is bimodal:
\[
\pi(\theta) \propto
\begin{cases}
\exp(-\nu |\theta|), & \text{if } \theta \in [-L, L] \\
\exp(-\nu (2L - |\theta|)), & \text{if } \theta \in [-2L, -L) \cup (L, 2L]
\end{cases}
\]
where $L \gg 1$ and $\nu > 0$. The parameters $L$ and $\nu$ control the problem's difficulty: large $L$ increases the separation between modes, while large $\nu$ makes each mode sharper and more concentrated.
\end{example}
\begin{figure}[h!]
\caption{Target $\pi$ defined on a 1-D circle, with $L=10, \nu =2$}
\end{figure}
A standard Metropolis Hastings sampler with a local $\delta$-ball random walk mixes poorly on this two–mode circle. As shown in Proposition (ref), its spectral gap collapses exponentially:
\[
\gamma(P)\ \le\ C_1 e^{-\nu (L-\delta)}.
\]
This indicates that the sampler becomes prohibitively slow, effectively getting trapped in one mode, whenever the modes are sharp (large $\nu$) or far apart (large $L$).
The slow mixing is caused by the low-probability barrier between the modes. To break this barrier, we design a proposal that respects the inherent symmetry of the state space. Let \[
s(\theta) =
\begin{cases}
\theta+2L, & \theta\in[-2L,0),\\[2pt]
\theta-2L, & \theta\in[0,2L).
\end{cases}
\]
be the antipodal shift and define the teleport kernel
\[
T(\theta,\cdot)=\tfrac12 \delta_\theta(\cdot)+\tfrac12 \delta_{s(\theta)}(\cdot), \quad\text{where }\delta_\theta(\cdot)\text{ denotes the Dirac measure at } \theta
\]
Let $P$ be one MH step with a symmetric $\delta$-ball proposal. The identification–aware transition is the two–stage composition $\widetilde P =\ PT $, i.e., at each iteration first apply $T$ (stay put with probability $1/2$ or jump to $s(\theta)$ with probability $1/2$), then perform one local MH step from the resulting point. It is worth noting that $\widetilde P$ is by itself reversible in this specific setting.
This small change to the proposal mechanism has a profound impact on performance. As proven in Proposition (ref) in Appendix (ref), we can obtain a uniform lower bound on the spectral gap:
\[
\gamma(\widetilde P)\ \ge\ C_2\ \min\{(\nu\delta)^2,\ 1\},
\]
so choosing $\delta \propto \nu^{-1}$ keeps $\gamma(\widetilde P)$ bounded away from zero independently of $L$ and $\nu$.
We now introduce a more general setup.
Let $\pi$ be a probability density on a connected state space
$\Omega\subseteq\mathbb{R}^n$. Fix $h>0$ with $m \geq 2$ modes. We assume that the support of $\pi$ can be decomposed into $m$ disjoint,
open components $A_1,\dots,A_m$ with
$\pi (\cup_{i=1}^m A_i )=1$ such that each $\pi$ has only one mode in each $A_i$.
Let $w_i:=\pi(A_i)$, $\pi_i(\cdot):=\pi(\cdot\cap A_i)/w_i$, $\mu_i\in A_i$ be the mode of $\pi$ restricted to $A_i$, and write
$d_i:= \operatorname{dist}(\mu_i, \partial A_i)=\inf_{z\in \partial A_i}\|\mu_i - z\|_2$, where $\partial A_i$ stands for the boundary of the closure of $A_i$, and $d_*:=\min_i d_i$.
\begin{assumption}
\begin{enumerate}[label=(\roman*)]
• There exist constants $c_i\ge1$ and
$\nu_i>0$ such that, for all $r\in(0,d_i]$,
\[
\pi_i \big(\{\theta\in A_i:\ \|\theta-\mu_i\|\ge r\}\big)
\le c_i e^{-\nu_i r}.
\]
Set $c_{\max}:=\max_i c_i$, $\nu_{\min}:=\min_i \nu_i$, and
$\nu_{\max}:=\max_i \nu_i$.
• Fix a step size $\delta>0$ and define
$A_i^{\mathrm{int}}:=\{\theta\in A_i:\ \operatorname{dist}(\theta,\partial A_i)\ge\delta\}$.
There exist $\varepsilon_1\in(0,1]$, $\varepsilon_2>0$, and measurable cores
$A_i^\circ\subset A_i$ with $\pi_i(A_i^\circ)\ge\varepsilon_1$ such that, for
every $i\neq j$ and all $\theta\in A_i^\circ$,
$$
T(\theta,A_j^{\mathrm{int}})\ \ge\ \varepsilon_2,
$$
where $T(\theta,\cdot)$ is a $\pi$-reversible teleport kernel used in the composition.
• For each $i$ there exists a convex set
$C_i\subset A_i$ with $\pi_i(C_i)\ge\alpha_i>0$ and
$\delta_0>0$ such that $C_i^{\mathrm{int}}
:=\{\theta\in C_i:\ \operatorname{dist}(\theta,\partial C_i)\ge\delta_0\}\neq\emptyset$,
and $\pi$ is continuous and bounded above/below on $C_i^{\mathrm{int}}$.
For some $n_0\in\mathbb N$ and $\eta_0>0$,
\[
\inf_{\theta\in A_i} P^{ n_0}(\theta,\cdot)\ \ge\ \eta_0 m_i(\cdot),
\qquad\text{with } m_i \text{ a probability measure supported on }C_i^{\mathrm{int}}.
\]
\end{enumerate}
\end{assumption}
Assumption (ref)(ref) is a mode-wise concentration condition: within each region $A_i$ it bounds the tail mass away from the mode $\mu_i$ up to the radius $d_i$. It is used to upper-bound the conductance of the RWM by showing that the $\delta$-boundary layer near $\partial A_i$ has exponentially small $\pi_i$-mass as $d_i$ grows, hence $\gamma(P)$ is exponentially small in $d_i$. Assumption (ref)(ref) states that a single teleport step places a fixed fraction of probability into the interior $A_j^{\mathrm{int}}$ of any other mode, uniformly over a core $A_i^\circ$. This guarantees the inter-mode communication for the composed kernel $\widetilde P=PT $. Assumption (ref)(ref) provides a within-mode “small set” $C_i\subset A_i$ on which a Doeblin–type minorization holds for some finite number of MH steps. This yields a uniform positive lower bound on the spectral gap of $P$ restricted to $A_i$ (see meyn2009markov).
An alternative condition for Assumptions (ref)(ref) and (ref) is log-concavity: if $\pi_i$ is log-concave on $A_i$, Lemma (ref) yields the mode-wise exponential tail bound required by (ref)(ref). If there exists a convex $C_i\subset A_i$ on which $\pi$ is log-concave, then for sufficiently small $\delta$, Lemma (ref) gives a positive spectral gap for the $\delta$-ball MH on $C_i$.
A widely used class of statistical models that meets these conditions is the $k$-mixture Gaussian $\pi(\theta)=\sum_{i=1}^k w_i \mathcal N(\theta;\mu_i,\Sigma_i)$ with well-separated means, take $A_i=\{\theta: w_i\varphi_i(\theta)=\max_j w_j\varphi_j(\theta)\}$. Separation guarantees $d_i>0$. On each $A_i$ the mixture behaves sub-Gaussian around $\mu_i$, so Assumption (ref)(ref) holds with $\nu_i$ comparable to $\lambda_{\max}(\Sigma_i)^{-1/2}$. Choosing $T(\theta,\cdot)=\pi(\cdot | A_j)$ and $\delta$ smaller than the interior margin of $A_j$ gives Assumptions (ref)(ref)–(ref). For Assumption (ref)(ref), take $C_i=\{\ (\theta-\mu_i)^\top\Sigma_i^{-1}(\theta-\mu_i)\le r_i^2\ \}\subset A_i$ with $r_i$ small. Then, $\pi$ is bounded above/below on $C_i^{\mathrm{int}}$ and the $\delta$-ball MH (with $\delta\le \delta_0$) satisfies the required minorization.
\begin{proposition}
Let $\pi$ be a probability density on $\Omega\subseteq\mathbb{R}^n$, and suppose Assumption (ref) holds with sets $A_1,\dots,A_m$, modes $\mu_i$, radii $d_i$, and the teleport kernel $T$.
Let $P$ be the RWM kernel with the uniform $\delta$--ball proposal,
\[
q(\theta,\cdot)=\mathrm{Unif}\bigl(B(\theta,\delta)\bigr),\qquad \delta>0,
\]
and let $\bar P := \frac12PT +\frac12 TP$ denote the reversible composed kernel. Then:
\begin{enumerate}
• For any $\delta>0$,
\[
\gamma(P) \le 2 c_{\max} \exp \big\{-\nu_{\min} (d_*-\delta)\big\},
\qquad d_*:=\min_{1\le i\le m} d_i.
\]
• There exists a constant $c_0>0$, depending only on $(\varepsilon_1,\varepsilon_2,n_0,m,\eta_0)$ (and not on the separations $d_i$), such that
\[
\gamma(\bar P)\ \ge\ c_0.
\]
\end{enumerate}
\end{proposition}
Refer to Appendix (ref) for a detailed proof.
Proposition (ref) shows that the plain $\delta$-ball RWM–MH can mix very slowly: its spectral gap admits an upper bound of order $\exp\{-\nu_{\min}(d_*-\delta)\}$, so it deteriorates exponentially as the modes become farther apart (large $d_*$) and/or the target is more sharply concentrated (large $\nu_{\min}$). In contrast, the IA–MH composition $ \bar P$ removes any dependence on inter-mode separation: under Assumption (ref) its gap is bounded below by a positive constant that does not involve the $d_i$. If Assumption (ref)(ref) and (ref) are replaced by log-concavity, the lower bound can depend on the within-mode “condition number” $\nu_{\min}/\nu_{\max}$. Choosing $\delta\propto 1/\nu_{\max}$ yields a gap that is uniform in the separations and degrades only as $\nu_{\max}$ grows relative to $\nu_{\min}$, not with the distance between modes.
For the Gaussian RWM case, we pursue an alternative set of assumptions: instead of working with Assumption (ref), we adopt a log-concavity framework and strengthen it to strong log-concavity. The corresponding result is stated in Proposition (ref) in Appendix (ref).
Proposition (ref) establishes that Hamiltonian Monte Carlo (HMC) can also experience significant bottlenecks in multi-modal settings, leading to an exponentially small spectral gap.
\begin{proposition}[Hamiltonian Monte Carlo]
Under the same setup as Proposition (ref), let $L_h$ decompose into $m$ disjoint open components $A_1, \dots, A_m$ with $\pi(\cup_i A_i)=1$. For each $i$, define $w_i = \pi(A_i)$, $\pi_i(\cdot) = \pi(\cdot \cap A_i)/w_i$, and let $\mu_i \in A_i$ be a mode of $\pi_i$. Set
$
d_i := \operatorname{dist}(\mu_i, \partial A_i)$ and $ d_* := \min_i d_i$.
Let Assumption (ref)(ref) hold.
In addition, assume the potential $U=-\log\pi$ is $L_s$–smooth on each $A_i$: $$\|\nabla U(x)-\nabla U(y)\|\le L_s\|x-y\|,\qquad \text{for all } x,y\in A_i$$.
Let $P$ be the standard HMC kernel with Gaussian momentum $p_0\sim\mathcal N(0,\sigma^2 I_n)$, $\ell$ leapfrog steps of size $\eta$, and the usual Metropolis accept/reject step. Then, for fixed $(\ell,\eta,\sigma^2)$,
$$
\gamma(P) \le C
\exp \Big(- \min\Big\{\tfrac12 \nu_{\min} d_*,\ \frac{c d_*^2}{\sigma^2(\ell\eta)^2}\Big\}\Big),
$$
for some constants $C,c>0$ that depend only on $(n,\ell,\eta,\sigma^2,L_s,c_{\max})$ but not on $d_*$.
\end{proposition}
The proof is provided in Appendix (ref).
Proposition (ref) shows that, in the multi-modal setup considered, HMC admits an exponentially small upper bound on its spectral gap as the inter–mode separation grows.
By contrast, as we show in Proposition (ref), the spectral gap of IA-RWM is bounded below by a constant independent of inter-mode distances, yielding provably faster mixing than standard HMC.
As for lower bounds, in general multimodal settings (without additional structure on the modes and barriers) there is no known universal, closed-form lower bound for the HMC spectral gap. Explicit lower bounds are available only for specific target families. Consequently, one cannot appeal to a general HMC lower bound to conclude that identification–aware HMC must mix faster than standard HMC.
As for IA–HMC, while a fully general spectral-gap lower bound is not yet available, the teleport mechanism directly targets the bottlenecks that hinder plain HMC, and empirical results indicate the same qualitative advantage.
Existing literature highlights the poor performance of HMC in multi-modal settings. For example, mangoubi2018does demonstrate that for certain multi-modal targets, HMC can perform even worse than random-walk Metropolis, as it struggles to cross between modes due to the geometry of Hamiltonian trajectories. These results align with the upper-bound analysis in Proposition (ref), where transitions between modes are shown to be exponentially unlikely. Similar arguments, such as those adapted from the conductance bounds in vishwanath2024repelling, also support the observation that HMC can be highly inefficient in multi-modal scenarios.
We leave establishing a general lower-bound theory for IA–HMC for future work.
\subsubsection{Set-identification}
Apart from multi-modality, another common source of identification failure in statistical models is set-identification. When the parameter is set-identified, there is typically a continuum of observationally equivalent parameter values forming a manifold in the parameter space.
Specifically, within an arbitrarily small neighborhood of a parameter, there always exists another parameter configuration that produces an identical data-generating process. In what follows, we compare the mixing times of the IA–RWM algorithm with those of the standard RWM in the presence of local unidentifiability. We begin by illustrating the concept with a concrete example.
\begin{example}
Consider the parameter space
\[
\Theta_D = X \times Y_D \subset \mathbb{R}^d,
\]
where \(d = m + r\). Here,
\(X \subset \mathbb{R}^m\) is a fixed, bounded set (with \(m \ge 1\)). Let \(Y_D = [-D,D]^r\).
Consider the target distribution \(\pi_D\) on \(\Theta_D\) such that $\pi_D$ is flat over the identified set, i.e., $\pi_D$ factors as
\[
\pi_D(\theta_x,\theta_y) = p(\theta_x) u_D(\theta_y), \quad (\theta_x,\theta_y)\in \Theta_D,
\]
where \(p(\theta_x)\) is a continuous density on \(X\), bounded away from \(0\) and \(\infty\).
\(u_D(\cdot)\) is the uniform density on \(Y_D\).
\end{example}
\begin{figure}[h!]
\caption{Target density \(\pi_D(\theta_x,\theta_y) = p(\theta_x) u_D(\theta_y)\), where \(D = 10\). Here, \(p(\theta_x) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{\theta_x^2}{2}\right)\) is the standard Gaussian restricted to \([-4,4]\), and \(u_D(\theta_y)\) is the uniform density on \([-10,10]\).}
\end{figure}
In this model, $\theta_x$ is assumed to be identified, while \(u_D\) is flat in the \(\theta_y\) direction, reflecting that \(\theta_y\) is unidentifiable. As \(D\to\infty\), the size of \(Y_D\) grows without bound.
For simplicity, we assume \(p(\theta_x)\) is the standard Gaussian density and use the standard RWM algorithm as a benchmark. Under the additional assumption that the proposal density \(q\) is sufficiently diffuse (so that, for a small number of steps, the transition probability \(q^n\left((\theta_x,\theta_y),(\theta_x',\theta_y)\right)\) is uniformly bounded below for all \(\theta_y\)), one can show that the RWM’s spectral gap \(\gamma(P)\) tends to 0 as the unidentifiable region expands. In contrast, the IA–RWM algorithm maintains a spectral gap \(\gamma( \bar P)\) that is uniformly bounded away from 0, which ensures faster convergence even as \(D\) grows. This result is formalized in Proposition (ref) and proved in Appendix (ref).
This idea extends beyond the rectangular setting for \(Y_D\). In differential geometry, one can often represent a family of smooth, compact, connected \(r\)-dimensional manifolds (known as the leaves of a foliation, or fibres) as the level sets of a smooth function, provided that certain regularity and topological conditions hold. Intuitively, think of the parameter space as stacked by smoothly curved surfaces rather than straight rectangular blocks. Each point belongs to one surface and one only. The surfaces never cross, and as you move the points, they shift gradually without sharp bends or gaps. Under these mild geometric conditions we can tag every surface with a single coordinate, the density depends only on that tag, and the identification–aware sampler can still update the uninformative directions by drawing a fresh point uniformly on the current surface. This is exactly the same strategy that works in the rectangular case. These assumptions are formalized as below.
\begin{assumption}
\begin{enumerate}[label=(\roman*)]
• There exists a $C^\infty$ map $\phi:\Theta\to \Phi\subset\mathbb{R}^{d-r}$ whose Jacobian
$D_\theta\phi$ has constant rank $d-r$ for every $\theta\in\Theta$.
• For every $u \in \Phi$, the fiber (i.e. observationally equivalent set)
\[
\mathcal{F}(u) := \{\theta \in \Theta : \phi(\theta) = u\}
\]
is compact and connected.
• For each $u\in\Phi$ there exists a smooth hypersurface $\Sigma(u)\subset\mathcal{F}(u)$ such that $\Sigma(u)$ separates $\mathcal{F}(u)$ into two subsets $\mathcal{F}(u)^{-}$ and $\mathcal{F}(u)^{+}$ with
$0<c_1\le \operatorname{Vol}_{\mathcal F}(\mathcal{F}(u)^{-})/\operatorname{Vol}_{\mathcal F}(\mathcal{F}(u))\le c_2<1$,
and $\Sigma(u)$ admits a tubular neighborhood of width at least $\varepsilon_{\min}>0$ inside $\mathcal{F}(u)$, i.e., there is a diffeomorphism
$\Sigma(u)\times(-\varepsilon_{\min},\varepsilon_{\min})\to\mathcal N_{\varepsilon_{\min}}(\Sigma(u))\subset\mathcal{F}(u)$.
• Assume
\[
\pi(d\theta) = f(u) w_u(\theta) \mu_u(d\theta) du,
\]
where $u=\phi(\theta)$, $0<\underline f\le f(u)\le \overline f<\infty$, and for each $u$ the function $w_u:\mathcal F(u)\to(0,\infty)$ satisfies
\[
0<w_{\min}\ \le\ w_u(\theta)\ \le\ w_{\max}<\infty \quad \text{for $\mu_u$–a.e.\ }\theta\in\mathcal F(u),
\]
and $w_u$ is locally Lipschitz on $(\mathcal F(u), d_{\mathcal F})$. That is, for every compact $K\subset \mathcal F(u)$ there exists $L_{u,K}<\infty$ such that
\[
|w_u(\theta)-w_u(\theta')|\ \le\ L_{u,K} d_{\mathcal F}(\theta,\theta')\qquad\text{for all }\theta,\theta'\in K,
\] where $d_{\mathcal F}(\theta,\theta')$ is the shortest-path distance within $\mathcal F(u)$ between $\theta$ and $\theta'$.
The conditional $\mu_u$ is the normalized $r$–dimensional Hausdorff measure on the fiber $\mathcal F(u)$:
\[
\mu_u(A) = \frac{\operatorname{Vol}_{\mathcal F}(A)}{\operatorname{Vol}_{\mathcal F}(\mathcal F(u))},
\qquad A\subset \mathcal F(u).
\]
• There exist $\rho>0$ and constants $0<c\le C<\infty$ such that, for every $\theta\in\Theta$, points with
$\|u-\phi(\theta)\|+\|s\|\le\rho$ admit a unique representation
\[
\theta=\Psi_\theta(u,s),\qquad u\in\mathbb{R}^{d-r},\ s\in\mathbb{R}^r,
\]
satisfying
\[
\phi(\Psi_\theta(u,s))=u, \quad \text{and }\]
\[c \|(u,s)-(u',s')\|\le \|\Psi_\theta(u,s)-\Psi_\theta(u',s')\|\le C \|(u,s)-(u',s')\|,
\]
and with Jacobian determinant of $\Psi_\theta$ bounded between $c$ and $C$.
• The proposal density $g$ is symmetric and translation-invariant, and satisfies:
(i) There exist $\delta>0$ and $c_g>0$ such that
\[
g(z)\ \ge\ c_g\qquad\text{for all }\ \|z\|\le \delta,
\]
and
(ii) its tail probability $\overline G(t):=\int_{\|z\|>t} g(z) dz$ has finite first moment,
\[
\int_{0}^{\infty} \overline G(t) dt\ <\ \infty.
\]
\end{enumerate}
\end{assumption}
The parameter vector $\phi(\theta) \in \mathbb{R}^{d-r}$ has a smaller dimension than $\theta$ and it corresponds to a
vector of reduced-form parameters commonly available in structural econometric models. See, for example, giacomini2021robust. Assumption (ref)(ref) ensures, by the Regular Level Set Theorem lee2012introduction, that each fiber $\mathcal{F}(u)$ is a smooth $r$-dimensional submanifold of $\Theta$.\footnote{A subset $M \subset \Theta$ is a \textit{smooth $r$-dimensional submanifold} if for every $x \in M$ there exists a neighborhood $U \subset \mathbb{R}^d$ of $x$ and a smooth map $F: U \to \mathbb{R}^{d-r}$ such that $M \cap U = \{y \in U : F(y) = 0\}$ and $DF(y)$ has full rank $d-r$ for all $y \in U$.}
It ensures that local neighborhoods of fibers behave regularly, avoiding pathological changes in fiber geometry. The assumption about the rank of $D_{\theta}\phi$ requires that the dimension of the $\phi$ corresponds to the dimension of minimally sufficient reduced-form parameters.
Assumption
(ref)(ref) rules out fibers that wander off to infinity or split into disjoint pieces, both of which can cause improper posteriors or poor mixing. In familiar SVAR settings, the rotation set (e.g., $SO(n)$) is compact, and after an ordering convention it is connected.
In Assumption
(ref) (ref), every fiber admits a smooth “mid-fiber” slice $\Sigma(u)$ that splits it into two parts of comparable size and has a uniform tubular neighborhood. For MCMC this means a small random-walk step cannot jump from one side of the fiber to the other unless the current point is close to the slice, which is exactly what lets us control crossing probabilities. The uniform neighborhood width is a standard geometric implication of smoothness and is made precise by the tubular neighborhood theorem lee2012introduction.
When $\pi(\theta)\propto L(y\mid\theta) p(\theta)$ and the likelihood is flat along each fiber $\mathcal F(u)$, Assumption (ref)(ref) matches the posterior disintegration with $f(u)\propto L(u) p(u)$ and $w_u(\theta)\propto p_u(\theta)$ (the prior’s conditional on $\mathcal F(u)$).
Thus $w_u$ encodes prior information along the observationally equivalent set: $w_u\equiv1$ corresponds to a uniform prior, while informative priors yield non-uniform $w_u$.
The bounds on $f$ and $w_u$ imply a uniform acceptance floor for moves that change $u$,
\[
\alpha(\theta,\theta')\ \ge\ \frac{\underline f w_{\min}}{\overline f w_{\max}} =: \alpha_0.
\]
Assumption (ref)(ref) provides uniform local coordinates $(u,s)$ near any parameter value, with controlled distortion of distances and volumes. In practice, this means we can vary $u$ and the fiber coordinate $s$ in a stable way everywhere in the parameter space, which is what we need for the geometric bounds used later.
Assumption (ref)(ref) requires the proposal $g$ to place positive mass on a fixed small ball and to have light tails. Both $\delta$--ball and Gaussian random walks satisfy this. These properties let us guarantee short moves occur with non-negligible probability and keep the contribution of very long jumps under control.
\begin{proposition}
Let Assumptions (ref)(ref)–(ref) hold.
For $u\in\Phi$, let $D(u):=\sup_{\theta,\theta'\in\mathcal F(u)} d_{\mathcal F}(\theta,\theta')$ and $D_{\max}:=\sup_{f(u)>0}D(u)$.
Let $P$ be the standard random–walk Metropolis kernel targeting $\pi$ (from (ref)) with proposal $g$ (from (ref)).
Define the teleport kernel $T$ by
\[
T(\theta,A) := \frac{\displaystyle\int_{A\cap \mathcal F(\phi(\theta))} w_{\phi(\theta)}(\xi)\, \mu_{\phi(\theta)}(d\xi)}{\displaystyle\int_{\mathcal F(\phi(\theta))} w_{\phi(\theta)}(\xi)\, \mu_{\phi(\theta)}(d\xi)},
\]
where $\mu_u$ is the normalized $r$–dimensional Hausdorff measure on $\mathcal F(u)$ from (ref), and set
\[
\bar P:=\tfrac12(PT )+\tfrac12(TP).
\]
Then:
\begin{enumerate}
• There exists $C<\infty$ independent of $D_{\max}$ such that $\gamma(P)\le C/D_{\max}$. In particular, $\gamma(P)\to0$ as $D_{\max}\to\infty$.
• There exist $n\in\mathbb N$ and $\varepsilon_0>0$, independent of $D_{\max}$, such that for all $\theta\in\Theta$ and measurable $A\subset\Theta$,
\[
\bar P^{\,n}(\theta,A)\ \ge\ \varepsilon_0\, \frac{\pi(A\cap T_0)}{\pi(T_0)},
\]
where $T_0\subset\Theta$ is a fixed measurable set with $\pi(T_0)>0$ independent of $D_{\max}$.
Consequently, $\bar P$ is uniformly ergodic and its $L_2(\pi)$ spectral gap is bounded below uniformly in $D_{\max}$. In particular, one may take
\[
\gamma(\bar P)\ \ge\ 1-\big(1-\varepsilon_0\big)^{1/n}\ >\ 0.
\]
\end{enumerate}
\end{proposition}
Proof can be found in Appendix (ref). Proposition (ref) formalizes that standard random-walk Metropolis mixes increasingly slowly as the fibers $F(u)$ grow in diameter, while the identification–aware RWM maintains uniformly fast mixing regardless of fiber size, because it refreshes uninformative directions via weighted draws along each fiber. Intuitively, IA–RWM avoids the bottleneck along nearly flat dimensions that traps standard RWM.
\section{Sampling Simulation}
In this section, we perform two simple simulation exercises to showcase the efficiency gain we achieved from employing the identification–aware step.
\begin{simulation}[Mixture Gaussian]
In this simulation exercise, I draw a sample of size 1,000 from the mixture Gaussian distribution $$X \sim p\mathcal{N}(\mu_1,\sigma_1^2)+(1-p)\mathcal{N}(\mu_2,\sigma_2^2),$$
where the true $(\mu_1, \mu_2,\sigma_1,\sigma_2, p) = (0, 20, 1, 5, 0.3)$.
Then, I sample the parameters $(\mu_1, \mu_2,\sigma_1,\sigma_2, p)$ based on the likelihood with 1,000 chains of length 100,000, randomize the initial guess of the parameters in each simulation, and use the random walk Markov chain as the baseline proposal distribution.
\end{simulation}
Figure (ref) plots the sample distribution of one out of the 1,000 chains. identification–aware Gaussian random-walk Metropolis-Hastings (RWM) is able to “teleport” across two modes while standard RWM is stuck in one of the two modes. In fact, based on 1,000 simulation exercises, the RWM consistently ends up trapped in one of the modes almost every time. In contrast, using that $K(\mu_1, \mu_2,\sigma_1,\sigma_2, p)=\{(\mu_1, \mu_2,\sigma_1,\sigma_2, p), (\mu_2, \mu_1,\sigma_2,\sigma_1, 1-p)\}$, the chain was able to explore both modes. Figure (ref) shows that the reported means of $(\mu_1,\sigma_1)$ are clustered near $(10, 3)$\footnote{Different scaling in the marginal distributions around each mode causes the posterior mean to deviate slightly from $(10, 3)$.} more tightly than those in Figure (ref), indicating that the identification–aware MCMC method reliably estimates the posterior means even when scaling varies between modes.
\begin{figure}[h!]
\begin{subfigure}[b]{0.5\textwidth}
\caption{Samples from RWM}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\caption{Samples from IA–RWM $\widetilde P$}
\end{subfigure}
\caption{Sampling distributions from a mixture Gaussian likelihood using RWM and IA–RWM $\widetilde P$}
\end{figure}
\begin{figure}[h!]
\begin{subfigure}[b]{0.5\textwidth}
\caption{Sample moments from RWM}
\end{subfigure}
\begin{subfigure}[b]{0.5\textwidth}
\caption{Sample Moments from IA–RWM}
\end{subfigure}
\caption{Comparison of sampled moments across \( 1,000 \) chains with length \( N = 1{,}000 \) }
\end{figure}
While our identification–aware MCMC method is highly effective at exploring multi-modal distributions, its gain in efficiency becomes less obvious when local identification fails. Specifically, in scenarios where maxima are not isolated but instead form continuous manifolds or more complex structures, the method's ability to thoroughly navigate the distribution is unclear.
\begin{simulation}[Conditional Gaussian]
In this simulation exercise, we sample two Gaussian distributions, one with two parameters and the other with seven parameters. They are both of the form
$$X \sim \mathcal{N}(\sum_{i}^k\mu_i, 1),$$
where $k =2$ and $10$, respectively.
Unlike the mixture Gaussian case, the observationally equivalent sets in this example will be $K(\mu_1,\dots, \mu_k)= \{\overline{\mu}_1,\dots, \overline{\mu}_k: \sum_{i}\overline{\mu}_i = \sum_{i} \mu_i \}$, an affine subspace with infinite elements.
\end{simulation}
In both experiments, a uniform prior is used. For
$k=2$, the true parameter values are set to $(\mu_1, \mu_2)=(0,2)$, and a sample of size 1,000 is generated. The naive Maximum Likelihood Estimator (MLE) is highly sensitive to the choice of initial values. For instance, initializing at $(0,0)$ yields estimates of $(\hat{\mu}_1, \hat{\mu}_2)= (0, 2)$.\footnote{Different starting points lead to varying MLE estimates, which is expected because any $(\mu_1,\mu_2)$ pair with the same sum results in identical likelihoods. This dependence on initial values arises solely from the implementation of the interior-point optimization method. Similar effects are observed in the movement of particles within SMC.} To see how well each algorithm explores the support, we bound each $\mu_i$ between $[-10, 10]$. The result of one simulation run with 10,000 iterations\footnote{For SMC, it has 10,000 particles with 10 iterations.} in Figure (ref) gives more credit to sequential Monte Carlo. Metropolis Hastings on average (across simulations) performs well, but it tends to be less stable and less evenly distributed on the global maxima.
\begin{figure}[h!]
\caption{Trace plot of posterior samples for \(\mu_1\) and \(\mu_2\) of size 10,000, $k=2$}
\end{figure}
Since standard SMC outperforms the other two classic algorithms in this example, we expanded the study to a higher dimensional scenario with $k=10$. This scaling allows for a more comprehensive evaluation of the algorithms' performance in higher-dimensional parameter spaces, where the number of particles and identification–aware MH (IA–RWM) samples remains relatively small compared to the dimensionality, limiting their ability to fully explore the space. In this scaled-up exercise, we compare only the performance of IA–RWM and SMC, utilizing the settings from herbst2014sequential. The true parameter values are set to $\mu_1=10$ and $\mu_i=0$ for $i\neq 1$. Intuitively, when dimension $n$ increases, the number of points needed to explore the full support grows exponentially. We ran SMC with 20 tempering stages and $10^8$ particles, effectively approaching the machine's hardware limit.
We sample the same amount of points with 10,000 Metropolis-Hastings sample, and 10,000 from $K(\mu_1,\dots, \mu_k)$ for each $(\mu_1,\dots, \mu_k)$, using the batch augmentation variant in (ref). Both methods were parallelized in MATLAB. On an Intel Xeon Gold 6246R CPU with 128 GB RAM, the SMC procedure took approximately 1,100 seconds to complete, whereas our method required about 700 seconds.
In Figure (ref), the spatial distribution of samples projected to the first two dimensions is depicted. SMC particles exhibit a clear clustering around the point $(10,0)$, which aligns closely with the initial particle distribution. This concentration suggests that SMC is strongly influenced by the starting values, potentially limiting its exploration of the parameter space. In contrast, IA–RWM samples are more evenly dispersed in the $[-10,10]^2 $ subspace, demonstrating a more thorough exploration and reduced dependence on initial conditions.
Figure (ref) illustrates the marginal kernel density estimates for the parameters $\mu_1$ and $\mu_2$. The SMC method shows a sharp concentration of $\mu_1$ around 10 and
$\mu_2$ around 0, reflecting the clustering observed in the scatter plot. This concentration near the starting values indicates a potential limitation in capturing the full posterior distribution's variability when local identification fails. In contrast, IA–RWM exhibits a more diffuse marginal density, approaching a uniform distribution while still reflecting a slight offset from the origin. This diffuse distribution underscores IA–RWM's capability to explore the parameter space more effectively, capturing a broader range of plausible parameter values.
Overall, the comparative analysis across the spatial distributions, marginal densities, and summary statistics underscores the enhanced performance and robustness of the IA–RWM method over the traditional SMC approach, particularly in scenarios where thorough exploration of high-dimensional, flat regions is essential.
\begin{figure}[h!]
\caption{2D projection of samples: SMC particles vs. IA–RWM samples in \((\mu_1, \mu_2)\) space}
\end{figure}
\begin{figure}[h!]
\caption{Marginal kernel density estimates for parameters \(\mu_1\) and \(\mu_2\): SMC vs. IA–RWM}
\end{figure}
\section{Identification-aware MCMC with informative prior}
The target posterior matters for teleportation solely through the prior conditional on the observationally equivalent set. Note that in earlier examples with a uniform prior, all observationally equivalent modes have the same posterior height, making the inferior performance of traditional samplers less apparent in the Bayesian case where the prior over observationally equivalent parameters is informative through the prior specification.
In this section, we explicitly incorporate the informative conditional prior into the teleportation step to ensure the resulting Markov chain maintains the posterior as its stationary distribution. We begin with an MA(1) example to show how the teleport step integrates into Bayesian samplers and improves posterior sampling efficiency.
\subsection{Moving Average Processes}
It is well known that an $\text{MA}(1)$ process,
\[
y_t = \epsilon_t + \theta \epsilon_{t-1},
\qquad \epsilon_t \sim \mathcal{N} \left(0,\sigma^2\right),
\]
exhibits an observational equivalence between $(\theta,\sigma)$ and $(\theta^{-1}, |\theta| \sigma)$. Throughout we analyze a single dataset generated at $(\theta,\sigma)=(0.5,1)$. The observationally equivalent point is $(\theta,\sigma)=(2,0.5)$, which induces the same likelihood.
Following plagborg2019bayesian, we consider independent priors on $(\theta,\sigma)$. The prior places $\theta\sim\mathcal{N}(1,0.5^2)$ and $\log \sigma\sim\mathcal{N}(0,0.25^2)$, which favors the vicinity of $(0.5,1)$ over $(2,0.5)$. For comparison, we also examine a likelihood-only (uniform prior) specification, in which the prior is effectively flat and the posterior coincides with the likelihood.\footnote{In the uniform case, the conditional “fiber” move between $(\theta,\sigma)$ and $(\theta^{-1},|\theta|\sigma)$ is a simple $1/2$ draw.}
We perform posterior inference in the transformed parameterization \((\theta, s=\log\sigma)\) using three samplers, each run for $50,000$ iterations: (i) a random-walk Metropolis (RWM; target acceptance \(0.234\)), (ii) the adaptive No-U-Turn Sampler (NUTS; target acceptance \(0.80\); hoffman2014no), and (iii) an identification–aware random-walk Metropolis (IA–RWM) that augments local RWM updates (target \(0.234\)) with a teleport move between observationally equivalent points.
NUTS is a Hamiltonian Monte Carlo method that adaptively adjusts both the leapfrog step size and the trajectory length, so the user need only supply gradients of the log posterior. It is well regarded for efficiently exploring the interior of a single mode, even in moderately high-dimensional settings. However, as noted by plagborg2019bayesian, performance can deteriorate under highly diffuse priors when the posterior is multi-modal, as exploration across modes may slow despite adaptation.
To assess sensitivity to the prior and initialization, we consider three setups: a uniform (likelihood-only) prior, an informative prior with the chain initialized at the true value $(0.5,1)$, and the same informative prior with the chain initialized at the observationally equivalent point $(2,0.5)$. For each configuration and sampler, we overlay in light gray the true marginal posterior obtained by direct grid integration of the joint log posterior and plot the sampler’s marginal density estimate from the draws.
\begin{figure}[h!]
\caption{$\theta$ posteriors. Light-gray curves are grid-based marginal posteriors; colored curves are KDEs from sampler draws. Vertical lines: true $\theta = 0.5$ (red dashed) and observationally equivalent $\theta = 2$ (black dotted).}
\end{figure}
\begin{figure}[h!]
\caption{$\sigma$ posteriors. Light-gray curves are grid-based marginal posteriors; colored curves are KDEs from sampler draws. Vertical lines: $\sigma = 1$ (red dashed) and observationally equivalent $\sigma = 0.5$ (black dotted).}
\end{figure}
From Figure (ref) and (ref) we can see that, under the informative prior, the posterior mass is strongly concentrated near $(0.5,1)$, with a smaller secondary mode near $(2,0.5)$.
When initialized near the dominant mode (row 2), both RWM and NUTS produce marginals that agree well with the gray truth curves.
However, when initialized at the observationally equivalent point (row 3), RWM and NUTS tend to get stuck in that local region for the entire sampling process and can under-represent the dominant mode, yielding marginals that deviate from the gray overlays and masking bimodality.
In contrast, IA–RWM, by explicitly traversing the observationally equivalent points, robustly mixes between the two regions and recovers the correct marginal shapes regardless of initialization.
Under the uniform (likelihood-only) specification (row 1), the mass allocated to the two observationally equivalent regions is more balanced. Here, too, IA–RWM delivers consistent bimodal marginals, while RWM/NUTS are more sensitive to their starting point.
Overall, the results show that set identification and multimodality can hinder standard samplers, whereas exploiting the observational-equivalence structure enables more effective posterior exploration.
\FloatBarrier
\subsection{Application: News Shocks and Business Cycles}
In this section we reassess the algorithms using the SVMA setup of plagborg2019bayesian. An $n$ variable SVMA writes the data as \[y_t=\sum_{\ell=0}^{q}\Theta_\ell \varepsilon_{t-\ell}, \qquad \varepsilon_t\sim \mathcal{N}(0,\Sigma). \]
We collect impulse responses in $\Theta=\{\Theta_{ij,\ell} : 1 \leq i,j \leq n\}$ and the shocks are assumed i.i.d. Gaussian with covariance matrix $\Sigma=\mathrm{diag}(\sigma_{1},\dots,\sigma_{n})$. Identification for this model has been studied extensively in lippi1994var and plagborg2019bayesian, and are briefly revisited in Appendix (ref). The application uses three quarterly U.S. series: TFP growth, real GDP growth, and the ex-ante real federal-funds rate, and three latent shocks: an unanticipated productivity shock, a technology news shock, and a monetary-policy shock. The sample spans 1954Q3 – 2007Q4 ($T=213$). Series are detrended using the Stock–Watson kernel smoother, and the real rate is the effective funds rate minus contemporaneous GDP-deflator inflation. We adopt a finite MA lag $q=16$: the Akaike Information Criterion suggests $q\approx 13$, while autocorrelation diagnostics support a slightly longer window, so $q=16$ is chosen to be conservative.
We place a multivariate Gaussian prior on the impulse responses. For each $(i,j)$ and horizon $0\le \ell\le q$,
\[
\Theta_{ij,\ell}\sim \mathcal N(\mu_{ij,\ell},\ \tau_{ij,\ell}^2),
\qquad
\mathrm{Corr}(\Theta_{ij,\ell+k},\Theta_{ij,\ell})=\rho_{ij}^{ k}\ \ (0\le \ell\le \ell+k\le q),
\]
and the vectors $\{\Theta_{ij,0:q}\}$ are a priori independent across $(i,j)$. Impact responses are normalized by fixing the own-shock impacts to unity, $\Theta_{jj,0}=1$, with zero prior variance $(\tau_{jj,0}^2=0)$. Shock standard deviations
$\sigma_1,\ldots,\sigma_n$ are mutually independent and independent of the IRFs, with
\[
\log \sigma_j \sim \mathcal N(\mu_j^\sigma,\ (\tau_j^\sigma)^2).
\]
Priors are centered on the log-linearized sticky-price DSGE model of sims2012news:
$\mu_{ij,\ell}$ equals the DSGE IRF, except for the news shock where the uncertain anticipation horizon is handled by setting the TFP mean to one-half of the DSGE impact and spreading that mass over $\ell=0,\ldots,6$. Variances $\tau_{ij,\ell}^2$ are calibrated so the DSGE IRFs lie inside 90% prior bands under reasonable parameter perturbations. Bands for news-shock IRFs are deliberately wide (often including zero), whereas monetary-policy IRFs are tighter to reflect stronger beliefs about their qualitative shape. Smoothness parameters use $\rho_{ij}=0.5$ for TFP responses and $\rho_{ij}=0.9$ for GDP and the real rate, encoding smoother behavior for output and rates and allowing spikier productivity. Shock-scale priors are deliberately vague, with $\mu_j^\sigma=\log 0.5$ and $\tau_j^\sigma=2$.
In his paper, plagborg2019bayesian uses a two-step heuristic to provide the No-U-Turn Sampler (NUTS) with a high-density starting value. First, he constructs a rough posterior-mode guess: the sample auto-covariance function of the data is computed, an invertible SVMA representation matching this auto-covariance is obtained, and a greedy search over all root flips of the characteristic polynomial is performed so that the candidate maximizes the Gaussian prior density.\footnote{This is essentially a discrete search among a subset of observationally equivalent points, as no rotation is performed.} This candidate mode is then blended with the prior mean along a convex grid. The weight that yields the highest posterior probability defines the initial parameter vector supplied to the sampler. The procedure ensures numerical stability (by starting from an invertible representation) and avoids low-density regions that would require a long burn-in. On the other hand, it starts in a high-probability region and a local Markov chain may still get trapped there if the target distribution is multi-modal.
We generate posterior samples of size 10,000 (thinned from 100,000, which is ten times the number in plagborg2019bayesian) with NUTS, IA–RWM (block-wise), and IA-NUTS separately, first using his original prior, and then with a bounded uniform prior. The identification–aware variant used here is the reversibilized teleport–local composition \(\bar P=\tfrac12(PT+TP)\) for both IA–RWM and IA–NUTS. Implementation details are provided in Appendix (ref).
Figure (ref) reports posterior results for $\log \sigma_i$ (marginal plots of $\Theta$ are in the Appendix) under the tighter prior. All samplers are initialized at the same point obtained by a greedy mode search. Because visualizing the full “true” posterior, let alone its marginals, is infeasible in this setting, we diagnose behavior using summary statistics and local mode analyses. On the trace plots and marginal densities, NUTS and IA-NUTS are more stable than IA–RWM and deliver very similar marginal shapes. Relative to NUTS, however, IA-NUTS yields (in our sample) a lower posterior mean, similar average log-posterior, a larger posterior mode, and greater variance (Table (ref)). The marginal for $\log \sigma_1$ and $\log \sigma_2$ suggests that IA-NUTS locates a distinct mode that coincides with the mode visited by IA–RWM. While the numerical values should not be read too literally, the table indicates that IA–RWM may have traversed a broader region of support, and that the target is at least bimodal (and possibly multimodal).
As an additional piece of evidence, IA–RWM’s draws frequently achieve high log-posterior values with non-trivial step sizes and the chain spends sustained time in a region that NUTS visits rarely (if at all). This persistent occupancy, together with relatively large log-posterior values, points to an alternative high-density region with non-negligible mass that is more readily uncovered by the IA-based procedures.
To further investigate, we performed a sample-based mode search that draws candidate seeds from each sampler and then conducts brief, box-constrained local optimization from condensed seed sets. This exercise consistently returned higher log-posterior maximizers from IA-NUTS or IA-RWM than from NUTS, and the top IA-NUTS modes were well separated from those reached by NUTS.
However, there are two caveats that could temper the interpretation of mode-search. First, none of the polished modes reproduces the prominent visual marginal mode of NUTS with very small $\log\sigma_2$ (below -2), and the optimizers can be sensitive to tuning parameters. The apparent marginal peak need not correspond to a distinct high-posterior maximizer in the full parameter space. It may reflect projection of a broader ridge or a different basin whose summit lies elsewhere in $\Theta$. Second, larger peak log-posterior values do not imply larger posterior mass. Our mode search certifies the existence of alternative high-density basins but does not quantify their volume.
Taken together, the evidence supports the practical advantage of IA-NUTS over standard NUTS in this problem: IA-NUTS more reliably discovers alternative basins with higher attained log-posterior than those found by NUTS under the same computational budget. At the same time, without being able to quantify the posterior mass of the explored regions, the extent of IA-RWM and IA-NUTS’s mixing remains unknown.
\begin{table}[htbp]
\caption{Summary statistics of 10,000 posterior draws}
\begin{tabular}{lccc}
\toprule
& NUTS & IA–RWM & IA--NUTS\\
\midrule
Log-posterior at mean & 75.82 & -312.66 & -582.326 \\
Avg log-posterior & 2.08 & 6.18 & 1.14 \\
Sum per-param SD & 23.70 & 25.98 & 36.49 \\
Highest log-posterior & 36.14 & 59.57 & 67.88 \\
\bottomrule
\end{tabular}
\end{table}
\begin{figure}[h!]
\caption{Trace and kernel-density plots of 10,000 posterior $\log \sigma$ draws, obtained by retaining every 10th sample after thinning.}
\end{figure}
Next, we replace the informative Gaussian prior with an independent, uniform prior on $(\Theta,\boldsymbol{\phi})$, where $\Theta=\{\Theta_{ij,\ell}\}$ and $\boldsymbol{\phi}=(\phi_1,\ldots,\phi_n)$ with $\phi_j=\log\sigma_j$. Specifically, $\Theta_{ij,\ell}\sim\text{Unif}[a_{ij,\ell},b_{ij,\ell}]$ and $\phi_j\sim\text{Unif}[L_j,U_j]$ independently across all indices (i.e., uniform on a hyper-rectangle), where $a_{ij}=b_{ij}=100, L_j=-8, U_j=5$ for all $i,j$. The sampler enforces these bound constraints via specular reflection at the boundaries.
Given the 153-dimensional parameter space and the flat directions created by set identification, and exacerbated under a uniform prior, the goal is to find a sampler that balances numerical stability with genuine state-space coverage. From Figure (ref), and Figure (ref) - (ref) in the Appendix, the IA-NUTS achieves that balance. Its trace plots are stationary yet display sustained movement across the support without prolonged residence near the parameter bounds, and its kernel densities are sharply peaked with credible shoulders and tails. Such shapes are expected when information is weak: many IRFs are locally close to zero under the normalization, the likelihood changes very little along broad manifolds, and the marginal projections of these manifolds concentrate probability near the origin while retaining non-negligible mass in the wings. IA-NUTS’s occasional global refresh moves reposition the chain across observationally equivalent regions, which helps prevent persistent max-depth saturation and reduces the artificial broadening that repeated boundary reflections can induce.
Baseline NUTS performs noticeably worse in this environment. The traces experience long excursions and show pronounced swings, and the corresponding marginals look overly broad and uneven across panels, suggesting the sampler is not settling into a stable geometry. In fact, baseline NUTS turns out to cost 3 times more gradient evaluations than IA-NUTS, suggesting frequent hits of maximum depth. By contrast, IA–RWM produces very persistent traces with tiny step-to-step movement, and its marginals collapse into narrow spikes near zero—patterns consistent with an overly conservative proposal that fails to traverse the flat directions. IA-NUTS strikes the better balance: its traces remain stable without sticking, and its marginals are tight where the data are informative yet still display believable shoulders and tails, indicating more faithful exploration of the weakly identified posterior.
\begin{figure}[h!]
\caption{Trace and density plots of 10,000 posterior $\log \sigma$ draws, obtained by retaining every 10th sample after thinning under the uniform prior.}
\end{figure}
\FloatBarrier
\section{Conclusions and discussion}
We proposed identification-aware sampling schemes and showed that they outperform the conventional RWM and HMC in terms of the speed of convergence.
The key idea of exploiting identified sets and introducing global moves along them in sampling algorithms readily extends to a broad class of Markov chain based methods, thereby expanding both their applicability and practical efficiency.
Moreover, although our exposition focuses on \(K(\theta)\) as the set of parameters sharing the same data generating process (i.e., an observationally equivalent set) or likelihood, this assumption can be relaxed. In practice, one may define \(K(\theta)\) as any subset of parameters yielding nearly identical or approximated likelihood (e.g., using Whittle approximations). When exact equivalence is replaced by approximate equivalence, one can introduce a suitable weighting or correction step to ensure that the sampler still targets the correct posterior distribution. This generalization is especially valuable in models or data scenarios where strict observational equivalences are difficult to characterize, but approximate regions of high posterior density can be identified. We leave the thorough treatment of this issue for future work.
\FloatBarrier