Robust parameter estimation using the ensemble Kalman filter

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.

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,$$

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 (2b), 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 (2a) 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 (2) 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,$$

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

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$$

in our notation, where $\alpha_{t}>0$ denotes the learning rate. We note that (2.1) shares with (2) the gain times innovation structure. However, while (2) approximates the Bayesian inference problem, formulation (2.1) 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 (2) also applies in simplified form to (2.1). We also point out that (2) 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}.$$

It is straightforward to combine this time discretisation with the Monte-Carlo approximation (2) 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 (2) and its appropriate time discretisation. Here we follow the Itô/Euler–Maruyama formulation of the data-driven term in (2),

$$\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 (2) 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,$$

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 (2.1):

$${\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.$$

We find that the learning rate $\alpha_{t}$ takes the role of the variance $\sigma_{t}$ in (2). 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 (2) replaces (2a) 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})$$

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),$$

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 (2) 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 (2) 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 (2) 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,$$

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],$$

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 (2.4) 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 (2) 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 (2), 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},$$

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 (4), 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,$$

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

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)}]$$

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.$$

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$$

to (64). Please note that it follows from (4) 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, (4) 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 (4) when taking the limit $\epsilon\to 0$.

Remark 4.5.

The analogous adaptation of (4) to the gradient descent formulation (2.4) 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).$$

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 (4), 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 (4) 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 (4). 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 (4.2f) 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 (4) 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 (4) 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 (4). 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 (4) 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 (2), 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,$$

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},$$

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

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 (4.2). 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 (4.2) behave exactly as the original equations (2) 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}$$

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 (4). 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 (4) 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 (4), 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 (4.2) and (4.2) 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.