###### Abstract

We consider the problem of parameter estimation for a class of continuous-time
state space models (SSMs). In particular, we explore the case of a partially
observed diffusion, with data also arriving according to a diffusion process.
Based upon a standard identity of the score function, we consider two particle
filter based methodologies to estimate the score function. Both methods rely
on an online estimation algorithm for the score function, as described, e.g.,
in [12], of $\mathcal{O}(N^{2})$ cost, with $N\in\mathbb{N}$ the number of
particles. The first approach employs a simple Euler discretization and
standard particle smoothers and is of cost
$\mathcal{O}(N^{2}+N\Delta_{l}^{-1})$ per unit time, where
$\Delta_{l}=2^{-l}$, $l\in\mathbb{N}_{0}$, is the time-discretization step.
The second approach is new and based upon a novel diffusion bridge
construction. It yields a new backward type Feynman-Kac formula in continuous-
time for the score function and is presented along with a particle method for
its approximation. Considering a time-discretization, the cost is
$\mathcal{O}(N^{2}\Delta_{l}^{-1})$ per unit time. To improve computational
costs, we then consider multilevel methodologies for the score function. We
illustrate our parameter estimation method via stochastic gradient approaches
in several numerical examples.
Keywords: Score Function, Parameter Estimation, Particle Filter, Diffusion
Bridges.
AMS subject classifications: 65C05, 65C35, 60G35, 60J60, 60J65, 60H10, 60H35,
91G60

Score-Based Parameter Estimation for a Class of
Continuous-Time State Space Models

BY ALEXANDROS BESKOS1, DAN CRISAN2, AJAY JASRA3, NIKOLAS KANTAS2 & HAMZA
RUZAYQAT3

1Department of Statistical Science, University College London, London, WC1E
6BT, UK. E-Mail<EMAIL_ADDRESS>
2Department of Mathematics, Imperial College London, London, SW7 2AZ, UK.
E-Mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
3Computer, Electrical and Mathematical Sciences and Engineering Division, King
Abdullah University of Science and Technology, Thuwal, 23955, KSA. E-Mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

We consider the problem of parameter estimation for continuous-time SSMs.
These are models comprising stochastic differential equations (SDEs)
describing a hidden dynamic state and their observations. Such models are
ubiquitous in a large number of practical applications in science,
engineering, finance, and economics, see [10] for an overview. Inference in
SSMs, also known as hidden Markov models, hinges upon computing conditional
probability distributions of the dynamic hidden state given the acquired
observations and unknown static parameters. This is referred to as the
stochastic filtering problem, which is in general intractable, but reliable
numerical approximations are routinely available [1, 10]. The problem of
inferring the unknown parameters is more challenging. In this paper we focus
on maximum likelihood inference and gradient methods that are performed in an
online manner. The offline or iterative case given a fixed batch of
observations can also be treated using our proposed methods. The approach
considered in this article is to make use of the gradient of the log-
likelihood, commonly referred to as the _score function_ , within a stochastic
gradient algorithm (see e.g. [12, 13]). Intrinsically, there are several
challenges arising with such an approach. Firstly, when one adopts a
continuous-time model and assumes access to arbitrarily high frequency
observations, one does not observe in practice truly in continuous-time,
therefore some sort of time-discretization is required. Secondly, in both
discrete-time and continuous-time formulations, there are very few cases when
the score function is analytically available. Both these issues imply that
numerical approximations are required.

We consider two different approaches for the numerical approximation of the
score function. The first is to simply time-discretize a representation of the
score function and then apply discrete-time numerical approximation schemes
[12, 13, 22]. The second, is to develop a numerical approximation scheme
directly on continuous-time path-space to estimate the score function and then
(necessarily) discretize the algorithm in time. The order of designing the
estimation method and time-discretization can be rather important. Often the
second approach is preferable in terms of both performance and robustness as
the discretization mesh vanishes, see e.g. [24, 23]. We then use the score
estimate for implementing recursive maximum likelihood, where the parameters
are updated at unit time intervals ([20, 26]). The particular choice of time
interval length is without loss of generality, and allows the score to
accumulate sufficient information from the observations before updating the
parameters.

In the first approach, we consider a well-known expression for the score
function, for instance as given in [7]. Given this formula, one can produce an
Euler discretized version of the score and work in discrete-time but with high
frequency. The score is an expectation of an additive functional of the hidden
state path conditioned on the available observations, which is commonly
referred to as the smoothing distribution. In this context, many well-known
particle smoothing schemes can now be adopted, such as the ones described in
[12, 22]. These latter approaches are simulation-based schemes whose
convergence is based upon the number of samples $N\in\mathbb{N}$. We prove
some technical results for the discretized problem, which together with the
work in [12] allow us conjecture that to obtain a mean square error (MSE) of
$\mathcal{O}(\epsilon^{2})$, for given $\epsilon>0$, we require a
computational cost of $\mathcal{O}(\epsilon^{-4})$ per unit time. The latter
derives from an algorithmic cost of $\mathcal{O}(N^{2}+N\Delta_{l}^{-1})$ per
unit time, where $\Delta_{l}=2^{-l}$ is the time-discretization. As we explain
later in the article, this corresponds to a best case scenario, due to the
intrinsic nature of the algorithm. In particular, we start with a continuous-
time formula and time-discretize it, but the deduced numerical algorithm can
be problematic in terms of computational complexity as $l\rightarrow\infty$.
Whilst in some scenarios one does not observe any issues, examples can be
found where the variance of the method can explode as $l$ grows ([35]), thus
putting into question the validity of the conjecture on the cost to achieve an
MSE of $\mathcal{O}(\epsilon^{2})$.

This motivates the introduction of our second approach, where we build upon a
change of measure technique proposed in [28]. This latter approach has so far
been used in very different contexts than the present paper, namely related to
discretely observed diffusions for Bayesian inference and Markov chain Monte
Carlo (MCMC) [33, 34] or smoothing for potentially non-linear observation
functions [21]. Our approach is a data augmentation scheme, whereby, at unit
time intervals, the end points for the hidden state are sampled and the path
is connected using diffusion bridges. Then, starting again from the formula
for the score function in [7], we will use this change of measure associated
to a diffusion bridge and its driving Brownian noise; see also [35] where a
related approach is used for a different class of problems. Based upon this
change of measure we develop a new backward type Feynman-Kac formula in
continuous-time. This new formula facilitates an adaptation of the method in
[12] in true continuous-time, albeit one cannot apply it in practice. We time-
discretize the algorithm and conjecture that to obtain a MSE of
$\mathcal{O}(\epsilon^{2})$, for given $\epsilon>0$, we need a cost of
$\mathcal{O}(\epsilon^{-6})$ per unit time, which derives from an algorithmic
cost of $\mathcal{O}(N^{2}\Delta_{l}^{-1})$ per unit time. We note however,
that this computational complexity will not explode with increasing $l$ as may
be the case in the first approach. To improve the cost required for a given
MSE, we develop a novel multilevel Monte Carlo extension that can, in some
cases, achieve a MSE of $\mathcal{O}(\epsilon^{2})$ at a cost per unit time of
$\mathcal{O}(\epsilon^{-4})$. We remark that although our MSE-cost statements
are based upon conjectures, they are verified numerically. Direct proofs of
these require substantial technical results that will be the topic of future
work.

###  1.1 Contributions and Organization

We conclude this introduction by emphasizing that our contributions are aimed
to deal with both continuous-time observations and hidden states. As mentioned
earlier this poses very particular challenges relative to earlier works that
deal with filtering and smoothing when discrete-time observations/models are
used as in [6, 14, 15, 21, 34, 27]. None of these works look at continuous-
time observations. Similarly online likelihood estimation of the parameters
using the score function is considered in [14, 15] only for the case of
discrete-time observations of hidden diffusions.

The contributions of this paper can be summarized as follows:

  * •

We investigate the efficiency and accuracy of two fundamentally different
numerical approximations of the score function on its own and when used for
the purpose of recursive maximum likelihood. Both methods rely on fairly
standard tools such as changes of measure, particle smoothing and Euler time-
discretization.

  * •

We provide a detailed discussion on the computational complexity of each
method. We illustrate the performance and computational cost for several
models in numerical examples that consider estimation of the score function
and parameter estimation.

  * •

The second approach is a novel method that operates directly on the path-
space. The approach improves performance and is robust to arbitrarily small
time-discretization at the expense of additional computational cost. The
latter is reduced using a new Multilevel Particle Filter; see [16, 19] for
some existing approaches.

This article is structured as follows. In Section 2 the basic problem is
formulated in continuous-time. In Section 3 we consider our first method for
online score estimation and explain the various features associated to it. In
Section 4 our second method for online score estimation is developed. In
Section 5 our numerical results are presented.

###  1.2 Notation

Let $(\mathsf{X},\mathcal{X})$ be a measurable space. We write
$\mathcal{B}_{b}(\mathsf{X})$ for the set of bounded measurable functions,
$\varphi:\mathsf{X}\rightarrow\mathbb{R}^{d}$, $d\in\mathbb{N}$, and
$\mathcal{C}(\mathsf{X})$ for the continuous ones. Let
$\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$;
$\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$ denotes the collection of real-
valued functions that are Lipschitz w.r.t. the Euclidean distance
$\|\cdot\|_{2}$ . That is,
$\varphi\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$ if there exists a
$C<+\infty$ such that for any $(x,y)\in\mathbb{R}^{2d}$,

$|\varphi(x)-\varphi(y)|\leq C\|x-y\|_{2}.$

$\mathcal{N}_{s}(\mu,\Sigma)$ denotes an $s$-dimensional Gaussian law of mean
$\mu$ and covariance $\Sigma$; if $s=1$ we omit subscript $s$. For a
vector/matrix $X$, $X^{*}$ denotes the transpose of $X$. For
$A\in\mathcal{X}$, $\delta_{A}(du)$ denotes the Dirac measure of $A$, and if
$A=\\{x\\}$ with $x\in\mathsf{X}$, we write $\delta_{x}(du)$. For a vector-
valued function in $d$-dimensions $\varphi(x)$ (resp. $d$-dimensional vector
$x$) we write the $i^{\textrm{th}}$-component, $1\leq i\leq d$, as
$\varphi^{(i)}(x)$ (resp. $x^{i}$). For a $d\times q$ matrix $x$ we write the
$(i,j)^{\textrm{th}}$-entry as $x^{(ij)}$. $\mathbb{N}=\\{1,2,\dots\\}$ and
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$.

##  2 Problem Formulation

###  2.1 Preliminaries

We consider the parameter space
$\theta\in\Theta\subset\mathbb{R}^{d_{\theta}}$, with $\Theta$ being compact,
$d_{\theta}\in\mathbb{N}$. The stochastic processes $\\{Y_{t}\\}_{t\geq 0}$,
$\\{X_{t}\\}_{t\geq 0}$ of interest are defined upon the probability triple
$(\Omega,\mathcal{F},\mathbb{P}_{\theta})$, with $Y_{t}\in\mathbb{R}^{d_{y}}$,
$X_{t}\in\mathbb{R}^{d_{x}}$, $d_{y},d_{x}\in\mathbb{N}$, initial conditions
$X_{0}=x_{*}\in\mathbb{R}^{d_{x}}$, $Y_{0}=y_{*}\in\mathbb{R}^{d_{y}}$, and
are determined as the solution of the system of SDEs:

$\displaystyle dY_{t}$ $\displaystyle=h_{\theta}(X_{t})dt+dB_{t};$ (1)
$\displaystyle dX_{t}$
$\displaystyle=b_{\theta}(X_{t})dt+\sigma(X_{t})dW_{t}.$ (2)

Here, for each $\theta\in\Theta$,
$h_{\theta}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{y}}$,
$b_{\theta}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$,
$\sigma:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times d_{x}}$, with
$\sigma$ being of full rank, and $\\{B_{t}\\}_{t\geq 0},\\{W_{t}\\}_{t\geq 0}$
are independent standard Brownian motions of dimension $d_{y}$, $d_{x}$
respectively.

To minimize technical difficulties, the following assumptions are made
throughout the paper:

  * (D1)

    * (i)

$\sigma$ is continuous, bounded; $a(x):=\sigma(x)\sigma(x)^{*}$ is uniformly
elliptic;

    * (ii)

for each $\theta$, $h_{\theta}$ and $b_{\theta}$ are bounded, measurable;
$h_{\theta}^{(i)}\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$, $1\leq
i\leq d_{y}$;

    * (iii)

the gradients
$\nabla_{\theta}h_{\theta}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{y}\times
d_{\theta}}$ and
$\nabla_{\theta}b_{\theta}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times
d_{\theta}}$ exist, and are continuous, bounded, measurable;
$\nabla_{\theta}h_{\theta}^{(ij)}\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$,
$1\leq i\leq d_{y}$, $1\leq j\leq d_{\theta}$;

    * (iv)

let $\phi_{\theta}(x)=(\nabla_{\theta}b_{\theta}(x))^{*}\,a(x)^{-1}\sigma(x)$;
for any $\theta$,
$\phi_{\theta}^{(ij)}\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$,
$1\leq i\leq d_{\theta}$, $1\leq j\leq d_{x}$.

We introduce the probability measure $\overline{\mathbb{P}}_{\theta}$, defined
via the Radon-Nikodym derivative:

$Z_{t,\theta}:=\frac{d\mathbb{P}_{\theta}}{d\overline{\mathbb{P}}_{\theta}}\Big{|}_{\mathcal{F}_{{}_{t}}}=\exp\Big{\\{}\int_{0}^{t}h_{\theta}(X_{s})^{*}dY_{s}-\tfrac{1}{2}\int_{0}^{t}h_{\theta}(X_{s})^{*}h_{\theta}(X_{s})ds\Big{\\}},$
(3)

with $\mathcal{F}_{t}=\sigma(\\{X_{s},Y_{s}\\}_{0\leq s\leq t})$. Henceforth,
$\overline{\mathbb{E}}_{\theta}$ denotes expectation w.r.t.
$\overline{\mathbb{P}}_{\theta}$, so that under
$\overline{\mathbb{P}}_{\theta}$, the process $\\{X_{t}\\}_{t\geq 0}$ follows
the dynamics in (2), whereas $\\{Y_{t}\\}_{t\geq 0}$ is a Brownian motion
independent of $\\{X_{t}\\}_{t\geq 0}$. We define, for
$\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$:

$\gamma_{t,\theta}(\varphi):=\overline{\mathbb{E}}_{\theta}\,\big{[}\,\varphi(X_{t})Z_{t,\theta}\,\big{|}\,\mathcal{Y}_{t}\,\big{]},$

