Sebastian Reich 11institutetext: University of Potsdam, Institute of
Mathematics, Potsdam, Germany
11email<EMAIL_ADDRESS>
# Robust parameter estimation using the ensemble Kalman filter

Sebastian Reich 11

###### Abstract

Standard maximum likelihood or Bayesian approaches to parameter estimation for
stochastic differential equations are not robust to perturbations in the
continuous-in-time data. In this paper, we give a rather elementary
explanation of this observation in the context of continuous-time parameter
estimation using an ensemble Kalman filter. We employ the frequentist
perspective to shed new light on three robust estimation techniques; namely
subsampling the data, rough path corrections, and data filtering. We
illustrate our findings through a simple numerical experiment.

###### keywords:

parameter estimation, stochastic differential equations, ensemble Kalman
filter, frequentist approach, rough path theory

##  1 Introduction

In this note, which is an extended version of [30], we consider the well-
studied problem of parameter estimation for stochastic differential equations
(SDEs) from continuous-time observations $X_{t}^{\dagger}$, $t\in[0,T]$ [26].
It is well-known that the corresponding maximum likelihood estimator does not
depend continuously on the observations $X_{t}^{\dagger}$, $t\in[0,T]$, which
can result in a systematic estimation bias [28, 15]. In other words, the
maximum likelihood estimator is not robust with respect to perturbations in
the observations. Here, we revisit this problem from the perspective of online
(time-continuous) parameter estimation [7, 12] using the popular ensemble
Kalman filter (EnKF) and its continuous-time ensemble Kalman–Bucy filter
(EnKBF) formulations [16, 11, 27]. As for the corresponding maximum likelihood
approaches, the EnKBF does not depend continuously on the incoming
observations $X_{t}^{\dagger}$, $t\geq 0$, with respect to the uniform norm
topology on the space of continuous functions. This fact has been first
investigated in [10] using rough path theory [17]. In particular, as already
demonstrated for the related maximum likelihood estimator in [15], rough path
theory allows one to specify an appropriately generalised topology which leads
to a continuous dependence of the EnKBF estimators on the observations. Here
we expand the analysis of [10] to a frequentist analysis of the EnKBF in the
spirit of [31], where the primary focus is on the expected behaviour of the
EnKBF estimators over all admissible observation paths. One recovers that the
discontinuous dependence of the EnKBF estimators on the driving observations
results in a systematic bias from a frequentist perspective. This is also a
well known fact for SDEs driven by multiplicative noise [24].

The proposed frequentist perspective naturally enables the study of known bias
correction methods, such as subsampling the data [28], a recently proposed
data filtering approach [2], as well as novel de-biasing approaches [10] in
the context of the EnKBF.

In order to facilitate a rather elementary mathematical analysis, we consider
only the very much simplified problem of parameter estimation for linear SDEs.
This restriction allows us to avoid certain technicalities from rough path
theory and enables a rather straightforward application of the numerical rough
path approach put forward in [14]. As a result we are able to demonstrate that
the popular approach of subsampling the data [3, 28, 6] can be well justified
from a frequentist perspective. The frequentist perspective also suggests a
rather natural approach to the estimation of the required correction term in
the case an EnKBF is implemented without subsampling.

We end this introductory paragraph with a reference to [1], which includes a
broad survey on alternative estimation techniques. We also point to [10] for
an in-depth discussion of rough path theory in connection to filtering and
parameter estimation.

The remainder of this paper is structured as follows. The problem setting and
the EnKBF are introduced in the subsequent Section 2. The frequentist
perspective and its implications on the specific implementations of an EnKBF
in the context of low and high frequency data assimilation are laid out in
Section 3. The importance of these considerations becomes transparent when
applying the EnKBF to perturbed data in Section 4. Here again, we restrict
attention to a rather simple model setting taken from [18] and also used in
[10]. As a result we build a clear connection between subsampling and the
necessity for a correction term in the case high frequency data is assimilated
directly. We also provide a discussion of the data filtering approach [2] in
the context of our simply model system. A brief numerical demonstration is
provided in Section 5, which is followed by a concluding remark in Section 6.

##  2 Ensemble Kalman parameter estimation

We consider the SDE parameter estimation problem

${\rm d}X_{t}=f(X_{t},\theta){\rm d}t+\gamma^{1/2}{\rm d}W_{t}$ (1)

subject to observations $X_{t}^{\dagger}$, $t\in[0,T]$, which arise from the
reference system

${\rm d}X_{t}^{\dagger}=f^{\dagger}(X_{t}^{\dagger}){\rm d}t+\gamma^{1/2}{\rm
d}W^{\dagger}_{t},$ (2)

where the unknown drift function $f^{\dagger}(x)$ typically satisfies
$f^{\dagger}(x)=f(x,\theta^{\dagger})$ and $\theta^{\dagger}$ denotes the true
parameter value. Here we assume for simplicity that the unknown parameter is
scalar-valued and that the state variable is $d$-dimensional with $d\geq 1$.
Furthermore, $W_{t}$ and $W_{t}^{\dagger}$ denote independent standard
$d$-dimensional Brownian motions and $\gamma>0$ is the (known) diffusion
constant.

Following the Bayesian paradigm, we treat the unknown parameter as a random
variable $\Theta$. Furthermore, we apply a sequential approach and update
$\Theta$ with the incoming data $X^{\dagger}_{t}$ as a function of time. Hence
we introduce the random variable $\Theta_{t}$ which obeys the Bayesian
posterior distribution given all observations $X_{\tau}^{\dagger}$,
$\tau\in[0,t]$, up to time $t>0$. Furthermore, instead of exactly solving the
time-continuous Bayesian inference problem as specified by the associated
Kushner–Stratonovitch equation [7, 27], we define the time evolution of
$\Theta_{t}$ by an application of the (deterministic) ensemble Kalman–Bucy
filter (EnKBF) mean-field equations [11, 27], which take the form

${\rm d}\Theta_{t}&=\gamma^{-1}\pi_{t}\left[(\theta-\pi_{t}[\theta])\otimes
f(X_{t}^{\dagger},\theta)\right]{\rm d}I_{t},\\\ {\rm d}I_{t}={\rm
d}X_{t}^{\dagger}-\frac{1}{2}\left(f(X_{t}^{\dagger},\Theta_{t})+\pi_{t}[f(X^{\dagger}_{t},\theta)]\right){\rm
d}t,$ (3)

where $\pi_{t}$ denotes the probability density function (PDF) of $\Theta_{t}$
and $\pi_{t}[g]$ the associated expectation value of a function $g(\theta)$.
The column vector $I_{t}$, defined by (3b), is called the innovation, while
the row vector

$K_{t}(\pi_{t})=\gamma^{-1}\pi_{t}\left[(\theta-\pi_{t}[\theta])\otimes
f(X_{t}^{\dagger},\theta)\right],$ (4)

premultiplying the innovation in (3a) is called the gain. Here the notation
$a\otimes b=ab^{\rm T}$, where $a,b$ can be any two column vectors, has been
used. The initial condition $\Theta_{0}\sim\pi_{0}$ is provided by the prior
PDF of the unknown parameter.

A Monte-Carlo implementation of the mean-field equations (3) leads to the
interacting particle system

${\rm
d}\Theta_{t}^{(i)}&=\gamma^{-1}\pi_{t}^{M}\left[(\theta-\pi_{t}^{M}[\theta])\otimes
f(X_{t}^{\dagger},\theta)\right]{\rm d}I_{t}^{(i)},\\\ {\rm d}I_{t}^{(i)}={\rm
d}X_{t}^{\dagger}-\frac{1}{2}\left(f(X_{t}^{\dagger},\Theta_{t}^{(i)})+\pi_{t}^{M}[f(X^{\dagger}_{t},\theta)]\right){\rm
d}t,$ (5)

$i=1,\ldots,M$, where expectations are now taken with respect to the empirical
measure. That is,

$\pi_{t}^{M}[g]=\frac{1}{M}\sum_{i=1}^{M}g(\Theta_{t}^{(i)})$ (6)

for given function $g(\theta)$, and all Monte-Carlo samples are driven by the
same (fixed) observations $X_{t}^{\dagger}$. The initial samples
$\Theta_{0}^{(i)}$, $i=1,\ldots,M$, are drawn identically and independently
from the prior distribution $\pi_{0}$.

We note in passing that there is also a stochastic variant of the innovation
process [27] defined by

${\rm d}I_{t}={\rm d}X_{t}^{\dagger}-f(X_{t}^{\dagger},\Theta_{t}){\rm
d}t-\gamma^{1/2}{\rm d}W_{t},$ (7)

which leads to the Monte-Carlo approximation

${\rm d}I_{t}^{(i)}={\rm
d}X_{t}^{\dagger}-f(X_{t}^{\dagger},\Theta_{t}^{(i)}){\rm d}t-\gamma^{1/2}{\rm
d}W_{t}^{(i)}$ (8)

of the innovation in (5).

######  Remark 2.1.

There is an intriguing connection to the stochastic gradient descent approach
to the estimation of $\theta^{\dagger}$, as proposed in [32], which is written
as

