###### Abstract

In this article we consider the application of multilevel Monte Carlo, for the
estimation of normalizing constants. In particular we will make use of the
filtering algorithm, the ensemble Kalman–Bucy filter (EnKBF), which is an
$N$-particle representation of the Kalman–Bucy filter (KBF). The EnKBF is of
interest as it coincides with the optimal filter in the continuous-linear
setting, i.e. the KBF. This motivates our particular setup in the linear
setting. The resulting methodology we will use is the multilevel ensemble
Kalman–Bucy filter (MLEnKBF). We provide an analysis based on deriving
$\mathbb{L}_{q}$-bounds for the normalizing constants using both the single-
level, and the multilevel algorithms, which is largely based on previous work
deriving the MLEnKBF [10]. Our results will be highlighted through numerical
results, where we firstly demonstrate the error-to-cost rates of the MLEnKBFs
comparing it to the EnKBF on a linear Gaussian model. Our analysis will be
specific to one variant of the MLEnKBF, whereas the numerics will be tested on
different variants. We also exploit this methodology for parameter estimation,
where we test this on the models arising in atmospheric sciences, such as the
stochastic Lorenz 63 and 96 model.

Keywords: Multilevel Monte Carlo, Filtering, Kalman–Bucy Filter, Normalizing
Constant,
Parameter estimation
AMS subject classifications: 60G35, 62F15, 65C05, 62M20

Multilevel Estimation of Normalization Constants
using Ensemble Kalman–Bucy Filters

HAMZA RUZAYQAT, NEIL K. CHADA & AJAY JASRA

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

Filtering [2, 13, 16] is the mathematical discipline concerned with the
conditional probability of an unobserved latent process, given sequentially
observed data. It can be found in a wide array of applications, most notably;
numerical weather prediction, mathematical finance, geophysical sciences and
more recently machine learning [4, 32, 33]. Mathematically, given a
$d_{x}$-dimensional unobserved signal process $\\{X_{t}\\}_{t\geq 0}$, and a
$d_{y}$-dimensional $\\{Y_{t}\\}_{t\geq 0}$ observed process, defined as

$\displaystyle dY_{t}$ $\displaystyle=h(X_{t})dt+dV_{t},$ (1.1) $\displaystyle
dX_{t}$ $\displaystyle=f(X_{t})dt+\sigma(X_{t})dW_{t},$ (1.2)

the aim of the filtering problem is to compute the following expectation
$\mathbb{E}[\varphi(X_{t})|\mathscr{F}_{t}]$, where
$\varphi:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}$ is an appropriately
integrable function and $\\{\mathscr{F}_{t}\\}_{t\geq 0}$ is the filtration
generated by the observed process (1.1). From (1.1) - (1.2), $V_{t}$ and
$W_{t}$ are independent $d_{y}$ and $d_{x}-$dimensional Brownian motions
respectively, with $h:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{y}}$ and
$f:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$ denoting potentially
nonlinear functions, and
$\sigma:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}\times d_{x}}$ acting as
a diffusion coefficient. Aside from computing the filtering distribution,
filtering can also be exploited to compute normalizing constants associated
with the filtering distribution [7, 8, 22, 28, 39], i.e. the marginal
likelihood, which is an important and useful computation in Bayesian
statistics. It is useful as it can be used for model comparison which is
commonly done through Bayes’ factor. These quantities are particularly useful
in mixture models such as time-series state space models, and hierarchical
models. Our motivation from this work is the estimation of the normalizing
constants, associated with the filtering distribution in the continuous
linear-Gaussian setting. It is well-known in this setting, that the optimal
filter is the Kalman–Bucy filter (KBF) [27]. However, the use of such a filter
can be challenging to work with, such as firstly actually simulating such
processes, or secondly the associated computational cost. As a result an
alternative to this filter is the ensemble Kalman–Bucy filter (EnKBF), which
is an $N$-particle representation of the KBF, or also as approximations of
conditional McKean–Vlasov-type diffusion processes. This filter has been
analyzed extensively of recent through various pieces of work which include,
but not limited to, understanding stability, deriving uniform propagation of
chaos bounds [5, 6, 17] and acquiring a multilevel estimator [10, 20, 21].
However in the context of estimating normalizing constants, recent work has
been done by Crisan et al. [12] where the authors provide a way to do so using
the EnKBF. However despite this, there is a computational burden associated
with such costs, for example with the EnKBF or other Monte-Carlo algorithms,
to attain a mean squared error (MSE) of a particular order. Therefore it would
be of interest to apply techniques, to reduce the cost associated of attaining
a particular MSE. This is the incentive behind the group of methods based on
multilevel Monte Carlo (MLMC).

MLMC is a numerical technique for stochastic computation, concerned with
reducing the computational complexity of Monte Carlo. Specifically it reduces
the cost of attaining an MSE of an order $\mathcal{O}(\epsilon^{2})$, for
$\epsilon>0$. The ideas of MLMC go back to the original work of Giles and
Hendrick [20, 21], which uses mesh refinements, and a telescoping sum
property, to reduce the cost. Since then it has been applied to numerous
applications, and disciplines, with filtering being one of them. Specifically
for Kalman-based filtering this has been applied to the discrete-EnKF, with
particular variants also analyzed, and applications arising in reservoir
modeling [11, 19, 23, 24]. However more recently there has been the extension
to the continuous-time setting, for the EnKBF [10], entitled MLEnKBF.

In this article we are interested in developing multilevel estimators related
to Kalman filtering, for the computation of normalizing constant. In
particular we will propose the use of various MLEnKBFs for the application of
normalizing constant estimation. Our motivation for this, as mentioned, is
that firstly applying ML techniques can reduce the cost, associated to
attaining an MSE of a particular oder, compared to the normalizing constant
estimator in [12]. Due to this we will make use of the methodology and proof
arguments discussed in [10]. Secondly as it coincides with KBF, it provides an
incentive to work in the linear setting, which can be more cost-effective in
terms of the cost of the algorithm, to other filtering techniques, which make
use of the MLMC, such as the multilevel particle filter and sequential Monte
Carlo sampler [14, 25, 26]. We emphasize with this paper, that we are not
focused on such a comparison.

###  1.1 Contributions

Our contributions of this manuscript are highlighted through the following
points:

  * •

We introduce and approximate a multilevel estimator for the computation of
normalizing constants associated with EnKBFs. This formulation is based on the
multilevel EnKBF introduced by Chada et al. [10], which looks to extend the
single-level estimator which was proposed by Crisan et al. [12].

  * •

Through our formulation we provide firstly a propagation of chaos result for
the normalizing constants associated with the EnKBF, which requires
appropriate $\mathbb{L}_{q}-$bounds, for $q\in[1,\infty)$. To achieve a MSE of
order $\mathcal{O}(\epsilon^{2})$, we require a cost of
$\mathcal{O}(\epsilon^{-3})$, for $\epsilon>0$. This is then extended to the
multilevel setting, where we prove a similar result for an ‘ideal’, or i.i.d.,
system of the MLEnKBF. In particular to achieve a MSE of order
$\mathcal{O}(\epsilon^{2})$, we require a cost of
$\mathcal{O}(\epsilon^{-2}|\log(\epsilon)|)$, for $\epsilon>0$. The analysis
is specific to the vanilla variant, which will form our main result of this
work and naturally follows from the analysis conducted in [10].

  * •

We verify the analysis derived in the paper for various numerical experiments.
We firstly demonstrate the rates attained on an Ornstein–Uhlenbeck process,
and provide parameter estimation on the stochastic Lorenz 63 and Lorenz 96
models. Furthermore we test other variants of the EnKBF, where the analysis
does not directly apply, such as the deterministic and deterministic-transport
variants, to see how they perform.

###  1.2 Outline

The outline of this paper is presented as follows. In Section 2 we review and
discuss the KBF which will motivate the introduction of the EnKBF and the
three variants we consider. This will lead onto Section 3 where we discuss
normalizing constants associated to the filtering distribution, and how one
can use the EnKBF and the MLEnKBF. Furthermore in this section we provide our
main result and demonstrate it on a toy linear example. In Section 4 we
provide an algorithm for parameter estimation using the MLEnKBF tested on
various models arising in atmospheric sciences. We conclude with some remarks
and suggest future directions of work in the Section 5. Finally the proofs of
our results are presented in the Appendix.

##  2 Model and Background

In this section we provide an overview of the necessary background material
required for the rest of the article. We begin by providing the continuous-
time filtering problem through the Kalman–Bucy filters (KBF). This will lead
onto a discussion of ensemble Kalman–Bucy filters (EnKBFs), which are
$N-$particle representation of the KBF. Finally we will discuss the concept of
multilevel Monte Carlo (MLMC), and review its extension to the EnKBF, referred
to as the MLEnKBF.

###  2.1 Kalman–Bucy Filters

Consider a linear-Gaussian filtering model of the following form

$\displaystyle dY_{t}$ $\displaystyle=CX_{t}dt+R^{1/2}dV_{t},$ (2.1)
$\displaystyle dX_{t}$ $\displaystyle=AX_{t}dt+Q^{1/2}dW_{t},$ (2.2)

where $(Y_{t},X_{t})\in\mathbb{R}^{d_{y}}\times\mathbb{R}^{d_{x}}$,
$(V_{t},W_{t})$ is a $(d_{y}+d_{x})-$dimensional standard Brownian motion, $A$
is a square $d_{x}\times d_{x}$ matrix, $C$ is a $d_{y}\times d_{x}$ matrix,
$Y_{0}=0$, $X_{0}\sim\mathcal{N}_{d_{x}}(\mathcal{M}_{0},\mathcal{P}_{0})$
($d_{x}-$dimensional Gaussian distribution, mean $\mathcal{M}_{0}$, covariance
matrix $\mathcal{P}_{0}$) and $R^{1/2},Q^{1/2}$ are square (of the appropriate
dimension) and symmetric and invertible matrices. It is well-known that,
letting $\\{\mathscr{F}_{t}\\}_{t\geq 0}$ be the filtration generated by the
observations, the conditional probability of $X_{t}$ given $\mathscr{F}_{t}$
is a Gaussian distribution with mean and covariance matrix

$\mathcal{M}_{t}:=\mathbb{E}[X_{t}|\mathscr{F}_{t}],\quad\mathcal{P}_{t}:=\mathbb{E}\Big{(}[X_{t}-\mathbb{E}(X_{t}|\mathscr{F}_{t})][X_{t}-\mathbb{E}(X_{t}|\mathscr{F}_{t})]^{\top}\Big{)},$

given by the Kalman–Bucy and Ricatti equations [17]

$\displaystyle d\mathcal{M}_{t}$
$\displaystyle=A\mathcal{M}_{t}dt+\mathcal{P}_{t}C^{\top}R^{-1}\Big{(}dY_{t}-C\mathcal{M}_{t}dt\Big{)},$
(2.3) $\displaystyle\partial_{t}\mathcal{P}_{t}$
$\displaystyle=\textrm{Ricc}(\mathcal{P}_{t}),$ (2.4)

where the Riccati drift term is defined as

$\textrm{Ricc}(G)=AG+GA^{\top}-GSG+Q,\quad\textrm{with}\ S:=C^{\top}R^{-1}C.$

A derivation of (2.3) - (2.4) can be found in [27]. KBF is viewed as the
$\mathbb{L}_{2}$-optimal state estimator for an Ornstein–Uhlenbeck process,
given the state is partially observed with linear and Gaussian assumptions. An
alternative approach is to consider a conditional McKean-Vlasov type diffusion
process. For this article we work with three different processes of the form

$\displaystyle d\overline{X}_{t}$ $\displaystyle=$ $\displaystyle
A\leavevmode\nobreak\ \overline{X}_{t}\leavevmode\nobreak\
dt\leavevmode\nobreak\ +\leavevmode\nobreak\ Q^{1/2}\leavevmode\nobreak\
d\overline{W}_{t}+{\cal P}_{\eta_{t}}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(C\overline{X}_{t}dt+R^{1/2}\leavevmode\nobreak\
d\overline{V}_{t}\right)\right],$ (2.5) $\displaystyle d\overline{X}_{t}$
$\displaystyle=$ $\displaystyle A\leavevmode\nobreak\
\overline{X}_{t}\leavevmode\nobreak\ dt\leavevmode\nobreak\
+\leavevmode\nobreak\ Q^{1/2}\leavevmode\nobreak\ d\overline{W}_{t}+{\cal
P}_{\eta_{t}}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(\frac{1}{2}C\left[\overline{X}_{t}+\eta_{t}(e)\right]dt\right)\right],$
(2.6) $\displaystyle d\overline{X}_{t}$ $\displaystyle=$ $\displaystyle
A\leavevmode\nobreak\ \overline{X}_{t}\leavevmode\nobreak\
dt\leavevmode\nobreak\ +\leavevmode\nobreak\ Q{\cal
P}_{\eta_{t}}^{-1}\left(\overline{X}_{t}-\eta_{t}(e)\right)\leavevmode\nobreak\
dt+{\cal P}_{\eta_{t}}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(\frac{1}{2}C\left[\overline{X}_{t}+\eta_{t}(e)\right]dt\right)\right],$
(2.7)

where $(\overline{V}_{t},\overline{W}_{t},\overline{X}_{0})$ are copies of the
process of $(V_{t},W_{t},X_{0})$ and covariance

$\mathcal{P}_{\eta_{t}}=\eta_{t}\Big{(}[e-\eta_{t}(e)][e-\eta_{t}(e)]^{\top}\Big{)},\quad\eta_{t}:=\mathrm{Law}(\overline{X}_{t}|\mathscr{F}_{t}),$

such that $\eta_{t}$ is the conditional law of $\overline{X}_{t}$ given
$\mathscr{F}_{t}$ and $e(x)=x$. We will explain the difference of each
diffusion process, in succeeding subsections. It is important to note that the
nonlinearity in (2.5)-(2.7) does not depend on the distribution of the state
$\mathrm{Law}(\overline{X}_{t})$ but on the conditional distribution
$\eta_{t}$, and $\mathcal{P}_{\eta_{t}}$ alone does not depend on
$\mathscr{F}_{t}$. These processes are commonly referred to as Kalman-Bucy
(nonlinear) diffusion processes. It is known that the conditional expectations
of the random states $\overline{X}_{t}$ and their conditional covariance
matrices $\mathcal{P}_{\eta_{t}}$, with respect to $\mathscr{F}_{t}$, satisfy
the Kalman–Bucy and the Riccati equations. In addition, for any
$t\in\mathbb{R}^{+}$

$\eta_{t}:=\mathrm{Law}(\overline{X}_{t}|\mathscr{F}_{t})=\mathrm{Law}({X}_{t}|\mathscr{F}_{t}).$

As a result, an alternative to recursively computing (2.3) - (2.4), is to
generate $N$ i.i.d. samples from any of (2.5)-(2.7) processes and apply a
Monte Carlo approximation, as mentioned which we now discuss.

###  2.2 Ensemble Kalman–Bucy Filters

Exact simulation from (2.5)-(2.7) is typically not possible, or feasible, as
one cannot compute $\mathcal{P}_{\eta_{t}}$ exactly. The ensemble Kalman–Bucy
filter (EnKBF) can be used to deal with this issue. The EnKBF coincides with
the mean-field particle interpretation of the Kalman-Bucy diffusion processes.
EnKBFs is an $N-$particle system simulated for the $i^{th}-$particle,
$i\in\\{1,\dots,N\\}$. The first variant of the EnKBF we consider is given as

$d\xi_{t}^{i}=A\leavevmode\nobreak\ \xi_{t}^{i}\leavevmode\nobreak\
dt\leavevmode\nobreak\ +\leavevmode\nobreak\ Q^{1/2}\leavevmode\nobreak\
d\overline{W}^{i}_{t}+P^{N}_{t}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(C\xi_{t}^{i}dt+R^{1/2}\leavevmode\nobreak\
d\overline{V}^{i}_{t}\right)\right],$ (2.8)

which is known as the vanilla EnKBF (VEnKBF). This is the standard EnKBF used
in theory and practice which contains perturbed observations, through the
brownian motion $\overline{V}_{t}$. It is a continuous-time derivation of the
EnKF [18]. The second EnKBF given, defined as,

$d\xi_{t}^{i}=A\leavevmode\nobreak\ \xi_{t}^{i}\leavevmode\nobreak\
dt\leavevmode\nobreak\ +\leavevmode\nobreak\ Q^{1/2}\leavevmode\nobreak\
d\overline{W}^{i}_{t}+P^{N}_{t}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(\frac{1}{2}C\left[\xi_{t}^{i}+\mathcal{M}_{t}\right]dt\right)\right],$
(2.9)

which is referred to as the deterministic EnKBF (DEnKBF), where unlike (2.8),
it contains no perturbed observations. This is the continuous-limiting object
of the deterministic EnKF defined in [36]. It is well-known in ensemble data
assimilation, that deterministic filters can perform better, partially due to
the containing exact observations with no noise. Our final variant of the
EnKBF is the deterministic transport EnKBF (DTEnKBF),

$d\xi_{t}^{i}=A\leavevmode\nobreak\ \xi_{t}^{i}\leavevmode\nobreak\
dt\leavevmode\nobreak\ +\leavevmode\nobreak\
QP_{t}^{-1}\left(\xi_{t}^{i}-\mathcal{M}_{t}\right)\leavevmode\nobreak\
dt+P^{N}_{t}C^{\top}R^{-1}\leavevmode\nobreak\
\left[dY_{t}-\left(\frac{1}{2}C\left[\xi_{t}^{i}+\mathcal{M}_{t}\right]dt\right)\right],$
(2.10)

where the modification is that it does not contain $\overline{V}_{t}$ and
$\overline{W}_{t}$, implying it is completely deterministic. This filter, is
motivated from the use of optimal transport methodologies within data
assimilation [34]. Such processes are discussed in further detail in [5]. For
all the variants of the EnKBF the sample mean and covariances are defined as

$\displaystyle{P}_{t}^{N}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}(\xi_{t}^{i}-m^{N}_{t})(\xi_{t}^{i}-m^{N}_{t})^{\top},$
$\displaystyle m^{N}_{t}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\xi_{t}^{i},$