where $\mathcal{Y}_{t}$ is the filtration generated by the process
$\\{Y_{s}\\}_{0\leq s\leq t}$. Our objective is to produce estimates of the
gradient of the score function $\nabla_{\theta}\log(\gamma_{T,\theta}(1))$.

######  Remark 2.1.

To connect the changes of measures with standard likelihood derivations,
notice that – via Girsanov’s theorem – $Z_{t,\theta}$ is the density (w.r.t.
to a Wiener measure) of the distribution of $\\{Y_{s}\\}_{0\leq s\leq t}$
conditionally on $\\{X_{s}\\}_{0\leq s\leq t}$. Then, $\gamma_{T,\theta}(1)$
integrates out $\\{X_{t}\\}_{0\leq t\leq T}$, thus corresponds to the marginal
density – i.e., the likelihood – of the observations $\\{Y_{t}\\}_{0\leq t\leq
T}$.

In our setting, the score function writes as (see e.g. [7]):

$\nabla_{\theta}\log(\gamma_{T,\theta}(1))=\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]},$
(4)

where we have defined:

$\displaystyle\lambda_{T,\theta}:=\int_{0}^{T}(\nabla_{\theta}b_{\theta}(X_{t}))^{*}$
$\displaystyle a(X_{t})^{-1}\sigma(X_{t})dW_{t}$
$\displaystyle+\int_{0}^{T}(\nabla_{\theta}h_{\theta}(X_{t}))^{*}dY_{t}-\int_{0}^{T}(\nabla_{\theta}h_{\theta}(X_{t}))^{*}h_{\theta}(X_{t})dt.$

For completeness, a derivation of (4) can be found in Appendix A. We remark
that one can derive a formula for the score function when $\sigma$ depends
upon $\theta$, which is given in Section 4. We will assume throughout that
$T\in\mathbb{N}$. Note also that an application of Bayes’ rule gives that,
almost surely:

$\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}=\mathbb{E}_{\theta}\,[\,\lambda_{T,\theta}\,|\,\mathcal{Y}_{T}\,],$
(5)

where $\mathbb{E}_{\theta}$ denotes expectation w.r.t. $\mathbb{P}_{\theta}$.

###  2.2 Parameter Estimation

In the offline case suppose one has obtained data $\\{Y_{t}\\}_{0\leq t\leq
T}$. Then it is possible to perform standard gradient descent using (5) and
updating $\theta$ iteratively:

$\displaystyle\theta^{m+1}=\theta^{m}+\alpha_{m}\,\mathbb{E}_{\theta^{m}}\,[\,\lambda_{T,\theta^{m}}\,|\,\mathcal{Y}_{T}\,],$
(6)

where $\alpha_{m}\in\mathbb{R}^{+}$, $m\in\mathbb{N}_{0}$, are decreasing
step-sizes. Instead here we will mainly focus on an online gradient estimation
procedure. To obtain this one can aim to maximize the following limiting
_average log-likelihood_ ,

$\mathcal{L}(\theta)=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}\log\gamma_{s,\theta}(1)ds.$

Let the filter be denoted as
$\pi_{t,\theta}(\varphi)=\mathbb{E}_{\theta}[\varphi(X_{t})|\mathcal{Y}_{s}]$
and using standard arguments (e.g. Lemma 3.29 p. 67 [1]) one can re-write
$\gamma_{t,\theta}(1)$ as

$\log\gamma_{t,\theta}(1)=\int_{0}^{t}\pi_{s,\theta}(h)^{T}dY_{s}-\frac{1}{2}\int_{0}^{t}\pi_{s,\theta}(h)^{T}\pi_{s,\theta}(h)ds.$

Under appropriate stability and regularity conditions for both $Y_{t}$ and
$\pi_{t,\theta}$ (see [30] for more details), then both $\mathcal{L}(\theta)$
and $\nabla_{\theta}\mathcal{L}(\theta)$ are ergodic averages. This means one
can implement stochastic gradient ascent using either
$\nabla_{\theta}\log\gamma_{t,\theta}(1)$ or
$\nabla_{\theta}\log\gamma_{t_{n},\theta}(1)-\nabla_{\theta}\log\gamma_{t_{n-1},\theta}(1)$
for any $t_{n-1}<t_{n}$ as estimates of $\mathcal{L}(\theta)$. Given an
initial $\theta_{0}\in\Theta$, as we obtain the observation path continuously
in time, we will update $\theta$ at times $T\in\mathbb{N}$ using the following
recursion:

$\displaystyle\theta_{T}=\theta_{T-1}+\alpha_{T}\Big{(}\nabla_{\theta}\log(\gamma_{T,\theta_{T-1}}(1))-\nabla_{\theta}\log(\gamma_{T-1,\theta_{T-1}}(1))\Big{)}$
(7)

where, for $T\in\mathbb{N}$, $\alpha_{T}\in\mathbb{R}^{+}$ is a collection of
step-sizes that satisfy $\sum_{T\in\mathbb{N}}\alpha_{T}=\infty$,
$\sum_{T\in\mathbb{N}}\alpha_{T}^{2}<\infty$ to ensure convergence of the
estimation as $T\rightarrow\infty$; see [3, 20] for details. This scheme can
provide an online estimate for the parameter vector as data arrive. Steps are
performed at $\mathcal{O}(1)$ times to ensure that enough information has
accumulated to update the parameter. The adoption of unit times is made only
for notational convenience. As both recursions (6) and (7) cannot be computed
exactly, we focus upon methodologies that approximate the score function
$\nabla_{\theta}\log(\gamma_{T,\theta}(1))$.

##  3 Direct Feynman-Kac Formulation

###  3.1 Discretized Model

In practice, we will have to work with a discretization of the model in
(1)-(2). We assume access to path of the data $\\{Y_{t}\\}_{0\leq t\leq T}$
which is available up-to an (almost) arbitrarily fine level of time
discretization. This would be a very finely discretized path, as accessing the
actual continuous path of observation is not possible; this point is discussed
later on. One could focus on a time-discretization of either side of (5),
however, as is conventional in the literature (e.g. [1, 19]) we focus on the
left hand side.

Let $l\in\mathbb{N}_{0}$ and consider an Euler-Maruyama time-discretization
with step-size $\Delta_{l}=2^{-l}$. That is, for
$k\in\\{1,2,\dots,T/\Delta_{l}\\}$:

$\displaystyle\widetilde{X}_{k\Delta_{l}}=\widetilde{X}_{(k-1)\Delta_{l}}+b_{\theta}(\widetilde{X}_{(k-1)\Delta_{l}})\Delta_{l}+\sigma(\widetilde{X}_{(k-1)\Delta_{l}})(W_{k\Delta_{l}}-W_{(k-1)\Delta_{l}}),\quad\widetilde{X}_{0}=x_{*}.$
(8)

Note that the Brownian motion in (8) is the same as in (2) under both
$\mathbb{P}_{\theta}$ and $\overline{\mathbb{P}}_{\theta}$. We set:

$\displaystyle\lambda_{T,\theta}^{l}(x_{0},$ $\displaystyle
x_{\Delta_{l}},\dots,x_{T}):=\sum_{k=0}^{T/\Delta_{l}-1}\Big{\\{}\,(\nabla_{\theta}b_{\theta}(x_{k\Delta_{l}}))^{*}a(x_{k\Delta_{l}})^{-1}\sigma(x_{k\Delta_{l}})(W_{(k+1)\Delta_{l}}-W_{k\Delta_{l}})$
$\displaystyle+(\nabla_{\theta}h_{\theta}(x_{k\Delta_{l}}))^{*}(Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}})-(\nabla_{\theta}h_{\theta}(x_{k\Delta_{l}}))^{*}h_{\theta}(x_{k\Delta_{l}})\Delta_{l}\,\Big{\\}}.$
(9)

We remark that $\lambda_{T,\theta}^{l}$ is a function also of the
observations, but this dependence is suppressed from the notation. For
$k\in\\{0,1,\dots,T/\Delta_{l}-1\\}$, we define:

$\displaystyle
g_{k,\theta}^{l}(x_{k\Delta_{l}}):=\exp\Big{\\{}h_{\theta}(x_{k\Delta_{l}})^{*}(y_{(k+1)\Delta_{l}}-y_{k\Delta_{l}})-\tfrac{\Delta_{l}}{2}h_{\theta}(x_{k\Delta_{l}})^{*}h_{\theta}(x_{k\Delta_{l}})\Big{\\}}.$

Note that:

$\displaystyle Z_{T,\theta}^{l}(x_{0},x_{\Delta_{l}},\dots,x_{T}):=$
$\displaystyle\prod_{k=0}^{T/\Delta_{l}-1}g_{k,\theta}^{l}(x_{k\Delta_{l}})$
$\displaystyle=$
$\displaystyle\exp\Big{\\{}\sum_{k=0}^{T/\Delta_{l}-1}\big{[}\,h_{\theta}(x_{k\Delta_{l}})^{*}(y_{(k+1)\Delta_{l}}-y_{k\Delta_{l}})-\tfrac{\Delta_{l}}{2}h_{\theta}(x_{k\Delta_{l}})^{*}h_{\theta}(x_{k\Delta_{l}})\,\big{]}\Big{\\}}$

is a time-discretization of $Z_{T,\theta}$. We thus obtain the discretized
approximation of the score function
$\nabla_{\theta}\log(\gamma_{T,\theta}(1))$:

$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1)):=\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}.$
(10)

We have the following result which establishes the convergence of our Euler
approximation. Below $\|\cdot\|_{2}$ is the $L_{2}-$norm for vectors. The
proof is given in Appendix B.

######  Theorem 3.1.

Assume (D(D1)). Then for any $(r,T)\in[1,\infty)\times\mathbb{N}$ there exists
a $C<+\infty$ such that for any $l\in\mathbb{N}_{0}$

$\mathbb{E}_{\theta}\left[\left\|\nabla_{\theta}\log(\gamma_{T,\theta}(1))-\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))\right\|_{2}^{r}\right]^{1/r}\leq
C\Delta_{l}^{1/2}.$

The result is fairly standard, but we note it is not a simple application of
results on SDEs in filtering (e.g. [25, 31]) and this is reflected in the
proof. The rate of convergence of the approximation will be very relevant for
some of our subsequent discussions.

###  3.2 Backward Feynman-Kac Model and Particle Smoothing

From herein the $\widetilde{X}$ notation is dropped for simplicity. Consider
the time interval $[k,k+1]$ and the $k$-th update of (7). We define the
discrete-time approximation (at level $l$) as:

$\displaystyle u_{k,l}=(x_{k+\Delta_{l}},\dots,x_{k+1})\in
E_{l}:=(\mathbb{R}^{d_{x}})^{\Delta_{l}^{-1}},\quad k\in\mathbb{N}_{0}.$

Recall (3). A discrete-time approximation of $p_{\theta}(\\{Y_{t}\\}_{k\leq
t\leq k+1}|\\{X_{t}\\}_{0\leq t\leq T})$ is:

$\displaystyle G_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle:=\prod_{p=0}^{\Delta_{l}^{-1}-1}g_{k+p,\theta}^{l}(x_{k+p\Delta_{l}}),$

where we set $u_{-1,l}=x_{*}$, for each $l\in\mathbb{N}_{0}$. We denote by
$m_{\theta}^{l}$ the Euler transition density induced by time-discretisation
(8), and then write the initial distribution and Markov transition kernel for
the discrete-time process with $k\in\mathbb{N}$ as follows:

$\displaystyle\eta_{0,\theta}^{l}(du_{0,l})$
$\displaystyle=\prod_{p=1}^{\Delta_{l}^{-1}}m_{\theta}^{l}(x_{(p-1)\Delta_{l}},x_{p\Delta_{l}})dx_{p\Delta_{l}};$
$\displaystyle\vspace{-5pt}M_{\theta}^{l}(u_{k-1,l},du_{k,l})$
$\displaystyle=\prod_{p=1}^{\Delta_{l}^{-1}}m_{\theta}^{l}(x_{k+(p-1)\Delta_{l}},x_{k+p\Delta_{l}})dx_{k+p\Delta_{l}}.$

######  Remark 3.1.

The definition of $G_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$,
$M_{\theta}^{l}(u_{k-1,l},du_{k,l})$ implies: i)
$G_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$ involves $u_{k-1,l}$ only via its very
last element, $x_{k}$; ii) the dynamics of $u_{k,l}$ conditionally on
$u_{k-1,l}$ depend only on the very last element, $x_{k}$, of $u_{k-1,l}$.

We can now state the discrete-time filtering distribution for
$k\in\mathbb{N}_{0}$:

$\pi_{k,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{k,l})\big{)}:=\frac{\big{(}\prod_{p=0}^{k}G_{p,\theta}^{l}(u_{p-1,l},u_{p,l})\big{)}\,\eta_{0,\theta}^{l}(du_{0,l})\prod_{p=1}^{k}M_{\theta}^{l}(u_{p-1,l},du_{p,l})}{\int_{E_{l}^{k+1}}\big{(}\prod_{p=0}^{k}G_{p,\theta}^{l}(u_{p-1,l},u_{p,l})\big{)}\,\eta_{0,\theta}^{l}(du_{0,l})\prod_{p=1}^{k}M_{\theta}^{l}(u_{p-1,l},du_{p,l})}.$
(11)

That is, $\pi_{k,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{k,l})\big{)}$ is a
discrete-time approximation of the filtering distribution:

$\pi_{k,\theta}\big{(}d(\\{X_{t}\\}_{0\leq t\leq
k})\big{)}:=\mathbb{P}_{\theta}(d\\{X_{t}\\}_{0\leq t\leq
k}|\\{Y_{t}\\}_{0\leq t\leq k}).$

Expression (11) corresponds to a standard Feynman-Kac model (see e.g. [11]),
thus one can approximate the involved filtering distributions via the
corresponding Monte Carlo methodology.

We develop a Monte Carlo method for the approximation of the discretised score
function in (10). This is accomplished by presenting a backward formula for
(10). We define for any $p\in\mathbb{N}_{0}$:

$\displaystyle f_{\theta}^{l}(x_{p\Delta_{l}},x_{(p+1)\Delta_{l}}):=$
$\displaystyle(\nabla_{\theta}b_{\theta}(x_{p\Delta_{l}}))^{*}a(x_{p\Delta_{l}})^{-1}\sigma(x_{p\Delta_{l}})(W_{(p+1)\Delta_{l}}-W_{p\Delta_{l}})$
$\displaystyle\quad+(\nabla_{\theta}h_{\theta}(x_{p\Delta_{l}}))^{*}(Y_{(p+1)\Delta_{l}}-Y_{p\Delta_{l}})-(\nabla_{\theta}h_{\theta}(x_{p\Delta_{l}}))^{*}h_{\theta}(x_{p\Delta_{l}})\Delta_{l}.$

