###### Abstract

In this article we consider the filtering problem associated to partially
observed diffusions, with observations following a marked point process. In
the model, the data form a point process with observation times that have its
intensity driven by a diffusion, with the associated marks also depending upon
the diffusion process. We assume that one must resort to time-discretizing the
diffusion process and develop particle and multilevel particle filters to
recursively approximate the filter. In particular, we prove that our
multilevel particle filter can achieve a mean square error (MSE) of
$\mathcal{O}(\epsilon^{2})$ ($\epsilon>0$ and arbitrary) with a cost of
$\mathcal{O}(\epsilon^{-2.5})$ versus using a particle filter which has a cost
of $\mathcal{O}(\epsilon^{-3})$ to achieve the same MSE. We then show how this
methodology can be extended to give unbiased (that is with no time-
discretization error) estimators of the filter, which are proved to have
finite variance and with high-probability have finite cost. Finally, we extend
our methodology to the problem of online static-parameter estimation.

Key words: Unbiased Methods, Multilevel Monte Carlo, Non-Linear Filtering,
Point Processes, Parameter Estimation.
MSC classes: 60G55, 60G35, 62M20, 62F30

Code available at: https://github.com/maabs/Multilevel-for-Diffusions-
Observed-via-Marked-Point-Processes
Corresponding author: Miguel Alvarez. E-mail:
<EMAIL_ADDRESS>
Unbiased and Multilevel Methods for a Class of Diffusions Partially Observed
via Marked Point Processes

BY MIGUEL ALVAREZ, AJAY JASRA & HAMZA RUZAYQAT

Applied Mathematics and Computational Science Program, Computer, Electrical
and Mathematical Sciences and Engineering Division, King Abdullah University
of Science and Technology, Thuwal, 23955, KSA. E-Mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

The class of partially observed diffusion processes can be found in wide
variety of real applications; see for instance the coverage in [2, 6]. In this
paper we consider observation models which constitute marked point processes
and in particular, when the points (observation times) are a non-homogeneous
Poisson process with intensity dependent upon the diffusion; this has
applications for instance in earthquake modeling see [12, 25, 29]. A related
problem is that of trying to also, simultaneously, estimate the static (time-
homogeneous) parameters of the model and this is also considered in this
article.

The filtering of partially observed diffusion processes is notoriously
challenging as even if the transition dynamics or density of the diffusion is
available, one must resort to numerical approximation of the filter. Such
approaches, at least for low-dimensional diffusions, are often focussed upon
the particle filter methodology; see [2, 6, 12, 21] for example. In many
cases, one is unable to exactly sample from the transition of the diffusion or
its transition density is unavailable even up-to a non-negative and unbiased
estimator (see [12] for diffusions where this is possible) and so often one
has to resort to using a time-discretization of the diffusion.

In the scenario where one considers approximating the time-discretized filter
as well as estimating the static parameters, there have been several works
including [1, 4, 8, 25, 21, 23, 24] to name but a few; as just noted these
methods will have bias in terms of the time-discretization. All of these
articles depend, in some way, on the particle filter method. This approach
will generate a collection of samples in parallel which undergo transitions
using the dynamics of the discretized diffusion and then are weighted and
resampled to approximate the filter; see [2, 6]. This has been further
extended in the works of [4, 21, 23, 24] to incorporate the multilevel Monte
Carlo (MLMC) method; see [13, 14, 16] for the initial articles and [18] for a
review when using MLMC in data problems. This technique relies upon simulating
approximations of the filters at multiple different levels of time-
discretizations. The approach in [21, 23, 24] is to use a coupling of the
time-discretized dynamics and of the resampling operation to give a type of
coupled particle filter. It has been shown to reduce the computational cost to
achieve a pre-specified mean square error (MSE) versus using an ordinary
particle filter. [4] consider the problem of on-line parameter estimation for
a type of partially observed diffusion process, using particle filters.

To the best of our knowledge, none of the above papers consider the
observation model that is under-study in this paper, with the exception of
[12, 25]. The work of [12] focusses upon using the exact simulation of
diffusions methodology to produce unbiased estimators of the filter. As is
well-known, such methodology rarely works in dimensions bigger than 1 and
often places significant restrictions on the diffusions that can be used,
albeit being very elegant and useful in the scenarios where it can be adopted.
The approach of this paper really only needs the strong error of the method to
fall at a fast enough rate; often this means only that the drift and diffusion
coefficients need Lipschitz and growth conditions; see [27]. The method of
[25], which provides unbiased inference (removes time discretization bias),
relies both on the exact simulation of the diffusion, which again is a
restricted class of problem, for which there is a transition density and the
well-known Poisson estimator [35], which in some examples can have a
substantial variance. In addition, the authors often have to truncate their
estimator to yield unbiased estimators as there is a (small) probability that
the estimator of the exponential of the integrated intensity is negative. [25]
also extend their ideas in the context of Markov chain Monte Carlo and batch
(fixed data set) static Bayesian parameter estimation; the approach in this
article can be extended to performing either a multilevel parameter estimation
using Euler-Maruyama scheme as in [22], or using Runge-Kutta schemes as in [7]
or unbiased Bayesian parameter estimation for diffusions driven by a Brownian
motion [5] or diffusions driven by Lévy processes [32], with the ideas of the
afore-mentioned references, although we do not do so. Whilst both of the
methods in [12, 25] clearly have several interesting contributions, they are
not as general as the approach that we will develop.

In this article we apply and analyze methodology for the filtering problem
associated to partially observed diffusions, with observations following a
marked point process. In particular, we apply the multilevel particle filter
(MLPF) [21] and unbiased particle filter (UPF) [20] adapted to these models.
The latter method is able to approximate the filter with no time-
discretization bias, whilst only working with time-discretized dynamics. We
also apply the ideas from [4] to provide a method which can estimate the
static (time-homogeneous) parameters of the model sequentially in time (online
static parameter estimation) via maximum likelihood approach. To summarize,
the contributions of this paper are as follows:

  * •

We develop the MLPF and UPF for the model under study.

  * •

We prove that the MLPF estimator can achieve a mean square error (MSE) of
$\mathcal{O}(\epsilon^{2})$, for arbitrary $\epsilon>0$, at a cost that is
$\mathcal{O}(\epsilon^{-2.5})$. If one uses a particle filter, then the cost
is $\mathcal{O}(\epsilon^{-3})$ to achieve the same MSE.

  * •

The afore-mentioned technical results also allow us to show that our UPF
estimator is unbiased, of finite variance and with high-probability (see e.g.
[31]) has finite cost.

  * •

We adapt the approach in [4] for online static parameter estimation.

  * •

We verify our findings by implementing the methodology on several examples.

This paper is structured as follows. In Section 2 we discuss the model and our
associated algorithms for filtering, multilevel filtering and unbiased
filtering as well as for online parameter estimation. In Section 3 we present
our mathematical results. In Section 4 we provide several numerical examples
that illustrate our algorithms and mathematical analysis. The proofs for our
mathematical results can be found in the appendix.

##  2 Modeling and Algorithms

###  2.1 Notations

Let $(\mathsf{X},\mathcal{X})$ be a measurable space. For
$\varphi:\mathsf{X}\rightarrow\mathbb{R}$ we write
$\mathcal{B}_{b}(\mathsf{X})$ as the collection of bounded measurable
functions. 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 with respect to (w.r.t.) $\|\cdot\|_{2}$
($\|\cdot\|_{p}$ denotes the $\mathbb{L}_{p}-$norm of a vector
$x\in\mathbb{R}^{d}$). 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}.$

We write $\|\varphi\|_{\textrm{Lip}}$ as the Lipschitz constant of a function
$\varphi\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d})$. For
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$, we write the supremum norm
$\|\varphi\|=\sup_{x\in\mathsf{X}}|\varphi(x)|$. $\mathcal{P}(\mathsf{X})$
denotes the collection of probability measures on $(\mathsf{X},\mathcal{X})$.
For a measure $\mu$ on $(\mathsf{X},\mathcal{X})$ and a function
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$, the notation
$\mu(\varphi)=\int_{\mathsf{X}}\varphi(x)\mu(dx)$ is used. $B(\mathbb{R}^{d})$
denote the Borel sets on $\mathbb{R}^{d}$. $dx$ is used to denote the Lebesgue
measure. If $K:\mathsf{X}\times\mathcal{X}\rightarrow[0,\infty)$ is a non-
negative operator and $\mu$ is a measure, we use the notations $\mu
K(dy)=\int_{\mathsf{X}}\mu(dx)K(x,dy)$ and for
$\varphi\in\mathcal{B}_{b}(\mathsf{X})$,
$K(\varphi)(x)=\int_{\mathsf{X}}\varphi(y)K(x,dy).$ For $A\in\mathcal{X}$, the
indicator function is written as $\mathbb{I}_{A}(x)$.
$\mathcal{N}_{s}(\mu,\Sigma)$ (resp. $\psi_{s}(x;\mu,\Sigma)$) denotes an
$s-$dimensional Gaussian distribution (density evaluated at
$x\in\mathbb{R}^{s}$) of mean $\mu$ and covariance $\Sigma$. If $s=1$ we omit
the subscript $s$. For a vector/matrix $X$, $X^{*}$ is used to denote 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 (resp.
$d-$dimensional vector), $\varphi(x)$ (resp. $x$) say, we write the
$i^{\textrm{th}}-$component ($i\in\\{1,\dots,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)}$.

###  2.2 Model

We consider data $Y_{s_{1}},Y_{s_{2}},\dots$,
$Y_{s_{k}}\in\mathsf{Y}\subseteq\mathbb{R}^{d_{y}}$, $k\in\mathbb{N}$, that
are observed at the event times $s_{1},s_{2},\dots$, $s_{k}\in\mathbb{R}^{+}$
of a non-homogenous Poisson process driven by a diffusion. More precisely,
consider a diffusion process

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t}$ (1)

where $X_{0}=x_{*}\in\mathbb{R}^{d_{x}}$ given,
$b:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$,
$\sigma:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times d_{x}}$ is non-
constant (except in some of our numerical examples) and $\\{W_{t}\\}_{t\geq
0}$ is a standard Brownian motion. At this stage we do not mention any static
parameters; this is considered in Section 2.5. To minimize certain technical
difficulties, the following assumption is made throughout the paper:

  * (D1)

We have:

    1. 1.

$\sigma^{(ij)}$ is bounded with
$\sigma^{(ij)}\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$, for all
$(i,j)\in\\{1,\dots,d_{x}\\}^{2}$ and

$a(x):=\sigma(x)\sigma(x)^{*}$

is uniformly elliptic for all $x\in\mathbb{R}^{d_{x}}$.

    2. 2.

$b^{(j)}$ are bounded and
$b^{(j)}\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$, for all
$j\in\\{1,\dots,d_{x}\\}$.

Note that the assumptions here and later are far from minimal. In general, as
stated in the introduction, it should be enough that the drift and diffusion
coefficients possess Lipschitz and growth conditions.

We then consider a non-homogeneous Poisson process with intensity function
$\lambda:\mathbb{R}^{d_{x}}\rightarrow[\underline{c},\overline{c}]$,
$0<\underline{c}<\overline{c}<+\infty$. That is to say that, conditional on
the process $\\{X_{r}\\}_{r\in[0,t]}$, the joint density of
$n_{t}\in\mathbb{N}$ event times $s_{1},\dots,s_{n_{t}}$ is precisely

