On Unbiased Score Estimation for Partially Observed Diffusions

Jeremy HengESSEC Business School<EMAIL_ADDRESS>Jeremie HoussineauUniversity of Warwick<EMAIL_ADDRESS>and Ajay JasraComputer, Electrical and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology<EMAIL_ADDRESS>
We consider the problem of statistical inference for a class of partially-observed diffusion processes,
with discretely-observed data and finite-dimensional parameters.
We construct unbiased estimators of the score function, i.e. the gradient of the log-likelihood function with respect to
parameters, with no time-discretization bias.
These estimators can be straightforwardly employed within stochastic gradient methods to perform maximum likelihood
estimation or Bayesian inference.
As our proposed methodology only requires access to a time-discretization scheme such as the Euler–Maruyama method,
it is applicable to a wide class of diffusion processes and observation models.
Our approach is based on a representation of the score as a smoothing expectation using Girsanov theorem, and
a novel adaptation of the randomization schemes developed in Mcleish, 2011, Rhee and Glynn, 2015, Jacob et al., 2020.
This allows one to remove the time-discretization bias and burn-in bias when computing smoothing expectations using the conditional particle filter of Andrieu et al., 2010.
Central to our approach is the development of new couplings of multiple conditional particle filters.
We prove under assumptions that our estimators are unbiased and have finite variance.
The methodology is illustrated on several challenging applications from population ecology and neuroscience.

Partially-observed diffusions, unbiased estimation, particle filters, coupling, stochastic gradient methods.

§ INTRODUCTION

We consider a diffusion process $(X_t)_{t\geq0}$ in $\mathbb{R}^d$, defined as the solution of
the stochastic differential equation (SDE)

\begin{align}\label{eq:intro_diff}
dX_t = a_{\theta}(X_t)dt + \sigma(X_t)dW_t,\quad X_0=x_{\star}\in\mathbb{R}^d,
\end{align}

where $(W_t)_{t\geq 0}$ is a standard Brownian motion in $\mathbb{R}^d$.
We will assume that the drift function $a:\Theta\times\mathbb{R}^d \rightarrow\mathbb{R}^d$
depends on a vector of unknown parameters $\theta\in\Theta\subseteq\mathbb{R}^{d_{\theta}}$,
but the symmetric diffusion coefficient $\sigma:\mathbb{R}^d \rightarrow\mathbb{R}^{d\times d}$ does not.
The general case where the diffusion coefficient is parameter dependent requires a different
treatment; we will discuss how to adapt our proposed methodology in Section <ref>.
The drift and diffusion coefficients are assumed to be regular enough for (<ref>)
to admit a (weakly) unique solution for all $t>0$; more precise regularity conditions needed in our
analysis will be stated.
We model observations $Y_{t_1},\dots,Y_{t_P}\in\mathbb{R}^{d_y}$ at a collection of time instances
$0 \leq t_1< \cdots < t_P = T$ as conditionally independent given the latent diffusion process
$X=(X_t)_{0\leq t\leq T}$, with observation density $g:\Theta\times\mathbb{R}^d\times\mathbb{R}^{d_y}\rightarrow\mathbb{R}^+$.
For notational ease, we assume that observations are given at unit times, i.e. $t_p = p$ for $p\in\{1,\ldots,P\}$ and $P=T$,
which covers the case of regularly observed data by a time rescaling.
Irregularly observed data can also be accommodated with minor modifications to our presentation
and considered for an application in Section <ref>.
We note that this partially observed setting is distinct from the discretely observed case
where the diffusion process is observed without error [Sørensen, 2004].
Extension of our methodology to continuous-time observation models is straightforward and illustrated
on an application in Section <ref>.

Given a realization $y_{1:T}=(y_{t})_{t=1}^T$ of the observation process, the marginal likelihood of our state
space model is

\begin{align}\label{eqn:marginal_likelihood}
p_{\theta}(y_{1:T}) = \mathbb{E}_{\theta}\left[\prod_{t=1}^T g_{\theta}(y_{t}|X_{t})\right],
\end{align}

where $\mathbb{E}_{\theta}$ denotes expectation w.r.t. the probability measure $\mathbb{P}_{\theta}$,
induced by the solution of (<ref>) on $\mathcal{C}_d([0,T])$, the space of continuous mappings from $[0,T]$ to $\mathbb{R}^d$.
To perform parameter inference using gradient-based algorithms, we shall consider approximations of
the score function $S(\theta)=\nabla_{\theta}\log p_{\theta}(y_{1:T})$, where $\nabla_{\theta}$ denotes the gradient
w.r.t. the parameter $\theta$.
If estimators $\widehat{S}(\theta)$ of $S(\theta)$ can be constructed,
one can approximate the maximum likelihood estimator (MLE) $\hat{\theta}\in \arg\max_{\theta\in\Theta} p_{\theta}(y_{1:T})$
with stochastic gradient ascent (SGA)

\begin{align}\label{eqn:SGA}
\theta_m = \theta_{m-1} + \varepsilon_m\widehat{S}(\theta_{m-1}), \quad m=1,2,\ldots,
\end{align}

where $(\varepsilon_m)_{m=1}^{\infty}$ is a sequence of learning rates. Similarly, in the Bayesian framework, one can
sample from the posterior distribution $p(d\theta|y_{1:T})\propto p(\theta)p_{\theta}(y_{1:T})d\theta$ using
the stochastic gradient Langevin dynamics (SGLD) [Welling and Teh, 2011]

\begin{align}\label{eqn:SGLD}
\theta_m = \theta_{m-1} + \frac{1}{2}\varepsilon_m\left(\nabla_{\theta}\log p(\theta_{m-1}) +
\widehat{S}(\theta_{m-1})\right) + \varepsilon_m^{1/2}\eta_m, \quad m=1,2,\ldots,
\end{align}

where $(\eta_m)_{m=1}^{\infty}$ is a sequence of independent standard Gaussian random vectors in $\mathbb{R}^{d_{\theta}}$.
In both stochastic gradient methods, constraints on parameters can be dealt with using a suitable transformation.
Under mild regularity conditions, we will see that Girsanov theorem can be applied to represent the score as a smoothing expectation

\begin{align}\label{eqn:score_function}
S(\theta) =
\mathbb{E}_{\theta}\left[ G_{\theta}(X) ~|~ y_{1:T} \right],
\end{align}

for some functional $G_{\theta}:\mathcal{C}_d([0,T])\rightarrow\mathbb{R}^{d_{\theta}}$.
For most problems of practical interest, statistical inference is a challenging task due to the intractability
of $p_{\theta}(y_{1:T})$ and $S(\theta)$.

The latter is typically the case for two reasons.
Firstly, as most diffusion processes do not have analytically tractable transition densities,
one has to resort to time-discretization.
Secondly, even if transition densities are available or a time-discretization scheme is employed,
the expectation over the latent process is usually intractable, and one has to rely on Monte Carlo approximations.
The following describes some existing approaches to these two issues, and how we will resolve them in this article.

To perform inference and state estimation without any time-discretization bias,
Beskos et al., 2006, Beskos et al., 2009, Fearnhead et al., 2008 adopted
the methodology developed in Beskos and Roberts, 2005, Beskos et al., 2006 that allows exact
simulation of diffusion sample paths; see also Pollock et al., 2016, Blanchet and Zhang, 2017, Fearnhead et al., 2017
for recent advances. Although these simulation techniques are elegant, they require a user to know
various properties of the diffusion process, which may not always be the case in practice.
An alternative approach to deal with the time-discretization bias, which we shall adopt,
is to employ the debiasing schemes proposed in
Mcleish, 2011, Rhee and Glynn, 2015.
These methods allow one to unbiasedly estimate an expectation of a functional w.r.t. an infinite-dimensional law
by randomizing over the level of the time-discretization in a type of multilevel Monte Carlo (MLMC) approach.
The use of unbiased estimators of the score function $S(\theta)$ within stochastic gradient methods
is particularly appealing. Under suitable assumptions and an appropriate choice of learning rates $(\varepsilon_m)_{m=1}^{\infty}$,
unbiased scores ensures convergence of the SGA iterates (<ref>) to a local maxima of the likelihood $p_{\theta}(y_{1:T})$
(see e.g. Kushner and Yin, 2003), and convergence of the SGLD and its ergodic averages
to the posterior distribution and posterior expectations, respectively [Teh et al., 2016].
In contrast, the use of biased gradient estimates
leads to asymptotic biases that have to be studied [Tadić and Doucet, 2017].
Moreover, these unbiased estimators can also attain better convergence rates
than Monte Carlo approaches based on the finest discretization level,
which can be made the same as MLMC under simple modifications [Vihola, 2018].

When employing such debiasing methods, we shall assume access to unbiased estimates
of differences of the score function associated to successive levels of discretization, i.e. $S_l(\theta) - S_{l-1}(\theta)$
where $S_l(\theta)$ denotes an approximation of $S(\theta)$ at discretization level $l\in\mathbb{N}$.
As $S_l(\theta)$ can be  written as an expectation w.r.t. a finite-dimensional smoothing distribution $\pi_{\theta}^l$,
biased but consistent approximations can be obtained using various particle smoothing algorithms
[Briers et al., 2010, Fearnhead et al., 2010, Douc et al., 2011] in the limit of the Monte Carlo sample size.
We will focus on the conditional particle filter (CPF) of Andrieu et al., 2010 which constructs
a Markov chain Monte Carlo (MCMC) algorithm targeting $\pi_{\theta}^l$ and produces biased
but consistent approximations of $S_l(\theta)$ in the limit of the number of MCMC iterations.
By employing recent advances in Jacob et al., 2020
that builds on earlier work by Glynn and Rhee, 2014, one can remove the MCMC burn-in bias and
estimate $S_l(\theta)$ unbiasedly by simulating a pair of CPF chains
that are coupled in such a way that they meet exactly after some random number of iterations.
For the debiasing schemes of Mcleish, 2011, Rhee and Glynn, 2015 to return estimators of $S(\theta)$
with finite variance, the variance of increments $S_l(\theta) - S_{l-1}(\theta)$ has to decrease sufficiently fast to zero with
the discretization level $l$.
Therefore, the scores in each increment have to be estimated in a dependent manner,
using appropriate couplings of pairs of CPF chains between successive discretization levels.
The construction of such couplings is our main methodological contribution.
This involves couplings of time-discretizations of the diffusion process
and couplings of resampling steps that are necessary to prevent weight degeneracy within the CPFs.
While the former can be achieved using common Brownian increments which is standard in MLMC
[Giles, 2008], the latter requires novel schemes to induce adequate dependencies
between multiple CPF chains as simple strategies based on common random variables will be inadequate.

Our main theoretical result is to establish, under assumptions, that the proposed methodology provides
unbiased estimators of the score function $S(\theta)$ with finite variance.

This article is structured as follows.
We introduce a representation of the score function as a smoothing expectation
in Section <ref> and its discrete-time approximation in Section <ref>.
Our proposed methodology to unbiasedly estimate the score function is presented
in Section <ref>. We establish various properties of our score estimators
under assumptions in Section <ref> and illustrate various aspects of our methodology
on three applications in Section <ref>.
The proof of our results are detailed in the appendix.
An R package to reproduce all numerical results can be found at

§ SCORE FUNCTIONS

§.§ Continuous-time representation

We begin by seeking a representation of the score function of the form in (<ref>).
For notational simplicity, we define $\Sigma(x)=\sigma(x)\sigma(x)^*$ and
$b_\theta(x) = \Sigma(x)^{-1}\sigma(x)^*a_{\theta}(x)$ for $x\in\mathbb{R}^d$, where $\sigma(x)^*$
denotes the transpose of $\sigma(x)$.
As the expectation in (<ref>) is w.r.t. $\mathbb{P}_{\theta}$, which
depends on the parameter $\theta$, we consider a change of measure to the law $\mathbb{Q}$ induced by
$dX_t = \sigma(X_t)dW_t$ with $X_0=x_{\star}$ on the space $\mathcal{C}_d([0,T])$.
Since $\mathbb{P}_{\theta}$ and $\mathbb{Q}$ are equivalent, by Girsanov theorem

\begin{align}\label{eqn:likelihood_girsanov}
p_{\theta}(y_{1:T}) = \mathbb{E}_{\mathbb{Q}}\left[\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(X)\prod_{t=1}^T g_{\theta}(y_{t}|X_{t})\right]
\end{align}

and the corresponding Radon–Nikodym derivative is

\begin{align}
\frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(X) =
\exp\left\lbrace -\frac{1}{2}\int_{0}^T \|b_{\theta}(X_t)\|_2^2dt +
\int_{0}^{T}b_{\theta}(X_t)^*\Sigma(X_t)^{-1}\sigma(X_t)^* dX_t\right\rbrace,
\end{align}

where $\|x\|_p$ denotes the $\mathbb{L}_p$-norm of a vector $x\in\mathbb{R}^d$.
We will assume throughout the article that $\theta\mapsto a_{\theta}(x)$ and $\theta\mapsto g_{\theta}(y|x)$ are
differentiable w.r.t. $\theta$ for each $(x,y)\in\mathbb{R}^d\times\mathbb{R}^{d_y}$.
Under mild regularity conditions, we may differentiate (<ref>) and represent
the score function as a smoothing expectation (<ref>) with the functional

\begin{align}\label{eqn:smoothing_functional}
G_{\theta}(X) = \nabla_{\theta} \log \frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(X) +
\sum_{t=1}^T\nabla_{\theta} \log g_{\theta}(y_{t}|X_{t}).
\end{align}

The first gradient term can be written as

\begin{align}
\nabla_{\theta} \log \frac{d\mathbb{P}_{\theta}}{d\mathbb{Q}}(X)
=-\int_{0}^{T}\left\{ \nabla_{\theta}a_{\theta}(X_{t})\right\} ^{*}\Sigma(X_{t})^{-1}a_{\theta}(X_{t})dt
+ \int_{0}^{T}\left\{ \nabla_{\theta}a_{\theta}(X_{t})\right\} ^{*}\Sigma(X_{t})^{-1}dX_{t},
\end{align}

where $\nabla_{\theta}a_{\theta}(x)\in\mathbb{R}^{d\times d_{\theta}}$ denotes the
Jacobian matrix w.r.t. $\theta$.

We can consider some generalization of the above representation.
Firstly, although we have assumed a deterministic initial condition of $X_0=x_{\star}$ to simplify our presentation,
random initializations can also be accommodated.
Assuming that $X_0\sim\mu_{\theta}$ is initialized from a distribution $\mu_{\theta}$
on $\mathbb{R}^d$ that admits a differentiable density w.r.t. the $d$-dimensional Lebesgue measure,
a term of $\nabla_{\theta}\log\mu_{\theta}(X_0)$ should be added to the expression in (<ref>).
This follows by conditioning on the value of $X_0$ and applying the above arguments.
Secondly, our methodology can also easily handle continuous-time observation models.
In this case, the term $\sum_{t=1}^T\nabla_{\theta} \log g_{\theta}(y_{t}|X_{t})$ in (<ref>)
would be replaced by the gradient of the log-conditional likelihood $\nabla_{\theta}\log p_{\theta}(y_{1:T}|X)$; see Section <ref> for an application where the observational model is given by  an inhomogenous Poisson point process.

When the diffusion coefficient also depends on unknown parameters, it may be possible to find an
invertible transformation $\Psi:\mathbb{R}^d\rightarrow\mathbb{R}^d$, such that the transformed process
$X_t = \Psi(Z_t)$ satisfies an SDE with a diffusion coefficient that is not parameter dependent.
We will illustrate two examples of this principle in Section <ref> using the Lamperti transformation
and a simple rescaling of each component of the diffusion process.
In more general cases where such transformations are not available, one has to seek a different representation
of the score function as the above change of measure is no longer applicable.
In Section <ref>, we consider such a representation and discuss how our proposed methodology
can be adapted.

§.§ Discrete-time approximation

As alluded to in the introduction, we will rely on time-discretizations of the diffusion process (<ref>).
We will employ a hierarchy of discretizations of the time interval $[0,T]$, indexed by a level parameter $l\in\mathbb{N}_0$ that determines the temporal resolution. Higher levels with finer time resolutions will require increased algorithmic cost.
For each level $l$, let $0= s_{0} < \cdots < s_{K_{l}} = T$ denote a dyadic uniform discretization of $[0,T]$,
defined as $s_k=k \Delta_l$ for $k\in\{0,1,\ldots,K_l\}$, where $\Delta_l=2^{-l}$ is the step-size and $K_l=2^l T$ is the
number of time steps.
Note that, by construction, $(s_k)_{k=0}^{K_l}$ contains the unit observation times $(t_p)_{p=1}^P$.

We consider the Euler–Maruyama scheme which defines a time-discretized process $X_{0:T}=(X_{s_k})_{k=0}^{K_l}$ according
to the following recursion for $k\in\{1,\ldots,K_{l}\}$:

\begin{align}\label{eq:disc_state}
X_{s_k} = X_{s_{k-1}} + a_{\theta}(X_{s_{k-1}}) \Delta_{l} + \sigma(X_{s_{k-1}})(W_{s_k} - W_{s_{k-1}}),
\quad X_0 = x_{\star}.
\end{align}

To simplify our notation, we omit the dependence of the time grid and
time-discretized process on the level parameter (until this distinction is necessary in Section <ref>),
and write the transition in (<ref>) as $X_{s_k}=F_{\theta}^l(X_{s_{k-1}},V_{s_k})$, where
$V_{s_k}=W_{s_k} - W_{s_{k-1}}$ denotes the Brownian increment.
Higher-order schemes such as the Milstein method could be employed but would be
difficult to implement for problems with dimension $d\geq 2$. Future work could
consider the antithetic truncated Milstein method of [Giles and Szpruch, 2014] for such settings.

In the following, we write $\mathcal{N}_d(\mu,\Sigma)$ to denote a $d$-dimensional
Gaussian distribution with mean $\mu\in\mathbb{R}^d$ and covariance $\Sigma\in\mathbb{R}^{d\times d}$, and
its density as $x\mapsto\mathcal{N}_d(x;\mu,\Sigma)$.
The notation $\mathcal{N}_d(0_d,I_d)$ refers to the standard Gaussian distribution, i.e.
with zero mean vector $0_d\in\mathbb{R}^d$ and identity covariance matrix $I_d\in\mathbb{R}^{d\times d}$.
Under time-discretization, the marginal likelihood of the resulting state space model is

\begin{align}
p_{\theta}^l(y_{1:T}) = \int_{\mathsf{X}^l}
\prod_{t=1}^T g_{\theta}(y_{t}|x_{t}) \delta_{x_{\star}}(dx_0)\prod_{k=1}^{K_l}f_{\theta}^l(dx_{s_k}|x_{s_{k-1}}),
\end{align}

where $\mathsf{X}^l=(\mathbb{R}^d)^{K_l+1}$ is the path space,
$\delta_{x_{\star}}(dx_0)$ refers to the Dirac measure at the deterministic initial condition $x_{\star}$,
and $f_{\theta}^l(dx_{s_k}|x_{s_{k-1}})=\mathcal{N}_d(x_{s_k}; x_{s_{k-1}} + a_{\theta}(x_{s_{k-1}}) \Delta_{l},\Delta_{l}\Sigma(x_{s_{k-1}}))dx_{s_k}$ denotes the Gaussian transition kernel corresponding to (<ref>).
Using Fisher's identity <cit.>, the score function
$S_l(\theta) = \nabla_{\theta}\log p_{\theta}^l(y_{1:T})$ at discretization level $l$
can be written as

\begin{align}\label{eqn:discretized_score}
S_l(\theta) = \int_{\mathsf{X}^l} G_{\theta}^l(x_{0:T}) \, \pi_{\theta}^l(dx_{0:T}),
\end{align}

where $G_{\theta}^l : \mathsf{X}^l \rightarrow \mathbb{R}^{d_{\theta}}$, defined as

\begin{align}\label{eqn:discrete_test_function}
= &-\sum_{k=1}^{K_l} \left\{ \nabla_{\theta}a_{\theta}(X_{s_{k-1}})\right\}^{*}\Sigma(X_{s_{k-1}})^{-1}
a_{\theta}(X_{s_{k-1}})\Delta_{l} \notag\\
& + \sum_{k=1}^{K_l} \left\{ \nabla_{\theta}a_{\theta}(X_{s_{k-1}})\right\}^{*}\Sigma(X_{s_{k-1}})^{-1}
(X_{s_{k}} - X_{s_{k-1}}) +
\sum_{t=1}^T\nabla_{\theta}\log g_{\theta}(y_{t}|X_{t})
\end{align}

can be seen as an approximation of $G_{\theta}(X)$ in (<ref>), and