and let:

$\displaystyle\Lambda_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle:=\sum_{p=0}^{\Delta_{l}^{-1}-1}f_{\theta}^{l}(x_{k+p\Delta_{l}},x_{k+(p+1)\Delta_{l}});$
$\displaystyle F_{T,\theta}^{l}(u_{0,l},\dots,u_{T-1,l})$
$\displaystyle:=\sum_{k=0}^{T-1}\Lambda_{k,\theta}^{l}(u_{k-1,l},u_{k,l})\,\,\Big{(}\equiv\lambda_{T,\theta}^{l}(x_{0},x_{\Delta_{l}},\dots,x_{T})\Big{)},$
(12)

for $\lambda_{T,\theta}^{l}(x_{0},x_{\Delta_{l}},\dots,x_{T})$ as defined in
(9) and used in the score function approximation (10). Thus:

$\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))=\int_{E_{l}^{T}}F_{T,\theta}^{l}(u_{0,l},\dots,u_{T-1,l})\,\mathbb{Q}_{T-1,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{T-1,l})\big{)},$
(13)

where $\mathbb{Q}_{T-1,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{T-1,l})\big{)}$ is
a time-discretisation of the smoothing law:

$\mathbb{Q}_{T-1,\theta}(d\\{X_{t}\\}_{0\leq t\leq
T}):=\mathbb{P}_{\theta}(d\\{X_{t}\\}_{0\leq t\leq T}|\\{Y_{t}\\}_{0\leq t\leq
T}).$

Now, by the time-reversal formula for hidden Markov models (see e.g. [12, 13])
one has:

$\displaystyle\mathbb{Q}_{T-1,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{T-1,l})\big{)}:=\pi_{T-1,\theta}^{l}(du_{T-1,l})\prod_{k=1}^{T-1}B_{k,\theta,\pi_{k-1,\theta}^{l}}^{l}(u_{k,l},du_{k-1,l}),$

for the backward Markov kernel:

$\displaystyle B_{k,\theta,\pi_{k-1,\theta}^{l}}^{l}(u_{k,l},$ $\displaystyle
du_{k-1,l}):=\frac{\pi_{k-1,\theta}^{l}(du_{k-1,l})\,G_{k,\theta}^{l}(u_{k-1,l},u_{k,l})m_{\theta}^{l}(u_{k-1,l},u_{k,l})}{\pi_{k-1,\theta}^{l}(G_{k,\theta}^{l}(\cdot,u_{k,l})m_{\theta}^{l}(\cdot,u_{k,l}))},$
(14)

under the standard notation:

$\pi_{k-1,\theta}^{l}(G_{k,\theta}^{l}(\cdot,u_{k,l})m_{\theta}^{l}(\cdot,u_{k,l}))=\int_{E_{l}}\pi_{k-1,\theta}^{l}(du_{k-1,l})G_{k,\theta}^{l}(u_{k-1,l},u_{k,l})m_{\theta}^{l}(u_{k-1,l},u_{k,l}).$

######  Remark 3.2.

The model structure gives important cancellations in (14), so that:

$B_{k,\theta,\pi_{k-1,\theta}^{l}}^{l}(u_{k,l},du_{k-1,l})\equiv\frac{\pi_{k-1,\theta}^{l}(du_{k-1,l})g_{k,\theta}^{l}(x_{k})m_{\theta}^{l}(x_{k},x_{k+\Delta_{l}})}{\int_{E_{l}}\pi_{k-1,\theta}^{l}(du_{k-1,l})g_{k,\theta}^{l}(x_{k})m_{\theta}^{l}(x_{k},x_{k+\Delta_{l}})}.$

The objective now is to approximate the right hand side of (13) using particle
approximations. Our online particle approximation of the gradient of the log-
likelihood in (13), for a given $l\in\mathbb{N}_{0}$ is presented in Algorithm
1. Our estimates are given in (15) and (17) in Algorithm 1. The approach is
the method introduced in [12, 13].

Algorithm 1 Online Score Function Estimation for a given $l\in\mathbb{N}_{0}$.

  1. 1.

For $i\in\\{1,\dots,N\\}$, sample $u_{0,l}^{i}$ i.i.d. from
$\eta_{0,\theta}^{l}(\cdot)$. The estimate of
$\nabla_{\theta}\log(\gamma_{1,\theta}^{l}(1))$ is:

$\widehat{\nabla_{\theta}\log(\gamma_{1,\theta}^{l}(1))}:=\frac{\sum_{i=1}^{N}G_{0,\theta}^{l}(x_{*},u_{0,l}^{i})\Lambda_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}{\sum_{i=1}^{N}G_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}.$
(15)

Set $k=1$, and for $i\in\\{1,\dots,N\\}$, $\check{u}_{-1,l}^{i}=x_{*}$.

  2. 2.

For $i\in\\{1,\dots,N\\}$ sample $\check{u}_{k-1,l}^{i}$ from:

$\sum_{i=1}^{N}\frac{G_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{i},u_{k-1,l}^{i})}{\sum_{j=1}^{N}G_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{j},u_{k-1,l}^{j})}\delta_{\\{u_{k-1,l}^{i}\\}}(\cdot).$

If $k=1$, for $i\in\\{1,\dots,N\\}$, set
$F_{k-1,\theta}^{l,N}(\check{u}_{0,l}^{i})=\Lambda_{0,\theta}^{l}(x_{*},\check{u}_{0,l}^{i})$.

  3. 3.

For $i\in\\{1,\dots,N\\}$, sample $u_{k,l}^{i}$ from
$M_{\theta}^{l}(\check{u}_{k-1,l}^{i},\cdot)$. For $i\in\\{1,\dots,N\\}$,
compute:

$F_{k,\theta}^{l,N}(u_{k,l}^{i})=\frac{\sum_{j=1}^{N}g_{k,\theta}^{l}(\check{x}_{k}^{j})m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})\\{F_{k-1,\theta}^{l,N}(\check{u}_{k-1,l}^{j})+\Lambda_{k,\theta}^{l}(\check{u}_{k-1,l}^{j},u_{k,l}^{i})\\}}{\sum_{j=1}^{N}g_{k,\theta}^{l}(\check{x}_{k}^{j})m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})}.$
(16)

The estimate of $\nabla_{\theta}\log(\gamma_{k+1,\theta}^{l}(1))$ is:

$\widehat{\nabla_{\theta}\log(\gamma_{k+1,\theta}^{l}(1))}:=\frac{\sum_{i=1}^{N}G_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})F_{k,\theta}^{l,N}(u_{k,l}^{i})}{\sum_{i=1}^{N}G_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})}.$
(17)

Set $k=k+1$ and return to the start of 2..

###  3.3 Discussion of Algorithm 1

There are several remarks worth making, before proceeding. Firstly, the cost
of this algorithm per unit time is $\mathcal{O}(N\Delta_{l}^{-1}+N^{2})$. In
detail, the cost of the particle filter is $\mathcal{O}(N\Delta_{l}^{-1})$.
The cost of calculating $F_{k,\theta}^{l,N}(u_{k,l}^{i})$, $1\leq i\leq N$, in
(16), is $\mathcal{O}(N^{2})$; $\Delta_{l}$ is not involved here due to
cancellations – see 3.2. The cost of (17) is $\mathcal{O}(N)$, given the
particle filter has already been executed. There are several implications of
this remark. Based upon the results in [12] and Theorem 3.1, in a sequel work
we prove, under appropriate assumptions, we will have the following MSE, for
$(k,N)\in\mathbb{N}^{2}$:

$\mathbb{E}_{\theta}\Big{[}\,\Big{\|}\,\widehat{\nabla_{\theta}\log(\gamma_{k,\theta}^{l}(1))}-\nabla_{\theta}\log(\gamma_{k,\theta}(1))\,\Big{\|}_{2}^{2}\,\Big{]}\leq
C\Big{(}\frac{1}{N}+\Delta_{l}\Big{)},$ (18)

where $C$ is a constant that does not depend on $N$ or $l$. To achieve an MSE
of $\mathcal{O}(\epsilon^{2})$ for some $\epsilon>0$ given, one sets $l$ so
that $\Delta_{l}=\epsilon^{2}$ (i.e. $l=\mathcal{O}(|\log(\epsilon)|)$) and
$N=\mathcal{O}(\epsilon^{-2})$. The cost per unit time of doing this is then
$\mathcal{O}(\epsilon^{-4})$. We note that to choose $l$ as specified, one has
to have access to an appropriately finely observed data set and this is
assumed throughout. Typically, one could use a multilevel Monte Carlo method,
as in [19], to reduce the cost to achieve an MSE of
$\mathcal{O}(\epsilon^{2})$. However, in this case as the $\mathcal{O}(N^{2})$
cost dominates and does not depend on $l$, one can easily check that such a
variance reduction method _will not improve_ our particle method. To
understand this, one can prove a bound on the MSE, for instance of the type
conjectured later in this article (37), and then try to minimize the cost, by
selecting the appropriate number of samples on each level to obtain a given
MSE. This latter problem leads to a constrained minimization problem that can
be solved using Lagrange multipliers (as in e.g. [8]), but one can show that
this yields that the order of the cost to achieve an MSE of
$\mathcal{O}(\epsilon^{2})$ is still $\mathcal{O}(\epsilon^{-4})$.

Secondly, it is important to note that the method in [22] can reduce the cost
of online score estimation to $\mathcal{O}(N\Delta_{l}^{-1})$ per unit time.
However, in order to do so, one requires that $m_{\theta}^{l}(x,x^{\prime})$
is uniformly lower-bounded in $(x,x^{\prime})$, which does not typically occur
for Euler-discretized diffusion densities. As a result, we only use the
approach shown in Algorithm 1. Note that [15], in a different but related
context, considers using unbiased and non-negative estimates of the transition
density in the approach in [22], but such estimates are not always available.

Thirdly and rather importantly, there is a potential issue related to the
construction of the algorithm. We have started with a continuous-time formula,
discretized it and applied what are essentially discrete-time methods for
smoothing of additive functionals. A serious caveat is that the algorithm is
not well-defined as $l\rightarrow\infty$, which is what we mean by saying it
has _no (Wiener) path-space formulation_. The source of the issue is related
to using approximations of the transition density
$m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})$ in (16), which can
degenerate when $l$ is high. This will result in increasing Monte Carlo
variance and computational cost and may mean that $C$ in (18) explodes
exponentially in $l$. We refer the reader to [35, Figure 1.1] for a numerical
example. This issue has manifested itself in MCMC schemes for inferring fully
observed SDEs (see e.g. [23]), but in the context of particle smoothing and
Algorithm 1 there are additional considerations. The resampling operation
introduces _discontinuities_. Often such terminology refers to the lack of
continuity of the transition density:

$\theta\rightarrow\widehat{\nabla_{\theta}\log(\gamma_{k,\theta}^{l}(1))},$

but here we are interested in the behavior of
$m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})$ when we combine
points $\check{x}_{k}^{j}$ and $x_{k+\Delta_{l}}^{i}$ that are intrinsically
not obtained in a continuous way as $\Delta_{l}$ diminishes. This issue has
not received attention in earlier numerical studies, but remains a concern. As
a result, we now consider defining an algorithm that is robust to the size of
the time-discretization mesh and hence has a path-space formulation.

##  4 Path-Space Feynman-Kac Formulation

###  4.1 Data Augmentation using Bridges

We begin this section with a review of the method in [28, 33]. For simplicity
we consider the case $t\in[0,1]$ and let
$\mathbf{X}:=\\{X_{t}\\}_{t\in[0,1]}$, and
$\mathbf{W}:=\\{W_{t}\\}_{t\in[0,1]}$. Let also $p_{\theta}(x,t;x^{\prime},1)$
denote the unknown transition density from time $t$ to $1$ associated to (2)
and let also $p_{\theta}(x,x^{\prime}):=p_{\theta}(x,0;x^{\prime},1)$. Suppose
one could sample from $p_{\theta}$ to obtain
$(x,x^{\prime})\in\mathbb{R}^{2d_{x}}$. Then one can interpolate these points
by using a bridge process which has a drift given by
$b_{\theta}(x)+a(x)\nabla_{x}\log{p}_{\theta}(x,t;x^{\prime},1)$, as we will
explain below. Let $\overline{\mathbb{P}}_{\theta,x,x^{\prime}}$ denote the
law of the solution of the SDE (2), on $[0,1]$, started at $x_{*}=x$ and
conditioned to hit $x^{\prime}$ at time $1$.

As $p_{\theta}$ is intractable in general, we consider a user-specified
auxiliary process $\\{\tilde{X}_{t}\\}_{t\in[0,1]}$ following:

$\displaystyle
d\tilde{X}_{t}=\tilde{b}_{\theta}(t,\tilde{X}_{t})dt+\tilde{\sigma}(t,\tilde{X}_{t})dW_{t},\quad
t\in[0,1],\quad\leavevmode\nobreak\ \tilde{X}_{0}=x,$ (19)

where for each parameter value $\theta\in\Theta$,
$\tilde{b}_{\theta}:[0,1]\times\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$
and $\tilde{\sigma}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times
d_{x}}$ is such that
$\tilde{a}(1,x^{\prime}):=\tilde{\sigma}(1,x^{\prime})\tilde{\sigma}(1,x^{\prime})^{*}\equiv
a(x^{\prime})$. Most importantly, (19) is chosen so that its transition
density is available. To avoid confusion – as the specification of process
(19) can involve parameter $\theta$ and a given ending position $x^{\prime}$ –
we note that the transition density of (19) from time $t$ to time $1$
corresponds to a mapping $z\to\tilde{p}_{\theta,x^{\prime}}(x,t;z,1)$. We also
use the notation
$\tilde{p}_{\theta,x^{\prime}}(x,z):=\tilde{p}_{\theta,x^{\prime}}(x,0;z,1)$.
One possible choice is to use an Ornstein-Uhlenbeck process (e.g. obtained
using linearizations or variational inference with (2) [28, Section 1.3]); see
also [28, Section 2.2] for technical conditions on
$\tilde{b}_{\theta},\tilde{a},\tilde{p}_{\theta,x^{\prime}}$. The main purpose
of $\\{\tilde{X}_{t}\\}_{t\in[0,1]}$ is to construct another process
$\\{X_{t}^{\circ}\\}_{t\in[0,1]}$ conditioned to hit a given $x^{\prime}$ at
$t=1$. The latter will form an importance proposal for
$\\{X_{t}\\}_{t\in[0,1]}$. Let:

$\displaystyle
dX^{\circ}_{t}=b_{\theta}^{\circ}(t,X^{\circ}_{t};x^{\prime})dt+\sigma(X^{\circ}_{t})dW_{t},\quad
t\in[0,1],\quad\leavevmode\nobreak\ X^{\circ}_{0}=x,$ (20)

where:

$b_{\theta}^{\circ}(t,x;x^{\prime})=b_{\theta}(x)+a(x)\nabla_{x}\log\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1),$

and denote by $\mathbb{P}^{\circ}_{\theta,x,x^{\prime}}$ the probability law
of the solution of (20). The SDE in (20) gives rise to a function:

$\mathbf{W}\rightarrow C_{\theta}(x,\mathbf{W},x^{\prime}),$ (21)

mapping the driving Wiener noise $\mathbf{W}$ to the solution of (20), so we
have effectively reparameterized the problem from $\mathbf{X}$ to
$(\mathbf{W},x^{\prime})$.

Now, following [28], the two measures
$\overline{\mathbb{P}}_{\theta,x,x^{\prime}}$ and
$\mathbb{P}^{\circ}_{\theta,x,x^{\prime}}$ are absolutely continuous w.r.t.
each other, with Radon-Nikodym derivative:

$\frac{d\overline{\mathbb{P}}_{\theta,x,x^{\prime}}}{d\mathbb{P}^{\circ}_{\theta,x,x^{\prime}}}(\mathbf{X})=\exp\Big{\\{}\int_{0}^{1}L_{\theta}(t,X_{t})dt\Big{\\}}\times\frac{\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})}{{p}_{\theta}(x,x^{\prime})},$
(22)

where:

$\displaystyle L_{\theta}$
$\displaystyle(t,x):=\big{(}b_{\theta}(x)-\tilde{b}_{\theta}(t,x)\big{)}^{*}\,\nabla_{x}\log\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1)$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!-\tfrac{1}{2}\textrm{Tr}\,\Big{\\{}\,\big{[}a(x)-\tilde{a}(t,x)\big{]}\big{[}-\nabla_{x}^{2}\log\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1)-\nabla_{x}\log\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1)\nabla_{x}\log\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1)^{*}\big{]}\,\Big{\\}}$

with $\textrm{Tr}(\cdot)$ denoting the trace of a squared matrix. Note that,
in the case when $\sigma=\sigma(x)$ is not a constant function, then,
typically, $x^{\prime}\rightarrow\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})$
will not integrate to $1$ and will give rise to a non-trivial distribution to
sample from. As the complete algorithm will require being able to sample from
the transition density, we rewrite:

$\frac{d\overline{\mathbb{P}}_{\theta,x,x^{\prime}}}{d\widetilde{\mathbb{P}}_{\theta,x,x^{\prime}}}(\mathbf{X})=\exp\Big{\\{}\int_{0}^{1}L_{\theta}(t,X_{t})dt\Big{\\}}\times\frac{\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})}{{p}_{\theta}(x,x^{\prime})\hat{{p}}_{\theta}(x,x^{\prime})}\times\hat{{p}}_{\theta}(x,x^{\prime}),$
(23)

where an arbitrary, tractable and easy to sample density
$\hat{{p}}_{\theta}(x,x^{\prime})$ is used to sample $x^{\prime}$.

###  4.2 Estimation of Score in Continuous-Time

We return to the expression of the score function in (4) and use the
alternative change of measure described above. Consider the processes:

$\displaystyle\mathbf{X}_{k}:=\\{X_{t}\\}_{t\in[k,k+1]},\quad\mathbf{Y}_{k}:=\\{Y_{t}\\}_{t\in[k,k+1]},\quad
k\in\mathbb{N}_{0}.$

We introduce the following notation:

$\displaystyle\Psi_{\theta}(\mathbf{X}_{k})$
$\displaystyle=\int_{k}^{k+1}L_{\theta}(t,X_{t})dt;$ $\displaystyle
J_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})$
$\displaystyle=\int_{k}^{k+1}h_{\theta}(X_{t})^{*}dY_{t}-\tfrac{1}{2}\int_{k}^{k+1}h_{\theta}(X_{t})^{*}h_{\theta}(X_{t})dt;$
$\displaystyle\Lambda_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})$
$\displaystyle=\int_{k}^{k+1}(\nabla_{\theta}b_{\theta}(X_{t}))^{*}a(X_{t})^{-1}(dX_{t}-b_{\theta}(X_{t})dt)$
$\displaystyle\qquad\qquad+\int_{k}^{k+1}(\nabla_{\theta}h_{\theta}(X_{t}))^{*}dY_{t}-\int_{k}^{k+1}(\nabla_{\theta}h_{\theta}(X_{t}))^{*}h_{\theta}(X_{t})dt;$
$\displaystyle\Phi_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})$
$\displaystyle=J_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})+\Psi_{\theta}(\mathbf{X}_{k})+\log\frac{\tilde{p}_{\theta,x_{k+1}}(x_{k},x_{k+1})}{\hat{{p}}_{\theta}(x_{k},x_{k+1})}.$

Note all the integrands can be computed point wise.

######  Remark 4.1.

The Wiener process in (21) is defined on the time interval $[0,1]$, thus so is
the transform $C_{\theta}(x,\mathbf{W},x^{\prime})$. In the derivations below,
one needs to calculate
$J_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})$ and
$\Lambda_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})$,
for $\mathbf{W}_{k}$’s that correspond to samples from the Wiener measure on
$[0,1]$. With some abuse of notation, it is to be understood that the path
$C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1})$ is ‘shifted’ from $[0,1]$ to
$[k,k+1]$, so all quantities below agree with the notation introduced above.
Also, the calculation of
$\Psi_{\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}))$ will be required,
but this should create no confusion.

Under these definitions, for any $T\in\mathbb{N}$ the score function in (4)
can be rewritten as:

$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}(1))\equiv\frac{\overline{\mathbb{E}}_{\theta}\Big{[}\,\big{(}\sum_{k=0}^{T-1}\Lambda_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})\,\big{)}\exp\big{\\{}\sum_{k=0}^{T-1}J_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}}{\overline{\mathbb{E}}_{\theta}\,\Big{[}\,\exp\big{\\{}\sum_{k=0}^{T-1}J_{k,\theta}(\mathbf{X}_{k},\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}}.$

Making use of the transform (21) and the density expression in (22), we can
equivalently write:

$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}(1))=$ (24)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\frac{\widetilde{\mathbb{E}}_{\theta}\,\Big{[}\,\big{(}\sum_{k=0}^{T-1}\Lambda_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\big{)}\exp\big{\\{}\sum_{k=0}^{T-1}\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}}{\widetilde{\mathbb{E}}_{\theta}\,\Big{[}\,\exp\big{\\{}\sum_{k=0}^{T-1}\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}},$

where, the expectation
$\widetilde{\mathbb{E}}_{\theta}\,[\,\cdot\,|\,\mathcal{Y}_{T}\,]$ is
considered under the probability measure:

$\displaystyle\widetilde{\mathbb{P}}_{\theta}\,\big{(}\,d(\mathbf{W}_{0},x_{1},\ldots,\mathbf{W}_{T-1},x_{T})\,\big{)}:=\bigotimes_{k=0}^{T-1}\big{[}\,\mathbb{W}(d\mathbf{W}_{k})\otimes\hat{p}_{\theta}(x_{k},x_{k+1})dx_{k+1}\,\big{]},$
(25)

independently of $\mathcal{Y}_{T}$; here, $\mathbb{W}$ is the standard Wiener
measure on $[0,1]$ and $x_{0}=x_{*}$.

######  Remark 4.2.

The approach that has been adopted here can also be used if $\sigma$ depends
upon $\theta$. Assuming the formula is well-defined, one would have a score
function with an expression of the type:

$\frac{\widetilde{\mathbb{E}}_{\theta}\,\Big{[}\,\big{(}\sum_{k=0}^{T-1}\Xi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1})\big{)}\exp\big{\\{}\sum_{k=0}^{T-1}\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}}{\widetilde{\mathbb{E}}_{\theta}\,\Big{[}\,\exp\big{\\{}\sum_{k=0}^{T-1}\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\big{\\}}\,\big{|}\,\mathcal{Y}_{T}\,\Big{]}},$

where:

$\Xi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1})=\\{\nabla_{\theta}\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k})\\}+\nabla_{\theta}\log\\{\hat{p}_{\theta}(x_{k},x_{k+1})\\},$

and with an appropriate modification of the approach to allow $\sigma$ to
depend on $\theta$. To keep consistency with the ideas in Section 3 we do not
consider the formula from herein, but remark that extension of the forthcoming
methodology to this case is straightforward.

######  Remark 4.3.

Recent advances in [5] have extended the construct of the auxiliary bridge
process – developed via (19), (20) herein – to the setting of _hypoelliptic_
SDEs. Though we do not pursue this direction here for the purpose of easing
the exposition, we remark that, given these new developments, one can now, in
principle, obtain score function estimates – thus, also carry out parameter
inference – in the hypoelliptic regime along the same steps we follow in the
current work.

The expression in (24), together with the (trivially) Markovian dynamics for
the process in (25) allow one to construct a backward Feynman-Kac type formula
as in (13). To better connect the approach here and that in Section 3, define
$u_{k}=(\mathbf{W}_{k},x_{k+1})$ for $k\in\mathbb{N}_{0}$, $u_{-1}=x_{*}$ and:

$\displaystyle\Lambda_{k,\theta}^{C}(u_{k-1},u_{k})$
$\displaystyle=\Lambda_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k});$
$\displaystyle\Phi_{k,\theta}^{C}(u_{k-1},u_{k})$
$\displaystyle=\Phi_{k,\theta}(C_{\theta}(x_{k},\mathbf{W}_{k},x_{k+1}),\mathbf{Y}_{k});$
$\displaystyle F_{T,\theta}(u_{0},\dots,u_{T-1})$
$\displaystyle=\sum_{k=0}^{T-1}\Lambda_{k,\theta}^{C}(u_{k-1},u_{k}).$

Superscript $C$ is motivated by the well-posedness of the formula in
_continuous-time_ path-space. Set:

$\pi_{k,\theta}\big{(}\,d(u_{0},\dots,u_{k})\,\big{)}:=\frac{\Big{(}\prod_{p=0}^{k}\exp\\{\Phi_{p,\theta}^{C}(u_{p-1},u_{p})\\}\Big{)}\widetilde{\mathbb{P}}_{\theta}\,\big{(}\,d(u_{0},\ldots,u_{k})\,\big{)}}{\int_{E^{k+1}}\Big{(}\prod_{p=0}^{k}\exp\\{\Phi_{p,\theta}^{C}(u_{p-1},u_{p})\\}\Big{)}\widetilde{\mathbb{P}}_{\theta}\,\big{(}\,d(u_{0},\ldots,u_{k})\,\big{)}},$

with $E=C([0,1],\mathbb{R}^{d_{x}})$. Then, we have the representation:

$\nabla_{\theta}\log(\gamma_{T,\theta}(1))=\int_{E^{T}}F_{T,\theta}(u_{0},\dots,u_{T-1})\,\mathbb{Q}_{T-1,\theta}\big{(}\,d(u_{0},\dots,u_{T-1})\,\big{)},$
(26)

where:

$\mathbb{Q}_{T-1,\theta}\big{(}\,d(u_{0},\dots,u_{T-1})\,\big{)}=\pi_{T-1,\theta}(du_{T-1})\prod_{k=1}^{T-1}B_{k,\theta,\pi_{k-1,\theta}}(u_{k},du_{k-1}),$

and:

$B_{k,\theta,\pi_{k-1,\theta}}(u_{k},du_{k-1}):=\frac{\pi_{k-1,\theta}(du_{k-1})\exp\\{\Phi_{k,\theta}^{C}(u_{k-1},u_{k})\\}\hat{p}_{\theta}(x_{k},x_{k+1})}{\pi_{k-1,\theta}(\exp\\{\Phi_{k,\theta}^{C}(\cdot,u_{k})\\}\hat{p}_{\theta}(\cdot,x_{k+1}))}.$

We remark that formula (26) is a type of backward Feynman-Kac formula in
continuous-time, which to our knowledge is new.

As in the case of Algorithm 1, an effective Monte Carlo approximation of such
a smoothing expectation (26) is given in Algorithm 2. The estimates of the
score function are given in equations (27)-(28) in Algorithm 2.

Algorithm 2 Online Score Function Estimation on Path-Space

  1. 1.

For $i\in\\{1,\dots,N\\}$, sample $u_{0}^{i}$ i.i.d. from
$\mathbb{W}(\cdot)\otimes\hat{p}_{\theta}(x_{*},\cdot)$. The estimate of
$\nabla_{\theta}\log(\gamma_{1,\theta}(1))$ is:

$\widehat{\nabla_{\theta}\log(\gamma_{1,\theta}(1))}:=\frac{\sum_{i=1}^{N}\exp\\{\Phi_{0,\theta}^{C}(x_{*},u_{0}^{i})\\}\Lambda_{0,\theta}^{C}(x_{*},u_{0}^{i})}{\sum_{i=1}^{N}\exp\\{\Phi_{0,\theta}^{C}(x_{*},u_{0}^{i})\\}}.$
(27)

Set $k=1$ and for $i\in\\{1,\dots,N\\}$, $\check{u}_{-1}^{i}=x_{*}$.

  2. 2.

For $i\in\\{1,\dots,N\\}$ sample $\check{u}_{k-1}^{i}$ from:

$\sum_{i=1}^{N}\frac{\exp\\{\Phi_{k-1,\theta}^{C}(\check{u}_{k-2}^{i},u_{k-1}^{i})\\}}{\sum_{j=1}^{N}\exp\\{\Phi_{k-1,\theta}^{C}(\check{u}_{k-2}^{j},u_{k-1}^{j})\\}}\delta_{\\{u_{k-1}^{i}\\}}(\cdot).$

If $k=1$, for $i\in\\{1,\dots,N\\}$, set
$\widetilde{F}_{k-1,\theta}^{N}(\check{u}_{0}^{i})=\Lambda_{0,\theta}^{C}(x_{*},\check{u}_{0}^{i})$.

  3. 3.

For $i\in\\{1,\dots,N\\}$, sample $u_{k}^{i}$ from
$\mathbb{W}(\cdot)\otimes\hat{p}_{\theta}^{l}(\check{x}_{k}^{i},\cdot)$. For
$i\in\\{1,\dots,N\\}$, compute:

$\widetilde{F}_{k,\theta}^{N}(u_{k}^{i})=\frac{\sum_{j=1}^{N}\exp\\{\Phi_{k,\theta}^{C}(\check{u}_{k-1}^{j},u_{k}^{i})\\}\hat{p}_{\theta}(\check{x}_{k}^{j},x_{k+1}^{i})\\{\widetilde{F}_{k-1,\theta}^{N}(\check{u}_{k-1}^{j})+\Lambda_{k,\theta}^{C}(\check{u}_{k-1}^{j},u_{k}^{i})\\}}{\sum_{j=1}^{N}\exp\\{\Phi_{k,\theta}^{C}(\check{u}_{k-1}^{j},u_{k}^{i})\\}\hat{p}_{\theta}(\check{x}_{k}^{j},x_{k+1}^{i})}.$

The estimate of $\nabla_{\theta}\log(\gamma_{k+1,\theta}(1))$ is:

$\widehat{\nabla_{\theta}\log(\gamma_{k+1,\theta}(1))}:=\frac{\sum_{i=1}^{N}\exp\\{\Phi_{k,\theta}^{C}(\check{u}_{k-1}^{i},u_{k}^{i})\\}\widetilde{F}_{k,\theta}^{N}(u_{k}^{i})}{\sum_{i=1}^{N}\exp\\{\Phi_{k,\theta}^{C}(\check{u}_{k-1}^{i},u_{k}^{i})\\}}.$
(28)

Set $k=k+1$ and return to the start of 2.

###  4.3 Time-Discretization

Whilst conceptually important, path-space valued Algorithm 2 can seldom be
implemented directly in practice; unbiased methods e.g. [4] may be possible,
but would be cumbersome. We develop a time-discretization procedure, in a
similar manner to that considered in Section 3.2.

We will discretize on the uniform grid with increment $\Delta_{l}=2^{-l}$. Let
$k\in\mathbb{N}_{0}$ and define:

$u_{k,l}=(z_{k+\Delta_{l}},z_{k+2\Delta_{l}},\dots,z_{k+1-\Delta_{l}},x_{k+1})\in(\mathbb{R}^{d_{x}})^{\Delta_{l}^{-1}}=E_{l},$

where $z_{k+\Delta_{l}},z_{k+2\Delta_{l}},\dots,z_{k+1-\Delta_{l}}$ will
represent increments of Brownian motion and $u_{-1,l}=x_{*}$. Define the
Markov kernel on $E_{l}$, for $k\in\mathbb{N}$:

$\widetilde{M}_{\theta}^{l}(u_{k-1,l},du_{k,l})=\Big{(}\prod_{s=1}^{\Delta_{l}^{-1}-1}\phi_{l}(z_{k+s\Delta_{l}})dz_{k+s\Delta_{l}}\Big{)}\hat{p}_{\theta}(x_{k},x_{k+1})dx_{k+1},$

where $\phi_{l}(z_{k+s\Delta_{l}})$ is the density associated to the
$\mathcal{N}_{d_{x}}(0,\Delta_{l}I_{d_{x}})$ distribution. We denote the
density of $\widetilde{M}_{\theta}^{l}$ as $\widetilde{Q}_{\theta}^{l}$. Set:

$\widetilde{\eta}_{0,\theta}^{l}(du_{0,l})=\Big{(}\prod_{s=1}^{\Delta_{l}^{-1}-1}\phi_{l}(z_{s\Delta_{l}})dz_{s\Delta_{l}}\Big{)}\hat{p}_{\theta}(x_{*},x_{1})dx_{1}.$

Now set, for
$(k,s)\in\\{0,1,\dots,T-1\\}\times\\{0,1,\dots,\Delta_{l}^{-1}-2\\}$:

$\displaystyle
X_{(s+1)\Delta_{l}+k}=X_{s\Delta_{l}+k}+b_{\theta}^{\circ}(s\Delta_{l},X_{s\Delta_{l}+k};x_{k+1})\Delta_{l}+\sigma(X_{s\Delta_{l}+k})Z_{(s+1)\Delta_{l}+k}.$
(29)

Define for $k\in\mathbb{N}_{0}$:

$\displaystyle\Psi_{\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle=\sum_{p=0}^{\Delta_{l}^{-1}-1}L_{\theta}(p\Delta_{l},x_{k+p\Delta_{l}})\Delta_{l};$
$\displaystyle J_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle=\sum_{p=0}^{\Delta_{l}^{-1}-1}h_{\theta}(x_{k+p\Delta_{l}})^{*}(Y_{k+(s+1)\Delta_{l}}-Y_{k+s\Delta_{l}})-\frac{1}{2}\sum_{p=0}^{\Delta_{l}^{-1}-1}h_{\theta}(x_{k+p\Delta_{l}})^{*}h_{\theta}(x_{k+p\Delta_{l}})\Delta_{l};$
$\displaystyle\Phi_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle=J_{k,\theta}^{l}(u_{k-1,l},u_{k,l})+\Psi_{\theta}^{l}(u_{k-1,l},u_{k,l})+\log\frac{\tilde{p}_{\theta,x_{k+1}}(x_{k},x_{k+1})}{\hat{{p}}_{\theta}(x_{k},x_{k+1})};$
$\displaystyle\widetilde{G}_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$
$\displaystyle=\exp\\{\Phi_{k,\theta}^{l}(u_{k-1,l},u_{k,l})\\}.$

Now, for $k\in\mathbb{N}_{0}$:

$\displaystyle\widetilde{\Lambda}_{k,\theta}^{l}($ $\displaystyle
u_{k-1,l},u_{k,l})=\sum_{p=0}^{\Delta_{l}^{-1}-1}(\nabla_{\theta}b_{\theta}(x_{k+p\Delta_{l}}))^{*}a(x_{k+p\Delta_{l}})^{-1}(x_{k+(p+1)\Delta_{l}}-x_{k+p\Delta_{l}}-b_{\theta}(x_{k+p\Delta_{l}})\Delta_{l})$
$\displaystyle+\sum_{p=0}^{\Delta_{l}^{-1}-1}(\nabla_{\theta}h_{\theta}(x_{k+p\Delta_{l}}))^{*}(Y_{k+(p+1)\Delta_{l}}-Y_{k+p\Delta_{l}})-\sum_{p=0}^{\Delta_{l}^{-1}-1}(\nabla_{\theta}h_{\theta}(x_{k+p\Delta_{l}}))^{*}h_{\theta}(x_{k+p\Delta_{l}})\Delta_{l}.$

Set:

$\widetilde{F}_{T,\theta}^{l}(u_{0,l},\dots,u_{T-1,l})=\sum_{k=0}^{T-1}\widetilde{\Lambda}_{k,\theta}^{l}(u_{k-1,l},u_{k,l}).$

Writing expectations w.r.t.
$\widetilde{\eta}_{\theta}^{l}(du_{0,l})\prod_{k=1}^{T-1}\widetilde{M}_{\theta}^{l}(u_{k-1,l},du_{k,l})$
as $\widetilde{\mathbb{E}}^{l}_{\theta}[\,\cdot\,|\mathcal{Y}_{T}]$, our
discretized approximation of $\nabla_{\theta}\log(\gamma_{T,\theta}(1))$ is:

$\nabla_{\theta}\log(\widetilde{\gamma}_{T,\theta}^{l}(1)):=\frac{\widetilde{\mathbb{E}}^{l}_{\theta}\,\big{[}\,\widetilde{F}_{T,\theta}^{l}(U_{0,l},\dots,U_{T-1,l})\prod_{k=0}^{T-1}\widetilde{G}_{k,\theta}^{l}(U_{k-1,l},U_{k,l})\,\big{|}\,\mathcal{Y}_{T}\,\big{]}}{\widetilde{\mathbb{E}}^{l}_{\theta}\,\big{[}\,\prod_{k=0}^{T-1}\widetilde{G}_{k,\theta}^{l}(U_{k-1,l},U_{k,l})\,\big{|}\,\mathcal{Y}_{T}\,\big{]}}.$

We note that, whilst terms
$\nabla_{\theta}\log(\widetilde{\gamma}_{T,\theta}^{l}(1))$,
$\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))$ should both converge to
$\nabla_{\theta}\log(\gamma_{T,\theta}(1))$, as $l\rightarrow\infty$, they
will in general be different for any fixed $l$.

One can also easily develop a discretized time reversal formula such as (13)
which will converge precisely to (26) as $l\rightarrow\infty$. Define, for
$k\in\mathbb{N}_{0}$:

$\widetilde{\pi}_{k,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{k,l})\big{)}:=\frac{\big{(}\prod_{p=0}^{k}\widetilde{G}_{p,\theta}^{l}(u_{p-1,l},u_{p,l})\big{)}\,\widetilde{\eta}_{0,\theta}^{l}(du_{0,l})\prod_{p=1}^{k}\widetilde{M}_{\theta}^{l}(u_{p-1,l},du_{p,l})}{\int_{E_{l}^{k+1}}\Big{(}\prod_{p=0}^{k}\widetilde{G}_{p,\theta}^{l}(u_{p-1,l},u_{p,l})\Big{)}\widetilde{\eta}_{0,\theta}^{l}(du_{0,l})\prod_{p=1}^{k}\widetilde{M}_{\theta}^{l}(u_{p-1,l},du_{p,l})}.$

Then, we have that:

$\nabla_{\theta}\log(\widetilde{\gamma}_{T,\theta}^{l}(1))=\int_{E_{l}^{T}}\widetilde{F}_{T,\theta}^{l}(u_{0,l},\dots,u_{T-1,l})\,\widetilde{\mathbb{Q}}_{T-1,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{T-1,l})\big{)},$
(30)

where we have defined:

$\displaystyle\widetilde{\mathbb{Q}}_{T-1,\theta}^{l}\big{(}d(u_{0,l},\dots,u_{T-1,l})\big{)}:=\widetilde{\pi}_{T-1,\theta}^{l}(du_{T-1,l})\prod_{k=1}^{T-1}\widetilde{B}_{k,\theta,\pi_{k-1,\theta}^{l}}^{l}(u_{k,l},du_{k-1,l})$

and:

$\displaystyle\widetilde{B}_{k,\theta,\widetilde{\pi}_{k-1,\theta}^{l}}^{l}(u_{k,l},$
$\displaystyle
du_{k-1,l}):=\frac{\widetilde{\pi}_{k-1,\theta}^{l}(du_{k-1,l})\widetilde{G}_{k,\theta}^{l}(u_{k-1,l},u_{k,l})\widetilde{Q}_{\theta}^{l}(u_{k-1,l},u_{k,l})}{\widetilde{\pi}_{k-1,\theta}^{l}(\widetilde{G}_{k,\theta}^{l}(\cdot,u_{k,l})\widetilde{Q}_{\theta}^{l}(\cdot,u_{k,l}))}.$

We remark that, due to the structure of the model:

$\widetilde{B}_{k,\theta,\widetilde{\pi}_{k-1,\theta}^{l}}^{l}(u_{k,l},du_{k-1,l})=\frac{\widetilde{\pi}_{k-1,\theta}^{l}(du_{k-1,l})\widetilde{G}_{k,\theta}^{l}(u_{k-1,l},u_{k,l})\hat{p}_{\theta}(x_{k},x_{k+1})}{\int_{E_{l}}\widetilde{\pi}_{k-1,\theta}^{l}(du_{k-1,l})\widetilde{G}_{k,\theta}^{l}(u_{k-1,l},u_{k,l})\hat{p}_{\theta}(x_{k},x_{k+1})}.$

######  Remark 4.4.

It is important to note that – in contrast to 3.2 – there is no cancellation
of terms of $\widetilde{G}_{k,\theta}^{l}(u_{k-1,l},u_{k,l})$ in the numerator
and denominator of this backward kernel. This is precisely due to recursion
(29) which leads to a path-dependence of the future coordinates of the
discretized bridge on the terminal position $x_{k+1}$.

###  4.4 Particle Approximation

Our online particle approximation of the gradient of the log-likelihood in
(30), for a given $l\in\mathbb{N}_{0}$ is presented in Algorithm 3. Our
estimates are given in (32) and (34) in Algorithm 3.

Algorithm 3 is simply the time-discretization of the procedure presented in
Algorithm 2. A number of remarks are again of interest. Firstly, the cost of
the algorithm per unit time is now $\mathcal{O}(N^{2}\Delta_{l}^{-1})$. The
increase in computational cost over Algorithm 1 is the fact that when
computing
$\widetilde{\Lambda}_{k,\theta}^{l}(\check{u}_{k-1,l}^{j},u_{k,l}^{i})$ in
(33), one must solve the recursion (29) for each
$(i,j)\in\\{1,\dots,N\\}^{2}$, which has a cost $\mathcal{O}(\Delta_{l}^{-1})$
and it is this cost that dominates. Secondly, following the discussion in
Section 3.3, we have proved in a companion work that, under appropriate
assumptions, the MSE for $(k,N)\in\mathbb{N}^{2}$:

$\mathbb{E}_{\theta}\Big{[}\,\Big{\|}\,\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{l}(1))}-\nabla_{\theta}\log(\gamma_{k,\theta}(1))\,\Big{\|}_{2}^{2}\,\Big{]}\leq
C\Big{(}\frac{1}{N}+\Delta_{l}\Big{)},$ (31)

for constant $C$ that does not depend on $N$, $l$. To achieve an MSE of
$\mathcal{O}(\epsilon^{2})$ for some $\epsilon>0$ given, one sets
$l=\mathcal{O}(|\log(\epsilon)|)$ and $N=\mathcal{O}(\epsilon^{-2})$. The cost
per unit time of doing this, is then $\mathcal{O}(\epsilon^{-6})$. This is
significantly worse than the approach in Algorithm 1, but we remark that when
discussing the cost of Algorithm 1, in the bound (18), we have assumed that
the constant $C$ does not depend upon $l$. However, in a sequel work we will
show that under assumptions that this afore-mentioned $C$ explodes
exponentially in $l$. Conversely, $C$ in (31) can be proved to be independent
of $l$, precisely due to the path-space development we have adopted. We
remark, however, that one can use an MLMC method to reduce this cost of
$\mathcal{O}(\epsilon^{-6})$ per unit time of Algorithm 3 and this algorithm
is presented in the next section.

Algorithm 3 Modified Online Score Function Estimation for a given
$l\in\mathbb{N}_{0}$.

  1. 1.

For $i\in\\{1,\dots,N\\}$, sample $u_{0,l}^{i}$ i.i.d. from
$\widetilde{\eta}_{0,\theta}^{l}(\cdot)$. The estimate of
$\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l}(1))$ is:

$\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l}(1))}:=\frac{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})\widetilde{\Lambda}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}.$
(32)

Set $k=1$ and for $i\in\\{1,\dots,N\\}$, $\check{u}_{-1,l}^{i}=x_{*}$.

  2. 2.

For $i\in\\{1,\dots,N\\}$ sample $\check{u}_{k-1,l}^{i}$ from:

$\sum_{i=1}^{N}\frac{\widetilde{G}_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{i},u_{k-1,l}^{i})}{\sum_{j=1}^{N}\widetilde{G}_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{j},u_{k-1,l}^{j})}\delta_{\\{u_{k-1,l}^{i}\\}}(\cdot).$

If $k=1$, for $i\in\\{1,\dots,N\\}$, set
$\widetilde{F}_{k-1,\theta}^{l,N}(\check{u}_{0,l}^{i})=\widetilde{\Lambda}_{0,\theta}^{l}(x_{*},\check{u}_{0,l}^{i})$.

  3. 3.

For $i\in\\{1,\dots,N\\}$, sample $u_{k,l}^{i}$ from
$\widetilde{M}_{\theta}^{l}(\check{u}_{k-1,l}^{i},\cdot)$. For
$i\in\\{1,\dots,N\\}$, compute:

$\widetilde{F}_{k,\theta}^{l,N}(u_{k,l}^{i})=\frac{\sum_{j=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{j},u_{k,l}^{i})\hat{p}_{\theta}(\check{x}_{k}^{j},x_{k+1}^{i})\\{\widetilde{F}_{k-1,\theta}^{l,N}(\check{u}_{k-1,l}^{j})+\widetilde{\Lambda}_{k,\theta}^{l}(\check{u}_{k-1,l}^{j},u_{k,l}^{i})\\}}{\sum_{j=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{j},u_{k,l}^{i})\hat{p}_{\theta}(\check{x}_{k}^{j},x_{k+1}^{i})}.$
(33)

The estimate of $\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l}(1))$
is:

$\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l}(1))}=\frac{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})\widetilde{F}_{k,\theta}^{l,N}(u_{k,l}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})}.$
(34)

Set $k=k+1$ and return to the start of 2..

###  4.5 Multilevel Particle Filter

We now present a new multilevel particle filter along with online estimation
of the score-function. We fix $l\in\mathbb{N}$ for now and for
$(k,s)\in\mathbb{N}_{0}\times\\{l,l-1\\}$ define:

$u_{k,s}=(z_{k+\Delta_{l},s},z_{k+2\Delta_{l},s},\dots,z_{k+1-\Delta_{l},s},x_{k+1,s})\in(\mathbb{R}^{d_{x}})^{\Delta_{s}^{-1}}=E_{s}.$

Algorithm 4 Coupled Online Score Function Estimation for a given
$l\in\mathbb{N}$.

  1. 1.

     * •

For $i\in\\{1,\dots,N\\}$, sample $u_{0,l}^{i}$ i.i.d. from
$\widetilde{\eta}_{0,\theta}^{l}(\cdot)$.

     * •

For $(i,p)\in\\{1,\dots,N\\}\times\\{1,\dots,\Delta_{l-1}^{-1}-1\\}$, set
$z_{p\Delta_{l-1},l-1}^{i}=z_{p\Delta_{l-1},l}^{i}+z_{p\Delta_{l-1}-\Delta_{l},l}^{i}$
and $x_{1,l-1}^{i}=x_{1,l}^{i}$.

The estimate of
$\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l}(1))-\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l-1}(1))$
is:

$\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l}(1))}-\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{1,\theta}^{l-1}(1))}:\vspace{-0.3cm}=$
$\frac{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})\widetilde{\Lambda}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l}(x_{*},u_{0,l}^{i})}-\frac{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l-1}(x_{*},u_{0,l-1}^{i})\widetilde{\Lambda}_{0,\theta}^{l-1}(x_{*},u_{0,l-1}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{0,\theta}^{l-1}(x_{*},u_{0,l-1}^{i})}.$
(35)

Set $k=1$ and for $i\in\\{1,\dots,N\\}$,
$\check{u}_{-1,l}^{i}=\check{u}_{-1,l-1}^{i}=x_{*}$.

  2. 2.

For $i\in\\{1,\dots,N\\}$, sample
$(\alpha_{l}(i),\alpha_{l-1}(i))\in\\{1,\dots,N\\}^{2}$ from a coupling of:

$\sum_{i=1}^{N}\frac{\widetilde{G}_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{\alpha_{l}(i)},u_{k-1,l}^{\alpha_{l}(i)})}{\sum_{j=1}^{N}\widetilde{G}_{k-1,\theta}^{l}(\check{u}_{k-2,l}^{j},u_{k-1,l}^{j})}\quad\quad\textrm{and}\quad\quad\sum_{i=1}^{N}\frac{\widetilde{G}_{k-1,\theta}^{l-1}(\check{u}_{k-2,l-1}^{\alpha_{l-1}(i)},u_{k-1,l}^{\alpha_{l-1}(i)})}{\sum_{j=1}^{N}\widetilde{G}_{k-1,\theta}^{l-1}(\check{u}_{k-2,l-1}^{j},u_{k-1,l-1}^{j})},$

and set
$(\check{u}_{k-1,l}^{i},\check{u}_{k-1,l-1}^{i})=(u_{k-1,l}^{\alpha_{l}(i)},u_{k-1,l-1}^{\alpha_{l-1}(i)})$.
If $k=1$, for $i\in\\{1,\dots,N\\}$, $s\in\\{l,l-1\\}$, set
$\widetilde{F}_{k-1,\theta}^{s,N}(\check{u}_{0,s}^{i})=\widetilde{\Lambda}_{0,\theta}^{s}(x_{*},\check{u}_{0,s}^{i})$.

  3. 3.

     * •

For $i\in\\{1,\dots,N\\}$, sample $(x_{k+1,l}^{i},x_{k+1,l-1}^{i})$ from a
coupling of $\hat{p}_{\theta}(\check{x}_{k,l}^{i},\cdot)$ and
$\hat{p}_{\theta}(\check{x}_{k,l-1}^{i},\cdot)$.

     * •

For $i\in\\{1,\dots,N\\}$ sample
$z_{k+\Delta_{l},l}^{i},\dots,z_{k+1-\Delta_{l},l}^{i}$ i.i.d. from
$\prod_{s=1}^{\Delta_{l}^{-1}-1}\phi_{l}(\,\cdot\,)$.

     * •

For $(i,p)\in\\{1,\dots,N\\}\times\\{1,\dots,\Delta_{l-1}^{-1}-1\\}$ set
$z_{k+p\Delta_{l-1},l-1}^{i}=z_{k+p\Delta_{l-1},l}^{i}+z_{k+p\Delta_{l-1}-\Delta_{l},l}^{i}$.

     * •

For $(i,s)\in\\{1,\dots,N\\}\times\\{l.l-1\\}$, compute:

$\widetilde{F}_{k,\theta}^{s,N}(u_{k,l}^{i})=\frac{\sum_{j=1}^{N}\widetilde{G}_{k,\theta}^{s}(\check{u}_{k-1,s}^{j},u_{k,s}^{i})\hat{p}_{\theta}(\check{x}_{k,s}^{j},x_{k+1,s}^{i})\\{\widetilde{F}_{k-1,\theta}^{s,N}(\check{u}_{k-1,s}^{j})+\widetilde{\Lambda}_{k,\theta}^{s}(\check{u}_{k-1,s}^{j},u_{k,s}^{i})\\}}{\sum_{j=1}^{N}\widetilde{G}_{k,\theta}^{s}(\check{u}_{k-1,s}^{j},u_{k,s}^{i})\hat{p}_{\theta}(\check{x}_{k,s}^{j},x_{k+1,s}^{i})}.$

The estimate of
$\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l}(1))-\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l-1}(1))$
is:

$\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l}(1))}-\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k+1,\theta}^{l-1}(1))}:=\vspace{-0.3cm}$
$\frac{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})\widetilde{F}_{k,\theta}^{l,N}(u_{k,l}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l}(\check{u}_{k-1,l}^{i},u_{k,l}^{i})}-\frac{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l-1}(\check{u}_{k-1,l-1}^{i},u_{k,l-1}^{i})\widetilde{F}_{k,\theta}^{l-1,N}(u_{k,{l-1}}^{i})}{\sum_{i=1}^{N}\widetilde{G}_{k,\theta}^{l-1}(\check{u}_{k-1,l-1}^{i},u_{k,l-1}^{i})}.$
(36)

Set $k=k+1$ and return to the start of 2..

We give the approach in Algorithm 4. Before explaining how one can use
Algorithm 4 to provide online estimates of the score function, several remarks
are required to continue. The first is related to the couplings mentioned in
Algorithm 4 point 2. and point 3. bullet 1. The coupling in point 2, requires
a way to resample the indices of the particles so that they have the correct
marginals. This topic has been investigated considerably in the literature,
see e.g. [18, 29], and techniques that have been adopted include sampling
maximal coupling, e.g. [16], or using the $\mathbb{L}_{2}$-Wasserstein optimal
coupling [2]; in general the latter is found to be better in terms of variance
reduction, but can only be implemented when $d_{x}=1$. We rely upon the
maximal coupling in this paper, which has a cost of $\mathcal{O}(N)$ per unit
time. For point 3. bullet 1, one again has a considerable degree of
flexibility. In this article we sample the maximal coupling which can be
achieved at a cost which is at most $\mathcal{O}(N)$ cost per-unit time using
the algorithm of [32]. The second main remark of interest is that the basic
filter that is sampled in Algorithm 4 is an entirely new coupled particle
filter for diffusions (i.e. different to [16, 19]). The utility of the
approach relative to [19] is of great interest, in the context of filtering.

Set $(l_{*},L)\in\mathbb{N}^{2}$ with $l_{*}<L$. The idea is to run Algorithm
4, independently, for $l\in\\{l_{*},\dots,L\\}$ each with $N_{l}$ particles
and, independently, Algorithm 3 for $l=l_{*}-1$ with $N_{l_{*}-1}$ particles.
We then consider the estimate, for $k\in\mathbb{N}$

$\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{L}(1))}_{ML}:=\sum_{l=l_{*}}^{L}\Big{\\{}\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{l}(1))}-\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{l-1}(1))}\Big{\\}}+\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{l_{*}-1}(1))},$

where the summands on the right hand side are defined in (35) and (36) and the
last term on the right hand side is as either (32) or (34) (depending on $k$).
Now, we show in an on-going companion work, under appropriate assumptions, one
has the following result for $(k,N_{l_{*}-1:L})\in\mathbb{N}^{L-l_{*}+2}$:

$\mathbb{E}_{\theta}\Big{[}\,\Big{\|}\,\widehat{\nabla_{\theta}\log(\widetilde{\gamma}_{k,\theta}^{L}(1))}_{ML}-\nabla_{\theta}\log(\gamma_{k,\theta}(1))\,\Big{\|}_{2}^{2}\,\Big{]}\leq
C\Big{(}\sum_{l=l_{*}-1}^{L}\frac{\Delta_{l}^{\beta}}{N_{l}}+\Delta_{L}\Big{)},$
(37)

for constant $C$ that does not depend on $N$, $l$; also, $\beta=1$ if $\sigma$
is a constant function, else $\beta=1/2$. Choose: i) $L$ so that
$\Delta_{L}=\mathcal{O}(\epsilon^{2})$, for $\epsilon>0$ given; ii) if
$\beta=1$, $N_{l}=\mathcal{O}(\epsilon^{-2}\Delta_{l}^{1/2+\rho})$ for some
$0<\rho<1/2$. These selections yield an MSE of $\mathcal{O}(\epsilon^{2})$ for
a cost of $\mathcal{O}(\epsilon^{-4})$. If $\beta=1/2$, one can set
$N_{l}=\mathcal{O}(\epsilon^{-2}\Delta_{l}^{1/2+\rho}\Delta_{L}^{-\rho})$ for
some $\rho>0$. This will yield an MSE of $\mathcal{O}(\epsilon^{2})$ for a
cost of $\mathcal{O}(\epsilon^{-4(1+\rho)})$. Such results are at least as
good as the method in Section 3.3, assuming that latter approach does not
collapse with $l$.

We remark that it is possible to produce an almost-surely unbiased estimator
of the score function, when $\theta$ is the true parameter, using a
combination of the multilevel method that has been developed here and the
approach in [17]. This is left for future work.

##  5 Numerical Results

In this section, we consider four models to investigate the various properties
of our algorithms. The score function is estimated using both Algorithms 1 and
3 for a fixed $\theta$. We will show, as expected, that they are equivalent
for a large number of particles $N$ and a high level of descritization $l$. We
then compare the cost of Algorithm 3 and its multilevel version Algorithm 4.
As an application of our methods, we use Algorithms 1 and 4 for parameter
estimation via stochastic gradient. The code is written in MATLAB and it can
be downloaded from https://github.com/ruzayqat/score_based_par_est.

We remark that we will not use Algorithm 3 for parameter estimation because it
is ‘slow’ compared to the algorithms as illustrated in the previous section
and in Figure 1. In Figure 1 we consider the Model 1, as described in the next
section, with $T=40$, $l=10$, $\theta=(-0.4,-0.5)$, $\kappa=2$ and
$x_{*}=0.2$. Figure 1 provides a comparison between the cost of Algorithms 1
and 3, which is the average machine time measured in seconds needed per each
simulation, versus the number of particles $N$. As predicted by our
theoretical conjectures, we see that the cost of Algorithm 1 is significantly
lower than that of Algorithm 3.

Figure 1: Comparison between the cost of Algorithms 1 and 3 per each
simulation versus the number of particles $N$. We run both algorithms on Model
1.

###  5.1 Models

In the following, parameters $(\kappa,\sigma)$ are fixed.

Model 1: Let $d_{x}=d_{y}=1$, $d_{\theta}=2$ and consider the following linear
SDEs:

$\displaystyle dX_{t}$ $\displaystyle=\theta_{1}X_{t}dt+\sigma dW_{t};$
$\displaystyle dY_{t}$ $\displaystyle=\theta_{2}(\kappa-X_{t})dt+dB_{t}.$

Model 2: Let $d_{x}=d_{y}=1$, $d_{\theta}=3$ and consider a nonlinear
diffusion process along with a linear diffusion process of observations:

$\displaystyle dX_{t}$
$\displaystyle=\big{(}\tfrac{\theta_{1}}{X_{t}}+\theta_{2}X_{t}\big{)}dt+\sigma
dW_{t};$ $\displaystyle dY_{t}$ $\displaystyle=\theta_{3}(\kappa-
X_{t})dt+dB_{t}.$

Model 3: Let $d_{x}=d_{y}=1$, $d_{\theta}=3$ and consider a nonlinear signal
along with a nonlinear diffusion process of observations. The first SDE is a
Cox-Ingersoll-Ross process after an 1-1 transform. Thus:

$\displaystyle dX_{t}$
$\displaystyle=\frac{1}{2}\big{(}\tfrac{\theta_{1}\theta_{2}-\sigma^{2}}{X_{t}}-\theta_{2}X_{t}\big{)}dt+\sigma
dW_{t};$ $\displaystyle dY_{t}$ $\displaystyle=\theta_{3}(\kappa-
X_{t}^{2})dt+dB_{t}.$

This model has a solution if and only if $\theta_{1}\theta_{2}>2\sigma^{2}$.

Model 4: Let $d_{x}=d_{y}=1$, $d_{\theta}=3$ and consider a type of Black-
Scholes model with a stochastic volatility:

$\displaystyle dX_{t}$
$\displaystyle=\theta_{1}(\theta_{2}-X_{t})dt+\sigma(X_{t})dW_{t};$
$\displaystyle dY_{t}$
$\displaystyle=(\theta_{3}-\frac{1}{2}X_{t}^{2})dt+dB_{t}.$

where $\sigma(X_{t})=\beta/\sqrt{X_{t}^{2}+1}$ and $\beta$ is fixed. In the
hidden process, $\theta_{1}$ and $\theta_{2}$ are the speed and level of mean
reversion and $\theta_{3}$ is a mean type level for the observation process,
We will apply our methodology (see Figure 8 later on) on the log mid-price of
Tesla Inc. stock in 2018. The dataset shown in Figure 2 represents the log
mid-price at every second during a trading day for a total of 250 trading
days.

Figure 2: The log of mid-price for every second of Tesla stock in 2018. The
dataset contains $5.85\times 10^{6}$ data points.

###  5.2 Simulation Results

In all our results data are generated from the model under the finest
discretization considered except in model 4, where we use a real data. In
Algorithm 3, we consider the auxiliary linear process following:

$\displaystyle d\tilde{X}_{t}=\sigma dW_{t}$

in models 1-3 and in model 4 it follows:

$\displaystyle d\tilde{X}_{t}=\sigma(x^{\prime})dW_{t}$

In models 1-3,
$\tilde{p}_{\theta,x^{\prime}}(x,t;x^{\prime},1)=\mathcal{N}(x^{\prime};x,(1-t)\sigma^{2})$,
hence
$\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})=\mathcal{N}(x^{\prime};x,\sigma^{2})$,
which is easy to sample $x^{\prime}$ from, and therefore,
$\hat{p}_{\theta}(x,x^{\prime})=\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})$.
But in model 4,
$\tilde{p}_{\theta,x^{\prime}}(x,x^{\prime})=\mathcal{N}(x^{\prime};x,\sigma^{2}(x^{\prime}))$
which is not easy to sample $x^{\prime}$ from. Therefore, we take
$\hat{p}_{\theta}(x,x^{\prime})=\mathcal{N}(x^{\prime};x,\sigma^{2}(x))$.

####  5.2.1 Estimation of the Score Function

For each model, we fix parameter $\theta$ and estimate the score function
using Algorithms 1 and 3. In Algorithm 1, $N\in\\{3000,7000,4000,5000\\}$ in
the 1st, 2nd, 3rd & 4th models, respectively. In Algorithm 3,
$N\in\\{1000,2000,1000,1500\\}$ in the 1st, 2nd, 3rd & 4th models,
respectively. In both algorithms, we set the discretization level to $l=10$.
In Models 1, 2 and 3, we set $\kappa=2$, $2.2$, $1.5$, $x_{*}=0.2$, $1$, $2$
and $\sigma=0.3$, $0.25$, $0.25$, respectively. While in model 4, we set
$x_{*}=1.3$ and $\beta=2$; $T=50$ for all 4 models ($T=50$ in model 4
corresponds to 14.22 hours of trading).

Figure 3 summarizes the results of 56 replications of estimates of the score
function for each model and for each unit time point. These simulations are
implemented in parallel using 8 CPUs. The figure illustrates that both
algorithms are equivalent for large $N$ and $l$ as one would expect.

(a)  (b)

(c)  (d)

Figure 3: Trajectories from the execution of Algorithms 1 and 3 for the
estimation of the score function in Models 1-4.

####  5.2.2 Cost Comparison

We now consider comparing the costs of Algorithms 3 and 4. We take $l_{*}$ to
be 7 in the models 1, 2 & 3 and 8 in the 4th model. The parameters of the
model are as in the previous section. The ground truth is computed at level 11
with $N=2000$ using Algorithm 3. We run 56 simulations of both algorithms for
$L\in\\{l_{*}-1,\cdots,10\\}$. For each $L$, the number of particles are
carefully chosen to give similar MSE values from both algorithms.
Particularly, the number of particles in Algorithm 3 is $N_{L}=\left\lfloor
C_{1}2^{L}\right\rfloor$ and for each level $l$ in Algorithm 3. In Figure 4
the number of particles is $N_{l}=\left\lfloor
C_{2}2^{L}(L-l_{*}+2)\Delta_{l}^{1/2+\rho}\right\rfloor$ (in models 1 to 3)
and $N_{l}=\left\lfloor
C_{2}2^{L}(L-l_{*}+2)\Delta_{l}^{1/2+\rho}\Delta_{L}^{-\rho}\right\rfloor$ (in
model 4), where $C_{1}$ and $C_{2}$ are constants. In we can observe the cost
against MSE curve, that appear to follow our conjectures over algorithmic
costs earlier in the article.

(a)  (b)

(c)  (d)

Figure 4: Cost per each simulation versus MSE on a log-log scale for
Algorithms 3 and 4. The dashed lines are for reference.

####  5.2.3 Parameter Estimation

We use Algorithms 1, 4 to estimate the parameters in each model. In Algorithm
1, the level of discretization, $l$, is 10 for models 1-3 and 9 for model 4,
and the number of particles, $N$, is 2,000 for models 1-3 and 2500 for model
4. In Algorithm 4, we use $l_{*}=7$, $L=10$ and the number of particles on
each level $l\in\\{l_{*}-1,\cdots,L\\}$ is
$N_{l}=2^{L}(L-l_{*}+2)\Delta_{l}^{1/2+\rho}$, where
$\rho\in\\{0.14,\,0.09,\,0.11\\}$ in Models 1, 2 and 3, respectively. In model
4, $l_{*}=8$, $L=9$ and the number of particles on each level
$l\in\\{l_{*}-1,\cdots,L\\}$ is $N_{l}=1.4\times
2^{L}(L-l_{*}+2)\Delta_{l}^{1/2+\rho}\Delta_{L}^{-\rho}$ where $\rho=0.1$.

Figure 5 considers Model 1. We fix $x_{*}=0.2$, $\sigma=0.3$, $\kappa=2$,
$T=20,000$. The parameter values used to generate the data are
$(\theta_{1}^{\star},\theta_{2}^{\star})=(-0.7,-0.5)$. For the stochastic
gradient algorithm, we used an initial value $(-0.05,-1.5)$ and step-size
$\alpha_{k}=k^{-0.85}$. Figure 6 considers Model 2. We fix $x_{*}=1.8$,
$\sigma=0.25$, $\kappa=2.2$, $T=20,000$. The parameter values used to generate
the data are
$(\theta_{1}^{\star},\theta_{2}^{\star},\theta_{3}^{\star})=(1.3,-0.5,0.18)$.
For the stochastic gradient algorithm, we used initial value $(0.8,-1,0.8)$
and step-size $\alpha_{k}=k^{-0.95}$. Figure 7 considers Model 3. We fix
$x_{*}=1.5$, $\sigma=0.25$, $\kappa=2$, $T=20,000$. The parameter values used
to generate the data are
$(\theta_{1}^{\star},\theta_{2}^{\star},\theta_{3}^{\star})=(2,1,0.45)$. For
the stochastic gradient algorithm, we used an initial value $(1.24,0.6,1.11)$
and step size $\alpha_{k}=k^{-0.9}$. Figure 8 considers Model 4 applied to the
data in Figure 2. We fix $x_{*}=1.3$, $\beta=2$, $T=11425$ (there is a
rescaling of the time parameter). For the stochastic gradient algorithm, we
used an initial value $(2.4,0.5,0.4)$ and step size $\alpha_{k}=k^{-0.82}$. In
all cases considered (Figures 5-8) our selected settings allow for an accurate
estimation of the parameter values over long time periods.

(a)  (b)

(c)  (d)

Figure 5: Trajectories from the execution of Algorithm 1 (left panel) and
Algorithm 4 (right panel) for the estimation of $(\theta_{1},\theta_{2})$ from
Model 1. The horizontal dashed lines in the plots show the true parameter
values $(\theta_{1}^{\star},\theta_{2}^{\star})=(-0.7,-0.5)$.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 6: Trajectories from the execution of Algorithm 1 (left panel) and
Algorithm 4 (right panel) for the estimation of
$(\theta_{1},\theta_{2},\theta_{3})$ of the second model. We used an initial
value $(0.8,-1,0.8)$. The horizontal dashed lines in the plots show the true
parameter values
$(\theta_{1}^{\star},\theta_{2}^{\star},\theta_{3}^{\star})=(1.3,-0.5,0.18)$.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 7: Trajectories from the execution of Algorithm 1 (left panel) and
Algorithm 4 (right panel) for the estimation of
$(\theta_{1},\theta_{2},\theta_{3})$ from Model 3. The horizontal dashed lines
in the plots show the true parameter values
$(\theta_{1}^{\star},\theta_{2}^{\star},\theta_{3}^{\star})=(2,1,0.45)$.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 8: Trajectories from the execution of Algorithm 1 (left panel) and
Algorithm 4 (right panel) for the estimation of
$(\theta_{1},\theta_{2},\theta_{3})$ from Model 4.

#### Acknowledgements

AJ & HR were supported by KAUST baseline funding. AB acknowledges support from
a Leverhulme Trust Prize. DC was partially supported by EU Synergy project
STUOD - DLV-856408. NK acknowledges funding by a JP Morgan A.I. Faculty award.

##  Appendix A Derivation of (4)

Recall that under $\overline{\mathbb{P}}_{\theta}$ the original processes
(2.1)-(2.2) have dynamics:

$\displaystyle dY_{t}=dB_{t};\quad
dX_{t}=b_{\theta}(X_{t})dt+\sigma(X_{t})dW_{t}.$

We consider the processes:

$\displaystyle dY_{t}=dB_{t};\quad dX_{t}=\sigma(X_{t})dW_{t},$

so that if $\mathbb{P}$ denotes their law, then we have the Radon-Nikodym
derivative:

$\displaystyle\frac{d\overline{\mathbb{P}}_{\theta}}{d\mathbb{P}}=\exp\Big{\\{}\int_{0}^{T}b_{\theta}(X_{s})^{*}a(X_{s})^{-1}dX_{s}-\tfrac{1}{2}\int_{0}^{T}b_{\theta}(X_{s})^{*}a(X_{s})^{-1}b_{\theta}(X_{s})ds\Big{\\}}.$

The log-likelihood is
$\log(\gamma_{T,\theta}(1))=\log\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}|\mathcal{Y}_{T}\right]$,
with $Z_{T,\theta}=d\mathbb{P}_{\theta}/d\overline{\mathbb{P}}_{\theta}$,
thus:

$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}(1))=\frac{1}{\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}|\mathcal{Y}_{T}\right]}\nabla_{\theta}\overline{\mathbb{E}}_{\theta}\big{[}d\mathbb{P}_{\theta}/d\overline{\mathbb{P}}_{\theta}\big{|}\mathcal{Y}_{T}\big{]}.$

For convenience we use the notation $\overline{\mathbb{E}}_{\theta,X}$,
$\mathbb{E}_{X}$ for the marginal expectations w.r.t. the original process $X$
and the $\theta$-free process $X$ defined above, respectively. Notice that we
can write:

$\displaystyle\nabla_{\theta}\overline{\mathbb{E}}_{\theta}\big{[}d\mathbb{P}_{\theta}/d\overline{\mathbb{P}}_{\theta}\big{|}\mathcal{Y}_{T}\big{]}$
$\displaystyle=\nabla_{\theta}\overline{\mathbb{E}}_{\theta,X}\big{[}d\mathbb{P}_{\theta}/d\overline{\mathbb{P}}_{\theta}\big{]}$
$\displaystyle=\nabla_{\theta}\mathbb{E}_{X}\big{[}(d\mathbb{P}_{\theta}/d\overline{\mathbb{P}}_{\theta})(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P})\big{]}$
$\displaystyle=\mathbb{E}_{X}\big{[}Z_{T,\theta}\nabla_{\theta}\log(Z_{T,\theta})(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P})\big{]}+\mathbb{E}_{X}\big{[}Z_{T,\theta}\nabla_{\theta}\log(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P})(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P})\big{]}$
$\displaystyle=\overline{\mathbb{E}}_{\theta}[(\nabla_{\theta}\log(Z_{T,\theta})+\nabla_{\theta}\log(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P}))Z_{T,\theta}].$

One can now verify that, for $\lambda_{T,\theta}$ as defined in the main text:

$\displaystyle\nabla_{\theta}\log(Z_{T,\theta})+\nabla_{\theta}\log(d\overline{\mathbb{P}}_{\theta}/d\mathbb{P})=\lambda_{T,\theta}.$

##  Appendix B $\mathbb{L}_{r}$ Bound for the Discretization Error

###  B.1 Formulation

We now consider proving a bound on ($\|\cdot\|_{2}$ is the $L_{2}-$norm for
vectors)

$\mathbb{E}\left[\left\|\nabla_{\theta}\log(\gamma_{T,\theta}(1))-\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))\right\|_{2}^{r}\right]^{1/r}$

where

$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}(1))$ $\displaystyle=$
$\displaystyle\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}$
$\displaystyle\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))$ $\displaystyle=$
$\displaystyle\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}.$

We begin by noting that for any random vector $X$ of dimension $d$ and finite
$r-$moments that

$\mathbb{E}[\|X\|_{2}^{r}]^{1/r}\leq\left(\sum_{i=1}^{d}\mathbb{E}[|X^{(i)}|^{r}]^{2/r}\right)^{1/2}.$

As a result, it will suffice to control for each
$i\in\\{1,\dots,d_{\theta}\\}$

$\mathbb{E}_{\theta}\left[\left|\nabla_{\theta}\log(\gamma_{T,\theta}(1))^{(i)}-\nabla_{\theta}\log(\gamma_{T,\theta}^{l}(1))^{(i)}\right|^{r}\right]^{1/r}.$

Throughout all of our proofs, $C$ is a deterministic constant whose value will
change upon each appearance. In addition we supress any dependencies on
$\theta$ below.

###  B.2 Technical Results

######  Lemma B.1.