$\left(\prod_{k=1}^{n_{t}}\lambda(x_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda(x_{r})dr\right\\}.$

The observations, conditional upon the event times $s_{1:n_{t}}$ and the
diffusion process $\\{X_{r}\\}_{r\in[0,t]}$, have joint density

$p(y_{s_{1}:s_{n_{t}}}|s_{1:n_{t}},\\{x_{r}\\}_{r\in[0,t]})=\prod_{k=1}^{n_{t}}g(x_{s_{k}},y_{s_{k}})$

where, for any $x\in\mathbb{R}^{d_{x}}$, $g(x,\cdot)$ is a probability density
on $\mathsf{Y}$.

Our objective is to do filtering at some finite collection of times. For now,
consider $t\in\mathbb{R}^{+}$ given, then, for
$\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$, we want to compute

$\pi_{t}(\varphi):=\frac{\mathbb{E}\left[\varphi(X_{t})\left(\prod_{k=1}^{n_{t}}g(X_{s_{k}},y_{s_{k}})\lambda(X_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda(X_{r})dr\right\\}\right]}{\mathbb{E}\left[\left(\prod_{k=1}^{n_{t}}g(X_{s_{k}},y_{s_{k}})\lambda(X_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda(X_{r})dr\right\\}\right]},$
(2)

where the expectation is w.r.t. the law of the diffusion process given in (1).
We will make the constraint that this computation is of interest at times
$t\in\mathbb{N}$.

###  2.3 Discretization

In practice, it is typically difficult (or impossible) to work directly with
(1) and (2). Often one works with a time discretization of (1), of which we
adopt the Euler-Maruyama approximation, with $\Delta_{l}=2^{-l}$,
$l\in\mathbb{N}_{0}$:

$\widetilde{X}_{(k+1)\Delta_{l}}=\widetilde{X}_{k\Delta_{l}}+b(\widetilde{X}_{k\Delta_{l}})+\sigma(\widetilde{X}_{k\Delta_{l}})[W_{(k+1)\Delta_{l}}-W_{k\Delta_{l}}]$
(3)

where $k\in\mathbb{N}_{0}$ and $\widetilde{X}_{0}=x_{*}$. We make the standard
extension that if $t\in[k\Delta_{l},(k+1)\Delta_{l})$ then

$\widetilde{X}_{t}=\widetilde{X}_{k\Delta_{l}}+(\widetilde{X}_{(k+1)\Delta_{l}}-\widetilde{X}_{k\Delta_{l}})\Delta_{l}^{-1}(t-k\Delta_{l}).$
(4)

This is clearly needed given the representation (2) that we seek to
approximate.

Given the Euler approximation (3), we can then consider the approximation, for
$t\in\mathbb{N}$, given by

$\pi_{t}^{l}(\varphi):=\frac{\mathbb{E}\left[\varphi(X_{t})\left(\prod_{k=1}^{n_{t}}g(\widetilde{X}_{s_{k}},y_{s_{k}})\lambda(\widetilde{X}_{s_{k}})\right)\exp\left\\{-\Delta_{l}\sum_{k=0}^{t\Delta_{l}^{-1}-1}\lambda(\widetilde{X}_{k\Delta_{l}})\right\\}\right]}{\mathbb{E}\left[\left(\prod_{k=1}^{n_{t}}g(\widetilde{X}_{s_{k}},y_{s_{k}})\lambda(\widetilde{X}_{s_{k}})\right)\exp\left\\{-\Delta_{l}\sum_{k=0}^{t\Delta_{l}^{-1}-1}\lambda(\widetilde{X}_{k\Delta_{l}})\right\\}\right]},$
(5)

where the expectation is taken w.r.t. the law associated to the Euler
approximation (3).

We introduce an additional assumption

  * (D2)

We have

    1. 1.

$\lambda\in\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$.

    2. 2.

For any $y\in\mathsf{Y}$,
$g(\cdot,y)\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})\cap\textrm{Lip}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$.

We have the following result.

######  Proposition 2.1.

Assume (D(D1)-(D2)). For
$(t,\varphi)\in\mathbb{N}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})\cap\textrm{\emph{Lip}}_{\|\cdot\|_{2}}(\mathbb{R}^{d_{x}})$
there exists a $C<+\infty$ such that for any $l\in\mathbb{N}_{0}$

$|\pi_{t}^{l}(\varphi)-\pi_{t}(\varphi)|\leq C\Delta_{l}.$

###### Proof.

One can see [2, 27] for justifications of this result as it is standard in the
literature. ∎

###  2.4 Filtering

We shall now follow [19] and present the filter in a recursive form. We will
then detail some algorithms for approximating the afore-mentioned filter.

####  2.4.1 The Filter Measure

Throughout the section $l\in\mathbb{N}_{0}$ is given and as $x_{*}$ is fixed,
it is removed from the notation where possible. This presentation closely
follows [19] and the purpose of these notations is to achieve a consistency
with the literature on particle filters. Let $p\in\mathbb{N}_{0}$ be given, we
use the notation

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

which denotes a path on a unit time $(p,p+1]$; we have removed the
$\widetilde{\cdot}$ notation for simplicity. We now set for
$p\in\mathbb{N}_{0}$

$G_{p}^{l}(u_{p-1}^{l},u_{p}^{l}):=\left(\prod_{k=n_{p}+1}^{n_{p+1}}g(x_{s_{k}},y_{s_{k}})\lambda(x_{s_{k}})\right)\exp\left\\{-\Delta_{l}\sum_{k=0}^{\Delta_{l}^{-1}-1}\lambda(x_{p+(k+1)\Delta_{l}})\right\\},$

where $n_{0}=0$, $u_{-1}^{l}=x_{*}$ and this is a function on $E_{l-1}\times
E_{l}$ (in fact it is a function on $\mathbb{R}^{d_{x}}\times E_{l}$ because
only $x_{p}^{l}$ from $u_{p-1}^{l}$ is used as we will see later) due to the
interpolation defined in (4). The density function of Euler-Maruyama
discretizations can be characterized by Gaussian densities; we denote by
$m^{l}$ the transition of (3). For $p\in\mathbb{N}$, the initial measures and
Markov kernels that we will need are

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

It is worth noting that the dependence of both
$G_{p}^{l}(u_{p-1}^{l},u_{p}^{l})$ and $M^{l}(u_{p-1}^{l},du_{p}^{l})$ on
$u_{p-1}^{l}$ is present only on its last element $x_{p}^{l}$, meaning that
for $M^{l}$ the transition is only dependent on $x_{p}^{l}$, we will use these
facts below.

The time-discretized filter _on path space_ (i.e. over time interval
$(p-1,p]$, $p\in\mathbb{N}$) can then be written, using the above notation,
for any $p\in\\{2,3,\dots\\}$ as

$\overline{\pi}_{p}^{l}(du_{p-1}^{l}):=\frac{\int_{E_{l}^{p-1}}\left\\{\prod_{k=0}^{p-1}G_{k}^{l}(u_{k-1}^{l},u_{k}^{l})\right\\}\eta_{0}^{l}(du_{0}^{l})\prod_{k=1}^{p-1}M^{l}(u_{k-1}^{l},du_{k}^{l})}{\int_{E_{l}^{p}}\left\\{\prod_{k=0}^{p-1}G_{k}^{l}(u_{k-1}^{l},u_{k}^{l})\right\\}\eta_{0}^{l}(du_{0}^{l})\prod_{k=1}^{p-1}M^{l}(u_{k-1}^{l},du_{k}^{l})}$
(6)

with

$\overline{\pi}_{1}^{l}(du_{0}^{l}):=\frac{G_{0}^{l}(u_{-1}^{l},u_{0}^{l})\leavevmode\nobreak\
\eta_{0}^{l}(du_{0}^{l})}{\int_{E_{l}}G_{0}^{l}(u_{-1}^{l},u_{0}^{l})\leavevmode\nobreak\
\eta_{0}^{l}(du_{0}^{l})}.$

This collection of probability measures will be of use later on in the
article. We note that, for any $p\in\mathbb{N}$, $l\in\mathbb{N}_{0}$, and
$\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$, we have that the expectation
of $\varphi$ w.r.t. the approximate filter distribution at time $p$ as

$\pi_{p}^{l}(\varphi)=\int_{E_{l}}\varphi(x_{p})\overline{\pi}_{p}^{l}(du_{p-1}^{l}).$

####  2.4.2 Particle and Coupled Particle Filter

We begin by describing the particle filter for approximating
$\pi_{p}^{l}(\varphi)$ for a given $l\in\mathbb{N}_{0}$, with $p\in\mathbb{N}$
and any $\varphi:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}$ that is
$\pi_{p}^{l}-$integrable. For a given $N\in\mathbb{N}$, the particle filter
generates a system of random variables on $(E_{l}^{N})^{n+1}$ at a time
$n\in\mathbb{N}_{0}$ according to the probability measure

$\mathbb{Q}(d(u_{0}^{l,1:N},\dots,u_{n}^{l,1:N}))=\left\\{\prod_{i=1}^{N}M^{l}(x_{*},du_{0}^{l,i})\right\\}\prod_{p=0}^{n-1}\prod_{i=1}^{N}\left\\{\sum_{j=1}^{N}\frac{G_{p}^{l}(u_{p-1}^{l,j},u_{p}^{l,j})}{\sum_{k=1}^{N}G_{p}^{l}(u_{p-1}^{l,k},u_{p}^{l,k})}M^{l}(u_{p}^{l,j},du_{p+1}^{l,i})\right\\}.$

An algorithmic description of the particle filter is given in Algorithm 1. For
$(t,\varphi)\in\mathbb{N}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$ one can
approximate the time-discretized filter $\pi_{t}^{l}(\varphi)$ corresponding
to step-size $\Delta_{l}$ via

$\pi_{t}^{l,N}(\varphi):=\frac{\sum_{i=1}^{N}G_{t-1}^{l}(u_{t-2}^{l,i},u_{t-1}^{l,i})\varphi(x_{t}^{l,i})}{\sum_{i=1}^{N}G_{t-1}^{l}(u_{t-2}^{l,i},u_{t-1}^{l,i})},$
(7)

which can easily be shown to converge (e.g. in probability as
$N\rightarrow\infty$) to $\pi_{t}^{l}(\varphi)$; see [6]. We define an
empirical measure at time $p-1$ and level $l$ which will be needed later on
through the following expectation. Let $\varphi\in\mathcal{B}_{b}(E_{l}^{2})$
and $p\in\mathbb{N}$, we define

$\eta_{p-1}^{l,N}(\varphi):=\frac{1}{N}\sum_{i=1}^{N}\varphi(u_{p-2}^{l,i},u_{p-1}^{l,i}).$
(8)

  1. 1.

Initialize: For $i\in\\{1,\dots,N\\}$, generate $u_{0}^{l,i}$ from
$M^{l}(x_{*},\cdot)$. Set $p=0$.

  2. 2.

Update: For $i\in\\{1,\dots,N\\}$, generate $u_{p+1}^{l,i}$ from

$\sum_{j=1}^{N}\frac{G_{p}^{l}(u_{p-1}^{l,j},u_{p}^{l,j})}{\sum_{k=1}^{N}G_{p}^{l}(u_{p-1}^{l,k},u_{p}^{l,k})}M^{l}(u_{p}^{l,j},du_{p+1}^{l,i}).$

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

Algorithm 1 Particle Filter.

We now consider the coupled particle filter (CPF) as developed in [21, 23]
(see also [24]). Contrary to the presentation in [19] we restrict ourselves to
an algorithmic rather than operator based description. The objective of the
coupled particle filter, for $l\in\mathbb{N}$ given and up-to a time
$t\in\mathbb{N}$, is to generate two clouds of particles $u_{t}^{l,1:N}$ and
$\bar{u}_{t}^{l-1,1:N}$ so that $u_{t}^{l,1:N}$ (resp.
$\bar{u}_{t}^{l-1,1:N}$) can be used to approximate $\pi_{t}^{l}(\varphi)$
(resp. $\pi_{t}^{l-1}(\varphi)$). Moreover, that there is a dependence between
these two clouds of particles which ensure that the variance of terms which
approximate the difference $\pi_{t}^{l}(\varphi)-\pi_{t}^{l-1}(\varphi)$ will
fall with $l$.

To describe the CPF in its most simple form, we need two algorithms: coupled
sampling and coupled resampling. We begin with the former, which is given in
Algorithm 2 and simply is the well-known synchronous coupling of Euler-
discretized diffusions. The method of coupled resampling, developed for
multilevel applications in [21] can be found in Algorithm 3. This resampling
algorithm maximizes the probability that two sampled indices are equal. As
noted in [17], Algorithm 3 is by no means optimal, but appears to be the most
used method in the literature.

We now give the CPF algorithm in Algorithm 4. This is simply a type of
particle filter that utilizes Algorithm 2 for sampling and Algorithm 3 for
resampling. The algorithm can estimate the difference
$\pi_{t}^{l}(\varphi)-\pi_{t}^{l-1}(\varphi)$, $t\in\mathbb{N}$, using the
expression