where
$\xi_{0}^{i}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d_{x}}(\mathcal{M}_{0},\mathcal{P}_{0})$.
Note that when $C=0$, (2.8) and (2.9) reduce to $N$–independent copies of an
Ornstein–Uhlenbeck process.

In practice, one will not have access to an entire trajectory of observations.
Thus numerically, one often works with a time discretization, such as an
Euler-type discretization. Let $\Delta_{l}=2^{-l}$ denote our level of
discretization, then we will generate the system for
$(i,k)\in\\{1,\dots,N\\}\times\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$ as

######  Remark 2.1.

As mentioned in the introduction, the function $f$ in (1.2) can be nonlinear.
If we assume in addition that $h$ is the identity function, then one can
modify equations (2.8)– (2.10) by replacing the term $A\leavevmode\nobreak\
\overline{X}_{t}$ with $f(\overline{X}_{t})$ and equations (2.11)– (2.13) by
replacing the term $A\leavevmode\nobreak\ \xi_{k\Delta_{l}}^{i}$ with
$f(\xi_{k\Delta_{l}}^{i})$.

######  Remark 2.2.

When mentioning perturbed observations, we are referring to the additional
noise in the innovation term. However, we emphasize it is not the observation
being perturbed, but rather the projected state. hence once can think of this
as a stochastic filter.

$\displaystyle\xi_{(k+1)\Delta_{l}}^{i}$ $\displaystyle=$
$\displaystyle\xi_{k\Delta_{l}}^{i}+A\xi_{k\Delta_{l}}^{i}\Delta_{l}+Q^{1/2}\big{\\{}\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}\big{\\}}+$
(2.11) $\displaystyle
P^{N}_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}\big{\\{}Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}\big{\\}}-\Big{[}C\xi_{k\Delta_{l}}^{i}\Delta_{l}+R^{1/2}\big{\\{}\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}\big{\\}}\Big{]}\Big{)},$
$\displaystyle\xi_{(k+1)\Delta_{l}}^{i}$ $\displaystyle=$
$\displaystyle\xi_{k\Delta_{l}}^{i}+A\xi_{k\Delta_{l}}^{i}\Delta_{l}+Q^{1/2}\big{\\{}\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}\big{\\}}+$
$\displaystyle
P^{N}_{k\Delta_{l}}C^{\top}R^{-1}\left(\big{\\{}Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}\big{\\}}-C\left(\dfrac{\xi_{k\Delta_{l}}^{i}+m^{N}_{k\Delta_{l}}}{2}\right)\Delta_{l}\right),$
(2.12) $\displaystyle\xi_{(k+1)\Delta_{l}}^{i}$ $\displaystyle=$
$\displaystyle\xi_{k\Delta_{l}}^{i}+A\xi_{k\Delta_{l}}^{i}\Delta_{l}+Q\left(P_{k\Delta_{l}}\right)^{-1}\left(\xi_{k\Delta_{l}}^{i}-m^{N}_{k\Delta_{l}}\right)\Delta_{l}+$
(2.13) $\displaystyle
P^{N}_{k\Delta_{l}}C^{\top}R^{-1}\left(\big{\\{}Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}\big{\\}}-C\left(\dfrac{\xi_{k\Delta_{l}}^{i}+m^{N}_{k\Delta_{l}}}{2}\right)\Delta_{l}\right),$

such that

$\displaystyle P_{k\Delta_{l}}^{N}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}(\xi_{k\Delta_{l}}^{i}-m^{N}_{k\Delta_{l}})(\xi_{k\Delta_{l}}^{i}-m_{k\Delta_{l}}^{N})^{\top},$
$\displaystyle m_{k\Delta_{l}}^{N}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\xi_{k\Delta_{l}}^{i},$

and
$\xi_{0}^{i}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d_{x}}(\mathcal{M}_{0},\mathcal{P}_{0})$.
For $l\in\mathbb{N}_{0}$ given, denote by $\eta_{t}^{N,l}$ as the
$N-$empirical probability measure of the particles
$(\xi_{t}^{1},\dots,\xi_{t}^{N})$, where
$t\in\\{0,\Delta_{l},2\Delta_{l},\dots\\}$. For
$\varphi:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$ we will use the
notation
$\eta_{t}^{N,l}(\varphi):=\tfrac{1}{N}\sum_{i=1}^{N}\varphi(\xi_{t}^{i})$.

######  Remark 2.3.

For ensemble based data assimilation, a common recurring assumption is that
the particles are i.i.d.. While this is true, through ensemble filtering
methodologies the assumption is practically false. This is based on the fact
that ensemble covariance is computed from all ensemble members, which
introduces a dependence, and in some cases the EnKF can be applied to
nonlinear functions, ensuring the ensemble is not normally distributed.

###  2.3 Multilevel EnKBFs

Let us define $\pi$ to be be a probability on a measurable space
$(\mathsf{X},\mathscr{X})$ and for $\pi-$integrable
$\varphi:\mathsf{X}\rightarrow\mathbb{R}$ consider the problem of estimating
$\pi(\varphi)=\mathbb{E}_{\pi}[\varphi(X)]$. Now let us assume that we only
have access to a sequence of approximations of $\pi$,
$\\{\pi_{l}\\}_{l\in\mathbb{N}_{0}}$, also each defined on
$(\mathsf{X},\mathscr{X})$ and we are now interested in estimating
$\pi_{l}(\varphi)$, such that
$\lim_{l\rightarrow\infty}|[\pi_{l}-\pi](\varphi)|=0$. Therefore one can use
the telescoping sum

$\pi_{L}(\varphi)=\pi_{0}(\varphi)+\sum^{L}_{l=1}[\pi_{l}-\pi_{l-1}](\varphi),$
(2.14)

as we know that the approximation error between $\pi$ and $\pi_{l}$ gets
smaller as $l\rightarrow\infty$. The idea of multilevel Monte Carlo (MLMC)
[20, 21] is to construct a coupled system, related to the telescoping sum,
such that the mean squared error can be reduced, relative to i.i.d. sampling
from $\pi_{L}$. Therefore our MLMC approximation of
$\mathbb{E}_{\pi_{L}}[\varphi(X)]$ is

$\pi_{L}^{ML}(\varphi):=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\varphi(X^{i,0})+\sum_{l=1}^{L}\frac{1}{N_{l}}\sum_{i=1}^{N_{l}}\\{\varphi(X^{i,l})-\varphi(\tilde{X}^{i,l-1})\\},$

where $N_{0}\in\mathbb{N}$ i.i.d. samples from $\pi_{0}$ as
$(X^{1,0},\dots,X^{N_{0},0})$ and for $l\in\\{1,\dots,L\\}$,
$N_{l}\in\mathbb{N}$ samples from a coupling of $(\pi_{l},\pi_{l-1})$ as
$((X^{1,l},\tilde{X}^{1,l-1}),\dots,(X^{N_{l},l},\tilde{X}^{N_{l},l-1}))$.
Then the MSE is then

$\mathbb{E}[(\pi_{L}^{ML}(\varphi)-\pi(\varphi))^{2}]=\underbrace{\mathbb{V}\textrm{ar}[\pi_{L}^{ML}(\varphi)]}_{\textrm{variance}}+[\underbrace{\pi_{L}(\varphi)-\pi(\varphi)}_{\textrm{bias}}]^{2}.$

It is important to note that the telescoping sum (2.14), despite its
simplicity and usefulness, is not an optimal choice to construct the
multilevel estimator. In particular Giles [21] makes this assumption, but
states that a potentially more robust weighted sum would be more optimal. The
work of Chada et al. [10] proposed the application of MLMC for the VEnKBF,
where we will now briefly review this. If we consider the discretized VEnKBF
(2.11), then, for $l\in\mathbb{N}$ and
$(i,k)\in\\{1,\dots,N\\}\times\mathbb{N}_{0}$ the ML adaption is given as

$\textbf{(F1)}\begin{cases}\begin{aligned}
\xi_{(k+1)\Delta_{l}}^{i,l}&=\xi_{k\Delta_{l}}^{i,l}+A\xi_{k\Delta_{l}}^{i,l}\Delta_{l}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]\\\
&+P_{k\Delta_{l}}^{N,l}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l}}-Y^{{i}}_{k\Delta_{l}}]-\Big{[}C\xi_{k\Delta_{l}}^{i,l}\Delta_{l}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}]\Big{]}\Big{)},\\\
\xi_{(k+1)\Delta_{l-1}}^{i,l-1}&=\xi_{k\Delta_{l-1}}^{i,l-1}+A\xi_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l-1}}^{i}-\overline{W}_{k\Delta_{l-1}}^{i}]\\\
&+P_{k\Delta_{l-1}}^{N,l-1}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l-1}}-Y^{{i}}_{k\Delta_{l-1}}]-\Big{[}C\xi_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l-1}}^{i}-\overline{V}_{k\Delta_{l-1}}^{i}]\Big{]}\Big{)}.\end{aligned}\end{cases}$

Similarly for the deterministic variant (2.12) we have

$\textbf{(F2)}\begin{cases}\begin{aligned}
\xi_{(k+1)\Delta_{l}}^{i,l}&=\xi_{k\Delta_{l}}^{i,l}+A\xi_{k\Delta_{l}}^{i,l}\Delta_{l}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]\\\
&+P_{k\Delta_{l}}^{N,l}C^{\top}R^{-1}\left([Y^{{i}}_{(k+1)\Delta_{l}}-Y^{{i}}_{k\Delta_{l}}]-C\left(\dfrac{\xi_{k\Delta_{l}}^{i,1}+m^{N,l}_{k\Delta_{l}}}{2}\right)\Delta_{l}\right),\\\
\xi_{(k+1)\Delta_{l-1}}^{i,l-1}&=\xi_{k\Delta_{l-1}}^{i,l-1}+A\xi_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l-1}}^{i}-\overline{W}_{k\Delta_{l-1}}^{i}]\\\
&+P_{k\Delta_{l-1}}^{N,l-1}C^{\top}R^{-1}\left([Y^{{i}}_{(k+1)\Delta_{l-1}}-Y^{{i}}_{k\Delta_{l-1}}]-C\left(\dfrac{\xi_{k\Delta_{l-1}}^{i,l-1}+m^{N,l-1}_{k\Delta_{l-1}}}{2}\right)\Delta_{l-1}\right),\end{aligned}\end{cases}$

and finally for the deterministic-transport variant (2.13)

$\textbf{(F3)}\begin{cases}\begin{aligned}
\xi_{(k+1)\Delta_{l}}^{i,l}&=\xi_{k\Delta_{l}}^{i,l}+A\xi_{k\Delta_{l}}^{i,l}\Delta_{l}+Q\left(P_{k\Delta_{l}}\right)^{-1}[\xi_{k\Delta_{l}}^{i,l}-m_{k\Delta_{l}}^{N,l}]\Delta_{l}\\\
&+P_{k\Delta_{l}}^{N,l}C^{\top}R^{-1}\left([Y^{{i}}_{(k+1)\Delta_{l}}-Y^{{i}}_{k\Delta_{l}}]-C\left(\dfrac{\xi_{k\Delta_{l}}^{i,l}+m^{N,l}_{k\Delta_{l}}}{2}\right)\Delta_{l}\right),\\\
\xi_{(k+1)\Delta_{l-1}}^{i,l-1}&=\xi_{k\Delta_{l-1}}^{i,l-1}+A\xi_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+Q\left(P_{k\Delta_{l-1}}\right)^{-1}[\xi_{k\Delta_{l-1}}^{i,l-1}-m_{k\Delta_{l-1}}^{N,l-1}]\Delta_{l-1}\\\
&+P_{k\Delta_{l-1}}^{N,l-1}C^{\top}R^{-1}\left([Y^{{i}}_{(k+1)\Delta_{l-1}}-Y^{{i}}_{k\Delta_{l-1}}]-C\left(\dfrac{\xi_{k\Delta_{l-1}}^{i-1}+m_{k\Delta_{l-1}}}{2}\right)\Delta_{l-1}\right),\end{aligned}\end{cases}$

where our sample covariances and means are defined accordingly as

$\displaystyle P_{k\Delta_{l}}^{N,l}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}(\xi_{k\Delta_{l}}^{i,l}-m_{k\Delta_{l}}^{N,l})(\xi_{k\Delta_{l}}^{i,l}-m_{k\Delta_{l}}^{N,l})^{\top},$
$\displaystyle m_{k\Delta_{l}}^{N,l}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\xi_{k\Delta_{l}}^{i,l},$
$\displaystyle P_{k\Delta_{l-1}}^{N,l-1}$
$\displaystyle=\frac{1}{N-1}\sum_{i=1}^{N}(\xi_{k\Delta_{l-1}}^{i,l-1}-m_{k\Delta_{l-1}}^{N,l-1})(\xi_{k\Delta_{l-1}}^{i,l-1}-m_{k\Delta_{l-1}}^{N,l-1})^{\top},$
$\displaystyle m_{k\Delta_{l-1}}^{N,l-1}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\xi_{k\Delta_{l-1}}^{i,l-1},$

and
$\xi_{0}^{i,l}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d_{x}}(\mathcal{M}_{0},\mathcal{P}_{0})$,
$\xi_{0}^{i,l-1}=\xi_{0}^{i,l}$. Then, one has the approximation of
$[\eta_{t}^{l}-\eta_{t}^{l-1}](\varphi)$, $t\in\mathbb{N}_{0}$,
$\varphi:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}$, given as

$[\eta_{t}^{N,l}-\eta_{t}^{N,l-1}](\varphi)=\frac{1}{N}\sum_{i=1}^{N}[\varphi(\xi_{t}^{i,l})-\varphi(\xi_{t}^{i,l-1})].$

Therefore for the multilevel estimation, one has an approximation for
$t\in\mathbb{N}_{0}$

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

Similarly with Remark 2.3, the i.i.d. assumption regarding each particle does
not also hold in the multilevel setting. Largely due to the correlation of the
particle in the multilevel setting.

######  Remark 2.4.

For the above forms of the MLEnKBF, only (F1) - (F2) were introduced and
numerically tested in [10]. Therefore this article aims to further test the
deterministic-transport variant of (F3).

##  3 Normalizing Constant Estimation

In this section we introduce the notion of normalizing constants (NC), which
is what we are concerned with estimating. We begin by recalling the NC
estimator using the EnKBF, aiming to show its cost associated to attain an MSE
of a particular order. We then present our NC estimator through the MLEnKBF.
This will lead onto our main result of the paper, which is the error-to-cost
ratio of our normalizing constant estimator, i.e. the MLEnKBF, which is
presented through Theorem 3.1.

###  3.1 Normalizing Constants

We begin by defining the normalizing constant associated with the filtering
distribution. This is defined as
$Z_{t}:=\frac{{\mathcal{L}}_{X_{0:t},Y_{0:t}}}{{\mathcal{L}}_{X_{0:t},W_{0:t}}}$
to be the density of ${\mathcal{L}}_{X_{0:t},Y_{0:t}}$, the law of the process
$(X,Y)$ and that of ${\mathcal{L}}_{X_{0:t},W_{0:t}}$, the law of the process
$(X,W)$. That is,

${\mathbb{E}}[f(X_{0:t})g(Y_{0:t})]={\mathbb{E}}[f(X_{0:t})g(W_{0:t})Z_{t}(X,Y)].$

One can show that (see Exercise 3.14 pp 56 in [2])

$Z_{t}(X,Y)=\exp{\left[\int_{0}^{t}\left[\langle
CX_{s},R^{-1}dY_{s}\rangle-\frac{1}{2}\langle
X_{s},SX_{s}\rangle\leavevmode\nobreak\ ds\right]\right]},$ (3.1)

recalling that $S:=C^{\top}R^{-1}C$. Now we let $\overline{Z}_{t}(Y)$ denote
the likelihood function defined by

$\overline{Z}_{t}(Y):=\mathbb{E}_{Y}\left(Z_{t}(X,Y)\right),$

where $\mathbb{E}_{Y}\left(\cdot\right)$ stands for the expectation w.r.t. the
signal process when the observation is fixed and independent of the signal.
From the work of Crisan et al. [12], the authors show that the normalizing
constant is given by

$\displaystyle\overline{Z}_{t}(Y)=\exp{\left[\int_{0}^{t}\left[\leavevmode\nobreak\
\langle
C\mathcal{M}_{s},R^{-1}dY_{s}\rangle-\frac{1}{2}\langle\mathcal{M}_{s},S\mathcal{M}_{s}\rangle\leavevmode\nobreak\
ds\right]\right]},$ (3.2)

and a sensible estimator of it is given as

$\displaystyle\overline{Z}_{t}^{N}(Y)=\exp{\left[\int_{0}^{t}\left[\leavevmode\nobreak\
\langle Cm_{s}^{N},R^{-1}dY_{s}\rangle-\frac{1}{2}\langle
m_{s}^{N},Sm_{s}^{N}\rangle\leavevmode\nobreak\ ds\right]\right]},$ (3.3)

which follows from replacing the conditional mean of the signal process with
the sample mean associated with the EnKBF. To briefly describe how (3.2) was
attained in [12], the authors used a change of measure rule, through
Girsanov’s Theorem, then an application of Bayes’ Theorem which can be
expressed through as the Kallianpur-Striebel formula defined on path space.
The derivation is a standard approach, which is discussed in Chapter 3 in [2].

In practice, one must time-discretize the EnKBF, and the normalizing constant
estimator (3.3) to yield for $t\in\mathbb{N}$

$\displaystyle\overline{Z}_{t}^{N,l}(Y)=\exp\Bigg{\\{}\sum_{k=0}^{t\Delta_{l}^{-1}-1}\left[\leavevmode\nobreak\
\langle
Cm_{k\Delta_{l}}^{N},R^{-1}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]\rangle-\frac{\Delta_{l}}{2}\langle
m_{k\Delta_{l}}^{N},Sm_{k\Delta_{l}}^{N}\rangle\right]\Bigg{\\}}.$ (3.4)

