# Particle-Based Inference for Continuous-Discrete State Space Models

Christopher Stanton
Department of Statistical Science, University College London
and
Alexandros Beskos
Department of Statistical Science, University College London

###### Abstract

This article develops a methodology allowing application of the complete
machinery of particle-based inference methods upon what we call the class of
_continuous-discrete State Space Models_ (CD-SSMs). Such models correspond to
a latent continuous-time Itô diffusion process which is observed with noise at
discrete time instances. Due to the continuous-time nature of the hidden
signal, standard Feynman-Kac formulations and their accompanying particle-
based approximations have to overcome several challenges, arising mainly due
to the following considerations: (i) finite-time transition densities of the
signal are typically intractable; (ii) ancestors of sampled signals are
determined w.p. 1, thus cannot be resampled; (iii) diffusivity parameters
given a sampled signal yield Dirac distributions. We overcome all above issues
by introducing a framework based on carefully designed proposals and
transformations thereof. That is, we obtain new expressions for the Feynman-
Kac model that accommodate the effects of a continuous-time signal and
overcome induced degeneracies. The constructed formulations will enable use of
the full range of particle-based algorithms for CD-SSMs: for
filtering/smoothing and parameter inference, whether online or offline. Our
framework is compatible with guided proposals in the filtering steps that are
essential for efficient algorithmic performance in the presence of informative
observations or in higher dimensions, and is applicable for a very general
class of CD-SSMs, including the case when the signal is modelled as a hypo-
elliptic diffusion. Our methods can be immediately incorporated to available
software packages for particle-based algorithms.

Keywords: Feynman-Kac Formula, Particle Gibbs, Hypo-Elliptic Diffusion, Guided
Proposal.

##  1 Introduction

This work considers the development of general particle-based Monte Carlo
methods for statistical inference within the broad class of what we call
_continuous-discrete State Space Models_ (CD-SSMs). For such a CD-SSM, the
latent process is a continuous-time diffusion process (Oksendal,, 2013)
defined as the solution of a Stochastic Differential Equation (SDE):

$dX(s)=b(s,X(s))ds+\sigma(s,X(s))dB(s),\qquad X(0)=x(0)\in\mathbb{R}^{d},$ (1)

for a $d_{w}$-dimensional Brownian motion $B=\\{B(s)\\}_{s\geq 0}$, a drift
function $b:[0,\infty)\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ and a diffusion
coefficient function
$\sigma:[0,\infty)\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d_{w}}$, with
$d\geq 1$, $1\leq d_{w}\leq d$. The signal process $X=\\{X(s)\\}_{s\geq 0}$ is
assumed to be indirectly observed at $T\geq 1$ discrete times
$0=s_{0}<s_{1}<\dots<s_{T}$, with data $Y_{1:T}=(Y_{1},\ldots,Y_{T})$ such
that:

$Y_{t}\,\big{|}\,\\{X(s)\\}_{s\geq 0},Y_{1:t-1},Y_{t+1:T}\sim
f_{t}(Y_{t}|X(s_{t})),\qquad t\in\\{1,\dots,T\\},$ (2)

independently over $t=1,\ldots,T$, for a conditional probability density
function (pdf) $f_{t}(\cdot|\cdot)$. We define $\mathcal{T}=\\{1,\ldots,T\\}$
and $\Delta_{t}=s_{t}-s_{t-1}$, $t\in\mathcal{T}$. For convenience, we often
write $X([u,v])$ for $\\{X(s)\\}_{s\in[u,v]}$, with $0\leq u<v\leq s_{T}$. We
will take under consideration the full range of particle-based algorithms
related with SSMs developed in the literature, including ones that involve
filtering, smoothing and parameter estimation (Chopin and Papaspiliopoulos,,
2020). For instance, the filtering problem, expressed within the CD-SSM
setting, considers inference about $X(s_{t})|Y_{1:t}$, for each $t\geq 1$, via
an online approach. The smoothing problem can be considered as either the
offline problem of inferring the joint distribution $X([0,s_{T}])|Y_{1:T}$, or
as the online problem of inferring $X([0,s_{t}])|Y_{1:t}$ for each
$t\in\mathcal{T}$. We note here that the name ‘continuous-discrete SSM’ we use
here resembles the terminology in Mider et al., (2021).

A standard Sequential Monte Carlo (SMC) algorithm (namely, a particle filter)
is immediately available for the approximation of the smoothing distribution
$X([0,s_{T}])|Y_{1:T}$ whereby iteratively for $t\in\mathcal{T}$:

  * (i)

particles are sampled in-between observation instances $s_{t-1}$, $s_{t}$,
according to some proposal dynamics
$M\big{(}dX([s_{t-1},s_{t}])\big{|}X(s_{t-1})\big{)}$;

  * (ii)

sampled particles representing $X([s_{t-1},s_{t}])$ are assigned un-normalised
weights according to the potential:

$\displaystyle\frac{dP(\cdot|X(s_{t-1}))}{dM(\cdot|X(s_{t-1}))}\big{(}X([s_{t-1},s_{t}])\big{)}\cdot
f_{t}(Y_{t}|X(s_{t}));$ (3)

  * (iii)

the weighted particles are resampled.

Here $P(\cdot|\cdot)$ is the probability law corresponding to the dynamics of
the signal SDE (1). Improved SMC algorithms provably overcome path degeneracy
characterising standard particle filters, by, e.g., updating the ancestors of
obtained particles (Godsill et al.,, 2004; Lindsten et al.,, 2014). Separate
particle-based algorithms are used in the setting when a parameter vector
$\theta\in\Theta\subseteq\mathbb{R}^{p}$, $p\geq 1$, must also be inferred.

A number of challenges arise due to the continuous-time nature of the signal
in the CD-SSM model class of interest in this work. The main ones are
summarised below:

  * (C.i)

First, sampling from $X(s_{t})|X(s_{t-1})$ is possible via approximate
numerical schemes, immediately enabling the use of ‘blind’ proposals in a
setting that ignores the continuous-time nature of the signal and only
considers its values at the discrete-time observation instances. Non-trivial
methods of superior performance do require access to the signal transition
density. However, for CD-SSMs, the distribution $P(X(s_{t})\in
dx^{\prime}|X(s_{t-1})=x)$ is typically intractable. As already implied by the
particle filter briefly described via steps (i)-(iii) in the previous
paragraph, a solution to this problem is to augment the state space by
considering the full continuous-time path of the SDE in-between observations.
Notice that expression (3) does not require a (finite-time) transition
density, and is analytically available via Girsanov’s theorem (Oksendal,,
2013).

  * (C.ii)

Then, a number of particle-based methods developed for smoothing and parameter
inference involve updating the ancestors of particles to overcome path
degeneracy, e.g. the Forward-Filtering, Backward-Sampling (FFBS) algorithm of
Godsill et al., (2004). However, owning to the continuous-time nature of the
signal, the ancestor of path $X([s_{t-1},s_{t}])$, $t\geq 2$, can only involve
particles $X([s_{t-2},s_{t-1}])$ for which:

$\displaystyle\texttt{end}(X([s_{t-2},s_{t-1}]))\equiv\texttt{start}(X([s_{t-1},s_{t}])).$
(4)

Thus, the overall method will be degenerate. From a probabilistic perspective,
this issue arises because the dominating measure of the path-valued proposal
for $X([s_{t-1},s_{t}])$ depends intrinsically on $X(s_{t-1})$. For the
algorithm to be allowed to select any ancestral particle
$X([s_{t-2},s_{t-1}])$, the path proposal for $X([s_{t-1},s_{t}])$ – or, more
precisely, of a transform of thereof – must have a density w.r.t. a reference
measure that will not depend on the previous state via $X(s_{t-1})$. This is
infeasible when working with $X([s_{t-1},s_{t}])$, but can be made possible
via careful transforms of the signal process. Indeed, such a
reparameterisation gives rise to a new CD-SSM formulation, for which both
target and proposals have densities w.r.t. the common reference measure
$\mathrm{Leb}\otimes\mathbb{W}_{t}$, with $\mathrm{Leb}(\cdot)$ denoting the
Lebesgue measure of dimension $d$ and $\mathbb{W}_{t}[\cdot]$ the distribution
of the standard $d_{w}$-dimensional Wiener process (to be fully specified in
the sequel). The construction has connections with the transforms of Chib et
al., (2004); Golightly and Wilkinson, (2008) used in the context of MCMC for
SDE models.

  * (C.iii)

The above two challenges are also present when one considers parameter
inference for CD-SSMs via application of particle-based MCMC methods (e.g.
particle MCMC, particle Gibbs (Andrieu et al.,, 2010)). Advanced versions of
particle Gibbs based on ancestral sampling improve mixing and have provably
robust properties w.r.t. $T$, see e.g. Godsill et al., (2004); Lindsten et
al., (2014). Implementation of such sampling approaches in the CD-SSM setting
immediately connects with the challenge raised in point (ii) above. In
addition, in the CD-SSM setting, the full conditional of
$\theta|\\{X(s)\\}_{s\in[0,s_{T}]}$ will be a Dirac measure for components of
$\theta$ appearing in the diffusion coefficient $\sigma(s,X(s))$. This is a
standard concern for data augmentation Bayesian methods for discretely
observed SDEs (Roberts and Stramer,, 2001).

In this contribution we take under consideration the full spectrum of
particle-based algorithms developed for standard SSMs, and produce a
methodology that overcomes challenges (C.i)-(C.iii) summarised above, thus
making such algorithms available for a large class of CD-SSMs. Our work
achieves that by developing careful transformations of the signal and
corresponding re-expressions of Feynman-Kac formulae for CD-SSMs.

We briefly summarise the main contributions of our work as follows.

  * •

The developed Feynman-Kac formulae and their Monte-Carlo approximation lead to
algorithms that bypass challenges (C.i)-(C.iii) for CD-SSMs. In brief, we
present two path proposals (labelled ‘forward’ and ‘backward’) and
corresponding transforms for CD-SSMs, and write down the deduced re-
expressions of the Feynman-Kac models. Though seemingly closely connected, we
illustrate that the produced forward method (FM) has subtle differences
compared to the backward one (BM). In particular, FM allows the application of
the developed methodology to the important class of hypo-elliptic SDEs.

  * •

Related to the above, we obtain re-expressions of the Feynman-Kac formulae,
appropriate for CD-SSMs. The full list of particle-based methods then becomes
immediately available for CD-SSMs. The relevant proposals and weights can be
seamlessly embedded within available probabilistic programming software. We
show a number of numerical applications, including FFBS and particle Gibbs
implementations.

We mention that instead of particle-based methods, standard MCMC algorithms
have been successfully applied for the class of CD-SSMs as defined herein, see
e.g. Chib et al., (2004); Golightly and Wilkinson, (2008) for some of the
first developments. If the logarithm of the posterior of parameters and latent
variable is differentiable, one can carry out inference for CD-SSMs with
standard MCMC (not particle-based) via use of available probabilistic
programming software, e.g. Stan (Carpenter et al.,, 2017). We stress that the
signal SDEs (elliptic and hypo-elliptic) upon our methodology develops in this
work are meant to serve as a blueprint for the treatment of more general
continuous-time model classes, including ones with discontinuous paths, i.e.
SDEs with jumps or discrete-valued pure jump-processes, when standard MCMC
methods are either inappropriate (e.g., derivative-driven algorithms,
including HMC used in Stan) or require complex model-specific development of
proposals over ‘blocks’ of signal paths, see e.g. Johannes and Polson, (2010)
for the case of SDEs with jumps and the references therein. Indicatively, one
of the authors has been involved in Yonekura and Beskos, (2022), where indeed
elliptic SDEs with jumps are considered, and a particle-based algorithm for
online parameter inference is proposed. Though that work lacks the general
scope of the present one, it does serve to illustrate the potential for the
full machinery of particle-based methods to be made available for general
continuous-time signals with discontinuities.

Another key advantage of particle-based methods is their scalability w.r.t.
the time index $T$. Karjalainen et al., (2023) prove that a particle Gibbs
method that makes use of FFBS to update the signal process allows for a number
of particles that scales as $\mathcal{O}(1)$ w.r.t. $T$, and provides mixing
time of $\mathcal{O}(\log T)$. Thus, the overall cost of the algorithm is
$\mathcal{O}(T\log T)$ w.r.t. $T$ (as $\mathcal{O}(T)$ is the computing cost
per step), with the number of particles selected independently of $T$. The
above contrast with reported results for standard MCMC, which point at ideal
costs of $\mathcal{O}(T^{2})$, $\mathcal{O}(T^{4/3})$ and
$\mathcal{O}(T^{5/4})$ for RWM, MALA and HMC, respectively (Gelman et al.,,
1997; Roberts and Rosenthal,, 1998; Beskos et al.,, 2013). Such robustness of
particle-based algorithms is achieved precisely due to the fact that careful
alteration of ancestors for particles within an FFBS (or similar) method
overcomes the path degeneracy of standard particle filter. Our contribution in
this work makes such ancestor alteration possible for the class CD-SSMs.

The remainder of the paper is structured as follows. In Section 2, we recall
the Feynman-Kac formulism of a SSM. In Section 3 we embed the CD-SSM within
the Feynman-Kac framework, in its original form, before consideration of path
transforms. In preparation of what follows, we describe forward and backward
decompositions of SDE paths. In Section 4 we discuss forward/backward path
proposals. In Section 5 we present the full methodology for CD-SSMs, by
describing the two path transforms and the corresponding Monte-Carlo methods,
FM and BM. We give concluding remarks and outline further research directions
in Section 6.