${\rm
d}\theta_{t}&=\frac{\alpha_{t}}{\gamma}\nabla_{\theta}f(X_{t}^{\dagger},\theta_{t}){\rm
d}\tilde{I}_{t},\\\ {\rm d}\tilde{I}_{t}={\rm
d}X_{t}^{\dagger}-f(X_{t}^{\dagger},\theta_{t}){\rm d}t$ (9)

in our notation, where $\alpha_{t}>0$ denotes the learning rate. We note that
(9) shares with (3) the gain times innovation structure. However, while (3)
approximates the Bayesian inference problem, formulation (9) treats the
parameter estimation problem from an optimisation perspective. Both
formulations share, however, the discontinuous dependence on the observation
path $X_{t}^{\dagger}$, and the proposed frequentist analysis of the EnKBF (3)
also applies in simplified form to (9). We also point out that (3) is affine
invariant [19] and does not require the computation of partial derivatives.

We now state a numerical implementation with step-size $\Delta t>0$ and denote
the resulting numerical approximations at $t_{n}=n\Delta t$ by
$\Theta_{n}\sim\pi_{n}$, $n\geq 1$. While a standard Euler–Maruyama
approximation could be applied, the following stable discrete-time mean-field
formulation of the EnKBF

$\Theta_{n+1}=\Theta_{n}+K_{n}\left\\{(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})-\frac{1}{2}\left(f(X_{t_{n}}^{\dagger},\Theta_{n})+\pi_{n}[f(X^{\dagger}_{t_{n}},\theta)]\right)\Delta
t\right\\}$ (10)

is inspired by [4] with Kalman gain

$K_{n}&=\pi_{n}\left[(\theta-\pi_{n}[\theta])\otimes
f(X_{t_{n}}^{\dagger},\theta)\right]\times\\\ \quad\left(\gamma+\Delta
t\pi_{n}\left[\left(f(X_{t_{n}}^{\dagger},\theta)-\pi_{n}[f(X_{t_{n}}^{\dagger},\theta)]\right)\otimes
f(X_{t_{n}}^{\dagger},\theta)\right]\right)^{-1}.$ (11)

It is straightforward to combine this time discretisation with the Monte-Carlo
approximation (5) in order to obtain a complete numerical implementation of
the EnKBF.

######  Remark 2.2.

The rough path analysis of the EnKBF presented in [10] is based on a
Stratonovich reformulation of (3) and its appropriate time discretisation.
Here we follow the Itô/Euler–Maruyama formulation of the data-driven term in
(3),

$\int_{0}^{T}g(X_{t}^{\dagger},t)\,{\rm d}X_{t}^{\dagger}=\lim_{\Delta t\to
0}\sum_{i=1}^{L}g(X_{t_{n}}^{\dagger},t_{n})(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})$
(12)

for any continuous function $g(x,t)$ and $\Delta t=T/L$, as it corresponds to
standard implementation of the EnKBF and is easier to analyse in the context
of this paper.

The EnKBF provides only an approximate solution to the Bayesian inference
problem for general nonlinear $f(x,\theta)$. However, it becomes exact in the
mean-field limit for affine drift functions $f(x,\theta)=\theta Ax+Bx+c$.

######  Example 2.3.

Consider the stochastic partial differential equation

$\partial_{t}u=-U\partial_{y}u+\rho\partial_{y}^{2}u+\dot{\mathcal{W}}$ (13)

over a periodic spatial domain $y\in[0,L)$, where $\mathcal{W}(t,y)$ denotes
space-time white noise, $U\in\mathbb{R}$, and $\rho>0$ are given parameters. A
standard finite-difference discretisation in space with $d$ grid points and
mesh-size $\Delta y$ leads to a linear system of SDEs of the form

${\rm d}{\bf u}_{t}=-(UD+\rho DD^{\rm T}){\bf u}_{t}{\rm d}t+\Delta
y^{-1/2}{\rm d}W_{t},$ (14)

where ${\bf u}_{t}\in\mathbb{R}^{d}$ denotes the vector of grid approximations
at time $t$, $D\in\mathbb{R}^{d\times d}$ a finite difference approximation of
the spatial derivative $\partial_{y}$, and $W_{t}$ the standard
$d$-dimensional Brownian motion. We can now set $X_{t}={\bf u}_{t}$,
$\gamma=\Delta y^{-1}$ and identify either $\theta=U$ or $\theta=\rho$ as the
unknown parameter in order to obtain an SDE of the form (1).

In this note, we further simplify our given inference problem to the case

$f(x,\theta)=\theta Ax\,,$ (15)

where $A\in\mathbb{R}^{d\times d}$ is a normal matrix with eigenvalues in the
left half plane. That is $\sigma(A)\subset\mathbb{C}_{-}$. The reference
parameter value is set to $\theta^{\dagger}=1$. Hence the SDE (2) possesses a
Gaussian invariant measure with mean zero and covariance matrix

$C=-\gamma(A+A^{\rm T})^{-1}.$ (16)

We assume from now on that the observations $X_{t}^{\dagger}$ are realisations
of (2) with initial condition $X_{0}^{\dagger}\sim{\rm N}(0,C)$.

Under these assumptions, the EnKBF (3) simplifies drastically, and we obtain

${\rm d}\Theta_{t}&=\frac{\sigma_{t}}{\gamma}(AX^{\dagger}_{t})^{\rm T}{\rm
d}I_{t},\\\ {\rm d}I_{t}={\rm
d}X_{t}^{\dagger}-\frac{1}{2}\left(\Theta_{t}+\pi_{t}[\theta]\right)AX_{t}^{\dagger}{\rm
d}t,$ (17)

with variance

$\sigma_{t}=\pi_{t}\left[(\theta-\pi_{t}[\theta])^{2}\right].$ (18)

######  Remark 2.4.

For completeness, we state the corresponding formulation for the stochastic
gradient descent approach (9):

${\rm d}\theta_{t}&=\frac{\alpha_{t}}{\gamma}(AX_{t}^{\dagger})^{\rm T}{\rm
d}\tilde{I}_{t},\\\ {\rm d}\tilde{I}_{t}={\rm
d}X_{t}^{\dagger}-\theta_{t}AX_{t}^{\dagger}{\rm d}t.$ (19)

We find that the learning rate $\alpha_{t}$ takes the role of the variance
$\sigma_{t}$ in (17). However, we emphasise again that the same pathwise
stochastic integrals arise from both formulations, and therefore, the same
robustness issue of the resulting estimators $\theta_{t}$, $t>0$, arises.

Similarly, the discrete-time mean-field EnKBF (10) reduces to

$\Theta_{n+1}=\Theta_{n}+K_{n}\left\\{(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})-\frac{1}{2}\left(\Theta_{n}+\pi_{n}[\theta]\right)AX_{t_{n}}^{\dagger}\Delta
t\right\\}$ (20)

with Kalman gain

$K_{n}=\sigma_{n}(AX_{t_{n}}^{\dagger})^{\rm T}\left(\gamma+\Delta
t\sigma_{n}(AX_{t_{n}}^{\dagger})^{\rm T}AX_{t_{n}}^{\dagger}\right)^{-1}\,.$
(21)

Furthermore, since $X_{t}^{\dagger}\sim{\rm N}(0,C)$,

$(AX_{t}^{\dagger})^{\rm T}AX_{t}^{\dagger}=(A^{\rm
T}A):(X_{t}^{\dagger}\otimes X_{t}^{\dagger})\approx(A^{\rm T}A):C$ (22)

for $d\gg 1$, and we may simplify the Kalman gain to

$K_{n}=\sigma_{n}\,(AX_{t_{n}}^{\dagger})^{\rm T}\left(\gamma+\Delta
t\sigma_{n}\,(A^{\rm T}A):C\right)^{-1}.$ (23)

Here we have used the notation $A:B=\mbox{tr}(A^{\rm T}B)$ to denote the
Frobenius inner product of two matrices $A,B\in\mathbb{R}^{d\times d}$. The
approximation (22) becomes exact in the limit $d\to\infty$, which we will
frequently assume in the following section. Please note that

$K_{n}=\frac{\sigma_{n}}{\gamma}\,(AX_{t_{n}}^{\dagger})^{\rm
T}+\mathcal{O}(\Delta t)$ (24)

under the stated assumptions.

######  Remark 2.5.

The Stratonovitch reformulation of (17) replaces (17a) by

${\rm d}\Theta_{t}=\frac{\sigma_{t}}{\gamma}\left\\{(AX^{\dagger}_{t})^{\rm
T}\circ{\rm d}I_{t}-\frac{\gamma}{2}\mbox{tr}\,(A)\,{\rm d}t\right\\}.$ (25)

The innovation $I_{t}$ remains as before. See Appendix B of [10] for more
details. An appropriate time discretisation of the innovation-driven term
replaces the Kalman gain (21) by

$K_{n+1/2}=\sigma_{n}(AX_{t_{n+1/2}}^{\dagger})^{\rm T}\left(\gamma+\Delta
t\sigma_{n}(AX_{t_{n+1/2}}^{\dagger})^{\rm
T}AX_{t_{n+1/2}}^{\dagger}\right)^{-1},$ (26)

where

$X_{t_{n+1/2}}^{\dagger}=\frac{1}{2}(X_{t_{n}}^{\dagger}+X_{t_{n+1}}^{\dagger})\,.$
(27)