$[\pi_{t}^{l}-\pi_{t}^{l-1}]^{N}(\varphi):=\frac{\sum_{i=1}^{N}G_{t-1}^{l}(u_{t-2}^{l,i},u_{t-1}^{l,i})\varphi(x_{t}^{l,i})}{\sum_{i=1}^{N}G_{t-1}^{l}(u_{t-2}^{l,i},u_{t-1}^{l,i})}-\frac{\sum_{i=1}^{N}G_{t-1}^{l-1}(\bar{u}_{t-2}^{l-1,i},\bar{u}_{t-1}^{l-1,i})\varphi(\bar{x}_{t}^{l-1,i})}{\sum_{i=1}^{N}G_{t-1}^{l-1}(\bar{u}_{t-2}^{l-1,i},\bar{u}_{t-1}^{l-1,i})}.$
(9)

We define two empirical measures at time $p-1$ and levels $l$, $l-1$ through
the following expectations. For
$(\varphi_{1},\varphi_{2})\in\mathcal{B}_{b}(E_{l}^{2})\times\mathcal{B}_{b}(E_{l-1}^{2})$,
$p\in\mathbb{N}$, we have

$\check{\eta}_{p-1}^{l,N}(\varphi_{1}):=\frac{1}{N}\sum_{i=1}^{N}\varphi_{1}(u_{t-2}^{l,i},u_{t-1}^{l,i})\quad\textrm{and}\quad\check{\bar{\eta}}_{p-1}^{l-1,N}(\varphi_{2}):=\frac{1}{N}\sum_{i=1}^{N}\varphi_{2}(\bar{u}_{t-2}^{l-1,i},\bar{u}_{t-1}^{l-1,i}).$
(10)

  1. 1.

Input: level $l$ and starting points $(x_{0}^{l},x_{0}^{l-1})$.

  2. 2.

Generate
$Z_{k}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d}(0,\Delta_{l}I_{d})$,
$k\in\\{1,2,\dots,\Delta_{l}^{-1}\\}$.

  3. 3.

Level $l$: For $k\in\\{0,1,2,\dots,\Delta_{l}^{-1}-1\\}$ with
$X_{0}^{l}=x_{0}^{l}$ generate

$X_{(k+1)\Delta_{l}}^{l}=X_{k\Delta_{l}}^{l}+b(X_{k\Delta_{l}}^{l})\Delta_{l}+\sigma(X_{k\Delta_{l}}^{l})Z_{k+1}.$

Set $U_{0}^{l}=(X_{\Delta_{l}}^{l},X_{2\Delta_{l}}^{l},\dots,X_{1}^{l})$.

  4. 4.

Level $l-1$: For $k\in\\{0,1,2,\dots,\Delta_{l-1}^{-1}-1\\}$ with
$X_{0}^{l-1}=x_{0}^{l-1}$ generate

$\displaystyle X_{(k+1)\Delta_{l-1}}^{l-1}$ $\displaystyle=$ $\displaystyle
X_{k\Delta_{l-1}}^{l-1}+b(X_{k\Delta_{l-1}}^{l-1})\Delta_{l-1}+\sigma(X_{k\Delta_{l-1}}^{l-1})\\{Z_{2(k+1)-1}+Z_{2(k+1)}\\}.$

Set
$\overline{U}_{0}^{l-1}=(X_{\Delta_{l-1}}^{l-1},X_{2\Delta_{l-1}}^{l-1},\dots,X_{1}^{l-1})$.

  5. 5.

Output: $(U_{0}^{l},\overline{U}_{0}^{l-1})$.

Algorithm 2 Coupled Euler Scheme on $[0,1]$.

  1. 1.

Input: $N\in\mathbb{N}$, two clouds of particles $(U_{1}^{1:N},U_{2}^{1:N})$
and their associated probabilities $(W_{1}^{1:N},W_{2}^{1:N})$.

  2. 2.

For $i\in\\{1,\dots,N\\}$, generate $R^{i}\sim\mathcal{U}_{[0,1]}$

     * •

If $R^{i}<\sum_{i=1}^{N}\min\\{W_{1}^{i},W_{2}^{i}\\}$, generate
$a^{i}\in\\{1,\dots,N\\}$ using the probability mass function

$\mathbb{P}(i)=\frac{\min\\{W_{1}^{i},W_{2}^{i}\\}}{\sum_{j=1}^{N}\min\\{W_{1}^{j},W_{2}^{j}\\}}$

and set $\tilde{U}_{j}^{i}=U_{j}^{a^{i}}$, $j\in\\{1,2\\}$.

     * •

Otherwise generate $(a_{1}^{i},a_{2}^{i})\in\\{1,\dots,N\\}^{2}$ using any
coupling of the probability mass functions:

$\mathbb{P}_{j}(i)=\frac{W_{j}^{i}-\min\\{W_{1}^{i},W_{2}^{i}\\}}{\sum_{k=1}^{N}[W_{j}^{k}-\min\\{W_{1}^{k},W_{2}^{k}\\}]}$

and set $\tilde{U}_{j}^{i}=U_{j}^{a_{j}^{i}}$, $j\in\\{1,2\\}$.

  3. 3.

Set: $U_{j}^{i}=\tilde{U}_{j}^{i}$, $(i,j)\in\\{1,\dots,N\\}\times\\{1,2\\}$.

  4. 4.

Output: $(U_{1}^{1:N},U_{2}^{1:N})$.

Algorithm 3 Maximal Coupling Resampling

  1. 1.

Input: $(l,N)\in\mathbb{N}^{2}$.

  2. 2.

Initialize: For $i\in\\{1,\dots,N\\}$, generate
$u_{0}^{l,i},\bar{u}_{0}^{l-1,i}$ by using Algorithm 2 with level $l$ and
starting points $(x_{*},x_{*})$. Set $p=1$.

  3. 3.

Iterate: For $i\in\\{1,\dots,N\\}$, compute the weights

$w_{p-1}^{l,i}=\frac{G_{p-1}^{l}(u_{p-2}^{l,i},u_{p-1}^{l,i})}{\sum_{j=1}^{N}G_{p-1}^{l}(u_{p-2}^{l,j},u_{p-1}^{l,j})}\quad\textrm{and}\quad\bar{w}_{p-1}^{l-1,i}=\frac{G_{p-1}^{l-1}(\bar{u}_{p-2}^{l-1,i},\bar{u}_{p-1}^{l-1,i})}{\sum_{j=1}^{N}G_{p-1}^{l-1}(\bar{u}_{p-2}^{l-1,j},\bar{u}_{p-1}^{l-1,j})}.$

Perform Algorithm 3 with inputs $N$,
$(u_{p-1}^{l,1:N},\bar{u}_{p-1}^{l-1,1:N})$ and
$(w_{p-1}^{l,1:N},\bar{w}_{p-1}^{l-1,1:N})$. For $i\in\\{1,\dots,N\\}$,
generate $u_{p}^{l,i},\bar{u}_{p}^{l-1,i}$ by using Algorithm 2 with level $l$
and starting points $(x_{p}^{l,i},\bar{x}_{p}^{l-1,i})$. Set $p=p+1$ and
return to the start of 3.

Algorithm 4 Coupled Particle Filter

####  2.4.3 Multilevel Particle Filter

We can now describe the MLPF using the PF and the CPF.

  1. 1.

Level 0: Run a PF as in Algorithm 1 with $N_{0}$ samples, independently of all
other levels.

  2. 2.

For each level $l\in\\{1,\dots,L\\}$: Run a CPF (to approximate the time-
discretized filters at levels $l$ and $l-1$) as in Algorithm 4 with $N_{l}$
samples, independently of all other levels.

An estimator of expectations $\pi_{t}^{L}(\varphi)$ for
$t\in\mathbb{N},\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$ with respect to
the time-discretized filter at the highest level $L$ is then given by

$\pi_{t}^{L,ML}(\varphi):=\pi_{t}^{0,N_{0}}(\varphi)+\sum_{l=1}^{L}[\pi_{t}^{l}-\pi_{t}^{l-1}]^{N_{l}}(\varphi),$
(11)

where $\pi_{t}^{0,N_{0}}(\varphi)$ is the estimator in (7) at level $l=0$ with
$N=N_{0}$ samples, and $[\pi_{t}^{l}-\pi_{t}^{l-1}]^{N_{l}}(\varphi)$ is the
estimator in (9) with $N=N_{l}$ samples.

####  2.4.4 Unbiased Particle Filter

The method of [20] can be used to remove the time-discretization bias and we
detail the method with as much brevity as possible; details can be found in
the afore-mentioned reference.

The basic method in [20] uses a double randomization technique which is based
upon the approaches in [30, 31] (see also [34]). Let
$\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$ be fixed but arbitrary and
note that the idea below can be used with little effort for many different
$\varphi$. This constitutes selecting a positive probability mass function
$\mathbb{P}_{L}$ on $\mathbb{N}_{0}$ and constructing a sequence of
independent random variables $(\Xi_{l})_{l\in\mathbb{N}_{0}}$ so that for any
$(l,t)\in\mathbb{N}^{2}$:

$\mathbb{E}[\Xi_{l}]=\pi_{t}^{l}(\varphi)-\pi_{t}^{l-1}(\varphi)$

with $\mathbb{E}[\Xi_{0}]=\pi_{t}^{0}(\varphi)$ and $\mathbb{E}$ is the
expectation w.r.t. the law associated to the simulation of the algorithm. Then
one can sample $L$ from $\mathbb{P}_{L}$ and then construct $\Xi_{L}$ to
obtain the estimator

$\widehat{\pi_{t}(\varphi)}=\frac{\Xi_{L}}{\mathbb{P}_{L}(L)},$ (12)

which is an unbiased estimator of $\pi_{t}(\varphi)$ and moreover, if one
selects $\mathbb{P}_{L}$ so that

$\sum_{l\in\mathbb{N}_{0}}\frac{\mathbb{E}[\Xi_{l}^{2}]}{\mathbb{P}_{L}(l)}<+\infty,$

then the variance of the estimator in (12) is finite. [20] provides a method
for constructing the sequence $(\Xi_{l})_{l\in\mathbb{N}_{0}}$ and we now
describe this.

We begin with the computation of $\Xi_{0}$ and the procedure that is needed is
detailed in Algorithm 5. The method requires a positive conditional
probability mass function $\mathbb{P}_{p}(\cdot|l)$ 111An independent
probability mass function is also possible; see Section 4.3. on
$\mathbb{N}_{0}$ and a sequence $(N_{p})_{p\geq 0}$ of non-decreasing positive
integers such that $N_{p}\uparrow\infty$ (e.g.
$N_{p}=N_{0}\leavevmode\nobreak\ 2^{p}$). To compute $\Xi_{0}$ we set, for
$(p,t,\varphi)\in\mathbb{N}_{0}\times\mathbb{N}\times\mathcal{B}_{b}(E_{0}^{2})$
given,

$\widetilde{\eta}_{t-1}^{0,N_{0:p}}(\varphi):=\sum_{q=0}^{p}\left(\frac{N_{q}-N_{q-1}}{N_{p}}\right)\eta_{t-1}^{0,N_{q}-N_{q-1}}(\varphi)$

with $N_{-1}=0$. Then set for $\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\widetilde{\pi}_{t}^{0,N_{0:p}}(\varphi):=\frac{\int_{E_{0}^{2}}\varphi(x_{t}^{0})\leavevmode\nobreak\
G_{t-1}^{0}(u_{t-2}^{0},u_{t-1}^{0})\leavevmode\nobreak\
\widetilde{\eta}_{t-1}^{0,N_{0:p}}\left(d(u_{t-2}^{0},u_{t-1}^{0})\right)}{\int_{E_{0}^{2}}G_{t-1}^{0}(u_{t-2}^{0},u_{t-1}^{0})\leavevmode\nobreak\
\widetilde{\eta}_{t-1}^{0,N_{0:p}}\left(d(u_{t-2}^{0},u_{t-1}^{0})\right)}.$

Then we have

$\Xi_{0}=\frac{1}{\mathbb{P}(p|l=0)}\left(\widetilde{\pi}_{t}^{0,N_{0:p}}(\varphi)-\widetilde{\pi}_{t}^{0,N_{0:p-1}}(\varphi)\right)$
(13)