Let $\overline{U}_{t}^{N,l}(Y)=\log(\overline{Z}_{t}^{N,l}(Y))$, where we now
consider the estimation of log-normalization constants. To enhance the
efficiency we consider a coupled ensemble Kalman–Bucy filter, as described in
Section 2. Let $l\in\mathbb{N}$ then we run the coupled system of the
different MLEnKBFs. Our multilevel estimator of log-normalizing constant given
in Algorithm 1.

######  Remark 3.1.

In order to use the analysis derived in [10], for the normalizing constant
estimation, our results will be specific to the vanilla MLEnKBF variant (F1).
This is important as we will use the results presented in [10], however, in
terms of the numerics, we will use all ML variants described for various
parameter estimation experiments.

###  3.2 Single-level EnKBF

To present a multilevel NC estimator using the EnKBF, we first require to
understand the single-level EnKBF. To aid our analysis we will consider the
i.i.d. particle system based upon the Euler discretization of EnKBF for (F1),
such that for $(i,k)\in\\{1,\dots,N\\}\times\mathbb{N}_{0}$

$\displaystyle\zeta_{(k+1)\Delta_{l}}^{i}$
$\displaystyle=(I+A\Delta_{l})\zeta_{k\Delta_{l}}^{i}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]$
(3.5)
$\displaystyle+P_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]-\Big{[}C\zeta_{k\Delta_{l}}^{i}\Delta_{l}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}]\Big{]}\Big{)},$

such that
$\zeta_{(k+1)\Delta_{l}}^{i}|\mathscr{F}_{(k+1)\Delta_{l}}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\mathcal{N}_{d_{x}}(m_{(k+1)\Delta_{l}},P_{(k+1)\Delta_{l}})$,
where the moments are defined as

$\displaystyle m_{(k+1)\Delta_{l}}$
$\displaystyle=m_{k\Delta_{l}}+Am_{k\Delta_{l}}\Delta_{l}+P_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]-Cm_{k\Delta_{l}}\Delta_{l}\Big{)},$
(3.6) $\displaystyle P_{(k+1)\Delta_{l}}$
$\displaystyle=P_{k\Delta_{l}}+\textrm{Ricc}(P_{k\Delta_{l}})\Delta_{l}+(A-P_{k\Delta_{l}}S)P_{k\Delta_{l}}(A^{\top}-SP_{k\Delta_{l}})\Delta_{l}^{2},$
(3.7)

which are satisfied by the Kalman–Bucy diffusion (2.5). We now consider the NC
estimator $\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}$, which instead
using the recursion of (3.5).

######  Remark 3.2.

Our reason for introducing the new estimator, based on (3.5) - (3.7), is that
the usual described estimator of $\overline{U}_{T}^{ML}(Y)$ will not hold, in
terms of the analysis, due to the recursion of the EnKBF, which is mentioned
and documented in [10]. However, as we will see, by using a number of
arguments one can show that both the multilevel estimators are close in
probability, as $N\rightarrow\infty$. Therefore our analysis will be specific
to ML NC estimator $\widehat{\overline{U}}_{T}^{ML}(Y)$ i.e.

$\widehat{\overline{U}}_{T}^{ML}(Y)=\widehat{\overline{U}}_{T}^{N_{l_{*}},l_{*}}(Y)+\sum_{l=l_{*}+1}^{L}\\{\widehat{\overline{U}}_{T}^{N_{l},l}(Y)-\widehat{\overline{U}}_{T}^{N_{l},l-1}(Y)\\},$

while we continue to use $\overline{U}_{T}^{ML}(Y)$ for the numerical
examples.

We now present our first result of the paper which is a propagation of chaos
result for the single-level vanilla EnKBF NC estimator. For a
$d_{x}-$dimensional vector $x$ we denote
$\|x\|_{2}=(\sum_{j=1}^{d_{x}}x(j)^{2})^{1/2}$, where $x(j)$ is the
$j^{th}-$element of $x$. We use the notation
$[{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}](Y)={\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}(Y)$.

######  Proposition 3.1.

For any
$(l,t,k_{1})\in\mathbb{N}_{0}\times\mathbb{R}^{+}\\{0,1,\dots,\Delta_{l}^{-1}\\}$
almost surely:

$\lim_{N\rightarrow\infty}[{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}](Y)=0.$

###### Proof.

The proof follows by using Proposition A.1 in the appendix and the
Marcinkiewicz–Zygmund inequality for i.i.d. random variables, along with a
standard first Borel–Cantelli lemma argument. ∎

Then using the same arguments, as described in [10], one can show, through the
Markov inequality that

$\mathbb{P}\left(\left|{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}(Y)\right|>\varepsilon\right)\leq\frac{\mathsf{C}}{\varepsilon^{2q}N^{q/2}}.$
(3.8)

for any $\varepsilon>0$ and $q>0$, where $\mathsf{C}$ is a constant that
depends on $(l,q,t,k_{1})$ but not $N$. Related to this in Proposition C.1, in
the appendix, we have shown that:

$\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}](Y)\Big{\|}_{2}^{2}\Big{]}\leq\mathsf{C}\Big{(}\frac{1}{N}+\Delta_{l}^{2}\Big{)},$

where $\mathsf{C}$ does not depend upon $l$ nor $N$. Thus to achieve an MSE of
order $\mathcal{O}(\epsilon^{2})$, for $\epsilon>0$, we require a cost of
$\mathcal{O}(\epsilon^{-3})$, under a suitable choice of $L$ and $N$, using
the single-level EnKBF [20].

###  3.3 Multilevel EnKBF

Now to analyze the multilevel estimator, we will aim to prove results for the
i.i.d. (ideal) coupled system for
$(i,k)\in\\{1,\dots,N\\}\times\mathbb{N}_{0}$:

$\displaystyle\zeta_{(k+1)\Delta_{l}}^{i,l}$
$\displaystyle=\zeta_{k\Delta_{l}}^{i,l}+A\zeta_{k\Delta_{l}}^{i,l}\Delta_{l}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]$
(3.9)
$\displaystyle+P^{N,l}_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l}}-Y^{{i}}_{k\Delta_{l}}]-\Big{[}C\zeta_{k\Delta_{l}}^{i,l}\Delta_{l}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}]\Big{]}\Big{)},$
$\displaystyle\zeta_{(k+1)\Delta_{l-1}}^{i,l-1}$
$\displaystyle=\zeta_{k\Delta_{l-1}}^{i,l-1}+A\zeta_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l-1}}^{i}-\overline{W}_{k\Delta_{l-1}}^{i}]$
(3.10)
$\displaystyle+P^{N,l-1}_{k\Delta_{l-1}}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l-1}}-Y^{{i}}_{k\Delta_{l-1}}]-\Big{[}C\zeta_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l-1}}^{i}-\overline{V}_{k\Delta_{l-1}}^{i}]\Big{]}\Big{)},$

Now we set

$\widehat{\overline{U}}_{t}^{ML}(Y):=\widehat{\overline{U}}_{t}^{N_{0},0}(Y)+\sum_{l-1}^{L}\big{[}\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}\big{]}(Y),$
(3.11)

Therefore using the same argument as in [10], one can establish, almost surely

$\lim_{\min_{l}N_{l}\rightarrow\infty}[{\overline{U}}_{t}^{ML}-\widehat{\overline{U}}_{t}^{ML}](Y)=0,\quad\mathbb{P}\left(\left|[{\overline{U}}_{t}^{ML}-\widehat{\overline{U}}_{t}^{ML}](Y)\right|>\varepsilon\right)\leq\frac{\mathsf{C}}{\varepsilon^{2q}}\left(\sum_{l=0}^{L}\frac{1}{N_{l}^{q/2}}\right),$

where $\mathsf{C}$ is a constant that can depend on $(L,q,t)$ but not
$N_{0:L}$.

We are now in a position to state our main result, concerned with using the
vanilla MLEnKBF, for the i.i.d. system, as a NC estimator.

######  Theorem 3.1.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T-1]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(L,N_{0:L})\in\mathbb{N}\times\\{2,3,\dots\\}^{L+1}$,

$\mathbb{E}\left[\left\|[\widehat{\overline{U}}_{t}^{ML}-{\overline{U}}_{t}](Y)\right\|_{2}^{2}\right]\leq\mathsf{C}\left(\sum_{l=0}^{L}\frac{\Delta_{l}}{N_{l}}+\sum_{l=1}^{L}\sum_{q=1,q\neq
l}^{L}\frac{\Delta_{l}\Delta_{q}}{N_{l}N_{q}}+\Delta_{L}^{2}\right).$

The above theorem translates as, in order to achieve an MSE of order
$\mathcal{O}(\epsilon^{2})$, for $\epsilon>0$, we have a cost of
$\mathcal{O}(\epsilon^{-2}\leavevmode\nobreak\ \log(\epsilon)^{2})$. This
implies a reduction in cost compared to the single-level NC estimator.

######  Remark 3.3.

It is important to note that the rates obtained in Proposition 3.1 and Theorem
3.1 for the normalizing constants, are identical to those derived in [10].
However despite this it should be emphasized as we will see in the appendix,
that the proofs are not trivial, and all of them do not directly follow
exactly.

Algorithm 1 (MLEnKBF-NC) Multilevel Estimation of Normalizing Constants

  1. 1.

Input: Target level $L\in\mathbb{N}$, start level $l_{*}\in\mathbb{N}$ such
that $l_{*}<L$, the number of particles on each level
$\\{N_{l}\\}_{l=l_{*}}^{L}$, the time parameter $T\in\mathbb{N}$ and initial
independent ensembles
$\Big{\\{}\\{\tilde{\xi}_{0}^{i,l_{*}}\\}_{i=1}^{N_{l_{*}}},\cdots,\\{\tilde{\xi}_{0}^{i,L}\\}_{i=1}^{N_{L}}\Big{\\}}$.

  2. 2.

Initialize: Set $l=l_{*}$. For
$(i,k)\in\\{1,\cdots,N_{l}\\}\times\\{0,\cdots,T\Delta_{l}^{-1}-1\\}$, set
$\\{\xi_{0}^{i,l}\\}_{i=1}^{N_{l}}=\\{\tilde{\xi}_{0}^{i,l}\\}_{i=1}^{N_{l}}$
and simulate any of the cases (2.11)-(2.13) to return
$\\{m^{N_{l},l}_{k\Delta_{l}}\\}_{k=0}^{T\Delta_{l}^{-1}-1}$. Then using
(3.4), return $\overline{Z}_{T}^{N_{l},l}(Y)$ or
$\overline{U}_{T}^{N_{l},l}(Y)$.

  3. 3.

Iterate: For $l\in\\{l_{*}+1,\cdots,L\\}$ and
$(i,k)\in\\{1,\cdots,N_{l}\\}\times\\{0,\cdots,T\Delta_{l}^{-1}-1\\}$, set
$\\{\xi_{0}^{i,l-1}\\}_{i=1}^{N_{l}}=\\{\xi_{0}^{i,l}\\}_{i=1}^{N_{l}}=\\{\tilde{\xi}_{0}^{i,l}\\}_{i=1}^{N_{l}}$,
and simulate any of the coupled ensembles (F1)-(F3) (the one which corresponds
to the case used in Step 2.) to return
$\\{m^{N_{l},l-1}_{k\Delta_{l-1}}\\}_{k=0}^{T\Delta_{l-1}^{-1}-1}$ and
$\\{m^{N_{l},l}_{k\Delta_{l}}\\}_{k=0}^{T\Delta_{l}^{-1}-1}$. Then using
(3.4), return $\overline{Z}_{T}^{N_{l},l-1}(Y)$ &
$\overline{Z}_{T}^{N_{l},l}(Y)$ or $\overline{U}_{T}^{N_{l},l-1}(Y)$ &
$\overline{U}_{T}^{N_{l},l}(Y)$.

  4. 4.

Output: Return the multilevel estimation of the normalizing constant:

$\displaystyle\overline{Z}_{T}^{ML}(Y)=\overline{Z}_{T}^{N_{l_{*}},l_{*}}(Y)+\sum_{l=l_{*}+1}^{L}\\{\overline{Z}_{T}^{N_{l},l}(Y)-\overline{Z}_{T}^{N_{l},l-1}(Y)\\}.$
(3.12)

or its logarithm:

$\displaystyle\overline{U}_{T}^{ML}(Y)=\overline{U}_{T}^{N_{l_{*}},l_{*}}(Y)+\sum_{l=l_{*}+1}^{L}\\{\overline{U}_{T}^{N_{l},l}(Y)-\overline{U}_{T}^{N_{l},l-1}(Y)\\}.$
(3.13)

##  4 Numerical simulations

In this section we provide various numerical experiments, which include both
verifying the ML rates obtained, and for parameter estimation on both linear
and nonlinear models. The former will be tested on an linear example, which
the later will include a toy linear Gaussian example, a stochastic Lorenz 63
model and Lorenz 96 model, which are common models that arise in atmospheric
sciences. Our parameter estimation will be conducted through a combination of
using recursive maximum likelihood estimation, and simultaneous perturbation
stochastic approximation, which is a gradient-free methodology.

###  4.1 Verification of multilevel rates

We now seek to verify our theory from Proposition 3.1 and Theorem 3.1 on an
Ornstein–Uhlenbeck process, taking the form of (2.1) - (2.2). We will compare
the error-to-cost rates of the MLEnKBF estimation versus the EnKBF estimation
of log-normalizing constant. This will be numerically tested using all
variants (F1) \- (F3), and their respective multilevel counterparts. We will
take $d_{x}=d_{y}=5$, $A=-0.8\leavevmode\nobreak\ Id$,
$\mathcal{M}_{0}=0.1\leavevmode\nobreak\ \textbf{1}$,
$\mathcal{P}_{0}=0.05\leavevmode\nobreak\ Id$, $R^{1/2}=2\leavevmode\nobreak\
Id$, $C$ to be a random matrix and $Q^{1/2}$ is a tri-diagonal matrix defined
as

$Q^{1/2}=\begin{bmatrix}2/3&1/3&&\ldots&0\\\ 1/3&2/3&\ddots&&\\\
&\ddots&\ddots&\ddots&\\\ \vdots&&&\ddots&1/3\\\
0&\ldots&&1/3&2/3\end{bmatrix},$

where $Id$ is the identity matrix of appropriate dimension and 1 is a vector
of ones. We computed the MSE for the target levels $L\in\\{7,8,9,10\\}$ using
416 simulations with $t=1$. In order to better estimate the quantity in (3.4),
we used $l_{*}=6$. The number of samples on each level is given by

$\displaystyle N_{l}=\left\lfloor C_{0}\leavevmode\nobreak\
2^{2L-l}\leavevmode\nobreak\ (L-l_{*}+1)\right\rfloor.$ (4.1)

The reference value is the mean of 208 simulations computed at a
discretization level $l=11$. For the NC estimator, using the EnKBF, the cost
of the estimator is $N_{L}\Delta_{L}^{-1}$. When using the MLEnKBF-NC, the
cost of the estimator is $\sum_{l=l_{*}}^{L}N_{l}\Delta_{l}^{-1}$.

Figure 1: MSE vs Cost plots on a log-log scale for the linear Gaussian model
(2.1) - (2.2), comparing EnKBF and MLEnKBF estimators of the normalizing
constant.

Our numerical results are presented in Figure 1, where we have plotted the ML
and single-level (SL) rates for each variant, within each subplot. Let us
first consider F(1), which is the vanilla MLEnKBF. The results obtained match
that from both Proposition 3.1 and Theorem 3.1, which suggest error-to-cost
rates of $\mathcal{O}(\epsilon^{-3})$ and
$\mathcal{O}(\epsilon^{-2}\log(\epsilon)^{2})$. This is indicated through the
slopes, which show that to attain an order of MSE, the MLEnKBF-NC is
computationally cheaper than that of the EnKBF-NC. Our next subplot concerns
F(2). Despite no existing theory for this case, or its SL counterpart, we see
that the results are similar to the vanilla variant. Specifically, rates given
for both slopes are similar, where there is little distinguishment. Finally
for the final variant of F(2), the results obtained are quite different.
Firstly we notice the slopes are more different to the other variants, but
also that the theory of the rate $\mathcal{O}(\epsilon^{-3})$ does not hold as
convincingly. This could be related to the lack of stochasticity within the
methodology. However from all subplots, one can see that by applying MLMC, one
does attain a lower cost for a particular order of MSE.

###  4.2 Parameter estimation

Let us firstly assume that the model (1.1) - (1.2) contains unknown parameters
$\theta\in\Theta\subseteq\mathbb{R}^{d_{\theta}}$. To account for this, we
rewrite the normalizing constant (3.2) with the additional subscript $\theta$
as $\overline{Z}_{t,\theta}(Y)$. To estimate these parameters we focus on
maximum likelihood inference and stochastic gradient methods that are
performed in an online manner. Mainly, we follow a recursive maximum
likelihood (RML) method, which has been proposed originally in [1], in [29]
for finite spaces, and in [3, 15, 35] in the context of sequential Monte Carlo
(SMC) approximations. Let

$U_{t,t+1,\theta}(Y):=\log\frac{\overline{Z}_{t+1,\theta}(Y)}{\overline{Z}_{t,\theta}(Y)}.$

RML relies on the following update scheme at any time $t\in\mathbb{N}$:

$\displaystyle\theta_{t+1}$
$\displaystyle=\theta_{t}+a_{t}\left(\nabla_{\theta}\log\overline{Z}_{t+1,\theta_{t}}(Y)-\nabla_{\theta}\log\overline{Z}_{t,\theta_{t}}(Y)\right)$
$\displaystyle=\theta_{t}+a_{t}\nabla_{\theta}\leavevmode\nobreak\
U_{t,t+1,\theta_{t}}(Y),$