Please note that a midpoint discretisation of the data-driven term in (25)
results in

$(AX_{t_{n+1/2}}^{\dagger})^{\rm
T}(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})&=(AX_{t_{n}}^{\dagger})^{\rm
T}(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})\,\,+\\\ \,\,\,\frac{1}{2}A^{\rm
T}:(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})\otimes(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})$
(28)

and that

$\frac{1}{2}A^{\rm
T}:(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})\otimes(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})\approx\frac{\Delta
t\,\gamma}{2}\mbox{tr}\,(A),$ (29)

which justifies the additional drift term in (25). A precise meaning of the
approximation in (29) will be given in Remark 3.2 below.

Alternatively, if one wishes to explicitly utilise the availability of
continuous-time data $X^{\dagger}_{t}$, one could apply the following variant
of (20):

$\Theta_{n+1}=\Theta_{n}+\frac{\sigma_{n}}{\gamma}\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm
T}{\rm
d}X_{t}^{\dagger}-\frac{1}{2}K_{n}AX_{t_{n}}^{\dagger}\left(\Theta_{n}+\pi_{n}[\theta]\right)\Delta
t,$ (30)

and following the Itô/Euler–Maruyama approximation (12), discretise the
integral with a small inner step-size $\Delta\tau=\Delta t/L$, $L\gg 1$; that
is,

$\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm T}{\rm
d}X_{t}^{\dagger}\approx\sum_{l=0}^{L-1}(AX_{\tau_{l}}^{\dagger})^{\rm
T}(X_{\tau_{l+1}}^{\dagger}-X_{\tau_{l}}^{\dagger})$ (31)

with $\tau_{l}=t_{n}+l\Delta\tau$. We note that

$\sum_{l=0}^{L-1}(AX_{\tau_{l}}^{\dagger})^{\rm
T}(X_{\tau_{l+1}}^{\dagger}-X_{\tau_{l}}^{\dagger})&=(AX_{t_{n}}^{\dagger})^{\rm
T}(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})\,\,+\\\ \,\,A^{\rm
T}:\left(\sum_{l=0}^{L-1}(X_{\tau_{l}}^{\dagger}-X_{t_{n}}^{\dagger})\otimes(X_{\tau_{l+1}}^{\dagger}-X_{\tau_{l}}^{\dagger})\right),$
(32)

which is at the heart of rough path analysis [14] and which we utilise in the
following section.

##  3 Frequentist analysis

It is well-known that the second-order contribution in (32) leads to a
discontinuous dependence of the integral on the observed $X_{t}^{\dagger}$ in
the uniform norm topology on the space of continuous functions. Rough path
theory fixes this problem by defining appropriately extended topologies and
has been extended to the EnKBF in [10]. In this section, we complement the
path-wise analysis from [10] by an analysis of the impact of second-order
contribution on the EnKBF (17) from a frequentist perspective, which analyses
the behaviour of EnKBF over all possible observations $X_{t}^{\dagger}$
subject to (2). In other words, one switches from a strong solution concept to
a weak one. While we assume that the observations satisfy (2), throughout this
section, we will analyse the impact of a perturbed observation process on the
EnKBF in Section 4.

We first derive evolution equations for the conditional mean and variance
under the assumption that $\Theta_{0}$ is Gaussian distributed with given
prior mean $m_{\rm prior}$ and variance $\sigma_{\rm prior}$. It follows
directly from (17) that the conditional mean $\mu_{t}=\pi_{t}[\theta]$, that
is the mean of $\Theta_{t}$, satisfies the SDE

${\rm d}\mu_{t}=\frac{\sigma_{t}}{\gamma}\left((AX_{t}^{\dagger})^{\rm T}{\rm
d}X^{\dagger}_{t}-\mu_{t}\,(A^{\rm T}A):(X_{t}^{\dagger}\otimes
X_{t}^{\dagger})\,{\rm d}t\right),$ (33)

which simplifies to

${\rm d}\mu_{t}=\frac{\sigma_{t}}{\gamma}\left((AX_{t}^{\dagger})^{\rm T}{\rm
d}X^{\dagger}_{t}-\mu_{t}\,(A^{\rm T}A):C\,{\rm d}t\right),$ (34)

under the approximation (22). The initial condition is $\mu_{0}=m_{\rm
prior}$. The evolution equation for the conditional variance, that is the
variance of $\Theta_{t}$, is given by

$\frac{\rm d}{{\rm d}t}\sigma_{t}=-\frac{\sigma_{t}^{2}}{\gamma}\,(A^{\rm
T}A):(X_{t}^{\dagger}\otimes X_{t}^{\dagger})$ (35)

with initial condition $\sigma_{0}=\sigma_{\rm prior}$ and which again reduces
to

$\frac{\rm d}{{\rm d}t}\sigma_{t}=-\frac{\sigma_{t}^{2}}{\gamma}\,(A^{\rm
T}A):C$ (36)

under the approximation (22).

We now perform a frequentist analysis of the estimator $\mu_{t}$ defined by
(34) and (36), that is, we perform a weak analysis of the SDE (34) in terms of
the first two moments of $\mu_{t}$ [31]. In the first step, we take the
expectation of (34) over all realisations $X_{t}^{\dagger}$ of the SDE (2),
which we denote by

$m_{t}:=\mathbb{E}^{\dagger}[\mu_{t}].$ (37)

The associated evolution equation is given by

$\frac{\rm d}{{\rm d}t}m_{t}=\frac{\sigma_{t}}{\gamma}\,(A^{\rm
T}A):\mathbb{E}^{\dagger}\left[X_{t}^{\dagger}\otimes
X_{t}^{\dagger}\right]-\frac{\sigma_{t}}{\gamma}\,(A^{\rm T}A):C\,m_{t},$ (38)

which reduces to

$\frac{\rm d}{{\rm d}t}m_{t}=\frac{\sigma_{t}}{\gamma}\,(A^{\rm
T}A):C\,(1-m_{t})=\sigma_{t}\,(A^{\rm T}A):(A+A^{\rm T})^{-1}\,(1-m_{t}).$
(39)

In the second step, we also look at the frequentist variance

$p_{t}:=\mathbb{E}^{\dagger}[(\mu_{t}-m_{t})^{2}].$ (40)

Using

${\rm d}(\mu_{t}-m_{t})&=\frac{\sigma_{t}}{\gamma}\left\\{(A^{\rm
T}A):\left(X_{t}^{\dagger}\otimes X_{t}^{\dagger}-C\right){\rm
d}t+\gamma^{1/2}(AX_{t}^{\dagger})^{\rm T}{\rm
d}W^{\dagger}_{t}\right\\}\,\,-\\\
\qquad\qquad\frac{\sigma_{t}}{\gamma}(A^{\rm T}A):C\,(\mu_{t}-m_{t}){\rm d}t,$
(41)

we obtain

$\frac{\rm d}{{\rm d}t}p_{t}&=-\frac{\sigma_{t}}{\gamma}\,(A^{\rm
T}A):C\left(2p_{t}-\sigma_{t}\right)\,\,+\\\
\qquad\qquad\frac{2\sigma_{t}}{\gamma}\,(A^{\rm
T}A):\mathbb{E}^{\dagger}\left[(X_{t}^{\dagger}\otimes
X_{t}^{\dagger}-C)\,(\mu_{t}-m_{t})\right],$ (42)

which we simplify to

$\frac{\rm d}{{\rm d}t}p_{t}=\frac{\sigma_{t}}{\gamma}\,(A^{\rm
T}A):C\left(\sigma_{t}-2p_{t}\right)=\sigma_{t}\,(A^{\rm T}A):(A+A^{\rm
T})^{-1}\left(\sigma_{t}-2p_{t}\right)$ (43)

under the approximation (22). The initial conditions are $m_{0}=m_{\rm prior}$
and $p_{0}=0$, respectively. We note that the differential equations (36) and
(43) are explicitly solvable. For example, it holds that

$\sigma_{t}=\frac{\sigma_{0}}{1+(A^{\rm T}A):(A^{\rm T}+A)^{-1}\,\sigma_{0}t}$
(44)

and one finds that $\sigma_{t}\sim 1/((A^{\rm T}A):(A^{\rm T}+A)^{-1}\,t)$ for
$t\gg 1$. It can also be shown that $p_{t}\leq\sigma_{t}$ for all $t\geq 0$.
Furthermore, this analysis suggests that the learning rate in the stochastic
gradient descent formulation (19) should be chosen as

$\alpha_{t}=\min\left\\{\bar{\alpha},\frac{1}{(A^{\rm T}A):(A^{\rm
T}+A)^{-1}\,t}\right\\},$ (45)

where $\bar{\alpha}>0$ denotes an initial learning rate; for example
$\bar{\alpha}=\sigma_{0}$.

We finally conduct a formal analysis of the ensemble Kalman filter time-
stepping (20) and demonstrate that the method is first-order accurate with
regard to the implied frequentist mean $m_{t}$. We recall (24) and conclude
from (20) that the implied update on the variance $\sigma_{n}$ satisfies