with the convention that $\widetilde{\pi}_{t}^{0,N_{0:-1}}(\varphi)=0$.

  1. 1.

Sample $P$ from $\mathbb{P}_{p}(\cdot|l)$.

  2. 2.

Run Algorithm 1 with $N_{0}$ samples until time $t-1$ and denote the empirical
measure (as in (8)) as $\eta_{t-1}^{l,N_{0}}(\cdot)$. Set $q=1$. If $p=0$
stop; otherwise go to the next step.

  3. 3.

Independently of all other random variables, run Algorithm 1 with
$N_{q}-N_{q-1}$ samples until time $t-1$ and denote the empirical measure (as
in (8)) as $\eta_{t-1}^{l,N_{q}-N_{q-1}}(\cdot)$. Set $q=q+1$. If $q=p+1$
stop; otherwise go to the start of 3.

Algorithm 5 Computing $\Xi_{0}$. Throughout
$(t,l)\in\mathbb{N}\times\mathbb{N}_{0}$ are given.

For the computation of $(\Xi_{l})_{l\in\mathbb{N}}$, the procedure is similar,
except one uses a CPF instead of a PF; see Algorithm 6. We set, for
$(p,t,\varphi_{1})\in\mathbb{N}_{0}\times\mathbb{N}\times\mathcal{B}_{b}(E_{l}^{2})$
given,

$\widetilde{\eta}_{t-1}^{l,N_{0:p}}(\varphi_{1}):=\sum_{q=0}^{p}\left(\frac{N_{q}-N_{q-1}}{N_{p}}\right)\check{\eta}_{t-1}^{l,N_{q}-N_{q-1}}(\varphi_{1})$

and, for
$(p,t,\varphi_{2})\in\mathbb{N}_{0}\times\mathbb{N}\times\mathcal{B}_{b}(E_{l-1}^{2})$
given,

$\widetilde{\bar{\eta}}_{t-1}^{l-1,N_{0:p}}(\varphi_{2}):=\sum_{q=0}^{p}\left(\frac{N_{q}-N_{q-1}}{N_{p}}\right)\check{\bar{\eta}}_{t-1}^{l-1,N_{q}-N_{q-1}}(\varphi_{2}).$

Then set for $\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\displaystyle[\widetilde{\pi}_{t}^{l}$
$\displaystyle-\widetilde{\pi}_{t}^{l-1}]^{N_{0:p}}(\varphi)=$
$\displaystyle\frac{\int_{E_{l}^{2}}\varphi(x_{t}^{l})G_{t-1}^{l}(u_{t-2}^{l},u_{t-1}^{l})\widetilde{\eta}_{t-1}^{l,N_{0:p}}\left(d(u_{t-2}^{l},u_{t-1}^{l})\right)}{\int_{E_{l}^{2}}\varphi(x_{t}^{l})G_{t-1}^{l}(u_{t-2}^{l},u_{t-1}^{l})\widetilde{\eta}_{t-1}^{l,N_{0:p}}\left(d(u_{t-2}^{l},u_{t-1}^{l})\right)}-\frac{\int_{E_{l-1}^{2}}\varphi(x_{t}^{l-1})G_{t-1}^{l-1}(\bar{u}_{t-2}^{l-1},\bar{u}_{t-1}^{l-1})\widetilde{\bar{\eta}}_{t-1}^{l-1,N_{0:p}}\left(d(\bar{u}_{t-2}^{l-1},\bar{u}_{t-1}^{l-1})\right)}{\int_{E_{l-1}^{2}}G_{t-1}^{l-1}(\bar{u}_{t-2}^{l-1},\bar{u}_{t-1}^{l-1})\widetilde{\bar{\eta}}_{t-1}^{l-1,N_{0:p}}\left(d(\bar{u}_{t-2}^{l-1},\bar{u}_{t-1}^{l-1})\right)}.$

Then we have

$\Xi_{l}=\frac{1}{\mathbb{P}(p|l)}\Big{(}[\widetilde{\pi}_{t}^{l}-\widetilde{\pi}_{t}^{l-1}]^{N_{0:p}}(\varphi)-[\widetilde{\pi}_{t}^{l}-\widetilde{\pi}_{t}^{l-1}]^{N_{0:p-1}}(\varphi)\Big{)}$
(14)

with the convention that
$[\widetilde{\pi}_{t}^{l}-\widetilde{\pi}_{t}^{l-1}]^{N_{0:-1}}(\varphi)=0$.

  1. 1.

Sample $P$ from $\mathbb{P}_{p}(\cdot|l)$.

  2. 2.

Run Algorithm 4 with $N_{0}$ samples until time $t-1$ and denote the empirical
measures (as in (10)) as $\check{\eta}_{t}^{l,N_{0}}(\cdot)$ and
$\check{\bar{\eta}}_{t-1}^{l-1,N_{0}}(\cdot)$. Set $q=1$. If $p=0$ stop;
otherwise go to the next step.

  3. 3.

Independently of all other random variables, run Algorithm 4 with
$N_{q}-N_{q-1}$ samples and denote the empirical measures (as in (10)) as
$\check{\eta}_{t-1}^{l,N_{q}-N_{q-1}}(\cdot)$ and
$\check{\bar{\eta}}_{t-1}^{l-1,N_{q}-N_{q-1}}(\cdot)$. Set $q=q+1$. If $q=p+1$
stop; otherwise go to the start of 3.

Algorithm 6 Computing $\Xi_{l}$. Throughout
$(t,l)\in\mathbb{N}\times\mathbb{N}$ are given.

Then to compute the unbiased estimator, one can simply use the expression in
(12), repeated $M\in\mathbb{N}$ times and in parallel; i.e. a Monte Carlo
estimator

$\displaystyle\widehat{\widehat{\pi_{t}(\varphi)}}:=\frac{1}{M}\sum_{i=1}^{M}\frac{\Xi_{L_{i}}}{\mathbb{P}_{L}(L_{i})},$
(15)

see [20, Algorithm 5] for further details. Note that the finite variance and
unbiasedness is to be established in Section 3.

###  2.5 Parameter Estimation

####  2.5.1 Model with Static Parameters and Score Function

We consider the model as in Section 2.2 except with some static parameters
$\theta\in\Theta$. For instance, the diffusion process $\\{X_{t}\\}_{t\geq 0}$
depends on a static parameter $\theta\in\Theta$ in the following manner

$dX_{t}=b_{\theta}(X_{t})dt+\sigma(X_{t})dW_{t}$ (16)

such that $X_{0}=x_{*}\in\mathbb{R}^{d_{x}}$ is given,
$b_{\theta}:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$
$\forall\theta\in\Theta$,
$\sigma:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times d_{x}}$ is non-
constant, and $\\{W_{t}\\}_{t\geq 0}$ is a standard Brownian motion. In the
remainder of the model (i.e. as in Section 2.2) we allow both $g$ and
$\lambda$ to depend on $\theta$ as well, using the subscript $\theta$ in the
notation from herein.

We define, for $\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$,
$t\in\mathbb{N}$