##  2 Feynman-Kac Model

Given measurable spaces $(\mathsf{X}_{t},\mathscr{X}_{t})$, $t\in\mathcal{T}$,
a Feynman-Kac model (Del Moral,, 2004, 2013) can be defined through an initial
law $M_{1}[dx_{1}]$, the Markov probability kernels $M_{t}[x_{t-1},dx_{t}]$,
$t\geq 2$, and the potential functions
$G_{1}:\mathsf{X}_{1}\rightarrow[0,\infty)$,
$G_{t}:\mathsf{X}_{t-1}\times\mathsf{X}_{t}\rightarrow[0,\infty)$, $t\geq 2$.
One then obtains, via a change of measure, the following probability law on
$\big{(}\mathsf{X}_{1}\times\cdots\times\mathsf{X}_{t},\mathscr{X}_{1}\otimes\cdots\otimes\mathscr{X}_{t}\big{)}$:

$\displaystyle\mathbb{Q}_{t}[dx_{1:t}]=\frac{1}{L_{t}}\Big{\\{}G_{1}(x_{1})\prod_{i=2}^{t}G_{i}(x_{i-1},x_{i})\Big{\\}}\times\mathbb{M}_{t}[dx_{1:t}],\qquad
t\in\mathcal{T},$ (5)

where we have set:

$\displaystyle\mathbb{M}_{t}[dx_{1:t}]=M_{1}[dx_{1}]\bigotimes_{i=2}^{t}M_{i}[x_{i-1},dx_{i}],$

and we have assumed that the normalising constant $L_{t}>0$ is well-defined,
that is:

$L_{t}=\mathbb{E}_{\mathbb{M}_{t}}\Big{[}G_{1}(X_{1})\prod_{i=2}^{t}G_{i}(X_{i-1},X_{i})\Big{]}<\infty.$

It is of interest to approximate integrals w.r.t. the distribution
$\mathbb{Q}_{t}[dx_{1:t}]$, In particular, $\mathbb{Q}_{t}[dx_{1:t}]$
corresponds to the _smoothing distribution_ of a SSM at time
$t\in\mathcal{T}$. We will also work with the _filtering distribution_ of a
SSM at $t\in\mathcal{T}$, i.e. the marginal of $\mathbb{Q}_{t}$ at $x_{t}$,
which we denote ${\mathbb{Q}}_{t}[dx_{t}]$. We use the notation
${\mathbb{Q}}_{t}[dx_{i:j}]$, $1\leq i\leq j\leq t$, to refer to general
marginal laws of $\mathbb{Q}_{t}[dx_{1:t}]$.

Successive filtering laws ${\mathbb{Q}}_{t}[dx_{t}]$ are related via the
following recursion:

$\displaystyle\mathbb{Q}_{t}[dx_{t-1:t}]$
$\displaystyle=const.\times\big{(}\mathbb{Q}_{t-1}[dx_{t-1}]\otimes
M_{t}[x_{t-1},dx_{t}]\big{)}\times G_{t}(x_{t-1},x_{t}).$ (6)

That is, starting from $\mathbb{Q}_{t-1}[dx_{t-1}]$, one obtains the one-step-
ahead predictive distribution, and the applies a change of measure, to arrive
at the next filtering law $\mathbb{Q}_{t}[dx_{t}]$. Here, we make use of the
convention $M_{1}[x_{0},dx_{1}]\equiv M_{1}[dx_{1}]$ and
$G_{1}(x_{0},x_{1})\equiv G_{1}(x_{1})$. The same convention is also applied
for other transition kernels, potentials and related quantities we use in the
sequel. Recursion (6) is analytically intractable except for the case of
finite state spaces, or when the Feynman-Kac model corresponds to a linear,
Gaussian SSM (LG-SSM). The particle filter (Gordon et al.,, 1993) provides an
effective Monte Carlo approximation of the recursion and, in particular, of
the filtering laws ${\mathbb{Q}}_{t}[dx_{t}]$, $t\in\mathcal{T}$.

The Feynman-Kac distribution ${\mathbb{Q}}_{T}[dx_{1:T}]$ is a Markov measure,
as it gives rise to the Markov chain $X_{1:T}\sim{\mathbb{Q}}_{T}[dx_{1:T}]$.
Assume that (with some abuse of notation)
$M_{t}[x_{t-1},dx_{t}]=M_{t}(x_{t}|x_{t-1})\nu{(dx_{t})}$, where
$M_{t}(\cdot|\cdot)$ is a transition density and $\nu(dx_{t})$ is a
$\sigma$-finite dominating measure, common for all $t\in\mathcal{T}$. FFBS
(Godsill et al.,, 2004) is an offline smoothing algorithm, used in this work
as representative of particle-based methods of improved efficiency due to
well-designed alterations in ancestors of particles during execution. Starting
with a forward pass of the particle filter that delivers Monte Carlo
approximations of the filtering laws $\mathbb{Q}_{t}[dx_{t}]$,
$t\in\mathcal{T}$, FFBS then also carries out a backward pass to generate an
approximate sample from the joint smoothing distribution
$\mathbb{Q}_{T}[dx_{1:T}]$. In brief, given current sample
$x_{t:T}\sim\mathbb{Q}_{T}[dx_{t:T}]$, $t\geq 2$, FFBS provides
$x_{t-1}|x_{t:T}\sim\mathbb{Q}_{T}[dx_{t-1}|x_{t:T}]$ via, first, use of the
formula:

$\displaystyle\mathbb{Q}_{T}[dx_{t-1}|x_{t:T}]$ $\displaystyle\propto
M_{t}(x_{t}|x_{t-1})\,G_{t}(x_{t-1},x_{t})\times\mathbb{Q}_{t-1}[dx_{t-1}].$

Then, the filtering distribution $\mathbb{Q}_{t-1}[dx_{t-1}]$ is replaced by
its Monte Carlo approximation provided during the forward pass, thus a proxy
for $\mathbb{Q}_{T}[dx_{t-1}|x_{t:T}]$ is obtained by re-weighting the
particles representing the filtering law $\mathbb{Q}_{t-1}[dx_{t-1}]$.

For concreteness, Algorithms 1, 2 provide pseudo-code for standard versions of
the particle filter and FFBS, respectively. We provide the following
explanations for the two pseudo-codes: $N\geq 1$ is the number of particles;
operations over $j$ are repeated independently for $j=1,\ldots,N$; operation
$\texttt{Resample}(W^{1:N})$, with $W^{j}\geq 0$ and $\sum_{j=1}^{N}W^{j}=1$,
applies multinomial (or some other type of) resampling; operation
$\mathcal{M}(W^{1:N})$ samples an index in $1,\ldots,N$ according to the
probabilities $W^{1:N}$.

Input: $M_{1}[dx_{1}]$, $M_{t}[x_{t-1},dx_{t}]$ and $G_{1}(x_{1})$,
$G_{t}(x_{t-1},x_{t})$, $t=2,\ldots,T$.

Output: Particles $X^{1:N}_{1:T}$, ancestors $A^{1:N}_{2:T}$, weights
$W_{1:T}^{1:N}$.

1 Sample: $X^{j}_{1}\sim M_{1}[dx_{1}]$;

2 Assign weights: $w^{j}_{1}\leftarrow G_{1}(X^{j}_{1})$;

3 Normalise weights:
$W^{j}_{1}\leftarrow\frac{w^{j}_{1}}{\sum_{k=1}^{N}w^{k}_{1}}$;

4 for  _$t=2,...,T$_ do

5 Sample: $A_{t}^{1:N}\sim\texttt{Resample}(W^{1:N}_{t-1})$;

6 Sample: $X_{t}^{j}\sim M_{t}[X_{t-1}^{A_{t}^{j}},dx_{t}]$;

7 Assign weights: $w_{t}^{j}\leftarrow G_{t}(X_{t-1}^{A_{t}^{j}},X_{t}^{j})$;

8 Normalise weights:
$W_{t}^{j}\leftarrow\frac{w_{t}^{j}}{\sum_{k=1}^{N}w_{t}^{k}}$;

9 end for

Algorithm 1 Particle Filter for Feynman-Kac model in (5).

Input: Full particle output, $X_{1:T}^{1:N}$, $A_{2:T}^{1:N}$,
$W_{1:T}^{1:N}$, from Algorithm 1.

Output: Approximate sample $(X_{1}^{B_{1}},\dots,X_{T}^{B_{T}})$ from
$\mathbb{Q}_{T}[dx_{1:T}]$.

1 Sample: $B_{T}\sim\mathcal{M}(W_{T}^{1:N})$;

2 for  _$t=T,...,2$_ do

3 Assign weights: $\hat{w}_{t-1}^{j}\leftarrow
M_{t}(X_{t}^{B_{t}}|X_{t-1}^{j})\,G_{t}(X_{t-1}^{j},X_{t}^{B_{t}})\cdot
W_{t-1}^{j}$;

4 Normalise weights:
$\hat{W}_{t-1}^{j}\leftarrow\frac{\hat{w}_{t-1}^{j}}{\sum_{k=1}^{N}\hat{w}_{t-1}^{k}}$;

5 Sample: $B_{t-1}\sim\mathcal{M}(\hat{W}_{t-1}^{1:N})$;

6 end for

Algorithm 2 FFBS for Feynman-Kac model in (5).

##  3 Continuous-Discrete SSM

###  3.1 Preliminaries

We return to the CD-SSM specified via (1)-(2). We introduce the path-valued
random variables $V_{t}=\\{V_{t}(s)\\}_{s\in[0,\Delta_{t}]}$, via use of the
time-shift $V_{t}(s)=X(s+s_{t-1})$. In addition, we consider the end point
$E_{t}=V_{t}(\Delta_{t})\equiv\texttt{end}(V_{t})$. We define the full latent
state at step $t\in\mathcal{T}$ as:

$\displaystyle X_{t}:=(V_{t},E_{t});\qquad
X_{t}\in\mathsf{X}_{t}:=C([0,\Delta_{t}],\mathbb{R}^{d})\times\mathbb{R}^{d}.$
(7)

We write realisations of $X_{t}$ as $x_{t}=(v_{t},e_{t})$.

The random variables $(X_{1:T},Y_{1:T})$ form a discrete-time SSM, with joint
distribution:

$P[X_{1:T}\in dx_{1:T},Y_{1:T}\in
dy_{1:T}]=\Big{[}\bigotimes_{t=1}^{T}f_{t}(y_{t}|e_{t})\,dy_{1:T}\Big{]}\otimes\mathbb{P}_{T}[dx_{1:T}],$
(8)

where $\mathbb{P}_{T}[dx_{1:T}]$ is the Markov measure on
$\mathsf{X}_{1}\times\cdots\times\mathsf{X}_{T}$ corresponding to the
distribution of time-shifted segments of the latent SDE model in (1) on the
time intervals $[0,\Delta_{t}]$, $t\in\mathcal{T}$. That is, we write:

$\displaystyle\mathbb{P}_{T}[dx_{1:T}]=P_{1}[dx_{1}]\bigotimes_{t=2}^{T}P_{t}[x_{t-1},dx_{t}],$

for accordingly defined kernels $P_{t}[x_{t-1},dx_{t}]$, $t\in\mathcal{T}$,
such that $P_{t}[x_{t-1},dx_{t}]\equiv P_{t}[e_{t-1},dx_{t}]$, since the
distribution of $X_{t}$ depends on $X_{t-1}$ only through the end-point
$E_{t-1}$. As done here, we will often replace argument $x_{t}=(v_{t},e_{t})$
with $e_{t}$ in probability expressions when these only depend on $e_{t}$
without further mention. From (8), we obtain the conditional law:

$P[X_{1:T}\in dx_{1:T}|Y_{1:T}\in
dy_{1:T}]=\frac{1}{f_{T}(y_{1:T})}\Big{[}\prod_{t=1}^{T}f_{t}(y_{t}|e_{t})\Big{]}\times\mathbb{P}_{T}[dx_{1:T}],$
(9)

with $f_{T}(\cdot)$ denoting the marginal pdf of $Y_{1:T}$.

A main contribution of the paper is to provide effective particle-based
algorithms for the approximation of the target smoothing measure (9), and of
related distributions arising when working with CD-SSMs (e.g. filtering
distributions, posterior laws of static parameters). In particular, we aim to
develop tailored Feynman-Kac models via careful transformations of the
variates involved in (9), so that important challenges arising due to the
continuous-time nature of the signal can be overcome – see challenges
(C.i)-(C.iii) stated in the Introduction. Once we appropriately define the
Feynman-Kac model, we will obtain access to the complete toolkit of state-of-
the-art particle-based methods for carrying out statistical inference for
SSMs.

###  3.2 Feynman-Kac Formulations

We begin by writing down the ‘forward decomposition’ for
$P_{t}[x_{t-1},dx_{t}]$, $t\in\mathcal{T}$:

$\displaystyle\textrm{F-Decomposition:}\qquad\qquad P_{t}[x_{t-1},dx_{t}]$
$\displaystyle=P_{t}[e_{t-1},dv_{t}]\otimes\delta_{\texttt{end}(v_{t})}(de_{t}).$
(10)