\begin{align}\label{eqn:smoothing_distribution}
\pi_{\theta}^l(dx_{0:T}) = p_{\theta}^l(y_{1:T})^{-1}
\prod_{t=1}^T g_{\theta}(y_{t}|x_{t}) \delta_{x_{\star}}(dx_0)\prod_{k=1}^{K_l}f_{\theta}^l(dx_{s_k}|x_{s_{k-1}})
\end{align}

is the smoothing distribution on $\mathsf{X}^l$.

In Section <ref>, under appropriate regularity conditions, we will study the rate at which the discrete-time
approximation $S_l(\theta)$ converges to the desired score $S(\theta)$ as the level $l$ tends to infinity.
The following section concerns numerical approximations of these score functions.

§ SCORE ESTIMATION

§.§ Unbiased estimation strategy

The goal of this section is to construct unbiased estimators of the score function $S(\theta)$ that admits the representation
in (<ref>) and (<ref>). For each $\theta\in\Theta$, convergence of the
time-discretized score function (<ref>) as $l\rightarrow\infty$ allows us to write

\begin{align}\label{eqn:score_as_limit}
S(\theta) = \lim_{L\rightarrow\infty}\sum_{l=0}^L I_l(\theta),
\end{align}

where $I_l(\theta) = S_l(\theta) - S_{l-1}(\theta)$ denotes the score increment at level $l\in\mathbb{N}_0$ (with $S_{-1}(\theta)=0$).
We will design algorithms that allow us to construct unbiased estimators
$\widehat{I}_l(\theta)$ of $I_l(\theta)$ independently for any $l\in\mathbb{N}_0$, i.e.

\begin{align}\label{eq:ub1}
\mathbb{E}\left[\widehat{I}_l(\theta)\right] = I_l(\theta),\quad l\in\mathbb{N}_0,
\end{align}

where $\mathbb{E}$ denotes expectation w.r.t. the law of our algorithms.
The key insight of the debiasing schemes proposed by Mcleish, 2011, Rhee and Glynn, 2015 is to
perform a random truncation of the highest level $L$ in (<ref>), so that unbiasedness
is retained when unbiased estimators of the increments are employed.

Let $(P_l)_{l=0}^{\infty}$ be a given probability mass function (PMF) with support on $\mathbb{N}_0$ and
define its cumulative tail probabilities as $\mathcal{P}_l = \sum_{k=l}^{\infty} P_k$ for $l\in\mathbb{N}_0$.
Suppose that

\begin{align}\label{eq:ub2}
\sum_{l=0}^\infty \mathcal{P}_l^{-1}\left\lbrace \mathrm{Var}\left[\widehat{I}_l(\theta)^j\right] +
\left(S_{l}(\theta)^j-S(\theta)^j\right)^2\right\rbrace < \infty
\end{align}

for all $j\in\{1,\ldots,d_{\theta}\}$, where $\mathrm{Var}$ denotes variance under $\mathbb{E}$ and
$x^j$ refers to the $j^{th}$-component of the vector $x\in\mathbb{R}^{d_{\theta}}$.
If we sample $L$ from $(P_l)_{l=0}^{\infty}$ independently of $(\widehat{I}_l(\theta))_{l=0}^{\infty}$,
it follows from Rhee and Glynn, 2015 that the independent sum estimator

\begin{align}\label{eq:ub3}
\widehat{S}(\theta) = \sum_{l=0}^L \frac{\widehat{I}_l(\theta)}{\mathcal{P}_l}
\end{align}

is an unbiased estimator of $S(\theta)$ with finite variance, i.e. $\mathbb{E}[\widehat{S}(\theta)]=S(\theta)$
and the entries of the covariance matrix $\mathrm{Var}[\widehat{S}(\theta)]$ are finite.
Inspection of (<ref>) reveals that we have to understand how fast $S_l(\theta)$ converges to $S(\theta)$
as $l\rightarrow\infty$, which is established in Theorem <ref>.
Moreover, for (<ref>) and (<ref>) to hold, we have to
compute an unbiased estimator of the score $S_0(\theta)$ at level $l=0$ with finite variance, and
construct unbiased estimators of the score increments $I_l(\theta)$, whose variance vanishes
sufficiently fast as $l\rightarrow\infty$ relative to the tails of $(P_l)_{l=0}^{\infty}$.
Developing a methodology that meets these two requirements will be the focus of
Sections <ref> and <ref> respectively.
Verifying that our proposed methodology satisfies these requirements constitutes our main theoretical result in Theorem <ref>.
Alternatives to the estimator in (<ref>) such as the single term estimator
of Rhee and Glynn, 2015 could also be considered here.

§.§ Unbiased estimation of time-discretized scores

This section considers unbiased estimation of the score $S_l(\theta)$ at discretization level $l\in\mathbb{N}_0$;
the case $l=0$ will be employed to construct the first summand $\widehat{I}_0(\theta)$ in (<ref>).
Our basic algorithmic building block is the CPF of Andrieu et al., 2010, which defines a Markov kernel
on the space of trajectories $X_{0:T}\in\mathsf{X}^l$ that admits
the smoothing distribution $\pi_{\theta}^l$ as its invariant distribution.
A detailed description for our application is given in Algorithm <ref>,
which has a complexity of $\mathcal{O}(NK_l)$.
The CPF involves simulating $N\geq 2$ trajectories under the time-discretized model dynamics (Steps 2a & 2b),
weighting samples according to the observation density (Step 2c), and resampling from the weighted particle approximation (Step 2d).
We will consider multinomial resampling, in which case, $\mathcal{R}(w_t^{1:N})$ refers to the categorical distribution
on $\{1,\ldots,N\}$ with probabilities $w_t^{1:N}=(w_t^n)_{n=1}^N$.
The main difference to a standard bootstrap particle filter (BPF) [Gordon et al., 1993] is that the input trajectory
$X_{0:T}^{\star}$ is conditioned to survive all resampling steps.
The algorithm outputs a trajectory $X_{0:T}^{\circ}$ by sampling the particle indexes
$(B_{s_k})_{k=0}^{K_l}$ after the terminal step (Steps 3a & 3b).
We will write $X_{0:T}^{\circ}\sim M_{\theta}^l(\cdot |X_{0:T}^{\star})$
to denote a trajectory generated by the CPF kernel at parameter $\theta\in\Theta$ and
discretization level $l\in\mathbb{N}_0$.

We will initialize the CPF Markov chain using the law of the time-discretized dynamics

\begin{align}\label{eqn:law_discrete_dynamics}
\nu_{\theta}^l(dx_{0:T})=\delta_{x_{\star}}(dx_0)\prod_{k=1}^{K_l}f_{\theta}^l(dx_{s_k}|x_{s_{k-1}}).
\end{align}

One could also consider the law of a trajectory sampled from a BPF,
which provides a good approximation of $\pi_{\theta}^l$ for sufficiently large $N$ <cit.>.
Under mild assumptions, the Markov chain $(X_{0:T}(i))_{i=0}^{\infty}$ generated by

\begin{align}\label{eqn:cpf_chain}
X_{0:T}(0)\sim\nu_{\theta}^l,\quad X_{0:T}(i)\sim M_{\theta}^l(\cdot |X_{0:T}(i-1)),
\end{align}

for $i\geq1$, is uniformly ergodic [Chopin and Singh, 2015, Lindsten et al., 2015, Andrieu et al., 2018].
Hence one can adopt a standard MCMC approach to approximate the score $S_l(\theta)$ by the average

\begin{align}
A_{l}^{b:I}(\theta) = \frac{1}{I-b+1}\sum_{i=b}^IG_{\theta}^l(X_{0:T}(i)),
\end{align}

for some fixed burn-in $0\leq b\leq I$, which is consistent as the number of iterations $I\rightarrow\infty$.
However, as the Markov chain is not started at stationarity, the MCMC average $A_{l}^{b:I}(\theta)$ is a biased estimator
for any finite $I\in\mathbb{N}$. Although this burn-in bias can be reduced by increasing $b$,
tuning it to control the bias is a difficult task in practice.

By building on the work of Glynn and Rhee, 2014, Jacob et al., 2020 showed how to obtain
unbiased estimators of smoothing expectations by simulating a pair of coupled CPF chains
$(X_{0:T}(i),\bar{X}_{0:T}(i))_{i=0}^{\infty}$ on the product space $\mathsf{Z}^l=\mathsf{X}^l\times\mathsf{X}^l$
with the same marginal law.
This is achieved using a coupling of two CPFs as described in Algorithm <ref>, which we will
refer to as the 2-CCPF.
Writing $(X_{0:T}^{\circ},\bar{X}_{0:T}^{\circ})\sim \bar{M}_{\theta}^l(\cdot|X_{0:T}^{\star},\bar{X}_{0:T}^{\star})$
to denote a pair of trajectories generated by the 2-CCPF kernel given $(X_{0:T}^{\star},\bar{X}_{0:T}^{\star})\in\mathsf{Z}^l$ as input,
marginally we have $X_{0:T}^{\circ}\sim M_{\theta}^l(\cdot |X_{0:T}^{\star})$ and
$\bar{X}_{0:T}^{\circ}\sim M_{\theta}^l(\cdot |\bar{X}_{0:T}^{\star})$.
The two main ingredients of this coupling are the use of common Brownian increments
(Steps 2a & 2b), and a coupling of the resampling
distributions $\mathcal{R}(w_{t}^{1:N})$ and $\mathcal{R}(\bar{w}_{t}^{1:N})$
denoted as $\bar{\mathcal{R}}(w_{t}^{1:N},\bar{w}_{t}^{1:N})$ (Steps 2e & 3a).
We refer readers to the references in Jacob et al., 2020 for various coupled resampling schemes.
Our focus is the maximal coupling [Chopin and Singh, 2015, Jasra et al., 2017] that maximizes the
probability of having identical ancestors at each step of the CPFs.
This is detailed in Algorithm <ref>, where we have suppressed the time dependence for notational simplicity.
We could also consider an improvement of Step 5
that samples from the residual distributions with common uniform random variables.

As the cost of implementing Algorithm <ref> is $\mathcal{O}(N)$, the overall cost of Algorithm
<ref> is still $\mathcal{O}(NK_l)$.

We initialize by sampling $(X_{0:T}(0),\bar{X}_{0:T}(0))$ from a coupling $\bar{\nu}_{\theta}^l$ with
$\nu_{\theta}^l$ as its marginals.
The choice of $\bar{\nu}_{\theta}^l$ could be explored but we will consider the independent coupling for simplicity.
The pair of CPF chains is then generated as

\begin{align}\label{eqn:coupled2_chains}
X_{0:T}(1)\sim M_{\theta}^l(\cdot |X_{0:T}(0)),\quad
(X_{0:T}(i+1),\bar{X}_{0:T}(i))\sim  \bar{M}_{\theta}^l(\cdot | X_{0:T}(i),\bar{X}_{0:T}(i-1)),
\end{align}

for $i\geq 1$. Marginally, $(X_{0:T}(i))_{i=0}^{\infty}$ and $(\bar{X}_{0:T}(i))_{i=0}^{\infty}$ have the same
law as the Markov chain generated by (<ref>).
It can be shown that each application of 2-CCPF allows the chains to meet with some positive probability
<cit.>, that depends on the number of trajectories $N$ and observations $T$
<cit.>. Moreover, the construction of 2-CCPF ensures that the chains are faithful,
i.e. $X_{0:T}(i)=\bar{X}_{0:T}(i-1)$ for all $i\geq \tau_{\theta}^l$, where $\tau_{\theta}^l=\inf\{i\geq 1: X_{0:T}(i)=\bar{X}_{0:T}(i-1)\}$
denotes the meeting time. Using the time-averaged estimator of Jacob et al., 2020,
an unbiased estimator of the score $S_l(\theta)$ is given by

\begin{align}\label{eqn:unbiased_discretized_score}
\widehat{S}_l(\theta) = A_{l}^{b:I}(\theta) + \sum_{i=b+1}^{\tau_{\theta}^l-1}\min\left(1,\frac{i-b}{I-b+1}\right)\left(G_{\theta}^l(X_{0:T}(i))
- G_{\theta}^l(\bar{X}_{0:T}(i-1))\right).
\end{align}

The second term corrects for the bias of the MCMC average $A_{l}^{b:I}(\theta)$ and is equal to zero if $b+1 > \tau_{\theta}^l-1$.
Under assumptions that will be stated in Section <ref>, $\widehat{S}_l(\theta)$
has finite variance and finite expected cost, for any choice of $N\geq2$ and $0\leq b\leq I$.
Assuming that 2-CCPF costs twice as much as CPF, the cost of computing $\widehat{S}_l(\theta)$
is $\max\{2\tau_{\theta}^l-1,I+\tau_{\theta}^l-1\}$ applications of the CPF kernel $M_{\theta}^l$.

Conditional particle filter (CPF) at parameter $\theta\in\Theta$ and discretization level $l\in\mathbb{N}_0$

Input: a trajectory $X_{0:T}^{\star}=(X_{s_k}^{\star})_{k=0}^{K_l}\in\mathsf{X}^{l}$.

For time step $k=0$

* (1a) Set $X_{0}^{n}=x_{\star}$ for $n\in\{1,\ldots,N\}$.

* (1b) Set $w_0^{n}=N^{-1}$ and $A_{0}^{n}=n$ for $n\in\{1,\ldots,N\}$.

For time step $k\in\{1,\ldots,K_l\}$

* (2a) Sample Brownian increment $V_{s_k}^n \sim \mathcal{N}_d(0_d,\Delta_{l}I_d)$ independently for $n\in\{1,\ldots,N-1\}$.

* (2b) Set $X_{s_k}^n =F_{\theta}^l(X_{s_{k-1}}^{A_{s_{k-1}}^{n}}, V_{s_k}^n)$
for $n\in\{1,\ldots,N-1\}$, and $X_{s_k}^N = X_{s_k}^{\star}$.

* If there is an observation at time $t=s_k\in\{1,\ldots,T\}$

* (2c) Compute normalized weight $w_t^n \propto g_{\theta}(y_{t}|X_{t}^n)$ for $n\in\{1,\ldots,N\}$.

* (2d) If $t < T$, sample ancestor $A_{t}^{n}\sim \mathcal{R}(w_{t}^{1:N})$ independently
for $n\in\{1,\ldots,N-1\}$ and set $A_{t}^N=N$.

* Else

* (2e) Set $A_{s_k}^{n}=n$ for $n\in\{1,\ldots,N\}$.

After the terminal step

* (3a) Sample particle index $B_{T}\sim \mathcal{R}(w_{T}^{1:N})$.

* (3b) Set particle index $B_{s_k} = A_{s_k}^{B_{s_{k+1}}}$ for $k\in\{0,1,\ldots,K_l-1\}$.

Output: a trajectory $X_{0:T}^{\circ}=(X_{s_k}^{B_{s_k}})_{k=0}^{K_l}\in\mathsf{X}^{l}$.

Two coupled CPF (2-CCPF) at parameter $\theta\in\Theta$ and discretization level $l\in\mathbb{N}_0$

Input: a pair of trajectories $(X_{0:T}^{\star},\bar{X}_{0:T}^{\star})=(X_{s_k}^{\star},\bar{X}_{s_k}^{\star})_{k=0}^{K_l}\in\mathsf{Z}^{l}$.

For time step $k=0$

* (1a) Set $X_{0}^{n}=x_{\star}$ and $\bar{X}_{0}^{n}=x_{\star}$ for $n\in\{1,\ldots,N\}$.

* (1b) Set $w_0^{n}=N^{-1}, \bar{w}_0^{n}=N^{-1}$ and $A_{0}^{n}=n, \bar{A}_{0}^{n}=n$ for $n\in\{1,\ldots,N\}$.

For time step $k\in\{1,\ldots,K_l\}$

* (2a) Sample Brownian increment $V_{s_k}^n \sim \mathcal{N}_d(0_d,\Delta_{l}I_d)$ independently for $n\in\{1,\ldots,N-1\}$.

* (2b) Set $X_{s_k}^n = F_{\theta}^l(X_{s_{k-1}}^{A_{s_{k-1}}^{n}},V_{s_k}^n)$ and
$\bar{X}_{s_k}^n = F_{\theta}^l(\bar{X}_{s_{k-1}}^{\bar{A}_{s_{k-1}}^{n}},V_{s_k}^n)$
for $n\in\{1,\ldots,N-1\}$.

* (2c) Set $X_{s_k}^N = X_{s_k}^{\star}$ and $\bar{X}_{s_k}^N = \bar{X}_{s_k}^{\star}$.

* If there is an observation at time $t=s_k\in\{1,\ldots,T\}$

* (2d) Compute normalized weights
$w_t^n \propto g_{\theta}(y_{t}|X_{t}^n)$ and
$\bar{w}_t^n \propto g_{\theta}(y_{t}|\bar{X}_{t}^n)$ for $n\in\{1,\ldots,N\}$.

* (2e) If $t < T$, sample ancestors $(A_{t}^{n},\bar{A}_{t}^{n})\sim \bar{\mathcal{R}}(w_{t}^{1:N},\bar{w}_{t}^{1:N})$ independently
for $n\in\{1,\ldots,N-1\}$ and set $A_{t}^N=N,\bar{A}_{t}^N=N$.

* Else

* (2f) Set $A_{s_k}^{n}=n$ and $\bar{A}_{s_k}^{n}=n$ for $n\in\{1,\ldots,N\}$.

After the terminal step

* (3a) Sample particle indexes $(B_{T},\bar{B}_{T})\sim \bar{\mathcal{R}}(w_{T}^{1:N},\bar{w}_{T}^{1:N})$.

* (3b) Set particle indexes $B_{s_k} = A_{s_k}^{B_{s_{k+1}}}$ and $\bar{B}_{s_k} = \bar{A}_{s_k}^{\bar{B}_{s_{k+1}}}$ for $k\in\{0,1,\ldots,K_l-1\}$.

Output: a pair of trajectories $(X_{0:T}^{\circ},\bar{X}_{0:T}^{\circ})=(X_{s_k}^{B_{s_k}},\bar{X}_{s_k}^{\bar{B}_{s_k}})_{k=0}^{K_l}\in\mathsf{Z}^l$.

Maximal coupling of two resampling distributions $\mathcal{R}(w^{1:N})$ and $\mathcal{R}(\bar{w}^{1:N})$

Input: normalized weights $w^{1:N}=(w^n)_{n=1}^N$ and $\bar{w}^{1:N}=(\bar{w}^n)_{n=1}^N$.

(1) Compute the overlap $o^n=\min\{w^n,\bar{w}^n\}$ for $n\in\{1,\ldots,N$}.

(2) Compute the mass of the overlap $\mu=\sum_{n=1}^No^n$ and normalize $O^n=o^n/\mu$ for $n\in\{1,\ldots,N\}$.

(3) Compute the residuals $r^n=(w^n-o^n)/(1-\mu)$ and $\bar{r}^n=(\bar{w}^n-o^n)/(1-\mu)$ for $n\in\{1,\ldots,N\}$.

With probability $\mu$

* (4) Sample $A\sim\mathcal{R}(O^{1:N})$ and set $\bar{A}=A$.

* (5) Sample $A\sim\mathcal{R}(r^{1:N})$ and $\bar{A}\sim\mathcal{R}(\bar{r}^{1:N})$ independently.

Output: indexes $(A,\bar{A})$.

§.§ Unbiased estimation of score increments

We now consider unbiased estimation of the score increment $I_l(\theta) = S_l(\theta) - S_{l-1}(\theta)$
at level $l\in\mathbb{N}$ to construct the term $\widehat{I}_l(\theta)$
in the independent sum estimator $\widehat{S}(\theta)$ in (<ref>).
A naive approach that employs the unbiased estimation framework in Section <ref>
to estimate the scores $S_{l-1}(\theta)$ and $S_l(\theta)$ independently, and take the difference to estimate
the increment $I_l(\theta)$, will satisfy the unbiasedness requirement in (<ref>).
However, as the variance of the increment will not decrease with the discretization level $l$ under independent estimation
of the successive scores, one cannot choose a PMF $(P_l)_{l=0}^{\infty}$ such that
the condition in (<ref>) holds.
This prompts an extension of the preceding framework that allows the scores in each increment to be estimated
in a dependent manner.