where $\\{a_{t}\\}_{t\in\mathbb{N}}$ is a sequence of positive real numbers
such that we assume the usual Robbins–Munro conditions, i.e.
$\sum_{t\in\mathbb{N}}a_{t}=\infty$ and
$\sum_{t\in\mathbb{N}}a_{t}^{2}<\infty$. Given an initial
$\theta_{0}\in\Theta$, this formula enables us to update $\theta$ online as we
obtain a new observation path in each unit time interval. Computing the
gradients in the above formula can be expensive (see e.g. [3]), therefore we
use a gradient-free method that is based on some type of finite differences
with simultaneous perturbation stochastic approximation (SPSA) [37, 38]. In a
standard finite difference approach, one perturbs $\theta_{t}$ in the positive
and negative directions of a unit vector $\textbf{e}_{k}$ (a vector of zeros
in all directions except $k$ it is 1). This means evaluating
$U_{t,t+1,\theta_{t}}(Y)$ $2d_{\theta}$–times.

Whereas in SPSA, we perturb $\theta_{t}$ with a magnitude of $b_{t}$ in the
positive and negative directions of a $d_{\theta}$-dimensional random vector
$\Psi_{t}$. The numbers $\\{b_{t}\\}_{t\in\mathbb{N}}$ are a sequence of
positive real numbers such that $b_{t}\to 0$,
$\sum_{t\in\mathbb{N}}a_{t}^{2}/b_{t}^{2}<\infty$, and for
$k\in\\{1\cdots,d_{\theta}\\}$, $\Psi_{t}(k)$ is sampled from a Bernoulli
distribution with success probability $1/2$ and support $\\{-1,1\\}$.
Therefore, this method requires only 2 evaluations of
$U_{t,t+1,\theta_{t}}(Y)$ to estimate the gradient. The MLEnKBF-NC estimator,
presented in Algorithm 1, will be used to estimate the ratio
$U_{t,t+1,\theta_{t}}(Y)$. In Algorithm 2 we illustrate how to implement these
approximations in order to estimate the model’s static parameters.

Algorithm 2 Parameter Estimation: using MLEnKBF-NC and RML-SPSA

  1. 1.

Input: Target level $L\in\mathbb{N}$, start level $l_{*}\in\mathbb{N}$ such
that $l_{*}<L$, the number of particles on each level
$\\{N_{l}\\}_{l=l_{*}}^{L}$, the number of iterations $M\in\mathbb{N}$,
initial $\theta_{0}\in\Theta$, step size sequences of positive real numbers
$\\{a_{t}\\}_{t\in\mathbb{N}}$, $\\{b_{t}\\}_{t\in\mathbb{N}}$ such that
$a_{t},b_{t}\to 0$, $\sum_{t\in\mathbb{N}}a_{t}=\infty$,
$\sum_{t\in\mathbb{N}}a_{t}^{2}/b_{t}^{2}<\infty$, and initial ensembles
$\\{\tilde{\xi}_{0}^{i}\\}_{i=1}^{N_{tot}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mathcal{N}(\mathcal{M}_{0},\mathcal{P}_{0})$, where
$N_{tot}=\sum_{l=l_{*}}^{L}N_{l}$.

  2. 2.

Iterate: For $t\in\\{0,\cdots,M-1\\}$:

     * •

Set
$\\{\xi_{0}^{i,l_{*}}\\}_{i=1}^{N_{l_{*}}}=\\{\tilde{\xi}_{0}^{i}\\}_{i=1}^{N_{l_{*}}}$,
$\cdots$,
$\\{\xi_{0}^{i,L}\\}_{i=1}^{N_{L}}=\\{\tilde{\xi}_{0}^{i}\\}_{i=N_{L-1}+1}^{N_{L}}$.

     * •

For $k\in\\{1,\cdots,d_{\theta}\\}$, sample $\Psi_{t}(k)$ from a Bernoulli
distribution with
success probability $1/2$ and support $\\{-1,1\\}$.

     * •

Set $\theta_{t}^{+}=\theta_{t}+b_{t+1}\Psi_{t}$ and
$\theta_{t}^{-}=\theta_{t}-b_{t+1}\Psi_{t}$.

     * •

Run Algorithm 1 twice, with $T=1$ and initial ensembles
$\Big{\\{}\\{\xi_{0}^{i,l_{*}}\\}_{i=1}^{N_{l_{*}}},\cdots,\\{\xi_{0}^{i,L}\\}_{i=1}^{N_{L}}\Big{\\}}$,
to generate the estimates $U_{t,t+1,\theta_{t}^{+}}(Y)$ and
$U_{t,t+1,\theta_{t}^{-}}(Y)$.

     * •

Set for $k\in\\{1,\cdots,d_{\theta}\\}$,

$\displaystyle\theta_{t+1}(k)$
$\displaystyle=\theta_{t}(k)+\frac{a_{t+1}}{2b_{t+1}\Psi_{t}(k)}\Big{[}U_{t,t+1,\theta_{t}^{+}}(Y)-U_{t,t+1,\theta_{t}^{-}}(Y)\Big{]}.$

     * •

Run the EnKBF up to time 1 under the new parameter $\theta_{t+1}$ with
discretization level $L$ and initial ensembles
$\\{\tilde{\xi}_{0}^{i}\\}_{i=1}^{N_{tot}}$ to return
$\\{\tilde{\xi}_{1}^{i}\\}_{i=1}^{N_{tot}}$.

     * •

Set
$\\{\tilde{\xi}_{0}^{i}\\}_{i=1}^{N_{tot}}=\\{\tilde{\xi}_{1}^{i}\\}_{i=1}^{N_{tot}}$.

We will now introduce our different models we work with. This includes a
linear and two nonlinear models, namely, the stochastic Lorenz 63 and Lorenz
96.

####  4.2.1 Linear Gaussian model

Our first numerical example is a linear Gaussian model with $d_{x}=d_{y}=2$,
based on (2.1) - (2.2). Specifically we choose the parameter values
$A=\theta_{1}Id$ and $Q^{1/2}=\theta_{2}\mathcal{Q}$ where $\mathcal{Q}$ is a
tri-diagnonal matrix defined as

$\mathcal{Q}=\begin{bmatrix}1&1/2&0&\ldots&0\\\ 1/2&1&1/2&&\vdots\\\
0&\ddots&\ddots&\ddots&0\\\ \vdots&&\ddots&\ddots&1/2\\\
0&\ldots&0&1/2&1\end{bmatrix}.$

We also take $C$ to be a uniform random matrix and
$R^{1/2}=0.556\leavevmode\nobreak\ Id$. The parameters of interest we aim to
estimate are $(\theta_{1},\theta_{2})\in\mathbb{R}^{2}$. We choose the target
level as $L=9$ and the start level as $l_{*}=7$, and specify an initial state
of $X_{0}\sim\mathcal{N}(4\mathbf{1},Id)$, where $\mathbf{1}$ is a vectors of
1’s. In all cases, (F1) \- (F3), we take

$\displaystyle N_{l}=\left\lfloor 0.04\leavevmode\nobreak\
2^{2L-l}\leavevmode\nobreak\ (L-l_{*}+1)\right\rfloor.$ (4.2)

The results are displayed in Figure 2. In order to reduce the variance between
the different simulations, we used the same Brownian increments (in the cases
(F1) and (F2)) that were needed to generate $U_{t,t+1,\theta_{t}^{+}}(Y)$ and
$U_{t,t+1,\theta_{t}^{-}}(Y)$ as well. We applied this trick only on the
linear model, but it can be extended to the nonlinear models below. We observe
from Figure 2 that all the variants learn of the true parameters
$(\theta_{1}^{*},\theta_{2}^{*})=(-2,1)$, and from the bottom panel,
$(\textbf{F3})$ has the smallest standard deviation.

######  Remark 4.1.

When testing (F3) in the linear example, there were difficulties in learning
$\theta_{2}$. This is because there is no natural coupling, hence why we had
to choose a specific step-size as described in Figure 2. However for the
nonlinear examples we did not experience this, where no specific modification
of the step-size was required.

Figure 2: Linear results: The outcomes of running Algorithm 2 for the
estimation of $(\theta_{1},\theta_{2})$ in the cases (F1) (top), (F2) (middle)
and (F3) (bottom). The black curve is the average of 6 independent runs and
the shaded area is the mean $\pm$ the standard deviation. The initial values
of the parameters are $(-1,2)$. The dashed lines represent the true parameters
values $(\theta^{*}_{1},\theta^{*}_{2})=(-2,1)$. In all cases we take
$b_{t}=t^{-0.1}$ for all $t\in\mathbb{N}$. In (F1) & (F2) cases, we set
$a_{t}=0.02$ when $t\leq 50$ and $a_{t}=t^{-0.75}$ (for $\theta_{1}$),
$a_{t}=t^{-0.82}$ (for $\theta_{2}$) when $t>50$. In (F3) case, we set
$a_{t}=0.02$ when $t\leq 500$ and $a_{t}=t^{-0.88}$ (for $\theta_{1}$),
$a_{t}=0.2\leavevmode\nobreak\ t^{-0.95}$ (for $\theta_{2}$) when $t>500$.

Figure 3: Lorenz 63 results: The outcomes of running Algorithm 2 for the
estimation of $(\theta_{1},\theta_{2},\theta_{3})$ in the cases (F1) (top),
(F2) (middle) and (F3) (bottom). The black curve is the average of 6
independent runs and the shaded area is the mean $\pm$ the standard deviation.
The initial values of the parameters are $(6,27,6.5)$. The dashed lines
represent the true parameters values
$(\theta^{*}_{1},\theta^{*}_{2},\theta_{3}^{*})=(10,28,8/3)$. The green
horizontal line represent the true parameters values
$(\theta^{*}_{1},\theta^{*}_{2},\theta_{3}^{*})=(10,28,8/3)$. In all cases, we
set $b_{t}=t^{-0.1}$ for all $t\in\mathbb{N}$, $a_{t}=0.01$ when $t\leq 100$
and $a_{t}=t^{-0.75}$ for $t>100$.

####  4.2.2 Stochastic Lorenz 63 model

Our next example is the Lorenz 63 model [30] with $d_{x}=d_{y}=3$, which is a
model for atmospheric convection. The model is based on three ordinary
differential equations, where now we have three parameters of interest to
estimate, i.e. $(\theta_{1},\theta_{2},\theta_{3})\in\mathbb{R}^{3}$. The
stochastic Lorenz 63 model is given as

$\displaystyle dX_{t}$ $\displaystyle=f(X_{t})dt+Q^{1/2}dW_{t},$
$\displaystyle dY_{t}$ $\displaystyle=CX_{t}dt+R^{1/2}dV_{t},$

such that

$\displaystyle f_{1}(X_{t})$ $\displaystyle=\theta_{1}(X_{t}(2)-X_{t}(1)),$
$\displaystyle f_{2}(X_{t})$
$\displaystyle=\theta_{2}X_{t}(1)-X_{t}(2)-X_{t}(1)X_{t}(3),$ $\displaystyle
f_{3}(X_{t})$ $\displaystyle=X_{t}(1)X_{t}(2)-\theta_{3}X_{t}(3),$

where $X_{t}(i)$ is the $i^{th}$ component of $X_{t}$. Furthermore we have
that $Q^{1/2}=Id$ and the variable $C$ is specified as

$C=\begin{cases}\frac{1}{2},\quad&\mathrm{if}\ i=j,\\\
\frac{1}{2},\quad&\mathrm{if}\ i=j-1,\quad j\in\\{1,2,3\\},\\\
0,\quad&\mathrm{otherwise},\end{cases}$

and $(R^{1/2})_{ij}=2\leavevmode\nobreak\
q(\frac{2}{5}\min\\{|i-j|,r_{2}-|i-j|\\}),\ i,j\in\\{1,2,3\\}$ such that

$q(x)=\begin{cases}1-\frac{3}{2}x+\frac{1}{2}x^{3},\quad&\mathrm{if}\ 0\leq
x\leq 1,\\\ 0,\quad&\mathrm{otherwise}.\end{cases}$

In Figure 3, we show the results for the parameters estimation of
$(\theta_{1},\theta_{2},\theta_{3})$ using the cases (F1) \- (F3). We set the
target level to be $L=9$, the start level $l_{*}=7$, and specify the initial
state $X_{0}\sim\mathcal{N}(\textbf{1},0.5\leavevmode\nobreak\ Id)$. The
number of samples on each level are the same as in (4.2). From Figure 3 it is
clear that all variants of the MLEnKBF perform similarly for this inference
problem. Interestingly we notice that the learning of the parameter
$\theta^{*}_{3}=8/3$ seems the most accurate, while taking the least amount of
time to reach the value of $\theta^{*}_{3}$. This is unlike the learning of
the other parameters, which require at least a time of $T=5000$ to get close
to $\theta^{*}_{1}$ and $T=3000$ to get close to $\theta^{*}_{2}$. Part of the
reason for this could be that $\theta_{3}$ acts as a coefficient only for
$X_{t}(3)$, unlike for $X_{t}(1)$ which depends on both $\theta_{1}$ and
$\theta_{2}$. In addition, we see that the learning of $\theta^{*}_{2}$ is
biased towards 27.6 in all variants.

Figure 4: Lorenz 96 results: The outcomes of running Algorithm 2 for the
estimation of $\theta$ in the cases (F1) (left), (F2) (middle) and (F3)
(right). The black curve is the average of 6 independent runs and the shaded
area is the mean $\pm$ the standard deviation. The initial value of $\theta$
is 10. The dashed line represents the true value $\theta^{*}=8$. In all cases,
we set $b_{t}=t^{-0.1}$ for all $t\in\mathbb{N}$, $a_{t}=0.03$ when $t\leq 50$
and $a_{t}=t^{-0.75}$ for $t>50$.

####  4.2.3 Stochastic Lorenz 96 model

Our final test model will be the Lorenz 96 model [31] with $d_{x}=d_{y}=40$,
which is a dynamical system designed to describe equatorial waves in
atmospheric science. The stochastic Lorenz 96 model takes the form

$\displaystyle dX_{t}$ $\displaystyle=f(X_{t})dt+Q^{1/2}dW_{t},$
$\displaystyle dY_{t}$ $\displaystyle=CX_{t}dt+R^{1/2}dV_{t},$

such that

$f_{i}(X_{t})=(X_{t}(i+1)-X_{t}(i-2))X_{t}(i-1)-X_{t}(i)+\theta,$

where again $X_{t}(i)$ is the $i^{th}$ component of $X_{t}$, and we assume
that $X_{t}(-1)=X_{t}(d_{x}-1)$, $X_{t}(0)=X_{t}(d_{x})$ and that
$X_{t}(d_{x}+1)=X_{t}(1)$. We specify our parameter values as
$Q^{1/2}=\sqrt{2}\leavevmode\nobreak\ Id$ and $R^{1/2}=0.5\leavevmode\nobreak\
Id$. The parameter $\theta$ is the external force in the system, while
$(X_{t}(i+1)-X_{t}(i-2))X_{t}(i-1)$ is the advection term and $-X_{t}(i)$ is
the damping term. In Figure 4, we show the results for the parameter
estimation of $\theta$ using the cases (F1)-(F3). We set the target level to
be $L=9$, and the start level $l_{*}=7$. In (F1) and (F2), we specify the
initial state as follows: We set $X_{0}(1)=8.01$ and $X_{0}(i)=8$ for $1<i\leq
d_{x}$. In (F3) case, to prevent the matrix $P_{0}^{N,l}$ from being equal to
zero, we set $X_{0}\sim\mathcal{N}(8\leavevmode\nobreak\
\textbf{1},0.05\leavevmode\nobreak\ Id)$. The number of samples on each level
are the same as in (4.2). As we observe from Figure 4, all MLEnKBF variants
learn the true value of $\theta^{*}=8$. The fluctuations again can be seen
more clearly in the first two variants, as they are stochastic, i.e. contain
either or both of $Q$ and $R$. This is not the case for (F3) whose dynamics
are completely deterministic. In general the results are more stable that
those conducted for the Lorenz 63 model, as the function
$f:\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{d_{x}}$, does not contain $\theta$
as a coefficient of the process $X_{t}$.

##  5 Conclusion

The purpose of this work was to apply MLMC strategies for normalizing constant
(NC) estimation. In particular our aim was to extend the work of [12], which
used EnKBF, to its multilevel counterpart which is the MLEnKBF [10]. As
stated, our motivation is primarily in the linear setting, for which the
optimal filter is the Kalman–Bucy filter, where we could provide propagation
of chaos-type results. This has the advantage of being more computationally
feasible than using other filtering strategies such as the particle filter or
SMC methods. From our results, we analyzed and provided various
$\mathbb{L}_{q}$ bounds associated to the i.i.d. normalizing constant
estimator. This was for both the single-level and multilevel cases. Numerical
experiments were conducted which showed the improvement, related to cost, of
adopting the multilevel methodology for the EnKBF. This was further tested on
various models for parameter estimation, where we compared different variants
of MLEnKBFs, introduced earlier on the stochastic Lorenz 63 and 96 models.

This work naturally leads to different fruitful and future directions of work
to consider. The first potential direction is to aim to provide the same
analysis for other MLEnKBFs, such as the deterministic and deterministic-
transport EnKBF [5]. This is of interest as from the numerical results, they
suggest alternative rates for (F3), different to that of the other variants.
One could also provide an unbiased estimation of the NC [39], which has
connections to MLMC [9, 40]. In order to do so, one would require a modified
multilevel analysis of the EnKBF, where one has uniform upper bounds, with
respect to levels $l=1,\ldots,L$. As stated, a comparison of both normalizing
estimators using both the MLEnKBF and MLPF would be interesting, but to make
such a comparison, one could exploit advanced methodologies, such as in [3],
and appropriate examples in low dimensions. Finally it would be of interest to
develop theory for the multilevel estimator $\overline{U}^{ML}_{t}$. We have
only tested this computationally which seems to match the rates attained for
the i.i.d. (ideal) ML estimator. To do so, one requires more sophisticated
mathematics to hold initially for the problem-setting in [10].

## Appendix

##  Appendix A Analysis for EnKBF NC Estimator

For the appendices, Appendix A will cover the propagation of chaos result,
which is required for the variance of the single-level EnKBF NC estimator. We
will then proceed to Appendix B which discusses various discretization biases
of the diffusion process related to both EnKBF and the NC estimator. Finally
our main theorem is proved in Appendix C. All of our results will be specific
to the vanilla variant of the ENKBF, F(1).

Before proceeding to our results, we will introduce the following assumptions
which will hold from herein, but not be added to any statements. For a square
matrix, $B$ say, we denote by $\mu(B)$ as the maximum eigenvalue of
$\textrm{Sym}(B)$.

  1. 1.

We have that $\mu(A)<0$.

  2. 2.

There exists a $\mathsf{C}<+\infty$ such that for any
$(k,l)\in\mathbb{N}_{0}^{2}$ we have that

$\max_{(j_{1},j_{2})\in\\{1,\dots,d_{x}\\}^{2}}|P_{k\Delta_{l}}(j_{1},j_{2})|\leq\mathsf{C}.$
(A.1)

We note that 1. is typically used in the time stability of the hidden
diffusion process $X_{t}$, see for instance [17]. In the case of 2. we expect
that it can be verified under 1., that $S=\mathsf{C}I$ with $\mathsf{C}$ a
positive constant and some controllability and observability assumptions (e.g.
[17, eq. (20)]). Under such assumptions, the Riccati equation has a solution
and moreover, by [17, Proposition 5.3] $\mathcal{P}_{t}$ is exponentially
stable w.r.t. the Frobenius norm; so that this type of bound exists in
continuous time.

Throughout the appendix we will make use of the $C_{q}-$inequality. For two
real-valued random variables $X$ and $Y$ defined on the same probability
space, with expectation operator $\mathbb{E}$, suppose that for some fixed
$q\in(0,\infty)$, $\mathbb{E}[|X|^{q}]$ and $\mathbb{E}[|Y|^{q}]$ are finite,
then the $C_{q}-$inequality is

$\mathbb{E}[|X+Y|^{q}]\leq\mathsf{C}_{q}\Big{(}\mathbb{E}[|X|^{q}]+\mathbb{E}[|Y|^{q}]\Big{)},$

where $\mathsf{C}_{q}=1$, if $q\in(0,1)$ and $\mathsf{C}_{q}=2^{q-1}$ for
$q\in[1,\infty)$.

In order to verify some of our claims for the analysis, we will rely on
various results derived in [10]. For convenience-sake we will state these
below, which are concerned with various $\mathbb{L}_{q}-$bounds, where

######  Lemma A.1.

For any $q\in[1,\infty)$ there exists a $\mathsf{C}<+\infty$ such that for any
$(k,l,j)\in\mathbb{N}_{0}^{2}\times\\{1,\dots,d_{y}\\}$:

$\mathbb{E}[|[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}](j)|^{q}]^{1/q}\leq\mathsf{C}\Delta_{l}^{1/2}.$