$\displaystyle\gamma_{t,\theta}(\varphi):=\mathbb{E}_{\theta}\left[\varphi(X_{t})\left(\prod_{k=1}^{n_{t}}g_{\theta}(X_{s_{k}},y_{s_{k}})\lambda_{\theta}(X_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda_{\theta}(X_{r})dr\right\\}\right]$

where the expectation $\mathbb{E}_{\theta}$ is w.r.t. law of the process (16).
Notice that for $\varphi(x)\equiv 1$, this represents the likelihood of the
data. Our objective is to perform an online (as the data arrive) static
parameter estimation using the score function which (under simple assumptions
which are not stated here), for $t\in\mathbb{N}$, is given by

$\displaystyle\nabla_{\theta}\log\left(\gamma_{t,\theta}(1)\right)=\frac{\mathbb{E}_{\theta}\left[\zeta_{t,\theta}\left(\prod_{k=1}^{n_{t}}g_{\theta}(X_{s_{k}},y_{s_{k}})\lambda_{\theta}(X_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda_{\theta}(X_{r})dr\right\\}\right]}{\mathbb{E}_{\theta}\left[\left(\prod_{k=1}^{n_{t}}g_{\theta}(X_{s_{k}},y_{s_{k}})\lambda_{\theta}(X_{s_{k}})\right)\exp\left\\{-\int_{0}^{t}\lambda_{\theta}(X_{r})dr\right\\}\right]},$
(17)

where

$\displaystyle\zeta_{t,\theta}:=\int_{0}^{t}\left(\nabla_{\theta}b_{\theta}\left(X_{r}\right)\right)^{*}a\left(X_{r}\right)^{-1}\sigma\left(X_{r}\right)dW_{r}+\sum_{k=1}^{n_{t}}\nabla_{\theta}\log\left(g_{\theta}(X_{s_{k}},y_{s_{k}})\lambda_{\theta}(X_{s_{k}})\right)-\int_{0}^{t}\nabla_{\theta}\lambda_{\theta}(X_{r})dr.$
(18)

An analogous derivation of (17) can be found in Appendix A of [4].

####  2.5.2 Stochastic Gradient Approach

In our parameter estimation strategy, we attempt to maximize the limiting
average log-likelihood

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

which, as well as its gradient $\nabla\mathcal{L}(\theta)$, can be shown to be
an ergodic average under appropriate stability and regularity conditions [33].
One can estimate the parameters using the gradient
$\nabla_{\theta}\log\left(\gamma_{t,\theta}(1)\right)$ and a stochastic
gradient ascent methodology as

$\displaystyle\theta_{m+1}=\theta_{m}+\alpha_{m}\nabla_{\theta}\log\left(\gamma_{m+1,\theta_{0:m}}(1)\right),$
(19)

where $m\in\mathbb{N}_{0}$, $\theta_{0}\in\mathbb{R}^{d_{\theta}}$ is given,
and $\\{\alpha_{m}\\}_{m\in\mathbb{N}}$ is a decreasing sequence of step-
sizes. The recursion (19) is not suitable for online computations since the
complexity of the score function is $\mathcal{O}(m)$. In practice, we use
online updates of the score function
$\nabla_{\theta}\log\left(\gamma_{m+1,\theta_{0:m}}(1)\right)$ that include
all updates of the parameters up to time $m-1$, where the filter in the
interval of time $(m,m+1]$ is updated using the parameter $\theta_{m}$, see
[28, 26] (more details in Section 4.4). This approach can be counterproductive
in the sense that the first parameters $\theta_{0}$,
$\theta_{1},\cdots,\theta_{q}$, $q<<m$ are always included in the score
function. In order to alleviate this problem we use a different estimator
$\nabla_{\theta}\log\left(\gamma_{c(m+1),\theta_{0:m}}(1)\right)-\nabla_{\theta}\log\left(\gamma_{cm,\theta_{0:m-1}}(1)\right)$,
$c\in\mathbb{N}$, that highlights the new observations. Its stochastic
gradient ascent equation is

$\displaystyle\theta_{m+1}=\theta_{m}+\alpha_{m}\left[\nabla_{\theta}\log\left(\gamma_{c(m+1),\theta_{0:m}}(1)\right)-\nabla_{\theta}\log\left(\gamma_{cm,\theta_{0:m-1}}(1)\right)\right]$
(20)

with the convention
$\nabla_{\theta}\log\left(\gamma_{0,\theta_{0:-1}}(1)\right)=0$. It can be
shown, for $\\{\alpha_{m}\\}_{m\in\mathbb{N}_{0}}$ such that $\alpha_{m}\to
0$, $\sum_{m\in\mathbb{N}_{0}}\alpha_{m}=\infty$ and
$\sum_{m\in\mathbb{N}_{0}}\alpha_{m}^{2}<\infty$, that $\theta_{m}$ in (19)
and (20) converges as $m\rightarrow\infty$, for proofs see [3, 28]. In the
following sections, we study the particle systems of discretized estimation of
the score function.

####  2.5.3 Discretized Score

In order to compute practical approximations of the score, we resort to the
same discretization scheme displayed in (3) and (5) where we replace the
functions $b,$ $g$ and $\lambda$ with its $\theta$ dependent counterparts. The
time discretization of $\zeta_{t,\theta}$ is written as

$\displaystyle\zeta_{t,\theta}^{l}\left(x_{0},x_{\Delta_{l}},\ldots,x_{t}\right):=$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}-1}\left\\{\left(\nabla_{\theta}b_{\theta}\left(x_{k\Delta_{l}}\right)\right)^{*}a\left(x_{k\Delta_{l}}\right)^{-1}\sigma\left(x_{k\Delta_{l}}\right)\left(W_{(k+1)\Delta_{l}}-W_{k\Delta_{l}}\right)\right.$
$\displaystyle+\sum_{k=1}^{n_{t}}\nabla_{\theta}\log\left(g_{\theta}(x_{s_{k}},y_{s_{k}})\lambda_{\theta}(x_{s_{k}})\right)-\sum_{k=0}^{t\Delta_{l}^{-1}-1}\nabla_{\theta}\lambda_{\theta}(x_{k\Delta_{l}})\Delta_{l}$

and recall the interpolation (4). Now, we approximate the score function as

$\displaystyle\nabla_{\theta}\log\left(\gamma_{t,\theta}^{l}(1)\right):=\frac{{\mathbb{E}}_{\theta}\left[\zeta_{t,\theta}^{l}\left(\widetilde{X}_{0},\widetilde{X}_{\Delta_{l}},\ldots,\widetilde{X}_{t}\right)\left(\prod_{i=1}^{n_{t}}g_{\theta}(\widetilde{X}_{s_{i}},y_{s_{i}})\lambda_{\theta}(\widetilde{X}_{s_{i}})\right)\exp\left\\{-\sum_{j=0}^{t\Delta_{l}^{-1}-1}\lambda_{\theta}(\widetilde{X}_{j\Delta_{l}})\Delta_{l}\right\\}\right]}{{\mathbb{E}}_{\theta}\left[\left(\prod_{i=1}^{n_{t}}g_{\theta}(\widetilde{X}_{s_{i}},y_{s_{i}})\lambda_{\theta}(\widetilde{X}_{s_{i}})\right)\exp\left\\{-\sum_{j=0}^{t\Delta_{l}^{-1}-1}\lambda_{\theta}(\widetilde{X}_{j\Delta_{l}})\Delta_{l}\right\\}\right]}.$

Following the proof of [4, Appendix B], one can establish convergence of
$\nabla_{\theta}\log\left(\gamma_{t,\theta}^{l}(1)\right)$ to
$\nabla_{\theta}\log\left(\gamma_{t,\theta}(1)\right)$ as
$l\rightarrow\infty$; we do not give the statement here as it is essentially
analogous to that in [4].

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

In the following, we define some objects to help us build the particle
smoother which will ultimately allow us to define an online algorithm for
estimating the score function. We shall suppose that $t\in\mathbb{N}$ is fixed
in the forthcoming description.

Now we can represent the discretized smoother for $p\in\mathbb{N}$ (see (6)
also)

$\displaystyle\overline{\pi}_{p,\theta}^{l,S}\big{(}d(u_{0}^{l},\dots,u_{p-1}^{l})\big{)}:=\frac{\big{(}\prod_{k=0}^{p-1}G_{k,\theta}^{l}(u_{k-1}^{l},u_{k}^{l})\big{)}\,\eta_{0,\theta}^{l}(du_{0}^{l})\prod_{k=1}^{p-1}M_{\theta}^{l}(u_{k-1}^{l},du_{k}^{l})}{\int_{E_{l}^{p}}\big{(}\prod_{k=0}^{p-1}G_{k,\theta}^{l}(u_{k-1}^{l},u_{k}^{l})\big{)}\,\eta_{0,\theta}^{l}(du_{0}^{l})\prod_{k=1}^{p-1}M_{\theta}^{l}(u_{k-1}^{l},du_{k}^{l})}.$
(21)

where we have added $\theta$ subscripts for the quantities $\eta_{0}^{l}$,
$M^{l}$ and $G_{k}^{l}$ as in Section 2.4.1. As in the previous subsection,
the Feynman-Kac structure of (21) allows a particle filter estimation. We aim
to estimate (17) online using backward smoothing; this is possible given the
additive structure of (18) in terms of the dependence of the variables
$u_{k}^{l}$ [11]. Such dependence will be shown in the following. For any
$p\in\mathbb{N}_{0}$, let $f_{p,\theta}^{l}:\mathbb{R}^{d_{x}\times
d_{x}}\rightarrow\mathbb{R}^{d_{\theta}}$ be defined as

$\displaystyle f_{p,\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+\sum_{k=n_{p\Delta_{l}}+1}^{n_{(p+1)\Delta_{l}}}\nabla_{\theta}\log\left(g_{\theta}({x}_{s_{k}},y_{s_{k}})\lambda_{\theta}({x}_{s_{k}})\right)-\nabla_{\theta}\lambda_{\theta}(x_{p\Delta_{l}})\Delta_{l}.$

Now we define for $p\in\\{0,\dots,t-1\\}$ so that

$\displaystyle\mu_{p,\theta}^{l}(u_{p-1}^{l},u_{p}^{l})$
$\displaystyle:=\sum_{k=0}^{\Delta_{l}^{-1}-1}f_{k,\theta}^{l}(x_{p+k\Delta_{l}},x_{p+(k+1)\Delta_{l}});$
$\displaystyle F_{t,\theta}^{l}(u_{0}^{l},\dots,u_{t-1}^{l})$
$\displaystyle:=\sum_{p=0}^{t-1}\mu_{p,\theta}^{l}(u_{p-1}^{l},u_{p}^{l})\,\,\Big{(}\equiv\zeta_{t,\theta}^{l}(x_{0},x_{\Delta_{l}},\dots,x_{t})\Big{)}.$
(22)

The discretized smoothing distribution can be written via the time-reversal
formula for hidden Markov models (see e.g. [4, Section 3.2]) as

$\displaystyle\overline{\pi}_{t,\theta}^{l,S}\big{(}d(u_{0}^{l},\dots,u_{t-1}^{l})\big{)}:=\overline{\pi}_{t,\theta}^{l}(du_{t-1}^{l})\prod_{k=1}^{t-1}B_{k,\theta,\overline{\pi}_{k,\theta}^{l}}^{l}(u_{k}^{l},du_{k-1}^{l}),$

where $\overline{\pi}_{t,\theta}^{l}$ is the path-wise filter as in (6) with
$\theta$ subscripts and the backward Markov kernel is defined as, where we
write the density of $M_{\theta}^{l}$ as $M_{\theta,d}^{l}$,

$\displaystyle B_{k,\theta,\overline{\pi}_{k,\theta}^{l}}^{l}(u_{k}^{l},$
$\displaystyle
du_{k-1}^{l}):=\frac{\overline{\pi}_{k,\theta}^{l}(du_{k-1}^{l})\,G_{k,\theta}^{l}(u_{k-1}^{l},u_{k}^{l})M_{\theta,d}^{l}(u_{k-1}^{l},u_{k}^{l})}{\overline{\pi}_{k,\theta}^{l}(G_{k,\theta}^{l}(\cdot,u_{k}^{l})M_{\theta,d}^{l}(\cdot,u_{k}^{l}))}$
(23)

with the notation

$\overline{\pi}_{k,\theta}^{l}(G_{k,\theta}^{l}(\cdot,u_{k}^{l})M_{\theta,d}^{l}(\cdot,u_{k}^{l})):=\int_{E_{l}}\overline{\pi}_{k,\theta}^{l}(du_{k-1}^{l})G_{k,\theta}^{l}(u_{k-1}^{l},u_{k}^{l})M_{\theta,d}^{l}(u_{k-1}^{l},u_{k}^{l}).$

Given the structure of both $G_{k,\theta}^{l}$ and $M_{\theta,d}^{l}$, instead
of depending on $u_{k}^{l}$, the backward kernel depends only on its first
element $x_{k+\Delta_{l}}$, and the cost of computing it does not depend on
the level $l$. For $k\in\left\\{0,1,\ldots,t\Delta_{l}^{-1}-1\right\\}$ we
define

$\overline{g}_{k,\theta}^{l}\left(x_{k\Delta_{l}},x_{(k+1)\Delta_{l}}\right):=\left(\prod_{i=n_{k\Delta_{l}}+1}^{n_{(k+1)\Delta_{l}}}g_{\theta}({x}_{s_{i}},y_{s_{i}})\lambda_{\theta}({x}_{s_{i}})\right)\exp\left\\{-\lambda_{\theta}(x_{k\Delta_{l}})\Delta_{l}\right\\},$

then, one has

$\displaystyle B_{k,\theta,\overline{\pi}_{k,\theta}^{l}}^{l}(u_{k}^{l},$
$\displaystyle
du_{k-1}^{l}):=\frac{\overline{\pi}_{k,\theta}^{l}(du_{k}^{l})\,\overline{g}_{k,\theta}^{l}(x_{k},x_{k+\Delta_{l}})m_{\theta}^{l}(x_{k},x_{k+\Delta_{l}})}{\overline{\pi}_{k,\theta}^{l}(\overline{g}_{k,\theta}^{l}(\cdot,x_{k+\Delta_{l}})m_{\theta}^{l}(\cdot,x_{k+\Delta_{l}}))},$
(24)

where $m_{\theta}^{l}$ is the transition of the Euler-Maruyama discretization
of equation (16). Now, we have the following representation of the (time-
discretized) score function

$\displaystyle\nabla_{\theta}\log(\gamma_{t,\theta}^{l}(1))=\int_{E_{l}^{t}}F_{t,\theta}^{l}(u_{0}^{l},\dots,u_{t-1}^{l})\,\overline{\pi}_{t-1,\theta}^{l,S}\big{(}d(u_{0}^{l},\dots,u_{t}^{l})\big{)}.$
(25)

We can exploit the structure of (22) to estimate (25), which is partially our
objective. The estimation is made online by computing particle estimators of
the backward kernel. For precise details of the online smoothing methods see
[11] for instance. The sequential estimation of the score function is detailed
in Algorithm 7.

Algorithm 7 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}(u_{-1}^{l,i},u_{0}^{l,i})\leavevmode\nobreak\
\mu_{0,\theta}^{l}(u_{-1}^{l,i},u_{0}^{l,i})}{\sum_{i=1}^{N}G_{0,\theta}^{l}(u_{-1}^{l,i},u_{0}^{l,i})}$

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

  2. 2.

(Resampling step) 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})=\mu_{0,\theta}^{l}(x_{*},\check{u}_{0}^{l,i})$.

  3. 3.

(Sampling step) 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}\overline{g}_{k,\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})\\{F_{k-1,\theta}^{l,N}(\check{u}_{k-1}^{l,j})+\mu_{k,\theta}^{l}(\check{u}_{k-1}^{l,j},u_{k}^{l,i})\\}}{\sum_{j=1}^{N}\overline{g}_{k,\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})m_{\theta}^{l}(\check{x}_{k}^{j},x_{k+\Delta_{l}}^{i})}.$
(26)

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})}.$
(27)

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

The cost of computing the estimation of the score function given in (27) is
the cost of the particle filter $\mathcal{O}(\Delta^{-1}_{l}N)$ plus the cost
of computing the term $F_{k,\theta}^{l,N}$ which is $\mathcal{O}(N^{2})$,
thus, the whole cost is $\mathcal{O}(\Delta^{-1}_{l}N+N^{2})$ per unit time.
As pointed out in [4], applying multilevel techniques will not improve the
complexity of the score function, therefore, we only apply a single-level
particle filter. We note that in [4] a further extension using diffusion
bridges is possible and could be used in the context of this paper; that
method has the advantage of a so-called path-space interpretation. One can
also consider the methods in [15], although it is unclear on its efficacy for
the class of models that we consider.

##  3 Mathematical Analysis

  * (D3)

For any $y\in\mathsf{Y}$ there exists a $0<C<+\infty$ such that
$\inf_{x\in\mathbb{R}^{d}}g(x,y)\geq C$.

A sketch proof of the below result can be found in the Appendix. We use the
notation
$[\pi_{t}^{l}-\pi_{t}^{l-1}](\varphi)=\pi_{t}^{l}(\varphi)-\pi_{t}^{l-1}(\varphi)$.

######  Proposition 3.1.