Our proposed methodology involves simulating two pairs of coupled CPF chains,
$(X_{0:T}^{l-1}(i),\bar{X}_{0:T}^{l-1}(i))_{i=0}^{\infty}$ on $\mathsf{Z}^{l-1}$
for discretization level $l-1$,
and $(X_{0:T}^{l}(i),\bar{X}_{0:T}^{l}(i))_{i=0}^{\infty}$ on $\mathsf{Z}^{l}$
for discretization level $l$.
This relies on the coupling of four CPFs detailed in Algorithm <ref>,
which will be referred to as the 4-CCPF.
In the algorithmic description, $(s_k)_{k=0}^{K_l}$ denotes a time discretization of $[0,T]$ at level $l$,
and a trajectory at the coarser level $l-1$ is written as $X_{0:T}^{l-1}=(X_{s_{2k}}^{l-1})_{k=0}^{K_{l-1}}$.
Given trajectories $(X_{0:T}^{l-1,\star},\bar{X}_{0:T}^{l-1,\star},X_{0:T}^{l,\star},\bar{X}_{0:T}^{l,\star})\in
\mathsf{Z}^{l-1}\times\mathsf{Z}^{l}$ as input,
we will write

\begin{align}\label{eqn:4CCPF}
\bar{M}_{\theta}^{l-1,l}(\cdot|X_{0:T}^{l-1,\star},\bar{X}_{0:T}^{l-1,\star},X_{0:T}^{l,\star},\bar{X}_{0:T}^{l,\star})
\end{align}

to denote the trajectories generated by 4-CCPF.
The 4-CCPF kernel $\bar{M}_{\theta}^{l-1,l}$ is a four-marginal coupling of the CPF kernels $M_{\theta}^{l-1}$
and $M_{\theta}^{l}$, in the sense that marginally, we have $X_{0:T}^{l-1,\circ}\sim M_{\theta}^{l-1}(\cdot |X_{0:T}^{l-1,\star})$ and
$\bar{X}_{0:T}^{l-1,\circ}\sim M_{\theta}^{l-1}(\cdot |\bar{X}_{0:T}^{l-1,\star})$ at level $l-1$, and
$X_{0:T}^{l,\circ}\sim M_{\theta}^l(\cdot |X_{0:T}^{l,\star})$ and
$\bar{X}_{0:T}^{l,\circ}\sim M_{\theta}^l(\cdot |\bar{X}_{0:T}^{l,\star})$ at level $l$.
The two main ingredients of 4-CCPF are the use of common Brownian increments within each level
(Steps 2b & 2e) and across levels (Steps 2a & 2d),
and an appropriate four-marginal coupling of the resampling
$\mathcal{R}(w_{t}^{l-1,1:N})$, $\mathcal{R}(\bar{w}_{t}^{l-1,1:N})$,
$\mathcal{R}(w_{t}^{l,1:N})$, $\mathcal{R}(\bar{w}_{t}^{l,1:N})$
denoted by
$\bar{\mathcal{R}}(w_{t}^{l-1,1:N},\bar{w}_{t}^{l-1,1:N},w_{t}^{l,1:N},\bar{w}_{t}^{l,1:N})$ (Steps 2h & 3a).
While the use of common Brownian increments is a standard choice in MLMC
[Giles, 2008], constructing coupled resampling schemes that induce sufficient dependencies
between the four CPF chains, for the variance of the estimated increment to decrease with the discretization level,
requires new algorithmic design.

Algorithm <ref> presents a coupled resampling scheme that will be the focus of our analysis in Section <ref>.
Here we suppress the time dependence for notational simplicity.
Given normalized weights $w^{l-1,1:N}=(w^{l-1,n})_{n=1}^N$, $\bar{w}^{l-1,1:N}=(\bar{w}^{l-1,n})_{n=1}^N$
at discretization level $l-1$ and $w^{l,1:N}=(w^{l,n})_{n=1}^N$, $\bar{w}^{l,1:N}=(\bar{w}^{l,n})_{n=1}^N$ at discretization level $l$,
the algorithm samples ancestor indexes $(A^{l-1},\bar{A}^{l-1},A^{l},\bar{A}^{l})$ from the maximal coupling of the maximal couplings
$\bar{\mathcal{R}}(w^{l-1,1:N},w^{l,1:N})$ and $\bar{\mathcal{R}}(\bar{w}^{l-1,1:N},\bar{w}^{l,1:N})$.
That is, amongst all possible four-marginal couplings, this scheme maximizes the probabilities of having
identical ancestors across levels, i.e. $A^{l-1}=A^l$ and $\bar{A}^{l-1}=\bar{A}^l$,
and identical pair of ancestors within the levels, i.e. $(A^{l-1},A^l)=(\bar{A}^{l-1},\bar{A}^l)$.
The cost of Algorithm <ref> is random as it employs rejection samplers [Thorisson, 2000]
in Steps 2b, 3b and 4, wherein $\mathcal{U}_{[0,1]}$ denotes the uniform distribution on $[0,1]$.
As the expected cost is $\mathcal{O}(N)$, the overall cost of Algorithm <ref> is $\mathcal{O}(NK_l)$ on average.
Note that a naive approach to sample from the desired coupled resampling scheme in place of
Algorithm <ref> would involve a deterministic but prohibitive cost of $\mathcal{O}(N^2)$.
In Step 4, we denote the respective PMF of
$\bar{\mathcal{R}}(w^{l-1,1:N},w^{l,1:N})$ and $\bar{\mathcal{R}}(\bar{w}^{l-1,1:N},\bar{w}^{l,1:N})$ as

\begin{align}\label{eqn:PMF_maximal}
&R^{l-1,l}(A,B)=\mathbb{I}_{\mathcal{D}}(A,B)o^{A} +
\frac{(w^{l-1,A}-o^A)(w^{l,B}-o^B)}{1-\mu}, \\
&\bar{R}^{l-1,l}(A,B)=\mathbb{I}_{\mathcal{D}}(A,B)\bar{o}^{A} +
\frac{(\bar{w}^{l-1,A}-\bar{o}^A)(\bar{w}^{l,B}-\bar{o}^B)}{1-\bar{\mu}}, \notag
\end{align}

for $(A,B)\in\{1,\ldots,N\}^2$, where $\mathbb{I}_{\mathcal{D}}(A,B)$ is the indicator function on the diagonal set
$\mathcal{D}=\{(A,B)\in\{1,\ldots,N\}^2 : A=B\}$,
$o^n=\min\{w^{l-1,n},w^{l,n}\}$ and $\bar{o}^n=\min\{\bar{w}^{l-1,n},\bar{w}^{l,n}\}$ for $n\in\{1,\ldots,N\}$
are the overlapping measures, and
$\mu=\sum_{n=1}^No^n$ and $\bar{\mu}=\sum_{n=1}^N\bar{o}^n$
are their corresponding mass.
From the expressions in (<ref>), one can check that the three cases considered in Algorithm
<ref> are necessary to ensure faithfulness of the pair of chains on each discretization level.
More precisely, if the input trajectories satisfy $X_{0:T}^{l-1,\star}=\bar{X}_{0:T}^{l-1,\star}$
and/or $X_{0:T}^{l,\star}=\bar{X}_{0:T}^{l,\star}$, then under the 4-CCPF (<ref>),
the output trajectories satisfy
$X_{0:T}^{l-1,\circ}=\bar{X}_{0:T}^{l-1,\circ}$ and/or $X_{0:T}^{l,\circ}=\bar{X}_{0:T}^{l,\circ}$
almost surely.

We note that the two-marginal couplings induced by the 4-CCPF kernel $\bar{M}_{\theta}^{l-1,l}$
on each discretization level are not the same as the 2-CCPF kernels $\bar{M}_{\theta}^{l-1}$ and $\bar{M}_{\theta}^l$.
Although it is not a requirement of our methodology,
this property would hold if we consider a modification of Algorithm <ref> that samples from
the maximal coupling of the maximal couplings
$\bar{\mathcal{R}}(w^{l-1,1:N},\bar{w}^{l-1,1:N})$ and $\bar{\mathcal{R}}(w^{l,1:N},\bar{w}^{l,1:N})$.
However, as with simple coupling strategies based on common random variables,
such a coupled resampling scheme does not induce adequate dependencies between the CPF chains across discretization levels.
This will be illustrated experimentally in Section <ref>.
To understand the rationale behind Algorithm <ref>, we observe that the
two-marginal couplings induced by the 4-CCPF kernel across discretization levels
are given by the multilevel CPF (ML-CPF) described in Algorithm <ref>.
That is, writing $M_{\theta}^{l-1,l}$ as the ML-CPF kernel,
which is a coupling of the CPF kernels $M_{\theta}^{l-1}$ and $M_{\theta}^{l}$, we have
$(X_{0:T}^{l-1,\circ},X_{0:T}^{l,\circ})\sim M_{\theta}^{l-1,l}(\cdot|X_{0:T}^{l-1,\star},X_{0:T}^{l,\star})$ and
$(\bar{X}_{0:T}^{l-1,\circ},\bar{X}_{0:T}^{l,\circ})\sim M_{\theta}^{l-1,l}(\cdot|\bar{X}_{0:T}^{l-1,\star},\bar{X}_{0:T}^{l,\star})$
under the 4-CCPF kernel in (<ref>).
The ML-CPF is similar to the multilevel particle filter of Jasra et al., 2017,
who proposed multilevel estimators of filtering expectations that are non-asymptotically biased
but consistent in the limit of our computational budget.
Even though our objective is markedly different, as we seek non-asymptotically unbiased and
(almost surely) finite cost estimators of the score function which is a smoothing expectation,
the connection to multilevel estimation alludes to better convergence rates than Monte Carlo
approaches based on the finest discretization level.

We now describe simulation of the two pairs of CPF chains
$(X_{0:T}^{l-1}(i),\bar{X}_{0:T}^{l-1}(i))_{i=0}^{\infty}$ and
$(X_{0:T}^{l}(i),\bar{X}_{0:T}^{l}(i))_{i=0}^{\infty}$ using ML-CPF and 4-CCPF.
We initialize $(X_{0:T}^{l-1}(0),\bar{X}_{0:T}^{l-1}(0),X_{0:T}^{l}(0),\bar{X}_{0:T}^l(0))$ from
a four-marginal coupling $\bar{\nu}_{\theta}^{l-1,l}$ that satisfies
$X_{0:T}^{l-1}(0),\bar{X}_{0:T}^{l-1}(0)\sim\nu_{\theta}^{l-1}$ and $X_{0:T}^{l}(0),\bar{X}_{0:T}^l(0)\sim\nu_{\theta}^{l}$.
For simplicity, we assume $\bar{\nu}_{\theta}^{l-1,l}$ is such that
each of the pairs across levels $(X_{0:T}^{l-1}(0),X_{0:T}^{l}(0))$ and $(\bar{X}_{0:T}^{l-1}(0),\bar{X}_{0:T}^{l}(0))$
independently follow a coupling of $\nu_{\theta}^{l-1}$ and $\nu_{\theta}^{l}$, denoted as $\nu_{\theta}^{l-1,l}$.
The choice of $\nu_{\theta}^{l-1,l}$ will be taken as the joint law of the time-discretized dynamics
under common Brownian increments in our analysis.

We then sample $(X_{0:T}^{l-1}(1),X_{0:T}^l(1))\sim M_{\theta}^{l-1,l}(\cdot | X_{0:T}^{l-1}(0),X_{0:T}^l(0))$ with ML-CPF,
and subsequently for $i\geq 1$, iteratively sample

\begin{align}\label{eqn:coupled4_chains}
\bar{M}_{\theta}^{l-1,l}(\cdot|X_{0:T}^{l-1}(i),\bar{X}_{0:T}^{l-1}(i-1),X_{0:T}^{l}(i),\bar{X}_{0:T}^{l}(i-1))
\end{align}

from 4-CCPF. Marginally, the single CPF chains have the same law as a Markov chain generated by (<ref>) at discretization level $l-1$ or $l$.
Since each application of 4-CCPF allows the pair of chains on each level to meet with some positive probability (see Lemma <ref> in  Appendix <ref>),
and by construction remain faithful thereafter, we define
the meeting time at level $l$ as $\tau_{\theta}^l=\inf\{i\geq 1: X_{0:T}^l(i)=\bar{X}_{0:T}^l(i-1)\}$
and the stopping time at level $l$ as $\bar{\tau}_{\theta}^l=\max\{\tau_{\theta}^{l-1},\tau_{\theta}^l\}$.
Note that the 4-CCPF collapses to the ML-CPF after the stopping time, i.e.
for $i>\bar{\tau}_{\theta}^l$, the transition in (<ref>) is equivalent to sampling
$(X_{0:T}^{l-1}(i+1),X_{0:T}^{l}(i+1))\sim M_{\theta}^{l-1,l}(\cdot|X_{0:T}^{l-1}(i),X_{0:T}^{l}(i))$ and
setting $\bar{X}_{0:T}^{l-1}(i)=X_{0:T}^{l-1}(i+1), \bar{X}_{0:T}^{l}(i)=X_{0:T}^{l}(i+1)$.
For any choice of burn-in $b\in\mathbb{N}_0$ and number of iterations $I\geq b$, we can compute unbiased estimators
$\widehat{S}_{l-1}(\theta)$ and $\widehat{S}_l(\theta)$
of the scores $S_{l-1}(\theta)$ and $S_l(\theta)$
using the time-averaged estimator in (<ref>) based on the pair of CPF chains
on level $l-1$ and $l$, respectively.
We can then obtain an unbiased estimator of the score increment $I_l(\theta)$ using the difference
$\widehat{I}_l(\theta)=\widehat{S}_l(\theta) - \widehat{S}_{l-1}(\theta)$,
which has finite variance and finite expected cost under the assumptions in Section <ref>.
The cost of computing $\widehat{I}_l(\theta)$ is[This assumes that
$\mathrm{Cost}(M_{\theta}^{l-1,l})=\mathrm{Cost}(M_{\theta}^{l-1}) + \mathrm{Cost}(M_{\theta}^{l})$,
$\mathrm{Cost}(\bar{M}_{\theta}^{l-1,l})=2\mathrm{Cost}(M_{\theta}^{l-1}) + 2\mathrm{Cost}(M_{\theta}^{l})$
and $\mathrm{Cost}(M_{\theta}^{l})=2\mathrm{Cost}(M_{\theta}^{l-1})$.]
$\max\{2\tau_{\theta}^{l-1}-1,I+\tau_{\theta}^{l-1}-1\}/2 + \max\{2\tau_{\theta}^l-1,I+\tau_{\theta}^l-1\}$ applications of the CPF kernel $M_{\theta}^l$.

Four coupled CPF (4-CCPF) at parameter $\theta\in\Theta$ and discretization levels $l-1$ and $l\in\mathbb{N}$.

Input: a pair of trajectories $(X_{0:T}^{l-1,\star},\bar{X}_{0:T}^{l-1,\star})=
and a pair of trajectories $(X_{0:T}^{l,\star},\bar{X}_{0:T}^{l,\star})=(X_{s_k}^{l,\star},\bar{X}_{s_k}^{l,\star})_{k=0}^{K_{l}}\in\mathsf{Z}^{l}$.

For time step $k=0$

* (1a) Set $X_{0}^{l-1,n}=x_{\star},\bar{X}_{0}^{l-1,n}=x_{\star}$ and
$X_{0}^{l,n}=x_{\star},\bar{X}_{0}^{l,n}=x_{\star}$ for $n\in\{1,\ldots,N\}$.

* (1b) Set $w_0^{l-1,n}=N^{-1},\bar{w}_0^{l-1,n}=N^{-1}$ and
$w_0^{l,n}=N^{-1},\bar{w}_0^{l,n}=N^{-1}$ for $n\in\{1,\ldots,N\}$.

* (1c) Set $A_{0}^{l-1,n}=n, \bar{A}_{0}^{l-1,n}=n$ and
$A_{0}^{l,n}=n,\bar{A}_{0}^{l,n}=n$ for $n\in\{1,\ldots,N\}$.

For time step $k\in\{1,\ldots,K_l\}$

* (2a) Sample Brownian increment $V_{s_k}^{l,n} \sim \mathcal{N}_d(0_d,\Delta_{l}I_d)$ at level $l$ independently for $n\in\{1,\ldots,N-1\}$.

* (2b) Set $X_{s_k}^{l,n} = F_{\theta}^l(X_{s_{k-1}}^{l,A_{s_{k-1}}^{l,n}},V_{s_k}^{l,n})$ and
$\bar{X}_{s_k}^{l,n} = F_{\theta}^l(\bar{X}_{s_{k-1}}^{l,\bar{A}_{s_{k-1}}^{l,n}}, V_{s_k}^{l,n})$
at level $l$ for $n\in\{1,\ldots,N-1\}$.

* (2c) Set $X_{s_k}^{l,N} = X_{s_k}^{l,\star}$ and $\bar{X}_{s_k}^{l,N} = \bar{X}_{s_k}^{l,\star}$.

* If $k\in\{2, 4, \ldots, K_l\}$

* (2d) Set Brownian increment $V_{s_k}^{l-1,n} =  V_{s_{k-1}}^{l,n} + V_{s_k}^{l,n}$ at level $l-1$ for $n\in\{1,\ldots,N-1\}$.

* (2e) Set $X_{s_k}^{l-1,n} = F_{\theta}^{l-1}(X_{s_{k-1}}^{l-1,A_{s_{k-1}}^{l-1,n}}, V_{s_k}^{l-1,n})$ and
$\bar{X}_{s_k}^{l-1,n} = F_{\theta}^{l-1}(\bar{X}_{s_{k-1}}^{l-1,\bar{A}_{s_{k-1}}^{l-1,n}}, V_{s_k}^{l-1,n})$
at level $l-1$ for $n\in\{1,\ldots,N-1\}$.

* (2f) Set $X_{s_k}^{l-1,N} = X_{s_k}^{l-1,\star}$ and $\bar{X}_{s_k}^{l-1,N} = \bar{X}_{s_k}^{l-1,\star}$.

* If there is an observation at time $t=s_k\in\{1,\ldots,T\}$

* (2g) Compute normalized weights
$w_t^{l-1,n} \propto g_{\theta}(y_{t}|X_{t}^{l-1,n}),
\bar{w}_t^{l-1,n} \propto g_{\theta}(y_{t}|\bar{X}_{t}^{l-1,n})$ and
$w_t^{l,n} \propto g_{\theta}(y_{t}|X_{t}^{l,n}),
\bar{w}_t^{l,n} \propto g_{\theta}(y_{t}|\bar{X}_{t}^{l,n})$ for $n\in\{1,\ldots,N\}$.

* (2h) If $t < T$, sample ancestors $(A_{t}^{l-1,n},\bar{A}_{t}^{l-1,n},A_{t}^{l,n},\bar{A}_{t}^{l,n})
\sim\bar{\mathcal{R}}(w_{t}^{l-1,1:N},\bar{w}_{t}^{l-1,1:N},w_{t}^{l,1:N},\bar{w}_{t}^{l,1:N})$ independently
for $n\in\{1,\ldots,N-1\}$ and set $A_{t}^{l-1,N}=N,\bar{A}_{t}^{l-1,N}=N,A_{t}^{l,N}=N,\bar{A}_{t}^{l,N}=N$.

* Else

* (2i) If $k\in\{2, 4, \ldots, K_l\}$, set $A_{s_k}^{l-1,n}=n$ and $\bar{A}_{s_k}^{l-1,n}=n$ at level $l-1$ for $n\in\{1,\ldots,N\}$.

* (2j) Set $A_{s_k}^{l,n}=n$ and $\bar{A}_{s_k}^{l,n}=n$ at level $l$ for $n\in\{1,\ldots,N\}$.

After the terminal step