######  Lemma A.2.

For any $(q,t,k,l)\in[1,\infty)\times\mathbb{N}_{0}^{3}$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(j,N)\in\\{1,\dots,d_{x}\\}\times\\{2,3,\dots\\}$:

$\mathbb{E}\Big{[}\Big{|}m_{t+k_{1}\Delta_{l}}^{N}(j)-m_{t+k_{1}\Delta_{l}}(j)\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}}{\sqrt{N}}.$

######  Lemma A.3.

For any $(q,k,l)\in(0,\infty)\times\mathbb{N}_{0}^{2}$ there exists a
$\mathsf{C}<+\infty$ such that for any $N\geq 2$ and $i\in\\{1,\dots,N\\}$:

$\mathbb{E}[|\xi_{k\Delta_{l}}^{i}(j)|^{q}]^{1/q}\leq\mathsf{C},$

where $\xi_{k\Delta_{l}}$ is defined through (2.11).

We now present our first result for the single-level EnKBF NC estimator, which
is presented as an $\mathbb{L}_{q}-$ error bound.

######  Proposition A.1.

For any $(q,t,k_{1},l)\in[1,\infty)\times\mathbb{N}_{0}^{3}$ there exists a
$\mathsf{C}<+\infty$ such that for any $N\in\\{2,3,\ldots\\}$ we have:

$\mathbb{E}\Big{[}\Big{|}[\overline{U}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)]\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}}{\sqrt{N}}.$

###### Proof.

Let us first consider
$\overline{U}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)$,
which for every $l\in\mathbb{N}_{0}$, we can decompose through a martingale
remainder-type decomposition,

$\overline{U}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)=M_{t+k_{1}\Delta_{l}}^{l}(Y)+R_{t+k_{1}\Delta_{l}}^{l},$
(A.2)

such that

$\displaystyle M_{t+k_{1}\Delta_{l}}^{l}(Y)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\langle
Cm_{k\Delta_{l}}^{N},R^{-1}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]\rangle-\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\langle
Cm_{k\Delta_{l}},R^{-1}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]\rangle,$
$\displaystyle R_{t+k_{1}\Delta_{l}}^{l}$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\langle
m_{k\Delta_{l}}^{N},Sm_{k\Delta_{l}}^{N}\rangle+\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\langle
m_{k\Delta_{l}},Sm_{k\Delta_{l}}\rangle.$

We can decompose the martingale term from (A.2) further through

$M_{t+k_{1}\Delta_{l}}^{l}(Y)=M_{t+k_{1}\Delta_{l}}^{l}(1)+R_{t+k_{1}\Delta_{l}}^{l}(1),$

where, by setting ${k=t\Delta_{l}^{-1}+k_{1}-1}$,

$\displaystyle M_{t+k_{1}\Delta_{l}}^{l}(1)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}C(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}R^{-1}(j_{1},j_{3})\times$
$\displaystyle\\{[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}](j_{3})-CX_{k\Delta_{l}}(j_{3})\Delta_{l}\\},$
$\displaystyle R_{t+k_{1}\Delta_{l}}^{l}(1)$ $\displaystyle=$
$\displaystyle\Delta_{l}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}C(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}\times$
$\displaystyle R^{-1}(j_{1},j_{3})CX_{k\Delta_{l}}(j_{3}).$

In order to proceed we construct a martingale associated with the term of
$M_{t}^{l}$. Let us first begin with the $M_{t}^{l}(1)$ term (A), where we
construct the filtration
$(\Omega,\mathscr{F},\mathscr{F}_{k\Delta_{l}},\mathbb{P})$ for our discrete-
time martingale $(M^{l}_{t}(1),\mathscr{F}_{k\Delta_{l}})$.

Then by using Hölder’s inequality

$\displaystyle\mathbb{E}[|M^{l}_{t+k_{1}\Delta_{l}}(1)|^{q}]^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{dy}C(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}R^{-1}(j_{1},j_{3})\times$
(A.5)
$\displaystyle\\{[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}](j_{3})-CX_{k\Delta_{l}}(j_{3})\Delta_{l}\\}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}C(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}R^{-1}(j_{1},j_{3})\Big{|}^{2q}\Big{]}^{1/2q}\times$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{3}=1}^{dy}\\{[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}](j_{3})-CX_{k\Delta_{l}}(j_{3})\Delta_{l}\\}\Big{|}^{2q}\Big{]}^{1/2q}$
$\displaystyle=:$ $\displaystyle T_{1}\times T_{2}.$

For $T_{1}$ we can apply the Minkowski inequality and Lemma A.2 to yield

$\displaystyle T_{1}$ $\displaystyle\leq$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}C(j_{1},j_{2})R^{-1}(j_{1},j_{3})\mathbb{E}\Big{[}\Big{|}\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}\Big{|}^{2q}\Big{]}^{1/2q}$
$\displaystyle\leq$ $\displaystyle\frac{\mathsf{C}}{\sqrt{N}}.$

For $T_{2}$ we know that the expression
$\\{[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}](j_{3})-CX_{k\Delta_{l}}(j_{3})\Delta_{l}\\}$
is a Brownian motion
increment, using the formulae (2.1) - (2.2). Therefore by using the
Burkholder–Davis–Gundy inequality, along with Minkowski, for $\tilde{q}=2q$,
we have

$\displaystyle T_{2}$ $\displaystyle\leq$
$\displaystyle\sum_{j_{3}=1}^{d_{y}}\mathbb{E}\Big{[}\Big{|}\sum^{t\Delta^{-1}_{l}+k_{1}-1}_{k=0}[V_{(k+1)\Delta_{l}}-V_{k\Delta_{l}}](j_{3})\Big{|}^{\tilde{q}}\Big{]}^{1/\tilde{q}}$
$\displaystyle\leq$
$\displaystyle\sum_{j_{3}=1}^{d_{y}}\sum^{t\Delta^{-1}_{l}+k_{1}-1}_{k=0}\mathsf{C}_{\tilde{q}}\mathbb{E}\Big{[}\Big{|}[V_{(k+1)\Delta_{l}}-V_{k\Delta_{l}}]^{2}(j_{3})\Big{|}^{\tilde{q}/2}\Big{]}^{1/\tilde{q}}$
$\displaystyle\leq$
$\displaystyle\sum_{j_{3}=1}^{d_{y}}\sum^{t\Delta^{-1}_{l}+k_{1}-1}_{k=0}\mathsf{C}_{\tilde{q}}\Big{(}\mathbb{E}\Big{[}\Big{|}[V_{(k+1)\Delta_{l}}-V_{k\Delta_{l}}](j_{3})\Big{|}^{\tilde{q}}\Big{]}^{2/\tilde{q}}\Big{)}^{1/2}$

Then using that fact that
$\mathbb{E}[|[V_{(k+1)\Delta_{l}}-V_{k\Delta_{l}}]|^{\tilde{q}}]=O(\Delta_{l}^{\tilde{q}/2})$,
and with the summation it is of order
$\mathcal{O}(\Delta^{(1/2-1/\tilde{q})})$, we can conclude $T_{2}$ is of order
$\mathcal{O}(1)$, as we have used the case of when $q\geq 1$.

For the $R_{t}^{l}(1)$ term, it follows similarly to $M_{t}^{l}(1)$, where we
require the use of Lemma A.2.

Again we make use of the Minkowski and Hölder inequality, and Lemma A.2,

$\displaystyle\mathbb{E}[|R^{l}_{t+k_{1}\Delta_{l}}(1)|^{q}]^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}C(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}R^{-1}(j_{1},j_{3})CX_{k\Delta_{l}}(j_{3})\Delta_{l}\Big{|}^{q}\Big{]}^{1/q}$
(A.6) $\displaystyle\leq$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum^{d_{x}}_{j_{4}=1}C(j_{1},j_{3})R^{-1}(j_{1},j_{3})\Big{(}\mathbb{E}\Big{[}\Big{|}m^{N}_{k\Delta
l}(j_{2})-m_{k\Delta l}(j_{2})\Big{|}^{2q}\Big{]}^{1/2q}\times$
$\displaystyle\mathbb{E}\Big{[}\Big{|}X_{k\Delta_{l}}(j_{4})\Delta_{l}\Big{|}^{2q}\Big{]}^{1/2q}\Big{)}$
$\displaystyle\leq$
$\displaystyle\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum^{d_{x}}_{j_{4}=1}\frac{\mathsf{C}}{\sqrt{N}}\mathbb{E}\Big{[}\Big{|}X_{k\Delta_{l}}(j_{4})\Delta_{l}\Big{|}^{2q}\Big{]}^{1/2q}.$

For the final term of (A.6) we can show it is of order $\mathcal{O}(1)$, using
the Cauchy–Schwarz
and Jensen’s inequality

$\displaystyle\mathbb{E}\Big{[}\Big{|}X_{k\Delta_{l}}(j_{4})\Delta_{l}\Big{|}^{2q}\Big{]}^{1/2q}$
$\displaystyle=$
$\displaystyle\Delta_{l}\mathbb{E}\Big{[}\Big{|}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}{\Delta_{l}^{-1}}X_{s}(j_{4})ds\Big{|}^{2q}\Big{]}^{1/2q}$
$\displaystyle\leq$
$\displaystyle\Delta_{l}\mathbb{E}\Big{[}\Big{(}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}\Delta^{-2}_{l}ds\Big{)}^{q}\Big{(}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}|X_{s}(j_{4})|^{2}ds\Big{)}^{q}\Big{]}^{1/2q}$
$\displaystyle=$
$\displaystyle\Delta^{1/2}_{l}\mathbb{E}\Big{[}\Big{(}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}|X_{s}(j_{4})|^{2}ds\Big{)}^{q}\Big{]}^{1/2q}$
$\displaystyle\leq$
$\displaystyle\Delta_{l}^{1/2}\mathbb{E}\Big{[}{\Delta_{l}^{q-1}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}}\Big{|}X_{s}(j_{4})\Big{|}^{2q}ds\Big{]}^{1/2q}$
$\displaystyle\leq$
$\displaystyle\Delta_{l}^{1-1/2q}\Big{(}\int^{(k+1)\Delta_{l}}_{k\Delta_{l}}{\mathbb{E}[|X_{s}(j_{4})|^{2q}]}ds\Big{)}^{1/2q}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\Delta_{l}^{1-1/2q}(\Delta_{l})^{1/2q}\quad\textrm{(By
\cite[cite]{[\@@bibref{}{DT18}{}{}]} Eq. 54)}$ $\displaystyle=$
$\displaystyle\mathsf{C}\Delta_{l},$

therefore combing this with the summation in (A.6), the above quantity is of
order $\mathcal{O}(1)$, resulting in
$\mathbb{E}[|R^{l}_{t+k_{1}\Delta_{l}}(1)|^{q}]^{1/q}\leq\frac{\mathsf{C}}{\sqrt{N}}$.
All that is left is the $R^{l}_{t+k_{1}\Delta_{l}}$ term. Before proceeding we
can express, or rewrite, the $R_{t+k_{1}\Delta_{l}}^{l}$ term as

$\displaystyle R_{t+k_{1}\Delta_{l}}^{l}$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\langle
m_{k\Delta_{l}}^{N},Sm_{k\Delta_{l}}^{N}\rangle+\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}-1}\langle
m_{k\Delta_{l}},Sm_{k\Delta_{l}}\rangle$ $\displaystyle=$
$\displaystyle-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\\{m_{k\Delta_{l}}^{N}(j_{1})S(j_{1},j_{2})m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{1})S(j_{1},j_{2})m_{k\Delta_{l}}(j_{2})\\}$
$\displaystyle=$
$\displaystyle-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\\{m_{k\Delta_{l}}^{N}(j_{1})-m_{k\Delta_{l}}(j_{1})\\}S(j_{1},j_{2})m_{k\Delta_{l}}^{N}(j_{2})$
$\displaystyle+$ $\displaystyle
m_{k\Delta_{l}}(j_{1})S(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}.$

By taking its associated $\mathbb{L}_{q}$-bound, from Minkowski’s inequality
and Lemma A.2 \- Lemma A.3, we have

$\displaystyle\mathbb{E}[|R_{t+k_{1}\Delta_{l}}^{l}|^{q}]^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\\{m_{k\Delta_{l}}^{N}(j_{1})-m_{k\Delta_{l}}(j_{1})\\}S(j_{1},j_{2})m_{k\Delta_{l}}^{N}(j_{2})$
$\displaystyle+$ $\displaystyle
m_{k\Delta_{l}}(j_{1})S(j_{1},j_{2})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle-\frac{\Delta_{l}}{2}\sum_{k=0}^{t\Delta_{l}^{-1}+k_{1}-1}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}S(j_{1},j_{2})\Big{(}\mathbb{E}\Big{[}\Big{|}\\{m_{k\Delta_{l}}^{N}(j_{1})-m_{k\Delta_{l}}(j_{1})\\}m_{k\Delta_{l}}^{N}(j_{2})\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle+$
$\displaystyle\mathbb{E}\Big{[}\Big{|}m_{k\Delta_{l}}(j_{1})\\{m_{k\Delta_{l}}^{N}(j_{2})-m_{k\Delta_{l}}(j_{2})\\}\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle\leq$ $\displaystyle\frac{\mathsf{C}}{\sqrt{N}}.$

Finally by using the Minkowski inequality, we can deduce that

$\displaystyle\mathbb{E}\Big{[}\Big{|}[\overline{U}_{t+k_{1}\Delta_{l}}^{N,l}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)]\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}M_{t+k_{1}\Delta_{l}}^{l}(1)+R_{t+k_{1}\Delta_{l}}^{l}(1)+R_{t+k_{1}\Delta_{l}}^{l}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}[|M^{l}_{t+k_{1}\Delta_{l}}(1)|^{q}]^{1/q}+\mathbb{E}[|R^{l}_{t+k_{1}\Delta_{l}}(1)|^{q}]^{1/q}+\mathbb{E}[|R^{l}_{t+k_{1}\Delta_{l}}|^{q}]^{1/q}$
$\displaystyle\leq$ $\displaystyle\frac{\mathsf{C}}{\sqrt{N}}.$

∎

##  Appendix B Analysis for Discretized Diffusion Process

In this appendix we consider deriving analysis for the discretized diffusion
process. This will include both the discretized i.i.d. particle system

$\displaystyle\zeta_{(k+1)\Delta_{l}}^{i}$
$\displaystyle=(I+A\Delta_{l})\zeta_{k\Delta_{l}}^{i}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]$
$\displaystyle+P_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]-\Big{[}C\zeta_{k\Delta_{l}}^{i}\Delta_{l}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}]\Big{]}\Big{)},$
(B.1)

and the discretized NC estimator. We recall, in the limit as
$N\rightarrow\infty$, the i.i.d. system coincides with discretized Kalman–Bucy
diffusion whose mean, for $(k,l)\in\mathbb{N}_{0}\times\mathbb{N}_{0}$, is
defined by

$m_{(k+1)\Delta_{l}}^{l}=m_{k\Delta_{l}}^{l}+Am^{l}_{k\Delta_{l}}\Delta_{l}+U_{k\Delta_{l}}^{l}\Big{(}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]-Cm_{k\Delta_{l}}^{l}\Delta_{l}\Big{)},$
(B.2)

We note in this appendix our results will use the notation $\overline{X}$ for
the Kalman–Bucy diffusion, to keep it consistent with [10]. However these
results also hold for the i.i.d. system (B.1). We require additional lemmas
from [10], which are discretization bias results for the discretized
Kalman–Bucy diffusion. We state these as follows. The notation of the
equations are modified for the multilevel, which we will discuss later.

######  Lemma B.1.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,j_{1},j_{2})\in\mathbb{N}_{0}\times\\{1,\dots,d_{x}\\}^{2}$:

$\Big{|}\mathcal{P}_{t}(j_{1},j_{2})-P_{\tau_{t}^{l}}^{l}(j_{1},j_{2})\Big{|}\leq\mathsf{C}\Delta_{l}.$

######  Lemma B.2.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T-1]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,j,k_{1})\in\mathbb{N}_{0}\times\\{1,\dots,d_{x}\\}\times\\{0,1,\dots,\Delta_{l}^{-1}\\}$:

$\mathbb{E}\Big{[}\Big{(}\overline{X}_{t+k_{1}\Delta_{l}}(j)-\overline{X}_{t+k_{1}\Delta_{l}}^{l}(j)\Big{)}^{2}\Big{]}\leq\mathsf{C}\Delta_{l}^{2}.$

We now present our first conditional bias result, which will be the weak error
of the Kalman–Bucy diffusion. This weak error will be analogous to the strong
error of Lemma B.2, which was not proved, or provided, in [10]. However this
result will be required for the current and succeeding appendix.

######  Lemma B.3.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T-1]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,k_{1})\in\mathbb{N}_{0}\times\\{1,\dots,d_{x}\\}\times\\{0,1,\dots,\Delta_{l}^{-1}\\}$:

$\Big{|}\mathbb{E}\Big{[}\overline{X}_{t+k_{1}\Delta_{l}}(j)-\overline{X}_{t+k_{1}\Delta_{l}}^{l}(j)\Big{]}\Big{|}\leq\mathsf{C}\Delta_{l}.$

###### Proof.

As before we can separate the above expression in different terms,

$\mathbb{E}\Big{[}\overline{X}_{t+k_{1}\Delta_{l}}(j)-\overline{X}_{t+k_{1}\Delta_{l}}^{l}(j)\Big{]}=T_{1}+T_{2}+T_{3},$

such that, for $\tau_{t}^{l}=[\tfrac{t}{\Delta_{l}}]\Delta_{l}$,
$t\in\mathbb{R}^{+}$, we have

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\int_{0}^{t+k_{1}\Delta_{l}}\Big{(}\sum_{j_{1}=1}^{d_{x}}A(j,j_{1})[\overline{X}_{s}(j_{1})-\overline{X}_{\tau_{s}^{l}}^{l}(j_{1})]+\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\mathcal{P}_{s}(j,j_{1})\hat{C}(j_{1},j_{2})[\overline{X}_{s}(j_{2})-\overline{X}_{\tau_{s}^{l}}^{l}(j_{2})]\Big{)}ds\Big{]},$
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\int_{0}^{t+k_{1}\Delta_{l}}\Big{(}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}[\mathcal{P}_{s}(j,j_{1})-P_{\tau_{s}^{l}}^{l}(j,j_{1})]\hat{C}(j_{1},j_{2})\overline{X}_{\tau_{s}^{l}}^{l}(j_{2})]\Big{)}ds\Big{]},$
$\displaystyle T_{3}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\sum^{d_{x}}_{j_{1}=1}\sum^{d_{y}}_{j_{2}=1}\int_{0}^{t+k_{1}\Delta_{l}}[\mathcal{P}_{s}(j,j_{1})-P_{\tau_{s}^{l}}^{l}(j,j_{1})]\tilde{C}(j_{1},j_{2})dY_{s}(j_{2})\Big{]},$

where $\hat{C}=\tilde{C}C$ with $\tilde{C}=C^{\top}R_{2}^{-1}$. Now let us
consider each individual term,
firstly for $T_{1}$, that (A.1) we have the following bound

$|T_{1}|\leq\mathsf{C}\int_{0}^{t+k_{1}\Delta_{l}}\max_{j\in\\{1,\dots,d_{x}\\}}\Big{|}\mathbb{E}\Big{[}\overline{X}_{s}(j)-\overline{X}_{\tau_{s}^{l}}^{l}(j)\Big{]}\Big{|}ds.$
(B.3)

For $T_{2}$, We can apply Lemma B.1, using the fact that
$\max_{j_{2}\in\\{1,\dots,d_{x}\\}}|\mathbb{E}[\overline{X}_{\tau_{s}^{l}}^{l}(j_{2})]|\leq\mathsf{C}$
we have

$|T_{2}|\leq\mathsf{C}\Delta_{l}.$ (B.4)

Similarly for $T_{3}$, we can use Lemma B.1 and Lemma A.1 which provides the
bound

$|T_{3}|\leq\mathsf{C}\Delta_{l}.$ (B.5)

Thus combining (B.3)-(B.5) leads to

$\max_{j\in\\{1,\dots,d_{x}\\}}\Big{|}\mathbb{E}\Big{[}\overline{X}_{t+k_{1}\Delta_{l}}(j)-\overline{X}_{t+k_{1}\Delta_{l}}^{l}(j)\Big{]}\Big{|}\leq\mathsf{C}\Delta_{l}+\mathsf{C}\int_{0}^{t+k_{1}\Delta_{l}}\max_{j\in\\{1,\dots,d_{x}\\}}\Big{|}\mathbb{E}\Big{[}\overline{X}_{s}(j)-\overline{X}_{\tau_{s}^{l}}^{l}(j)\Big{]}\Big{|}ds.$

Finally by applying Grönwall’s lemma, leads to the desired result. ∎

We now proceed with our result of the discretized NC estimator, which is the
strong error, through the following lemma.

######  Lemma B.4.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T-1]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,k_{1})\in\mathbb{N}_{0}\times\\{0,1,\dots,\Delta_{l}^{-1}\\}$:

$\mathbb{E}\Big{[}\Big{(}[\overline{U}_{t+k_{1}\Delta_{l}}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)]\Big{)}^{2}\Big{]}\leq\mathsf{C}\Delta^{2}_{l}.$

###### Proof.

Let us first recall that,

$\displaystyle\overline{U}^{l}_{t+k_{1}\Delta_{l}}(Y)$ $\displaystyle=$
$\displaystyle\sum^{t\Delta^{-1}_{l}+k_{1}-1}_{k=0}\langle
Cm_{k\Delta_{l}},R^{-1}[Y_{(k+1)\Delta_{l}}-Y_{k\Delta_{l}}]\rangle-\frac{\Delta_{l}}{2}\sum^{t\Delta^{-1}_{l}+k_{1}-1}_{k=0}\langle
m_{k\Delta_{l}},Sm_{k\Delta_{l}}\rangle,$
$\displaystyle\overline{U}_{t+k_{1}\Delta_{l}}(Y)$ $\displaystyle=$
$\displaystyle\int^{t+k_{1}\Delta_{l}}_{0}\Big{[}\langle
Cm_{s},R^{-1}dY_{s}\rangle-\frac{1}{2}\langle m_{s},Sm_{s}\rangle ds\Big{]}.$

In order to proceed we again consider a martingale-remainder type
decomposition. Therefore by setting
$\tau_{t}^{l}=[\tfrac{t}{\Delta_{l}}]\Delta_{l}$, $t\in\mathbb{R}^{+}$, and
expanding on the angle brackets, we have

$\displaystyle M_{t+k_{1}\Delta_{l}}(1)$ $\displaystyle=$
$\displaystyle\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\int^{t+k_{1}\Delta_{l}}_{0}R^{-1}(j_{1},j_{3})C(j_{1},j_{2})\Big{[}m_{s}(j_{2})-m^{l}_{\tau_{s}^{l}}(j_{2})\Big{]}R^{1/2}(j_{1},j_{3})dV_{s}(j_{3}),$
$\displaystyle R_{t+k_{1}\Delta_{l}}$ $\displaystyle=$
$\displaystyle\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l}}_{0}\frac{S(j_{1},j_{2})}{2}\Big{[}m_{s}(j_{1})m_{s}(j_{2})-m^{l}_{\tau_{s}^{l}}(j_{1})m^{l}_{\tau_{s}^{l}}(j_{2})\Big{]}ds,$
$\displaystyle R_{t+k_{1}\Delta_{l}}(1)$ $\displaystyle=$
$\displaystyle\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\int^{t+k_{1}\Delta_{l}}_{0}R^{-1}(j_{1},j_{3})C(j_{1},j_{2})\Big{[}m_{s}(j_{1})-m^{l}_{\tau_{s}^{l}}(j_{1})\Big{]}CX_{s}(j_{3})ds,$

where we have used the formula for the observational process (2.1), combined
both remainder terms into one, and taken the scaled Brownian motion $V_{t}$.
Let us first consider the remainder term of $R_{t+k_{1}\Delta_{l}}(1)$.
Through Jensen’s inequality we have

$\displaystyle\mathbb{E}[R_{t+k_{1}\Delta_{l}}(1)]^{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\int^{t+k_{1}\Delta_{l}}_{0}R^{-1}(j_{1},j_{3})C(j_{1},j_{2})\Big{[}m_{s}(j_{2})-m^{l}_{\tau_{s}^{l}}(j_{2})\Big{]}CX_{s}(j_{3})ds\Big{]}^{2}$
$\displaystyle\leq$
$\displaystyle\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum_{j_{4}=1}^{d_{x}}R^{-2}(j_{1},j_{3})C^{2}(j_{1},j_{2})\int^{t+k_{1}\Delta_{l}}_{0}\mathbb{E}\Big{[}[m_{s}(j_{2})-m^{l}_{\tau_{s}}(j_{2})]X_{s}(j_{4})\Big{]}^{2}ds.$

Then by using $d^{2}_{x}d^{2}_{y}$ applications of the $C_{2}$-inequality we
get

$\displaystyle\mathbb{E}[R_{t+k_{1}\Delta_{l}}(1)]^{2}$ $\displaystyle\leq$
$\displaystyle\mathsf{C}\Big{(}\int^{t+k_{1}\Delta_{l}}_{0}\max_{j_{2}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}[m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})]^{2}\max_{j_{4}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}[X_{s}(j_{4})]^{2}ds\Big{)}.$

We know that
$\max_{j_{4}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}[X_{s}(j_{4})]^{2}\leq\mathsf{C}$,
therefore all we need is to bound
$|m_{s}(j_{1})-m^{l}_{\tau^{l}_{s}}(j_{1})|^{2}$. Therefore by using the fact
that $\mathbb{E}[\overline{X}_{s}(j)]=m_{s}(j)$ and
$\mathbb{E}[\overline{X}^{l}_{\tau^{l}_{s}}(j)]=m^{l}_{\tau^{l}_{s}}(j)$, we
can use the weak error, i.e. Lemma B.3, to conclude that
$\mathbb{E}[R_{t+k_{1}\Delta_{l}}(1)]^{2}\leq\mathsf{C}\Delta^{2}_{l}$.

Now to proceed with $R_{t+k_{1}\Delta_{l}}(Y)$, we can split the difference of
the mean term

$\displaystyle
m_{s}(j_{1})m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{1})m^{l}_{\tau^{l}_{s}}(j_{2})$
$\displaystyle=$ $\displaystyle
m_{s}(j_{1})m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{1})m_{s}(j_{2})+m^{l}_{\tau^{l}_{s}}(j_{1})m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{1})m^{l}_{\tau^{l}_{s}}(j_{2})$
$\displaystyle=$
$\displaystyle\\{m_{s}(j_{1})-m^{l}_{\tau^{l}_{s}}(j_{1})\\}m_{s}(j_{2})+m^{l}_{\tau^{l}_{s}}(j_{1})\\{m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})\\}.$
(B.6)

Therefore one can substitute (B.6) into $R_{t}(Y)$, and by Jensen’s and the
$C_{2}$-inequality, results in

$\displaystyle\mathbb{E}[R_{t+k_{1}\Delta_{l}}]^{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l}}_{0}\frac{S(j_{1},j_{2})}{2}\Big{[}\\{m_{s}(j_{1})-m^{l}_{\tau^{l}_{s}}(j_{1})\\}m_{s}(j_{2})+m^{l}_{\tau^{l}_{s}}(j_{1})\\{m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})\\}\Big{]}ds\Big{]}^{2}$
(B.7) $\displaystyle\leq$
$\displaystyle\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\frac{S^{2}(j_{1},j_{2})}{2}\int^{t+k_{1}\Delta_{l}}_{0}\mathbb{E}\Big{[}\\{m_{s}(j_{1})-m^{l}_{\tau^{l}_{s}}(j_{1})\\}m_{s}(j_{2})\Big{]}^{2}ds$
$\displaystyle+$
$\displaystyle\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\frac{S^{2}(j_{1},j_{2})}{2}\int^{t+k_{1}\Delta_{l}}_{0}\mathbb{E}\Big{[}m^{l}_{\tau^{l}_{s}}(j_{1})\\{m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})\\}\Big{]}^{2}ds.$

Therefore by again taking the $C_{2}$-inequality $d^{2}_{x}$ times, and noting
that
$\max_{j_{2}\in\\{1,\dots,d_{x}\\}}m_{s}(j_{2})=\max_{j_{2}\in\\{1,\dots,d_{x}\\}}\mathbb{E}[\overline{X}_{s}(j_{2})]\leq\mathsf{C}$,
from Lemma A.2 and similarly done for $m^{l}_{\tau^{l}_{s}}(j_{1})$, leads to

$\displaystyle\mathbb{E}[R_{t+k_{1}\Delta_{l}}]^{2}\leq\mathsf{C}\Big{(}\int^{t+k_{1}\Delta_{l}}_{0}\max_{j_{1}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}\Big{[}m_{s}(j_{1})-m^{l}_{\tau^{l}_{s}}(j_{1})\Big{]}^{2}+\max_{j_{2}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}\Big{[}m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})\Big{]}^{2}ds\Big{)}.$

Through the same substitution as before, and using the weak error, i.e. Lemma
B.3,
we have that
$\mathbb{E}[R_{t+k_{1}\Delta_{l}}]^{2}\leq\mathsf{C}\Delta^{2}_{l}$.

Lastly we have the martingale term $M_{t+k_{1}\Delta_{l}}(1)$. As before we
can apply Jensens inequality

$\displaystyle\mathbb{E}[M_{t+k_{1}\Delta_{l}}(1)]^{2}$ $\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\int^{t+k_{1}\Delta_{l}}_{0}R^{-1}(j_{1},j_{3})C(j_{1},j_{2})\Big{[}m_{s}(j_{2})-m^{l}_{\tau_{s}^{l}}(j_{2})\Big{]}R^{1/2}(j_{1},j_{3})dV_{s}(j_{3})\Big{]}^{2}$
$\displaystyle\leq$
$\displaystyle\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}{R}^{-1}(j_{1},j_{3})\
C^{2}(j_{1},j_{2})\int^{t+k_{1}\Delta_{l}}_{0}\mathbb{E}\Big{[}m_{s}(j_{2})-m^{l}_{\tau_{s}^{l}}(j_{2})dV_{s}(j_{3})\Big{]}^{2}.$

Then by using the Ito isometry and $d_{x}d^{2}_{y}$ applications of the
$C_{2}$-inequality, we have

$\displaystyle\mathbb{E}[M_{t+k_{1}\Delta_{l}}(1)]^{2}$ $\displaystyle\leq$
$\displaystyle\mathsf{C}\int^{t+k_{1}\Delta_{l}}_{0}\max_{j_{2}\in\\{1,\ldots,d_{x}\\}}\mathbb{E}\Big{[}m_{s}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{2})\Big{]}^{2}ds,$

and, as before, by using Lemma B.3 as done for $R_{t+k_{1}\Delta_{l}}$, we can
conclude that
$\mathbb{E}[M_{t+k_{1}\Delta_{l}}(1)]^{2}\leq\mathsf{C}\Delta^{2}_{l}$.
Therefore by combining all terms and a further application of the
$C_{2}$-inequality three times

$\displaystyle\mathbb{E}\Big{[}\Big{(}[\overline{U}_{t+k_{1}\Delta_{l}}(Y)-\overline{U}_{t+k_{1}\Delta_{l}}^{l}(Y)]\Big{)}^{2}\Big{]}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}[M_{t+k_{1}\Delta_{l}}(1)]^{2}+\mathbb{E}[R_{t+k_{1}\Delta_{l}}(1)]^{2}+\mathbb{E}[R_{t+k_{1}\Delta_{l}}]^{2}$
$\displaystyle\leq$ $\displaystyle\mathsf{C}\Delta^{2}_{l}.$

∎

##  Appendix C Analysis for i.i.d. MLEnKBF NC Estimator

We now discuss the analysis, related to the variance, of both the NC estimator
using the EnKBF and the i.i.d. MLEnKBF. This will lead onto the proof of our
main result, presented as Theorem 3.1. We note that in our notations, we
extend the case of the discretized EnKBF, to the discretized MLEnKBF, by
adding superscripts $l$ as above. Specifically the analysis now considers the
i.i.d. couple particle system

$\displaystyle\zeta_{(k+1)\Delta_{l}}^{i,l}$
$\displaystyle=\zeta_{k\Delta_{l}}^{i,l}+A\zeta_{k\Delta_{l}}^{i,l}\Delta_{l}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l}}^{i}-\overline{W}_{k\Delta_{l}}^{i}]$
$\displaystyle+P^{N,l}_{k\Delta_{l}}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l}}-Y^{{i}}_{k\Delta_{l}}]-\Big{[}C\zeta_{k\Delta_{l}}^{i,l}\Delta_{l}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l}}^{i}-\overline{V}_{k\Delta_{l}}^{i}]\Big{]}\Big{)},$
$\displaystyle\zeta_{(k+1)\Delta_{l-1}}^{i,l-1}$
$\displaystyle=\zeta_{k\Delta_{l-1}}^{i,l-1}+A\zeta_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+Q^{1/2}[\overline{W}_{(k+1)\Delta_{l-1}}^{i}-\overline{W}_{k\Delta_{l-1}}^{i}]$
$\displaystyle+P^{N,l-1}_{k\Delta_{l-1}}C^{\top}R^{-1}\Big{(}[Y^{{i}}_{(k+1)\Delta_{l-1}}-Y^{{i}}_{k\Delta_{l-1}}]-\Big{[}C\zeta_{k\Delta_{l-1}}^{i,l-1}\Delta_{l-1}+R^{1/2}[\overline{V}_{(k+1)\Delta_{l-1}}^{i}-\overline{V}_{k\Delta_{l-1}}^{i}]\Big{]}\Big{)},$