Assume (D(D1)-(D3)). Then for any $t\in\mathbb{N}$ there exists a $C<+\infty$
such that for any
$(l,N,\varphi)\in\mathbb{N}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})\cap\textrm{\emph{Lip}}(\mathbb{R}^{d_{x}})$

$\mathbb{E}\left[\Big{(}[\pi_{t}^{l}-\pi_{t}^{l-1}]^{N}(\varphi)-[\pi_{t}^{l}-\pi_{t}^{l-1}](\varphi)\Big{)}^{2}\right]\leq\frac{C(\|\varphi\|+\|\varphi\|_{\textrm{\emph{Lip}}})^{2}\Delta_{l}^{1/2}}{N}.$

The following result can be established via Lemma A.2 and the proofs for
Proposition 3.3.

######  Proposition 3.2.

Assume (D(D1)-(D3)). Then for any $t\in\mathbb{N}$ there exists a $C<+\infty$
such that for any
$(l,N,\varphi)\in\mathbb{N}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})\cap\textrm{\emph{Lip}}(\mathbb{R}^{d_{x}})$

$\Big{|}\mathbb{E}\left[[\pi_{t}^{l}-\pi_{t}^{l-1}]^{N}(\varphi)-[\pi_{t}^{l}-\pi_{t}^{l-1}](\varphi)\right]\Big{|}\leq\frac{C(\|\varphi\|+\|\varphi\|_{\textrm{\emph{Lip}}})\Delta_{l}^{1/4}}{N}.$

Given the above two results, one can show that for the MLPF estimator in (11),
that the associated MSE (when centering by $\pi_{t}(\varphi)$) is, for
$\epsilon>0$ given, of $\mathcal{O}(\epsilon^{2})$ with a computational effort
of $\mathcal{O}(\epsilon^{-2.5})$, with a selection of
$N_{l}=\mathcal{O}(\epsilon^{-2.5}\Delta_{l}^{3/4})$ with $L$ chosen so that
$\mathcal{O}(\Delta_{L})=\mathcal{O}(\epsilon)$. This can be inferred by using
standard arguments such as in [21].

In addition, one has the following result which follows directly from the
above results and the proofs in [20, Theorem 2].

######  Proposition 3.3.

Assume (D(D1)-(D3)). Then for any
$(t,\varphi)\in\mathbb{N}\cap\mathcal{B}_{b}(\mathbb{R}^{d_{x}})\cap\textrm{\emph{Lip}}(\mathbb{R}^{d_{x}})$
there exist choices of $\mathbb{P}(p|l)$ and $\mathbb{P}_{L}(l)$ so that (12)
with $\Xi_{L}$ (resp. $\Xi_{0}$) as in (14) (resp. (13)) is an unbiased and
finite variance estimator of $\eta_{t}(\varphi)$.

We note that the choice of $\mathbb{P}(p|l)$ and $\mathbb{P}_{L}(l)$ are
exactly as [20] so we do not discuss it here. A discussion of the cost is the
same as [20, Section 3.2] and hence omitted.

##  4 Numerical Results

###  4.1 Models

In the following, we introduce four different diffusion models (through all
the numerical examples in this section and the following we set $d_{x}=1$).
Let $X_{0}=x^{*}$ be the initial point at time $t=0$, the diffusion processes
that will be associated with the hidden Markov process are given through the
following stochastic differential equations (SDEs)

  * •

Ornstein-Uhlenbeck (OU): Let $\nu,$ $\sigma\in\mathbb{R}^{+}$. For $t\geq 0$,
an Ornstein-Uhlenbeck process is defined as

$\displaystyle dX_{t}=-\nu X_{t}dt+\sigma dW_{t}.$ (28)

  * •

Langevin process: Let $\pi$ be a probability density function; the overdamped
Langevin SDE is defined as

$\displaystyle dX_{t}=-\nabla\log\left(\pi(X_{t})\right)dt+dW_{t}.$

The asymptotic distribution of this process is precisely $\pi$. In the
numerical experiments we use Student’s t-distribution with $\nu=10$ degrees of
freedom, thus, $\nabla\log(\pi(x))=-(\nu+1)x(x^{2}+\nu)^{-1}$. Therefore, we
have

$\displaystyle dX_{t}=(\nu+1)X_{t}(X_{t}^{2}+\nu)^{-1}dt+dW_{t}.$ (29)

  * •

Nonlinear diffusion term (NLDT): The SDE of this process is defined as

$\displaystyle dX_{t}=\frac{1}{\sqrt{1+X_{t}^{2}}}dW_{t}.$ (30)

  * •

Geometric Brownian Motion (GBM): Let $\nu\in\mathbb{R}$ and
$\sigma\in\mathbb{R}^{+}$, the GBM is defined by the SDE

$\displaystyle dX_{t}=\nu X_{t}dt+\sigma X_{t}dW_{t},$ (31)

with $x^{*}\geq 0$.

The likelihood $g(x,\cdot)$ is chosen to be the normal distribution with
variance $\Sigma\in\mathbb{R}^{+}$ and mean $x\in\mathbb{R}$. The intensity
function is taken as $\lambda(x)=a|x|$, where $a\in\mathbb{R}^{+}$. The
observations and the times at which they are observed are generated from
running one realization of each of the SDEs above. The parameters in the
models above are chosen so that two criteria are met: the first one is the
stability of the particle filter for a feasible level of discretization and a
certain number of particles, and the second is the strength of the coupling,
i.e., how small is the constant $C$ in Propositions 3.1 compared to the
variance of the single-level particle filter. The code of the simulations is
written in Python and it can be downloaded from
https://github.com/maabs/Multilevel-for-Diffusions-Observed-via-Marked-Point-
Processes.

###  4.2 Multilevel Particle Filter

In Figure 1, we plot the computational cost of the MLPF and compare it with
the cost of the single-level PF. We can clearly see the predicted error-to-
cost rates of each algorithm and the reduction in cost when using the MLPF
compared to using the single-level PF. We choose the final time $t=100$. The
number of particles that we need in order to attain the desired MSE
$\varepsilon^{2}$ in the MLPF is
$N_{l}=\mathcal{O}(\Delta_{l}^{3/4}\varepsilon^{-(2.5)})$, for
$l\in\\{l_{0},l_{0}+1,\cdots,L\\}$, where $l_{0}\in\mathbb{N}_{0}$.
$L=\mathcal{O}(\log(\varepsilon^{-1}))\in\mathbb{N}$, $L>l_{0}$, depends on
the desired level of error.

Figure 1: Computional complexity of the multilevel particle filter and the
single-level particle filter. On the horizontal line we set MSE
$=\varepsilon^{2}$.

###  4.3 Unbiased Particle Filter

Similar to the MLPF, we plot the cost vs. MSE of the unbiased estimator in
(15) at the final time $t=100$. We use
$\mathbb{P}_{L}(l)\propto\log(l+2)^{2}(l+1)\Delta^{\frac{1}{2}}_{l}$ as the
probability mass distribution of the discretization level $l$, where
$l\in\\{l_{0},l_{0}+1,\cdots,L_{\text{Trunc}}\\}$. For the number of particles
we use a geometric sequence $N_{p}=N_{0}2^{p}$, where the randomization of the
particles (in terms of $p$) follows the probability mass distribution
$\mathbb{P}_{P}(p)\propto\log(p+2)^{2}(p+1)/N_{p}$ for
$p\in\\{0,1,\cdots,P_{\text{Trunc}}\\}$. $L_{\text{Trunc}}$ and
$P_{\text{Trunc}}$ are the truncation levels of our now unbiased algorithm; we
choose them large enough to ensure the bias is negligible w.r.t. the variance
of the estimator. For more details on these choices, see [20]. We use
different values of $L_{\text{Trunc}}$, $P_{\text{Trunc}}$ and $N_{0}$
depending on the diffusion process, for the OU, Langevin and NLDT the values
are $L_{\text{Trunc}}=10$, $P_{\text{Trunc}}=11$ and $N_{0}=5$; for the GBM we
have $L_{\text{Trunc}}=10$, $P_{\text{Trunc}}=5$ and $N_{0}=100$; the value of
$l_{0}=0$ is the same for all the processes.

The computational complexity of the unbiased estimator in (15) is plotted in
Figure 2. Due to the random nature of the estimator, we can observe different
rates depending on the number of realizations $M$. In this Section, we are
interested in the rates rather than the comparison to the single-level PF or
MLPF; this is because one of the perks of the unbiased estimator comes from
its embarrassingly parallel characteristic, thus making the unbiased estimator
viable depending on the parallelization capabilities of the hardware. A
comparison between the unbiased and multilevel particle filters of partially
observed diffusions at regular times can be found in [20].

Figure 2: Computational complexity of the unbiased particle filters.

###  4.4 Parameter Estimation

We use Algorithm 7 along with (20) to estimate
$\nabla_{\theta}\log\left(\gamma_{c(k+1),\theta_{0:k}}(1)\right)$ and the
parameters $\theta$ for each iteration $k$, specifically, for each $k\geq 0$
we replace $\theta$ by $\theta_{k}$ in steps 2 and 3 of Algorithm 7. In this
Section we estimate the parameters for the four different hidden diffusions
given above. Also we estimate the parameter
$\theta_{\lambda}\in\mathbb{R}^{+}$ in the intensity function
$\lambda_{\theta}(x)=\theta_{\lambda}|x|$, and the and variance parameter
$\theta_{\Sigma}\in\mathbb{R}^{+}$ in the likelihood function
$g_{\theta}(\cdot,y)$, which is defined as a Gaussian density with mean $y$
and variance $\theta_{\Sigma}$. The hidden processes and their drift terms
are:

  * •

Ornstein-Uhlenbeck process in (28), with $b_{\theta}(x)=-\theta_{b}$ with
$\theta_{b}\in\mathbb{R}^{+}$. We will estimate
$\theta=(\theta_{b},\theta_{\lambda},\theta_{\Sigma})$.

  * •

Langevin process in (29) with $b_{\theta}(x)=(\nu+1)x(\nu+x^{2})^{-1}$,
$\nu\in\mathbb{R}^{+}$. We will estimate
$\theta=(\theta_{\lambda},\theta_{\Sigma})$.

  * •

The nonlinear diffusion process in (30) with $b_{\theta}(x)=0$. We will
estimate $\theta=(\theta_{\lambda},\theta_{\Sigma})$.

  * •

The GBM process in (31) with $b_{\theta}(x)=\theta_{b}x$ with
$\theta_{b}\in\mathbb{R}$. We will estimate
$\theta=(\theta_{b},\theta_{\lambda},\theta_{\Sigma})$.

In order to control the trade-off between the number of new observations
included in the estimation of the score function and the resampling rate of
the particle filter, we introduce the hyperparameter $c$ in (20). Different
step-sizes are used depending on the parameter we estimate, i.e., we modify
slightly (20) changing $\alpha_{m}$ from a scalar to a positive definite
diagonal matrix in $\mathbb{R}^{d_{\theta}\times d_{\theta}}$, with diagonal
terms $\alpha_{m}^{(i)}=\alpha_{0}^{(i)}(m+1)^{-\beta},$ $\beta\in(0.5,1]$,
$\alpha_{0}^{(i)}\in\mathbb{R}^{+}$, $i\in\\{1,2,3\\}$, $m\in\mathbb{N}_{0}$.
Plots of the evolving values of $\theta_{m}$ are provided in Figure 3. For the
OU process we choose the true parameters to be
$(\overline{\theta}_{\lambda},\overline{\theta}_{\Sigma},\overline{\theta}_{b})=(3.5,1,0.98)$
and the initial guesses
$({\theta}_{\lambda,0},{\theta}_{\Sigma,0},{\theta}_{b,0})=(1.5,1.5,0.48)$.
For the Langevin process the true parameters are
$(\overline{\theta}_{\lambda},\overline{\theta}_{\Sigma})=(1,1)$ and the
initial values $({\theta}_{\lambda,0},{\theta}_{\Sigma,0})=(2,2.5)$. For the
nonlinear diffusion process the true parameters are
$(\overline{\theta}_{\lambda},\overline{\theta}_{\Sigma})=(0.222,1)$ and the
initial values $({\theta}_{\lambda,0},{\theta}_{\Sigma,0})=(2.222,2)$.
Finally, for the GBM the true parameters are
$(\overline{\theta}_{\lambda},\overline{\theta}_{\Sigma},\overline{\theta}_{b})=(0.5,1,0.015)$
and $({\theta}_{\lambda,0},{\theta}_{\Sigma,0},{\theta}_{b,0})=(2.5,2,1.015)$.
We can see in the figures the relatively fast convergence to the true
parameters in all models considered.