* (3a) Sample particle indexes $(B_{T}^{l-1},\bar{B}_{T}^{l-1},B_{T}^{l},\bar{B}_{T}^{l})\sim
\bar{\mathcal{R}}(w_{T}^{l-1,1:N},\bar{w}_{T}^{l-1,1:N},

* (3b) Set particle indexes $B_{s_{2k}}^{l-1} = A_{s_{2k}}^{l-1,B_{s_{2(k+1)}}^{l-1}}$ and
$\bar{B}_{s_{2k}}^{l-1} = \bar{A}_{s_{2k}}^{l-1,\bar{B}_{s_{2(k+1)}}^{l-1}}$
at level $l-1$ for $k\in\{0,1,\ldots,K_{l-1}-1\}$.

* (3c) Set particle indexes $B_{s_k}^l = A_{s_k}^{l,B_{s_{k+1}}^l}$ and $\bar{B}_{s_k}^l = \bar{A}_{s_k}^{l,\bar{B}_{s_{k+1}}^l}$ at level $l$ for $k\in\{0,1,\ldots,K_l-1\}$.

Output: a pair of trajectories $(X_{0:T}^{l-1,\circ},\bar{X}_{0:T}^{l-1,\circ})=
and a pair of trajectories $(X_{0:T}^{l,\circ},\bar{X}_{0:T}^{l,\circ})=(X_{s_k}^{l,B_{s_k}^l},\bar{X}_{s_k}^{l,\bar{B}_{s_k}^l})_{k=0}^{K_l}\in\mathsf{Z}^{l}$.

Maximal coupling of the maximal couplings $\bar{\mathcal{R}}(w^{l-1,1:N},w^{l,1:N})$
and $\bar{\mathcal{R}}(\bar{w}^{l-1,1:N},\bar{w}^{l,1:N})$

Input: normalized weights $w^{l-1,1:N}=(w^{l-1,n})_{n=1}^N$, $\bar{w}^{l-1,1:N}=(\bar{w}^{l-1,n})_{n=1}^N$ at level $l-1$ and
$w^{l,1:N}=(w^{l,n})_{n=1}^N$, $\bar{w}^{l,1:N}=(\bar{w}^{l,n})_{n=1}^N$ at level $l$.

(1) Sample $(A^{l-1},A^l)\sim \bar{\mathcal{R}}(w^{l-1,1:N},w^{l,1:N})$ using Algorithm <ref>.

If normalized weights at level $l-1$ are identical and normalized weights at level $l$ are non-identical

* (2a) Set $\bar{A}^{l-1}=A^{l-1}$.

* (2b) With probability $\bar{w}^{l,\bar{A}^{l-1}}/\bar{w}^{l-1,\bar{A}^{l-1}}$, set $\bar{A}^l=\bar{A}^{l-1}$;
otherwise sample $A\sim\mathcal{R}(\bar{w}^{l,1:N})$ and $U\sim\mathcal{U}_{[0,1]}$
until $U>\bar{w}^{l-1,A}/\bar{w}^{l,A}$, and set $\bar{A}^l=A$.

If normalized weights at level $l-1$ are non-identical and normalized weights at level $l$ are identical

* (3a) Set $\bar{A}^{l}=A^{l}$.

* (3b) With probability $\bar{w}^{l-1,\bar{A}^{l}}/\bar{w}^{l,\bar{A}^{l}}$, set $\bar{A}^{l-1}=\bar{A}^{l}$;
otherwise sample $A\sim\mathcal{R}(\bar{w}^{l-1,1:N})$ and $U\sim\mathcal{U}_{[0,1]}$
until $U>\bar{w}^{l,A}/\bar{w}^{l-1,A}$, and set $\bar{A}^{l-1}=A$.

* (4) With probability $\bar{R}^{l-1,l}(A^{l-1},A^l)/R^{l-1,l}(A^{l-1},A^l)$, set $(\bar{A}^{l-1},\bar{A}^l)=(A^{l-1},A^l)$;
otherwise sample $(A,B)\sim \bar{\mathcal{R}}(\bar{w}^{l-1,1:N},\bar{w}^{l,1:N})$ and $U\sim\mathcal{U}_{[0,1]}$
until $U>R^{l-1,l}(A,B)/\bar{R}^{l-1,l}(A,B)$, and set $(\bar{A}^{l-1},\bar{A}^l)=(A,B)$.

Output: indexes $(A^{l-1},\bar{A}^{l-1},A^{l},\bar{A}^{l})$.

Multilevel CPF (ML-CPF) at parameter $\theta\in\Theta$ and discretization levels $l-1$ and $l\in\mathbb{N}$.

Input: a trajectory $X_{0:T}^{l-1,\star}=
and a trajectory $X_{0:T}^{l,\star}=(X_{s_k}^{l,\star})_{k=0}^{K_{l}}\in\mathsf{X}^{l}$.

For time step $k=0$

* (1a) Set $X_{0}^{l-1,n}=x_{\star}$ and
$X_{0}^{l,n}=x_{\star}$ for $n\in\{1,\ldots,N\}$.

* (1b) Set $w_0^{l-1,n}=N^{-1}, w_0^{l,n}=N^{-1}$ and
$A_{0}^{l-1,n}=n, A_{0}^{l,n}=n$ for $n\in\{1,\ldots,N\}$.

For time step $k\in\{1,\ldots,K_l\}$

* (2a) Sample Brownian increment $V_{s_k}^{l,n} \sim \mathcal{N}_d(0_d,\Delta_{l}I_d)$ at level $l$ independently for $n\in\{1,\ldots,N-1\}$.

* (2b) Set $X_{s_k}^{l,n} = F_{\theta}^l(X_{s_{k-1}}^{l,A_{s_{k-1}}^{l,n}},V_{s_k}^{l,n})$
at level $l$ for $n\in\{1,\ldots,N-1\}$, and $X_{s_k}^{l,N} = X_{s_k}^{l,\star}$.

* If $k\in\{2, 4, \ldots, K_l\}$

* (2c) Set Brownian increment $V_{s_k}^{l-1,n} =  V_{s_{k-1}}^{l,n} + V_{s_k}^{l,n}$ at level $l-1$ for $n\in\{1,\ldots,N-1\}$.

* (2d) Set $X_{s_k}^{l-1,n} = F_{\theta}^{l-1}(X_{s_{k-1}}^{l-1,A_{s_{k-1}}^{l-1,n}}, V_{s_k}^{l-1,n})$
at level $l-1$ for $n\in\{1,\ldots,N-1\}$, and $X_{s_k}^{l-1,N} = X_{s_k}^{l-1,\star}$.

* If there is an observation at time $t=s_k\in\{1,\ldots,T\}$

* (2e) Compute normalized weights
$w_t^{l-1,n} \propto g_{\theta}(y_{t}|X_{t}^{l-1,n})$ and
$w_t^{l,n} \propto g_{\theta}(y_{t}|X_{t}^{l,n})$ for $n\in\{1,\ldots,N\}$.

* (2f) If $t < T$, sample ancestors $(A_{t}^{l-1,n},A_{t}^{l,n})
\sim\bar{\mathcal{R}}(w_{t}^{l-1,1:N},w_{t}^{l,1:N})$ independently
for $n\in\{1,\ldots,N-1\}$ and set $A_{t}^{l-1,N}=N,A_{t}^{l,N}=N$.

* Else

* (2i) If $k\in\{2, 4, \ldots, K_l\}$, set $A_{s_k}^{l-1,n}=n$ at level $l-1$ for $n\in\{1,\ldots,N\}$.

* (2j) Set $A_{s_k}^{l,n}=n$ at level $l$ for $n\in\{1,\ldots,N\}$.

After the terminal step

* (3a) Sample particle indexes $(B_{T}^{l-1},B_{T}^{l})\sim \bar{\mathcal{R}}(w_{T}^{l-1,1:N}, w_{T}^{l,1:N})$.

* (3b) Set particle indexes $B_{s_{2k}}^{l-1} = A_{s_{2k}}^{l-1,B_{s_{2(k+1)}}^{l-1}}$
at level $l-1$ for $k\in\{0,1,\ldots,K_{l-1}-1\}$.

* (3c) Set particle indexes $B_{s_k}^l = A_{s_k}^{l,B_{s_{k+1}}^l}$ at level $l$ for $k\in\{0,1,\ldots,K_l-1\}$.

Output: a trajectory $X_{0:T}^{l-1,\circ}=(X_{s_{2k}}^{l-1,B_{s_{2k}}^{l-1}})_{k=0}^{K_{l-1}}\in\mathsf{X}^{l-1}$
and a trajectory $X_{0:T}^{l,\circ}=(X_{s_k}^{l,B_{s_k}^l})_{k=0}^{K_l}\in\mathsf{X}^{l}$.

§.§ Summary of proposed methodology and choice of tuning parameters

We consolidate the algorithms presented in this section by
summarizing our proposed methodology to unbiasedly estimate the score function $S(\theta)$ below.

Input: PMF $(P_l)_{l=0}^{\infty}$, number of particles $N$, burn-in $b$ and number of iterations $I$.

* Sample highest discretization level $L$ from $(P_l)_{l=0}^{\infty}$.

* Simulate coupled CPF chains $(X_{0:T}(i),\bar{X}_{0:T}(i))$ using CPF (Algorithm <ref>) and
2-CCPF (Algorithm <ref>) with $N$ particles at discretization level $0$ as described in (<ref>)
for iteration $i=0,1,\ldots,\max\{I, \tau_{\theta}^0\}$, and
compute unbiased estimator $\widehat{I}_0(\theta)$ of ${I}_0(\theta)$ using time-averaged estimator
(<ref>) with burn-in $b$ and $I$ iterations.

* Independently for $l = 1,\ldots,L$, simulate two pairs of coupled CPF chains $(X_{0:T}^{l-1}(i),\bar{X}_{0:T}^{l-1}(i))$
and $(X_{0:T}^{l}(i),\bar{X}_{0:T}^{l}(i))$ using ML-CPF (Algorithm <ref>) and
4-CCPF (Algorithm <ref>) with $N$ particles at discretization levels $l-1$ and $l$ as described in (<ref>)
for iteration $i=0,1,\ldots,\max\{I, \bar{\tau}_{\theta}^l\}$, and
compute unbiased estimators $\widehat{S}_{l-1}(\theta)$ and $\widehat{S}_l(\theta)$ of the scores
$S_{l-1}(\theta)$ and $S_l(\theta)$ using time-averaged estimator
(<ref>) with burn-in $b$ and $I$ iterations.
Compute unbiased estimator of the score increment $I_l(\theta)$ with
$\widehat{I}_l(\theta)=\widehat{S}_l(\theta) - \widehat{S}_{l-1}(\theta)$.

* Compute unbiased estimator of the score function $S(\theta)$ using
independent sum estimator $\widehat{S}(\theta) = \sum_{l=0}^L \widehat{I}_l(\theta)/\mathcal{P}_l$,
where $\mathcal{P}_l = \sum_{k=l}^{\infty} P_k$.

Output: unbiased estimator $\widehat{S}(\theta)$ of the score function $S(\theta)$.

We will establish unbiasedness and finite variance properties of $\widehat{S}(\theta)$ in Section <ref>.
The cost of the above procedure is $c(\theta)=\sum_{l=0}^Lc_l(\theta)$, where
$c_l(\theta)$ is the cost of computing $\widehat{I}_l(\theta)$.
From Sections <ref> and <ref>, we have
$c_l(\theta)=a_\theta^l\times\mathrm{Cost}(M_{\theta}^l)$, with
$a_\theta^l=\max(2\tau_{\theta}^0-1,I+\tau_{\theta}^0-1)$ for level $l=0$, and
$a_\theta^l=\max(2\tau_{\theta}^{l-1}-1,I+\tau_{\theta}^{l-1}-1)/2 + \max(2\tau_{\theta}^l-1,I+\tau_{\theta}^l-1)$
for level $l\in\mathbb{N}$.
We take the view here that the cost per application of the CPF kernel $M_{\theta}^l$ is
$\mathrm{Cost}(M_{\theta}^l)=N K_l = N 2^l T$. Hence the expected cost of computing $\widehat{S}(\theta)$ is

\begin{align}\label{eqn:expected_cost}
\mathbb{E}\left[c(\theta)\right] = \sum_{l=0}^{\infty}\mathbb{E}\left[c_l(\theta)\right] \mathcal{P}_l
= NT\sum_{l=0}^{\infty} \mathbb{E}\left[a_{\theta}^l\right] 2^l \mathcal{P}_l.
\end{align}

We will see that one cannot find a PMF $(P_l)_{l=0}^{\infty}$ so that the variance and expected cost of $\widehat{S}(\theta)$
are both finite.
We defer further discussions and the selection of the distribution of $L$ to Section <ref>, and
assume for now that we have a given PMF $(P_l)_{l=0}^{\infty}$ that ensures finite variance
but infinite expected cost. In this regime, we will refrain from discussions of asymptotic efficiency
in the sense of Glynn and Whitt, 1992.

We now discuss the choice of tuning parameters and various algorithmic considerations.
In the above description, the choice of $(N, b, I)$ could be level-dependent but optimizing
these tuning parameters is outside the scope of this work.
Following the discussion in <cit.> and the empirical findings in Jacob et al., 2020,
we will scale the number of particles $N$ linearly with the number of observations $T$.
Although the variance of $\widehat{I}_l(\theta)$ decreases as we increase the burn-in
parameter $b\in\mathbb{N}_0$, setting $b$ too large would be inefficient.
Jacob et al., 2020, Jacob et al., 2020 proposed choosing $b$ according to the distribution of the meeting time.
In our context, as the stopping time $\bar{\tau}_{\theta}^l$ typically decreases as the level $l$ increases,
a conservative strategy is to select $b$ based on the stopping time of a low discretization level,
which can be simulated by running ML-CPF and 4-CCPF as in Step 3.
We will illustrate this numerically in Section <ref> and experiment with
various choices of $b$.

After selecting $b$, one can choose the number of iterations $I\geq b$ to further reduce the variance of
$\widehat{I}_l(\theta)$, and hence that of $\widehat{S}(\theta)$, at a cost (<ref>)
that grows linearly with $I$.
On the other hand, when employing score estimators within a stochastic gradient method,
taking large values of $I$ to obtain low variance gradient estimators would be inefficient.
Choosing the tuning parameters $(N, b, I)$ to maximize the efficiency of the resulting
stochastic gradient method is a highly non-trivial problem, and could be the topic of future work.

Given a choice of $(N, b, I)$, one can produce $R\in\mathbb{N}$ independent replicates
$\widehat{S}(\theta)_r = \sum_{l=0}^{L_r} \widehat{I}_l(\theta)_r/\mathcal{P}_l, r\in\{1,\ldots,R\}$,
of $\widehat{S}(\theta)$ in parallel,
and compute the average $\bar{S}(\theta)=R^{-1}\sum_{r=1}^R\widehat{S}(\theta)_r$ to approximate $S(\theta)$.
To see the connection to MLMC, we follow <cit.> by noting that $\bar{S}(\theta)$
has the same distribution as the random variable

\begin{align}\label{eqn:MLMC_connection}
\sum_{l=0}^{\infty}\frac{1}{\mathbb{E}\left[R_l\right]}\sum_{r=1}^{R_l}\widehat{I}_l(\theta)_r,
\end{align}

where $R_l=\sum_{r=1}^R\mathbb{I}(L_r\geq l)$ has expectation $\mathbb{E}[R_l] = R \mathcal{P}_l$.
Vihola, 2018 proposed new unbiased estimators with lower variance than $\bar{S}(\theta)$
by sampling the random variables $(L_r)_{r=1}^R$ that define $(R_l)_{l=0}^{\infty}$ in (<ref>) using stratification.
As the number of replicates $R\rightarrow\infty$, these improved estimators also attain the same limiting variance
as the idealized MLMC estimator

\begin{align}
\widetilde{S}(\theta) = \sum_{l=0}^{\infty}\frac{1}{\tilde{R}_l}\sum_{r=1}^{\tilde{R}_l}\widehat{I}_l(\theta)_r,
\end{align}

where $\tilde{R}_l = \lfloor R\mathcal{P}_l\rfloor$ is allocated using the PMF $(P_l)_{l=0}^{\infty}$.
From <cit.>, this asymptotic variance is given by
\sum_{l=0}^{\infty} \mathrm{Var}[\widehat{I}_l(\theta)^j] /\mathcal{P}_l
for all $j\in\{1,\ldots,d_{\theta}\}$.

§ ANALYSIS

This section is concerned with the theoretical validity of our approach. We first introduce some notation
needed to state the assumptions which we will rely on.
Let $(\mathsf{E},\mathcal{E})$ be an arbitrary measurable space.
We write $\mathcal{B}_b(\mathsf{E})$ as the collection of real-valued, bounded and measurable functions on $\mathsf{E}$.
For real-valued $\varphi:\mathsf{E}\rightarrow\mathbb{R}$ (and vector-valued $\varphi:\mathsf{E}\rightarrow\mathbb{R}^d$),
let $\mathcal{C}^j(\mathsf{E})$ (and $\mathcal{C}_d^j(\mathsf{E})$) denote the collection of $j\in\mathbb{N}$ times continuously differentiable functions,
and $\mathcal{C}(\mathsf{E})$ (and $\mathcal{C}_d(\mathsf{E})$) for the collection of continuous functions.
We write $\varphi\in\textrm{Lip}_{\|\cdot\|_2}(\mathbb{R}^{d})$ if the real-valued function
$\varphi:\mathbb{R}^d\rightarrow\mathbb{R}$ is Lipschitz w.r.t. the $\mathbb{L}_2$-norm $\|\cdot\|_2$, i.e.
if there exists a constant $C<\infty$ such that $|\varphi(x)-\varphi(y)| \leq C\|x-y\|_2$ for all $x,y\in\mathbb{R}^{d}$,
and $\|\varphi\|_{\textrm{Lip}}$ as the Lipschitz constant.

We recall the definitions of $\Sigma(x)=\sigma(x)\sigma(x)^*$ and
$b_\theta(x) = \Sigma(x)^{-1}\sigma(x)^*a_{\theta}(x)$ for $x\in\mathbb{R}^d$
as they appear in the following.

The drift function $a:\Theta\times\mathbb{R}^d \rightarrow\mathbb{R}^d$ and diffusion coefficient
$\sigma:\mathbb{R}^d \rightarrow\mathbb{R}^{d\times d}$ satisfy:

* (Smoothness) For any $\theta\in\Theta$, $a_{\theta}^j\in \mathcal{C}^2(\mathbb{R}^d)$
for all $j\in\{1,\ldots,d\}$ components of $a_{\theta}$, and $\sigma^{j,k} \in \mathcal{C}^2(\mathbb{R}^d)$
for all $(j,k)\in\{1,\ldots, d\}$ components of $\sigma$.
Also, for any $x\in\mathbb{R}^d$, $\theta\mapsto a_{\theta}^j(x)\in\mathcal{C}(\Theta)$ for all $j\in\{1,\dots,d\}$.

* (Uniform ellipticity) $\Sigma(x)$ is uniformly positive definite for all $x\in \mathbb{R}^d$.

* (Globally Lipschitz) For any $\theta\in\Theta$, there exists a constant $C<\infty$ such that
$|a_{\theta}^j(x)-a_{\theta}^j(x')|+|\sigma^{j,k}(x)-\sigma^{j,k}(x')| \leq C \|x-x'\|_2$ for all
$(x,x') \in \mathbb{R}^d\times\mathbb{R}^d$ and $(j,k)\in\{1,\dots,d\}^2$.

The drift function $a:\Theta\times\mathbb{R}^d \rightarrow\mathbb{R}^d$, diffusion coefficient
$\sigma:\mathbb{R}^d \rightarrow\mathbb{R}^{d\times d}$ and
observation density $g:\Theta\times\mathbb{R}^d\times\mathbb{R}^{d_y}\rightarrow\mathbb{R}^+$ satisfy:

* the inverse of $x\mapsto\Sigma(x)$ satisfies
$[\Sigma^{-1}]^{j,k}\in\mathcal{B}_b(\mathbb{R}^d)\cap\textrm{\emph{Lip}}_{\|\cdot\|_2}(\mathbb{R}^d)$ for all $(j,k)\in\{1,\dots,d\}^2$.

* For any $\theta\in\Theta$, $a_\theta^j\in\mathcal{B}_b(\mathbb{R}^d)$ for all $j\in\{1,\dots,d\}$, and
$\sigma^{j,k}\in\mathcal{B}_b(\mathbb{R}^d)$ for all $(j,k)\in\{1,\dots,d\}^2$.

* For any $\theta\in\Theta$, there exists $0<\underline{C}<\overline{C}<\infty$ such that
$\underline{C}\leq g_{\theta}(y|x)\leq \overline{C}$ for all $(x,y)\in \mathbb{R}^d\times\mathbb{R}^{d_y}$.
In addition, for any $(\theta,y)\in\Theta\times\mathbb{R}^{d_y}$, we have $g_{\theta}(y|\cdot)\in\textrm{\emph{Lip}}_{\|\cdot\|_2}(\mathbb{R}^d)$.

* For any $(\theta,y)\in\Theta\times\mathbb{R}^{d_y}$, $[\nabla_{\theta}\log g_{\theta}(y|\cdot)]^j
\in\mathcal{B}_b(\mathbb{R}^d)\cap\textrm{\emph{Lip}}_{\|\cdot\|_2}(\mathbb{R}^d)$ for all $j\in\{1,\dots,d_{\theta}\}$.

* For any $\theta\in\Theta$, $[\nabla_{\theta} [b_{\theta}]^j]^k, [\nabla_{\theta} (b_{\theta}^j)^2]^k
\in\mathcal{B}_b(\mathbb{R}^d)\cap\textrm{\emph{Lip}}_{\|\cdot\|_2}(\mathbb{R}^d)$ for all

Assumptions <ref> and <ref> should be understood as sufficient conditions
to verify the validity of our proposed methodology, and are not necessary for its implementation.
Although some of these assumptions are strong, they have been adopted to simplify the exposition
of our analysis, as is common in theoretical works on particle filtering.
Some assumptions can be relaxed at the expense of more involved and lengthy
technical arguments.

We first give an intermediate result on the convergence of the time-discretized score functions $S_l(\theta)$
defined in (<ref>).

Under Assumptions <ref> and <ref>, for any $(T,\theta)\in\mathbb{N}\times\Theta$, there exists a constant $C<\infty$ such that
for any $l\in\mathbb{N}_0$, $\|S_l(\theta) - S(\theta) \|_1 \leq C\Delta_l^{1/2}$.

We note that the upper-bound can be sharpened to $\mathcal{O}(\Delta_l)$.
The proof which involves studying the time-discretization of diffusions
can be found in Appendix <ref>.
The following is our main result that establishes basic properties of our score estimator.

Under Assumptions <ref> and <ref>, for any number of particles $N\geq 2$,
burn-in $b\in\mathbb{N}_0$ and number of iterations $I\geq b$,
there exists a choice of PMF $(P_l)_{l=0}^{\infty}$ such that for any $\theta\in\Theta$,
the score estimator $\widehat{S}(\theta)$ in (<ref>) is unbiased and has finite variance.

The proof of Theorem <ref> in Appendix <ref> involves analyzing
various aspects of the 4-CCPF chains (Appendix <ref> & <ref>)
and its initialization (Appendix <ref>), followed by the
properties of our score estimators that are constructed using these coupled Markov chains (Appendix <ref>).

It follows from the proof that the left-hand side of (<ref>) is upper-bounded by

\begin{align}\label{eqn:upperbound_criterion}
\end{align}

where $C(\theta)<\infty$ is a parameter-dependent constant and
$\phi\in(0,1/2)$ is a constant determined in our analysis that does not depend on $l$.
Hence any choice of PMF $(P_l)_{l=0}^{\infty}$ such that the sum (<ref>)
is finite would be valid; e.g. $P_l\propto \Delta_l^{2\phi\alpha}$ for any $\alpha\in(0,1)$.
Due to the technical complexity of the problem and algorithms under consideration,
the rate in (<ref>) is not sharp.
We conjecture that the correct rate corresponds to having $\phi=1/4$ and
a better rate of $\phi=1/2$ can be obtained in the case of constant diffusion coefficient $\sigma$.

Using Lemma <ref> in Appendix <ref>, we can upper-bound the expected cost (<ref>)

\begin{align}\label{eqn:upperbound_cost}
C(\theta,T,N,b,I)NT\sum_{l=0}^{\infty} 2^l \mathcal{P}_l,
\end{align}

where $C(\theta,T,N,b,I)<\infty$ is another constant that is independent of $l$.
As we are in a setting where $\phi\leq 1/2$, there is no choice of PMF $(P_l)_{l=0}^{\infty}$
that can keep both (<ref>) and (<ref>) finite
This is a consequence of employing the Euler–Maruyama discretization scheme (<ref>)
and the choice of coupled resampling scheme in Algorithm <ref>;
future work could consider the antithetic truncated Milstein scheme of Giles and Szpruch, 2014
and improved coupled resampling schemes such as Ballesio et al., 2020.
Our numerical implementations will follow the approach of Rhee and Glynn, 2015
by choosing the PMFs $P_l\propto \Delta_l^{1/2}l(\log_2(1+l))^2$
and $P_l\propto \Delta_l l (\log_2(1+l))^2$ in the case of non-constant and constant diffusion coefficients, respectively.
These choices yield score estimators with finite variance but infinite expected cost,
which is a sensible compromise when employing these gradients within stochastic gradient methods.
Moreover, under these choices, the unbiased estimators also achieve
computational complexities that are similar to MLMC estimators [Giles, 2008].

§ APPLICATIONS

§.§ Ornstein–Uhlenbeck process

We consider an Ornstein–Uhlenbeck process $X=(X_t)_{0\leq t\leq T}$ in $\mathbb{R}$, defined by the SDE

\begin{align}\label{eqn:OU_SDE}
dX_t = \theta_1(\theta_2-X_t)dt + \sigma dW_t,\quad X_0=0.
\end{align}

The parameter $\theta_{1}> 0$ can be interpreted as the speed of the mean reversion to the
long-run equilibrium value $\theta_{2}\in\mathbb{R}$.
This corresponds to (<ref>) with initial condition $x_{\star}=0$, linear drift function $a_{\theta}(x)=\theta_1(\theta_2-x)$
and constant diffusion coefficient $\sigma(x)=\sigma>0$ for $x\in\mathbb{R}$.
We assume that the process is observed at unit times with Gaussian measurement errors,
i.e. $Y_t|X\sim g_{\theta}(\cdot|X_t)=\mathcal{N}(X_t,\theta_3)$ independently for $t\in\{1,\ldots,T\}$ and some $\theta_3>0$.
We will generate observations $y_{1:T}$ by simulating from the model with
parameter $\theta=(\theta_1,\theta_2,\theta_3)=(2,7,1)$.
This setting is considered as it is possible to compute the score function (<ref>)
exactly using Kalman smoothing; see Appendix <ref> for details and
model-specific expressions to evaluate (<ref>).

Figure <ref> illustrates how the distribution of the
stopping time $\bar{\tau}_{\theta}^l$ varies with the discretization level $l$ on a simulated dataset with $T=25$ observations.
We took $l=3$ as the lowest discretization level as lower levels lead to numerically unstable trajectories.
As alluded earlier, the coupled resampling scheme proposed in Algorithm <ref> (referred to as “maximal”)
leads to smaller and more stable stopping times for large enough levels.
As alternatives, we consider a modification (“other maximal”) that ensures 4-CCPF admits 2-CCPF as marginals on each level
and a scheme that uses common uniform random variables (“common uniforms”).
While the schemes based on maximal couplings have comparable stopping times, the approach based on common uniform
variables gives rise to significantly larger stopping times.
Figure <ref> reveals that these two alternative coupled resampling schemes
do not induce sufficient dependencies between the four CPF chains.
As the variance of the estimated increment does not decrease
with the discretization level, this precludes their use within our unbiased score estimation framework.

To show the impact of the choice of $b$ and $I$, we considered three types of estimators corresponding to
having $b=0$ and $I=b$ (“naive”); $b=90\%\textrm{-quantile}(\bar{\tau}_{\theta}^l)$ at level $l=3$ and $I=b$ (“simple”);
and $b=90\%\textrm{-quantile}(\bar{\tau}_{\theta}^l)$ and $I= 10b$ (“time-averaged”), where
$90\%\textrm{-quantile}(\bar{\tau}_{\theta}^l)$ denotes the 90% sample quantile of the stopping time at level $l$.
The benefits of increasing $b$ and $I$ in terms of variance reduction are consistent with findings in Jacob et al., 2020, Jacob et al., 2020.
Under our proposed coupling, all three choices yield estimators of increments whose variance decrease exponentially with the level,
which agrees with our theoretical results (see Lemma <ref> in Appendix <ref>).
Hence we can employ any of these estimators within the score estimation methodology outlined
in Section <ref>.
Figure <ref> displays the resulting squared error $\| \widehat{S}(\theta) - S(\theta) \|_2^2$
and cost of $100$ independent replicates. This plot suggests having $N=128$ particles is sufficient in the case of $T=25$ observations.
The choice of $b$ and $I$ also allows a tradeoff between error and cost.
As we increase the number of observations $T$, Figure <ref> shows it is important
to scale the number of particles $N$ linearly with $T$ to obtain stable and non-exponentially increasing stopping times.
Lastly, Figures <ref> and <ref> concern the averaging of independent replicates of the score estimator
$(\widehat{S}(\theta)_r)_{r=1}^R$. Figure <ref> shows that the average $\bar{S}(\theta)=R^{-1}\sum_{r=1}^R\widehat{S}(\theta)_r$
satisfies the standard Monte Carlo rate as $R\rightarrow\infty$, which is consistent with its properties in Theorem <ref>,
at a linear cost in $R$ as illustrated in Figure <ref>.

Boxplots of stopping time $\bar{\tau}_{\theta}^l$ against discretization level $l$

Variance of score increment estimator summed over components
$\sum_{j=1}^{d_{\theta}}\mathrm{Var}[\widehat{I}_l(\theta)^j]$ against discretization level $l$

Squared error $\| \widehat{S}(\theta) - S(\theta) \|_2^2$ against cost with $N\in\{128,512\}$ particles

Boxplots of stopping time $\bar{\tau}_{\theta}^l$ against number of observations $T$ for discretization level $l=6$ with $N=128$ fixed or linearly increasing number of particles

Mean squared error $\mathrm{MSE}(\theta) = \mathbb{E}[\| \bar{S}(\theta) - S(\theta) \|_2^2]$ against number of averaged replicates $R$

Boxplots of cost against number of averaged replicates $R$

Behaviour of different coupling methods and score estimators at the data generating parameter $\theta=(2,7,1)$ of the Ornstein–Uhlenbeck model in Section <ref>. $T=25$ observations and $N=128$ particles were employed unless stated otherwise. These plots are based on $100$ independent repetitions.

§.§ Logistic diffusion model for population dynamics of red kangaroos

Next we consider an application from population ecology to model the dynamics
of a population of red kangaroos (Macropus rufus) in New South Wales, Australia.
Figure <ref> displays data $y_{t_1},\ldots,y_{t_P}\in\mathbb{N}_{0}^{2}$
from Caughley et al., 1987, which are double transect counts on $P=41$ occasions at
irregular times $(t_p)_{p=1}^P$ between $1973$ to $1984$.
The latent population size $Z=(Z_t)_{t_1\leq t\leq t_P}$ is assumed to follow a logistic diffusion process
with environmental variance [Dennis and Costantino, 1988, Knape and De Valpine, 2012]

\begin{align}\label{eqn:logistic_diffusion}
\end{align}

where $\mathcal{LN}$ denotes the log-Normal distribution.
The parameters $\theta_{1}\in\mathbb{R}$ and $\theta_{2}>0$ can be seen as coefficients describing how the growth rate depends on the population size.
As the parameter $\theta_{3}>0$ appears in the diffusion coefficient of (<ref>),
we apply the Lamperti transformation $X_t = \Psi(Z_t)=\log(Z_t)/\theta_3$.
By Itô's lemma, the transformed process $X=(X_t)_{t_1\leq t \leq t_P}$ satisfies the SDE (<ref>)
with random initialization $X_{t_1}\sim\mu_{\theta}=\mathcal{N}(5/\theta_3,10^2/\theta_3^{2})$,
drift function $a_{\theta}(x)=\theta_1/\theta_3-(\theta_2/\theta_3)\exp(\theta_3x)$ and
unit diffusion coefficient $\sigma(x)=1$ for $x\in\mathbb{R}$.
The observations $(Y_{t_p})_{p=1}^P$ are modelled as conditionally independent given $X$ and
negative Binomial distributed, i.e.
the observation density at time $t\in\{t_1,\ldots,t_P\}$ is
where $\theta_{4}>0$.
We will use a parameterization of the negative Binomial distribution that is common in ecology,
$\mathcal{NB}(y;r,\mu)=\frac{\Gamma(y+r)}{\Gamma(r)y!}(\frac{r}{r+\mu})^{r}(\frac{\mu}{r+\mu})^{y}$ for $y\in\mathbb{N}_{0}$, where $r>0$ is the dispersion parameter and $\mu>0$ is the mean parameter. The $d_{\theta}=4$ unknown parameters to be inferred are

Application of our methodology requires some minor modifications.
As the initial distribution $\mu_{\theta}$ depends on $\theta_{3}$,
the representation of score functions in (<ref>) and (<ref>) require
adding $\nabla_{\theta}\log\mu_{\theta}(X_{t_1})$ to (<ref>) and (<ref>);
see Appendix <ref> for model-specific expressions.
To deal with irregular observation times $(t_p)_{p=1}^P$,
we set the step-size at discretization level zero as the size of the smallest time interval, i.e. $\Delta_{0}=\min_{p=2,\ldots,P}t_{p}-t_{p-1}$.
Higher levels $l\in\mathbb{N}$ will employ $\Delta_{l}=\Delta_{0}2^{-l}$.
For level $l\in\mathbb{N}_{0}$, the first time interval $[t_1,t_2]$ is discretized using $\Delta_{l}$ sequentially,
i.e. we set $s_{k}=t_1+k\Delta_{l}$ for $k\in\{0,\ldots,m_{l,1}\}$ with
$m_{l,1}=\left\lfloor (t_{2}-t_{1})/\Delta_{l}\right\rfloor$,
and $s_{k}=t_{2}$ for $k=m_{l,1}+1$ if $(t_{2}-t_{1})/\Delta_{l}\notin\mathbb{N}$.
The subsequent time intervals are then discretized in the same manner.

Figure <ref> illustrates how the median and the 90% quantile of the
stopping time $\bar{\tau}_{\theta}^l$ vary with the discretization level $l$, the impact of
the number of particles $N$ and the benefits of employing adaptive resampling.
As before, the coupled resampling scheme in Algorithm <ref>
results in stopping times that are smaller for higher discretization levels,
with less variability over levels as the number of particles increases.
Moreover, resampling only when the effective sample size is less than $N/2$ allows us to induce
more dependencies between the multiple CPF chains at lower discretization levels.
Using $N=256$ particles and adaptive resampling, Figure <ref> examines
the rate at which the variance of the estimated increment decreases with the discretization level.
Here we consider the “naive” and “simple” estimators described in Section <ref>,
with a burn-in of $b=90\%\textrm{-quantile}(\bar{\tau}_{\theta}^l)$ at level $l=3$,
and omit the more costly “time-averaged” estimator.
From the plot, both type of estimators have similar rate of decay and are valid choices
in our score estimation methodology.
Using the “simple” estimator, Figure <ref> verifies that
the average of $R$ independent replicates of the resulting score estimator $\bar{S}(\theta)$
satisfies the standard Monte Carlo rate as $R\rightarrow\infty$.

Lastly, we perform Bayesian parameter inference by employing our score estimators
within the SGLD framework [Welling and Teh, 2011]. We rely on logarithmic transformations
to deal with positivity parameter constraints, and specify the prior distribution
$(\theta_1, \log \theta_2, \log \theta_3, \log \theta_4)\sim \mathcal{N}_{d_{\theta}}(\mu_0,\Sigma_0)$,
where $\mu_0=(0, -1, -1, -1)$ and $\Sigma_0=\textrm{diag}(5^2, 2^2, 2^2, 2^2)$.
As $\log \theta_3$ has a significantly different scale compared to the other parameters,
we let the learning rate in (<ref>) be component-dependent
by taking $\varepsilon_m=\textrm{diag}((100+m)^{-0.6}(10^{-2}, 10^{-2}, 10^{-4}, 10^{-2}))$ at iteration $m\geq 1$.
The algorithmic settings used to produce score estimators are the same as in Figure <ref>
with $R=1$ realization.
Figure <ref> shows the trace plot over $7500$ iterations of the resulting SGLD algorithm
for each parameter.

Double transect counts

Sample median and 90% quantile of stopping time $\bar{\tau}_{\theta}^l$ against discretization level $l$

Variance of score increment estimator summed over components
$\sum_{j=1}^{d_{\theta}}\mathrm{Var}[\widehat{I}_l(\theta)^j]$ against discretization level $l$

Mean squared error $\mathrm{MSE}(\theta)^j = \mathbb{E}[|\bar{S}(\theta)^j - S(\theta)^j |^2]$ against number of averaged replicates $R$
for component $j\in\{1,\ldots,d_{\theta}\}$

SGLD for each component.

Behaviour at parameter $\theta=(2.397, 4.429 \times 10^{-3}, 0.840, 17.631)$ of the logistic diffusion model in Section <ref>. $N=256$ particles were employed unless stated otherwise. These plots are based on $100$ independent repetitions.

§.§ Neural network model for grid cells in the medial entorhinal cortex

As our final application, we consider a neural network model for single neurons to analyze grid
cells spike data[<https://www.ntnu.edu/kavli/research/grid-cell-data>]
recorded in the medial entorhinal cortex of rats that were running
on a linear track [Hafting et al., 2008]. The neural states $Z_{t}=(Z_{t}^{1},Z_{t}^{2})$
of two grid cells that were simultaneously recorded is assumed to follow

\begin{align}\label{eqn:nn_model}
dZ_{t}^{1} & =\left(\alpha_{1}\tanh(\beta_{1}Z_{t}^{2}+\gamma_{1})-\delta_{1}Z_{t}^{1}\right)dt+\sigma_{1}dW_{t}^{1},\\
dZ_{t}^{2} & =\left(\alpha_{2}\tanh(\beta_{2}Z_{t}^{1}+\gamma_{2})-\delta_{2}Z_{t}^{2}\right)dt+\sigma_{2}dW_{t}^{2},\notag
\end{align}

for $t\in[0,T]$, where $(\alpha_{1},\alpha_{2})\in\mathbb{R}^{2}$
controls the amplitude, $(\beta_{1},\beta_{2})\in\mathbb{R}^{2}$
describes the connectivity between the cells, $(\gamma_{1},\gamma_{2})\in\mathbb{R}^{2}$
are baseline levels, $(\delta_{1},\delta_{2})\in(0,\infty)^{2}$
determines the strength of the mean reversion towards the origin.
We assume $Z_{0}=(0,0)$ at the beginning of the experiment.
This diffusion is motivated
by an example in [Kappen and Ruiz, 2016], and modified for our purposes.
To infer the unknown diffusivity parameters $(\sigma_{1},\sigma_{2})\in(0,\infty)^2$,
we consider the transformation
which rescales each component of the diffusion. By Itô's formula,
the transformed process $X=(X_t)_{0\leq t\leq T}$ satisfies the diffusion model (<ref>)
with initialization $x_{\star}=(0,0)$, drift function

\begin{align}
\end{array}\right)=\left(\begin{array}{c}
\alpha_{1}\tanh(\beta_{1}\sigma_{2}x^{2}+\gamma_{1})/\sigma_{1}-\delta_{1}x^{1}\\
\alpha_{2}\tanh(\beta_{2}\sigma_{1}x^{1}+\gamma_{2})/\sigma_{2}-\delta_{2}x^{2}
\end{array}\right),
\end{align}

and diffusion coefficient $\sigma(x)=I_{2}$ for $x=(x^1,x^2)\in\mathbb{R}^2$.

The experimental data over a duration of $T=20$ seconds contains time stamps
in $[0,T]$ when a spike at one of the two cells is recorded using
tetrodes. Following [Brown, 2005], we adopt an inhomogenous Poisson
point process to model these times.
Let $t_p = p T 2^{-6}$ for $p\in\{0,1,\ldots,P\}$ denote a dyadic uniform discretization of $[0,T]$
into $P=2^6$ time intervals.
Given the latent process $X=(X_t)_{0\leq t\leq T}$, the number of spikes $Y_{t_{p}}^{i}$ occurring
in each time interval $[t_{p-1},t_{p}]$ at cell $i=1,2$
is assumed to be conditionally independent of the other time intervals and
the activity in the other cell, and follow a Poisson distribution
with rate $\int_{t_{p-1}}^{t_{p}}\lambda_{i}(X_{t}^{i})dt$. The intensity
function for grid cell $i=1,2$ is modelled as $\lambda_{i}(X_{t}^{i})=\exp(\kappa_{i}+X_{t}^{i})$,
where $\kappa_{i}\in\mathbb{R}$ represents a baseline level.
The observed counts $y_{t_p}=(y_{t_{p}}^{1},y_{t_{p}}^{2})$ for interval $p\in\{1,\ldots,P\}$,
computed from the experimental data, are displayed in Figure <ref>.
The conditional likelihood of the observation model is
$p_{\theta}(y_{t_{1}},\ldots,y_{t_{P}}|X) = \prod_{p=1}^{P}g_{\theta}(y_{t_p}| (X_t)_{t_{p-1}\leq t\leq t_{p}})$
with the intractable observation density

\begin{align}
g_{\theta}(y_{t_p}| (X_t)_{t_{p-1}\leq t\leq t_{p}}) =
\prod_{i=1}^2\mathcal{P}oi\left(y_{t_p}^i;\int_{t_{p-1}}^{t_p} \lambda_i(X^i_t)dt\right),
\end{align}

where $\mathcal{P}oi(y;\lambda)=\lambda^y\exp(-\lambda)/y!$ for $y\in\mathbb{N}_0$
denotes the PMF of a Poisson distribution with rate $\lambda>0$.
To approximate the conditional likelihood, at level $l \geq 6$, we discretize the time interval $[0,T]$
in a similar manner using $s_k=k \Delta_l$ for $k\in\{0,1,\ldots,K_l\}$, where $\Delta_l=T2^{-l}$
is the step-size and $K_l=2^l$ is the number of time steps.
Under the time-discretized process $X_{0:T}=(X_{s_k})_{k=0}^{K_l}$,
the resulting approximation of the conditional likelihood is
$p^l_{\theta}(y_{t_{1}},\ldots,y_{t_{P}}| X_{0:T}) = \prod_{p=1}^P g^l_{\theta} (y_{t_p} | (X_t)_{t_{p-1}\leq t\leq t_{p}})$
with the corresponding observation density

\begin{align}\label{eqn:NN_obs_density}
g^l_{\theta} (y_{t_p} | (X_t)_{t_{p-1}\leq t\leq t_{p}}) =
\prod_{i=1}^2\mathcal{P}oi\left(y_{t_p}^i; \Delta_l\sum_{t: t_{p-1}\leq t\leq t_p}\lambda_i(X^i_t)\right).
\end{align}

By using these level-dependent observation densities (<ref>) in Section <ref>,
our proposed methodology can then be applied.
There are $d_{\theta}=12$ parameters $\theta=(\theta_1,\theta_2)$ to be inferred, where
$\theta_i = (\alpha_{i},\beta_{i},\gamma_{i},\delta_{i},\sigma_{i},\kappa_{i})$ denote the parameters
associated to cell $i=1,2$.
We refer the reader to Appendix <ref> for model-specific expressions
to evaluate (<ref>).

We consider an extension of the proposed method based on the conditional ancestor sampling particle filter
(CASPF) [Lindsten et al., 2014] as the basic algorithmic building block.
As CASPF has better mixing properties than CPF, we observe smaller stopping times in Figure <ref>
for lower discretization levels.
Figure <ref> verifies the validity of using both MCMC algorithms and
“naive” and “simple” estimators (as described in Section <ref>) within our methodology.
For “simple” estimators, the burn-in was taken as $b=90\%\textrm{-quantile}(\bar{\tau}_{\theta}^l)$ at level $l=11$.
The rate at which the variance of the estimated increment decays with the discretization level
is similar across algorithms and estimators.

Lastly, we combine our score estimators and the SGA scheme in (<ref>) to perform
maximum likelihood estimation.
The score estimation relies on the CASPF algorithm and the “simple” estimator with a burn-in of $b=100$.
Positivity parameter constraints are dealt with using logarithmic transformations and
a constant learning rate of $\varepsilon_m = 10^{-3}$ is employed.
Figure <ref> illustrates how the distribution of the Polyak–Ruppert average
evolves over the iterations, estimated using independent runs of SGA.
We note that only $87$ out of $100$ runs were considered, as there were
$13$ instances of the variance of the score estimator driving the SGA algorithm to
regions of the parameter space where the stopping times are prohibitively large,
causing the SGA to stall.
This pathological behaviour is due to poor mixing properties of the
underlying CASPF algorithm at very unlikely regions of the parameter space.
To improve the MCMC algorithm, one has to simulate particle dynamics
that incorporate information from the observation sequence $y_{t_{1}},\ldots,y_{t_{P}}$,
which could be investigated in future work.
The parameter estimates of $\beta_1$ and $\beta_2$ support the use of a joint model
(<ref>) for both grid cells, and indicates that these cells are positively dependent.

Counts $y_{t_p}=(y_{t_{p}}^{1},y_{t_{p}}^{2})$ on time intervals of duration $T2^{-6}=0.3125$ second

Sample median and 90% quantile of stopping time $\bar{\tau}_{\theta}^l$ against discretization level $l$

Variance of score increment estimator summed over components
$\sum_{j=1}^{d_{\theta}}\mathrm{Var}[\widehat{I}_l(\theta)^j]$ against discretization level $l$

Mean ($\pm$ one standard deviation) of Polyak–Ruppert average for
parameter estimates of $\theta_1$ over $87$ runs of SGA

Behaviour at parameter $\theta=(1,\dots,1)$ of the neural network model in Section <ref>.
The algorithmic settings involve $N=256$ particles and adaptive resampling.
These plots are based on $1000$ independent repetitions unless stated otherwise.

§ PARAMETER DEPENDENCE IN DIFFUSION COEFFICIENT

We end this article by considering how one can extend the ideas in this article to accommodate
the case where the diffusion coefficient also depends on the parameter $\theta \in \Theta$.
In this case, we have the SDE

\begin{align}
dX_t = a_{\theta}(X_t)dt + \sigma_{\theta}(X_t)dW_t,\quad X_0=x_{\star}\in\mathbb{R}^d,
\end{align}

where $\sigma:\Theta\times\mathbb{R}^d\rightarrow\mathbb{R}^{d\times d}$ is assumed to be such that
$\theta\mapsto\sigma_{\theta}$ is invertible and satisfies the conditions in Assumption <ref> uniformly in $\theta\in\Theta$.
Moreover, we shall suppose that all the conditions in Schauer et al., 2017 hold. For $k\in\mathbb{N}_0$ and $t\in[k,k+1]$,
consider the diffusion bridge

\begin{align}\label{eqn:diffusion_bridge}
dX_t = a_{\theta}^{\circ}(X_t)dt + \sigma_{\theta}(X_t)dW_t,\quad X_k=x_{k},~X_{k+1}=x_{k+1},
\end{align}

where the drift function $a_{\theta}^{\circ}$ is described in Schauer et al., 2017.
Given a Brownian path $\bm{W}_k=(W_t)_{k\leq t\leq k+1}$, we denote the path-wise solution of the diffusion bridge
as $F_{\theta,k}(\bm{W}_k,x_k,x_{k+1})$.
Furthermore, for $G:\Theta\times\mathbb{R}^d\rightarrow\mathbb{R}$ given in <cit.> and a process
$\bm{Z}_k=(Z_t)_{k\leq t \leq k+1}$, we define the functional

\begin{align}
H_{\theta,k}(\bm{Z}_k) = \int_{k}^{k+1} G_{\theta}(Z_t)dt.
\end{align}

Using the change of measure in Schauer et al., 2017 along with the approach in Beskos et al., 2021
and Yonekura and Beskos, 2020, one can write the marginal likelihood of observations $y_{1:T}=(y_{t})_{t=1}^T$ as

\begin{align}\label{eqn:general_marginal_likelihood}
p_{\theta}(y_{1:T}) = \widetilde{\mathbb{E}}_{\theta}\left[\prod_{t=1}^T g_{\theta}(y_t|X_t)
\exp\Big\{\sum_{t=0}^{T-1} H_{\theta,t}(F_{\theta,t}(\bm{W}_t,X_t,X_{t+1}))\Big\}\right].
\end{align}

In the above, $\widetilde{\mathbb{E}}_{\theta}$ denotes expectation w.r.t. the probability measure $\widetilde{\mathbb{P}}_{\theta}$
defined as

\begin{align}
\widetilde{\mathbb{P}}_{\theta}(d(x_1,\dots,x_T,\bm{W}_0,\dots,\bm{W}_{T-1}))
= \prod_{t=1}^T \bigg\{\widetilde{p}_{\theta}(x_{t-1},x_t)dx_t\bigg\}\widetilde{\mathbb{W}}(d(\bm{W}_0,\dots,\bm{W}_{T-1})),
\end{align}

where $\widetilde{p}_{\theta}(x_{t-1},x_t)$ is the transition density of an auxiliary process on a unit time interval
as constructed in [Schauer et al., 2017] that is known and can be sampled, and
$\widetilde{\mathbb{W}}(d(\bm{W}_0,\dots,\bm{W}_{T-1})) = \bigotimes_{k=0}^{T-1} \mathbb{W}(d\bm{W}_{k})$
is given by the Wiener measure $\mathbb{W}$.

Under regularity conditions, one can differentiate (<ref>) and represent the score function as

\begin{align}
S(\theta) =
\breve{\mathbb{E}}_{\theta}\left[
\sum_{t=1}^T \nabla_{\theta}\log g_{\theta}(y_t|X_t) +
\sum_{t=1}^T \nabla_{\theta}\log\widetilde{p}_{\theta}(X_{t-1},X_t) +
\sum_{t=1}^{T-1} \nabla_{\theta} H_{\theta,t}(F_{\theta,t}(\mathbf{W}_t,X_t,X_{t+1}))
\right],
\end{align}

where $\breve{\mathbb{E}}_{\theta}$ denotes expectation w.r.t. the probability measure

\begin{align}
&\breve{\mathbb{P}}_{\theta}(d(x_1,\dots,x_T,\bm{W}_0,\dots,\bm{W}_{T-1})) \notag\\
=~&p_{\theta}(y_{1:T})^{-1}\prod_{t=1}^T g_{\theta}(y_t|x_t)
\exp\Big\{
\sum_{t=0}^{T-1} H_{\theta,t}(F_{\theta,t}(\mathbf{W}_t,x_t,x_{t+1}))\Big\}
\widetilde{\mathbb{P}}_{\theta}(d(x_1,\dots,x_T,\bm{W}_0,\dots,\bm{W}_{T-1})).
\end{align}

Practical implementation will require a discretized approximation of $\nabla_{\theta} H_{\theta,p}(F_{\theta,p}(\cdot))$,
which involves a gradient w.r.t. $\theta$ of a path-wise solution of the diffusion bridge (<ref>).
Although Euler-type approximations can be obtained, the resulting bias in the sense of Theorem <ref>
is significantly more complicated to analyze and is thus left as future work.
We stress that only small modifications to our proposed methodology is necessary
to obtain unbiased estimators of the score function in this case.
A similar approach is considered in [Beskos et al., 2021] for a class of continuous-time models.
Alternatively, one could also consider using Malliavan techniques (see e.g. Fournié et al., 1999),
instead of the ideas described here.

§.§.§ Acknowledgements

Ajay Jasra was supported by KAUST baseline funding.
Jeremy Heng was funded by CY Initiative of Excellence (grant “Investissements
d'Avenir” ANR-16-IDEX-0008).

[Andrieu et al., 2010]
C. Andrieu, A. Doucet, and R. Holenstein.
Particle Markov chain Monte Carlo methods.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 720 (3):0 269–342, 2010.

[Andrieu et al., 2018]
C. Andrieu, A. Lee, and M. Vihola.
Uniform ergodicity of the iterated conditional SMC and geometric
ergodicity of particle Gibbs samplers.
Bernoulli, 240 (2):0 842–872, 2018.

[Ballesio et al., 2020]
M. Ballesio, A. Jasra, E. von Schwerin, and R. Tempone.
A Wasserstein coupled particle filter for multilevel estimation.
arXiv preprint arXiv:2004.03981, 2020.

[Beskos and Roberts, 2005]
A. Beskos and G. O. Roberts.
Exact simulation of diffusions.
The Annals of Applied Probability, 150 (4):0
2422–2444, 2005.

[Beskos et al., 2006]
A. Beskos, O. Papaspiliopoulos, and G. O. Roberts.
Retrospective exact simulation of diffusion sample paths with
Bernoulli, 120 (6):0 1077–1098,

[Beskos et al., 2006]
A. Beskos, O. Papaspiliopoulos, G. O. Roberts, and P. Fearnhead.
Exact and computationally efficient likelihood-based estimation for
discretely observed diffusion processes (with discussion).
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 680 (3):0 333–382, 2006b.

[Beskos et al., 2009]
A. Beskos, O. Papaspiliopoulos, and G. Roberts.
Monte Carlo maximum likelihood estimation for discretely observed
diffusion processes.
The Annals of Statistics, 370 (1):0 223–245,

[Beskos et al., 2021]
A. Beskos, D. Crisan, A. Jasra, N. Kantas, and H. Ruzayqat.
Score-based parameter estimation for a class of continuous-time state
space models.
SIAM Journal on Scientific Computing (to appear,
arXiv:2008.07803), 2021.

[Blanchet and Zhang, 2017]
J. Blanchet and F. Zhang.
Exact simulation for multivariate Itô diffusions.
arXiv preprint arXiv:1706.05124, 2017.

[Briers et al., 2010]
M. Briers, A. Doucet, and S. Maskell.
Smoothing algorithms for state–space models.
Annals of the Institute of Statistical Mathematics,
620 (1):0 61, 2010.

[Brown, 2005]
E. N. Brown.
The theory of point processes for neural systems.
In C. Chow, B. Gutkin, D. Hansel, C. Meunier, and J. Dalibard,
editors, Methods and Models in Neurophysics, pages 691–726. Elsevier,

[Cappé et al., 2006]
O. Cappé, E. Moulines, and T. Rydén.
Inference in hidden Markov models.
Springer Science & Business Media, 2006.

[Caughley et al., 1987]
G. Caughley, N. Shepherd, and J. Short.
Kangaroos: their ecology and management in the sheep rangelands
of Australia.
Cambridge University Press, 1987.

[Chopin and Singh, 2015]
N. Chopin and S. S. Singh.
On particle Gibbs sampling.
Bernoulli, 210 (3):0 1855–1883, 2015.

[Crisan, 2011]
D. Crisan.
Discretizing the continuous time filtering problem. Order of
The Oxford Handbook of Nonlinear Filtering, pages 572–597,

[Dennis and Costantino, 1988]
B. Dennis and R. F. Costantino.
Analysis of steady-state populations with the gamma abundance model:
application to Tribolium.
Ecology, 690 (4):0 1200–1213, 1988.

[Douc et al., 2011]
R. Douc, A. Garivier, E. Moulines, and J. Olsson.
Sequential Monte Carlo smoothing for general state space hidden
Markov models.
The Annals of Applied Probability, 210 (6):0
2109–2145, 2011.

[Fearnhead et al., 2008]
P. Fearnhead, O. Papaspiliopoulos, and G. O. Roberts.
Particle filters for partially observed diffusions.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 700 (4):0 755–777, 2008.

[Fearnhead et al., 2010]
P. Fearnhead, D. Wyncoll, and J. Tawn.
A sequential smoothing algorithm with linear computational cost.
Biometrika, 970 (2):0 447–464, 2010.

[Fearnhead et al., 2017]
P. Fearnhead, K. Latuszynski, G. O. Roberts, and G. Sermaidis.
Continuous-time importance sampling: Monte Carlo methods which
avoid time-discretisation error.
arXiv preprint arXiv:1712.06201, 2017.

[Fournié et al., 1999]
E. Fournié, J.-M. Lasry, J. Lebuchoux, P.-L. Lions, and N. Touzi.
Applications of Malliavin calculus to Monte Carlo methods in
Finance and Stochastics, 30 (4):0 391–412,

[Giles, 2008]
M. B. Giles.
Multilevel Monte Carlo path simulation.
Operations research, 560 (3):0 607–617, 2008.

[Giles and Szpruch, 2014]
M. B. Giles and L. Szpruch.
Antithetic multilevel Monte Carlo estimation for
multi-dimensional SDEs without Lévy area simulation.
The Annals of Applied Probability, 240 (4):0
1585–1620, 2014.

[Glynn and Rhee, 2014]
P. W. Glynn and C.-h. Rhee.
Exact estimation for Markov chain equilibrium expectations.
Journal of Applied Probability, 510 (A):0
377–389, 2014.

[Glynn and Whitt, 1992]
P. W. Glynn and W. Whitt.
The asymptotic efficiency of simulation estimators.
Operations research, 400 (3):0 505–520, 1992.

[Gordon et al., 1993]
N. J. Gordon, D. J. Salmond, and A. F. Smith.
Novel approach to nonlinear/non-Gaussian Bayesian state
In IEE proceedings F (radar and signal processing), volume
140, pages 107–113. IET, 1993.

[Hafting et al., 2008]
T. Hafting, M. Fyhn, T. Bonnevie, M.-B. Moser, and E. I. Moser.
Hippocampus-independent phase precession in entorhinal grid cells.
Nature, 4530 (7199):0 1248–1252, 2008.

[Heng et al., 2021]
J. Heng, A. Jasra, K. J. Law, and A. Tarakanov.
On unbiased estimation for discretized models.
arXiv preprint arXiv:2102.12230, 2021.

[Ikeda and Watanabe, 2014]
N. Ikeda and S. Watanabe.
Stochastic differential equations and diffusion processes.
Elsevier, 2014.

[Jacob et al., 2020]
P. E. Jacob, F. Lindsten, and T. B. Schön.
Smoothing with couplings of conditional particle filters.
Journal of the American Statistical Association, 1150
(530):0 721–729, 2020a.

[Jacob et al., 2020]
P. E. Jacob, J. O'Leary, and Y. F. Atchadé.
Unbiased Markov chain Monte Carlo methods with couplings.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 820 (3):0 543–600, 2020b.

[Jasra and Yu, 2020]
A. Jasra and F. Yu.
Central limit theorems for coupled particle filters.
Advances in Applied Probability, 520 (3):0
942–1001, 2020.

[Jasra et al., 2017]
A. Jasra, K. Kamatani, K. J. Law, and Y. Zhou.
Multilevel particle filters.
SIAM Journal on Numerical Analysis, 550 (6):0
3068–3096, 2017.

[Kappen and Ruiz, 2016]
H. J. Kappen and H. C. Ruiz.
Adaptive importance sampling for control and inference.
Journal of Statistical Physics, 1620 (5):0
1244–1266, 2016.

[Kloeden and Platen, 2013]
P. E. Kloeden and E. Platen.
Numerical solution of stochastic differential equations,
volume 23.
Springer Science & Business Media, 2013.

[Knape and De Valpine, 2012]
J. Knape and P. De Valpine.
Fitting complex population models by combining particle filters with
Markov chain Monte Carlo.
Ecology, 930 (2):0 256–263, 2012.

[Kushner and Yin, 2003]
H. Kushner and G. G. Yin.
Stochastic approximation and recursive algorithms and
applications, volume 35.
Springer Science & Business Media, 2003.

[Lee et al., 2020]
A. Lee, S. S. Singh, and M. Vihola.
Coupled conditional backward sampling particle filter.
Annals of Statistics, 480 (5):0 3066–3089,

[Lindsten et al., 2014]
F. Lindsten, M. I. Jordan, and T. B. Schon.
Particle Gibbs with ancestor sampling.
Journal of Machine Learning Research, 15:0 2145–2184,

[Lindsten et al., 2015]
F. Lindsten, R. Douc, and E. Moulines.
Uniform ergodicity of the particle Gibbs sampler.
Scandinavian Journal of Statistics, 420 (3):0
775–797, 2015.

[Mcleish, 2011]
D. Mcleish.
A general method for debiasing a Monte Carlo.
Monte Carlo Methods and Applications, 17:0 301–315,

[Pollock et al., 2016]
M. Pollock, A. M. Johansen, and G. O. Roberts.
On the exact and $\varepsilon$-strong simulation of (jump)
Bernoulli, 220 (2):0 794–856, 2016.

[Rhee and Glynn, 2015]
C.-H. Rhee and P. W. Glynn.
Unbiased estimation with square root convergence for SDE models.
Operations Research, 630 (5):0 1026–1043,

[Rogers and Williams, 2000]
L. C. G. Rogers and D. Williams.
Diffusions, Markov processes and martingales: Volume 2, Itô
calculus, volume 2.
Cambridge university press, 2000.

[Schauer et al., 2017]
M. Schauer, F. Van Der Meulen, and H. Van Zanten.
Guided proposals for simulating multi-dimensional diffusion bridges.
Bernoulli, 230 (4A):0 2917–2950, 2017.

[Sørensen, 2004]
H. Sørensen.
Parametric inference for diffusion processes observed at discrete
points in time: a survey.
International Statistical Review, 720 (3):0
337–354, 2004.

[Tadić and Doucet, 2017]
V. B. Tadić and A. Doucet.
Asymptotic bias of stochastic gradient search.
The Annals of Applied Probability, 270 (6):0
3255–3304, 2017.

[Teh et al., 2016]
Y. W. Teh, A. H. Thiery, and S. J. Vollmer.
Consistency and fluctuations for stochastic gradient Langevin
Journal of Machine Learning Research, 17, 2016.

[Thorisson, 2000]
H. Thorisson.
Coupling, stationarity, and regeneration.
Springer New York, 2000.

[Vihola, 2018]
M. Vihola.
Unbiased estimators and multilevel Monte Carlo.
Operations Research, 660 (2):0 448–462, 2018.

[Welling and Teh, 2011]
M. Welling and Y. W. Teh.
Bayesian learning via stochastic gradient Langevin dynamics.
In Proceedings of the 28th international conference on machine
learning (ICML-11), pages 681–688. Citeseer, 2011.

[Yonekura and Beskos, 2020]
S. Yonekura and A. Beskos.
Online smoothing for discretely observed jump diffusions: Application
to parameter estimation and model selection.
Technical Report, UCL, 2020.

§ THEORETICAL ANALYSIS

§.§ Introduction and preliminaries

Section <ref> provides some results on time-discretization of diffusions, which are needed for the proofs associated to
Theorem <ref> as well as the 4-CCPF (Algorithm <ref>).
Our main technical arguments associated to Theorem <ref> are given in Section <ref>,
followed by several remarks about the proofs and discussions of alternative strategies.
This section of the appendix is intended to be read in the order in which it is presented.
Some familiarity with the approach in [Jasra et al., 2017] is also useful.

Note that our results concerning $\mathbb{L}_r$-norms are stated for $r\in[1,\infty)$; and can be extended to the case $r\in(0,1)$ by Hölder's inequality. We will use this fact without further elaboration. Throughout our arguments, $C$ will represent a finite constant whose value may change from line to line, but does not depend upon the discretization level. Any other dependencies in the various parameters considered will be made explicit in the statement of our results.

§.§ Results on time-discretized diffusion processes

In this section, we consider two diffusion process $X = (X_t)_{t\geq 0}$ and $X^{\star}=(X_t^{\star})_{t\geq 0}$ on the filtered probability space $(\Omega,\mathscr{F},\{\mathscr{F}_t\}_{t\geq 0},\mathbb{P}_{\theta})$ following (<ref>), with the respective initial conditions $X_0=x\in\mathbb{R}^d$ and $X_0^{\star}=x_{\star}\in\mathbb{R}^d$, and driven by the same Brownian motion.
We will consider Euler discretizations (<ref>) of $(X_t)_{t\geq 0}$ and $(X_t^{\star})_{t\geq 0}$ at some given level $l$, denoted as $\widetilde{X}_{0:T}$ and $\widetilde{X}_{0:T}^{\star}$, driven by the same Brownian motion and with the initial conditions $\widetilde{X}_0=x$ and $\widetilde{X}_0^{\star}=x_{\star}$.
The expectation operator for the described processes is written as $\mathbb{E}_{\theta}$.
Although many proof strategies in this section follow Crisan, 2011, we note that many of these arguments will be used in Section <ref> to study coupled CPFs.

In addition to the previously defined terms $b_\theta(x) = \Sigma(x)^{-1}\sigma(x)^*a_{\theta}(x)$ and $\Sigma(x)=\sigma(x)\sigma(x)^*$,
we introduce the function $\rho_{\theta}(x) = b_{\theta}(x)^*\Sigma(x)^{-1}\sigma(x)^*$ which allows us
to rewrite (<ref>) and (<ref>) as

\begin{align}
G_{\theta}(X) &= -\frac{1}{2}\int_{0}^{T}\nabla_{\theta} \| b_{\theta}(X_t) \|_2^2 dt
+ \int_{0}^{T}\nabla_{\theta}\rho_{\theta}(X_{t})dX_{t} +
\sum_{t=1}^T\nabla_{\theta} \log g_{\theta}(y_{t}|X_{t}),\\
= &-\frac{1}{2}\sum_{k=1}^{K_l} \nabla_{\theta} \| b_{\theta}(X_{s_{k-1}}) \|_2^2 \Delta_l
+ \sum_{k=1}^{K_l} \nabla_{\theta}\rho_{\theta}(X_{s_{k-1}})(X_{s_{k}} - X_{s_{k-1}}) +
\sum_{t=1}^T\nabla_{\theta}\log g_{\theta}(y_{t}|X_{t})\label{eqn:rewrite_discrete_test_function}.
\end{align}

For notational convenience, we define the $d\times 1$ vector of derivatives of $\rho_{\theta}$ as

\kappa_{\theta,i}(x)^* = \Big(\frac{\partial}{\partial \theta_i}[\rho_{\theta}(x)]^1,\dots,
\frac{\partial}{\partial \theta_i}[\rho_{\theta}(x)]^d\Big),

for any $(i,x)\in\{1,\dots,d_{\theta}\}\times\mathbb{R}^d$,

and the conditional likelihood given states $x_1,\dots,x_T \in \mathbb{R}^d$ as $\varphi_{\theta}(x_1,\dots,x_T) = \prod_{t = 1}^T g_{\theta}(y_t | x_t)$.
We now give the proof of Theorem <ref> followed by several technical lemmata
that are required to establish the theorem.

We consider the proof for any given component $i\in\{1,\dots,d_{\theta}\}$ and decompose
the error of the score function (<ref>) at level $l \in \mathbb{N}_0$ as

\begin{equation}\label{eq:theo_main_eq}
[S_l(\theta) - S(\theta)]^i = T_1 + T_2
\end{equation}

\begin{align*}
T_1 & = \frac{\mathbb{E}_{\theta}[\varphi_{\theta}(\widetilde{X}_1^{\star},\dots,\widetilde{X}_T^{\star})[G_{\theta}^l(\widetilde{X}_{0:T}^{\star})]^i]}
\mathbb{E}_{\theta}[\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})]
\mathbb{E}_{\theta}[\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})] -
\mathbb{E}_{\theta}[\varphi_{\theta}(\widetilde{X}_1^{\star},\dots,\widetilde{X}_T^{\star})]
\Big), \\
T_2 & = \frac{1}{\mathbb{E}_{\theta}[\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})]}
\Big(
\mathbb{E}_{\theta}[\varphi_{\theta}(\widetilde{X}_1^{\star},\dots,\widetilde{X}_T^{\star})[G_{\theta}^l(\widetilde{X}_{0:T}^{\star})]^i] - %\\ & &
\mathbb{E}_{\theta}[\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})[G_{\theta}(X^{\star})]^i]
\Big).
\end{align*}