$\sigma_{n+1}=\sigma_{n}-\frac{\sigma_{n}^{2}}{\gamma}\,(A^{\rm T}A):C\Delta
t+\mathcal{O}(\Delta t^{2}),$ (46)

which provides a first-order approximation to (36).

We next analyse the evolution equation (34) for the conditional mean $\mu_{t}$
and its numerical approximation

$\mu_{n+1}=\mu_{n}+K_{n}\left\\{(X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger})-\mu_{n}AX_{t_{n}}^{\dagger}\Delta
t\right\\}$ (47)

arising from (20). Here we follow [14] in order to analyse the impact of the
data $X_{t}^{\dagger}$ on the estimator. An in-depth theoretical treatment can
be found in [10].

Comparing (47) to (34) and utilising (24), we find that the key quantity of
interest is

$J^{\dagger}_{t_{n},t_{n+1}}:=\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm
T}{\rm d}X_{t}^{\dagger},$ (48)

which we can rewrite as

$J^{\dagger}_{t_{n},t_{n+1}}=A^{\rm T}:(X^{\dagger}_{t_{n}}\otimes
X^{\dagger}_{t_{n},t_{n+1}})+A^{\rm
T}:\mathbb{X}_{t_{n},t_{n+1}}^{\dagger}\,.$ (49)

Here, motivated by (32) and following standard rough path notation, we have
used

$X^{\dagger}_{t_{n},t_{n+1}}:=X_{t_{n+1}}^{\dagger}-X_{t_{n}}^{\dagger}$ (50)

and the second-order iterated Itô integral

$\mathbb{X}_{t_{n},t_{n+1}}^{\dagger}:=\int_{t_{n}}^{t_{n+1}}(X^{\dagger}_{t}-X^{\dagger}_{t_{n}})\otimes{\rm
d}X_{t}^{\dagger}.$ (51)

The difference between the integral (48) and its corresponding approximation
in (47) is provided by $A^{\rm T}:\mathbb{X}_{t_{n},t_{n+1}}^{\dagger}$ plus
higher-order terms arising from (24). The iterated integral
$\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}$ becomes a random variable from the
frequentist perspective. Taking note of (2), we find that the drift,
$f(x)=Ax$, contributes with terms of order $\mathcal{O}(\Delta t^{2})$ to
$\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}$ and the expected value of
$\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}$ therefore satisfies

$\mathbb{E}^{\dagger}[\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}]=\mathcal{O}(\Delta
t^{2}),$ (52)

since $\mathbb{E}^{\dagger}[W^{\dagger}_{t_{n},\tau}]=0$ for $\tau>t_{n}$, and

$\mathbb{E}^{\dagger}[\mathbb{W}_{t_{n},t_{n+1}}^{\dagger}]=\frac{1}{2}\mathbb{E}^{\dagger}[W^{\dagger}_{t_{n},t_{n+1}}\otimes
W^{\dagger}_{t_{n},t_{n+1}}-[W^{\dagger}_{t_{n}},W^{\dagger}_{t_{n},t_{n+1}}]]-\frac{\Delta
t}{2}I=0,$ (53)

where we have introduced the commutator

$[W_{t_{n}}^{\dagger},W_{t_{n},t_{n+1}}^{\dagger}]:=W^{\dagger}_{t_{n}}\otimes
W^{\dagger}_{t_{n},t_{n+1}}-W^{\dagger}_{t_{n},t_{n+1}}\otimes
W_{t_{n}}^{\dagger}.$ (54)

Hence we find that, while (47) is not a first-order (strong) approximation of
the SDE (34), the approximation becomes first-order in $m_{t}$ when averaged
over realisations $X_{t}^{\dagger}$ of the SDE (2). More precisely, one
obtains

$\mathbb{E}^{\dagger}[J^{\dagger}_{t_{n},t_{n+1}}]=(A^{\rm T}A):C\Delta
t+\mathcal{O}(\Delta t^{2}).$ (55)

We note that the modified scheme (30) leads to the same time evolution in the
variance $\sigma_{n}$ while the update in $\mu_{n}$ is changed to

$\mu_{n+1}=\mu_{n}+\frac{\sigma_{n}}{\gamma}\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm
T}{\rm d}X_{t}^{\dagger}-K_{n}AX_{t_{n}}^{\dagger}\mu_{n}\Delta t.$ (56)

This modification results in a more accurate evolution in the conditional mean
$\mu_{n}$, but because of (52) it does not impact to leading order the
evolution of the underlying frequentist mean,
$m_{n}=\mathbb{E}^{\dagger}[\mu_{n}]$. We summarise our findings in the
following proposition.

######  Proposition 3.1.

The discrete-time EnKBF implementations (20) and (30) both provide first-order
approximations to the time evolution of the frequentist mean, $m_{t}$, and the
frequentist variance, $p_{t}$. In other words, both methods converge weakly
with order one.

We also note that the frequentist uncertainty is essentially data-independent
and depends only on the time window $[0,T]$ over which the data gets observed.
Hence, for fixed observation interval $[0,T]$, it makes sense to choose the
step-size $\Delta t$ such that the discretisation error (bias) remains on the
same order of magnitude as $p_{T}^{1/2}\approx\sigma_{T}^{1/2}$. Selecting a
much smaller step-size would not significantly reduce the frequentist
estimation error in the conditional estimator $\mu_{T}$.

######  Remark 3.2.

We can now give a precise reformulation of the approximation (29):

$\frac{1}{2}\mathbb{E}^{\dagger}\left[A^{\rm
T}:(X_{t_{n},t_{n+1}}^{\dagger}\otimes
X_{t_{n},t_{n+1}}^{\dagger})\right]=\frac{\Delta
t\,\gamma}{2}\mbox{tr}\,(A)+\mathcal{O}(\Delta t^{2}),$ (57)

which is at the heart of the Stratonovich formulation (25) of the EnKFB [10].

##  4 Multi-scale data

We now have all the material in place to study the dependency of the EnKBF
estimator on a set of observations $X_{t}^{(\epsilon)}$, $\epsilon>0$, which
approach the theoretical $X_{t}^{\dagger}$ with respect to the uniform norm
topology on the space of continuous functions as $\epsilon\to 0$. Since the
second-order contribution in (32), that is (51), does not depend continuously
on such perturbations, we demonstrate in this section that a systematic bias
arises in the EnKBF. Furthermore, we show how the bias can be eliminated
either via subsampling the data, which effectively amounts to ignoring these
second-order contributions, or via an appropriate correction term, which
ensures a continuous dependence on observations $X_{t}^{(\epsilon)}$ with
respect to the uniform norm topology. More specifically, we investigate the
impact of a possible discrepancy between the SDE model (1), for which we aim
to estimate the parameter $\theta$, and the data generating SDE (2). We
therefore replace (2) by the the following two-scale SDE [18]:

${\rm d}X^{(\epsilon)}_{t}&=AX^{(\epsilon)}_{t}\,{\rm
d}t+\frac{\gamma^{1/2}}{\epsilon}MP^{(\epsilon)}_{t}\,{\rm d}t,\\\ {\rm
d}P^{(\epsilon)}_{t}=-\frac{1}{\epsilon}MP^{(\epsilon)}_{t}\,{\rm d}t+{\rm
d}W^{\dagger}_{t},$ (58)

where

$M=\left(\begin{array}[]{cc}1&\beta\\\ -\beta&1\end{array}\right),$ (59)

$\beta=2$ and $\epsilon=0.01$. The dimension of state space is $d=2$
throughout this section. While we restrict here to the simple two-scale model
(58), similar scenarios can arise from deterministic fast-slow systems [25,
8].

The associated EnKBF mean-field equations in the parameter $\Theta_{t}$, which
we now denote by $\Theta_{t}^{(\epsilon)}$ in order to explicitly record its
dependence on the scale parameter $\epsilon\ll 1$, become

${\rm
d}\Theta_{t}^{(\epsilon)}&=\frac{\sigma_{t}^{(\epsilon)}}{\gamma}(AX^{(\epsilon)}_{t})^{\rm
T}{\rm d}I_{t}^{(\epsilon)},\\\ {\rm d}I_{t}^{(\epsilon)}={\rm
d}X_{t}^{(\epsilon)}-\frac{1}{2}\left(\Theta_{t}^{(\epsilon)}+\pi_{t}^{(\epsilon)}[\theta]\right)AX_{t}^{(\epsilon)}{\rm
d}t,$ (60)

with variance

$\sigma_{t}^{(\epsilon)}=\pi_{t}^{(\epsilon)}\left[(\theta-\pi_{t}^{(\epsilon)}[\theta])^{2}\right]$
(61)

and $\Theta^{\epsilon}_{t}\sim\pi_{t}^{(\epsilon)}$. The discrete-time mean-
field EnKBF (20) turns into

$\Theta_{n+1}^{(\epsilon)}=\Theta_{n}^{(\epsilon)}+K_{n}^{(\epsilon)}\left\\{\left(X_{t_{n+1}}^{(\epsilon)}-X_{t_{n}}^{(\epsilon)}\right)-\frac{1}{2}\left(\Theta_{n}^{(\epsilon)}+\pi_{n}^{(\epsilon)}[\theta]\right)AX_{t_{n}}^{(\epsilon)}\Delta
t\right\\}$ (62)

with Kalman gain