Assume (D(D1)). Then for any $(r,T)\in[1,\infty)\times\mathbb{N}$ there exists
a $C<+\infty$ such that

$\overline{\mathbb{E}}_{\theta}\left[\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]^{r}}\right]\leq
C.$

###### Proof.

By the conditional Jensen inequality

$\overline{\mathbb{E}}_{\theta}\left[\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]^{r}}\right]\leq\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}^{-r}\right].$

It is now simple to use the properties of the process under study to deduce
the result. ∎

######  Remark B.1.

A standard result is that
$\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}^{r}\right]\leq C$ for any
fixed $r\in\mathbb{R}$.

######  Lemma B.2.

Assume (D(D1)). Then for any $(r,T)\in[1,\infty)\times\mathbb{N}$ there exists
a $C<+\infty$ such that for any $l\in\mathbb{N}_{0}$

$\overline{\mathbb{E}}_{\theta}\left[\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]^{r}}\right]\leq
C.$

###### Proof.

By the conditional Jensen inequality

$\overline{\mathbb{E}}_{\theta}\left[\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]^{r}}\right]\leq\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{-r}\right].$

The result now follows by the arguments stated in [19, eq. (20)-(21)]. ∎

######  Remark B.2.

By the arguments stated in [19, eq. (20)-(21)]
$\overline{\mathbb{E}}_{\theta}\left[Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{r}\right]\leq
C$ for any fixed $r\in\mathbb{R}$ and $C$ does not depend on $l$.

######  Lemma B.3.

Assume (D(D1)). Then for any $(r,T)\in[1,\infty)\times\mathbb{N}$ there exists
a $C<+\infty$ such that for any
$(l,i)\in\mathbb{N}_{0}\times\\{1,\dots,d_{\theta}\\}$

$\overline{\mathbb{E}}_{\theta}\left[\left|\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]\right|^{r}\right]\leq
C.$

###### Proof.

Using the conditional Jensen inequality, Cauchy-Schwarz and Remark B.2 it
suffices to bound

$\overline{\mathbb{E}}_{\theta}\left[\left|\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\right|^{2r}\right].$

Via the $C_{2r}-inequality$ one can then simply focus on the 3 terms

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\sum_{k=0}^{T/\Delta_{l}-1}\,\left(\frac{\partial}{\partial\theta_{i}}b_{\theta}(x_{k\Delta_{l}})\right)^{*}a(x_{k\Delta_{l}})^{-1}\sigma(x_{k\Delta_{l}})(W_{(k+1)\Delta_{l}}-W_{k\Delta_{l}})\right|^{2r}\right]$
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\sum_{k=0}^{T/\Delta_{l}-1}\left(\frac{\partial}{\partial\theta_{i}}h_{\theta}(x_{k\Delta_{l}})\right)^{*}(Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}})\right|^{2r}\right]$
$\displaystyle T_{3}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\sum_{k=0}^{T/\Delta_{l}-1}\left(\frac{\partial}{\partial\theta_{i}}h_{\theta}(x_{k\Delta_{l}})\right)^{*}h_{\theta}(x_{k\Delta_{l}})\Delta_{l}\right|^{2r}\right].$

To bound $T_{1}$ and $T_{2}$ one can simply apply the Burkholder-Gundy-Davis
inequality and combine this with the boundedness of the terms which are
functions of $x$; this is a standard argument in the literature. The bound on
$T_{3}$ is immediate by the boundedness of the summands. This concludes the
proof. ∎

######  Remark B.3.

It is more-or-less the same argument as in the proof of Lemma B.3 to deduce
that
$\overline{\mathbb{E}}_{\theta}\left[\left|\lambda_{T,\theta}^{(i)}\right|^{r}\right]\leq
C$.

######  Lemma B.4.

Assume (D(D1)). Then for any $(r,T)\in[1,\infty)\times\mathbb{N}$ there exists
a $C<+\infty$ such that for any
$(l,i)\in\mathbb{N}_{0}\times\\{1,\dots,d_{\theta}\\}$

$\overline{\mathbb{E}}_{\theta}\left[\left|\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{(i)}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]-\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]\right|^{r}\right]^{1/r}\leq
C\Delta_{l}^{1/2}.$

###### Proof.

Applying conditional Jensen and Minkowski we have the upper-bound

$\overline{\mathbb{E}}_{\theta}\left[\left|\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{(i)}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]-\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]\right|^{r}\right]^{1/r}\leq
T_{1}+T_{2}$

where

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\left\\{\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}-\lambda_{T,\theta}^{(i)}\right\\}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\right|^{r}\right]^{1/r}$
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\left\\{Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})-Z_{T,\theta}\right\\}\lambda_{T,\theta}^{(i)}\right|^{r}\right]^{1/r}.$

For $T_{1}$ one can use Cauchy-Schwarz and the result in Remark B.2 to deduce
the upper-bound

$T_{1}\leq
C\overline{\mathbb{E}}_{\theta}\left[\left|\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}-\lambda_{T,\theta}^{(i)}\right|^{2r}\right]^{1/(2r)}$

Then the term on the R.H.S. can be dealt with by using standard results in the
discretization of Riemann-integrals coupled with the Burkholder-Gundy-Davis
inequality, Lipschitz properties of the various functions and Euler
discretizations. As the results are almost identical to the calculations in
[9, pp.589] they are omitted. That is, one can deduce that

$T_{1}\leq C\Delta_{l}^{1/2}.$

For $T_{2}$, again, using Cauchy-Schwarz and the result in Remark B.3 we have

$T_{2}\leq
C\overline{\mathbb{E}}_{\theta}\left[\left|Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})-Z_{T,\theta}\right|^{2r}\right]^{1/(2r)}.$

Then by [19, Lemma A.5.]

$T_{2}\leq C\Delta_{l}^{1/2}.$

The end of the proof is now clear. ∎

###  B.3 Proof of Theorem 3.1

###### Proof.

We need only to bound

$T:=\mathbb{E}_{\theta}\left[\left|\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{(i)}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}-\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}\right|^{r}\right]^{1/r}.$

Now we have by Cauchy-Schwarz

$T\leq T_{1}T_{2}$

where

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}[Z_{T,\theta}^{2r}]^{1/(2r)}$
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\left[\left|\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{(i)}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}-\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}\right|^{2r}\right]^{1/(2r)}.$

By the result in Remark B.1 $T_{1}\leq C$, so we need only consider $T_{2}$.
We have by using a standard decomposition and the Minkowski inequality that

$T_{2}\leq T_{3}+T_{4}$

where

$\displaystyle T_{3}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\Bigg{[}\Bigg{|}\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}\times$
$\displaystyle\left\\{\overline{\mathbb{E}}_{\theta}[Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})|\mathcal{Y}_{t}]-\overline{\mathbb{E}}_{\theta}[Z_{T,\theta}|\mathcal{Y}_{T}]\right\\}\Bigg{|}^{2r}\Bigg{]}^{1/(2r)}$
$\displaystyle T_{4}$ $\displaystyle=$
$\displaystyle\overline{\mathbb{E}}_{\theta}\Bigg{[}\Bigg{|}\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}\Big{\\{}\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{(i)}Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]-$
$\displaystyle\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]\Big{\\}}\Bigg{|}^{2r}\Bigg{]}^{1/(2r)}.$

For $T_{3}$ one can apply Cauchy-Schwarz and [19, Lemma A.5.] we have the
upper-bound

$T_{3}\leq
C\Delta_{l}^{1/2}\overline{\mathbb{E}}_{\theta}\Bigg{[}\Bigg{|}\frac{\overline{\mathbb{E}}_{\theta}\,[\,\lambda_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})^{(i)}\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}^{l}(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\dots,\widetilde{X}_{T})\,|\,\mathcal{Y}_{T}\,]}\Bigg{|}^{4r}\Bigg{]}^{1/(4r)}.$

Now for the expectation on the R.H.S. one can use the Hölder inequality along
with Lemmata B.1-B.3 to deduce that

$T_{3}\leq C\Delta_{l}^{1/2}.$

For $T_{4}$ using Cauchy-Schwarz and Lemma B.4 we have the upper-bound

$T_{4}\leq
C\Delta_{l}^{1/2}\overline{\mathbb{E}}_{\theta}\Bigg{[}\Bigg{|}\frac{1}{\overline{\mathbb{E}}_{\theta}\,[\,Z_{T,\theta}\,|\,\mathcal{Y}_{T}\,]}\Bigg{|}^{4r}\Bigg{]}^{1/(4r)}.$

Applying Lemma B.1 we have

$T_{4}\leq C\Delta_{l}^{1/2}$

from which we conclude. ∎

## References

  * [1] Bain, A. & Crisan, D. (2009). _Fundamentals of Stochastic Filtering_. Springer: New York.
  * [2] Ballesio, M., Jasra, A., von Schwerin, E. & Tempone, R. (2020). A Wasserstein coupled particle filter for multilevel estimation. arXiv:2004.03981.
  * [3] Benveniste, A., Métivier, M. & Priouret, P. (1990). _Adaptive Algorithms and Stochastic Approximation_. New York: Springer-Verlag.
  * [4] Beskos, A., Papaspiliopoulos, O., Roberts, G., Fearnhead, P. (2006). Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes (with discussion). _J. R. Statist. Soc. Ser. B_ , 68, 333-382.
  * [5] Bierkens, J., Van Der Meulen, F., Schauer, M. (2020). Simulation of elliptic and hypo-elliptic conditional diffusions, _Advances in Applied Probability_ , 52, 173-212.
  * [6] Botha, I., Kohn, R., & Drovandi, C. (2020). Particle methods for stochastic differential equation mixed effects models. _Bayes. Anal._ (to appear).
  * [7] Campillo, F. & Le Gland, F. (1989). Maximum likelihood estimation for partially observed diffusions: Direct Maximization vs The EM algorithm. _Stoch. Proc. Appl._ , 33, 245–274.
  * [8] Cliffe, K. A., Giles, M. B., Scheichl, R., & Teckentrup, A. L. (2011). Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. _Comput. Vis. Sci._ , 14, 3–15.
  * [9] Crisan, D. (2011). Discretizing the continuous-time filtering problem: Order of Convergence. In [10], 572–597.
  * [10] Crisan, D., & Rozovskii, B. (2011)._The Oxford Handbook of Nonlinear Filtering_. Oxford University Press.
  * [11] Del Moral, P. (2004)._Feynman-Kac Formulae_. Springer.
  * [12] Del Moral, P., Doucet, A., & Singh S. S. (2010). A backward particle interpretation of Feynman-Kac formuale. _M2AN_ , 44, 947–975.
  * [13] Del Moral, P., Doucet, A., & Singh S. S. (2010). Forward smoothing using sequential Monte Carlo, arXiv:1012.5390
  * [14] Etienne, M. P., Gloaguen, P., Corff, S. L., & Olsson, J. (2020). Backward importance sampling for partially observed diffusion processes. arXiv:2002.05438.
  * [15] Gloaguen, P., Etienne, M. P. & Le Corff, S. (2018). Online sequential Monte Carlo smoother for partially observed diffusion processes. _EURASIP J. Adv. Sig. Proc_ , article 9.
  * [16] Jasra, A., Kamatani, K., Law K. J. H. & Zhou, Y. (2017). Multilevel particle filters. _SIAM J. Numer. Anal._ , 55, 3068-3096.
  * [17] Jasra, A., Law, K. J. H., & Yu, F. (2020). Unbiased filtering of a class of partially observed diffusions. arXiv:2002.03747.
  * [18] Jasra, A., & Yu, F. (2020). Central limit theorems for coupled particle filters. _Adv. Appl. Probab._ 52, 942–1001.
  * [19] Jasra, A., Yu, F. & Heng, J. (2020). Multilevel particle filters for the nonlinear filtering problem in continuous time. _Stat. Comp._ 30, 1381–1402.
  * [20] Le Gland, F. and Mevel, M. (1997). Recursive identification in hidden Markov models. _Proc. 36th IEEE Conf. Dec. Contr._ , 3468-3473.
  * [21] Mider, M., Schauer, M. & van der Meulen, F. (2020). Continuous-discrete smoothing of diffusions. arXiv:1712.03807.
  * [22] Olsson, J. & Westerborn, J. (2017). Efficient particle-based online smoothing in general hidden Markov models: The PaRIS algorithm. _Bernoulli_ , 23, 1951-1996.
  * [23] Papaspiliopoulos, O., Roberts, G. O., & Stramer, O. (2013). Data augmentation for diffusions. _J. Comp. Graph. Stat._ , 22, 665-688.
  * [24] Papaspiliopoulos, O. & Roberts, G. (2012). Importance sampling techniques for estimation of diffusion models. _Stat. Meth. Stoch. Diff. Eq._ , 124, 311-340.
  * [25] Picard, J. (1984). Approximation of nonlinear filtering problems and order of convergence. In _Filtering and control of random processes_ , 219-236, Springer, Berlin, Heidelberg.
  * [26] Poyiadjis, G., Doucet, A., & Singh, S. S. (2011). Particle approximations of the score and observed information matrix in state space models with application to parameter estimation. Biometrika, 98, 65-80.
  * [27] Särkkä, S., & Sottinen, T. (2008). Application of Girsanov theorem to particle filtering of discretely observed continuous-time non-linear systems. _Bayes. Anal._ , 3, 555-584.
  * [28] Schauer, M., van der Meulen, F. & van Zanten, H. (2017). Guided proposals for simulating multi-dimensional diffusion bridges. _Bernoulli_ , 23, 2917–2950.
  * [29] Sen, D., Thiery, A., Jasra, A. (2018). On coupling particle filters. _Statist. Comp._ , 28, 461-475.
  * [30] Surace, S. C., & Pfister, J. P. (2018). Online Maximum-Likelihood Estimation of the Parameters of Partially Observed Diffusion Processes. _IEEE Transactions on Automatic Control_ , 64(7), 2814-2829.
  * [31] Talay, D. (1984). Efficient numerical schemes for the approximation of expectations of functionals of the solution of a SDE, and applications. In _Filtering and control of random processes_ , 294-313, Springer, Berlin, Heidelberg.
  * [32] Thorisson, H. (2000). _Coupling, stationarity, and regeneration_. Springer:New York.
  * [33] van der Meulen, F., & Schauer, M. (2017). Bayesian estimation of discretely observed multi-dimensional diffusion processes using guided proposals. _Elec. J. Stat._ , 11, 2358-2396.
  * [34] Whitaker, G. A., Golightly, A., Boys, R. J., & Sherlock, C. (2017). Improved bridge constructs for stochastic differential equations. _Stat. Comp._ , 27, 885-900.
  * [35] Yonekura, S. & Beskos, A. (2020). Online smoothing for diffusion processes observed with noise. arXiv: 200312247.