Thus our objective is to provide bounds on the quantities $T_1$ and $T_2$ to conclude the proof.

For $T_1$, using Assumption <ref>, one has the upper-bound

T_1 \leq C \sum_{t=1}^T \mathbb{E}_{\theta}[\|\widetilde{X}_t^{\star}-X_t^{\star}\|_2],

then by using results on the convergence of Euler approximations (e.g. [Kloeden and Platen, 2013]), for $r>0$

\begin{equation}\label{eq:eul_conv_strong}
\mathbb{E}_{\theta}\big[\|\widetilde{X}^{\star}_t-X^{\star}_t\|_2^r\big]^{1/r} \leq C\Delta_l^{1/2}
\end{equation}

one has

\begin{equation}\label{eq:t1_bias_bd}
T_1 \leq C \Delta_l^{1/2}.
\end{equation}

Note that using standard results on weak errors for diffusions one can improve this upper-bound to $T_1 \leq C \Delta_l$.

For $T_2$, using Assumption <ref>, we have $T_2\leq C(T_3+T_4)$ where

\begin{align*}
T_3 & = \mathbb{E}_{\theta}\big[\{\varphi_{\theta}(\widetilde{X}_1^{\star},\dots,\widetilde{X}_T^{\star})-\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})\}
[G_{\theta}^l(\widetilde{X}_{0:T}^{\star})]^i\big], \\
T_4 & = \mathbb{E}_{\theta}\big[\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})\{
\end{align*}

For $T_3$, using Cauchy-Schwarz, we have the upper-bound

T_3 \leq \mathbb{E}_{\theta}\big[\{\varphi_{\theta}(\widetilde{X}_1^{\star},\dots,\widetilde{X}_T^{\star})-\varphi_{\theta}(X_1^{\star},\dots,X_T^{\star})\}^2\big]^{1/2}
\mathbb{E}_{\theta}\big[\{[G_{\theta}^l(\widetilde{X}_{0:T}^{\star})]^i\}^2\big]^{1/2}.

As the second term is bounded by $C$, we consider only the first. We have the upper-bound

\begin{equation}\label{eq:t3_bias_bd}
T_3 \leq C\sum_{t=1}^T \mathbb{E}_{\theta}\big[\|\widetilde{X}_t^{\star}-X_t^{\star}\|_2^2\big]^{1/2} \leq C\Delta_l^{1/2}.
\end{equation}

For $T_4$, noting that $\varphi_{\theta}$ is a bounded function under Assumption <ref>,
applying Lemma <ref> allows one to conclude that $T_4 \leq C\Delta_l^{1/2}$.
Therefore, using $T_2\leq C(T_3+T_4)$ along with (<ref>), we have

\begin{equation}\label{eq:t2_bias_bd}
T_2 \leq C \Delta_l^{1/2}.
\end{equation}

Combining (<ref>) and (<ref>) with (<ref>) allows one to conclude the proof.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta,i)\in \mathbb{N}\times[1,\infty)\times\Theta\times\{1,\dots,d_{\theta}\}$, there exists a constant $C<\infty$ such that for any $(l,x)\in\mathbb{N}_0\times\mathbb{R}^d$