Figure 3: Estimated parameters in terms of the iterations for the four
different processes.

#### Acknowledgements

All authors were supported by KAUST baseline funding.

##  Appendix A Mathematical Proofs

The proof of Proposition 3.3 is virtually identical to that of [19, Theorem
4.1]. The only issue is that there is minor mistake in the proof of [19, Lemma
A.1.] and we correct that here. As the rest of the proof is then basically as
[19, Theorem 4.1] we do not repeat the technical proofs. The reason that the
proofs are so similar, is that our hidden process is identical to that of [19]
and the potential functions $\mathbf{G}_{t}^{l}$ (defined below) are uniformly
in $l$ upper and lower bounded by deterministic constants as well having
similar functional form to the counter-parts in [19]. To that end, we give a
similar (and corrected) result to [19, Lemma A.1.] and leave the remainder of
the proof as a simple exercise in reading and adapting the results in [19].

###  A.1 Notations

To assist moving between the proofs of this paper and that of [19] we adopt a
similar notations to that paper, which are different from the main text.

For $p\in\mathbb{N}_{0}$ set

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

For technical reasons, this notation differs slightly from that in Section
2.4.1.

For $\varphi\in\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$, we define, for any $l\geq
0$, $\boldsymbol{\varphi}^{l}:E_{l}\rightarrow\mathbb{R}$

$\boldsymbol{\varphi}^{l}(x_{0},x_{\Delta_{l}},\dots,x_{1}):=\varphi(x_{1}).$

Set, with $p\in\mathbb{N}_{0}$

$\mathbf{G}_{p}^{l}(u_{p}^{l}):=\left(\prod_{k=n_{p}+1}^{n_{p+1}}g(x_{s_{k}},y_{s_{k}})\lambda(x_{s_{k}})\right)\exp\left\\{-\Delta_{l}\sum_{k=0}^{\Delta_{l}^{-1}-1}\lambda(x_{p+(k+1)\Delta_{l}})\right\\}.$

Denote by $M^{l}:\mathbb{R}^{d_{x}}\rightarrow\mathcal{P}(E_{l})$ the joint
Markov transition of $(x_{0},x_{\Delta_{l}},\dots,x_{1})$ defined via the
Euler discretization and a Dirac on a point $x\in\mathbb{R}^{d_{x}}$: for
$(x,\varphi)\in\mathbb{R}^{d_{x}}\times\mathcal{B}_{b}(E_{l})$,

$M^{l}(\varphi)(x):=\int_{E_{l}}\varphi(x_{0},x_{\Delta_{l}},\dots,x_{1})\delta_{x}(dx_{0})\Big{[}\prod_{k=1}^{\Delta_{l}^{-1}}m(x_{(k-1)\Delta_{l}},x_{k\Delta_{l}})\Big{]}d(x_{\Delta_{l}},\dots,x_{1}).$

For $p\in\mathbb{N}$, we use the notation

$\mu(\mathbf{G}_{p-1}^{l}\mathbf{M}^{l}(\varphi))=\int_{E_{l}}\mu(d(x_{p-1},x_{p-1+\Delta_{l}},\dots,x_{p}))\mathbf{G}_{p-1}^{l}(x_{p-1},x_{p-1+\Delta_{l}},\dots,x_{p-\Delta_{l}},x_{p})M^{l}(\varphi)(x_{p})$

where $(\mu,\varphi)\in\mathcal{P}(E_{l})\times\mathcal{B}_{b}(E_{l})$. For
$p\in\mathbb{N}$, define the operator
$\Phi_{p}^{l}:\mathcal{P}(E_{l})\rightarrow\mathcal{P}(E_{l})$ with
$(\mu,\varphi)\in\mathcal{P}(E_{l})\times\mathcal{B}_{b}(E_{l})$ as:

$\Phi_{p}^{l}(\mu)(\varphi):=\frac{\mu(\mathbf{G}_{p-1}^{l}\mathbf{M}^{l}(\varphi))}{\mu(\mathbf{G}_{p-1}^{l})}.$

Now, define, for $(p,\varphi)\in\mathbb{N}_{0}\times\mathcal{B}_{b}(E_{l})$,

$\eta_{p}^{l}(\varphi):=\int_{\mathbb{R}^{d_{x}}}\eta_{p}^{l}(dx)\int_{E_{l}}M^{l}(x,du)\varphi(u).$

Then one can establish that for $p\in\mathbb{N}$

$\eta_{p}^{l}(\varphi)=\Phi_{p}^{l}(\eta_{p-1}^{l})(\varphi).$

Moreover, for
$(t,\varphi)\in\mathbb{N}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\pi_{t}^{l}(\varphi)=\frac{\eta_{t-1}^{l}(\mathbf{G}_{t-1}^{l}\boldsymbol{\varphi}^{l})}{\eta_{t-1}^{l}(\mathbf{G}_{t-1}^{l})}$

which is the time discretized filter.

Some operators are now defined. Let $(l,p,n)\in\mathbb{N}_{0}^{3}$, $n>p$,
$(x,\varphi)\in E_{l}\times\mathcal{B}_{b}(E_{l})$

$\mathbf{Q}_{p,n}^{l}(\varphi)(u_{p}):=\int\varphi(u_{n})\Big{(}\prod_{q=p}^{n-1}\mathbf{G}_{q}^{l}(u_{q})\Big{)}\prod_{q=p+1}^{n}M^{l}(u_{q-1},du_{q}).$

where we use the convention
$\mathbf{Q}_{p,p}^{l}(\varphi)(u_{p})=\varphi(u_{p})$. In addition, for
$(l,p,n)\in\mathbb{N}_{0}^{3}$, $n>p$, $(x,\varphi)\in
E_{l}\times\mathcal{B}_{b}(E_{l})$:

$\mathbf{D}_{p,n}^{l}(\varphi)(u_{p}):=\frac{\mathbf{Q}_{p,n}^{l}(\varphi-\eta_{n}^{l}(\varphi))(u_{p})}{\pi_{p}^{l}(\mathbf{Q}_{p,n}^{l}(1))}$

where
$\mathbf{D}_{p,p}^{l}(\varphi)(u_{p})=\varphi(u_{p})-\eta_{p}^{l}(\varphi)$.
Throughout our arguments, $C$ is a finite constant whose value may change from
line to line, but does not depend upon $l$ nor $N$. The particular
dependencies of a given constant will be clear from the statement of a given
result.

Set, for $l\in\mathbb{N}$,
$(n,p,\varphi)\in\mathbb{N}_{0}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$,
$p<n$

$\displaystyle
T_{p,n}^{l}(\varphi):=\mathbb{E}[(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\|\varphi\|^{2}\mathbb{E}[(\mathbf{G}_{p}^{l}(U_{p}^{l,1})-\mathbf{G}_{p}^{l-1}(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\Delta_{l}^{2}$

and if $p=n$

$T_{p,n}^{l}(\varphi):=\mathbb{E}[(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\Delta_{l}^{2}.$

We remark that using these modified operators in the analysis of the
multilevel and unbiased methods do not change the algorithms and are just an
alternative representation. In the notation just introduced, we use
$\check{\eta}_{t}^{l,N}$ and $\check{\bar{\eta}}_{t}^{l-1,N}$ to denote the
$N-$empirical measures from the coupled particle filter at levels $l$ and
$l-1$ associated to the samples $u_{n}^{l,1:N}$ and $\bar{u}_{n}^{l-1,1:N}$.

###  A.2 Technical Results

######  Lemma A.1.

Assume (D(D1)-(D3)). Then for any $n\in\mathbb{N}$ there exists a $C<+\infty$
such that for any
$(l,N,\varphi)\in\mathbb{N}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\mathbb{E}\left[\left(\check{\eta}_{n}^{l,N}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\check{\bar{\eta}}_{n}^{l-1,N}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})-\\{\eta_{n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\eta_{n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})\\}\right)^{2}\right]\leq\frac{C}{N}\sum_{p=0}^{n}T_{p,n}^{l}(\varphi).$

###### Proof.

As in [19, Lemma A.1.] we apply the following standard Martingale plus
remainder decomposition [9, Lemma 6.3] followed by the $C_{2}-$inequality
multiple times to yield the upper-bound

$\mathbb{E}\left[\left(\check{\eta}_{n}^{l,N}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\check{\bar{\eta}}_{n}^{l-1,N}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})-\\{\eta_{n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\eta_{n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})\\}\right)^{2}\right]\leq
C\Big{(}\sum_{p=0}^{n}\mathbb{E}[T_{1}(p)^{2}]+\sum_{p=0}^{n-1}\mathbb{E}[T_{2}(p)^{2}]\Big{)}$

where

$\displaystyle T_{1}(p)$ $\displaystyle:=$
$\displaystyle(\check{\eta}_{p}^{l,N}-\Phi_{p}^{l}(\check{\eta}_{p-1}^{l,N}))(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))-(\check{\bar{\eta}}_{p}^{l-1,N}-\Phi_{p}^{l-1}(\check{\bar{\eta}}_{p-1}^{l-1,N}))(\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1}))$
$\displaystyle T_{2}(p)$ $\displaystyle:=$
$\displaystyle\left(\frac{\check{\eta}_{p}^{l,N}(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))}{\check{\eta}_{p}^{l,N}(\mathbf{G}_{p}^{l})}[\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})-\frac{\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1}))}{\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})}[\eta_{p}^{l-1}-\check{\bar{\eta}}_{p}^{l-1,N}](\mathbf{G}_{p}^{l-1})\right).$

We have to control the terms $T_{1}(p)$, $p\in\\{0,1,\dots,n\\}$ and
$T_{2}(p)$, $p\in\\{0,1,\dots,n-1\\}$ in an appropriate way.

For $T_{1}(p)$ applying the conditional Marcinkiewicz-Zygmund inequality
followed by Jensen’s inequality

$\mathbb{E}[T_{1}(p)^{2}]\leq\frac{C}{N}\overline{\mathbb{E}}[(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}.$
(32)

For $T_{2}(p)$ we have

$T_{2}(p)=T_{3}(p)+T_{4}(p)+T_{5}(p),$

where

$\displaystyle T_{3}(p)$ $\displaystyle:=$
$\displaystyle[\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})\frac{\check{\eta}_{p}^{l,N}(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))}{\check{\eta}_{p}^{l,N}(\mathbf{G}_{p}^{l})\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})}\Big{\\{}\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})-\check{\eta}_{p}^{l,N}(\mathbf{G}_{p}^{l})\Big{\\}}$
(33) $\displaystyle T_{4}(p)$ $\displaystyle:=$
$\displaystyle[\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})\frac{1}{\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})}\Big{\\{}\check{\eta}_{p}^{l,N}(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))-\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1}))\Big{\\}}$
(34) $\displaystyle T_{5}(p)$ $\displaystyle:=$
$\displaystyle\frac{\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1}))}{\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})}\left([\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})-[\eta_{p}^{l-1}-\check{\bar{\eta}}_{p}^{l-1,N}](\mathbf{G}_{p}^{l-1})\right).$
(35)

Using the $C_{2}-$inequality, we need only to bound the second moment of each
of the terms $T_{3}(p)$, $T_{4}(p)$ and $T_{5}(p)$ to conclude the proof.

For $T_{3}(p)$ applying Cauchy-Schwarz twice and using the uniform in $l$
lower-bounds on $\mathbf{G}_{p}^{l}$ and $\mathbf{G}_{p}^{l-1}$ yields that

$\mathbb{E}[T_{3}(p)^{2}]\leq
C\mathbb{E}[|[\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})|^{8}]^{1/4}\mathbb{E}[|\check{\eta}_{p}^{l,N}(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))|^{8}]^{1/4}\mathbb{E}[|\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{G}_{p}^{l-1})-\check{\eta}_{p}^{l,N}(\mathbf{G}_{p}^{l})|^{4}]^{1/2}.$