The kernel $P_{t}[e_{t-1},dv_{t}]$ corresponds to the law of the time-shifted
process $V_{t}=\\{X_{t}(s)\\}_{s\in[0,\Delta_{t}]}$. Note that $X_{t}(s)$ is
an diffusion process with drift function $b_{t}(s,x)=b(s_{t-1}+s,x)$ and
diffusion $\sigma_{t}(s,v)=\sigma(s_{t-1}+s,v)$. I.e., $V_{t}$ is obtained as
the solution of the following SDE on $[0,\Delta_{t}]$:

$\displaystyle P_{t}[e_{t-1},dv_{t}]:\qquad\quad
dX_{t}(s)=b_{t}(s,X_{t}(s))ds+\sigma_{t}(s,X_{t}(s))dB(s),\quad
X_{t}(0)=e_{t-1}.$ (11)

Expression (10) involves the distribution of the path $V_{t}$ given the
starting point $E_{t-1}$, along with that of the ending point $E_{t}$, the
latter being fully determined by $V_{t}$ in the construction shown so far.

An alternative ‘backward decomposition’ of $P_{t}[x_{t-1},dx_{t}]$ can be
constructed by instead considering first the distribution of the ending point
$E_{t}$ given the starting point $E_{t-1}$ and then the distribution of
$V_{t}|E_{t-1},E_{t}$, that latter variate being a _diffusion bridge_
conditioned to start at $E_{t-1}$ and end at $E_{t}$. In a probabilistic
language, this second construction provides the following expression for
$P_{t}[x_{t-1},dx_{t}]$, $t\in\mathcal{T}$:

$\displaystyle\textrm{B-Decomposition:}\qquad P_{t}[x_{t-1},dx_{t}]$
$\displaystyle=P_{t}[e_{t-1},de_{t}]\otimes P_{t}[(e_{t-1},e_{t}),dv_{t}]$
$\displaystyle\qquad\qquad=\big{(}p_{t}(e_{t}|e_{t-1})de_{t}\big{)}\otimes
P_{t}[(e_{t-1},e_{t}),dv_{t}].$ (12)

Here, $p_{t}(\cdot|e_{t-1})$ is the (typically intractable) pdf of
$X_{t}(\Delta_{t})|X_{t}(0)=e_{t-1}$ and $P_{t}[(e_{t-1},e_{t}),dv_{t}]$ is
the law of the bridge of the signal SDE (1) between times $s_{t-1}$ and
$s_{t}$, pinned at the starting point $e_{t-1}$ and the ending point $e_{t}$.
Under regularity conditions, such a bridge is known to correspond to the path-
valued process defined as the solution of the SDE:

$\displaystyle P_{t}[(e_{t-1},e_{t}),dv_{t}]:\qquad dX_{t}(s)$
$\displaystyle=b_{t}(s,X_{t}(s);e_{t})ds+\sigma_{t}(s,X_{t}(s))dB_{s},\quad
X_{t}(0)=e_{t-1},$ (13)

for the following drift function ($\nabla_{\\!x}$ denotes the vector of
partial derivatives w.r.t. argument $x$):

$\displaystyle b_{t}(s,x;e_{t})$
$\displaystyle=b_{t}(s,x)+\sigma_{t}(s,x)\sigma^{\top}_{t}(s,x)\nabla_{\\!x}\log(p_{\Delta_{t},e_{t}}(s,x)),$
(14)

where we have made use of the following function related with the SDE
transition density:

$\displaystyle p_{\Delta_{t},e_{t}}(s,x)=P[X_{t}(\Delta_{t})\in
de_{t}|X_{t}(s)=x]/\mathrm{Leb}(de_{t}),\qquad(s,x)\in[0,\Delta_{t}]\times\mathbb{R}^{d}.$
(15)

The bridge SDE determined via (13)-(14) can be obtained via use of Doob’s
$h$-transform, see e.g. Mider et al., (2021) for a recent exposition.

###  3.3 Bootstap Particle Filter

Given the above sequence $P_{t}[x_{t-1},dx_{t}]$, $t\in\mathcal{T}$, the
targeted conditional law (9) of the CD-SSM has been written in the form of the
Feynman-Kac formula (5), with Markov kernels $M_{t}\equiv P_{t}$ and weights
given by $G_{t}(x_{t-1},x_{t})\equiv f_{t}(y_{t}|e_{t})$. For such a
specification, the probability law $\mathbb{Q}_{t}[dx_{1:t}]$ in (5) coincides
with the conditional distribution of $X_{1:t}|Y_{1:t}$, as given in (9) for
$t=T$.

We stress that, as discussed in the Introduction, the above construction can
provide a non-degenerate particle filter that targets our CD-SSM via execution
of Algorithm 1. However, attempts to use this construct within FFBS in
Algorithm 2, and related smoothing-based algorithms, will lead to a degenerate
methodology. Even in the context of a particle filter, greater flexibility is
desirable and can be made possible with the selection of ‘guided’ proposal
kernels. We present such a guided formalism in the next section, where we
target the conditional distribution $X_{1:t}|Y_{1:t}$, by making use of
proposal kernels that are non-blind to the data.

##  4 Guided Feynman-Kac Formulae for CD-SSMs

For a ‘guided’ Feynman-Kac formalism of a SSM (see e.g. Chopin and
Papaspiliopoulos, (2020)) we can assume the selection of a Markov measure
$\mathbb{M}_{T}[dx_{1:T}]$ in our Feynman-Kac formulation that does not
coincide with the prior dynamics of the latent signal,
$\mathbb{P}_{T}[dx_{1:T}]$. That is, the Feynman-Kac expression will now
involve the following weights:

$\displaystyle
G_{t}(x_{t-1},x_{t})=\frac{dP_{t}[x_{t-1},\cdot]}{dM_{t}[x_{t-1},\cdot]}(x_{t})\times
f_{t}(y_{t}|x(s_{t})).$ (16)

Given that absolute continuity $P_{t}[x_{t-1},dx_{t}]\ll
M_{t}[x_{t-1},dx_{t}]$ holds, the weights in (16) will be well-defined.

We outline two different approaches for the selection of guided proposals that
take the data point $y_{t}$ into account. We will refer to our first construct
as the ‘forward proposal’, as it can be thought as building upon the forward
decomposition (10) of $P_{t}[x_{t-1},dx_{t}]$. We call the second the
‘backward proposal’, as it is based on the corresponding backward
decomposition (3.2). Though seemingly similar, it will turn out that these two
constructs have subtle differences. In particular, the backward proposal will
allow for the derivation of a Feynman-Kac formula (and access to all
accompanying particle-based algorithms) in the case of hypo-elliptic SDEs,
when the diffusion coefficient is degenerate. That is not possible for the
forward construct. The forward construct resembles several standard approaches
followed in the literature in the context of MCMC algorithms and particle
filters. The backward construct has not appeared before in the literature, to
the best of our knowledge.

######  Remark 1.

In the context of the forward proposal, and the forward method (FM) that we
will later on obtain via use of such a proposal, we will henceforth be
assuming that $d_{w}\equiv d$, and that $\sigma_{t}(s,x)\in\mathbb{R}^{d\times
d}$ is invertible for all $(s,x)\in[0,\Delta_{t}]\times\mathbb{R}^{d}$. That
is, FM will only be applicable to elliptic SDEs. In contrast, we allow full
generality, with $d_{w}\leq d$ in the setting of the backward proposal, and
the backward method (BM), so that BM will also cover the important class of
hypo-elliptic SDEs.

###  4.1 Forward Proposal

####  4.1.1 Construction of Forward Proposal

We define the forward proposal $M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ as
follows:

$\displaystyle\textrm{F-Proposal:}\qquad\qquad\quad
M_{t}^{\rightarrow}[x_{t-1},dx_{t}]=M_{t}^{\rightarrow}[e_{t-1},dv_{t}]\otimes\delta_{\texttt{end}(v_{t})}(de_{t}),$
(17)

where we assume that $P_{t}[e_{t-1},dv_{t}]\ll
M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$. The proposal kernel
$M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$ corresponds to the solution of an SDE
with the same diffusion coefficient as the target $P_{t}[e_{t-1},dv_{t}]$ in
(11), and a drift function given as
$b^{\rightarrow}_{t}(\cdot,\cdot;y_{t}):[0,\Delta_{t}]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$.
In particular, the drift function of the proposal SDE aims to take into
account information from the observation $y_{t}$. We thus define the proposal
SDE and its corresponding law as follows:

$\displaystyle M_{t}^{\rightarrow}[e_{t-1},dv_{t}]:\qquad\qquad dX_{t}(s)$
$\displaystyle=b^{\rightarrow}_{t}(s,X_{t}(s);y_{t})ds+\sigma_{t}(s,X_{t}(s))dB(s),$
(18) $\displaystyle X_{t}(0)$ $\displaystyle=e_{t-1}.$

We set $\Sigma_{t}=\sigma_{t}\sigma_{t}^{\top}\in\mathbb{R}^{d\times d}$.
Recall that we assume that $\Sigma_{t}=\Sigma_{t}(s,x)$ is invertible for all
$(s,x)\in[0,\Delta_{t}]\times\mathbb{R}^{d}$. Girsanov’s theorem (Oksendal,,
2013) states that, under regularity conditions, we have the following Radon-
Nikodym derivative:

$\displaystyle\frac{dP_{t}[e_{t-1},\cdot]}{dM_{t}^{\rightarrow}[e_{t-1},\cdot]}(v_{t})$
$\displaystyle\qquad\qquad=\exp\Bigg{[}\int_{0}^{\Delta_{t}}\big{\\{}(b_{t}-b^{\rightarrow}_{t})^{\top}\Sigma_{t}^{-1}\big{\\}}(s,v_{t}(s);y_{t})dv_{t}(s)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\tfrac{1}{2}\int_{0}^{\Delta_{t}}\big{\\{}(b_{t}-b^{\rightarrow}_{t})^{\top}\Sigma_{t}^{-1}(b_{t}+b^{\rightarrow}_{t})\big{\\}}(s,v_{t}(s);y_{t})ds\Bigg{]}.$
(19)

In brief, sufficient regularity conditions involve global Lipschitz continuity
and linear or slower growth for the drifts and diffusion coefficient, together
with Novikov’s condition (Oksendal,, 2013). The likelihood (19) can be
evaluated up-to time-discretisation.

Upon choosing a tractable drift function $b^{\rightarrow}_{t}(s,x;y_{t})$, one
can sample from $M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$. The corresponding
potential $G_{t}^{\rightarrow}(x_{t-1},x_{t})$ of the generated Feynmac-Kac
model writes as:

$\displaystyle\textrm{F-Weight:}\qquad\qquad\qquad
G_{t}^{\rightarrow}(x_{t-1},x_{t})=\frac{dP_{t}[x_{t-1},\cdot]}{dM_{t}^{\rightarrow}[x_{t-1},\cdot]}(x_{t})\times
f_{t}(y_{t}|x(s_{t})),$ (20)

where the involved Radon-Nikodym derivative is the one given in (19).

######  Remark 2.

The proposal $M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ is defined via the forward
decomposition (17), in agreement with a similar expression for the signal
kernel $P_{t}[x_{t-1},dx_{t}]$ in (10). Notice though that
$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ also admits a backward-type
decomposition, as one can write:

$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]=M_{t}^{\rightarrow}[e_{t-1},dx_{t}]=\big{(}m_{t}^{\rightarrow}(e_{t}|e_{t-1})de_{t}\big{)}\otimes
M_{t}^{\rightarrow}[(e_{t-1},e_{t}),dv_{t}].$ (21)

In (21), $m_{t}^{\rightarrow}(\cdot|\cdot)$ is the, typically intractable,
transition density from time $0$ to $\Delta_{t}$ of the proposal SDE (18), and
$M_{t}^{\rightarrow}[(e_{t-1},e_{t}),dv_{t}]$ is the distribution of the
bridge of SDE (18). Expression (21) will be relevant when presenting the
complete Forward Method (FM) in the sequel.

####  4.1.2 Choice of Forward Proposals

We give some guidance over the choice of the drift function
$b^{\rightarrow}_{t}(s,x;y_{t})$ in the proposal SDE (18) that takes the
observation $y_{t}$ into account. Sottinen and Särkkä, (2008) develop proposal
SDEs for the case of state-independent diffusion coefficient by using
expressions from a continuous-discrete extended Kalman Filter to form the
drift function. Alternative proposal SDEs can be obtained by adapting ideas
for diffusion bridges (in settings of observations without noise) to the case
of noisy observations (Delyon and Hu,, 2006; Llopis et al.,, 2018).

We present in more detail approaches based on the ‘optimal’ proposal SDE and
its approximation, following the recent developments in Schauer et al.,
(2017); Mider et al., (2021) – the relevant methodology provides a flexible
framework that provides proposal SDEs for a broad class of diffusion models.
The optimal proposal minimises the variance of the incremental weights (Doucet
et al.,, 2000) and is obtained by conditioning the signal dynamics on the
observation $y_{t}$. That is:

$\displaystyle M_{t}^{\rightarrow,\textrm{opt}}[e_{t-1},dv_{t}]\propto
P_{t}[e_{t-1},dv_{t}]\times f_{t}(y_{t}|e_{t}).$

Under conditions, the distribution
$M_{t}^{\rightarrow,\textrm{opt}}[e_{t-1},dv_{t}]$ coincides with the law of
an SDE with the following drift function (see Mider et al., (2021) for a
proof):

$b^{\rightarrow,\textrm{opt}}_{t}(s,x;y_{t})=b_{t}(s,x)+\Sigma_{t}(s,x)\nabla_{x}(\log\rho_{\Delta_{t},y_{t}}(s,x)),$
(22)

