The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
82,945 characters
Exact Rejection Sampling for Non-Gaussian State Space Models
\title{Exact Rejection Sampling for Non-Gaussian State Space Models}
\author{Joshua C.C. Chan \\
{\small Purdue University}
}
\date{August 2026}
\maketitle
\onehalfspacing
\begin{abstract}
\noindent Rejection sampling requires a proposal that dominates the target by a known constant, generally unavailable for non-Gaussian state space models. We construct such a proposal for the latent state path, yielding independent exact smoothing draws and an unbiased likelihood estimator whose relative variance is at most $1/p-1$ per draw at acceptance probability $p$. The method covers scalar states with affine Gaussian dynamics and log-concave observation densities, including multivariate observations. Transition twisting makes the log target-to-proposal ratio separable, and tangent-line twists make each term nonpositive, producing an attained, sharp dominating constant. With a companding node placement, the accumulated envelope error is $O(T/G^2)$ for a sample of length $T$ with $G$ nodes per date, so $G\propto\sqrt{T}$ keeps acceptance bounded away from zero; for stochastic volatility, the required conditions hold almost surely. A simpler mode-centered grid shows the same scaling empirically. At $T=2{,}000$, acceptance is $75\%$, versus roughly $10^{-16}$ for the Gaussian envelope.
\smallskip
\noindent \textbf{JEL classification:} C11, C15, C32
\smallskip
\noindent \textbf{Keywords:} rejection sampling, exact simulation, log-concavity, simulation smoothing, state space model, stochastic volatility, duration model
\end{abstract}
\thispagestyle{empty}
\newpage
\section{Introduction} \label{s:intro}
Latent state variables are central to many models in macroeconomics and finance. Stochastic volatility, unobserved components, count and duration models, among others, require simulation from the conditional distribution of a latent state path. When the model is linear and Gaussian, the path can be drawn exactly and cheaply using standard simulation smoothers.\footnote{These include the smoothers of \citet{CK94}, \citet{dS95}, \citet{DK02}, \citet{CJ09} and \citet{MMP11}.} Outside that class, state simulation typically relies on Metropolis--Hastings updates, importance sampling, particle methods, or approximations that replace the measurement density by a mixture of normals so that linear Gaussian machinery applies \citep{KSC98, OCSN07}. These methods can be highly accurate, but generally do not provide independent exact draws of the entire state path along with a computable global bound on the target-to-proposal ratio.
This paper develops an exact rejection sampler for a class of non-Gaussian state space models with a scalar latent state, affine Gaussian dynamics and a finite, continuously differentiable, concave observation log-density. The observation may be multivariate. Throughout, ``exact'' refers to the conditional state draw at a fixed parameter vector: an accepted path is distributed exactly from the joint smoothing distribution of the intended model. There is no discretization or approximation-model error, and the state draw requires neither burn-in nor a Markov-chain correction. The exactness claim concerns the algorithm in exact arithmetic; the implementation is subject to floating-point roundoff.
The main obstacle to rejection sampling the entire path is to construct a proposal that dominates the target up to a \emph{known} constant. For a $T$-dimensional state path, no practically useful constant has generally been available. A Gaussian envelope centered at the posterior mode provides a computable bound, but its acceptance probability decays exponentially with $T$; for the observation densities considered here, changing the Gaussian covariance does not in general resolve the problem. Existing non-Gaussian simulation smoothers instead construct accurate approximations to the smoothing distribution. Their target-to-proposal ratio may be bounded, but the value of that bound is not computed, whereas rejection sampling requires it.
Our construction makes the global bound tractable in two steps. First, twisting the transition densities makes the log target-to-proposal ratio, apart from an additive constant, a \emph{sum of univariate functions} of the states. Its $T$-dimensional supremum therefore separates into $T$ scalar suprema. Second, we take each twist to be the minimum of tangent lines to a concave backward function. The resulting tangent envelope dominates that function and touches it at the nodes, so every scalar supremum is exactly zero. At the same time, each twisted transition remains a finite mixture of truncated normals with closed-form normalizing constants and draws. The global dominating constant is then simply the normalizing constant of the first twisted density, computed as a by-product of the backward recursion.
The paper makes two main contributions. The first is an exact simulation smoother for the model class described above. We prove that the tangent-twisted proposal dominates the unnormalized smoothing density everywhere on $\mathbb{R}^T$, with a dominating constant that is attained and is the smallest admissible for the proposal. Rejection sampling therefore produces independent draws from the exact joint smoothing distribution. To our knowledge, this is the first practical rejection sampler for the joint smoothing distribution of a continuous, unbounded state process with Gaussian transitions for which the global dominating constant is both computable and sharp.
We also characterize how the computational cost of exactness scales with the path length. In nondegenerate cases, a fixed per-period envelope error accumulates along the path and causes acceptance to deteriorate exponentially in $T$. For the tangent construction, the per-period error is instead of order $G^{-2}$ in the number of nodes $G$, so the accumulated error, and hence $-\log$ of the acceptance probability, is of order $T/G^2$. Increasing $G$ at rate $\sqrt{T}$ therefore keeps the acceptance probability bounded away from zero. Theorem~\ref{thm:scaling} establishes this result under uniform curvature and proposal-tail conditions using a companding rule for node placement; for the stochastic volatility model, these conditions hold almost surely under the data generating process. The same square-root rate is recommended for the discretization filter of \citet{Farmer2021}, where it instead controls the error from replacing the continuous state process by a finite Markov chain.
The second contribution is likelihood evaluation with importance weights bounded by the same constant. The construction gives a nonnegative unbiased likelihood estimator whose variance is finite at every parameter value, with relative variance bounded in terms of the rejection acceptance probability rather than diagnosed from realized weights. Finite variance cannot be taken for granted for Gaussian proposals in this literature: \citet*{KoopmanShephardCreal2009} argue that it should be tested rather than assumed, and for the Gaussian approximation at the mode the relevant test fails throughout the stochastic volatility experiments and at all six empirical fits reported below. The bound also provides an exact reference for evaluating approximate smoothers, a computable uniform-ergodicity bound for the associated independence sampler, and a building block for pseudo-marginal and marginal likelihood estimation, as used for large vector autoregressions with stochastic volatility by \citet{chan23JE} and \citet*{CYZ26}.
The Monte Carlo results show that these guarantees remain useful for long state paths. For the stochastic volatility model, the rejection sampler accepts $75\%$ of proposals at $T=2{,}000$, compared with fewer than one in $10^{16}$ for the Gaussian envelope, with similarly high acceptance across seven other observation densities. Taking $G$ proportional to $\sqrt{T}$ keeps acceptance approximately constant across a sixteenfold increase in sample size, at about $0.87$ for the practical grid used in the implementation; the companding grid of Theorem~\ref{thm:scaling} shows the same scaling at a uniformly lower level. Under a common computing-time budget, the likelihood estimator is two to three thousand times more efficient than a bootstrap particle filter and comparable to the efficient importance sampler of \citet{RichardZhang2007}. Against importance sampling from the Gaussian approximation at the mode no such ratio is well defined, because those weights have infinite variance in every design we consider. What distinguishes the tangent estimator is therefore not the margin but that it attains this efficiency with weights bounded by a constant the algorithm computes. Exact draws also provide a direct benchmark for quantifying approximation error in standard state smoothers.
We illustrate the method with the stochastic conditional duration model of \citet{BauwensVeredas2004}. Marginal likelihoods strongly favor gamma over exponential errors in both series considered. Imposing exponential errors forces the latent state to absorb dispersion that would otherwise be captured by the error distribution, increasing the estimated state innovation variance several-fold.
The construction connects to several strands of the literature. The twisted transition densities are Doob $h$-transforms of the type used in sequential Monte Carlo \citep{WhiteleyLee2014, GuarnieroJohansenLee2017, HengBishopDeligiannidisDoucet2020}. The backward recursion is closely related to the efficient importance sampling method of \citet{RichardZhang2007} and the HESSIAN method of \citet{McCausland2012}, and more broadly to the Gaussian importance smoothers of \citet{DanielssonRichard1993}, \citet{SP97} and \citet{DurbinKoopman1997}. These methods construct accurate global proposals for importance sampling or Metropolis--Hastings, but do not supply the value of a global bound over the entire state path. The tangent construction is the path-space analogue of the piecewise linear upper hull used for univariate log-concave densities by \citet{Devroye1986} and, with adaptive refinement of the node set, by \citet{GilksWild1992}.
The rest of the paper is organized as follows. Section~\ref{s:method} establishes exactness of the tangent-twisted rejection sampler, and Section~\ref{s:implementation} presents the implementation and scaling theory. Section~\ref{s:uses} develops uses beyond exact state simulation, and Section~\ref{s:related} relates the construction to existing smoothers. Sections~\ref{s:MC} and~\ref{s:app} report the Monte Carlo evidence and the application, and Section~\ref{s:conclusion} concludes.
\section{Exact Rejection Sampling for the State Path} \label{s:method}
This section sets out the modeling framework and the assumption maintained throughout, and shows why the obvious Gaussian envelope cannot work. It then builds the proposal in two steps: twisting the transition densities makes the log target-to-proposal ratio separable across dates, and tangent-line twists make every term of that sum nonpositive. The main result is Theorem~\ref{thm:exact}: the algorithm is exact and its dominating constant is attained and smallest admissible.
\subsection{The Modeling Framework} \label{ss:model}
Let $\mathbf{h}=(h_1,\ldots,h_T)'$ denote a scalar latent state sequence and $\mathbf{y}=(\mathbf{y}_1',\ldots,\mathbf{y}_T')'$ the observations, where each $\mathbf{y}_t$ may be of any dimension. The observations are conditionally independent given the states, and the state evolves as a Markov chain with Gaussian transitions,
\begin{equation} \label{eq:model}
\log p(\mathbf{y}_t \,|\, h_t = u) = \ell_t(u), \qquad
h_1 \sim \mathcal{N}(\mu_1,\omega_1), \qquad
(h_t \,|\, h_{t-1}) \sim \mathcal{N}\left(a_t(h_{t-1}),\omega_t\right),
\end{equation}
for $t=2,\ldots,T$. Write $p_1(\cdot)$ for the density of $h_1$ and $p_t(\cdot \,|\, u)$ for the transition density. The functions $\ell_t$ and $a_t$ depend on static parameters, which are suppressed from the conditioning set to keep the notation uncluttered. The unnormalized posterior of the state path and its normalizing constant are
\begin{equation} \label{eq:target}
\gamma(\mathbf{h}) = p_1(h_1)\text{e}^{\ell_1(h_1)}\prod_{t=2}^{T}p_t(h_t \,|\, h_{t-1})\text{e}^{\ell_t(h_t)},
\qquad Z = \int_{\mathbb{R}^T}\gamma(\mathbf{h})\,\text{d}\mathbf{h},
\end{equation}
so that the target is $p(\mathbf{h} \,|\, \mathbf{y}) = \gamma(\mathbf{h})/Z$. The objective is to obtain independent draws from $p(\mathbf{h} \,|\, \mathbf{y})$ exactly.
The leading example is the stochastic volatility model $y_t = \text{e}^{h_t/2}\varepsilon_t$ with $\varepsilon_t\sim\mathcal{N}(0,1)$ and $h_t = \mu+\phi(h_{t-1}-\mu)+\sigma_h\eta_t$, $\eta_t\sim\mathcal{N}(0,1)$, where $|\phi|<1$ and the state is initialized from its stationary distribution. In the notation of \eqref{eq:model}, $a_t(u) = \mu+\phi(u-\mu)$, $\omega_t=\sigma^2_h$, $\mu_1=\mu$, $\omega_1 = \sigma^2_h/(1-\phi^2)$ and $\ell_t(u) = -\frac{1}{2}(u+y_t^2\text{e}^{-u})-\frac{1}{2}\log(2\pi)$. The restriction $|\phi|<1$ is what makes $\omega_1$ positive and finite; a nonstationary state is covered by \eqref{eq:model} under any other initialization with $\omega_1>0$.
We maintain two conditions throughout. The first restricts the state equation and is used only to establish a log-concavity property in Lemma~\ref{lem:concave}; the second restricts the observation densities.
\begin{assumption} \label{as:model}
The following hold.
\begin{enumerate}
\item[(i)] For each $t\geqslant 2$ the conditional mean $a_t(u) = \alpha_t+\beta_tu$ is affine in $u$ and the conditional variance $\omega_t>0$ does not depend on the state; and $\omega_1>0$.
\item[(ii)] For each $t$ the observation log-density $\ell_t:\mathbb{R}\rightarrow\mathbb{R}$ is finite, continuously differentiable and concave.
\end{enumerate}
\end{assumption}
Assumption~\ref{as:model}(ii) accommodates a wide range of observation densities. For the stochastic volatility model, $\ell_t''(u)=-\frac{1}{2}y_t^2\text{e}^{-u}\leqslant 0$. For the stochastic volatility in mean model $y_t = \alpha\text{e}^{h_t}+\text{e}^{h_t/2}\varepsilon_t$ of \citet{KoopmanHolUspensky2002}, $\ell_t''(u) = -\frac{1}{2}y_t^2\text{e}^{-u}-\frac{1}{2}\alpha^2\text{e}^{u}\leqslant0$. For Student's $t$ measurement error with $\nu$ degrees of freedom, $\ell_t''(u) = -\frac{1}{2}(\nu+1)z/(1+z)^2 \leqslant0$, where $z=y_t^2\text{e}^{-u}/\nu$, with equality only if $y_t=0$. For count data with a latent log-intensity, $y_t$ Poisson with mean $\text{e}^{h_t}$, $\ell_t''(u)=-\text{e}^{u}<0$. For a dynamic probit model with $\mathbb P(y_t=1)=\Phi(h_t)$, concavity of $\ell_t$ follows from the log-concavity of the normal distribution function. For the stochastic conditional duration model $y_t = \text{e}^{h_t}\varepsilon_t$ of \citet{BauwensVeredas2004} with exponential errors, $\ell_t(u) = -u-y_t\text{e}^{-u}$ and $\ell_t''(u) = -y_t\text{e}^{-u}<0$, and the fixed-shape gamma and Weibull versions have the same linear-minus-exponential form.
A further class deserves separate mention, because it shows that the restriction to a scalar state places no restriction on the dimension of the observation. In the common stochastic volatility model of \citet{CCM16}, $\mathbf{y}_t = \mathbf{B}'\mathbf{x}_t+\text{e}^{h_t/2}\bm{\varepsilon}_t$ with $\bm{\varepsilon}_t\sim\mathcal{N}(\mathbf{0},\bm{\Sigma})$, the observation log-density conditional on $\mathbf{B}$ and $\bm{\Sigma}$ is
\[
\ell_t(u) = -\frac{n}{2}u-\frac{1}{2}\text{e}^{-u}(\mathbf{y}_t-\mathbf{B}'\mathbf{x}_t)'\bm{\Sigma}^{-1}(\mathbf{y}_t-\mathbf{B}'\mathbf{x}_t)+\text{constant},
\]
which is concave in $u$ for every $n$. The same holds for a one-factor dynamic count model in which $y_{it}$ is Poisson with mean $\text{e}^{\alpha_i+\lambda_ih_t}$ for $i=1,\ldots,n$, since $\ell_t$ is then a sum of $n$ concave functions of $u$.
Assumption~\ref{as:model}(ii) requires concavity but not strict concavity. This matters in practice: in the stochastic volatility model an observation recorded as exactly zero gives $\ell_t(u)=-u/2-\frac{1}{2}\log(2\pi)$, which is affine rather than strictly concave.
\subsection{Why the Obvious Envelopes Fail} \label{ss:gaussian}
Before developing the proposal it is useful to see precisely why the obvious candidate fails. Let $\mathbf{Q}$ denote the prior precision matrix of $\mathbf{h}$ implied by \eqref{eq:model}, which is tridiagonal, and let $\mathbf{m}$ denote the posterior mode, with $t$th element $m_t$. Since each $\ell_t$ is concave, the tangent line at $m_t$ dominates it, and hence $\gamma(\mathbf{h}) \leqslant \text{e}^{M_G}q_G(\mathbf{h})$, where $M_G = \sum_{t=1}^{T}\ell_t(m_t)-\tfrac{1}{2}(\mathbf{m}-\bm{\alpha})'\mathbf{Q}(\mathbf{m}-\bm{\alpha})$, $\bm{\alpha}$ collects the prior means, and $q_G$ is the $\mathcal{N}(\mathbf{m},\mathbf{Q}^{-1})$ density including its normalizing constant.
The bound is attained at $\mathbf{h}=\mathbf{m}$, so the dominating constant is exact and the acceptance probability is $Z\text{e}^{-M_G}$, for which a Laplace approximation gives $\text{e}^{-\kappa}$ with $\kappa=\frac{1}{2}(\log|\mathbf{K}|-\log|\mathbf{Q}|)\geqslant0$ and $\mathbf{K} = \mathbf{Q}+\text{diag}\{-\ell_1''(m_1),\ldots,-\ell_T''(m_T)\}$. When the curvature $-\ell_t''(m_t)$ is bounded away from zero for a positive fraction of the sample and the eigenvalues of $\mathbf{Q}$ are bounded above uniformly in $T$, as they are under Assumption~\ref{as:model}(i) with fixed state-equation parameters, $\kappa$ grows linearly in $T$ and the acceptance probability decays exponentially.
Whether this deficiency can be repaired by changing the Gaussian covariance matrix depends on the tails of the observation density. Along the ray $\mathbf{h}=\mathbf{m}+r s_j\mathbf{e}_j$ for a coordinate $j$ and a direction $s_j\in\{-1,1\}$, if $\ell_j(m_j+rs_j) = o(r^2)$ as $r\rightarrow\infty$, then $\log\gamma \sim -\frac{1}{2}r^2\mathbf{Q}_{jj}$, and boundedness of $\log\gamma-\log q$ for a Gaussian proposal with precision $\mathbf{Q}_q$ requires $(\mathbf{Q}_q)_{jj} \leqslant \mathbf{Q}_{jj}$: no Gaussian envelope can then add diagonal curvature in coordinate $j$, which is exactly the curvature that the observation density contributes and that $\log|\mathbf{K}|-\log|\mathbf{Q}|$ measures. The tail condition holds for the standard stochastic volatility model, where $\ell_t(u)\rightarrow -u/2$ as $u\rightarrow\infty$, and for the Student's $t$, Poisson and dynamic probit specifications discussed earlier.
A second route to a bounded weight avoids concavity altogether: the defensive importance sampling of \citet{Hesterberg1995} mixes the proposal with the state prior $p_0$, and taking $q_{\text{def}} = (1-\lambda)p_0+\lambda q$ gives $\gamma(\mathbf{h})/q_{\text{def}}(\mathbf{h}) \leqslant (1-\lambda)^{-1}\exp\{\sum_t\sup_u\ell_t(u)\}$ from $T$ scalar maximizations. Two difficulties arise. The bound is finite only if every $\sup_u\ell_t(u)$ is finite, which Assumption~\ref{as:model} does not require: the affine $\ell_t$ produced by an observation recorded as exactly zero is unbounded above. When all $T$ suprema are finite, the bound is still too loose, because the implied acceptance probability is $(1-\lambda)\exp[-\sum_t\{\sup_u\ell_t(u)-\log p(\mathbf{y}_t \,|\, \mathbf{y}_1,\ldots,\mathbf{y}_{t-1})\}]$, in which every summand is nonnegative and vanishes only if $\ell_t$ is flat where the predictive state distribution puts mass. Whenever the state matters, in the sense that the average of these gaps is bounded away from zero, the acceptance probability again decays exponentially in $T$.
\subsection{Separability and Tangent Twisting} \label{ss:separability}
The proposal is built by twisting the transition densities of the state equation: each is multiplied by a positive function of the current state and renormalized, so that the result remains a Markov chain that can be simulated forward in one pass. Let $\psi_1,\ldots,\psi_T$ be functions on $\mathbb{R}$ for which the integrals below are finite, and define
\[
\widehat C_1 = \int_{\mathbb{R}}p_1(v)\text{e}^{\psi_1(v)}\,\text{d} v,
\qquad
\widehat C_t(u) = \int_{\mathbb{R}}p_t(v \,|\, u)\text{e}^{\psi_t(v)}\,\text{d} v, \quad t=2,\ldots,T.
\]
The twisted proposal is the Markov chain with initial and transition densities
\begin{equation} \label{eq:proposal}
q_1(h_1) = \frac{p_1(h_1)\text{e}^{\psi_1(h_1)}}{\widehat C_1},
\qquad
q_t(h_t \,|\, h_{t-1}) = \frac{p_t(h_t \,|\, h_{t-1})\text{e}^{\psi_t(h_t)}}{\widehat C_t(h_{t-1})},
\end{equation}
and joint density $q(\mathbf{h}) = q_1(h_1)\prod_{t=2}^{T}q_t(h_t \,|\, h_{t-1})$. Define the backward functions
\begin{equation} \label{eq:f}
f_T = \ell_T, \qquad f_t = \ell_t + \log\widehat C_{t+1}, \quad t=1,\ldots,T-1,
\end{equation}
and the residuals $d_t = f_t-\psi_t$.
\begin{lemma}[Separability] \label{lem:sep}
For every $\mathbf{h}\in\mathbb{R}^T$,
\begin{equation} \label{eq:ratio}
\frac{\gamma(\mathbf{h})}{q(\mathbf{h})} = \widehat C_1\exp\left\{\sum_{t=1}^{T}d_t(h_t)\right\},
\end{equation}
and consequently $\sup_{\mathbf{h}\in\mathbb{R}^T}\sum_{t=1}^{T}d_t(h_t) = \sum_{t=1}^{T}\sup_{u\in\mathbb{R}}d_t(u)$.
\end{lemma}
\begin{proof}
Dividing \eqref{eq:target} by the joint density implied by \eqref{eq:proposal}, every transition density cancels and
\begin{align*}
\frac{\gamma(\mathbf{h})}{q(\mathbf{h})}
&= \widehat C_1\prod_{t=1}^{T}\text{e}^{\ell_t(h_t)-\psi_t(h_t)}\prod_{t=2}^{T}\widehat C_t(h_{t-1})\\
&= \widehat C_1\exp\left\{\sum_{t=1}^{T}\left[\ell_t(h_t)-\psi_t(h_t)\right]+\sum_{t=1}^{T-1}\log \widehat C_{t+1}(h_t)\right\}.
\end{align*}
Collecting the terms that involve $h_t$ and using \eqref{eq:f} gives \eqref{eq:ratio}. The second claim follows because the summands depend on disjoint coordinates of $\mathbf{h}$.
\end{proof}
The identity \eqref{eq:ratio} is the standard weight decomposition for twisted proposals, implicit in \citet{RichardZhang2007}. What matters here is its consequence: a supremum over the $T$-dimensional path is a sum of $T$ scalar suprema, whatever the twisting functions. The optimal twist $\psi_t=f_t$ makes every residual vanish, so that $\widehat C_1=Z$ and every proposal is accepted, but it is not available in closed form. We instead choose $\psi_t$ so that each residual is nonpositive while each $\widehat C_t$ remains analytic. Specifically, for each $t$ let $V_t=\{v_{t,1}<v_{t,2}<\cdots<v_{t,G_t}\}$ be a finite nonempty set of nodes, so that $G_t\geqslant1$, and define
\begin{equation} \label{eq:tangent}
\psi_t(u) = \min_{1\leqslant j\leqslant G_t}\left\{f_t(v_{t,j})+f_t'(v_{t,j})(u-v_{t,j})\right\}.
\end{equation}
For a fixed concave function, \eqref{eq:tangent} is the upper hull of univariate adaptive rejection sampling \citep{GilksWild1992}; what is new here is that it is applied to functions the construction itself generates. The definition is recursive: $\psi_T$ is constructed from $f_T=\ell_T$, then $\widehat C_T$ determines $f_{T-1}$ through \eqref{eq:f}, which determines $\psi_{T-1}$, and so on down to $t=1$. Each step integrates the twisted density of the following period, so the construction is well posed only if concavity and differentiability are preserved under that integration; Lemma~\ref{lem:concave} establishes that they are, and records the domination property on which Theorem~\ref{thm:exact} rests. Its proof is in \ref{app:proofs}.
\begin{lemma}[Concavity and domination] \label{lem:concave}
Under Assumption~\ref{as:model}, for every $t=1,\ldots,T$ the following hold.
\begin{enumerate}
\item[(i)] $f_t$ is finite, concave and continuously differentiable on $\mathbb{R}$, so that \eqref{eq:tangent} is well defined.
\item[(ii)] $\psi_t$ is concave and piecewise linear, with $\psi_t(u)\geqslant f_t(u)$ for every $u\in\mathbb{R}$ and $\psi_t(v)=f_t(v)$ for every $v\in V_t$.
\item[(iii)] For $t\geqslant 2$, $0<\widehat C_t(u)<\infty$ for every $u$, and $\log\widehat C_t$ is concave and continuously differentiable; and $0<\widehat C_1<\infty$.
\end{enumerate}
\end{lemma}
Part~(ii) is also what makes the twisted kernel tractable. Where $\psi_t$ is affine, multiplying the Gaussian transition density by $\text{e}^{\psi_t}$ leaves a normal density with the same variance and a shifted mean, truncated to that piece. Hence $q_t(\cdot \,|\, h_{t-1})$ is a mixture of truncated normals, and Section~\ref{ss:pieces} obtains its weights and $\widehat C_t$ in closed form.
\subsection{Exactness and Bounded Weights} \label{ss:exact}
The two lemmas now combine: each residual $d_t=f_t-\psi_t$ is nonpositive by Lemma~\ref{lem:concave}(ii) and vanishes on the nodes, so the sum in \eqref{eq:ratio} is at most zero and the bound is attained.
\begin{theorem}[Exactness and the smallest dominating constant] \label{thm:exact}
Let Assumption~\ref{as:model} hold and let $\psi_t$ be given by \eqref{eq:tangent}. Then the following hold.
\begin{enumerate}
\item[(i)] $q$ is a probability density on $\mathbb{R}^T$.
\item[(ii)] $\gamma(\mathbf{h})\leqslant \widehat C_1 q(\mathbf{h})$ for every $\mathbf{h}\in\mathbb{R}^T$, with equality whenever $h_t\in V_t$ for every $t$; hence $\widehat C_1 = \sup_{\mathbf{h}}\gamma(\mathbf{h})/q(\mathbf{h})$, and no smaller dominating constant is admissible for this proposal.
\item[(iii)] Let $\mathbf{h}\sim q$ and $U\sim\mathcal{U}(0,1)$ be independent, and accept $\mathbf{h}$ if
\begin{equation} \label{eq:accept}
\log U \leqslant \sum_{t=1}^{T}\left[f_t(h_t)-\psi_t(h_t)\right].
\end{equation}
Then the conditional distribution of $\mathbf{h}$ given acceptance is exactly $p(\mathbf{h} \,|\, \mathbf{y})$, and a proposal is accepted with probability $Z/\widehat C_1$.
\end{enumerate}
\end{theorem}
The proof is in \ref{app:proofs}. Theorem~\ref{thm:exact} supplies what has been missing in this setting: a dominating constant that is known. It requires no numerical maximization, no safety factor and no search over the tails, because the supremum of each residual is zero by construction and is attained at the nodes. The theorem holds for any finite nonempty node sets. Exactness therefore does not depend on the accuracy of the posterior mode: poor placement lowers the acceptance probability but cannot invalidate the bound.
Theorem~\ref{thm:exact} also bounds the importance weights, which is what makes the twisted proposal useful beyond rejection sampling; the proof of the corollary is in \ref{app:proofs}.
\begin{corollary}[Bounded importance weights] \label{cor:is}
Let the conditions of Theorem~\ref{thm:exact} hold and let $w(\mathbf{h}) = \gamma(\mathbf{h})/q(\mathbf{h})$. Then $0<w(\mathbf{h})\leqslant\widehat C_1$ for every $\mathbf{h}$, $\mathbb E_q\{w(\mathbf{h})\}=Z$, and
\[
\frac{\text{Var}_q\{w(\mathbf{h})\}}{\left[\mathbb E_q\{w(\mathbf{h})\}\right]^2} \ \leqslant\ \frac{\widehat C_1}{Z}-1 \ =\ \frac{1}{p}-1,
\]
where $p$ is the acceptance probability of Theorem~\ref{thm:exact}. Hence for $\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(M)}$ drawn independently from $q$, the estimator $\widehat Z = M^{-1}\sum_{m}w(\mathbf{h}^{(m)})$ is unbiased for $Z$ with relative variance at most $(1/p-1)/M$.
\end{corollary}
\section{Implementation and Scaling} \label{s:implementation}
The construction of Section~\ref{s:method} is useful only if the twisted densities can be sampled and their normalizing constants evaluated. This section shows that both are available in closed form, because a piecewise linear twist leaves each twisted transition a mixture of truncated normals. It then sets out how the nodes are chosen, how many of them the sample size calls for, and states the algorithm, closing with the numerical precautions that an implementation requires.
\subsection{Closed-Form Expressions} \label{ss:pieces}
Write $g_{t,j}=f_t'(v_{t,j})$ for the slope of the $j$th tangent line and $\alpha_{t,j}=f_t(v_{t,j})-g_{t,j}v_{t,j}$ for its intercept, so that the line is $u\mapsto \alpha_{t,j}+g_{t,j}u$. The following lemma records the interval on which each line is active; the intersection points and the interval rule are equations (1) and (2) of \citet{GilksWild1992}, restated for a general strictly concave $f_t$ and with the interlacing $v_{t,j}<z_{t,j}<v_{t,j+1}$ established rather than assumed. The proof is in \ref{app:proofs}.
\begin{lemma}[Tangent intervals] \label{lem:intervals}
Suppose in addition that $f_t$ is strictly concave. Then $g_{t,1}>g_{t,2}>\cdots>g_{t,G_t}$, the intersection points
\[
z_{t,j} = \frac{\alpha_{t,j+1}-\alpha_{t,j}}{g_{t,j}-g_{t,j+1}}, \qquad j=1,\ldots,G_t-1,
\]
satisfy $v_{t,j}<z_{t,j}<v_{t,j+1}$ and hence $z_{t,1}<z_{t,2}<\cdots<z_{t,G_t-1}$, and
$\psi_t(v) = \alpha_{t,j}+g_{t,j}v$ for $v\in(z_{t,j-1},z_{t,j}]$, with the conventions $z_{t,0}=-\infty$ and $z_{t,G_t}=+\infty$.
\end{lemma}
When $f_t$ is concave but not strictly so, distinct nodes can share a tangent slope and the intersection $z_{t,j}$ is undefined. Equal slopes at two nodes imply that $f_t$ is affine between them, so the two tangent lines coincide rather than merely running parallel. The implementation therefore merges coincident lines and keeps one representative, then retains only the lines that appear on the lower envelope of the remaining family, discarding any that is nowhere the minimum, and forms intersections between adjacent distinct slopes. No slope is altered, so every retained line remains a tangent and $\psi_t\geqslant f_t$ is preserved exactly. When $f_t$ is affine a single line survives and the twisted transition reduces to a single untruncated normal, which is the case used as an implementation check in Section~\ref{s:MC}. After merging coincident lines and removing inactive lines, let $\mathcal{A}_t$ index the retained distinct lines. In the formulas below, sums over mixture components are understood to run over $j\in\mathcal{A}_t$. Under strict concavity, $\mathcal{A}_t=\{1,\ldots,G_t\}$.
On each of these intervals the twisted density is normal. Completing the square,
\[
p_t(v \,|\, u)\text{e}^{\psi_t(v)}
= \exp\left\{\alpha_{t,j}+g_{t,j}a_t(u)+\tfrac{1}{2}g_{t,j}^2\omega_t\right\}
\mathcal{N}\left(v;b_{t,j}(u),\omega_t\right),
\]
where $b_{t,j}(u)=a_t(u)+g_{t,j}\omega_t$ and $\mathcal{N}(v;b,\omega)$ denotes the normal density with mean $b$ and variance $\omega$ evaluated at $v$. Integrating over the $j$th interval gives the weight
\begin{equation} \label{eq:weights}
w_{t,j}(u) = \exp\left\{\alpha_{t,j}+g_{t,j}a_t(u)+\tfrac{1}{2}g_{t,j}^2\omega_t\right\}
\left[\Phi\left(\zeta^{+}_{t,j}\right)-\Phi\left(\zeta^{-}_{t,j}\right)\right],
\end{equation}
where $\zeta^{+}_{t,j} = (z_{t,j}-b_{t,j}(u))/\sqrt{\omega_t}$,
$\zeta^{-}_{t,j} = (z_{t,j-1}-b_{t,j}(u))/\sqrt{\omega_t},$ and $\Phi$ is the standard normal distribution function. The normalizing constant is $\widehat C_t(u)=\sum_{j\in\mathcal{A}_t}w_{t,j}(u)$, so that $q_t(\cdot \,|\, u)$ is a mixture of $|\mathcal{A}_t|\leqslant G_t$ normal densities that share the variance $\omega_t$ and are each truncated to one interval. A draw is obtained by selecting the $j$th component with probability $w_{t,j}(u)/\widehat C_t(u)$ and then drawing from $\mathcal{N}(b_{t,j}(u),\omega_t)$ truncated to $(z_{t,j-1},z_{t,j}]$ by inversion. Both steps are exact and together require $O(G_t)$ evaluations of $\Phi$.
The same expressions deliver the initial density: $q_1$ is obtained from \eqref{eq:weights} on replacing $a_t(u)$ by $\mu_1$ and $\omega_t$ by $\omega_1$, and $\widehat C_1$ is the corresponding sum of weights. The dominating constant of Theorem~\ref{thm:exact} is therefore computed as a by-product of the backward recursion, at the cost of a single additional evaluation.
Constructing $\psi_t$ requires $f_t$ and its first derivative at the nodes. The first is available from \eqref{eq:weights} through \eqref{eq:f}. The second is available at no additional cost. Differentiation under the integral sign, justified in the proof of Lemma~\ref{lem:concave}, gives
\begin{equation} \label{eq:Cprime}
\widehat C_t'(u) = \frac{\beta_t}{\omega_t}\int_{\mathbb{R}}\left(v-a_t(u)\right)p_t(v \,|\, u)\text{e}^{\psi_t(v)}\,\text{d} v,
\end{equation}
and dividing by $\widehat C_t(u)$ gives
\begin{equation} \label{eq:derivative}
\frac{\text{d}}{\text{d} u}\log\widehat C_t(u) = \frac{\beta_t}{\omega_t}\left\{\mathcal{M}_t(u)-a_t(u)\right\},
\qquad \mathcal{M}_t(u) = \mathbb E_{q_t(\cdot \,|\, u)}\left[h_t\right],
\end{equation}
so that $f_{t-1}'(u) = \ell_{t-1}'(u)+\beta_t\{\mathcal{M}_t(u)-a_t(u)\}/\omega_t$. The conditional mean of the twisted transition is a weighted average of truncated normal means,
\[
\mathcal{M}_t(u) = \frac{\sum_{j\in\mathcal{A}_t}w_{t,j}(u)\mu_{t,j}(u)}{\sum_{j\in\mathcal{A}_t}w_{t,j}(u)},
\qquad
\mu_{t,j}(u) = b_{t,j}(u)+\sqrt{\omega_t}\,
\frac{\varphi\left(\zeta^{-}_{t,j}\right)-\varphi\left(\zeta^{+}_{t,j}\right)}
{\Phi\left(\zeta^{+}_{t,j}\right)-\Phi\left(\zeta^{-}_{t,j}\right)},
\]
where $\varphi$ is the standard normal density (not to be confused with the persistence parameter $\phi$ of the stochastic volatility model). Every quantity appearing here is already formed in computing $\widehat C_t$, so the derivative recursion adds only $O(G_t)$ elementary operations per node.
\subsection{Choosing the Nodes} \label{ss:nodes}
The nodes affect only efficiency, so they should be placed where the smoothed states have mass. Our default implementation sets $v_{t,j}=m_t+w_js_t$, where $m_t$ is the $t$th element of the posterior mode of $\mathbf{h}$, $s_t$ is the corresponding standard deviation under the Gaussian approximation at the mode, and $w_1<\cdots<w_G$ is a fixed grid on $[-4,6]$. The mode is computed by Newton's method, with each iteration requiring an $O(T)$ band solve. The standard deviations are the square roots of diagonal elements of the inverse of a tridiagonal precision matrix and are obtained in $O(T)$ operations using the selected inversion recursion of \citet*{TakahashiFaganChin1973}, rather than by forming the inverse at $O(T^2)$. The grid is asymmetric because in the stochastic volatility and duration models the observation log-density decays much faster to the left of the mode than to the right. The favorable orientation is model specific---it is reversed for Poisson counts and depends on $y_t$ for dynamic probit.
We refer to this mode-centered, fixed-range construction as the \emph{practical grid}. The scaling theorem below is proved for a different \emph{theoretical grid}, obtained from a companding rule chosen to permit uniform control of both curvature and proposal tails. This distinction concerns efficiency, not exactness: Theorem~\ref{thm:exact} applies to any finite nonempty node sets, including both grids. We use the simpler practical grid as the default implementation.
\subsection{How Many Nodes: Bounds and Scaling} \label{ss:scaling}
How many nodes are needed has two answers. We first give a finite-sample bound that applies to arbitrary node sets, including the practical grid of Section~\ref{ss:nodes}. We then introduce a particular companding grid for which that bound can be strengthened to a uniform $T/G^2$ scaling result. Write $S_T=-\sum_{t=1}^{T}d_t(h_t)$ for the accumulated envelope residual of a proposed path. Since $d_t\leqslant0$, $S_T\geqslant0$, and the acceptance probability of Theorem~\ref{thm:exact} is $\mathbb E_q(\text{e}^{-S_T})$. Jensen's inequality gives
$ \mathbb E_q(\text{e}^{-S_T})\geqslant \exp\{-\mathbb E_q(S_T)\},$ and the expected residual admits the following finite-sample bound.
\begin{proposition}[Acceptance bound] \label{prop:rate}
Let the conditions of Theorem~\ref{thm:exact} hold, suppose $G_t\geqslant2$ for every $t$, let $I_t$ denote the interval spanned by $V_t$, let $\Delta_t$ be the largest spacing between adjacent nodes of $V_t$, and suppose $f_t'$ is Lipschitz on $I_t$ with constant $L_t$, which holds in particular if $f_t$ is twice differentiable there with $-f_t''\leqslant L_t$. Then
\[
0\leqslant \psi_t(u)-f_t(u)\leqslant \frac{L_t\Delta_t^2}{8} \quad \text{for } u\in I_t,
\]
and consequently
\[
\mathbb P(\text{accept}) \geqslant \exp\left\{-\sum_{t=1}^{T}\left[\frac{L_t\Delta_t^2}{8}+\mathbb E_q(R_t)\right]\right\},
\]
where $R_t = \{\psi_t(h_t)-f_t(h_t)\}\mathbf{1}\{h_t\notin I_t\}$ is the contribution from outside the node range.
\end{proposition}
The proof is in~\ref{app:proofs}. The proposition separates the two sources of envelope error. Inside the node range, the error is quadratic in the largest node spacing. If curvature is uniformly bounded and $\Delta_t=O(G^{-1})$, this contribution accumulates at rate $T/G^2$. The remaining term is the contribution from outside the node range. Thus, if the node ranges also provide sufficiently uniform tail coverage, taking $G$ proportional to $\sqrt{T}$ keeps acceptance bounded away from zero. At that rate the backward pass costs $O(TG^2)=O(T^2)$ and each proposal costs $O(TG)=O(T^{3/2})$, the latter being the marginal cost of an additional draw once the backward pass has been constructed.
Proposition~\ref{prop:rate} remains valid for the practical grid, but it does not impose the uniform control of the tail terms $\mathbb E_q(R_t)$ that an asymptotic claim requires. The following assumption and theorem instead construct a companding grid for which both the curvature and tail contributions can be controlled uniformly. The result is stated in centered coordinates. For deterministic or data-dependent centers $c_{t,T}$, write $\widetilde f_t(x)=f_t(c_{t,T}+x)$ and $X_t = h_t-c_{t,T}$ for the centered state. Both parts of the assumption are imposed over the single family of companding node sets constructed in Theorem~\ref{thm:scaling} from the constants $b$ and $W$ appearing below.
\begin{assumption}[Curvature and tail control] \label{as:scaling}
Each $\ell_t$ is twice continuously differentiable, and there exist constants $b>0$ and $\varsigma>2b/3$, a continuous function $W\geqslant1$ with $|\log W(x)-\log W(y)|\leqslant L_W|x-y|$ for all $x,y$, centers $c_{t,T}$, and finite nonnegative coefficients $A_{t,T}$ and $B_{t,T}$ with $\sup_TT^{-1}\sum_{t=1}^TA_{t,T}B_{t,T}<\infty$, such that for every $t$ and $T$, uniformly over the node sets of Theorem~\ref{thm:scaling}: (i) $0\leqslant-\widetilde f_t''(x)\leqslant A_{t,T}\,W(x)$ for every $x$; and (ii) $\mathbb E_q\big[\text{e}^{\varsigma X_t^2}\big]\leqslant B_{t,T}$. The centers and the coefficients may be data dependent, in which case they are measurable functions of the data, defined together with $b$, $\varsigma$ and $W$ on a common probability-one event, every expectation under $q$ is read conditionally on the realized data, and the average condition on $A_{t,T}B_{t,T}$ is required to hold almost surely on that event.
\end{assumption}
Part (i) can be checked directly from the observation density. The backward recursion cannot generate uncontrolled curvature on its own: Lemma~\ref{lem:curvature} in \ref{app:proofs} shows that $0\leqslant-f_t''\leqslant-\ell_t''+\beta_{t+1}^2/\omega_{t+1}$, whatever the number and placement of the nodes at later dates. Part (ii) is the substantive condition, because it restricts an object generated by the algorithm rather than a primitive feature of the model. It requires the centered proposal marginals to be sub-Gaussian, uniformly in $t$, $T$ and $G$, at an exponent $\varsigma$ strictly above $2b/3$; the proof of Theorem~\ref{thm:scaling} shows that this is exactly what the approximation bound consumes. It is not automatic: the marginals of $q$ are outputs of the backward recursion, and the bound can fail for otherwise admissible node sets. The quantification is joint---the assumption asserts the existence of $b$, $\varsigma$ and $W$ such that the associated grids satisfy both (i) and (ii). In \ref{app:svscaling}, we construct such a triple for the stochastic volatility model and verify both conditions for the resulting grids.
\begin{theorem}[Scaling of the acceptance probability] \label{thm:scaling}
Let Assumptions~\ref{as:model} and \ref{as:scaling} hold, and at each date place $G\geqslant2$ nodes at $v_{t,j} = c_{t,T}+x_{j,G}$, where $x_{j,G} = \Lambda^{-1}\{(j-\frac{1}{2})/G\}$ and $\Lambda$ is the distribution function of the density proportional to $\{W(x)\text{e}^{-bx^2}\}^{1/3}$. Then there is a constant $C$, depending only on $b$, $\varsigma$, $L_W$ and $W$, such that
\[
-\log\mathbb P(\text{accept}) \leqslant \frac{C}{G^{2}}\sum_{t=1}^{T}A_{t,T}B_{t,T}.
\]
\end{theorem}
The proof is in~\ref{app:proofs}. Under the average bound of Assumption~\ref{as:scaling}, the right-hand side is $O(T/G^2)$. Hence $G$ proportional to $\sqrt{T}$ keeps the acceptance probability bounded away from zero, $G/\sqrt{T}\rightarrow\infty$ drives it to one, and $G$ proportional to $\sqrt{T/\log T}$ prevents it from decaying faster than polynomially. These conclusions apply to the companding grid in the theorem. The Monte Carlo results in Section~\ref{s:MC} show that the practical grid of Section~\ref{ss:nodes} displays the same $T/G^2$ scaling empirically, while delivering higher acceptance in the designs considered there. For the stochastic volatility model, the assumptions underlying the theoretical grid can be verified rather than imposed.
\begin{corollary} \label{cor:svscaling}
Consider the stochastic volatility model of Section~\ref{ss:model} with $\sigma_h^2>0$ and $|\phi|<1$, the state initialized from its stationary distribution. With centers $c_{t,T} = \log(1+y_t^2)$ and envelope $W(x) = 1+\text{e}^{|x|}$, there are deterministic constants $b>0$ and $\varsigma>2b/3$, depending only on the model parameters, such that Assumption~\ref{as:scaling} holds almost surely under the data generating process. Consequently, for almost every data path there is a finite random constant $C_y$ for which
\[
-\log\mathbb P(\text{accept})\leqslant C_y\,\frac{T}{G^{2}}
\]
at the node sets of Theorem~\ref{thm:scaling}, simultaneously for every $T$ and every $G\geqslant2$.
\end{corollary}
The square-root rate has a counterpart in the discretization filter of \citet{Farmer2021}, which replaces the continuous state process by a finite Markov chain and, for a scalar state, recommends growing the number of chain states in proportion to $\sqrt{T}$. The two prescriptions control different objects. There, refinement drives the error of approximating the continuous-state model to zero fast enough for the approximate likelihood to support valid inference; here, exactness holds at every finite node set by Theorem~\ref{thm:exact}, and refinement governs only the acceptance probability.
\subsection{The Algorithm} \label{ss:algorithm}
Everything the algorithm requires is now in place: the closed-form weights and derivatives of Section~\ref{ss:pieces} and the nodes of Section~\ref{ss:nodes}. A single backward pass forms the twists $\psi_1,\ldots,\psi_T$ and delivers the dominating constant $\widehat C_1$ as a by-product, after which each proposal is one forward pass, accepted or rejected by \eqref{eq:accept}. It is summarized in Algorithm~\ref{alg:main}.
\begin{algorithm}[ht]
\caption{Tangent-twisted rejection sampling of the state path.}
\label{alg:main}
\begin{enumerate}
\item \textit{Nodes.} Compute the posterior mode $m_t$ and the scale $s_t$, and set $v_{t,j}=m_t+w_js_t$.
\item \textit{Backward pass.} Set $f_T(v_{T,j})=\ell_T(v_{T,j})$ and $f_T'(v_{T,j})=\ell_T'(v_{T,j})$. Then for $t=T,T-1,\ldots,2$:
\begin{enumerate}
\item[(a)] form $\psi_t$ from $\{v_{t,j},f_t(v_{t,j}),f_t'(v_{t,j})\}$, that is the slopes $g_{t,j}$, intercepts $\alpha_{t,j}$ and intersections $z_{t,j}$, using Lemma~\ref{lem:intervals} under strict concavity and the merging-and-pruning rule following that lemma otherwise;
\item[(b)] evaluate $\widehat C_t(v_{t-1,j})$ and $\mathcal{M}_t(v_{t-1,j})$ from \eqref{eq:weights} for $j=1,\ldots,G_{t-1}$;
\item[(c)] set $f_{t-1}(v_{t-1,j})=\ell_{t-1}(v_{t-1,j})+\log\widehat C_t(v_{t-1,j})$ and obtain $f_{t-1}'(v_{t-1,j})$ from \eqref{eq:derivative}.
\end{enumerate}
Finally form $\psi_1$ and compute $\widehat C_1$.
\item \textit{Proposal.} Draw $h_1\sim q_1$ and set $L=0$. Then for $t=2,\ldots,T$: compute the weights $w_{t,j}(h_{t-1})$ and $\widehat C_t(h_{t-1})$ from \eqref{eq:weights}; update
$L\leftarrow L+\ell_{t-1}(h_{t-1})+\log \widehat C_t(h_{t-1})-\psi_{t-1}(h_{t-1})$; and draw $h_t$ from the piecewise normal density $q_t(\cdot \,|\, h_{t-1})$. Finally update $L\leftarrow L+\ell_T(h_T)-\psi_T(h_T)$.
\item \textit{Accept-reject.} Draw $U\sim\mathcal{U}(0,1)$. Return $\mathbf{h}$ if $\log U\leqslant L$; otherwise return to step 3.
\end{enumerate}
\end{algorithm}
Step 1 costs $O(T)$ operations per Newton iteration, because the precision matrix is tridiagonal and each iteration is therefore a band solve, as in \citet{chan17}. Step 2 evaluates $G$ weights at each of $G$ nodes for each of $T$ periods and therefore requires $O(TG^2)$ evaluations of $\Phi$, where $G$ denotes the common number of nodes. Each proposal in step 3 requires $O(TG)$ evaluations, and the expected number of proposals per accepted draw is $\widehat C_1/Z$. Steps 1 and 2 depend on the data and on the model parameters, but not on the number of draws required, so their cost is amortized when many draws are taken at a fixed parameter vector.
Theorem~\ref{thm:exact} guarantees that the right-hand side of \eqref{eq:accept} is nonpositive, but floating-point error can produce small positive values. Three precautions are therefore important. The mixture weights \eqref{eq:weights} are accumulated by log-sum-exp rather than by direct exponentiation and summation. Differences of normal distribution functions in \eqref{eq:weights} are evaluated using log distribution and survival functions, according to the signs of their arguments, since direct subtraction loses precision when both arguments lie far in the same tail. Finally, we monitor the largest value of $\sum_t d_t(h_t)$ over all proposals; a value materially above machine precision indicates a numerical or coding error and is investigated rather than silently truncated.
These diagnostics are used throughout the Monte Carlo study and empirical application. On deterministic grids containing the nodes, tangent intersections, interval midpoints and far-tail points, the residual never exceeds $4.5\times10^{-13}$, even in the demanding designs with $T=4{,}000$, $\phi=0.995$ or $\sigma_h^2=0.5$; positive values occur only at the level of double-precision roundoff. Direct exponentiation and subtraction of normal probabilities are numerically unstable in these designs, so the log-domain calculations described above are used throughout.
\section{Uses Beyond Exact State Simulation} \label{s:uses}
Beyond producing independent exact state draws, the tangent-twisted proposal has two uses that are particularly important for this paper. First, its bounded importance weights give a nonnegative unbiased likelihood estimator with guaranteed finite variance. Second, exact draws provide a benchmark for assessing approximate state smoothers. We develop these two uses first and then record several additional consequences of the same bound. The backward pass of Algorithm~\ref{alg:main} depends on the data and on the model parameters, so it must be rebuilt whenever the parameter vector changes, at a cost of $O(TG^2)$. Uses that vary the parameters are therefore most attractive when several state draws or likelihood evaluations are taken at each parameter value, or when computations can be parallelized.
Let $\bm{\theta}$ collect the parameters of the state and observation equations, and write $\gamma_{\bm{\theta}}$, $q_{\bm{\theta}}$, $\widehat C_1(\bm{\theta})$ and $\ell_{t,\bm{\theta}}$ for the objects of Sections~\ref{s:method} and \ref{s:implementation} evaluated at $\bm{\theta}$. The normalizing constant
\[
Z(\bm{\theta})
=
\int_{\mathbb{R}^T}\gamma_{\bm{\theta}}(\mathbf{h})\,\text{d}\mathbf{h}
=
p(\mathbf{y}\,|\,\bm{\theta})
\]
is the observed-data likelihood, and $w_{\bm{\theta}}(\mathbf{h}) =\gamma_{\bm{\theta}}(\mathbf{h})/q_{\bm{\theta}}(\mathbf{h})$ satisfies $ 0<w_{\bm{\theta}}(\mathbf{h})\leqslant\widehat C_1(\bm{\theta})$ and $ \mathbb E_{q_{\bm{\theta}}} \{w_{\bm{\theta}}(\mathbf{h})\}
= Z(\bm{\theta})$ by Theorem~\ref{thm:exact} and Corollary~\ref{cor:is}. Write
$p_{\bm{\theta}}=Z(\bm{\theta})/\widehat C_1(\bm{\theta})$ for the corresponding
rejection acceptance probability, reserving $p(\bm{\theta})$ for the prior
density.
\subsection{Likelihood Evaluation} \label{ss:likeuse}
For $\mathbf{h}^{(1)},\ldots,\mathbf{h}^{(M)}$ drawn independently from $q_{\bm{\theta}}$, define
\begin{equation} \label{eq:Zhat}
\widehat Z_M(\bm{\theta})
=
\frac{1}{M}
\sum_{m=1}^{M}
w_{\bm{\theta}}\big(\mathbf{h}^{(m)}\big).
\end{equation}
The estimator is nonnegative and unbiased for $Z(\bm{\theta})$.
Corollary~\ref{cor:is} further gives
\[
\frac{\text{Var}\{\widehat Z_M(\bm{\theta})\}}
{Z(\bm{\theta})^2}
\leqslant
\frac{1/p_{\bm{\theta}}-1}{M}.
\]
Thus the likelihood estimator has finite variance at every parameter value, with a bound determined by a quantity that has a direct algorithmic interpretation. For example, an acceptance probability of $0.75$ implies relative variance at most $1/(3M)$. The bound is the exact relative variance of the cruder estimator $\widehat C_1(\bm{\theta})M^{-1}\sum_m\mathbf{1}\{\text{accept}^{(m)}\}$ constructed from the rejection indicators alone. Estimator~\eqref{eq:Zhat} is its Rao--Blackwellization \citep{CasellaRobert1996}, so the continuous weights should be retained rather than discarded at the rejection step.
Estimator~\eqref{eq:Zhat} can be used for simulated likelihood evaluation,
likelihood ratios and comparison of parameter values, in the manner of the
integrated likelihood estimators of \citet{CJ09}. The likelihood itself is
estimated without bias, although $\log\widehat Z_M(\bm{\theta})$ is not an
unbiased estimator of the log likelihood. The same estimator also supports the marginal likelihood calculation used in the empirical application in Section~\ref{s:app}. Let $p(\bm{\theta})$ be the prior density and let $g$ be an importance density whose support contains that of $Z(\bm{\theta})p(\bm{\theta})$. If $\bm{\theta}^{(1)},\ldots,\bm{\theta}^{(N)}$ are drawn from $g$, and an independent estimate $\widehat Z_M(\bm{\theta}^{(n)})$ is constructed at each draw, then
\begin{equation} \label{eq:marglike}
\widehat p(\mathbf{y})
=
\frac{1}{N}
\sum_{n=1}^{N}
\frac{p(\bm{\theta}^{(n)})}
{g(\bm{\theta}^{(n)})}
\widehat Z_M(\bm{\theta}^{(n)})
\end{equation}
is unbiased for $p(\mathbf{y})$. Bounded state weights rule out an infinite variance arising from integration over the high-dimensional state path at a fixed $\bm{\theta}$; the remaining tail requirement concerns the lower-dimensional importance density $g$. \ref{app:data} gives a sufficient condition for the choice used in the application.
More generally, because $\widehat Z_M(\bm{\theta})$ is nonnegative and unbiased, it can also be used in pseudo-marginal parameter samplers \citep{Beaumont2003,AndrieuRoberts2009,ADH10}; \citet{LiuPlagborgMoller2023} illustrate the reach of this idea in economics, basing full-information Bayesian inference for heterogeneous agent models with unobserved aggregate states on a numerically unbiased likelihood estimator.
\subsection{Benchmarking State Approximations} \label{ss:benchmark}
Exact simulation is especially useful for benchmarking approximate state smoothers. Assessing a particle smoother, a Laplace or variational approximation, or a Metropolis--Hastings state sampler ordinarily requires a reference that is itself approximate, so the comparison measures the difference between two approximation errors. Accepted paths from Algorithm~\ref{alg:main} instead provide draws from the intended smoothing distribution, leaving only ordinary Monte Carlo error. Posterior expectations of integrable functions can therefore be estimated from accepted paths without burn-in, thinning or autocorrelation correction. Alternatively, all proposal draws can be retained and used by ordinary or self-normalized importance sampling. The bounded weights prevent a small number of paths from dominating these estimates.
The rejection bound also gives a direct measure of the quality of the proposal
itself. Writing $\pi_{\bm{\theta}}(\mathbf{h}) = p(\mathbf{h}\,|\,\mathbf{y},\bm{\theta}),$ we have
$\pi_{\bm{\theta}}(\mathbf{h})/q_{\bm{\theta}}(\mathbf{h}) \leqslant 1/p_{\bm{\theta}},$ with equality attained. Write $D_{\infty}$, $\mathrm{KL}$ and $\chi^{2}$ for the R\'{e}nyi divergence of order infinity, the Kullback--Leibler divergence and the chi-squared divergence. We have $D_{\infty} (\pi_{\bm{\theta}}\Vert q_{\bm{\theta}}) = -\log p_{\bm{\theta}},$ $ \mathrm{KL} (\pi_{\bm{\theta}}\Vert q_{\bm{\theta}}) \leqslant -\log p_{\bm{\theta}},$ and $ \chi^2 (\pi_{\bm{\theta}}\Vert q_{\bm{\theta}}) \leqslant 1/p_{\bm{\theta}}-1.$ Section~\ref{s:MC} uses exact simulation in precisely this benchmarking role.
\subsection{Further Consequences} \label{ss:furtheruses}
Several other consequences follow from the same bound. Algorithm~\ref{alg:main} can be used as an exact blocked state update in a posterior sampler, drawing $\mathbf{h}$ from $p(\mathbf{h}\,|\,\mathbf{y},\bm{\theta})$ before updating $\bm{\theta}$. This removes approximation and within-block Markov-chain error from the state update, although dependence remains through the parameter draws. Likewise, an independence sampler based on $q_{\bm{\theta}}$ satisfies
$P_{\bm{\theta}}(\mathbf{h},A) \geqslant p_{\bm{\theta}}\pi_{\bm{\theta}}(A)$ and is therefore uniformly ergodic \citep{tierney94,MengersenTweedie1996}. These observations are theoretical consequences rather than recommended implementations, since each draws only one state path per parameter value.
Finally, although $p_{\bm{\theta}}$ is unknown, it is itself an acceptance
probability and can therefore be estimated from bounded random variables.
A lower confidence bound for $p_{\bm{\theta}}$ immediately gives corresponding
upper bounds on $1/p_{\bm{\theta}}-1$ and $-\log p_{\bm{\theta}}$. This is
qualitatively different from diagnosing the variance of unbounded importance
weights, which may be driven by rare events absent from a finite simulation.
\section{Relation to Existing State Simulation Methods} \label{s:related}
The proposal developed here shares its architecture with several existing simulation smoothers, but differs in what it asks of the approximation error. In twisted methods, each transition density is multiplied by a function of the current state and renormalized, giving an importance density of the form \eqref{eq:proposal}. Lemma~\ref{lem:sep} is a property of this architecture: apart from an additive constant, the log target-to-proposal ratio is a sum of univariate functions of the states. Efficient importance sampling chooses Gaussian twists by backward regressions on simulated draws \citep{RichardZhang2007, LMRD2013}; the HESSIAN method of \citet{McCausland2012}, confined like the present paper to univariate states, matches derivatives of the target conditional log density through fifth order; and \citet{ScharthKohn2016} combine global twisting with particle resampling. A related family instead approximates the smoothing density around its mode \citep{SP97, DurbinKoopman1997, JK08, MMP11, chan17}. These methods seek a target-to-proposal ratio that varies little; the tangent construction instead controls its sign.
That distinction is a trade-off. The construction here matches only the level and first derivative of the backward function at each node, giving an error of order $\Delta^2$ in the node spacing, coarser than a fifth-order fit or a variance-minimizing regression. In exchange, every $d_t$ is nonpositive, so Lemma~\ref{lem:sep} turns the global dominating constant into a computable quantity. Importance sampling requires only support coverage and integrable weights; bounded weights, as in \citet{McCausland2012}, are a stronger sufficient condition that also guarantees finite variance. Rejection sampling requires more: the value of a finite dominating constant must be known. A period-by-period bound does not solve this problem, because multiplying the $T$ local bounds can lead to exponential deterioration when their average log gap is bounded away from zero. The tangent construction instead makes every coordinatewise residual supremum zero. Its univariate antecedents are \citet{Devroye1986} and \citet{GilksWild1992}; separability carries the envelope to the full state path.
The ensemble rejection sampler of \citet*{DeligiannidisDoucetRubenthaler2020} reaches exactness by a different route and covers a broader class of models, including nonlinear and non-Gaussian transitions. It performs rejection sampling on an extended particle space and obtains its cost guarantee from two-sided bounds on incremental weights. Under their Proposition~4, keeping acceptance bounded away from zero requires an ensemble size of order $T$, implying
an $O(T^3)$ cost per exact draw. Those bounds do not apply in our setting because the state space is unbounded and the incremental weights have infimum zero; Theorem~\ref{thm:scaling} instead exploits concavity and Gaussian transitions to obtain the $T/G^2$ bound on the global envelope error.
Another class of methods approximates the measurement density rather than twisting the transitions. The auxiliary mixture sampler of \citet*{KSC98} for stochastic volatility replaces the measurement density by a finite mixture of normals, after which standard Gaussian simulation smoothers apply; extensions include leverage \citep{OCSN07}, stochastic volatility in mean \citep{HCO25}, and count and binomial models \citep{FruhwirthSchnatterWagner2006,FruhwirthSchnatterFruhwirthHeldRue2009}. These methods draw exactly from an approximating model, with any remaining discrepancy corrected, if desired, by reweighting. The construction here instead draws directly from the intended model under the log-concavity condition in Assumption~\ref{as:model}(ii).
The signed residual also distinguishes the method as an importance sampler. Corollary~\ref{cor:is} gives the explicit relative-variance bound $1/p-1$, which neither efficient importance sampling nor the HESSIAN method provides. Finite variance cannot in general be assumed for Gaussian proposals \citep{KoopmanShephardCreal2009}. For a Gaussian importance density with precision $\mathbf{Q}_q$, the second moment is infinite whenever $\mathbf{d}'(2\mathbf{Q}-\mathbf{Q}_q)\mathbf{d}<0$ for a nonnegative direction $\mathbf{d}$ satisfying the tail condition of Section~\ref{ss:gaussian}. When $\mathbf{Q}_q-\mathbf{Q}$ is diagonal and the state coefficients $\beta_t$ are nonnegative, this reduces to an eigenvalue check on a tridiagonal matrix; for the Gaussian approximation at the mode, $\mathbf{Q}_q$ is the matrix $\mathbf{K}$ of Section~\ref{ss:gaussian} and the difference is diagonal by construction. The criterion fails in all stochastic volatility experiments and all six empirical fits below. The tangent weights are bounded by construction.
\section{Monte Carlo Evidence} \label{s:MC}
This section examines three aspects of the proposed method. First, we assess whether exact rejection sampling remains practical as the sample size grows and across different models. Second, we compare the resulting likelihood estimator with existing simulation-based likelihood methods. Third, we use exact simulation to benchmark commonly used approximate smoothers.
For the stochastic volatility design, $\mu=0$, $\phi=0.95$ and $\sigma^2_h=0.01$ unless stated otherwise. Acceptance probabilities are estimated by the Rao--Blackwellized form of the realized acceptance rate, $\widehat p = M^{-1}\sum_{m=1}^{M}\text{e}^{-S_T^{(m)}}$, where $S_T=\sum_{t=1}^{T}\{\psi_t(h_t)-f_t(h_t)\}\geqslant0$ is the total envelope residual of a proposed path, as in Section~\ref{ss:scaling}, and $S_T^{(m)}$ is its value at the $m$th of $M$ proposals. All computations are in $\mathrm{M}\mathrm{{\scriptstyle ATLAB}}$ on an otherwise idle desktop machine. The comparison in Section~\ref{ss:MClike} is calibrated by wall clock: the number of draws or particles each method receives is read off a measured cost model. The exact efficiency ratios are therefore hardware-specific; the common-budget protocol, rather than the numerical ratios, is the reproducible object. We run four checks to validate the implementation, each targeting a property the construction guarantees in theory.\footnote{First, on a dense grid, the residual $d_t$ never exceeds $9.1\times10^{-13}$ for any of the eight models in Table~\ref{tab:robust}, so the envelope is respected to numerical precision. Second, the likelihood estimates agree with deterministic quadrature for $T\leqslant50$. Third, when $y_t=0$ for every $t$, the observation log density is affine, the smoothing distribution is Gaussian, and the acceptance probability is $1.0$. Finally, at $T=60$, smoothed moments and quantiles from $50{,}000$ draws agree to within $0.01$ with those from the Gaussian-envelope sampler of Section~\ref{ss:gaussian}, which is exact and computationally feasible at that sample size.}
\subsection{Acceptance Performance} \label{ss:MCscaling}
Figure~\ref{fig:scaling} shows that the tangent construction delivers high acceptance probabilities even for long state paths. Panel~(a) compares it with the Gaussian envelope of Section~\ref{ss:gaussian}. For the Gaussian envelope, the acceptance probability decays approximately exponentially with $T$, falling to $\text{e}^{-73.7}$ at $T=4{,}000$.\footnote{Probabilities this small are computed as follows: the two envelopes dominate the same unnormalized target, so the Gaussian acceptance probability equals $\widehat p\,\widehat C_1\text{e}^{-M_G}$ exactly, where $M_G$ is the Gaussian dominating constant of Section~\ref{ss:gaussian}, and only the factor $\widehat p$ is estimated. The same identity produces the Gaussian column of Table~\ref{tab:robust}.} With the tangent envelope and a fixed number of nodes, acceptance also declines with $T$, at a rate about one hundred times slower: at $T=4{,}000$ it is still $0.56$ with $G=81$. When instead $G$ is increased in proportion to $\sqrt{T}$, acceptance remains essentially constant at about $0.87$.
Panels~(b) and (c) show the source of this stability. The acceptance curves separate when plotted against $G$, but nearly coincide when plotted against $G/\sqrt{T}$. Thus the finite-sample behavior closely matches the $T/G^2$ envelope-error rate of Theorem~\ref{thm:scaling}: increasing the number of nodes at rate $\sqrt{T}$ prevents acceptance from deteriorating as the state dimension grows. Regressing $\log\mathbb E_q(S_T)$ on $\log T$ and $\log G$ gives exponents of $1.00$ and $-2.03$, and the companding grid of Theorem~\ref{thm:scaling} gives $0.97$ and $-2.01$ at a uniformly lower level of acceptance.
\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{fig1_scaling.pdf}
\caption{Acceptance, sample size and node refinement for the stochastic volatility model. Panel (a) reports $-\log(\text{acceptance})$. Panels (b) and (c) report acceptance against $G$ and $G/\sqrt{T}$, respectively. Each point in panels (b) and (c) is based on 600 proposals.}
\label{fig:scaling}
\end{figure}
The high acceptance rates are not specific to stochastic volatility. Table~\ref{tab:robust} reports results for eight models at $T=500$, sharing the Gaussian AR(1) state equation and differing only in the observation density. Acceptance for the proposed method ranges from $0.811$ to $0.931$. In contrast, acceptance for the Gaussian envelope ranges from $1.1\times10^{-4}$ to $4.9\times10^{-66}$. The results are also stable across independently generated data sets: the standard deviation of the acceptance probability over twenty replications, reported in Table~\ref{tab:robust}, is below $0.004$ in seven of the eight designs, and $0.021$ for the Poisson. Hence the practical advantage of the tangent construction extends across the observation densities considered, rather than being specific to a particular stochastic volatility sample.
\begin{table}[H]
\caption{Acceptance probabilities across eight models, $T=500$ and $G=81$. ``Across data sets'' is the standard deviation of acceptance over twenty independent data sets. The Student's $t$, duration and Poisson designs use $\sigma^2_h=0.05$, with $\mu=2$ for the Poisson, and the dynamic probit design $\phi=0.98$ and $\sigma^2_h=0.1$. The duration models have unit-mean exponential, gamma (shape $2$) and Weibull (shape $1.5$) errors. Every design uses the default node range $[-4,6]$ of Section~\ref{ss:nodes}, even where its orientation is unfavorable.}
\label{tab:robust}
\centering
\setlength{\tabcolsep}{4.5pt}
\begin{tabular}{lccc}
\hline\hline
Observation density & Tangent envelope & Across data sets & Gaussian envelope \\ \hline
Stochastic volatility & 0.931 & 0.0010 & $1.1\times10^{-4}$ \\
\rowcolor{lightgray}
Stochastic volatility in mean & 0.918 & 0.0013 & $3.5\times10^{-6}$ \\
Student's $t$ errors & 0.896 & 0.0006 & $4.3\times10^{-10}$ \\
\rowcolor{lightgray}
Dynamic probit & 0.855 & 0.0034 & $1.2\times10^{-19}$ \\
Duration, exponential errors & 0.867 & 0.0004 & $3.2\times10^{-20}$ \\
\rowcolor{lightgray}
Duration, gamma errors & 0.851 & 0.0004 & $2.6\times10^{-30}$ \\
Duration, Weibull errors & 0.850 & 0.0004 & $7.7\times10^{-32}$ \\
\rowcolor{lightgray}
Poisson counts & 0.811 & 0.0206 & $4.9\times10^{-66}$ \\
\hline\hline
\end{tabular}
\end{table}
Acceptance also remains high when the stochastic volatility design is made more difficult: it is $0.885$ at $\phi=0.98$ and $0.855$ at $\phi=0.995$ as the state approaches a unit root, and $0.869$ and $0.835$ at $\sigma_h^2=0.1$ and $0.5$ as the innovation variance grows, while setting fifty observations exactly to zero, or inflating twenty observations by a factor of twenty, leaves it essentially unchanged.
\subsection{Likelihood Evaluation} \label{ss:MClike}
We next assess the tangent-twisted proposal for likelihood evaluation. We compare its unbiased likelihood estimator with efficient importance sampling (EIS), a Laplace--Gaussian importance sampler based on the Gaussian approximation at the posterior mode, and a bootstrap particle filter. The experiment uses $T=500$, $1{,}000$, $2{,}000$ and $4{,}000$. For each $T$, we generate five independent series and conduct four independent likelihood evaluations per method and series. Efficiency ratios are computed within each series and then averaged across series.
For likelihood evaluation we fix $G=121$ at every sample size rather than use the acceptance-based $G\propto\sqrt{T}$ rule of Section~\ref{ss:scaling}, since the relevant criterion is weight variance per unit of computing time. Figure~\ref{fig:like} compares the methods under a common wall-clock budget equal to the cost of one complete tangent-twisted evaluation with two hundred draws at a new parameter vector, including method-specific setup.\footnote{For EIS, setup consists of five backward regression sweeps using one thousand draws each; the Gaussian methods include the posterior-mode calculation. Comparator draw or particle counts are chosen to exhaust the same wall-clock budget. For the tangent-twisted and EIS estimators, the quantities in Figure~\ref{fig:like} are constructed from the relative variance $\mathrm{rv}_w$ of a single importance weight, estimated from separate long runs; panel~(a) uses the lognormal moment-matched value $\{\log(1+\mathrm{rv}_w/N)\}^{1/2}$. The particle-filter relative variance is estimated from independent replicated likelihood evaluations, since its likelihood estimate is a product of sequential averages rather than an average of independent path weights.}
\begin{figure}[H]
\centering
\includegraphics[width=0.92\textwidth]{fig3_likelihood.pdf}
\caption{Likelihood estimation for the stochastic volatility model at a fixed parameter vector under a common computing-time budget, with $G=121$. Each point averages over five independent series and four evaluations per method and series. The budget equals the cost of a tangent-twisted evaluation with two hundred draws, including setup; comparator draw or particle counts match it. Panel~(a) reports the implied standard deviation of $\log\widehat Z$, and panel~(b) the relative variance of $\widehat Z$ multiplied by computing time. Lower values indicate greater efficiency. The Laplace--Gaussian importance sampler is omitted because its weights have infinite variance in every design.}
\label{fig:like}
\end{figure}
The tangent-twisted estimator is substantially more efficient than the particle filter, by a factor of $1{,}900$ to $3{,}500$, and is comparable to EIS. The tangent proposal is more costly per draw because each transition is a mixture over $G$ pieces, but its lower weight variance largely offsets that cost, with neither method uniformly dominating across series and sample sizes. What distinguishes the tangent estimator is therefore not a systematic efficiency advantage over EIS, but its guarantee: its weights are bounded by a constant the algorithm computes, so Corollary~\ref{cor:is} ensures finite relative variance by construction. Consistent with this bound, the largest tangent weight exceeds the mean by less than $1.5\%$ at every sample size, compared with a factor of $55$ for EIS at $T=4{,}000$.
The Laplace--Gaussian comparison illustrates why this distinction matters. The finite-variance criterion of Section~\ref{s:related} fails for all twenty data sets, so neither quantity in Figure~\ref{fig:like} is defined for this proposal.\footnote{The smallest eigenvalue of $2\mathbf{Q}-\mathbf{K}$ ranges from $-0.41$ to $-0.26$. At $T=4{,}000$, the largest observed Laplace--Gaussian weight exceeds its mean by more than two orders of magnitude, against $1.013$ for the tangent proposal, and it grows with the number of draws, as the maximum of an infinite-variance quantity must. The same finite-variance criterion fails at all six empirical fits in Section~\ref{s:app}.} The tangent weights, by contrast, are bounded by construction.
\subsection{Benchmarking Approximate Smoothers} \label{ss:MCaccuracy}
A separate use of the method is to benchmark approximate state smoothers. Such assessments ordinarily compare one approximation with another; the tangent construction instead provides a reference for the intended smoothing distribution. We compare the Gaussian approximation at the posterior mode, EIS with self-normalized importance weights, and the genealogy smoother associated with a bootstrap particle filter. Each simulation-based method uses $8{,}000$ draws or particles, while the Gaussian approximation is evaluated in closed form. For greater precision, the reference uses $8{,}000$ tangent draws with self-normalized weights, whose effective sample size equals the nominal one to three decimal places.\footnote{\label{fn:refcheck}The reference is checked in three ways. First, the bound of Section~\ref{s:uses} gives $\chi^2(\pi_{\bm{\theta}}\Vert q_{\bm{\theta}})\leqslant1/p-1\leqslant0.339$ for $T\leqslant2{,}000$. Second, the largest realized weight exceeds the mean by at most $2.6\%$. Finally, the reference agrees with $4{,}000$ exact draws from Algorithm~\ref{alg:main}, with the largest absolute discrepancy in the smoothed mean ranging from $0.058$ to $0.070$ reference standard deviations. For comparison, two independent reference replicates differ by $0.045$ to $0.060$ under the same maximum-over-$t$ metric.} For each $t$, we examine the posterior mean and the $5\%$ quantile of $h_t$, with errors normalized by the reference marginal posterior standard deviation.
Table~\ref{tab:accuracy} shows three distinct patterns. EIS is extremely accurate: its root-mean-square error in the smoothed mean is $0.014$ to $0.017$ reference standard deviations, and its largest absolute error is $0.045$ to $0.058$, comparable to the Monte Carlo variation between independent reference replicates. Thus the case for the tangent construction is not that EIS gives materially inaccurate state estimates, but that it supplies exact draws and a computable global bound on the importance weights.
The Gaussian approximation exhibits a small but persistent approximation bias. Its root-mean-square error is about $0.08$ posterior standard deviations for the smoothed mean and $0.09$ for the lower quantile, with little change as $T$ increases. Because the Gaussian functionals are evaluated in closed form, these discrepancies reflect approximation rather than simulation error.
The genealogy smoother behaves differently: its error increases with $T$ and is concentrated near the beginning of the path, where particle ancestry has collapsed. At $T=2{,}000$, only about one percent of the original particles remain distinct over the first five percent of dates, and errors in the lower tail reach $0.72$ posterior standard deviations. The terminal filtering weights remain well behaved, with effective sample size at least $0.98$ of the nominal size, so the deterioration comes from genealogical collapse rather than poor filtering. This illustrates the benchmarking value of the exact sampler: approximation errors that are otherwise difficult to separate become directly measurable.
\begin{table}[H]
\caption{Error in the smoothed posterior mean and in the $5\%$ quantile of $h_t$, relative to the reference, in units of the reference marginal posterior standard deviation. RMSE is the root-mean-square error and Max the largest absolute error, both over $t$.}
\label{tab:accuracy}
\centering
\begin{tabular}{lrrrrrr}
\hline\hline
& \multicolumn{2}{c}{$T=500$} & \multicolumn{2}{c}{$T=1{,}000$} & \multicolumn{2}{c}{$T=2{,}000$} \\
& RMSE & Max & RMSE & Max & RMSE & Max \\ \hline
\multicolumn{7}{l}{\textit{Posterior mean}} \\
\rowcolor{lightgray}
Gaussian approximation & 0.075 & 0.109 & 0.079 & 0.116 & 0.078 & 0.123 \\
Efficient importance sampling & 0.016 & 0.046 & 0.014 & 0.045 & 0.017 & 0.058 \\
\rowcolor{lightgray}
Genealogy smoother & 0.051 & 0.158 & 0.082 & 0.240 & 0.094 & 0.394 \\
\multicolumn{7}{l}{\textit{$5\%$ quantile}} \\
\rowcolor{lightgray}
Gaussian approximation & 0.092 & 0.176 & 0.093 & 0.154 & 0.095 & 0.173 \\
Efficient importance sampling & 0.033 & 0.090 & 0.031 & 0.110 & 0.033 & 0.111 \\
\rowcolor{lightgray}
Genealogy smoother & 0.114 & 0.349 & 0.156 & 0.588 & 0.194 & 0.720 \\
\hline\hline
\end{tabular}
\end{table}
\section{Application: Stochastic Conditional Duration} \label{s:app}
We illustrate the method using the stochastic conditional duration model of \citet{BauwensVeredas2004}, providing an application outside stochastic volatility in which the choice of conditional duration distribution is an economically relevant model comparison. \citet{BauwensGalli2009} develop likelihood evaluation for this model using efficient importance sampling.
\subsection{Model Specifications and Data}
The stochastic conditional duration model is $y_t=\text{e}^{h_t}\varepsilon_t,$
where $\text{e}^{h_t}$ is the latent expected duration, $h_t$ follows the affine Gaussian process in \eqref{eq:model}, and $\varepsilon_t$ is a positive unit-mean error. We compare exponential, gamma and Weibull errors. The unit-mean normalization identifies the level of the state. Writing $u$ for the state, the corresponding observation log-densities are
\[
\ell_t(u)=-u-y_t\text{e}^{-u}, \qquad
\ell_t(u)=-au-ay_t\text{e}^{-u}, \qquad
\ell_t(u)=-ku-(y_t/\lambda)^k\text{e}^{-ku},
\]
for the exponential, gamma with shape $a$, and Weibull with shape $k$, respectively, where $\lambda=1/\Gamma(1+1/k)$, up to terms independent of $u$. All three are concave, so Algorithm~\ref{alg:main} applies directly, and the gamma and Weibull families nest the exponential at unit shape.
We use Binance trade records for LINKUSDT and ALGOUSDT and construct seasonally adjusted volume durations for 15 May 2024; \ref{app:data} gives the data sources and construction details. Table~\ref{tab:appdata} shows substantial persistence and overdispersion in both series. The standard deviation exceeds the mean in each, with the greater dispersion for ALGOUSDT.
\begin{table}[H]
\caption{Adjusted volume durations, 15 May 2024. The seasonal range is the ratio of the largest to the smallest value of the periodic factor $\varphi_t$. Autocorrelations are of $\log y_t$.}
\label{tab:appdata}
\centering
\begin{tabular}{lrrrrrr}
\hline\hline
Pair & $T$ & Mean (s) & sd/mean & $\rho_1$ & $\rho_{20}$ & Seasonal range \\ \hline
LINKUSDT & 2{,}711 & 33.5 & 1.18 & 0.52 & 0.10 & 2.5 \\
\rowcolor{lightgray}
ALGOUSDT & 2{,}151 & 44.9 & 1.48 & 0.59 & 0.29 & 4.1 \\
\hline\hline
\end{tabular}
\end{table}
\subsection{Estimates and Model Comparison} \label{ss:appresults}
Each model is estimated by maximizing the unbiased likelihood estimator \eqref{eq:Zhat} using common random numbers; \ref{app:data} gives the implementation details. Table~\ref{tab:appfit} reports the estimates. Acceptance ranges from $0.78$ to $0.88$ across the six specifications. At all six estimates, the Gaussian approximation fails the finite-variance criterion of Section~\ref{s:related}, whereas the tangent weights remain bounded.
The model ranking is the same in both series: gamma is preferred to Weibull, which is preferred to exponential. The estimated gamma shapes, $0.642$ and $0.455$, are both well below the exponential benchmark of one and decline as the dispersion in Table~\ref{tab:appdata} rises. The error specification also materially affects the estimated state process: imposing exponential errors produces larger innovation variances and lower persistence, because variation that can instead be captured by the error distribution is forced into the latent state.
\begin{table}[H]
\caption{Simulated maximum likelihood estimates and the Bayesian information criterion (BIC). The exponential family has no shape parameter. Acceptance is the probability of Theorem~\ref{thm:exact} evaluated at the estimates with $G=188$.}
\label{tab:appfit}
\centering
\begin{tabular}{llrrrrrr}
\hline\hline
Pair & Family & $\mu$ & $\phi$ & $\sigma^2_h$ & Shape & Acceptance & BIC \\ \hline
LINKUSDT & Gamma & 3.21 & 0.991 & 0.015 & 0.642 & 0.856 & 22{,}741 \\
\rowcolor{lightgray}
& Weibull & 3.22 & 0.988 & 0.023 & 0.813 & 0.856 & 22{,}970 \\
& Exponential & 3.14 & 0.978 & 0.046 &---& 0.852 & 23{,}101 \\
\rowcolor{lightgray}
ALGOUSDT & Gamma & 3.21 & 0.981 & 0.072 & 0.455 & 0.883 & 17{,}167 \\
& Weibull & 3.24 & 0.969 & 0.175 & 0.606 & 0.880 & 17{,}367 \\
\rowcolor{lightgray}
& Exponential & 2.53 & 0.830 & 2.081 &---& 0.777 & 17{,}651 \\
\hline\hline
\end{tabular}
\end{table}
Because \eqref{eq:Zhat} is nonnegative and unbiased with bounded state weights, it can also be used in the marginal likelihood estimator \eqref{eq:marglike}. We use the priors
$\mu\sim\mathcal{N}(0,10)$, $(\phi+1)/2\sim\mathcal{B}(20,1.5)$, $\sigma_h^2\sim\mathcal{IG}(2.5,0.1)$,
and a standard normal prior for the log shape parameter. The parameter importance density is a heavy-tailed Student's $t$ centered at the simulated maximum likelihood estimate; \ref{app:data} gives the remaining details. Table~\ref{tab:appml} confirms the BIC ranking by a wide margin. Relative to gamma, the log marginal likelihood is lower by $101$--$114$ units for Weibull and by $178$--$249$ units for exponential. These differences are far larger than the numerical standard errors, which are at most $0.06$.
\begin{table}[H]
\caption{Log marginal likelihoods from \eqref{eq:marglike}, with numerical standard errors from $2{,}000$ parameter draws.}
\label{tab:appml}
\centering
\begin{tabular}{llrr}
\hline\hline
Pair & Family & $\log\widehat p(\mathbf{y})$ & NSE \\ \hline
LINKUSDT & Gamma & $-11{,}367.9$ & 0.06 \\
\rowcolor{lightgray}
& Weibull & $-11{,}481.7$ & 0.05 \\
& Exponential & $-11{,}546.4$ & 0.03 \\
\rowcolor{lightgray}
ALGOUSDT & Gamma & $-8{,}579.6$ & 0.05 \\
& Weibull & $-8{,}680.3$ & 0.05 \\
\rowcolor{lightgray}
& Exponential & $-8{,}828.9$ & 0.04 \\
\hline\hline
\end{tabular}
\end{table}
\section{Concluding Remarks and Future Research} \label{s:conclusion}
This paper develops an exact rejection sampler for non-Gaussian state space models with a scalar latent state, affine Gaussian dynamics, and a concave observation log-density. Backward twisting makes the log target-to-proposal ratio separable across dates, while tangent envelopes make every coordinatewise residual nonpositive and attain zero at the nodes. The resulting global dominating constant is therefore computable, attained, and smallest admissible for the proposal. Rejection sampling then produces independent draws from the intended joint smoothing distribution without burn-in, Markov-chain correction, or approximation-model error.
Exactness would be of limited practical value if acceptance deteriorated rapidly with the length of the state path. For the companding grid, the accumulated envelope error is \(O(T/G^2)\), so increasing the number of nodes at rate \(G\propto\sqrt{T}\) keeps the acceptance probability bounded away from zero. For stochastic volatility, the conditions underlying this result hold almost surely, while the simpler mode-centered grid used in practice displays the same scaling empirically. Importantly, exactness holds for every finite node set: refinement governs computational efficiency rather than the validity of accepted draws. The same proposal also yields a nonnegative unbiased likelihood estimator with bounded importance weights and therefore finite variance. In the numerical experiments, it is competitive with efficient importance sampling and provides an exact benchmark for assessing approximate state smoothers.
The main restriction of the present construction is the scalar latent state. For vector-valued states, the domination argument continues to apply, but the affine regions of the tangent envelope become polyhedra, making Gaussian probability evaluation and truncated-normal simulation the principal computational challenges. Developing practical methods for low-dimensional vector states is therefore a natural next step. Another direction is adaptive node placement that reduces the cost of the backward recursion while preserving the global bound. More broadly, the results suggest that one-sided approximations can be useful in latent-variable computation when a computable guarantee is more valuable than minimizing an unsigned approximation error.
\bigskip
\singlespacing
\bibliographystyle{econometrica}
\bibliography{var,exact-AR}
\onehalfspacing
\pagebreak