The first two terms on the R.H.S. are, uniformly in $l$, $\mathcal{O}(N^{-1})$
(see e.g. the proofs in [19, Appendix A.5]) and the last term is upper-bounded
by
$\mathbb{E}[(\mathbf{G}_{p}^{l}(U_{p}^{l,1})-\mathbf{G}_{p}^{l-1}(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}$,
which gives

$\mathbb{E}[T_{3}(p)^{2}]\leq
C\frac{\|\varphi\|^{2}}{N^{2}}\mathbb{E}[(\mathbf{G}_{p}^{l}(U_{p}^{l,1})-\mathbf{G}_{p}^{l-1}(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}.$

For $T_{4}(p)$ applying Cauchy-Schwarz and using the uniform in $l$ lower-
bound on $G_{p}^{l-1}$ gives

$\mathbb{E}[T_{4}(p)^{2}]\leq
C\mathbb{E}[|[\eta_{p}^{l}-\check{\eta}_{p}^{l,N}](\mathbf{G}_{p}^{l})|^{4}]^{1/2}\mathbb{E}[|\check{\eta}_{p}^{l,N}(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l}))-\check{\bar{\eta}}_{p}^{l-1,N}(\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1}))|^{4}]^{1/2},$

thus, it easily follows that

$\mathbb{E}[T_{4}(p)^{2}]\leq\frac{C}{N}\mathbb{E}[(\mathbf{D}_{p,n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}.$

One can use similar arguments to the above to establish that

$\mathbb{E}[T_{5}(p)^{2}]\leq
C\frac{\|\varphi\|^{2}}{N}\left(\mathbb{E}[(\mathbf{G}_{p}^{l}(U_{p}^{l,1})-\mathbf{G}_{p}^{l-1}(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\Delta_{l}^{2}\right)$

where the $\Delta_{l}^{2}$ term is from the weak error of the filter/predictor
(see Proposition 2.1). Therefore, we have shown that

$\mathbb{E}[T_{2}(p)^{2}]\leq CT_{p,n}^{l}(\varphi).$ (36)

Combining (32) and (36) the proof can easily be concluded. ∎

######  Remark A.1.

Under (D(D1)-(D3)) one can also prove that: for any $(n,q)\in\mathbb{N}^{2}$,
$n\leq q$ there exists a $C<+\infty$ such that for any
$(l,N,\varphi)\in\mathbb{N}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\mathbb{E}\left[\left(\check{\eta}_{n}^{l,N}(\mathbf{D}_{n,q}^{l}(\mathbf{G}_{q}^{l}\boldsymbol{\varphi}^{l}))-\check{\bar{\eta}}_{n}^{l-1,N}(\mathbf{D}_{n,q}^{l-1}(\mathbf{G}_{q}^{l-1}\boldsymbol{\varphi}^{l-1}))-\\{\eta_{n}^{l}(\mathbf{D}_{n,q}^{l}(\mathbf{G}_{q}^{l}\boldsymbol{\varphi}^{l}))-\eta_{n}^{l-1}(\mathbf{D}_{n,q}^{l-1}(\mathbf{G}_{q}^{l-1}\boldsymbol{\varphi}^{l-1}))\\}\right)^{2}\right]\leq\frac{C}{N}\sum_{p=0}^{n}T_{p,n,q}^{l}(\varphi),$

where for $p<n$

$\displaystyle
T_{p,n,q}^{l}(\varphi):=\mathbb{E}[(\mathbf{D}_{p,q}^{l}(\mathbf{G}_{q}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,q}^{l-1}(\mathbf{G}_{q}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\|\varphi\|^{2}\mathbb{E}[(\mathbf{G}_{p}^{l}(U_{p}^{l,1})-\mathbf{G}_{p}^{l-1}(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\Delta_{l}^{2}$

and if $p=n$

$T_{p,n,q}^{l}(\varphi):=\mathbb{E}[(\mathbf{D}_{p,q}^{l}(\mathbf{G}_{q}^{l}\boldsymbol{\varphi}^{l})(U_{p}^{l,1})-\mathbf{D}_{p,q}^{l-1}(\mathbf{G}_{q}^{l-1}\boldsymbol{\varphi}^{l-1})(\bar{U}_{p}^{l-1,1}))^{4}]^{1/2}+\Delta_{l}^{2}.$

######  Lemma A.2.

Assume (D(D1)-(D3)). Then for any $n\in\mathbb{N}$ there exists a $C<+\infty$
such that for any
$(l,N,\varphi)\in\mathbb{N}^{2}\times\mathcal{B}_{b}(\mathbb{R}^{d_{x}})$

$\left|\mathbb{E}\left[\check{\eta}_{n}^{l,N}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\check{\bar{\eta}}_{n}^{l-1,N}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})-\\{\eta_{n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\eta_{n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})\\}\right]\right|\leq$
$\frac{C}{N}\sum_{p=0}^{n}\left[\left(\sum_{q=0}^{p}T_{q,p}^{l}(\varphi)\right)^{1/2}+\left(\sum_{p=0}^{n}T_{p,n,q}^{l}(\varphi)\right)^{1/2}\right].$

###### Proof.

Using the Martingale plus remainder decomposition we have that

$\left|\mathbb{E}\left[\check{\eta}_{n}^{l,N}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\check{\bar{\eta}}_{n}^{l-1,N}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})-\\{\eta_{n}^{l}(\mathbf{G}_{n}^{l}\boldsymbol{\varphi}^{l})-\eta_{n}^{l-1}(\mathbf{G}_{n}^{l-1}\boldsymbol{\varphi}^{l-1})\\}\right]\right|\leq\sum_{j=3}^{5}\mathbb{E}[T_{j}(p)]$

where the terms $T_{3:5}(p)$ are given in (33)-(35). The proof essentially
follows that of Lemma A.1 except that one will use Lemma A.1 and Remark A.1
where it is appropriate; as the proof is fairly trivial, therefore, it is
omitted for brevity. ∎

## References

  * [1] Archibald, R., Bao, F. & Tu, X. (2019). A direct filter method for parameter estimation. _J. Comput. Phys._ , 398, pp. 108871.
  * [2] Bain, A. & Crisan, D. (2009). _Fundamentals of Stochastic Filtering_. Springer: New York.
  * [3] Benveniste, A., Métivier, M. & Priouret, P. (1990). _Adaptive Algorithms and Stochastic Approximation_. New York: Springer-Verlag.
  * [4] Beskos, A., Crisan, D., Jasra, A., Kantas, N., & Ruzayqat, H. (2021). Score-based parameter estimation for a class of continuous-time state space models. _SIAM J. Sci. Comp._ , 43, pp. A2555–A2580.
  * [5] Chada, N., Franks, J., Jasra A., Law K., & Vihola M. (2021). Unbiased inference for discretely observed hidden Markov model diffusions. _SIAM/ASA JUQ_ , 9, pp. 763–787.
  * [6] Del Moral, P. (2004). _Feynman-Kac Formulae_. Springer: New York.
  * [7] Del Moral, P., Hu, S., Jasra, A., Ruzayqat, H. & Wang, X. (2023). Bayesian Parameter Inference for Partially Observed Diffusions using Multilevel Stochastic Runge-Kutta Methods. _ArXiv preprint_ arXiv: 2309.13557.
  * [8] Del Moral, P., Jacod, J., & Protter, P. (2001). The Monte Carlo method for filtering with discrete-time observations. _Probab. Theory Rel. Fields_ , 120, pp. 346–368.
  * [9] Del Moral, P., Doucet, A. & Jasra, A. (2012). On adaptive resampling procedures for sequential Monte Carlo methods. Bernoulli, 18, pp. 252–272.
  * [10] Del Moral, P., Doucet, A. & Singh S. S. (2010). A backward particle interpretation of Feynman-Kac formuale. _M2AN_ , 44, pp. 947–975.
  * [11] Del Moral, P., Doucet, A., & Singh S. S. (2010). Forward smoothing using sequential Monte Carlo. _ArXiv preprint_ arXiv:1012.5390.
  * [12] Fearnhead, P., Papaspiliopoulos, O. & Roberts, G. O. (2008). Particle filters for partially observed diffusions. _J. R. Stat. Soc. Ser. B_ , 70, pp. 755–777.
  * [13] Giles, M. B. (2008). Multilevel Monte Carlo path simulation. _Op. Res._ , 56, pp. 607–617.
  * [14] Giles, M. B. (2015) Multilevel Monte Carlo methods. _Acta Numerica_ , 24, pp. 259–328.
  * [15] Hai-Dang, D. & Chopin, N. (2022). On the complexity of backward smoothing algorithms. _ArXiv preprint_ arXiv:2207.00976.
  * [16] Heinrich, S. (2001). Multilevel Monte Carlo methods. In _Large-Scale Scientific Computing_ , (eds. S. Margenov, J. Wasniewski & P. Yalamov), Springer: Berlin.
  * [17] Jasra, A. & Yu, F. (2020). Central limit theorems for coupled particle filters. _Adv. Appl. Probab._ , 52, pp. 942–1001.
  * [18] Jasra, A., Law K. J. H. & Suciu, C. (2020). Advanced Multilevel Monte Carlo. _Intl. Stat. Rev._ , 88, pp. 548–579.
  * [19] Jasra, A., Yu, F. & Heng, J. (2020). Multilevel particle filters for the non-linear filtering problem in continuous time. _Stat. Comp._ , 30, pp. 1381–1402.
  * [20] Jasra, A., Law, K. J. H. & Yu, F. (2022). Unbiased filtering of a class of partially observed diffusions. _Adv. Appl. Probab._ , 54, pp. 661–687.
  * [21] Jasra, A., Kamatani, K., Law K. J. H. & Zhou, Y. (2017). Multilevel particle filters. _SIAM J. Numer. Anal._ , 55, pp. 3068-3096.
  * [22] Jasra, A. , Kamatani K., Law, K. & Zhou, Y. (2018). Bayesian Static Parameter Estimation for Partially Observed Diffusions via Multilevel Monte Carlo. _SIAM J. Sci. Comp._ , 40, pp. A887–A902.
  * [23] Jasra, A., Kamatani, K., Osei, P. P. & Zhou, Y. (2018). Multilevel particle filters: normalizing constant estimation. _Statist. Comp._ , 28, pp. 47–60.
  * [24] Jasra, A., Ballesio, M., Von Schwerin, E. & Tempone, R. (2023). A coupled particle filter for multilevel estimation. _Stoch. Anal. Appl._ , 41, pp. 820–859.
  * [25] Jin, R., Singh, S. S. & Chopin, N. (2023). De-biasing particle filtering for a continuous time hidden Markov model with a Cox process observation model. _Stat. Sin._ , (to appear).
  * [26] Kantas, N., Doucet, A., Singh, S.S., Maciejowski, J. & Chopin, N. (2015). On particle methods for parameter estimation in state-space models. _Statistical science_ 30(3), pp. 328–351.
  * [27] Kloeden, P. & Platten, E. (1999). _Numerical Solution of Stochastic Differential Equations_. Third Edition. Springer: New York.
  * [28] LeGland, F. & Mevel, L., (1997). Recursive estimation in hidden Markov models. _Proceedings of the 36th IEEE Conference on Decision and Control_ , 4, pp. 3468–3473.
  * [29] Martin, J. M., Jasra, A. & McCoy, E. (2013). Inference for a class of partially observed point process models. _Ann. Inst. Stat. Math._ , 65, pp. 413–437.
  * [30] McLeish, D. (2011). A general method for debiasing a Monte Carlo estimator. _Monte Carlo Meth. Appl._ , 17, pp. 301–315.
  * [31] Rhee, C. H. & Glynn, P. (2015). Unbiased estimation with square root convergence for SDE models. _Op. Res._ , 63, pp. 1026–1043.
  * [32] Ruzayqat, H. & Jasra, A. (2022). Unbiased Parameter Inference for a Class of Partially Observed Lévy-Process Models. Found. data sci., 4 (2), pp. 299–322.
  * [33] 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), pp. 2814–2829.
  * [34] Vihola, M. (2018). Unbiased estimators and multilevel Monte Carlo. _Oper. Res._ , 66(2) pp.448–462.
  * [35] Wagner, W. (1988). Unbiased multi-step estimators for the Monte Carlo evaluation of certain functional integrals. _J. Comp. Phys._ , 79, pp. 336–352.