where we make use of the density function:

$\rho_{\Delta_{t},y_{t}}(s,x)=\int_{\mathbb{R}^{d}}f_{t}(y_{t}|x^{\prime})p_{\Delta_{t},x^{\prime}}(s,x)dx^{\prime},\qquad(s,x)\in[0,\Delta_{t}]\times\mathbb{R}^{d},$
(23)

for $p_{\Delta_{t},x^{\prime}}(s,x)$ as defined in (15).

The expression for $\rho_{t,y}$ is typically intractable, as it involves an
integral w.r.t. the transition density of the signal SDE. Following the
approach of Mider et al., (2021), a reasonable choice can be obtained by
replacing $\rho_{\Delta_{t},y_{t}}(s,x)$ in (22) with an appropriate proxy
$\tilde{\rho}_{\Delta_{t},y_{t}}(s,x)$. Thus, one makes use of the following
drift function:

$b^{\rightarrow}_{t}(s,x;y_{t})=b_{t}(s,x)+\Sigma_{t}(s,x)\nabla_{x}(\log\tilde{\rho}_{t,y_{t}}(s,x)).$
(24)

Function $\tilde{\rho}_{\Delta_{t},y_{t}}$ can be obtained by replacing
$p_{\Delta_{t},x^{\prime}}(s,x)$ in the integral (23) with a corresponding
proxy $\tilde{p}_{\Delta_{t},x^{\prime}}(s,x)$. If
$x^{\prime}\mapsto\tilde{p}_{\Delta_{t},x^{\prime}}(s,x)$ is Gaussian and the
observation operator for $y_{t}$ is linear (the latter being a very common
case in applications), then the integral (23) is tractable. In Mider et al.,
(2021), $\tilde{p}_{\Delta_{t},x^{\prime}}(s,x)$ is chosen via use of the
transition density of a member within a class of linear SDEs.

######  Remark 3.

Consider the following linear SDE:

$d\tilde{X}_{t}(s)=\big{\\{}\tilde{b}_{0}(s)+\tilde{b}_{1}(s)\tilde{X}_{t}(s)\big{\\}}ds+\tilde{\sigma}(s)dB(s),$
(25)

for $\tilde{b}_{0}:[0,\Delta_{t}]\rightarrow\mathbb{R}^{d}$,
$\tilde{b}_{1}:[0,\Delta_{t}]\rightarrow\mathbb{R}^{d\times d}$ and
$\tilde{\sigma}:[0,\Delta_{t}]\rightarrow\mathbb{R}^{{d\times d_{w}}}$, with
$d_{w}\leq d$. Such a linear SDE has an analytical solution that writes as
follows, for $s\in[0,\Delta_{t}]$ and $\tilde{X}_{t}(s)=x\in\mathbb{R}^{d}$:

$\displaystyle\tilde{X}_{t}(\Delta_{t})=\Phi(\Delta_{t},s)\,x+\int_{s}^{\Delta_{t}}\Phi(\Delta_{t},\tau)\tilde{b}_{0}(\tau)\,d\tau+\int_{s}^{\Delta_{t}}\Phi(\Delta_{t},\tau)\tilde{\sigma}(\tau)\,dB(\tau),$

where, for $s\in[0,\Delta_{t}]$, we set $\Phi(s,u)=\Phi(s)\Phi(u)^{-1}$, with
$\Phi(\cdot)$ obtained from the solution of the ODE specified as
$d\Phi(t)=\tilde{b}_{1}(t)\Phi(t)dt$, $\Phi(0)=\mathrm{I}_{d}$. Thus, the
conditional law of $\tilde{X}_{t}(\Delta_{t})|\tilde{X}_{t}(s)=x$ corresponds
to a Gaussian distribution $\mathcal{N}(\mu(x),\mathrm{V})$, with:

$\displaystyle\mu(x)=\Phi(\Delta_{t},s)x+\int_{s}^{\Delta_{t}}\Phi(\Delta_{t},\tau)\tilde{b}_{0}(\tau)\,d\tau;\qquad\mathrm{V}=\int_{s}^{\Delta_{t}}\Phi(\Delta_{t},\tau)(\tilde{\sigma}\tilde{\sigma}^{\top})(\tau)\Phi(\Delta_{t},\tau)^{\top}\,d\tau.$

One aims to choose $\tilde{b}_{0}(s)$, $\tilde{b}_{1}(s)$, and
$\tilde{\sigma}(s)$ so that the linear SDE (25) approximates the dynamics
$P_{t}[e_{t-1},dv_{t}]$ of the signal SDE. A strategy for specifying (25) can
involve linearisation of the drift of the signal SDE (11), or the following
simpler approaches (Mider et al.,, 2021):

$\displaystyle\tilde{b}_{0}(s)=b_{t}(0,e_{t-1}),\,\,\tilde{b}_{1}(s)=0,\,\,\tilde{\sigma}(s)=\sigma_{t}(0,e_{t-1});\qquad\tilde{b}_{0}(s)=0,\,\,\tilde{b}_{1}(s)=0,\,\,\tilde{\sigma}(s)=\mathrm{I}_{d}.$

###  4.2 Backward Proposal

####  4.2.1 Construction of Backward Proposal

Recall that we now allow for $d_{w}\leq d$, as we cover both the elliptic and
hypo-elliptic case for the signal SDE in (1).

In this direction, we define the backward proposal
$M^{\leftarrow}_{t}[x_{t-1},dx_{t}]$ as follows:

$\displaystyle\textrm{B-Proposal:}\qquad\quad
M_{t}^{\leftarrow}[x_{t-1},dx_{t}]=\big{(}m_{t}^{\leftarrow}(e_{t}|e_{t-1})de_{t}\big{)}\otimes
M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}].$ (26)

The density $m_{t}^{\leftarrow}(e_{t}|e_{t-1})$ is user-specified, under the
requirements that one can evaluate it and sample from it. Then,
$M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]$ is a proposal for the target
bridge law $P_{t}[(e_{t-1},e_{t}),dv_{t}]$ in (13) connecting $e_{t-1}$ and
$e_{t}$, chosen so that $P_{t}[(e_{t-1},e_{t}),dv_{t}]\ll
M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]$. In such a case, the potential
$G_{t}^{\leftarrow}(x_{t-1},x_{t})$ for the induced Feynman-Kac model is
expressed in the form of the following product:

$\displaystyle\textrm{B-Weight:}\qquad G_{t}^{\leftarrow}(x_{t-1},x_{t})$
$\displaystyle=\frac{dP_{t}[e_{t-1},\cdot]}{dM_{t}^{\leftarrow}[e_{t-1},\cdot]}(v_{t})$
$\displaystyle\qquad\qquad=\frac{p_{t}(e_{t}|e_{t-1})}{m_{t}^{\leftarrow}(e_{t}|e_{t-1})}\times\frac{dP_{t}[(e_{t-1},e_{t}),\cdot]}{dM^{\leftarrow}_{t}[(e_{t-1},e_{t}),\cdot]}(v_{t})\times
f_{t}(y_{t}|e_{t}).$ (27)

####  4.2.2 Choice of Backward Proposals

We provide some details on the specification of the backward proposal so that
$y_{t}$ is taken into account.

Choice of $m_{t}^{\leftarrow}(e_{t}|e_{t-1})$ – Elliptic and Hypo-Elliptic
Case:

The optimal proposal (Doucet et al.,, 2000), as also discussed in the forward
setting of Section 4.1, denoted here as $m_{t}^{\leftarrow,\textrm{opt}}$, is
obtained from the true conditional distribution of $(e_{t},v_{t})$ given
$e_{t-1}$ and $y_{t}$. Thus, in the current backward setting, we can write
that:

$M_{t}^{\leftarrow,\textrm{opt}}[e_{t-1},dx_{t}]=\big{(}m_{t}^{\leftarrow,\textrm{opt}}(e_{t}|e_{t-1})de_{t}\big{)}\otimes
P_{t}[(e_{t-1},e_{t}),dv_{t}],$ (28)

where the optimal proposal for the end point $e_{t}$, is given by:

$m_{t}^{\leftarrow,\textrm{opt}}(e_{t}|e_{t-1})\propto
f_{t}(y_{t}|e_{t})p_{t}(e_{t}|e_{t-1}).$ (29)

The expression for $m_{t}^{\leftarrow,\textrm{opt}}$ involves the typically
intractable transition density of the SDE model for the latent signal. From
this point onwards, one can follow an approach similar to the one set out
within Section 4.1, where the true transition density needed within the drift
of the optimal SDE proposal is replaced by a tractable proxy. Thus, in the
current setting, to obtain an alternative to
$m_{t}^{\leftarrow,\textrm{opt}}(e_{t}|e_{t-1})$ in (29), one can replace
$p_{t}(e_{t}|e_{t-1})$ with a proxy, $\tilde{p}_{t}(e_{t}|e_{t-1})$, so that:

$m_{t}^{\leftarrow}(e_{t}|e_{t-1})\propto
f_{t}(y_{t}|e_{t})\tilde{p}_{t}(e_{t}|e_{t-1}).$ (30)

In the scenario when $e_{t}\mapsto\tilde{p}_{t}(e_{t}|e_{t-1})$ in a Gaussian
pdf and the observation operator providing $y_{t}$ is linear in $e_{t}$, then
one obtains a Gaussian density for $m^{\leftarrow}_{t}(e_{t}|e_{t-1})$.

Proceeding as in the discussion for the forward proposal in Section 4.1, one
can construct proxies $\tilde{p}_{t}(e_{t}|e_{t-1})$ that are Gaussian
conditionally on $e_{t-1}$, by utilising the transition density of the linear
SDE (25). We consider here in detail SDE (25) as specified by a linearisation
of the drift $b_{t}(s,x;y_{t})$ of the signal SDE (11), as its transition
density will be non-degenerate even for hypo-elliptic SDEs. That is, applying
a first-order Taylor expansion around $e_{t-1}$ upon the $x$-component of the
drift function $b_{t}(s,x)$, one obtains the following linear proxy of the
true drift:

$\displaystyle b_{t}(0,e_{t-1})+(D_{x}b_{t})(0,e_{t-1})(x-e_{t-1}),$

where $(D_{x}b_{t})(s,x^{\prime})\in\mathbb{R}^{d\times d}$ with
$(D_{x}b_{t})(s,x^{\prime})_{i,j}=(\partial_{x_{j}}b_{t}^{i})(s,x^{\prime})$,
$1\leq i,j\leq d$. Thus, in the notation of (25), we have:

$\displaystyle\tilde{b}_{0}(s)=b_{t}(0,e_{t-1})-(D_{x}b_{t})(0,e_{t-1})e_{t-1};\quad\tilde{b}_{1}(s)=(D_{x}b_{t})(0,e_{t-1})x;\quad\tilde{\sigma}(s)=\sigma_{t}(0,e_{t-1}).$

The transition density $\tilde{p}_{t}(e_{t}|e_{t-1})$ can now be obtained via
the formulae in Remark 3. In combination with a linear observation operator,
standard calculations will yield a Gaussian density for
$m_{t}^{\leftarrow}(e_{t}|e_{t-1})$ in (30). We note that simpler choices for
$\tilde{p}_{t}(e_{t}|e_{t-1})$ can also be considered. These can range from a
standard Euler scheme in the elliptic case, to schemes tailored to particular
classes of hypo-elliptic SDEs, see for example Pokern et al., (2009);
Ditlevsen and Samson, (2019).

Choice of $M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]$:

We now proceed to the choice of the proposal
$M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$ for the bridge. Recall that the
true diffusion bridge is given by (13)-(14), with a drift that involves the
intractable transition density of the signal SDE. A general expression for a
proposal bridge SDE can be follows:

$\displaystyle M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]:\qquad\quad
dX_{t}(s)$
$\displaystyle=b^{\leftarrow,br}_{t}(s,X_{t}(s);e_{t})ds+\sigma_{t}(s,X_{t}(s))dB(s),$
(31) $\displaystyle X_{t}(0)$ $\displaystyle=e_{t-1}.$

The drift $b^{\leftarrow,br}_{t}(s,x;e_{t})$ appearing above must have a
particular structure, as it must force paths to terminate at point $e_{t}$ at
time $\Delta_{t}$, while preserving absolute continuity w.r.t. the law of the
bridge of the latent SDE (1) in the CD-SSM. The works in Delyon and Hu,
(2006); Schauer et al., (2017); Bierkens et al., (2020) provide detailed
analytical studies into the theme of diffusion bridges. The latter of these
(Bierkens et al.,, 2020) covers both elliptic and hypoelliptic SDEs. We
briefly summarise some standard choices for SDE bridges that one can sample
from, and for which analytical expressions of Radon-Nikodym derivatives w.r.t.
target SDE bridges have been obtained. We stress that further options can be
directly integrated within the framework of our methodology.

Choice of $M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$ – Elliptic Case, Delyon
and Hu, (2006):

In the context of elliptic SDEs, with $d_{w}=d$ and
$\sigma_{t}(s,x)\in\mathbb{R}^{d\times d}$ assumed positive-definite, a simple
choice used in the literature is the Delyon-Hu bridge (Delyon and Hu,, 2006),
when:

$\displaystyle b^{\leftarrow,br}_{t}(s,x;e_{t})=\frac{e_{t}-x}{\Delta_{t}-s},$

in which case the bridge SDE (31) writes as follows:

$\displaystyle M_{t}^{\leftarrow,\textrm{DH}}[(e_{t-1},e_{t}),dv_{t}]:\qquad
dX_{t}(s)=\frac{e_{t}-X_{t}(s)}{\Delta_{t}-s}ds+\sigma_{t}(s,X_{t}(s))dB(s),\quad
X_{t}(0)=e_{t-1}.$ (32)

Under regularity conditions, Delyon and Hu, (2006) obtain the following Radon-
Nikodym derivative:

$\displaystyle\frac{dP_{t}[(e_{t-1},e_{t}),\cdot]}{dM^{\leftarrow,\textrm{DH}}_{t}[(e_{t-1},e_{t}),\cdot]}(v_{t})$
$\displaystyle\qquad\qquad\qquad=\frac{\mathcal{N}(e_{t};e_{t-1},\Sigma_{t}(0,e_{t-1}))}{p_{t}(e_{t}|e_{t-1})}\sqrt{\frac{\det(\Sigma_{t}(0,e_{t-1}))}{\det(\Sigma_{t}(\Delta_{t},e_{t}))}}\times\exp\big{\\{}\phi_{t}(v_{t})\big{\\}},$
(33)

where we have defined the following functional:

$\displaystyle\phi_{t}(v_{t})$
$\displaystyle=\int_{0}^{\Delta_{t}}(b_{t}^{\top}\Sigma_{t}^{-1})(s,v_{t}(s))dv_{t}(s)$
$\displaystyle\qquad-\tfrac{1}{2}\int_{0}^{\Delta_{t}}(b_{t}^{\top}\Sigma_{t}^{-1}b_{t})(s,v_{t}(s))ds$
$\displaystyle\qquad\qquad-\tfrac{1}{2}\int_{0}^{\Delta_{t}}\frac{e_{t}-v_{t}(s)}{\Delta_{t}-s}\diamond
d\big{\\{}\Sigma_{t}^{-1}(s,v_{t}(s))(e_{t}-v_{t}(s))\big{\\}}.$

In brief, the integral with ‘$\diamond$’ is defined by taking sums on sub-
intervals that partition $[0,\Delta_{t}]$ and replacing the integrand in each
summand with its right-side value, before taking the limit for vanishing
length of the sub-intervals (when Itô integrals use left-side values). We note
that the normalising constant in the above Radon-Nikodym derivative is omitted
in Delyon and Hu, (2006), but can be obtained via use of Bayes’ theorem, see
e.g. (van der Meulen and Schauer,, 2017, Section 4.1). Use of (4.2.2) within
(4.2.1) provides an analytical expression of the potential
$G_{t}^{\leftarrow}(x_{t-1},x_{t})$ for the induced Feynman-Kac model. The
intractable transition density $p_{t}(e_{t}|e_{t-1})$ appearing in (4.2.1)
will cancel out when obtaining $G^{\leftarrow}_{t}(x_{t-1},x_{t})$ as it also
appears as a denominator in the density expression for the diffusion bridges
in (4.2.2). All remaining terms in the formula for $G_{t}(x_{t-1},x_{t})$ can
be analytically evaluated, up-to time-discretisation. Note that the above
choice of proposal bridge does not take into account the drift of the signal
SDE.

Choice of $M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$ – Elliptic Case,
Schauer et al., (2017):

Schauer et al., (2017) develop a class of ‘guided proposals’ for diffusion
bridges, in the setting of elliptic SDEs, that are of the following form:

$\displaystyle
M_{t}^{\leftarrow,\textrm{G}}[(e_{t-1},e_{t}),dv_{t}]:\qquad\qquad dX_{t}(s)$
$\displaystyle=b_{t}^{\mathrm{G}}(s,X_{t}(s);e_{t})ds+\sigma_{t}(s,X_{t}(s))dB(s),$
(34) $\displaystyle X_{t}(0)$ $\displaystyle=e_{t-1},$

where we have defined:

$\displaystyle b_{t}^{\mathrm{G}}(s,x;e_{t})$
$\displaystyle=b_{t}(s,x)+\Sigma_{t}(s,x)\nabla_{x}\log(\tilde{p}_{\Delta_{t},e_{t}}(s,x)),$
(35)

for a transition density $\tilde{p}_{\Delta_{t},e_{t}}(s,x)$, used as a proxy
to the true $p_{\Delta_{t},e_{t}}(s,x)$. In Schauer et al., (2017),
$\tilde{p}_{\Delta_{t},e_{t}}(s,x)$ is chosen to correspond to the transitions
of a linear SDE as given in (25), with constituent functions
$\tilde{b}_{0}(s)$, $\tilde{b}_{1}(s)$ and $\tilde{\sigma}(s)$, and
corresponding mean $\mu(x)$, variance $V$ and matrix $\Phi(\Delta_{s},s)$,
involved in expression of the Gaussian transition density of the linear SDE,
as specified in Remark 3. We define:

$\displaystyle\tilde{\Sigma}(s)$
$\displaystyle=\tilde{\sigma}(s)\tilde{\sigma}(s)^{\top};$
$\displaystyle\tilde{r}(s,x)$
$\displaystyle=\nabla_{x}\log\tilde{p}_{\Delta_{t},e_{t}}(s,x)\equiv\Phi(\Delta_{t},s)^{\top}V^{-1}(e_{t}-\mu(x));$
$\displaystyle\tilde{H}(s)$
$\displaystyle=-(\nabla_{x}\nabla_{x}^{\top})\log\tilde{p}_{\Delta_{t},e_{t}}(s,x)\equiv\Phi(\Delta_{t},s)^{\top}V^{-1}\Phi(\Delta_{t},s).$

Schauer et al., (2017) show that, under the main ‘matching’ requirement that:

$\displaystyle\tilde{\Sigma}(\Delta_{t})\equiv\Sigma_{t}(\Delta_{t},e_{t}),$
(36)

and given that additional regularity conditions are satisfied, then the law
for the class of SDEs defined in (34)-(35), with
$\tilde{p}_{\Delta_{t},e_{t}}(s,x)$ obtained via the linear SDE (25),
dominates the law of the diffusion bridge $P_{t}[(e_{t-1},e_{t}),dv_{t}]$,
with a Radon-Nikodym derivative:

$\frac{dP_{t}[(e_{t-1},e_{t}),\cdot]}{dM_{t}^{\leftarrow,\mathrm{G}}[(e_{t-1},e_{t}),\cdot]}(v_{t})=\frac{p^{\mathrm{G}}_{t}(e_{t}|e_{t-1})}{p_{t}(e_{t}|e_{t-1})}\exp\Big{\\{}\int_{0}^{\Delta_{t}}\phi(s,v_{t}(s))ds\Big{\\}},$
(37)

with
$p^{\mathrm{G}}_{t}(de_{t}|e_{t-1})=\mathcal{N}(e_{t-1},\Delta_{t}\mathrm{I}_{d})$,
where we have defined the function:

$\displaystyle\phi(s,x)=$
$\displaystyle(b_{t}(s,x)-\tilde{b}(s,x))^{\top}\tilde{r}(s,x)$
$\displaystyle\qquad-\tfrac{1}{2}\operatorname{trace}\Big{(}\big{(}\Sigma_{t}(s,x)-\tilde{\Sigma}(s)\big{)}\big{(}\tilde{H}(s)-\tilde{r}(s,x)\tilde{r}(s,x)^{\top}\big{)}\Big{)}.$

Specific choices of $\tilde{b}_{0}(s)$, $\tilde{b}_{1}(s)$, and
$\tilde{\sigma}(s)$ are discussed in van der Meulen and Schauer, (2017). E.g.,
a simple linear SDE satisfying the matching condition (36) is deduced when:

$\displaystyle
d\tilde{X}_{t}(s)=\sigma_{t}(\Delta_{t},e_{t})dB(s),\quad\tilde{X}_{t}(0)=e_{t-1},$

which yields the transition density
$p^{\mathrm{G}}_{t}(e_{t}|e_{t-1})=\mathcal{N}(e_{t};e_{t-1},\Delta_{t}\Sigma_{t}(\Delta_{t},e_{t}))$
to be used in (37). Note, that the drift of the guided proposal in this case
writes as:

$\displaystyle b_{t}^{\mathrm{G}}(s,x;e_{t})$
$\displaystyle=b_{t}(s,x)+\Sigma_{t}(s,x)\Sigma_{t}(\Delta_{t},e_{t})^{-1}\frac{e_{t}-x}{\Delta_{t}-s}.$

Choice of $M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$ – Hypo-Elliptic Case,
Bierkens et al., (2020):

The methodology developed in Bierkens et al., (2020) follows along the lines
of the elliptic case. That is, a guided proposal is constructed as in
(34)-(35), with a proxy $\tilde{p}_{\Delta_{t},e_{t}}(s,x)$ of the true
transition density $p_{\Delta_{t},e_{t}}(s,x)$ used within the drift function.
The tractable density $\tilde{p}_{\Delta_{t},e_{t}}(s,x)$ derives from the
dynamics of the linear SDE (25), however in the hypo-elliptic setting extra
care is required on the specification of the drift
$\tilde{b}_{0}(s)+\tilde{b}_{1}(s)x$ to yield absolute continuity between the
constructed proposal and the bridge of the signal SDE. ‘Matching’ restrictions
must now be placed both upon $\tilde{\sigma}(s)$ and the linear drift
$\tilde{b}_{0}(s)+\tilde{b}_{1}(s)x$. Under conditions (some of which we
briefly discuss below) one has absolute continuity between target and proposal
bridge, with a density given as in (37), with terms as defined therein. Note,
the Radon-Nikodym derivative (37) is well-defined for coefficient matrices
$\sigma_{t}(s,x),\tilde{\sigma}(s)\in\mathbb{R}^{d\times d_{w}}$, $d_{w}\leq
d$, even if these are non-rectangular.

##  5 Reparameterised Feynman-Kac Formulae

Assume now that we have chosen a kernel $M_{t}[x_{t-1},dx_{t}]$ within the
guided Feynman-Kac formula framework, either by using a forward proposal
$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ or a backward one
$M_{t}^{\leftarrow}[x_{t-1},dx_{t}]$ as in (17) or (26), respectively. The
corresponding potentials are then given by
$G_{t}^{\rightarrow}(x_{t-1},x_{t})$ in (20) or
$G_{t}^{\leftarrow}(x_{t-1},x_{t})$ in (4.2.1), so that the induced Feynman-
Kac model delivers the conditional distributions of the signal
$X_{1:t}|Y_{1:t}$, $t\in\mathcal{T}$, for the CD-SSM of interest. One can
immediately apply particle-based methods to approximate the filtering
distributions, using e.g. the standard version of a particle filter as
described in Algorithm 1. In contrast, and as already explained in the
Introduction, due to the continuous-time nature of the signal, the derivations
so far provided will lead to degenerate particle-based smoothing and parameter
inference methods, as will be the case for FFBS described in Algorithm 2. We
repeat here that in a probabilistic language this matter relates directly with
the fact that the proposal kernel $M_{t}[x_{t-1},dx_{t}]$, as defined so far,
does not have a density w.r.t. a common reference measure that will not depend
on the starting point $e_{t-1}$. The critical practical consequence of such a
property is that the ancestor $x_{t-1}$ of a particle $x_{t}$ cannot be
altered within any of the many algorithms that involve a smoothing
approximation procedure.

To overcome the above issue, we now construct new Feynman-Kac formulae, via a
transform of the models we have so far defined. That is, we carefully set-up
new proposal kernels $\bar{M}_{t}[z_{t-1},dz_{t}]$ and potentials
$\bar{G}_{t}(z_{t-1},z_{t})$, on corresponding measurable spaces
$(\mathsf{Z}_{t},\mathscr{Z}_{t})$, so that we overcome challenges
(C.i)-(C.iii) for CD-SSMs stated in the Introduction. Such constructions will
lead to a new sequence of target distributions
$\bar{\mathbb{Q}}_{t}[dz_{1:t}]$ on
$\big{(}\mathsf{Z}_{1}\times\cdots\times\mathsf{Z}_{t},\mathscr{Z}_{1}\otimes\cdots\otimes\mathscr{Z}_{t}\big{)}$,
through the standard Feynmac-Kac expression:

$\bar{\mathbb{Q}}_{t}[dz_{1:t}]=\frac{1}{\bar{L}_{t}}\Big{\\{}\bar{G}_{1}(z_{1})\prod_{i=2}^{t}\bar{G}_{i}(z_{i-1},z_{i})\Big{\\}}\times\bar{\mathbb{M}}_{t}[dz_{1:t}],\quad
t\in\mathcal{T}.$ (38)

We will sometimes refer to the Feynman-Kac model given by (38) as the
‘transformed Feynman-Kac formula/model’ of a CD-SSM. The key elements of this
construction (yet to be defined) are summarised as follows:

  * i.

We can sample from the kernels $\bar{M}_{t}[z_{t-1},dz_{t}]$, up-to time-
discretisation.

  * ii.

The kernels $\bar{M}_{t}[z_{t-1},dz_{t}]$ have a tractable density w.r.t. a
common reference measure that does not depend on the ancestor $z_{t-1}$. The
densities can be evaluated, up-to time-discretisation.

  * iii.

The potentials $\bar{G}_{t}(z_{t-1},z_{t})$ can be evaluated, up-to time-
discretisation.

  * iv.