\mathbb{E}_{\theta}\Big[\big|[G_{\theta}^l(\widetilde{X}_{0:T})]^i-[G_{\theta}(X)]^i\big|^r\Big]^{1/r} \leq C\Delta_l^{1/2},

with $\widetilde{X}_0 = X_0 = x$.

We have that

\begin{equation}\label{eq:main_lem_main_eq}
\mathbb{E}_{\theta}\Big[\big|[G_{\theta}^l(\widetilde{X}_{0:T})]^i-[G_{\theta}(X)]^i\big|^r\Big] \leq C(T_1+T_2)
\end{equation}

\begin{align*}
T_1 & = \mathbb{E}_{\theta}\Big[\big|[G_{\theta}^l(\widetilde{X}_{0:T})]^i-[G^l_{\theta}(X_{0:T})]^i\big|^r\Big], \\
T_2 & = \mathbb{E}_{\theta}\Big[\big|[G_{\theta}^l(X_{0:T})]^i-[G_{\theta}(X)]^i\big|^r\Big],
\end{align*}

where $X_{0:T} = (X_{s_k})_{k=0}^{K_l}$ are the states of the process $(X_t)_{t\geq 0}$ at the discretization times of the process $\widetilde{X}_{0:T}$.
From (<ref>), we have that $T_1 \leq C\sum_{j=3}^6 T_j$, where