within the NC estimator.

We will use the fact that the i.i.d. system coincides with the Kalman–Bucy
Diffusion $\overline{X}_{t}$, in the limit of $N\rightarrow\infty$. This
implies the mean and covariance are defined through the Kalman–Bucy filter and
the Ricatti equations, which allows us to use the results from Appendix B, for
the process $\zeta^{i}_{k\Delta_{l}}$.

###  C.1 MSE bound on EnKBF NC estimator

Here we use the notation: for a $d_{x}-$dimensional vector $x$ denote
$\|x\|_{2}=(\sum_{j=1}^{d_{x}}x(j)^{2})^{1/2}$.

######  Proposition C.1.

For any $T\in\mathbb{N}$ fixed and $t\in[0,T-1]$ there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,N,k_{1})\in\mathbb{N}_{0}\times\\{2,3,\dots\\}\times\\{0,1,\dots,\Delta_{l}^{-1}\\}$:

$\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-\overline{U}_{t+k_{1}\Delta_{l}}](Y)\Big{\|}_{2}^{2}\Big{]}\leq\mathsf{C}\Big{(}\frac{1}{N}+\Delta_{l}^{2}\Big{)}.$

###### Proof.

Using the $C_{2}-$inequality one has

$\displaystyle\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}](Y)\Big{\|}_{2}^{2}\Big{]}$
$\displaystyle\leq\mathsf{C}\Big{(}\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}](Y)\Big{\|}_{2}^{2}\Big{]}+\mathbb{E}\Big{[}\Big{\|}[{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}-{\overline{U}}_{t+k_{1}\Delta_{l}}](Y)\Big{\|}_{2}^{2}\Big{]}\Big{)}.$
(C.1)

The first term on the R.H.S. can be controlled by standard results for i.i.d.
sampling (recall that
$\zeta_{t+k_{1}\Delta_{l}}^{i}|\mathscr{F}_{t+k_{1}\Delta_{l}}$ are i.i.d.
Gaussian with mean $m_{t+k_{1}\Delta_{l}}$ and covariance
$P_{t+k_{1}\Delta_{l}}$), that is

$\mathbb{E}\Big{[}\Big{\|}[{\widehat{\overline{U}}}_{t+k_{1}\Delta_{l}}^{N,l}-{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}](Y)\Big{\|}_{2}^{2}\Big{]}\leq\frac{\mathsf{C}}{N}.$
(C.2)

The formula in (C.2) can be proved by using the formulae for the NC
estimators, in the usual integral form, and through a simple application of
the general Minkowski inequality.
Note that it is crucial that (A.1) holds, otherwise the upper-bound can
explode as a function of $l$. For the right-most term on the R.H.S. of (C.1)
by Jensen’s inequality and Lemma B.4:

$\mathbb{E}\Big{[}\Big{\|}[{\overline{U}}_{t+k_{1}\Delta_{l}}^{l}-{\overline{U}}_{t+k_{1}\Delta_{l}}](Y)\Big{\|}_{2}^{2}\Big{]}\leq\mathsf{C}\Delta_{l}^{2}.$
(C.3)

So the proof can be concluded by combining (C.1), (C.2) and (C.3).

∎

###  C.2 Variance of i.i.d. MLEnKBF NC estimator

######  Proposition C.2.

For any $(t,q)\in\mathbb{N}_{0}\times[1,\infty)$, there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,N,k_{1})\in\mathbb{N}\times\\{2,3,\dots\\}\times\\{0,1,\dots,\Delta_{l-1}^{-1}-1\\}$:

$\mathbb{E}\Big{[}\Big{|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{N,l}-\widehat{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{N,l-1}]-[{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{l}-{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{l-1}]\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta^{1/2}_{l}}{\sqrt{N}}.$
(C.4)

As before, in order to proceed we will make use of a martingale-remainder
decomposition.
Recall that, for level $s\in\\{l-1,l\\}$, we have

$\displaystyle\overline{U}_{t+k_{1}\Delta_{l-1}}^{N,s}(Y)-\overline{U}_{t+k_{1}\Delta_{l-1}}^{s}(Y)$
$\displaystyle=M_{t+k_{1}\Delta_{l-1}}^{s}(Y)+R_{t+k_{1}\Delta_{l-1}}^{s}$
$\displaystyle=M^{s}_{t+k_{1}\Delta_{l-1}}(1)+R^{s}_{t+k_{1}\Delta_{l-1}}(1)+R^{s}_{t+k_{1}\Delta_{l-1}},$

Therefore substituting into the LHS of (C.4), and using Minkowski’s
inequality, results in

$\displaystyle\mathbb{E}\Big{[}\Big{|}[\widehat{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{N,l}-\widehat{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{N,l-1}]-[{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{l}-{\overline{U}}_{t+k_{1}\Delta_{l-1}}^{l-1}]\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=\mathbb{E}\Big{[}\Big{|}(R_{t+k_{1}\Delta_{l-1}}^{l}(1)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(1))+(M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1))$
$\displaystyle+(R_{t+k_{1}\Delta_{l-1}}^{l}(Y)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(Y))\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq\mathbb{E}\Big{[}\Big{|}(R_{t+k_{1}\Delta_{l-1}}^{l}(1)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(1))\Big{|}^{q}\Big{]}^{1/q}+\mathbb{E}\Big{[}\Big{|}(M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1))\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle+\mathbb{E}\Big{[}\Big{|}(R_{t+k_{1}\Delta_{l-1}}^{l}(Y)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(Y))\Big{|}^{q}\Big{]}^{1/q}.$

Therefore in order to prove Proposition C.2 we will split it into three
lemmas, which are stated and proved below.

######  Lemma C.1.

For any $(t,q)\in\mathbb{N}_{0}\times[1,\infty)$, there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,N,k_{1})\in\mathbb{N}\times\\{2,3,\dots\\}\times\\{0,1,\dots,\Delta_{l-1}^{-1}-1\\}$:

$\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l-1}}^{l}(1)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta^{1/2}_{l}}{\sqrt{N}}.$
(C.5)

###### Proof.

$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l-1}}^{l}(1)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum^{d_{y}}_{j_{1}=1}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}C(j_{1},j_{2})\Big{(}m^{N,l}_{k\Delta_{l-1}}(j_{2})-m^{l}_{k\Delta_{l-1}}(j_{2})\Big{)}R^{-1}(j_{1},j_{3})CX_{k\Delta_{l-1}}(j_{3})\Delta_{l}$
$\displaystyle-$
$\displaystyle\sum^{d_{y}}_{j_{1}=1}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}C(j_{1},j_{2})\Big{(}m^{N,l-1}_{k\Delta_{l-1}}(j_{2})-m^{l-1}_{k\Delta_{l-1}}(j_{2})\Big{)}R^{-1}(j_{1},j_{3})CX_{k\Delta_{l-1}}(j_{3})\Delta_{l-1}\Big{|}^{q}\Big{]}^{1/q}.$

Then through the generalized Minkowski inequality

$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t}^{l}(1)-R_{t}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum_{j_{4}=1}^{d_{x}}C(j_{1},j_{2})C(j_{3},j_{4})R^{-1}(j_{1},j_{3})\times$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})X_{\tau_{s}^{l}}(j_{4})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})X_{\tau_{s}^{l-1}}(j_{4})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\Big{(}{\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}X_{\tau^{l}_{s}}(j_{4})ds\Big{|}^{q}\Big{]}^{1/q}}$
$\displaystyle+$
$\displaystyle{\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{(}X_{\tau^{l}_{s}}(j_{4})-X_{\tau^{l-1}_{s}}(j_{4})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}}\Big{)}$
$\displaystyle=:$ $\displaystyle T_{1}+T_{2},$

where we have used $\tau^{s}_{t}=\lfloor\frac{t}{\Delta_{s}}\rfloor\Delta_{s}$
for $t\in\mathbb{R}^{+}$, and

$\overline{m}^{N,s}_{k\Delta_{l}}={m}^{N,s}_{k\Delta_{l}}-m^{l}_{k\Delta_{s}},\quad
s\in\\{l-1,l\\}.$ (C.6)

For $T_{1}$ we can express it as

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})+\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}X_{\tau^{l}_{s}}(j_{4})ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\mathbb{E}\Big{[}\Big{|}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}X_{\tau^{l}_{s}}(j_{4})\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle+$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}X_{\tau^{l}_{s}}(j_{4})\Big{|}^{q}\Big{]}^{1/q}\Big{)}ds,$
$\displaystyle=:$ $\displaystyle T_{3}+T_{4}.$

again by using the generalized Minkowski inequality and Jensen’s inequality.

For $T_{3}$ we can apply the Marcinkiewicz–Zygmund and Hölder inequalities,
and using the fact that means can be expressed as the expectations of (3.9) -
(3.10)

$\displaystyle T_{3}$ $\displaystyle=$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l-1}}_{0}\mathbb{E}\Big{[}\Big{|}\frac{1}{N}\sum^{N}_{i=1}\Big{[}\Big{(}{\zeta}^{i,l}_{\tau^{l}_{s}}(j_{2})-{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}-\Big{(}m^{l}_{\tau^{l}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}\Big{]}X_{\tau^{l}_{s}}(j_{4})\Big{|}^{q}\Big{]}^{1/q}ds$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l-1}}_{0}\mathbb{E}\Big{[}\Big{|}\frac{1}{N}\sum^{N}_{i=1}\Big{[}\Big{(}{\zeta}^{i,l}_{\tau^{l}_{s}}(j_{2})-{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}-\Big{(}m^{l}_{\tau^{l}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}\Big{]}\Big{|}^{2q}\Big{]}^{1/2q}\times$
$\displaystyle\mathbb{E}\Big{[}\Big{|}X_{\tau^{l}_{s}}(j_{4})\Big{|}^{2q}\Big{]}^{1/2q}ds$
$\displaystyle\leq$
$\displaystyle\frac{\mathsf{C}_{q}}{\sqrt{N}}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{4}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l-1}}_{0}\mathbb{E}\Big{[}\Big{|}[{\zeta}^{l}_{\tau^{l}_{s}}(j_{2})-{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{2})]\Big{|}^{2q}\Big{]}^{1/2q}\mathbb{E}\Big{[}\Big{|}X_{\tau^{l}_{s}}(j_{4})\Big{|}^{2q}\Big{]}^{1/2q}ds.$

The process $\mathbb{E}[|X_{\tau^{l}_{s}}]^{2q}]$ is of order $\mathcal{O}(1)$
and and the recursion is of order $\mathcal{O}({\Delta_{l}})$, through the
strong error Lemma B.2. For $T_{4}$ we know it is sufficiently small, which is
of order $\mathcal{O}(\Delta_{l})$.

For $T_{2}$, we use the definition of the discretized diffusion process

$X_{\tau^{l}_{s}}-X_{\tau^{l-1}_{s}}=\int^{\tau^{l}_{s}}_{\tau^{l-1}_{s}}AX_{u}du+Q^{1/2}\Big{(}W_{\tau^{l}_{s}}-W_{\tau^{l-1}_{s}}\Big{)}.$
(C.7)

As before, we know the difference of the Brownian motion increment
$\mathbb{E}\Big{[}\Big{|}W_{\tau^{l}_{s}}-W_{\tau^{l-1}_{s}}\Big{|}^{q}\Big{]}$
is of order $\mathcal{O}(\Delta^{1/2}_{l})$, and, as before,
$\mathbb{E}[|X_{u}|^{q}]^{1/q}\leq\mathsf{C}$. Therefore the integral term of
(C.7) is of order $\mathcal{O}(\Delta^{1/2}_{l})$. Finally for
$\mathbb{E}[|\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}|^{q}]^{1/q}$, as it is of
order $\mathcal{O}(N^{-\frac{1}{2}})$, therefore, combining all terms, we can
deduce from that

$\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l-1}}^{l}(1)-R_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta^{1/2}_{l}}{\sqrt{N}}.$

∎

######  Lemma C.2.

For any $(t,q)\in\mathbb{N}_{0}\times[1,\infty)$, there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,N,k_{1})\in\mathbb{N}\times\\{2,3,\dots\\}\times\\{0,1,\dots,\Delta_{l-1}^{-1}-1\\}$:

$\mathbb{E}\Big{[}\Big{|}\Big{(}M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta_{l}^{1/2}}{\sqrt{N}}.$
(C.8)

###### Proof.

As before, we set $\tau^{s}_{t}=\lfloor\frac{t}{\Delta_{s}}\rfloor\Delta_{s}$
for $t\in\mathbb{R}^{+}$, and make use of (C.6)

$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum^{d_{y}}_{j_{1}=1}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\Big{(}\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}C(j_{1},j_{2})\Big{(}\overline{m}^{N,l}_{k\Delta_{l-1}}(j_{2})\Big{)}R^{-1}(j_{1},j_{3})\times$
$\displaystyle\Big{(}[Y_{(k+1)\Delta_{l-1}}-Y_{k\Delta_{l-1}}](j_{3})-CX_{k\Delta_{l-1}}(j_{3})\Delta_{l}\Big{)}$
$\displaystyle-$
$\displaystyle\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}C(j_{1},j_{2})\Big{(}\overline{m}^{N,l-1}_{k\Delta_{l-1}}(j_{2})\Big{)}R^{-1}(j_{1},j_{3})\Big{(}[Y_{(k+1)\Delta_{l-1}}-Y_{k\Delta_{l-1}}](j_{3})-CX_{k\Delta_{l-1}}(j_{3})\Delta_{l-1}\Big{)}\Big{)}\Big{|}^{q}\Big{]}^{1/q}$

Then by using generalized Minkowski and Jensen’s inequality

$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\sum_{j_{1}=1}^{d_{y}}\sum_{j_{2}=1}^{d_{x}}\sum^{d_{y}}_{j_{3}=1}C(j_{1},j_{2})R^{-1}(j_{1},j_{3})\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}dY_{s}(j_{3})$
$\displaystyle-$
$\displaystyle\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})CX_{\tau_{s}^{l}}(j_{3})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})CX_{\tau_{s}^{l-1}}(j_{3})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum_{j_{4}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}X_{\tau^{l}_{s}}(j_{4})ds+\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}dV_{s}(j_{3})$
$\displaystyle-$
$\displaystyle\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})X_{\tau_{s}^{l}}(j_{4})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})X_{\tau_{s}^{l-1}}(j_{4})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{2}=1}^{d_{x}}\sum_{j_{3}=1}^{d_{y}}\sum_{j_{4}=1}^{d_{x}}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\mathbb{E}\Big{[}\Big{|}-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{(}X_{\tau_{s}^{l}}(j_{4})-X_{\tau_{s}^{l-1}}(j_{4})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle+$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}dV_{s}(j_{3})\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle=:$ $\displaystyle T_{2}+T_{3}.$

For $T_{2}$, we know it follows the same analysis as Lemma C.1. Therefore we
can conclude that

$T_{2}\leq\frac{\mathsf{C}\Delta^{1/2}_{l}}{\sqrt{N}}.$

For $T_{3}$ by using the Burkholder–Davis–Gundy and Hölder inequality, as done
previously from (A.5), and we use the bound from Lemma C.1, for the
$\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}-\overline{m}^{N,l-1}_{\tau^{l-1}_{s}}\Big{)}$
term. Therefore this implies that

$\mathbb{E}\Big{[}\Big{|}\Big{(}M_{t+k_{1}\Delta_{l-1}}^{l}(1)-M_{t+k_{1}\Delta_{l-1}}^{l-1}(1)\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta_{l}^{1/2}}{\sqrt{N}}.$

∎

######  Lemma C.3.

For any $(t,q)\in\mathbb{N}_{0}\times[1,\infty)$, there exists a
$\mathsf{C}<+\infty$ such that for any
$(l,N,k_{1})\in\mathbb{N}\times\\{2,3,\dots\\}\times\\{0,1,\dots,\Delta_{l-1}^{-1}-1\\}$:

$\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l-1}}^{l}-R_{t+k_{1}\Delta_{l-1}}^{l-1}\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta_{l}}{\sqrt{N}}.$
(C.9)

###### Proof.

Again we let $\tau^{s}_{t}=\lfloor\frac{t}{\Delta_{s}}\rfloor\Delta_{s}$ for
$t\in\mathbb{R}^{+}$, and using the generalized Minkowski inequality

$\displaystyle\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l}}^{l}-R_{t+k_{1}\Delta_{l}}^{l-1}\Big{)}\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=$
$\displaystyle\Big{(}\mathbb{E}\Big{[}\Big{|}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\frac{S(j_{1},j_{2})}{2}\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}\Big{(}m^{N,l}_{k\Delta_{l-1}}(j_{1})m^{N,l}_{k\Delta_{l-1}}(j_{2})-m^{l}_{k\Delta_{l-1}}(j_{1})m^{l}_{k\Delta_{l-1}}(j_{2})\Big{)}$
$\displaystyle-$ $\displaystyle\sum^{t\Delta^{-1}_{l-1}+k_{1}-1}_{k=0}\
\Big{(}m^{N,l-1}_{k\Delta_{l-1}}(j_{1})m^{N,l-1}_{k\Delta_{l-1}}(j_{2})-m^{l}_{k\Delta_{l-1}}(j_{1})m^{l}_{k\Delta_{l-1}}(j_{2})\Big{)}\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle=$
$\displaystyle\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\frac{S(j_{1},j_{2})}{2}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}m^{N,l}_{\tau^{l}_{s}}(j_{1})m^{N,l}_{\tau^{l}_{s}}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{1})m^{l}_{\tau^{l}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}m^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle+$
$\displaystyle\Big{(}m^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}m^{N,l-1}_{\tau^{l-1}_{s}}(j_{1})m^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}m^{N,l}_{\tau^{l}_{s}}(j_{1})m^{N,l}_{\tau^{l}_{s}}(j_{2})-m^{l}_{\tau^{l}_{s}}(j_{1})m^{l}_{\tau^{l}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}m^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle+$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\Big{(}m^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}m^{N,l-1}_{\tau^{l-1}_{s}}(j_{1})m^{N,l-1}_{\tau^{l-1}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})m^{l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle=:$ $\displaystyle T_{1}+T_{2}.$

For $T_{1}$ we make use of (C.6) and with the generalized Minkowski inequality

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}m^{N,l}_{\tau^{l}_{s}}(j_{1})m^{N,l}_{\tau^{l}_{s}}(j_{2})-m_{\tau^{l}_{s}}(j_{1})m_{\tau^{l}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}m^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})-m_{\tau^{l-1}_{s}}(j_{1})m_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle=$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{1})m^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle+$
$\displaystyle\Big{(}m_{\tau^{l}_{s}}(j_{1})\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-m_{\tau^{l-1}_{s}}(j_{1})\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{1})m^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{1})m^{N,l}_{\tau^{l-1}_{s}}(j_{2})ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle+$
$\displaystyle\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}m^{l}_{\tau^{l}_{s}}(j_{1})\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-m^{l}_{\tau^{l-1}_{s}}(j_{1})\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})ds\Big{|}^{q}\Big{]}^{1/q}$
$\displaystyle=:$ $\displaystyle T_{3}+T_{4}.$

For $T_{3}$ we can use the difference of mean trick, as in (B.6),

$T_{3}=\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{1})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}m^{N,l}_{\tau^{l}_{s}}(j_{2})+\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{1})\Big{(}{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q}.$

We know $\mathbb{E}[|m^{N,l}_{\tau^{l}_{s}}|^{q}]\leq\mathsf{C}$, and
$\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{1})$ is of oder
$\mathcal{O}(N^{-\frac{1}{2}})$. The last bracket term is of order
$\mathcal{O}(\Delta_{l})$, arising from the strong error. The first bracket
term is the same, that appears in Lemma C.1, which is of order
$\mathcal{O}(\frac{\Delta_{l}}{\sqrt{N}})$.

Similarly, $T_{4}$ can be expressed as

$T_{4}=\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\overline{m}^{N,l}_{\tau^{l}_{s}}(j_{2})-\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}m^{l}_{\tau^{l}_{s}}(j_{1})+\overline{m}^{N,l}_{\tau^{l-1}_{s}}(j_{2})\Big{(}{m}^{l}_{\tau^{l}_{s}}(j_{1})-{m}^{l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}ds\Big{|}^{q}\Big{]}^{1/q},$

which contains the same bounds as $T_{3}$, therefore implying that

$T_{1}\leq\frac{\mathsf{C}\Delta_{l}}{\sqrt{N}}.$

For $T_{2}$, we can rewrite it as

$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}$
$\displaystyle-$
$\displaystyle\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{1})]\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{2})]+\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{1})]\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{2})]ds\Big{|}^{q}\Big{]}^{1/q}\Big{)}.$

Now to proceed, we again use
$\overline{m}^{N,s}_{\tau^{l-1}_{s}}={m}^{N,s}_{\tau^{l-1}_{s}}-m_{\tau^{l-1}_{s}}$,
for $s\in\\{l-1,l\\}$, which we can rewrite, in terms of its expectation, as

$\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j)=\frac{1}{N}\sum^{N}_{i=1}{{\zeta}}^{i,s}_{\tau^{l-1}_{s}}(j)-\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j)].$

Therefore the integral terms of $T_{2}$, for $s\in\\{l-1,l\\}$, become

$\displaystyle\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{2})\Big{)}-\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{1})]\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{2})]$
$\displaystyle=\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{{\zeta}}^{i,s}_{\tau^{l-1}_{s}}(j_{1})+\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{1})]\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{2})+\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{2})]\Big{)}-\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{1})]\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{2})]$
$\displaystyle=\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{2})\Big{)}+\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,s}(j_{1})\Big{)}\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{2})]+\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,s}_{\tau^{l-1}_{s}}(j_{2})\Big{)}\mathbb{E}[{\zeta}^{s}_{\tau^{l-1}_{s}}(j_{1})].$

By substituting this into $T_{2}$, we have

$T_{2}=\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}T_{5}+T_{6}+T_{7}\
ds\Big{|}^{q}\Big{]}^{1/q}\Big{)},$

where

$\displaystyle T_{5}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}-\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)},$
$\displaystyle T_{6}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{2})]-\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{2})],$
$\displaystyle T_{7}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})\Big{)}\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{1})]-\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{1})].$

For $T_{6}$ and $T_{7}$, they can be expressed as

$\displaystyle T_{6}+T_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})\Big{(}\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{2})]-\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{2})]\Big{)}+\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{2})]\frac{1}{N}\sum^{n}_{i=1}\Big{(}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})-\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l-1}(j_{1})\Big{)}$
$\displaystyle+$
$\displaystyle\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})\Big{(}\mathbb{E}[{\zeta}^{l}_{\tau^{l-1}_{s}}(j_{1})]-\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{1})]\Big{)}+\mathbb{E}[{\zeta}^{l-1}_{\tau^{l-1}_{s}}(j_{1})]\frac{1}{N}\sum^{n}_{i=1}\Big{(}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})-\overline{\zeta}^{i,l-1}(j_{2})\Big{)}.$

For $T_{5}$, we can rewrite it as

$\displaystyle
T_{5}=\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})\Big{(}\frac{1}{N}\sum^{N}_{i=1}\Big{(}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})-\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\Big{)}+\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{1})\Big{)}\frac{1}{N}\sum^{N}_{i=1}\Big{(}\overline{\zeta}^{i,l}_{\tau^{l-1}_{s}}(j_{2})-\overline{\zeta}^{i,l-1}_{\tau^{l-1}_{s}}(j_{2})\Big{)}.$

Using these expressions, and through the generalized Minkowski and Jensen’s
inequality, $T_{2}$ can be simplified as

$\displaystyle T_{2}$ $\displaystyle\leq$
$\displaystyle\mathsf{C}\sum_{j_{1}=1}^{d_{x}}\sum_{j_{2}=1}^{d_{x}}\Big{(}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})\Big{(}\frac{1}{N}\sum^{N}_{i=1}\Big{(}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{1})-\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l-1}(j_{1})\Big{)}\Big{)}ds\Big{|}^{p}\Big{]}^{1/p}$
$\displaystyle+$
$\displaystyle\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\Big{(}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{1})\Big{(}\frac{1}{N}\sum^{N}_{i=1}\Big{(}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})-\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l-1}(j_{2})\Big{)}\Big{)}ds\Big{|}^{p}\Big{]}^{1/p}$
$\displaystyle+$
$\displaystyle\Big{|}\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l}(j_{2})]-\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l-1}(j_{2})]\Big{|}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{1})ds\Big{|}^{p}\Big{]}^{1/p}$
$\displaystyle+$
$\displaystyle\Big{|}\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l-1}(j_{2})]\Big{|}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\frac{1}{N}\sum^{n}_{i=1}\Big{(}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{1})-\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l-1}(j_{1})\Big{)}ds\Big{|}^{p}\Big{]}^{1/p}$
$\displaystyle+$
$\displaystyle\Big{|}\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l}(j_{1})]-\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l-1}(j_{1})]\Big{|}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\frac{1}{N}\sum^{N}_{i=1}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})ds\Big{|}^{p}\Big{]}^{1/p}$
$\displaystyle+$
$\displaystyle\Big{|}\mathbb{E}[\zeta_{\tau^{l-1}_{s}}^{l-1}(j_{1})]\Big{|}\mathbb{E}\Big{[}\Big{|}\int^{t+k_{1}\Delta_{l-1}}_{0}\frac{1}{N}\sum^{n}_{i=1}\Big{(}\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l}(j_{2})-\overline{\zeta}_{\tau^{l-1}_{s}}^{i,l-1}(j_{2})\Big{)}ds\Big{|}^{p}\Big{]}^{1/p}\Big{)}.$

Then by using the results from the strong and weak error of the diffusion
processes, in Lemma B.2 and B.3, we reach the following bound of

$\mathbb{E}\Big{[}\Big{|}\Big{(}R_{t+k_{1}\Delta_{l-1}}^{l}-R_{t+k_{1}\Delta_{l-1}}^{l-1}\Big{)}\Big{|}^{q}\Big{]}^{1/q}\leq\frac{\mathsf{C}\Delta_{l}}{\sqrt{N}}.$

∎

Therefore by combining all the results from Lemma C.1 \- Lemma C.3, leads to
the desired result of (C.4) in Proposition C.2.

###  C.3 Proof of Theorem 3.1

###### Proof.

Noting (2.15) one has

$[\widehat{\overline{U}}_{t}^{ML}-{\overline{U}}_{t}](Y)=[\widehat{\overline{U}}_{t}^{N_{0},0}-{\overline{U}}_{t}^{0}](Y)+\sum_{l=1}^{L}[\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}-{\overline{U}}_{t}^{l}+{\overline{U}}_{t}^{l-1}](Y)+[{\overline{U}}_{t}^{L}-{\overline{U}}_{t}](Y).$

Thus, by using three applications of the $C_{2}-$inequality we have

$\displaystyle\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t}^{ML}-{\overline{U}}_{t}](Y)\Big{\|}_{2}^{2}\Big{]}$
$\displaystyle\leq$
$\displaystyle\mathsf{C}\Big{(}\mathbb{E}\Big{[}\Big{\|}[\widehat{\overline{U}}_{t}^{N_{0},0}-{\overline{U}}_{t}^{0}](Y)\Big{\|}_{2}^{2}\Big{]}$
$\displaystyle+$
$\displaystyle\mathbb{E}\Big{[}\Big{\|}\sum_{l=1}^{L}[\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}-{\overline{U}}_{t}^{l}+{\overline{U}}_{t}^{l-1}](Y)\Big{\|}_{2}^{2}\Big{]}+\mathbb{E}\Big{[}\Big{\|}[{\overline{U}}_{t}^{L}-{\overline{U}}_{t}](Y)\Big{\|}_{2}^{2}\Big{]}\Big{)}.$

For the first term on the R.H.S. one can use (C.2) and for the last term on
the R.H.S., we can use (C.3). For the middle term, one has

$\displaystyle\Big{\|}\sum_{l=1}^{L}[\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}-{\overline{U}}_{t}^{l}+{\overline{U}}_{t}^{l-1}](Y)\Big{\|}_{2}^{2}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{L}\sum_{j=1}^{d_{x}}\Big{(}[\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}-{\overline{U}}_{t}^{l}+{\overline{U}}_{t}^{l-1}](Y)^{2}\Big{)}(j)^{2}$
$\displaystyle+$
$\displaystyle\sum_{l=1}^{L}\sum_{q=1}^{L}\mathbb{I}_{D^{c}}(l,q)\sum_{j=1}^{d_{x}}\Big{(}[\widehat{\overline{U}}_{t}^{N_{l},l}-\widehat{\overline{U}}_{t}^{N_{l},l-1}-{\overline{U}}_{t}^{l}+{\overline{U}}_{t}^{l-1}](Y)\Big{)}(j)\times$
$\displaystyle\Big{(}[\widehat{\overline{U}}_{t}^{N_{q},q}-\widehat{\overline{U}}_{t}^{N_{q},q-1}-{\overline{U}}_{t}^{q}+{\overline{U}}_{t}^{q-1}](Y)\Big{)}(j).$

Then using a combination of the independence of the coupled particle systems
along with Proposition C.2 the proof can be concluded.

∎

## Acknowledgments

This work was supported by KAUST baseline funding.

## References

  * [1] A. Arapostathis and S. I. Marcus. Analysis of an identification algorithm arising in the adaptive estimation of Markov chains.  Math Control Signals Syst 3:1–29, 1990.
  * [2] A. Bain and D. Crisan.  Fundamentals of Stochastic Filtering.  Springer, New York, 2009.
  * [3] A. Beskos, D. Crisan, A. Jasra, N. Kantas and H. Ruzayqat. Score-based parameter estimation for a class of continuous-time state space models.  SIAM J. Sci. Comp. 43(4), A2555-A2580, 2021.
  * [4] R. Bhar.  Stochastic Filtering with Applications in Finance.  World Scientific, 2010.
  * [5] A. N. Bishop and P. Del Moral.  On the mathematical theory of ensemble (linear-gaussian) Kalman–Bucy filtering. arXiv preprint arXiv:2006.08843, 2020.
  * [6] A. N. Bishop and P. Del Moral. On the stability of Kalman–Bucy diffusion processes. SIAM J. Control Optim., 55(6), 4015-4047, 2017.
  * [7] N. Brosse, A. Durmus and E. Moulines. Normalizing constants of log-concave densities.  Elect. J. Stat.,12, 851–889, 2018.
  * [8] F. Cerou, P. Del Moral and A. Guyader. A non-asymptotic variance theorem for unnormalized Feynman–Kac particle models.  Ann. Inst. Henri Poincare, 4, 629–649, 2011.
  * [9] N. K. Chada, J. Franks, A. Jasra, K. J. H. Law and M. Vihola. Unbiased inference for discretely observed hidden Markov model diffusions.  SIAM/ASA J. Uncertainty Quantification, 9(2), 763–787, 2021.
  * [10] N. K. Chada, A. Jasra and F. Yu.  Multilevel ensemble Kalman–Bucy filters.  SIAM/ASA J. Uncertainty Quantification, (to appear), 2022.
  * [11] A. Chernov, H. Hoel, K. J. H. Law, F. Nobile and R. Tempone.  Multilevel ensemble Kalman filtering for spatio-temporal processes.  Numer. Math., 1–55, 2021.
  * [12] D. Crisan, P. Del Moral, A. Jasra and H. Ruzayqat. Log-normalization constant estimation using the ensemble Kalman-Bucy filter with application to high-dimensional models.  Adv. Appl. Probab., (to appear), 2022.
  * [13] D. Crisan and B. Rozovskii.  The Oxford Handbook of Nonlinear Filtering. Oxford University Press, Oxford, 2011.
  * [14] P. Del Moral, A. Jasra, K. Law, and Y. Zhou. Multilevel sequential Monte Carlo samplers for normalizing constants.  TOMACS, 27, 2017.
  * [15] P. Del Moral, A. Doucet and S. S. Singh. Forward smoothing using sequential Monte Carlo. arXiv preprint arxiv:1012.5390, 2010.
  * [16] P. Del Moral. Feynman-Kac Formulae. Springer-Verlag, New York, 2004.
  * [17] P Del Moral, and J. Tugaut. On the stability and the uniform propagation of chaos properties of ensemble Kalman–Bucy filters. Ann. Appl. Probab., 28, 790–850, 2018.
  * [18] G. Evensen. Data Assimilation: The Ensemble Kalman Filter. Springer Science & Business Media, 2nd edition, 2009.
  * [19] K. Fossum, T. Mannseth and A. S. Stordal.  Assessment of multilevel ensemble-based data assimilation for reservoir history matching.  Computational Geosciences, 24, 217–239, 2020.
  * [20] M. B. Giles. Multilevel Monte Carlo path simulation.  Op. Res., 56, 607–617, 2008.
  * [21] M. B. Giles.  Multilevel Monte Carlo methods.  Acta Numerica, 24, 259–328, 2015.
  * [22] A. Gelman and X.-L. Meng. Simulating normalizing constants: From importance sampling to bridge sampling to path sampling.  Statistical science, 163–185, 1998.
  * [23] H. Hoel, K. J. H. Law and R. Tempone.  Multilevel ensemble Kalman filtering.  SIAM J. Numer. Anal., 54(3), 1813–1839, 2016.
  * [24] H. Hoel, G. Shaimerdenova and R. Tempone. Multi-index ensemble Kalman filtering  arXiv preprint arXiv:2104.07263, 2021.
  * [25] A. Jasra, K. Kamatani, K. J. H. Law and Y. Zhou.  Multilevel particle filters.  SIAM J. Numer. Anal., 55(6), 3068–3096, 2017.
  * [26] A. Jasra, K. Kamatani, P. P. Osei and Y. Zhou. Multilevel particle filters: normalizing constant estimation.  Stat Comput, 28:47–60, 2018.
  * [27] A. Jazwinski.  Stochastic processes and filtering theory,  vol.x 63. Academic Pr, 1970.
  * [28] S. Kostov and N. Whiteley. An algorithm for approximating the second moment of the normalizing constant estimate from a particle filter.  Methodol Comput Appl Probab, 19:799–818, 2017.
  * [29] F. Le Gland and M. Mevel. Recursive identification in hidden Markov models. Proc. 36th IEEE Conf. Dec. Contr., 3468–3473, 1997.
  * [30] E. N. Lorenz.  Deterministic nonperiodic flow.  Journal of the Atmospheric Sciences, 20(2): 130–141, 1963.
  * [31] E. N. Lorenz.  Predictability: A problem partly solved.  Proc. ECMWF Seminar on predictability, 1, 1–18, 1996.
  * [32] A. Majda and X. Wang.  Non-linear Dynamics and Statistical Theories for Basic Geophysical Flows, Cambridge University Press,  2006\.
  * [33] D. Oliver, A. C. Reynolds and N. Liu.  Inverse Theory for Petroleum Reservoir Characterization and History Matching, Cambridge University Press, 1st edn,  2008\.
  * [34] S. Reich. Data assimilation: the Schrödinger perspective.  Acta Numerica, 28, 635–711, 2019.
  * [35] G. Poyiadjis, A. Doucet and S. S. Singh. Particle approximations of the score and observed information matrix in state space models with application to parameter estimation. Biometrika 98: 65–80, 2011.
  * [36] P. Sakov and P.R. Oke. A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters.  Tellus A., 60(2), 361–371, 2008.
  * [37] J. C. Spall. Introduction to stochastic search and optimization, 1st edn. Wiley, New York, 2003.
  * [38] J. C. Spall. Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.  IEEE Trans. Auto. Control., 37, 332–341, 1992.
  * [39] M. Rischard, P. E. Jacob and N. Pillai. Unbiased estimation of log normalizing constants with applications to Bayesian cross-validation. arXiv preprint arxiv:1810.01382, 2018.
  * [40] M. Vihola. Unbiased estimators and multilevel Monte Carlo.  Op. Res., 66, 448–462, 2018.