The transform is obtained via maps ${\mathsf{Z}}_{t}\mapsto{\mathsf{X}}_{t}$,
$t\in\mathcal{T}$. The image of $\bar{\mathbb{Q}}_{t}[dz_{1:t}]$ under such
maps gives back $\mathbb{Q}_{t}[dx_{1:t}]$ from the original Feynman-Kac
model.

Given such a transformed Feynman-Kac model, one will be able to implement the
full spectrum of particle-based inference methods for the class of CD-SSMs.

We now provide the details for the construction of the transformed model (38).
Following the development of the forward and backward proposals in Section 4,
we will also present two different definitions for
$\bar{\mathbb{Q}}_{t}[dz_{1:t}]$. We refer to the first as the ‘forward
transform’, as it is applicable to the case when the original Feynman-Kac
model builds upon forward proposals, as described in Section 4.1. This
reparameterisation affords high flexibility in the choice of corresponding
proposals, but will only be applicable when the diffusion coefficient
$\sigma_{t}(s,x)$ is invertible. We refer to the second construct as the
‘backward transform’, as it is applicable to for the Feynman-Kac models that
make use of the backward proposals in Section 4.2. This direction does not
require invertibility of $\sigma_{t}(s,x)$ and will be applicable to a broader
class of signal SDEs for the CD-SSM, including hypoelliptic SDEs.

###  5.1 Defining $\bar{M}^{\rightarrow}_{t}$: Forward Transform

We consider the setting of the forward proposal
$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ in Section 4.1. Recall from (17) that
this kernel derives from the law of the proposal SDE in (18), extended with
the end point. We will define a new kernel
$\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]$, obtained from the forward
proposal $M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ via use of an invertible map,
the latter defined via an auxiliary bridge SDE.

The construction of this map will relate with the backwards decomposition of
$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ discussed in Remark 2 and expressed in
(21) in terms of the transition density $m_{t}^{\rightarrow}(e_{t}|e_{t-1})$
for the ending point and a connecting SDE bridge
$M_{t}^{\rightarrow}[(e_{t-1},e_{t}),dv_{t}]$. Indeed, we define an auxiliary
bridge SDE with distribution $R_{t}[(e_{t-1},e_{t}),dv_{t}]$ that dominates
the law of the bridge of the forward proposal, that is,
$M_{t}^{\rightarrow}[(e_{t-1},e_{t}),dv_{t}]\ll
R_{t}[(e_{t-1},e_{t}),dv_{t}]$. The auxiliary bridge SDE writes as follows:

$\displaystyle R_{t}[(e_{t-1},e_{t}),dv_{t}]:\qquad
dX_{t}(s)=a_{t}^{br}(s,X_{t}(s);y_{t})+\sigma_{t}(s,X_{t}(s))dB(s),\quad
X_{t}(0)=e_{t-1},$ (39)

for a drift function
$a_{t}^{br}(\cdot,\cdot;y_{t}):[0,\Delta_{t}]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$.
This drift must be such so that the solution of (39) is forced to hit the
ending point $e_{t}$. Following the earlier discussion about the Delyon-Hu
bridge, a simple choice of drift is the following:

$a_{t}^{br,\mathrm{DH}}(s,x;y_{t})=\frac{e_{t}-x}{\Delta_{t}-s}.$ (40)

We denote the corresponding bridge distribution for this choice as
$R_{t}^{\mathrm{DH}}[(e_{t-1},e_{t}),dv_{t}]$.

######  Remark 4.

The bridge SDE law $R_{t}^{\mathrm{DH}}[(e_{t-1},e_{t}),dv_{t}]$ coincides
with $M^{\leftarrow,\textrm{DH}}_{t}[(e_{t-1},e_{t}),dv_{t}]$ defined in the
context of setting up a backward proposal in Section 4.2. We note though that
this bridge SDE will now take up a very different role compared to the one in
Section 4.2. That is, here bridge SDE (39) will provide the means for a
transform of paths of the proposal SDE (18). In contrast, in Section 4.2,
$M^{\leftarrow,\textrm{DH}}_{t}[(e_{t-1},e_{t}),dv_{t}]$ formed part of the
proposal.

Following Delyon and Hu, (2006), we obtain a density function which is as the
one in (4.2.2), with the difference that we now work with the bridge
$M_{t}^{\rightarrow}[(e_{t-1},e_{t}),dv_{t}]$ of the proposal SDE (18) rather
than the bridge of the signal SDE, $P_{t}[(e_{t-1},e_{t}),dv_{t}]$. That is,
under the drift choice (40), and regularity conditions, we have:

$\frac{dM_{t}^{\rightarrow}[(e_{t-1},e_{t}),\cdot]}{dR_{t}^{\mathrm{DH}}[(e_{t-1},e_{t}),\cdot]}(v_{t})=\frac{\mathcal{N}(e_{t};e_{t-1},\Sigma_{t}(0,e_{t-1}))}{m_{t}^{\rightarrow}(e_{t}|e_{t-1})}\sqrt{\frac{\det(\Sigma_{t}(0,e_{t-1}))}{\det(\Sigma_{t}(\Delta_{t},e_{t}))}}\times\exp\big{\\{}\phi_{t}(v_{t})\big{\\}},$
(41)

for the functional:

$\displaystyle\phi_{t}(v_{t})=$ $\displaystyle\int_{0}^{\Delta
s}((b_{t}^{\rightarrow})^{\top}\Sigma_{t}^{-1})(s,v_{t}(s))dv_{t}(s)$
$\displaystyle\qquad-\tfrac{1}{2}\int_{0}^{\Delta
s}((b_{t}^{\rightarrow})^{\top}\Sigma_{t}^{-1}b_{t}^{\rightarrow})(s,v_{t}(s))ds$
$\displaystyle\qquad\qquad-\tfrac{1}{2}\int_{0}^{\Delta
s}\frac{e_{t}-v_{t}(s)}{\Delta_{t}-s}\diamond
d\big{\\{}\Sigma_{t}^{-1}(s,v_{t}(s))(e_{t}-v_{t}(s))\big{\\}}.$

It is now important to consider the auxiliary bridge SDE (39) as the means for
obtaining a map
$\mathfrak{F}_{t}(\cdot;e_{t-1},e_{t}):\mathsf{U}_{t}\to\mathsf{X}_{t}(e_{t})$
from the driving noise $B_{t}=\\{B(s)\\}_{s\in[0,\Delta_{t}]}$ to the solution
$X_{t}=\\{X_{t}(s)\\}_{s\in[0,\Delta_{t}]}$ of (39), so that:

$\displaystyle B_{t}\mapsto X_{t}=\mathfrak{F}_{t}(B_{t};e_{t-1},e_{t}).$ (42)

We have defined spaces $\mathsf{U}_{t}$ and $\mathsf{X}_{t}(e_{t})$, so that
the former corresponds to the set of realisations of a standard
$d$-dimensional Brownian motion on $[0,\Delta_{t}]$ and the latter to
realisations of the solution of (39), thus emphasising that paths must
terminate at $e_{t}$ in the case of $\mathsf{X}_{t}(e_{t})$. Recall that in
the setting of the forward approach we only consider $d_{w}=d$.

######  Assumption 1.

The diffusion coefficient matrix $\sigma_{t}(s,x)\in\mathbb{R}^{d\times d}$ is
invertible for all $(s,x)\in[0,\Delta_{t}]\times\mathbb{R}^{d}$,
$t\in\mathcal{T}$.

Under Assumption 1, we can make use of the inverse of map $\mathfrak{F}_{t}$
to transform paths from the proposal SDE (18). That is, for $v_{t}\sim
M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$ we set $e_{t}=\texttt{end}(v_{t})$ and
define the map $\mathfrak{F}_{t}^{-1}:\mathsf{X}_{t}(e_{t})\to\mathsf{U}_{t}$,
given $e_{t-1}$, $e_{t}$, as follows:

$\displaystyle u_{t}:=\mathfrak{F}_{t}^{-1}(v_{t};e_{t-1},e_{t}).$ (43)

Note that $u_{t}$ is, in practice, obtained via recursive application of the
following expression:

$\displaystyle
du_{t}(s)=\sigma_{t}(s,v_{t}(s))^{-1}\big{\\{}dv_{t}(s)-a_{t}^{br}(s,v_{t}(s))ds\big{\\}}.$
(44)

We now have all necessary ingredients for the definition of the transformed
proposal kernel $\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]$. In particular, we
have that:

$\displaystyle
z_{t}=(u_{t},e_{t})\in\mathsf{Z}_{t}=\mathsf{U}_{t}\times\mathbb{R}^{d},\qquad
t\in\mathcal{T}.$

Then, given $z_{t-1}=(u_{t-1},e_{t-1})$, we have that the variate
$z_{t}=(u_{t},e_{t})\sim\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]$ is obtained
as follows:

  * i.

Sample $v_{t}\sim M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$, with the kernel
defined via the proposal SDE (18) and set $e_{t}=\texttt{end}(v_{t})$.

  * ii.

Retrieve $u_{t}=\mathfrak{F}_{t}^{-1}(v_{t};e_{t-1},e_{t})$ via use of (44).

From the above specification of the reparameterised kernel
$\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]$, one can obtain the following
expression for its backward decomposition:

$\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]=\big{(}m_{t}^{\rightarrow}(e_{t}|e_{t-1})de_{t}\big{)}\otimes\bar{M}^{\rightarrow}_{t}[(e_{t-1},e_{t}),du_{t}].$
(45)

Recall that $m_{t}^{\rightarrow}(e_{t}|e_{t-1})$ is the transition density of
the forward proposal SDE (18). Then, the kernel
$\bar{M}^{\rightarrow}_{t}[(e_{t-1},e_{t}),du_{t}]$ is defined as the
distribution of $u_{t}\sim\mathfrak{F}^{-1}_{t}(v_{t};e_{t-1},e_{t})$, for the
path-variable $v_{t}\sim M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$ and for
$e_{t}=\texttt{end}(v_{t})$. This definition of $u_{t}$ has connections with
the ‘innovation process’ used in the context of MCMC algorithms, e.g. in Chib
et al., (2004); Golightly and Wilkinson, (2008).

Let $\mathbb{W}_{t}$ denote the law of a standard $d$-dimensional Brownian
motion on $[0,\Delta_{t}]$. Notice that application of the map
$\mathfrak{F}_{t}(\cdot;e_{t-1},e_{t})$ to a Brownian motion $B_{t}$ (resp. to
$u_{t}\sim\bar{M}^{\rightarrow}_{t}[(e_{t-1},e_{t}),du_{t}]$) will yield a
random path distributed according to the solution of the bridge SDE (39)
(resp. distributed according to
$M^{\rightarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]$). Thus, we immediately obtain
the following density, due to the 1–1 transformation of the random variates
involved:

$\frac{d\bar{M}^{\rightarrow}_{t}[(e_{t-1},e_{t}),\cdot]}{d\mathbb{W}_{t}}(u_{t})=\frac{dM_{t}^{\rightarrow}[(e_{t-1},e_{t}),\cdot]}{dR_{t}[(e_{t-1},e_{t}),\cdot]}\big{(}\mathfrak{F}_{t}(u_{t};e_{t-1},e_{t})\big{)}.$
(46)

For the case when the Delyon-Hu drift function (40) is used in the
derivations, the density on the right-hand-side of (46) is as provided in
(41). Continuing from (46), the complete proposal kernel
$\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]$ admits the following Radon-Nikodym
derivative (recall that $z_{t}=(e_{t},u_{t})$):

$\displaystyle\bar{M}^{\rightarrow}_{t}(z_{t}|z_{t-1})$
$\displaystyle=\frac{\bar{M}^{\rightarrow}_{t}[z_{t-1},dz_{t}]}{\mathrm{Leb}(de_{t})\otimes\mathbb{W}_{t}[du_{t}]}$
$\displaystyle=m_{t}^{\rightarrow}(e_{t}|e_{t-1})\times\frac{dM^{\rightarrow}_{t}[(e_{t-1},e_{t}),\cdot]}{dR_{t}[(e_{t-1},e_{t}),\cdot]}\big{(}\mathfrak{F}_{t}(u_{t};e_{t-1},e_{t})\big{)}.$
(47)

Notice that the transition density $m_{t}^{\rightarrow}(e_{t}|e_{t-1})$ of the
proposal SDE (18) cancels out when taking the product in the right-hand-side
of (47), in the case
$R_{t}[(e_{t-1},e_{t}),dv_{t}]=R_{t}^{\mathrm{DH}}[(e_{t-1},e_{t}),dv_{t}]$,
thus the provided density for $\bar{M}_{t}^{\rightarrow}[z_{t-1},dz_{t}]$ is
analytically tractable.

###  5.2 Defining $\bar{M}_{t}^{\leftarrow}$: Backward Transform

We present a method for the construction of a transform
$\bar{M}_{t}^{\leftarrow}[z_{t-1},dz_{t}]$ of the backward proposal
$M^{\leftarrow}_{t}[x_{t-1},dx_{t}]$, the latter specified via the backward
decomposition in (26), consisted of the kernels
$m_{t}^{\leftarrow}(e_{t}|e_{t-1})$ and
$M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$.