\begin{align*}
T_3 & = \mathbb{E}_{\theta}\Big[\Big|\sum_{t=1}^T \Big\{\big[\nabla_{\theta} \log g_{\theta}(y_t|\widetilde{X}_t)\big]^i -
\big[\nabla_{\theta} \log g_{\theta}(y_t|X_t)\big]^i\Big\}\Big|^r\Big], \\
T_4 & = \Delta_l^r\mathbb{E}_{\theta}\Big[\Big|\sum_{k=1}^{K_l}\Big\{\big[\nabla_{\theta}\|b_{\theta}(\widetilde{X}_{s_{k-1}})\|_2^2\big]^i - \big[\nabla_{\theta}\|b_{\theta}(X_{s_{k-1}})\|_2^2\big]^i\Big\}\Big|^r\Big], \\
T_5 & = \mathbb{E}_{\theta}\Big[\Big| \sum_{k=1}^{K_l} \{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^* - \kappa_{\theta,i}(X_{s_{k-1}})^*\} [\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}]\Big|^r\Big], \\
T_6 & = \mathbb{E}_{\theta}\Big[\Big| \sum_{k=1}^{K_l} \kappa_{\theta,i}(X_{s_{k-1}})^* [(\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}) - (X_{s_k}-X_{s_{k-1}}) ]\Big|^r\Big].
\end{align*}

The term $T_3$ can be treated in almost the same manner as $T_1$ in the proof of Theorem <ref>, i.e.
using a similar argument to the proof of the bound on $T_1$ in Theorem <ref>, one can deduce that

\begin{equation}\label{eq:t3_bd_lem-b1}
T_3 \leq  C\Delta_l^{r/2}.
\end{equation}

For $T_4$, using the fact that $\partial/\partial\theta_i [b_{\theta}^2]^j\in \textrm{Lip}_{\|\cdot\|_2}(\mathbb{R}^d)$ for any $(i,j)\in\{1,\dots,d_{\theta}\}\times\{1,\dots,d\}$, we have by first applying Minkowski's inequality

T_4 \leq  C\Delta_l^r\Big(\sum_{k=1}^{K_l}\mathbb{E}_{\theta}\big[\|\widetilde{X}_{s_{k-1}}-X_{s_{k-1}}\|_2^r\big]^{1/r}\Big)^r.

Then using (<ref>), it follows that

\begin{equation}\label{eq:t4_bd_lem-b1}
T_4 \leq  C\Delta_l^{r/2}.
\end{equation}

The terms $T_5$ and $T_6$ are bounded in Lemmata <ref>-<ref>, so combining (<ref>), (<ref>) and the afore-mentioned lemmata with $T_1 \leq C\sum_{j=3}^6 T_j$ yields

T_1 \leq  C\Delta_l^{r/2}.

By Lemma <ref>, $T_2 \leq  C\Delta_l^{r/2}$ and thus by (<ref>) the proof is concluded.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta)\in \mathbb{N}\times[1,\infty)\times\Theta$, there exists a constant $C<\infty$ such that for any $(l,x)\in\mathbb{N}_0\times\mathbb{R}^d$

\mathbb{E}_{\theta}\Big[ \big\| G_{\theta}^l(\widetilde{X}_{0:T})-G_{\theta}(X) \big\|_2^r \Big]^{1/r} \leq C\Delta_l^{1/2}.

By Minkowski's inequality

\mathbb{E}_{\theta}\Big[ \big\| G_{\theta}^l(\widetilde{X}_{0:T})-G_{\theta}(X) \big\|_2^r \Big]^{1/r} \leq
\bigg(\sum_{i=1}^{d_{\theta}} \mathbb{E}_{\theta}\Big[ \big| [G_{\theta}^l(\widetilde{X}_{0:T})]^i - [G_{\theta}(X)]^i \big|^r \Big]^{2/r} \bigg)^{1/2}

so the proof follows by Lemma <ref>.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta,i)\in \mathbb{N}\times[1,\infty)\times\Theta\times\{1,\dots,d_{\theta}\}$, there exists a constant $C<\infty$ such that for any $(l,x)\in\mathbb{N}_0\times\mathbb{R}^d$

\mathbb{E}_{\theta}\bigg[\Big|
\sum_{k=1}^{K_l}
\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}
[\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}]\Big|^r\bigg] \leq C\Delta_l^{r/2},

with $\widetilde{X}_0 = X_0 = x$.

We have the decomposition

\sum_{k=1}^{K_l}
\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}
[\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}] = M_{K_l} + R_{K_l},

\begin{align*}
M_{K_l} & = \sum_{k=1}^{K_l}\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}\sigma(\widetilde{X}_{s_{k-1}})[W_{s_k}-W_{s_{k-1}}],\\
R_{K_l} & = \Delta_l\sum_{k=1}^{K_l}
\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}a_{\theta}(\widetilde{X}_{s_{k-1}}).
\end{align*}

Thus by the $C_r-$inequality,

\begin{equation}\label{eq:diff_mart_prf_lem6}
\mathbb{E}_{\theta}\bigg[\Big|
\sum_{k=1}^{K_l}
\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}
[\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}] \Big|^r\bigg] \leq C\big( \mathbb{E}_{\theta}[|M_{K_l}|^r]+\mathbb{E}_{\theta}[|R_{K_l}|^r]\big).
\end{equation}

We will bound the two terms on the R.H.S. of (<ref>) individually.

Bound for $\mathbb{E}_{\theta}[|M_{K_l}|^r]$. If we define $M_{u} =  \sum_{k=0}^{u-1}\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}\sigma(\widetilde{X}_{s_{k-1}})[W_{s_k}-W_{s_{k-1}}]$ for any $u \in \mathbb{N}_0$,
then $(M_u,\mathscr{F}_{u\Delta_l})_{u\in \mathbb{N}_0}$ is a martingale. It follows from this fact and from an application of the Burkholder-Gundy-Davis (BGD) inequality that

\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C
\mathbb{E}_{\theta}\bigg[\Big|\sum_{k=1}^{K_l}
\kappa_{\theta,i}(X_{s_{k-1}})^*\}\sigma(\widetilde{X}_{s_{k-1}})[W_{s_k}-W_{s_{k-1}}])^2
\Big|^{r/2}\bigg],

from which Minkowski's inequality yields

\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C
\bigg(
\sum_{k=1}^{K_l}  \mathbb{E}_{\theta}\Big[\big|
\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}\sigma(\widetilde{X}_{s_{k-1}})[W_{s_k}-W_{s_{k-1}}]\big|^r\Big]^{2/r}
\bigg)^{r/2}.

Using the $C_r-$inequality $d^2$ times, we obtain the bound

\begin{multline*}
\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq \\
C\Bigg( \sum_{k=1}^{K_l} \bigg( \sum_{(m,j)\in\{1,\dots,d\}^2}
\mathbb{E}_{\theta}\Big[\big|\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}}) - \kappa_{\theta,i}(X_{s_{k-1}})\}^m
\sigma(\widetilde{X}_{s_{k-1}})^{m,j}[W_{s_k}-W_{s_{k-1}}]^j\big|^r\Big]
\bigg)^{2/r}\Bigg)^{r/2}.
\end{multline*}

Using the fact that $\sigma^{m,j}\in\mathcal{B}_b(\mathbb{R}^d)$ along with the Cauchy-Schwarz inequality yields

\begin{multline}
\label{eq:diff_mart_prf_lem1}
\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq \\
C\Bigg( \sum_{k=1}^{K_l}\bigg(\sum_{(m,j)\in\{1,\dots,d\}^2}
\mathbb{E}_{\theta}\Big[ \big| \{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})-
\kappa_{\theta,i}(X_{s_{k-1}})\}^m \big|^{2r} \Big]^{1/2}
\mathbb{E}_{\theta}\Big[ \big|[W_{s_k}-W_{s_{k-1}}]^j\big|^{2r} \Big]^{1/2}
\bigg)^{2/r} \Bigg)^{r/2}.
\end{multline}

Since it holds that $[\kappa_{\theta,i}]^m\in\textrm{Lip}_{\|\cdot\|}(\mathbb{R}^d)$, it follows from the same type of inequality as (<ref>) that

\begin{equation}\label{eq:diff_mart_prf_lem2}
\mathbb{E}_{\theta}\Big[ \big| \{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}}) - \kappa_{\theta,i}(X_{s_{k-1}})\}^m \big|^{2r} \Big]^{1/2} \leq C\Delta_l^{r/2},
\end{equation}

and by standard properties of Brownian motion, we obtain

\begin{equation}\label{eq:diff_mart_prf_lem3}
\mathbb{E}_{\theta}\Big[ \big| [W_{s_k}-W_{s_{k-1}}
]^j \big|^{2r} \Big]^{1/2}
\leq C\Delta_l^{r/2}.
\end{equation}

Combining (<ref>) with (<ref>) and (<ref>) yields the upper-bound

\begin{equation}\label{eq:diff_mart_prf_lem4}
\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C\Delta_l^{r/2}.
\end{equation}

Bound for $\mathbb{E}_{\theta}[|R_{K_l}|^r]$. We have the upper-bound by Minkowski's inequality

\mathbb{E}_{\theta}[|R_{K_l}|^r] \leq \Delta_l^r\bigg(\sum_{k=1}^{K_l}
\mathbb{E}_{\theta}\Big[\big|\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})^*-
\kappa_{\theta,i}(X_{s_{k-1}})^*\}a_{\theta}(\widetilde{X}_{s_{k-1}}) \big|^r \Big]^{1/r}
\bigg)^r.

Then applying the $C_r-$inequality $d$ times and using the assumption that $a_{\theta}^j\in\mathcal{B}_b(\mathbb{R}^d)$ we have the upper-bound

\mathbb{E}_{\theta}[|R_{K_l}|^r] \leq C\Delta_l^r\bigg(\sum_{k=1}^{K_l}
\sum_{j=1}^d
\mathbb{E}_{\theta}\Big[\big|\{\kappa_{\theta,i}(\widetilde{X}_{s_{k-1}})-
\kappa_{\theta,i}(X_{s_{k-1}})\}^j \big|^{r} \Big]^{1/r} \bigg)^r.

Using the same argument to obtain (<ref>), we have

\begin{equation}\label{eq:diff_mart_prf_lem5}
\mathbb{E}_{\theta}[|R_{K_l}|^r] \leq C\Delta_l^r\Delta_l^{-r/2} =C\Delta_l^{r/2}.
\end{equation}

Combining (<ref>) and (<ref>) with (<ref>) allows us to conclude.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta,i)\in \mathbb{N}\times[1,\infty)\times\Theta\times\{1,\dots,d_{\theta}\}$, there exists a constant $C<\infty$ such that for any $(l,x)\in\mathbb{N}_0\times\mathbb{R}^d$

\mathbb{E}_{\theta}\bigg[\Big| \sum_{k=1}^{K_l} \kappa_{\theta,i}(X_{s_{k-1}})^* \big[(\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}})- (X_{s_k}-X_{s_{k-1}}) \big]\Big|^r\bigg] \leq C\Delta_l^{r/2},