$K_{n}^{(\epsilon)}=\sigma_{n}^{(\epsilon)}(AX_{t_{n}}^{(\epsilon)})^{\rm
T}\left(\gamma+\Delta t\sigma_{n}^{(\epsilon)}(AX_{t_{n}}^{(\epsilon)})^{\rm
T}AX_{t_{n}}^{(\epsilon)}\right)^{-1}\,.$ (63)

We also consider the appropriately modified scheme (30):

$\Theta_{n+1}^{(\epsilon)}=\Theta_{n}^{(\epsilon)}+\frac{\sigma_{n}^{(\epsilon)}}{\gamma}\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm
T}{\rm
d}X_{t}^{(\epsilon)}-\frac{1}{2}K_{n}^{(\epsilon)}AX_{t_{n}}^{(\epsilon)}\left(\Theta_{n}^{(\epsilon)}+\pi_{n}^{(\epsilon)}[\theta]\right)\Delta
t.$ (64)

In order to understand the impact of the modified data generating process on
the two mean-field EnKBF formulations (62) and (64), respectively, we follow
[18] and investigate the difference between $X^{(\epsilon)}_{t}$ and
$X^{\dagger}_{t}$:

${\rm
d}(X^{(\epsilon)}_{t}-X^{\dagger}_{t})&=A(X^{(\epsilon)}_{t}-X^{\dagger}_{t}){\rm
d}t+\frac{\gamma^{1/2}}{\epsilon}MP_{t}^{(\epsilon)}{\rm d}t-\gamma^{1/2}{\rm
d}W_{t}^{\dagger}\\\ =A(X_{t}^{(\epsilon)}-X_{t}^{\dagger}){\rm
d}t-\gamma^{1/2}{\rm d}P_{t}^{(\epsilon)}.$ (65)

When $P^{(\epsilon)}_{t}$ is stationary, it is Gaussian with mean zero and
covariance

$\mathbb{E}_{\rm stat}\left[P_{t}^{(\epsilon)}\otimes
P_{t}^{(\epsilon)}\right]=\epsilon\,(M+M^{\rm T})^{-1}=\frac{\epsilon}{2}I.$
(66)

Hence $P^{(\epsilon)}_{t}\rightarrow 0$ as $\epsilon\rightarrow 0$ and also

$X^{(\epsilon)}_{t}\rightarrow X^{\dagger}_{t}$ (67)

in $L^{2}$ uniformly in $t$, provided $\sigma(A)\subset\mathbb{C}_{-}$ and
$X^{(\epsilon)}_{0}=X^{\dagger}_{0}$. This is illustrated in Figure 1.

Figure 1: SDE driven by mathematical vs. physical Brownian motion
($\epsilon=0.01$). The top panel displays both $X_{t}^{\dagger}$ (blue) and
$X_{t}^{(\epsilon)}$ (red) over the long time interval $t\in[0,10]$, while the
lower panel provides a zoomed in perspective over the interval $t\in[0,1]$.

In order to investigate the problem further, we study the integral

$J^{(\epsilon)}_{t_{n},t_{n+1}}:=\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm
T}{\rm d}X_{t}^{(\epsilon)}$ (68)

and its relation to (48). As for (48), we can rewrite (68) as

$J^{(\epsilon)}_{t_{n},t_{n+1}}=A^{\rm T}:(X^{(\epsilon)}_{t_{n}}\otimes
X^{(\epsilon)}_{t_{n},t_{n+1}})+A^{\rm
T}:\mathbb{X}_{t_{n},t_{n+1}}^{(\epsilon)}.$ (69)

We now investigate the limit of the second-order iterated integral

$\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}&=\int_{t_{n}}^{t_{n+1}}X^{(\epsilon)}_{t_{n},t}\otimes{\rm
d}X_{t}^{(\epsilon)}\\\ =\frac{1}{2}X_{t_{n},t_{n+1}}^{(\epsilon)}\otimes
X_{t_{n},t_{n+1}}^{(\epsilon)}-\frac{1}{2}\int_{t_{n}}^{t_{n+1}}[X_{t_{n},t}^{(\epsilon)},{\rm
d}X_{t}^{(\epsilon)}]$ (70)

as $\epsilon\to 0$ [18]. Here $[.,.]$ denotes the commutator defined by (54).

######  Proposition 4.1.

The second-order iterated integral $\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}$
satisfies

$\lim_{\epsilon\to
0}\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}=\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}+\frac{\Delta
t\,\gamma}{2}M$ (71)

######  Proof 4.2.

The proof follows [18] and can be summarised as follows:

$\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}&=\int_{t_{n}}^{t_{n+1}}X^{(\epsilon)}_{t_{n},t}\otimes{\rm
d}X^{(\epsilon)}_{t}\\\
\rightarrow\int_{t_{n}}^{t_{n+1}}X^{\dagger}_{t_{n},t}\otimes{\rm
d}X^{\dagger}_{t}-\gamma^{1/2}\int_{t_{n}}^{t_{n+1}}X^{(\epsilon)}_{t_{n},t}\otimes{\rm
d}P^{(\epsilon)}_{t}\\\
=\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}-\gamma^{1/2}X^{(\epsilon)}_{t_{n},t_{n+1}}\otimes
P^{(\epsilon)}_{t_{n+1}}+\gamma^{1/2}\int_{t_{n}}^{t_{n+1}}{\rm
d}X^{(\epsilon)}_{t}\otimes P^{(\epsilon)}_{t}\\\
\rightarrow\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}+\gamma^{1/2}\int_{t_{n}}^{t_{n+1}}\left\\{AX_{t}^{(\epsilon)}+\frac{\gamma^{1/2}}{\epsilon}MP^{(\epsilon)}_{t}\right\\}\otimes
P^{(\epsilon)}_{t}{\rm d}t\\\
\rightarrow\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}+\frac{\Delta
t\,\gamma}{\epsilon}M\,\mathbb{E}_{\rm
stat}\left[P_{t_{n}}^{(\epsilon)}\otimes P_{t_{n}}^{(\epsilon)}\right]\\\
=\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}+\frac{\Delta t\,\gamma}{2}M.$ (72)

As discussed in detail in [10] already, Proposition 4.1 implies that the
scheme (64) does not, in general, converge to the scheme (64) as $\epsilon\to
0$ since

$J^{\dagger}_{t_{n},t_{n+1}}=\lim_{\epsilon\to
0}J^{(\epsilon)}_{t_{n},t_{n+1}}-\frac{\Delta t\,\gamma}{2}A^{\rm T}:M\,.$
(73)

This observation suggests the following modification

$\Theta_{n+1}^{(\epsilon)}&=\Theta_{n}^{(\epsilon)}+\frac{\sigma_{n}^{(\epsilon)}}{\gamma}\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm
T}{\rm d}X_{t}^{(\epsilon)}-\frac{\Delta t}{2}\sigma_{n}^{(\epsilon)}\,A^{\rm
T}:M\,\,-\\\
\qquad\qquad\frac{1}{2}K_{n}^{(\epsilon)}AX_{t_{n}}^{(\epsilon)}\left(\Theta_{n}^{(\epsilon)}+\pi_{n}^{(\epsilon)}[\theta]\right)\Delta
t$ (74)

to (64). Please note that it follows from (70) that

$\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm T}{\rm
d}X_{t}^{(\epsilon)}=A^{\rm T}:\left(X_{t_{n+1/2}}^{(\epsilon)}\otimes
X_{t_{n},t_{n+1}}^{(\epsilon)}-\frac{1}{2}\int_{t_{n}}^{t_{n+1}}[X_{t_{n},t}^{(\epsilon)},{\rm
d}X_{t}^{(\epsilon)}]\right).$ (75)

######  Proposition 4.3.

The discrete-time EnKBF (62) converges to (20) for fixed $\Delta t$ as
$\epsilon\to 0$. Similarly, (74) converges to (30) under the same limit.

######  Proof 4.4.

The first statement follows from $\sigma_{n}^{(\epsilon)}=\sigma_{n}$, the
limiting behaviour (67), and

$\lim_{\epsilon\to 0}K_{n}^{(\epsilon)}=K_{n}.$ (76)

The second statement additionally requires (73) to be substituted into (74)
when taking the limit $\epsilon\to 0$.

######  Remark 4.5.

The analogous adaptation of (74) to the gradient descent formulation (19) with
$X_{t}^{\dagger}$ replaced by $X_{t}^{(\epsilon)}$ becomes

$\theta_{n+1}^{(\epsilon)}&=\theta_{n}^{(\epsilon)}+\frac{\alpha_{t_{n}}}{\gamma}\left(\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm
T}{\rm d}X_{t}^{(\epsilon)}-\frac{\gamma\Delta t}{2}A^{\rm T}:M\,\,-\right.\\\
\qquad\qquad\left.\theta_{n}^{(\epsilon)}(AX_{t_{n}}^{(\epsilon)})^{\rm
T}AX_{t_{n}}^{(\epsilon)}\Delta t\right).$ (77)

Alternatively, subsampling the data can be applied which leads to the simpler
formulation

$\theta_{n+1}^{(\epsilon)}=\theta_{n}^{(\epsilon)}+\frac{\alpha_{t_{n}}}{\gamma}(AX_{t_{n}}^{(\epsilon)})^{\rm
T}\left((X_{t_{n+1}}^{(\epsilon)}-X_{t_{n}}^{(\epsilon)})-\theta_{n}^{(\epsilon)}AX_{t_{n}}^{(\epsilon)}\Delta
t\right).$ (78)

######  Remark 4.6.

A two-scale SDE, closely related to (58), has been investigated in [9] in
terms of the time integrated autocorrelation function of $P_{t}^{(\epsilon)}$
and modified stochastic integrals. In our case, the modified quadrature rule,
here denoted by $\diamond$, has to satisfy

$\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm T}\diamond{\rm
d}X_{t}^{\dagger}=\lim_{\epsilon\to
0}\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm T}{\rm
d}X_{t}^{(\epsilon)},$ (79)

and it is therefore related to the standard Itô integral via

$\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm T}\diamond{\rm
d}X_{t}^{\dagger}=\int_{t_{n}}^{t_{n+1}}(AX_{t}^{\dagger})^{\rm T}{\rm
d}X_{t}^{\dagger}+\frac{\Delta t\gamma}{2}A^{\rm T}:M.$ (80)

Hence $M$ playes the role of the integrated autocorrelation function of
$P_{t}^{(\epsilon)}$ in our approach. We note that the modified quadrature
rule reduces to the standard Stratonovitch integral if either $\beta=0$ in
(59) or $A$ is symmetric. While the results from [9] could, therefore, also be
used as a starting point for discussing the induced estimation bias, practical
implementations would still require knowledge of the integrated
autocorrelation function of $P_{t}^{(\epsilon)}$ or, equivalently, the
estimation of $M$ in addition to observing $X_{t}^{(\epsilon)}$. We address
this aspect next.

###  4.1 Numerical implementation

The numerical implementation of (74) requires an estimator for the generally
unknown $M$ in (73). This task is challenging as we only have access to
$X_{t}^{(\epsilon)}$ without any explicit knowledge of the underlying
generating process (58). While the estimator proposed in [10] is based on the
idea of subsampling the data, the frequentist perspective taken in this note
suggests the alternative estimator $M_{\rm est}$ defined by

$\frac{\Delta t\,\gamma}{2}M_{\rm
est}=\mathbb{E}^{\dagger}[\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}],$ (81)

which follows from (72f) and (52). That is,
$\mathbb{E}^{\dagger}[\mathbb{X}^{\dagger}_{t_{n},t_{n+1}}]=\mathcal{O}(\Delta
t^{2})$ for $\Delta t$ sufficiently small. Note that second-order iterated
integral $X_{t_{n},t_{n+1}}^{(\epsilon)}$ satisfies (70) and is therefore easy
to compute. In practice, the frequentist expectation value can be replaced by
an approximation along a given single observation path $X^{(\epsilon)}_{t}$,
$t\in[0,T]$, under the assumption of ergodicity.

An appropriate choice of the outer or sub-sampling step-size $\Delta t$ [28]
constitutes an important aspect for the practical implementation of the EnKBF
formulation (62) for finite values of $\epsilon>0$ [27]. Consistency of the
second-order iterated integrals [14] implies

$\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+2}}=\mathbb{X}^{(\epsilon)}_{t_{n},t_{n+1}}+\mathbb{X}^{(\epsilon)}_{t_{n+1},t_{n+2}}+X^{(\epsilon)}_{t_{n},t_{n+1}}\otimes
X^{(\epsilon)}_{t_{n+1},t_{n+2}}.$ (82)

A sensible choice of $\Delta t$ is dictated by

$\mathbb{E}^{\dagger}\left[X^{(\epsilon)}_{t_{n},t_{n+1}}\otimes
X^{(\epsilon)}_{t_{n+1},t_{n+2}}\right]=\mathcal{O}(\Delta t^{2})\,,$ (83)

that is, the sub-sampled data $X_{t_{n}}^{(\epsilon)}$ behaves to leading
order like solution increments from the reference model (2) at scale $\Delta
t$ independent of the specific value of $\epsilon$. Note that, on the other
hand,

$\mathbb{E}^{\dagger}\left[X^{(\epsilon)}_{\tau_{l},\tau_{l+1}}\otimes
X^{(\epsilon)}_{\tau_{l+1},\tau_{l+2}}\right]=\mathcal{O}(\epsilon^{-1}\Delta\tau^{2})$
(84)

for an inner step-size $\Delta\tau\sim\epsilon$. In other words, a suitable
step-size $\Delta t>0$ can be defined by making

$h(\Delta t):=\Delta
t^{-2}\left\|\mathbb{E}^{\dagger}\left[X^{(\epsilon)}_{t_{n},t_{n+1}}\otimes
X^{(\epsilon)}_{t_{n+1},t_{n+2}}\right]\right\|$ (85)

as small as possible while still guaranteeing an accurate numerical
approximation in (62).

######  Remark 4.7.

The choice of the outer time step $\Delta t$ is less critical for the EnKBF
formulation (74) since it does not rely on sub-sampling the data and is robust
with regard to perturbations in the data provided the appropriate $M$ is
explicitly available or has been estimated from the available data using (81).
Furthermore, if $A$ is symmetric, then it follows from (75) and the skew-
symmetry of the commutator $[.,.]$ that

$\int_{t_{n}}^{t_{n+1}}(AX_{t}^{(\epsilon)})^{\rm T}{\rm
d}X_{t}^{(\epsilon)}=A:\left(X^{(\epsilon)}_{t_{n+1/2}}\otimes
X^{(\epsilon)}_{t_{n},t_{n+1}}\right),$ (86)

which can be used in (74). The same simplification arises when $M$ is
symmetric. This insight is at the heart of the geometric rough path approach
followed in [10] and which starts from the Stratonovich formulation (25) of
the EnKBF. See also [29] on the convergence of Wong–Zakai approximations for
stochastic differential equations. In all other cases, a more refined
numerical approximation of the data-driven integral in (74) is necessary; such
as, for example, (31). For that reason, we rely on the Itô/Euler–Maruyama
interpretation of (68) in this note instead, that is the approximation (12).

###  4.2 Filtered data

We finally discuss a recently proposed [2] robust modification to the
parameter estimation problem in the light of the mean-field EnKBF equations
considered in this paper. The essential idea is to filter the observation
paths $X_{t}^{\dagger}$, $t\geq 0$, via

${\rm d}Z_{t}^{\dagger}=\frac{1}{\delta}(X_{t}^{\dagger}-Z_{t}^{\dagger}){\rm
d}t+\delta_{\rm noise}\sqrt{2}{\rm d}V_{t}^{\dagger},$ (87)

where $\delta>0$ is a sufficiently small parameter and $V_{t}^{\dagger}$
denotes independent Brownian motion with $\delta_{\rm noise}=1$ (noise added)
or $\delta_{\rm noise}=0$ (no noise added). Extending the methodology proposed
in [2] to the mean-field EnKFB equations (17), we now consider

${\rm d}\Theta_{t}&=\frac{\sigma_{t}}{\gamma}(AZ^{\dagger}_{t})^{\rm T}{\rm
d}I_{t},\\\ {\rm d}I_{t}={\rm
d}X_{t}^{\dagger}-\frac{1}{2}\left(\Theta_{t}+\pi_{t}[\theta]\right)AX_{t}^{\dagger}{\rm
d}t,$ (88)

with the variance $\sigma_{t}$ defined as before.

Let us first investigate the long-time behaviour of the extended data
generating system

${\rm d}X_{t}^{\dagger}&=AX_{t}^{\dagger}{\rm d}t+\gamma^{1/2}{\rm
d}W_{t}^{\dagger},\\\ {\rm
d}Z_{t}^{\dagger}=\frac{1}{\delta}(X_{t}^{\dagger}-Z_{t}^{\dagger}){\rm
d}t+\delta_{\rm noise}\sqrt{2}{\rm d}V_{t}^{\dagger},$ (89)

in some detail. Its stationary distribution is Gaussian with mean
$m_{\infty}^{x}=m_{\infty}^{z}=0$. The stationary covariance matrices satisfy
the relations

$0&=A\Sigma_{\infty}^{xx}+\Sigma_{\infty}^{xx}A^{\rm T}+\gamma I,\\\
0=\Sigma_{\infty}^{xz}+\Sigma_{\infty}^{zx}-2(\Sigma_{\infty}^{zz}-\delta_{\rm
noise}\delta I),\\\
0=A\Sigma_{\infty}^{xz}+\frac{1}{\delta}(\Sigma_{\infty}^{xx}-\Sigma_{\infty}^{xz}).$
(90)

We note that $\Sigma_{\infty}^{xx}=C$ with the matrix $C$ defined in (16) and
that the symmetric part of $\Sigma_{\infty}^{zx}$ and $\Sigma_{\infty}^{xz}$,
respectively, are equivalent to $\Sigma_{\infty}^{zz}-\delta_{\rm noise}\delta
I$. Hence, following (22), we again make the crucial approximation

$(A^{\rm T}A):Z_{t}^{\dagger}\otimes X_{t}^{\dagger}\approx(A^{\rm
T}A):\Sigma_{\infty}^{zx}=(A^{\rm T}A):(\Sigma_{\infty}^{zz}-\delta_{\rm
noise}\delta I)$ (91)

for $d\gg 1$. Let us therefore introduce the shorthand

$\tilde{C}=\Sigma_{\infty}^{zz}-\delta_{\rm noise}\delta I.$ (92)

We also note that

$\tilde{C}=C+\mathcal{O}(\delta).$ (93)

The frequentist analysis from Section 3 delivers

${\rm d}\mu_{t}=\frac{\sigma_{t}}{\gamma}\left((AZ_{t}^{\dagger})^{\rm T}{\rm
d}X_{t}^{\dagger}-\mu_{t}(A^{\rm T}A):\tilde{C}\,{\rm d}t\right)$ (94)

for the conditional mean and

$\frac{\rm d}{{\rm d}t}\sigma_{t}=-\frac{\sigma_{t}^{2}}{\gamma}(A^{\rm
T}A):\tilde{C}$ (95)

for the conditional variance of the random variable $\Theta_{t}$, as defined
by the modified mean-field evolution equations (88). Furthermore, we find that
$\mu_{t}$ still provides an asymptotically unbiased estimator since
$m_{t}=\mathbb{E}^{\dagger}[\mu_{t}]$ satisfies

$\frac{\rm d}{{\rm d}t}m_{t}=\frac{\sigma_{t}}{\gamma}(A^{\rm
T}A):\tilde{C}\,(1-m_{t}).$ (96)

Similarly, the variance, $p_{t}$ of the estimator $\mu_{t}$ satisfies

$\frac{\rm d}{{\rm d}t}p_{t}=\frac{\sigma_{t}}{\gamma}(A^{\rm
T}A):\tilde{C}\,(\sigma_{t}-2p_{t}).$ (97)

In summary, we find that the modified EnKBF mean-field equations (88) behave
exactly as the original equations (17) with the only difference that the
stationary covariance matrix $C=\Sigma_{\infty}^{xx}$ is replaced everywhere
by (92).

Again following [2], given multi-scale observations $X_{t}^{(\epsilon)}$,
$t\geq 0$, we define associated filtered $Z_{t}^{(\epsilon)}$ via

$\frac{\rm d}{{\rm
d}t}Z_{t}^{(\epsilon)}=\frac{1}{\delta}(X_{t}^{(\epsilon)}-Z_{t}^{(\epsilon)})+\delta_{\rm
noise}\sqrt{2}{\rm d}V_{t}^{\dagger}.$ (98)

The intriguing observation is that

$\tilde{J}_{t_{n},t_{n+1}}^{(\epsilon)}:=\int_{t_{n}}^{t_{n+1}}(AZ_{t}^{(\epsilon)})^{\rm
T}{\rm d}X_{t}^{(\epsilon)}$ (99)

converges to

$\tilde{J}^{\dagger}_{t_{n},t_{n+1}}:=\int_{t_{n}}^{t_{n+1}}(AZ_{t}^{\dagger})^{\rm
T}{\rm d}X_{t}^{\dagger}$ (100)

as $\epsilon\to 0$. This simply follows from the fact that both integrals can
be interpreted as standard Riemann–Stieltjes integrals and convergence of
$X_{t}^{(\epsilon)}\to X_{t}^{\dagger}$ and $Z_{t}^{(\epsilon)}\to
Z_{t}^{\dagger}$ as $\epsilon\to 0$ in the standard topology of continuous
functions is sufficient to conclude the convergence of the associated
integrals
$\tilde{J}_{t_{n},t_{n+1}}^{(\epsilon)}\to\tilde{J}_{t_{n},t_{n+1}}^{\dagger}$.
In other words, the extended EnKBF formulation

${\rm
d}\Theta_{t}^{(\epsilon)}&=\frac{\sigma_{t}^{(\epsilon)}}{\gamma}(AZ^{(\epsilon)}_{t})^{\rm
T}{\rm d}I_{t}^{(\epsilon)},\\\ {\rm d}I_{t}^{(\epsilon)}={\rm
d}X_{t}^{(\epsilon)}-\frac{1}{2}\left(\Theta_{t}^{(\epsilon)}+\pi_{t}^{(\epsilon)}[\theta]\right)AX_{t}^{(\epsilon)}{\rm
d}t,\\\ {\rm
d}Z_{t}^{(\epsilon)}=\frac{1}{\delta}(X_{t}^{(\epsilon)}-Z_{t}^{(\epsilon)}){\rm
d}t+\delta_{\rm noise}\sqrt{2}V_{t}^{\dagger}$ (101)

will converge to the correct parameter value $\theta^{\dagger}=1$ as
$t\to\infty$ in the limit $\epsilon\to 0$, that is,

$\lim_{t\to\infty}\lim_{\epsilon\to
0}\Theta_{t}^{(\epsilon)}=\theta^{\dagger}.$ (102)

This statement is in line with the results from [2]. The intriguing point is
that the data filtering approach does not require knowledge of the rough path
correction term implied by (73) while still delivering an unbiased estimator.

##  5 Numerical example

We consider the linear SDE (2) with $\gamma=1$ and

$A=\frac{-1}{2}\left(\begin{array}[]{cc}1&-1\\\ 1&1\end{array}\right).$ (103)

We find that $C=I$ and $A^{\rm T}A=1/2I$. Hence $(A^{\rm T}A):C=1$, and the
posterior variance simply satisfies $\sigma_{t}=\sigma_{0}/(1+\sigma_{0}t)$
according to (44). We set $m_{\rm prior}=0$ and $\sigma_{\rm prior}=4$ for the
Gaussian prior distribution of $\Theta_{0}$, and the observation interval is
$[0,T]$ with $T=6$. We find that $\sigma_{T}=0.16$. Solving (39) for given
$\sigma_{t}$ with initial condition $m_{0}=0$ yields

$m_{t}=1-\frac{\sigma_{t}}{\sigma_{0}}$ (104)

and $m_{T}=0.96$. The corresponding curves are displayed in red in Figure 2.

Figure 2: a)–b): frequentist mean, $m_{t}$ and variance, $p_{t}$, from EnKBF
implementation (20) with step-size $\Delta t=0.06$; c)–d): same results from
EnKBF implementation (30) with inner time-step $\Delta\tau=\Delta t/600$. We
also display the curves arising for $\sigma_{t}$ and $m_{t}$ from the standard
Kalman theory using the approximation (22). Note that the posterior variance,
$\sigma_{t}$, should provide an upper bound on the frequentist uncertainty
$p_{t}$.

We implement the EnKBF schemes (20) and (30) with $t_{n}=n\,\Delta t$. The
inner time-step is $\Delta\tau=10^{-4}$ while $\Delta t=0.06$, that is,
$L=600$. We repeat the experiment $N=10^{4}$ times and compare the outcome
with the predicted mean value of $m_{T}=0.96$ and the posterior variance of
$\sigma_{T}=0.16$ in Figure 2. The differences in the computed time evolutions
of $m_{t}$ and $p_{t}$ are rather minor and support the idea that it is not
necessary to assimilate continuous-time data beyond $\Delta t$. We also find
that the simple prediction (104), based on standard Kalman filter theory, is
not very accurate for this low-dimensional problem ($d=2$). The corresponding
approximation for $\sigma_{t}$ provides, however, a good upper bound for
$p_{t}$.

Figure 3: Same experimental setting as in Figure 2 but with the data now
generated from the multi-scale SDE (58). Again, subsampling the data in
intervals of $\Delta t=0.06$ and high-frequency assimilation with step-size
$\Delta\tau=10^{-4}$ lead to very similar results in terms of their
frequentist means and variances.

We now replace the data generating SDE model (2) by the multi-scale
formulation (58) with $\epsilon=0.01$ and $\beta=2$. This parameter choice
agrees with the one used in [10]. We again find that assimilating the data at
the slow time-scale $\Delta t=0.06$ leads to very similar results obtained
from an assimilation at the fast time-scale $\Delta\tau=10^{-4}$ with the
EnKBF formulation (74), provided the correction term resulting from the
second-order iterated integral (73) is included. See Figure 3. We also
verified numerically that $\Delta t=0.06$ constitutes a nearly optimal step-
size in the sense of making (85) sufficiently small while maintaining
numerical accuracy. For example, reducing the outer step-size to $\Delta
t=0.02$ leads to $h(0.02)-h(0.06)\approx 10$ in (85).

We finally implement the data filtering approaches (88) and (101) with
$\delta=0.1$ using the true signal $X_{t}^{\dagger}$ and its multi-scale
representation $X_{t}^{(\epsilon)}$, respectively. The numerical
implementation with a step-size $\Delta t=\Delta\tau=10^{-4}$ resulted in
approximations $\mu_{t_{n}}$ which converged to the true parameter value
$\theta^{\dagger}=1$ as $t\to T=6$ without the need for including further
corrections terms; as expected from the results in Section 4.2.

##  6 Conclusions

In this follow-up note to [10], we have investigated the impact of
subsampling, data filtering, and high-frequency data assimilation on the
corresponding conditional mean estimators, $\mu_{t}$, both for data generated
from the standard SDE model and a modified multi-scale SDE. A frequentist
analysis supports the basic finding that all three approaches lead to
comparable results provided that the systematic biases due to different
second-order iterated integrals are properly accounted for. While the EnKBF is
relatively easy to analyse and a full rough path approach can be avoided,
extending these results to the nonlinear feedback particle filter [27, 10]
will prove more challenging. Extensions to systems without a strong scale
separation [5, 33] and applications to geophysical fluid dynamics [23, 13] are
also of interest. In this context, the approximation quality of the proposed
estimator (81) and the choice of the step-size $\Delta t$ following (85) (and
potentially $\Delta\tau$) will be of particular interest. Finally, while we
have investigated the univariate parameter estimation problem, a semi-
parametric parametrisation of the drift term $f$ in (1), such as random
feature maps [22], lead to high-dimensional parameter estimation problems and
their statistics [20, 21]. This provides another fertile direction for future
research.

#### Acknowledgements.

SR has been partially funded by Deutsche Forschungsgemeinschaft (DFG) -
Project-ID 318763901 - SFB1294 and Project-ID 235221301 - SFB1114. He would
also like to thank Nikolas Nüsken for many fruitful discussions on the subject
of this paper.

## References

  * Abdulle et al. [2021] A. Abdulle, G. Garegnani, G. A. Pavliotis, A. M. Stuart, and A. Zanoni.  Drift estimation of multiscale diffusions based on filtered data.  _Foundations of Computational Mathematics_ , published online 2021/10/13:in press, 2021.  10.1007/s10208-021-09541-9.
  * Abdulle et al. [2023] A. Abdulle, G. Garegnani, G. Pavliotis, A. Stuart, and A. Zanoni.  Drift estimation of multiscale diffusions based on filtered data.  _Foundations of Computational Mathematics_ , 23:33–84, 2023\.  10.1007/s10208-021-09541-9.
  * Ait-Sahalia et al. [2005] Y. Ait-Sahalia, P. A. Mykland, and L. Zhang.  How often to sample a continuous-time process in the presence of market microstructure noise.  _The Review of Financial Studies_ , 18:351–416, 2005.
  * Amezcua et al. [2014] J. Amezcua, E. Kalnay, K. Ide, and S. Reich.  Ensemble transform Kalman-Bucy filters.  _Q.J.R. Meteor. Soc._ , 140:995–1004, 2014.
  * Arnold [2001] L. Arnold.  Hasselmann’s program revisited: The analysis of stochasticity in deterministic climate models.  In _Stochastic Climate Models_ , pages 141–158. Birkhäuser Basel, 2001.  10.1007/978-3-0348-8287-3.
  * Azencott et al. [2013] R. Azencott, A. Beri, A. Jain, and I. Timofeyev.  Sub-sampling and parametric estimation for multiscale dynamics.  _Communications in Mathematical Sciences_ , 11:939–970, 2013.
  * Bain and Crisan [2009] A. Bain and D. Crisan.  _Fundamentals of Stochastic Filtering_ , volume 60 of _Stoch. Model. Appl. Probab._ Springer, New York, 2009.  10.1007/978-0-387-76896-0.
  * Bálint and Melbourne [2018] P. Bálint and I. Melbourne.  Statistical properties for flows with unbounded roof function, including the Lorenz attractor.  _Journal of Statistical Physics_ , 172:1101–1126, 2018\.  10.1007/s10955-018-2093-y.
  * Bo and Celani [2013] S. Bo and A. Celani.  White-noise limit of nonwhite nonequilibrium processes.  _Physical Review E_ , 88:062150, 2013.  10.1103/PhysRevE.88.062150.
  * Coghi et al. [2023] M. Coghi, T. Nilssen, N. Nüsken, and S. Reich.  Rough McKean–Vlasov dynamics for robust ensemble Kalman filtering.  _Ann. Appl. Probab._ , 33 (6B):5693–5752, 2023.
  * Cotter and Reich [2013] C. Cotter and S. Reich.  Ensemble filter techniques for intermittent data assimilation.  _Radon Ser. Comput. Appl. Math._ , 13:91–134, 2013.  10.1515/9783110282269.91.
  * Crisan et al. [2013] D. Crisan, J. Diehl, P. K. Friz, H. Oberhauser, et al.  Robust filtering: correlated noise and multidimensional observation.  _The Annals of Applied Probability_ , 23:2139–2160, 2013\.
  * Culina et al. [2011] J. Culina, S. Kravtsov, and A. H. Monahan.  Stochastic parameterization schemes for use in realistic climate models.  _Journal of the Atmospheric Sciences_ , 68:284 – 299, 2011\.  10.1175/2010JAS3509.1.
  * Davie [2008] A. M. Davie.  Differential equations driven by rough paths: An approach via discrete approximation.  _Applied Mathematics Research eXpress_ , 2008, 2008.  10.1093/amrx/abm009.  abm009.
  * Diehl et al. [2016] J. Diehl, P. Friz, and H. Mai.  Pathwise stability of likelihood estimators for diffusion via rough paths.  _The Annals of Applied Probability_ , 26:2169–2192, 2016\.  10.1214/15-AAP1143.
  * Evensen [2009] G. Evensen.  _Data assimilation_.  Springer-Verlag, Berlin, second edition, 2009.  ISBN 978-3-642-03710-8.  10.1007/978-3-642-03711-5.
  * Friz and Hairer [2020] P. Friz and M. Hairer.  _A course on rough paths_.  Springer-Verlag, 2020.
  * Friz et al. [2015] P. Friz, P. Gassiat, and T. Lyons.  Physical Brownian motion in a magnetic field as a rough path.  _Transactions of the American Mathematical Society_ , 367:7939–7955, 2015.
  * Garbuno-Inigo et al. [2020] A. Garbuno-Inigo, N. Nüsken, and S. Reich.  Affine invariant interacting Langevin dynamics for Bayesian inference.  _SIAM J. Appl. Dyn. Syst._ , 19:1633–1658, 2020.  10.1137/19M1304891.
  * Ghosal and van der Vaart [2017] S. Ghosal and A. van der Vaart.  _Fundamentals of Nonparametric Bayesian Inference_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2017.  10.1017/9781139029834.
  * Giné and Nickl [2016] E. Giné and R. Nickl.  _Mathematical Foundations of Infinite-Dimensional Statistical Models_.  Cambridge University Press, Cambridge, 2016.  10.1017/CBO9781107337862.
  * Gottwald and Reich [2021] G. A. Gottwald and S. Reich.  Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation.  _Physica D: Nonlinear Phenomena_ , 423:132911, 2021.  ISSN 0167-2789.  https://doi.org/10.1016/j.physd.2021.132911.
  * Hasselmann [1976] K. Hasselmann.  Stochastic climate models Part I. Theory.  _Tellus_ , 28:473–485, 1976.  10.1111/j.2153-3490.1976.tb00696.x.
  * Ikeda and Watanabe [1989] N. Ikeda and S. Watanabe.  _Stochastic differential equations and diffusion pro- cesses_.  North Holland Publishing Company, Amsterdam-New York, 2nd edition, 1989\.
  * Kelly and Melbourne [2017] D. Kelly and I. Melbourne.  Deterministic homogenization for fast-slow systems with chaotic noise.  _Journal of Functional Analysis_ , 272:4063–4102, 2017\.  10.1016/j.jfa.2017.01.015.
  * Kutoyants [2013] Y. A. Kutoyants.  _Statistical inference for ergodic diffusion processes_.  Springer Science & Business Media, 2013.
  * Nüsken et al. [2019] N. Nüsken, S. Reich, and P. J. Rozdeba.  State and parameter estimation from observed signal increments.  _Entropy_ , 21(5):505, 2019.  10.3390/e21050505.
  * Papavasiliou et al. [2009] A. Papavasiliou, G. Pavliotis, and A. Stuart.  Maximum likelihood estimation for multiscale diffusions.  _Stochastic Processes and their Applications_ , 19:3173–3210, 2009.
  * Pathiraja [2020] S. Pathiraja.  $L^{2}$ convergence of smooth approximations of stochastic differential equations with unbounded coefficients, 2020.  arXiv:2011.13009.
  * Reich [2023] S. Reich.  Frequentist perspective on robust parameter estimation using the ensemble Kalman filter.  In B. Chapron, D. Crisan, D. Holm, E. Mémin, and A. Radomska, editors, _Stochastic Transport in Upper Ocean Dynamics. STUOD 2021. Mathematics of Planet Earth_ , volume 10, pages 237–258, Cham., 2023. Springer.  10.1007/978-3-031-18988-3_15.
  * Reich and Rozdeba [2020] S. Reich and P. Rozdeba.  Posterior contraction rates for non-parametric state and drift estimation.  _Foundation of Data Science_ , 2:333–349, 2020.  10.3934/fods.2020016.
  * Sirignano and Spiliopoulos [2017] J. Sirignano and K. Spiliopoulos.  Stochastic gradient descent in continuous time.  _SIAM J. Financial Math._ , 8:933–961, 2017.  10.1137/17M1126825.
  * Wouters and Gottwald [2019] J. Wouters and G. A. Gottwald.  Stochastic model reduction for slow-fast systems with moderate time scale separation.  _Multiscale Modeling & Simulation_, 17:1172–1188, 2019\.