Nothing needs to be done about $m_{t}^{\leftarrow}(e_{t}|e_{t-1})$. Moving to
kernel $M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$ corresponding to the
proposal bridge SDE (31), we define the map
$\mathfrak{H}_{t}(\cdot;e_{t-1},e_{t}):\mathsf{U}_{t}\rightarrow\mathsf{X}_{t}(e_{t})$
that projects the driving noise $B_{t}=\\{B(s)\\}_{s\in[0,\Delta_{t}]}$ of
(31) onto the solution of (31). We now denote the solution of the bridge SDE
(31) as
$X^{\leftarrow,br}_{t}=\\{X(s)^{\leftarrow,br}\\}_{s\in[0,\Delta_{t}]}$ for
clarity. That is, we have:

$\displaystyle B_{t}\mapsto
X^{\leftarrow,br}_{t}=\mathfrak{H}_{t}(B_{t};e_{t-1},e_{t}).$

The definition resembles the one of $\mathfrak{F}_{t}(\cdot;\cdot)$ in (42),
used in the case of the forward reparameterisation. For $v_{t}\sim
M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$, the backward reparameterisation
considers instead the driving Wiener noise, $u_{t}$, in the construction of
$v_{t}$ via $v_{t}=\mathfrak{H}(u_{t};e_{t-1},e_{t})$. In this case: (i) one
simply has that $u_{t}\sim\mathbb{W}_{t}$; (ii) no inversion of
$u_{t}\mapsto\mathfrak{H}(u_{t};e_{t-1},e_{t})$ is required. This second
characteristic of the backward approach is important, as it sidesteps the
requirement for invertibility of the diffusion coefficient in the SDE (1) for
the latent signal, thus allowing application of the backward method to hypo-
elliptic models.

Overall, we define $\bar{M}^{\leftarrow}_{t}[z_{t-1},dz_{t}]$ as follows, for
$z_{t}=(e_{t},u_{t})$ given $z_{t-1}=(e_{t-1},u_{t-1})$:

$\bar{M}_{t}^{\leftarrow}[z_{t-1},dz_{t}]=\big{(}m_{t}^{\leftarrow}(e_{t}|e_{t-1})de_{t}\big{)}\otimes\mathbb{W}_{t}[dz_{t}].$
(48)

The definition of $\bar{M}_{t}^{\leftarrow}[z_{t-1},dz_{t}]$ immediately
implies that one can sample from this kernel, by independently generating a
proposed Brownian motion path $z_{t}\sim\mathbb{W}_{t}$ and a proposed end-
point $e_{t}\sim m_{t}^{\leftarrow}(e_{t}|e_{t-1})$. Further, this proposal
kernel has a tractable density given by:

$\bar{M}_{t}^{\leftarrow}(z_{t}|z_{t-1})=\frac{\bar{M}_{t}^{\leftarrow}[z_{t-1},dz_{t}]}{\mathrm{Leb}(de_{t})\otimes\mathbb{W}_{t}[du_{t}]}=m_{t}(v_{t}|v_{t-1}).$
(49)

######  Remark 5.

We note a main difference between backward and forward proposals,
$M_{t}^{\rightarrow}[x_{t-1},dx_{t}]$ and
$M_{t}^{\leftarrow}[x_{t-1},dx_{t}]$, that ultimately leads to the method
based on $M_{t}^{\leftarrow}[x_{t-1},dx_{t}]$ being applicable for a non-
invertible diffusion coefficient $\sigma_{t}(s,x)$. Kernel
$M^{\leftarrow}_{t}[x_{t-1},dx_{t}]$ is constructed in (26) via direct
consideration of the backward decomposition of the target
$P_{t}[x_{t-1},dx_{t}]$ in (3.2). That is, one selects _tractable_ laws
$m_{t}^{\leftarrow}(e_{t}|e_{t-1})$ and
$M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$, both for the proposed ending
point $e_{t}$ and the proposed bridge SDE. Then, the transform
$(e_{t},v_{t})\mapsto(e_{t},u_{t})$, with $(e_{t},u_{t})$ admitting a density
w.r.t. the reference measure
$\mathrm{Leb}(de_{t})\otimes\mathbb{W}_{t}[du_{t}]$, is obtained by working
with the noise used to sample $v_{t}\sim
M_{t}^{\leftarrow}[(e_{t-1},e_{t}),dv_{t}]$. Kernel
$M^{\rightarrow}_{t}[x_{t-1},dx_{t}]$ in (17) relates to the forward
decomposition of $P_{t}[x_{t-1},dx_{t}]$ in (10), thus one must now work with
$v_{t}\sim M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$. However, use of the driving
noise for $v_{t}$ cannot lead to non-degenerate particle-based smoothing
methods. Instead, one must rely to the auxiliary bridge SDE, of law
$R_{t}[(e_{t-1},e_{t}),dv_{t}]$, which serves to transform the signal
$(e_{t},v_{t})\mapsto(e_{t},u_{t})$, via a map that indeed requires the
inverse of $\sigma_{t}(s,x)$.

###  5.3 Potentials $\bar{G}_{t}^{\rightarrow}$, $\bar{G}_{t}^{\leftarrow}$

We have now defined two classes of tractable kernels
$\bar{M}_{t}^{\rightarrow}[z_{t-1},dz_{t}]$,
$\bar{M}^{\leftarrow}_{t}[z_{t-1},dz_{t}]$ that we can sample from, and which
have analytical densities, given in (47) and (49) w.r.t. the reference measure
$\mathrm{Leb}(de_{t})\otimes\mathbb{W}_{t}[du_{t}]$. Critically, the reference
measure does not involve the end-point $e_{t-1}$. We now complete our
construction of the transformed Feynman-Kac models, by determining the
corresponding potentials, denoted $\bar{G}_{t}^{\rightarrow}(z_{t-1},z_{t})$
and $\bar{G}_{t}^{\leftarrow}(z_{t-1},z_{t})$.

####  5.3.1 Forward Method, $\bar{G}_{t}^{\rightarrow}$

Recall that we work with the dynamics of the proposal SDE (18) to obtain the
path-valued variate $v_{t}\sim M_{t}^{\rightarrow}[e_{t-1},dv_{t}]$. Upon
taking $u_{t}=\mathfrak{F}_{t}^{-1}(v_{t};e_{t-1},e_{t})$ with
$e_{t}=\texttt{end}(v_{t})$, one obtains
$z_{t}=(e_{t},v_{t})\sim\bar{M}_{t}^{\rightarrow}[z_{t-1},dz_{t}]$. The
transform:

$\displaystyle x_{t}=(e_{t},v_{t})\mapsto z_{t}=\big{(}e_{t},u_{t});\qquad
u_{t}=\mathfrak{F}_{t}^{-1}(v_{t};e_{t-1},e_{t}),$

when applied on the target signal variate $x_{t}\sim P_{t}[x_{t-1},dx_{t}]$
itself will similarly give rise to a probability law
$\bar{P}_{t}^{\rightarrow}[z_{t-1},dz_{t}]$ on the $z$-space. Thus, the
potential $\bar{G}_{t}^{\rightarrow}(z_{t-1},z_{t})$ on the $z$-space is equal
to:

$\displaystyle\bar{G}_{t}^{\rightarrow}(z_{t-1},z_{t})=\frac{d\bar{P}_{t}^{\rightarrow}[z_{t-1},\cdot]}{d\bar{M}^{\rightarrow}_{t}[z_{t-1},\cdot]}(z_{t})\times
f_{t}(y_{t}|e_{t}),$ (50)

for the Radon-Nikodym derivative:

$\displaystyle\frac{d\bar{P}_{t}^{\rightarrow}[z_{t-1},\cdot]}{d\bar{M}_{t}^{\rightarrow}[z_{t-1},\cdot]}(z_{t})$
$\displaystyle=\frac{p_{t}(e_{t}|e_{t-1})}{m_{t}^{\rightarrow}(e_{t}|e_{t-1})}\times\frac{d\bar{P}_{t}[(e_{t-1},e_{t}),\cdot]}{d\bar{M}^{\rightarrow}_{t}[(e_{t-1},e_{t}),\cdot]}(u_{t})$
$\displaystyle=\frac{p_{t}(e_{t}|e_{t-1})}{m_{t}^{\rightarrow}(e_{t}|e_{t-1})}\times\frac{dP_{t}[(e_{t-1},e_{t}),\cdot]}{dM^{\rightarrow}_{t}[(e_{t-1},e_{t}),\cdot]}\big{(}\mathfrak{F}_{t}(u_{t};e_{t-1},e_{t})\big{)}$
$\displaystyle\equiv\frac{dP_{t}[e_{t-1},\cdot]}{dM^{\rightarrow}_{t}[e_{t-1},\cdot]}\big{(}\mathfrak{F}_{t}(u_{t};e_{t-1},e_{t})\big{)}.$
(51)

The expression in the last line of (51) can be analytically evaluated up-to
time-discretisation via use of the Girsanov density in (19).

####  5.3.2 Backward Method, $\bar{G}_{t}^{\leftarrow}$

In this case, we work with the transform:

$\displaystyle z_{t}=\big{(}e_{t},u_{t})\mapsto x_{t}=(e_{t},v_{t});\qquad
v_{t}=\mathfrak{H}_{t}(u_{t};e_{t-1},e_{t}).$

When considering the target signal path $x_{t}\sim P_{t}[x_{t-1},dx_{t}]$, we
have that $x_{t}=(e_{t},v_{t})$ with laws $e_{t}\sim
p_{t}(e_{t}|e_{t-1})de_{t}$ and $v_{t}\sim P_{t}[(e_{t-1},e_{t}),dv_{t}]$.
Then, $u_{t}\sim\bar{P}_{t}^{\leftarrow}[(e_{t-1},e_{t}),du_{t}]$, with this
latter distribution being defined via the requirement
$v_{t}=\mathfrak{H}_{t}(u_{t};e_{t-1},e_{t})$ must follow
$P_{t}[(e_{t-1},e_{t}),dv_{t}]$.

Thus, the potential on the $z$-space writes as follows:

$\displaystyle\bar{G}_{t}^{\leftarrow}(z_{t-1},z_{t})$
$\displaystyle=\frac{d\bar{P}_{t}^{\leftarrow}[z_{t-1},\cdot]}{d\bar{M}^{\leftarrow}_{t}[z_{t-1},\cdot]}(z_{t})\times
f_{t}(y_{t}|e_{t})$
$\displaystyle=\frac{p_{t}(e_{t}|e_{t-1})}{m_{t}^{\leftarrow}(e_{t}|e_{t-1})}\times\frac{d\bar{P}_{t}^{\leftarrow}[(e_{t-1},e_{t}),\cdot]}{d\bar{M}^{\leftarrow}_{t}[(e_{t-1},e_{t}),\cdot]}(u_{t})\times
f_{t}(y_{t}|e_{t})$
$\displaystyle=\frac{p_{t}(e_{t}|e_{t-1})}{m_{t}^{\leftarrow}(e_{t}|e_{t-1})}\times\frac{dP_{t}[(e_{t-1},e_{t}),\cdot]}{dM^{\leftarrow}_{t}[(e_{t-1},e_{t}),\cdot]}\big{(}\mathfrak{H}_{t}(u_{t};e_{t-1},e_{t})\big{)}\times
f_{t}(y_{t}|e_{t}).$ (52)

The Radon-Nikodym derivative appearing on the right-hand-side of the last line
in the above expression is given in (4.2.2) for proposal
$M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]=M^{\leftarrow,\mathrm{DH}}_{t}[(e_{t-1},e_{t}),dv_{t}]$,
and in (37) for the case when
$M^{\leftarrow}_{t}[(e_{t-1},e_{t}),dv_{t}]=M_{t}^{\leftarrow,\mathrm{G}}[(e_{t-1},e_{t}),dv_{t}]$.

###  5.4 Complete Forward and Backward Methods

We provide in Algorithms 3 and 4 pseudo-codes for the particle filter and FFBS
based on the the Forward Method. Then, Algorithms 5 and 6 provide similar
pseudo-codes for the Backward Method.

Input: CD-SSM in (1)-(2).

Output: Particles $z^{1:N}_{1:T}=(e_{t}^{1:N},u_{t}^{1:N})_{t\in[1:T]}$,
$x^{1:N}_{1:T}=(e_{t}^{1:N},v_{t}^{1:N})_{t\in[1:T]}$, ancestors
$A^{1:N}_{2:T}$, weights $W_{1:T}^{1:N}$.

1 for  _$t=1,...,T$_ do

2 If $t=1$, set $A_{t}^{j}\leftarrow j$; else
$A_{t}^{j}\sim\texttt{Resample}(W^{1:N}_{t-1})$;

3 Sample $v^{j}_{t}$ from the proposal SDE (18) with starting point
$e_{t-1}^{A_{t}^{j}}$ ;

4 Set $e_{t-1}^{j}=\texttt{start}(v_{t}^{j})$ and
$e_{t}^{j}=\texttt{end}(v_{t}^{j})$;

5 Apply transform:
$u^{j}_{t}=\mathfrak{F}_{t}^{-1}(v^{j}_{t};e_{t-1}^{j},e_{t}^{j})$;

6 Collect particles: $z^{j}_{t}=(e_{t}^{j},u^{j}_{t})$,
$x^{j}_{t}=(e_{t}^{j},v^{j}_{t})$ ;

7 Assign weights:
$w_{t}^{j}=\bar{G}_{t}^{\rightarrow}(z_{t-1}^{A_{t}^{j}},z_{t}^{j})$ using the
expressions in (50)-(51) for $\bar{G}^{\rightarrow}_{t}(\cdot,\cdot)$;

8 Normalise weights: $W_{t}^{j}=\frac{w_{t}^{j}}{\sum_{k=1}^{N}w_{t}^{k}}$;

9 end for

Algorithm 3 Particle Filter for CD-SSM in (1)-(2) – under Forward Method

Input: Output $z^{1:N}_{1:T}=(e_{t}^{1:N},u_{t}^{1:N})_{t\in[1:T]}$,
$A_{2:T}^{1:N}$, $W_{1:T}^{1:N}$, from Algorithm 3.

Output: Approximate sample
$\big{(}x_{1}^{B_{1}}=(e_{1},v_{1})^{B_{1}},\ldots,x_{T}^{B_{T}}=(e_{T},v_{T})^{B_{T}}\big{)}$
from the smoothing distribution $P\,[\,X_{1:T}\in
dx_{1:T}\,|\,Y_{1:T}=y_{1:T}\,]$.

1 Sample $B_{T}\sim\mathcal{M}(W_{T}^{1:N})$;

2 for  _$t=T,...,2$_ do

3 Assign weights:
$\hat{w}_{t-1}^{j}\leftarrow\bar{M}_{t}^{\rightarrow}(z_{t}^{B_{t}}|z_{t-1}^{j})\,\bar{G}^{\rightarrow}_{t}(z_{t-1}^{j},z_{t}^{B_{t}})\cdot
W_{t-1}^{j}$ using expression (47) for
$\bar{M}^{\rightarrow}_{t}(\cdot|\cdot)$;

4 Normalise weights:
$\hat{W}_{t-1}^{j}=\frac{\hat{w}_{t-1}^{j}}{\sum_{k=1}^{N}\hat{w}_{t-1}^{k}}$;

5 $B_{t-1}\sim\mathcal{M}(\hat{W}_{t}^{1:N})$;

6 end for

Return
$x_{1}^{B_{1}}=\big{(}e_{1}^{B_{1}},\mathfrak{F}_{1}(u^{B_{1}}_{1};e_{0},e_{1}^{B_{1}})\big{)}$,
$x_{2}^{B_{2}}=\big{(}e_{2}^{B_{2}},\mathfrak{F}_{2}(u^{B_{2}}_{2};e_{1}^{B_{1}},e_{2}^{B_{2}})\big{)}$,
…,
$x_{T}^{B_{T}}=\big{(}e_{T}^{B_{T}},\mathfrak{F}_{T}(u^{B_{T}}_{T};e_{T-1}^{B_{T-1}},e_{T}^{B_{T}})\big{)}$.

Algorithm 4 FFBS for CD-SSM in (1)-(2) – under Forward Method

Input: CD-SSM in (1)-(2).

Output: Particles $z^{1:N}_{1:T}=(e_{t}^{1:N},u_{t}^{1:N})_{t\in[1:T]}$,
$x^{1:N}_{1:T}=(e_{t}^{1:N},v_{t}^{1:N})_{t\in[1:T]}$, ancestors
$A^{1:N}_{2:T}$, weights $W_{1:T}^{1:N}$.

1 for  _$t=1,...,T$_ do

2 If $t=1$ set $A_{t}^{j}\leftarrow j$; else
$A_{t}^{j}\sim\texttt{Resample}(W^{1:N}_{t-1})$;

3 Sample $e^{j}_{t}\sim m^{\leftarrow}_{t}(e_{t}|e_{t-1}^{A_{t}^{j}})$ and
$u_{t}^{j}\sim\mathbb{W}_{t}[du_{t}]$ ;

4 Obtain the path:
$v_{t}^{j}\leftarrow\mathfrak{h}_{t}(u_{t}^{j};e_{t-1}^{A_{t}^{j}},e_{t}^{j})$
;

5 Collect particles: $z^{j}_{t}=(e_{t}^{j},u^{j}_{t})$,
$x^{j}_{t}=(e_{t}^{j},v^{j}_{t})$;

6 Assign weights:
$w_{t}^{j}=\bar{G}_{t}^{\leftarrow}(z_{t-1}^{A_{t}^{j}},z_{t}^{j})$ using the
expression in (52) for $\bar{G}^{\leftarrow}_{t}(\cdot,\cdot)$;

7 Normalise weights: $W_{t}^{j}=\frac{w_{t}^{j}}{\sum_{k=1}^{N}w_{t}^{k}}$;

8 end for

Algorithm 5 Particle Filter for CD-SSM in (1)-(2) – under Backward Method

Input: Full output, $z^{1:N}_{1:T}=(e_{t}^{1:N},u_{t}^{1:N})_{t\in[1:T]}$,
$A_{2:T}^{1:N}$, $W_{1:T}^{1:N}$, from Algorithm 5.

Output: Approximate sample
$\big{(}x_{1}^{B_{1}}=(e_{1},v_{1})^{B_{1}},\ldots,x_{T}^{B_{T}}=(e_{T},u_{T})^{B_{T}}\big{)}$
from the smoothing distribution $P\,[\,X_{1:T}\in
dx_{1:T}\,|\,Y_{1:T}=y_{1:T}\,]$.

1 Sample $B_{T}\sim\mathcal{M}(W_{T}^{1:N})$;

2 for  _$t=T,...,2$_ do

3 Assign weights:
$\hat{w}_{t-1}^{j}\leftarrow\bar{M}_{t}^{\leftarrow}(z_{t}^{B_{t}}|z_{t-1}^{j})\,\bar{G}^{\leftarrow}_{t}(z_{t-1}^{j},z_{t}^{B_{t}})\cdot
W_{t-1}^{j}$ using the expression in (49) for
$\bar{M}^{\leftarrow}_{t}(\cdot|\cdot)$;

4 Normalise weights:
$\hat{W}_{t-1}^{j}=\frac{\hat{w}_{t-1}^{j}}{\sum_{k=1}^{N}\hat{w}_{t-1}^{k}}$;

5 $B_{t-1}\sim\mathcal{M}(\hat{W}_{t}^{1:N})$;

6 end for

Return
$x_{1}^{B_{1}}=\big{(}e_{1}^{B_{1}},\mathfrak{H}_{1}(u^{B_{1}}_{1};e_{0},e_{1}^{B_{1}})\big{)}$,
$x_{2}^{B_{2}}=\big{(}e_{2}^{B_{2}},\mathfrak{H}_{2}(u^{B_{2}}_{2};e_{1}^{B_{1}},e_{2}^{B_{2}})\big{)}$,
…,
$x_{T}^{B_{T}}=\big{(}e_{T}^{B_{T}},\mathfrak{H}_{T}(u^{B_{T}}_{T};e_{T-1}^{B_{T-1}},e_{T}^{B_{T}})\big{)}$.

Algorithm 6 FFBS for CD-SSM in (1)-(2) – under Backward Method

######  Remark 6.

A key point in the implementation of Algorithm 6 is the calculation of the
backward weights $\hat{w}_{t-1}^{j}$ for the current particle
$z_{t}^{B_{t}}=(e_{t}^{B_{t}},u_{t}^{B_{t}})$. Implementation of the
methodology has allowed for all particles
$z_{t-1}^{j}=(e_{t-1}^{j},u_{t-1}^{j})$ to be given positive backward weights.

##  6 Conclusions

We have presented a comprehensive methodology over particle-based algorithms
for carrying out effective statistical inference for CD-SSMs. By allowing for
careful alternation of ancestors within algorithmic executions, our particle-
based methods for CD-SSMs overcome path degeneracy characterising standard
approaches. Further advances in the development of guided proposals for the
filtering problem for CD-SSMs can be directly incorporated in our framework,
see e.g. Mider et al., (2021); Chopin et al., (2023) for some recent
contributions in this direction. We have focused on continuous-time stochastic
signals defined via SDEs, and have allowed the latter to be members of the
elliptic or hypoelliptic class of diffusion processes. We aim to extend our
approach to signals determined by more general classes of continuous-time
processes, including SDEs with with jumps and discrete-values pure-jump
processes.

## References

  * Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.
  * Beskos et al., (2013) Beskos, A., Pillai, N., Roberts, G., Sanz-Serna, J.-M., and Stuart, A. (2013).  Optimal tuning of the hybrid Monte Carlo algorithm.  Bernoulli, 19(5A):1501–1534.
  * Bierkens et al., (2020) Bierkens, J., Van Der Meulen, F., and Schauer, M. (2020).  Simulation of elliptic and hypo-elliptic conditional diffusions.  Advances in Applied Probability, 52(1):173–212.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1).
  * Chib et al., (2004) Chib, S., Pitt, M. K., and Shephard, N. (2004).  Likelihood based inference for diffusion driven models.
  * Chopin et al., (2023) Chopin, N., Fulop, A., Heng, J., and Thiery, A. H. (2023).  Computational Doob h-transforms for online filtering of discretely observed diffusions.  In International Conference on Machine Learning, pages 5904–5923. PMLR.
  * Chopin and Papaspiliopoulos, (2020) Chopin, N. and Papaspiliopoulos, O. (2020).  An introduction to sequential Monte Carlo, volume 4.  Springer.
  * Del Moral, (2004) Del Moral, P. (2004).  Feynman-Kac formulae: genealogical and interacting particle systems with applications.  Springer.
  * Del Moral, (2013) Del Moral, P. (2013).  Mean field simulation for Monte Carlo integration.  Monographs on Statistics and Applied Probability, 126:26.
  * Delyon and Hu, (2006) Delyon, B. and Hu, Y. (2006).  Simulation of conditioned diffusion and application to parameter estimation.  Stochastic Processes and their Applications, 116(11):1660–1675.
  * Ditlevsen and Samson, (2019) Ditlevsen, S. and Samson, A. (2019).  Hypoelliptic diffusions: filtering and inference from complete and partial observations.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 81(2):361–384.
  * Doucet et al., (2000) Doucet, A., Godsill, S., and Andrieu, C. (2000).  On sequential Monte Carlo sampling methods for Bayesian filtering.  Statistics and Computing, 10:197–208.
  * Gelman et al., (1997) Gelman, A., Gilks, W. R., and Roberts, G. O. (1997).  Weak convergence and optimal scaling of random walk Metropolis algorithms.  The Annals of Applied Probability, 7(1):110–120.
  * Godsill et al., (2004) Godsill, S. J., Doucet, A., and West, M. (2004).  Monte Carlo smoothing for nonlinear time series.  Journal of the American Statistical Association, 99(465):156–168.
  * Golightly and Wilkinson, (2008) Golightly, A. and Wilkinson, D. J. (2008).  Bayesian inference for nonlinear multivariate diffusion models observed with error.  Computational Statistics & Data Analysis, 52(3):1674–1693.
  * Gordon et al., (1993) Gordon, N. J., Salmond, D. J., and Smith, A. F. (1993).  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  In IEE proceedings F (radar and signal processing), volume 140, pages 107–113. IET.
  * Johannes and Polson, (2010) Johannes, M. and Polson, N. (2010).  MCMC methods for continuous-time financial econometrics.  In Handbook of Financial Econometrics: Applications, pages 1–72. Elsevier.
  * Karjalainen et al., (2023) Karjalainen, J., Lee, A., Singh, S. S., and Vihola, M. (2023).  Mixing time of the conditional backward sampling particle filter.  arXiv preprint arXiv:2312.17572.
  * Lindsten et al., (2014) Lindsten, F., Jordan, M. I., and Schon, T. B. (2014).  Particle Gibbs with ancestor sampling.  Journal of Machine Learning Research, 15:2145–2184.
  * Llopis et al., (2018) Llopis, F. P., Kantas, N., Beskos, A., and Jasra, A. (2018).  Particle filtering for stochastic Navier–Stokes signal observed with linear additive noise.  SIAM Journal on Scientific Computing, 40(3):A1544–A1565.
  * Mider et al., (2021) Mider, M., Schauer, M., and Van der Meulen, F. (2021).  Continuous-discrete smoothing of diffusions.  Electronic Journal of Statistics, 15(2):4295–4342.
  * Oksendal, (2013) Oksendal, B. (2013).  Stochastic Differential Equations: An introduction with applications.  Springer Science & Business Media.
  * Pokern et al., (2009) Pokern, Y., Stuart, A. M., and Wiberg, P. (2009).  Parameter estimation for partially observed hypoelliptic diffusions.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 71(1):49–73.
  * Roberts and Rosenthal, (1998) Roberts, G. O. and Rosenthal, J. S. (1998).  Optimal scaling of discrete approximations to Langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268.
  * Roberts and Stramer, (2001) Roberts, G. O. and Stramer, O. (2001).  On inference for partially observed nonlinear diffusion models using the Metropolis–Hastings algorithm.  Biometrika, 88(3):603–621.
  * Schauer et al., (2017) Schauer, M., Van Der Meulen, F., and Van Zanten, H. (2017).  Guided proposals for simulating multi-dimensional diffusion bridges.  Bernoulli, 23(4A):2917–2950.
  * Sottinen and Särkkä, (2008) Sottinen, T. and Särkkä, S. (2008).  Application of Girsanov theorem to particle filtering of discretely observed continuous-time non-linear systems.  Bayesian Analysis, 3(3):555–584.
  * van der Meulen and Schauer, (2017) van der Meulen, F. and Schauer, M. (2017).  Bayesian estimation of discretely observed multi-dimensional diffusion processes using guided proposals.  Electronic Journal of Statistics, 11(1):2358–2396.
  * Yonekura and Beskos, (2022) Yonekura, S. and Beskos, A. (2022).  Online smoothing for diffusion processes observed with noise.  Journal of Computational and Graphical Statistics, 31(4):1344–1360.