with $\widetilde{X}_0 = X_0 = x$.

We have the decomposition

\sum_{k=1}^{K_l} \kappa_{\theta,i}(X_{s_{k-1}})^* \big[(\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}})- (X_{s_k}-X_{s_{k-1}}) \big] = M_{K_l} + R_{K_l},

\begin{align*}
M_{K_l} & = \sum_{k=1}^{K_l}\sum_{(m,j)\in\{1,\dots,d\}^2}\kappa_{\theta,i}(X_{s_{k-1}})^m \int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j} dW_s^j,\\
R_{K_l} & = \sum_{k=1}^{K_l}\sum_{j=1}^d \kappa_{\theta,i}(X_{s_{k-1}})^j \int_{s_{k-1}}^{s_k} \big[a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_s) \big]^j ds.
\end{align*}

Thus by the $C_r-$inequality,

\begin{equation}\label{eq:diff_mart_prf_lem7}
\mathbb{E}_{\theta}\bigg[ \Big|
\sum_{k=1}^{K_l} \kappa_{\theta,i}(X_{s_{k-1}})^*
\big[ (\widetilde{X}_{s_k}-\widetilde{X}_{s_{k-1}}) - (X_{s_k}-X_{s_{k-1}}) \big]
\Big|^r \bigg] \leq C\big( \mathbb{E}_{\theta}[|M_{K_l}|^r]+\mathbb{E}_{\theta}[|R_{K_l}|^r] \big).
\end{equation}

We will bound the two terms on the R.H.S. of (<ref>) individually.

Bound for $\mathbb{E}_{\theta}[|M_{K_l}|^r]$. By applying the $C_r-$inequality $d^2$ times we have the upper-bound

\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C \sum_{(m,j)\in\{1,\dots,d\}^2} \mathbb{E}_{\theta}\bigg[\Big|
\sum_{k=1}^{K_l} \kappa_{\theta,i}(X_{s_{k-1}})^m \int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j \Big|^r\bigg].

For any $(m,j) \in \{1,\dots,d\}^2$, we define $M_u^{m,j} = \sum_{k=0}^{u-1} \kappa_{\theta,i}(X_{s_{k-1}})^m \int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j$ for $u \in \mathbb{N}_0$.
As $(M_u^{m,j},\mathscr{F}_{u\Delta_l})_{u\in \mathbb{N}_0}$ is a martingale, applying the BGD inequality yields

\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C \sum_{(m,j)\in\{1,\dots,d\}^2} \mathbb{E}_{\theta} \bigg[ \Big| \sum_{k=1}^{K_l} \Big(\kappa_{\theta,i}(X_{s_{k-1}})^m \int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j\Big)^2 \Big|^{r/2}\bigg].

Applying Minkowski's inequality and using the fact that $[\kappa_{\theta,i}]^m\in\mathcal{B}_b(\mathbb{R}^d)$, we obtain

\begin{equation}\label{eq:diff_mart_prf_lem8}
\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C \sum_{(m,j)\in\{1,\dots,d\}^2} \Bigg( \sum_{k=1}^{K_l}
\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j
\Big|^{r}\bigg]^{2/r} \Bigg)^{r/2}.
\end{equation}

Now we deal with the expectation on the R.H.S. of (<ref>). Using the martingale property of the stochastic integral, it follows from applying the BGD inequality again that

\begin{align*}
\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j
\Big|^{r}\bigg] & \leq
C\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k} \big\{ [\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j} \big\}^2 ds
\Big|^{r/2}\bigg] \\
& \leq C\Delta_l^{r/2-1}\mathbb{E}_{\theta}\Big[\int_{s_{k-1}}^{s_k} \big|[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j} \big|^r ds\Big],
\end{align*}

where we have used Jensen's inequality in the second line. Using the $C_r-$inequality, we have

\begin{multline*}
\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j
\Big|^{r}\bigg] \leq \\
\Big\{ \mathbb{E}_{\theta}\Big[ \big| [\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_{{s_{k-1}}})]^{m,j} \big|^r \Big]
+ \mathbb{E}_{\theta}\Big[ \big| [\sigma(X_{{s_{k-1}}})-\sigma(X_s)]^{m,j} \big|^r \Big]
\Big\}ds.
\end{multline*}

Using the fact that $\sigma^{m,j}\in\textrm{Lip}_{\|\cdot\|_2}(\mathbb{R}^d)$, we then obtain

\begin{multline*}
\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k}[\sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s)]^{m,j}dW_s^j
\Big|^{r}\bigg] \leq \\
C\Delta_l^{r/2-1}\int_{s_{k-1}}^{s_k} \Big\{\mathbb{E}_{\theta}
\big[ \|
\widetilde{X}_{s_{k-1}}-X_{{s_{k-1}}}\|_2^r \big] + \mathbb{E}_{\theta}\big[ \|X_{{s_{k-1}}}-X_s\|_2^r \big] \Big\}ds.
\end{multline*}

By the property (<ref>) for $r>0$ (e.g. [Ikeda and Watanabe, 2014]), it follows that

\begin{equation}\label{eq:diff_cont_l_r}
\sup_{(t,s)\in[0,T]^2}\mathbb{E}_{\theta} \big[ \|X_t-X_s\|_2^r \big] \leq C|t-s|^{r/2},
\end{equation}

hence we obtain the upper-bound

\begin{equation}\label{eq:diff_mart_prf_lem9}
\mathbb{E}_{\theta}\bigg[\Big|\int_{s_{k-1}}^{s_k} \big[ \sigma(\widetilde{X}_{s_{k-1}})-\sigma(X_s) \big]^{m,j} dW_s^j
\Big|^{r}\bigg] \leq C\Delta_l^{r/2-1}\Delta_l^{1+r/2}= C\Delta_l^r.
\end{equation}

Combining (<ref>) with (<ref>) gives

\begin{equation}\label{eq:diff_mart_prf_lem10}
\mathbb{E}_{\theta}[|M_{K_l}|^r] \leq C\Delta_l^{r/2}.
\end{equation}

Bound for $\mathbb{E}_{\theta}[|R_{K_l}|^r]$. Using Minkowski's inequality followed by Jensen's inequality,
we obtain the following bounds

\begin{align*}
\mathbb{E}_{\theta}[|R_{K_l}|^r] & \leq
\Bigg( \sum_{k=1}^{K_l} \sum_{j=1}^d
\mathbb{E}_{\theta} \bigg[ \Big| \kappa_{\theta,i}(X_{s_{k-1}})^j
\int_{s_{k-1}}^{s_k}[a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_s)]^jds \Big|^r \bigg]^{1/r} \Bigg)^r\\
& \leq
\Bigg(\sum_{k=1}^{K_l}\sum_{j=1}^d \Delta_l^{1-1/r} \mathbb{E}_{\theta}\Big[
\big| \kappa_{\theta,i}(X_{s_{k-1}})^j \big|^r
\int_{s_{k-1}}^{s_k}\big| [a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_s)]^j \big|^rds
\Big]^{1/r} \Bigg)^r.
\end{align*}

Using the assumption $[\kappa_{\theta,i}]^j\in\mathcal{B}_b(\mathbb{R}^d)$, the decomposition

[a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_s)]^j =
[a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_{s_{k-1}})]^j + [a_{\theta}(X_{s_{k-1}})-a_{\theta}(X_s)]^j,

and the $C_r-$inequality, we have

\begin{align*}
& \mathbb{E}_{\theta}[|R_{K_l}|^r] \leq \\
& \Bigg(\sum_{k=1}^{K_l} \sum_{j=1}^d C\Delta_l^{1-1/r} \bigg( \int_{s_{k-1}}^{s_k} \mathbb{E}_{\theta}\Big[ \big| [a_{\theta}(\widetilde{X}_{s_{k-1}})-a_{\theta}(X_{s_{k-1}})]^j \big|^r \Big] ds + \int_{s_{k-1}}^{s_k} \mathbb{E}_{\theta}\Big[ \big| [a_{\theta}(X_{s_{k-1}})-a_{\theta}(X_s)]^j \big|^r \Big] ds \bigg)^{1/r} \Bigg)^r.
\end{align*}

Since we have assumed that $a_{\theta}^j\in\textrm{Lip}_{\|\cdot\|_2}(\mathbb{R}^d)$, one has

\begin{equation*}
\mathbb{E}_{\theta}[|R_{K_l}|^r] \leq \Bigg(\sum_{k=1}^{K_l}\sum_{j=1}^d C\Delta_l^{1-1/r} \bigg(
\int_{s_{k-1}}^{s_k} \mathbb{E}_{\theta} \big[ \|\widetilde{X}_{s_{k-1}}-X_{s_{k-1}}\|_2^r \big] ds + \int_{s_{k-1}}^{s_k} \mathbb{E}_{\theta} \big[ \|X_{s_{k-1}}-X_s\|_2^r \big] ds \bigg)^{1/r} \Bigg)^r.
\end{equation*}

Using (<ref>)  and (<ref>), we therefore obtain

\begin{equation}\label{eq:diff_mart_prf_lem11}
\mathbb{E}_{\theta}[|R_{K_l}|^r] \leq C\Delta_l^{r/2}.
\end{equation}

Combining (<ref>) and (<ref>) with (<ref>) allows us to conclude.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta,i)\in \mathbb{N}\times[1,\infty)\times\Theta\times\{1,\dots,d_{\theta}\}$, there exists a constant $C<\infty$ such that for any $(l,x)\in\mathbb{N}_0\times\mathbb{R}^d$

\mathbb{E}_{\theta}\Big[\big|[G_{\theta}^l(X_{0:T})]^i - [G_{\theta}(X)]^i\big|^r \Big] \leq C\Delta_l^{r/2},

with $\widetilde{X}_0 = X_0 = x$.

We have the decomposition

\big| [G_{\theta}^l(X_{0:T})]^i - [G_{\theta}(X)]^i \big| = R_{K_l}^{(1)} + M_{K_l}
+ R_{K_l}^{(2)},

\begin{align*}
R_{K_l}^{(1)} & = \sum_{k=1}^{K_l}\int_{s_{k-1}}^{s_k}
\Big([\nabla_{\theta}\|b_{\theta}(X_{s_{k-1}})\|_2^2]^i - [\nabla_{\theta}\|b_{\theta}(X_{s})\|_2^2]^i\Big)
ds, \\
M_{K_l} & = \sum_{k=1}^{K_l}\sum_{(m,j)\in\{1,\dots,d\}^2}
\int_{s_{k-1}}^{s_k}[\kappa_{\theta,i}(X_{s_{k-1}})^m\sigma(\widetilde{X}_{s_{k-1}})-
\kappa_{\theta,i}(X_{s})^m
\sigma(X_s)]^{m,j}dW_s^j, \\
R_{K_l}^{(2)} & = \sum_{k=1}^{K_l}\sum_{j=1}^d\int_{s_{k-1}}^{s_k}
\Big(\kappa_{\theta,i}(X_{s_{k-1}})^ja_{\theta}(X_{s_{k-1}})^j -
\kappa_{\theta,i}(X_{s})^ja_{\theta}(X_{s})^j\Big)ds.
\end{align*}

Thus, by the $C_r-$inequality, we have

\mathbb{E}_{\theta}\Big[ \big| [G_{\theta}^l(X_{0:T})]^i - [G_{\theta}(X)]^i \big|^r \Big] \leq
C\Big( \mathbb{E}_{\theta}\big[ |M_{K_l}|^r \big]+\mathbb{E}_{\theta}\big[ \big| R_{K_l}^{(1)} \big|^r \big]+\mathbb{E}_{\theta}\big[ \big| R_{K_l}^{(2)}\big|^r \big]
\Big).

In order to prove that $\mathbb{E}_{\theta}[|M_{K_l}|^r]\leq C\Delta_l^{r/2}$, one can rely on very similar arguments to (<ref>) in the proof of Lemma <ref>. Likewise, for $m\in\{1,2\}$ one can prove that $\mathbb{E}_{\theta}[|R_{K_l}^{(m)}|^r]\leq C\Delta_l^{r/2}$ using similar arguments to (<ref>) in the proof of Lemma <ref>.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta)\in \mathbb{N}\times[1,\infty)\times\Theta$, there exists a constant $C<\infty$ such that for any $(x,x_{\star})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$

\mathbb{E}_{\theta}\big[\| G_{\theta}(X)-G_{\theta}(X^{\star}) \|_2^r \big]^{1/r} \leq C\|x-x_{\star}\|_2,

with $X_0 = x$ and $X_0^{\star} = x_{\star}$.

This proof follows a similar type of arguments to those considered in the proofs of Lemmata <ref>-<ref>. The main difference is that one must use the following result (which can be deduced from <cit.> and Grönwall's inequality)

\begin{equation*}%\label{eq:cont_exp_diff}
\sup_{t\in[0,T]}\mathbb{E}_{\theta}[\|X_t-X_t^{\star}\|_2^{2r}]^{\frac{1}{2r}} \leq C\|x-x_{\star}\|_2
\end{equation*}

instead of using Euler convergence. Given that the proofs of Lemmata <ref>-<ref> are repetitive, these arguments are omitted.

Under Assumptions <ref> and <ref>, for any $(T,r,\theta)\in \mathbb{N}\times[1,\infty)\times\Theta$, there exists a constant $C<\infty$ such that for any $(l,x,x_{\star})\in\mathbb{N}\times\mathbb{R}^{d}\times\mathbb{R}^{d}$

\mathbb{E}_{\theta}\big[ \|G_{\theta}^l(\widetilde{X}_{0:T})- G_{\theta}^{l-1}(\widehat{X}^{\star}_{0:T}) \|_2^r \big]^{1/r}
\leq C\big(\Delta_l^{1/2} + \|x - x_{\star}\|_2\big),

with $\widetilde{X}_0 = x$ and $\widetilde{X}_0^{\star} = x_{\star}$,
where $\widehat{X}^{\star}_{0:T} = (\widetilde{X}^{\star}_{s_{2k}})_{k=0}^{K_{l-1}}$ are
the sequence of states in $\widetilde{X}^{\star}_{0:T}=(\widetilde{X}^{\star}_{s_k})_{k=0}^{K_l}$ at the time steps corresponding to level $l-1$.

The expectation in the statement of the lemma is upper-bounded by $\sum_{j=1}^3 T_j$, where

\begin{align*}
T_1 & = \mathbb{E}_{\theta}\big[\|G_{\theta}^l(\widetilde{X}_{0:T})-G_{\theta}(X)\|_2^r\big]^{1/r}, \\
T_2 & = \mathbb{E}_{\theta}\big[\|G_{\theta}(X) - G_{\theta}(X^{\star})\|_2^r\big]^{1/r}, \\
T_3 & = \mathbb{E}_{\theta}\big[\|G_{\theta}^{l-1}
\end{align*}

For $T_1$ and $T_3$, one can apply Corollary <ref>; for $T_2$, we use Lemma <ref>.
This allows us to conclude the result.

We now introduce two additional functions which will be useful in the following section.
For $(t,l)\in\{1,\dots,T\}\times\mathbb{N}_0$, we define
$G_{\cdot,t}^l:\Theta\times(\mathbb{R}^d)^{2^l t+1}\rightarrow\mathbb{R}^{d_{\theta}}$ as

G_{\theta,t}^l(X_{0:t}) =  -\frac{1}{2}
\sum_{k=1}^{2^lt} \nabla_{\theta}\|b_{\theta}(X_{s_{k-1}})\|_2^2\Delta_l
+ \sum_{k=1}^{2^lt} \nabla_{\theta}\rho_{\theta}(X_{s_{k-1}})
(X_{s_k}-X_{s_{k-1}})+\sum_{p=1}^t \nabla_{\theta} \log g_{\theta}(y_p|X_p),

and $G_{\cdot,t-1:t}^l:\Theta \times \mathbb{R}^d \times (\mathbb{R}^d)^{2^l} \rightarrow \mathbb{R}^{d_{\theta}}$ as

\begin{align*}
&G_{\theta,t-1:t}^l(X_{t-1},X_{t-1+\Delta_l:t}) \\
&=  -\frac{1}{2}
\sum_{k=0}^{2^l-1} \nabla_{\theta}\|b_{\theta}(X_{t-1+k\Delta_l})\|_2^2\Delta_l
+ \sum_{k=0}^{2^l-1} \nabla_{\theta}\rho_{\theta}(X_{t-1+k\Delta_l})
(X_{t-1+(k+1)\Delta_l}-X_{t-1+k\Delta_l}) + \nabla_{\theta} \log g_{\theta}(y_t|X_t).
\end{align*}

From (<ref>), we note that $G_{\theta,T}^l(X_{0:T})=G_{\theta}^l(X_{0:T})$.
The following remarks will facilitate our proofs in the next section.

One can easily extend Lemma <ref> as follows. Under Assumptions <ref> and <ref>,
for any $(t,r,\theta)\in \{1,\dots,T\}\times[1,\infty)\times\Theta$, there exists a constant $C<\infty$ such that for any

\mathbb{E}_{\theta}\Big[ \big\|G_{\theta,t}^l(\widetilde{X}_{0:t})-
\big\|_2^r \Big]^{1/r} \leq C\big(\Delta_l^{1/2} + \|x-x_{\star}\|_2\big),

with $\widetilde{X}_0 = x$ and $\widetilde{X}_0^{\star} = x_{\star}$,
where $\widehat{X}^{\star}_{0:t} = (\widetilde{X}^{\star}_{s_{2k}})_{k=0}^{2^{l-1}t}$.

One can also deduce the following result. Under Assumptions <ref> and <ref>,
for any $(t,r,\theta)\in \{1,\dots,T\}\times[1,\infty)\times\Theta$, there exists a constant $C<\infty$
such that for any $(l,x,x_{\star})\in\mathbb{N}\times\mathbb{R}^{d}\times\mathbb{R}^{d}$

\mathbb{E}_{\theta}\Big[ \big\|G_{\theta,t-1:t}^l(x,\widetilde{X}_{t-1+\Delta_l:t})-
\big\|_2^r \Big| \widetilde{X}_{t-1}=x,\widetilde{X}_{t-1}^{\star}=x_{\star} \Big]^{1/r}
\leq  C\big(\Delta_l^{1/2} + \|x-x_{\star}\|_2 \big),

where $\widehat{X}_{t-1+\Delta_{l-1}:t}^{\star} = (\widetilde{X}^{\star}_{t-1+\Delta_{l-1}})_{k=1}^{2^{l-1}}$.

§.§ Results on coupled conditional particle filters

We begin with some definitions associated to Algorithm <ref>.

For any $(i,t,s)\in\{1,\dots,N\}\times\mathbb{N}_0\times\{l-1,l\}$,
we will write $A_t^{s}(i) = A_t^{s,i}$ and $\bar{A}_t^{s}(i) = \bar{A}_t^{s,i}$.
Using this notation, for any $(t,l)\in\{0,\dots,T-1\}\times\mathbb{N}$, we define

\begin{align*}
& \mathsf{S}_t^l=\{i\in\{1,\ldots, N\}: A_{t}^l(i)=A_{t}^{l-1}(i),A_{t-1}^l\circ A_{t}^l(i)=A_{t-1}^{l-1}\circ A_{t}^{l-1}(i),\dots,%\\
A_{0}^l\circ\cdots\circ A_{t}^l(i)=A_{0}^{l-1}\circ\cdots\circ A_{t}^{l-1}(i)\},\\
& \bar{\mathsf{S}}_t^l=\{i\in\{1,\ldots, N\}: \bar{A}_{t}^l(i)=\bar{A}_{t}^{l-1}(i),\bar{A}_{t-1}^l\circ \bar{A}_{t}^l(i)=\bar{A}_{t-1}^{l-1}\circ \bar{A}_{t}^{l-1}(i),\dots,%\\
\bar{A}_{0}^l\circ\cdots\circ \bar{A}_{t}^l(i)=\bar{A}_{0}^{l-1}\circ\cdots\circ \bar{A}_{t}^{l-1}(i)\}.
\end{align*}

For $(l,\beta,C)\in\mathbb{N}\times\mathbb{R}^+\times\mathbb{R}^+$, we also introduce the following sets

\begin{align*}
\mathsf{B}^l_{\beta,C} & = \{(x_{0:T},\bar{x}_{0:T})\in\mathsf{X}^l\times\mathsf{X}^{l-1}:%(\mathbb{R}^d)^n\times (\mathbb{R}^d)^n:
\|x_t-\bar{x}_t\|_2 \leq C\Delta_l^{\beta}, t\in\{1,\dots,T\}\},\\
\mathsf{G}^l_{\beta,C} & = \{(x_{0:T},\bar{x}_{0:T})\in\mathsf{X}^l\times\mathsf{X}^{l-1}: