# Coupled Stochastic-Statistical Equations for Filtering Multiscale Turbulent
Systems

Di Qia and Jian-Guo Liub

###### Abstract

We present a new strategy for filtering high-dimensional multiscale systems
characterized by high-order non-Gaussian statistics using observations from
leading-order moments. A closed stochastic-statistical modeling framework
suitable for systematic theoretical analysis and efficient numerical
simulations is designed. Optimal filtering solutions are derived based on the
explicit coupling structures of stochastic and statistical equations subject
to linear operators, which satisfy an infinite-dimensional Kalman-Bucy filter
with conditional Gaussian dynamics. To facilitate practical implementation, we
develop a finite-dimensional stochastic filter model that approximates the
optimal filter solution. We prove that this approximating filter effectively
captures key non-Gaussian features, demonstrating consistent statistics with
the optimal filter first in its analysis step update, then at the long-time
limit guaranteeing stable convergence to the optimal filter. Finally, we build
a practical ensemble filter algorithm based on the approximating filtering
model, which enables accurate recovery of the true model statistics. The
proposed modeling and filtering strategies are applicable to a wide range
challenging problems in science and engineering, particularly for statistical
prediction and uncertainty quantification of multiscale turbulent states.

aDepartment of Mathematics, Purdue University, 150 North University Street,
West Lafayette, IN 47907, USA

bDepartment of Mathematics and Department of Physics, Duke University, Durham,
NC 27708, USA

##  1 Introduction

Complex turbulent phenomena are widely observed in various science and
engineering systems [51, 14, 42]. Such systems are characterized by a wide
spectrum of nonlinearly coupled spatiotemporal scales and high degrees of
inherent internal instability [37, 21]. A probabilistic formulation containing
highly non-Gaussian statistics is required to quantify uncertainties in the
high-dimensional turbulent state [44, 38]. Traditional ensemble approaches
using a particle system representation to approximate the probability
evolution quickly become computationally prohibitive since a sufficiently
large sample size is necessary to capture the extreme non-Gaussian outliers
even for relatively low dimensional systems [35, 55]. As a result, rigorous
analysis often becomes intractable and direct numerical simulations are likely
to be expensive and inaccurate.

Filtering strategies [50, 33, 52] have long been used for finding the optimal
probability estimates of a stochastic state with large uncertainty based on
partial and noisy observation data. Applications to improve the probability
forecast can be found in various practical dynamical systems [23, 29, 3, 5].
In predicting nonlinear turbulent signals, ensemble Kalman filters [18, 26,
19] as well as the related particle methods [13, 15, 7] provide an effective
tool for state and parameter estimations. Filtering theories [58, 8, 45] and
corresponding numerical solutions [32, 4, 9] for general nonlinear systems
have been investigated through different approaches. Despite wide applications
[17, 40, 12], difficulties persist for accurate statistical forecast of
turbulent states especially when non-Gaussian features are present in the
target probability distribution. Conventional ensemble-based approaches often
suffer inherent difficulties in estimating the crucial higher-order moment
statistics and maintaining stable prediction with finite number of particles
[53, 24, 54]. On the other hand, in many situations observations of the
statistical states, such as the mean and covariance, are easier to extract
from various sources of data [39]. Therefore, a promising research direction
is to propose new ensemble filtering models to recover the highly non-Gaussian
probability distributions using statistical observations from the leading-
order mean and covariance [1].

###  1.1 General problem setup

We start with a general mathematical formulation [38] modeling a high-
dimensional stochastic state variable $u_{t}\in\mathbb{R}^{d}$ involving
nonlinear multiscale interactions satisfying the following stochastic
differential equation (SDE)

$\frac{\mathrm{d}u_{t}}{\mathrm{d}t}=\Lambda
u_{t}+B\left(u_{t},u_{t}\right)+F_{t}+\sigma_{t}\dot{W}_{t}.$ (1.1)

On the right hand side of the above equation, the linear operator,
$\Lambda:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, represents linear
dispersion and dissipation effects. The nonlinear effect in the dynamical
system is introduced via a bilinear quadratic operator,
$B:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. The system is
subject to external forcing effects that are decomposed into a deterministic
component, $F_{t}$, and a stochastic component represented by a Gaussian white
noise $W_{t}\in\mathbb{R}^{s}$ with $\sigma_{t}\in\mathbb{R}^{d\times s}$. The
model emphasizes the important role of quadratic interactions through
$B\left(u,u\right)$. This typical structure is inherited from a finite-
dimensional truncation of the continuous system, for example, a spectral
projection of the nonlinear advection in fluid model [34]. Many realistic
systems with wide applications [27, 48, 22] can be categorized in the general
dynamical equation (1.1).

The evolution of the model state $u_{t}$ depends on the sensitivity to the
randomness in initial conditions and external stochastic effects, which will
be further amplified in time by the inherent internal instability due to the
nonlinear coupling term [37, 49]. Assuming that the stochastic solution
satisfies a continuous probability density function (PDF), the time evolution
of the PDF $p_{t}$ is given by the associated Fokker-Planck equation (FPE)

$\frac{\partial p_{t}}{\partial
t}=\mathcal{L}_{\mathrm{FP}}p_{t}\coloneqq-\nabla_{u}\cdot\left[\left(\Lambda
u+B\left(u,u\right)+F_{t}\right)p_{t}\right]+\frac{1}{2}\nabla_{u}\cdot\left[\nabla_{u}\cdot\left(\sigma_{t}\sigma_{t}^{\intercal}p_{t}\right)\right],\;p_{t=0}=p_{0},$
(1.2)

where $\mathcal{L}_{\mathrm{FP}}$ represents the Fokker-Planck operator with
$\nabla\cdot\left(\nabla\cdot
A\right)=\sum_{k,l}\frac{\partial^{2}A_{kl}}{\partial u_{k}\partial u_{l}}$.
However, it remains a challenging task for directly solving the FPE (1.2) as a
high-dimensional PDE. As an alternative approach, ensemble forecast by
tracking the Monte-Carlo solutions estimates the essential statistics through
empirical averages among a group of samples drawn i.i.d. from the initial
distribution $u^{\left(i\right)}\left(0\right)\sim p_{0}$ at the starting time
$t=0$. The PDF solution $p_{t}\left(u\right)$ and the associated statistical
expectation of any test function $\varphi\left(u\right)$ at each time instant
$t>0$ are then approximated by the empirical ensemble representation

$p_{t}\left(u\right)\simeq
p_{t}^{N}\left(u\right)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\delta\left(u-u_{t}^{\left(i\right)}\right),\quad\mathbb{E}_{p_{t}}\varphi\left(u\right)\simeq\mathbb{E}_{p_{t}^{N}}\varphi\left(u\right)=\frac{1}{N}\sum_{i=1}^{N}\varphi\left(u_{t}^{\left(i\right)}\right),$
(1.3)

In practice, large errors will be introduced in the above empirical
estimations since only a finite sample approximation is available in modeling
the probability distribution and statistics in a high dimensional system.

It is expected that the model errors in the predicted PDF of the stochastic
states from the finite ensemble simulation can be effectively corrected
through the available observation data. In designing new filtering strategies,
we propose to use statistical observations from mean and covariance to improve
the accuracy and stability in the finite ensemble forecast using the finite
ensemble approximation. However, the general formulation (1.1) as well as the
associated FPE (1.2) becomes inconvenient to use since all the multiscale
stochastic processes in the equation are mixed together containing all high-
order statistics. The main goal of this paper is thus to develop a systematic
modeling framework with strategies to accurately capture the (potentially
highly non-Gaussian) PDF $p_{t}$ with the help of statistical measurements in
the leading moments.

###  1.2 Overview of the paper

In this paper, we study nonlinear filtering of the general multiscale
turbulent system (1.1) according to the two-step filtering procedure. In the
first forecast step, a new modeling framework is designed combining
statistical equations (2.3) for the leading moments and stochastic equations
(2.5) for capturing higher-order statistics. In the second analysis step, we
propose the observation process from the statistical equations of the leading-
order mean and covariance to improve the probability prediction of the
stochastic model. The statistical equations for observation processes (3.1b)
then become linearly dependent on the PDF of stochastic state, while the
evolution of the PDF is also subject to conditional linear dynamics (3.1a).
Therefore, optimal filtering equations can be explicitly derived based on the
linearity in both signal and observation processes. Finally, the expensive
optimal filter equations are approximated by an efficient filtering model
demonstrating equivalent statistics according to the most important
observation functions (3.4).

The new multiscale nonlinear filtering model is constructed under the
following step-by-step procedure, which will finally be combined to build a
practical ensemble filtering strategy able to recover non-Gaussian PDFs:

  * •

First, we propose a precise coupled stochastic-statistical formulation (2.7)
as the forecast model for multiscale turbulence: the stochastic dynamics will
serve as the signal process in filtering including high-order non-Gaussian
features, while the reinforced statistical equations provide the observation
process;

  * •

Second, the optimal filtering problem is formulated based on the particular
coupling structure in the multiscale stochastic-statistical model: the optimal
filter solution (3.12) is given based on the PDF of the stochastic state,
combined with the mean and covariance as a natural choice of the observed
state;

  * •

Third, an effective statistical filtering model (3.16) is developed
approximating the optimal filter solution in key leading-order statistics: an
equivalent McKean-Vlasov SDE containing higher moments feedbacks is adapted
from the optimal filter solution;

  * •

Last, discrete ensemble filtering schemes (5.2) are constructed as a particle
approximation of the statistical filtering model: the filter SDE is
approximated by interacting particles and the statistical moments are computed
by empirical ensemble averages.

The coupled stochastic-statistical model (2.7) by itself can serve as an
effective tool for statistical forecasts and uncertainty quantification.
Efficient computational algorithms have been developed [47, 48] that
demonstrate high skill in capturing crucial non-Gaussian phenomena such as
extreme events and fat-tailed PDFs. Then combined with the observation data,
the resulting filtering McKean-Vlasov SDE (3.16) is only implicitly dependent
on the probability distribution, which can be computed directly from the
statistical equations. This enables efficient computational strategies to
effectively improve the accuracy and stability in capturing high-order non-
Gaussian features based on only observation from the lower moments. We
illustrate the main steps in building effective models for capturing
probability distributions in Fig. 1.1.

Predicting Probability Distributions of Multiscale Turbulent States (1.1) A
Coupled Stochastic-Statistical Closure Modeling Framework (2.7) Optimal
Filtering Solution in Infinite-Dimensional Equations (3.12) state
decomposition filtering density with statistical observations A Statistical
Filtering Model (3.16) approximating Optimal Filter Moments statistical
consistency in observation functions ensemble-based approach Figure 1.1:
Diagram illustrating main ideas in constructing the modeling framework

Our main focus of this paper is to provide a detailed description of the
forecast and analysis step of the filtering method so that various
applications can directly follow based on this highly adaptive general
framework. We provide detailed discussion in each component of the proposed
model and methods: i) we show in Proposition 2 that the highly tractable
closure model (2.7) demonstrates consistent statistics as the original system
(1.1) which plays a fundamental role in the construction of the filtering
methods; ii) the filter solution (3.12) is shown optimal in the mean square
sense by exploiting the conditional Gaussian structure of the forward
equation, and the approximation model is shown to recover the same key
statistics during the analysis step update in Theorem 7; iii) the long-time
convergence in statistics to the optimal filter is demonstrated in Theorem 9
concerning the complete filtering procedure using the statistical filtering
model (3.16); and iv) the ensemble approximation with discrete time step is
shown to recover to the true statistics at the large ensemble limit.

The structure of the paper is organized as follows: In Section 2, we first set
up a statistically consistent formulation for the general multiscale system
(1.1) that is suitable for the construction of the statistical filtering
models. Section 3 shows the major ideas of finding the optimal filter solution
and constructing approximate filtering models. Long-time convergence and
stability of the approximation filter model is then discussed in Section 4.
Combining each component of the ideas, the final ensemble statistical
filtering model is developed in Section 5 Finally, a summary of this paper is
given in Section 6.

##  2 A statistically consistent modeling framework for multiscale dynamical
systems

In this section, we first introduce a new formulation for the original system
(1.1) using an explicit macroscopic and microscopic decomposition of the
multiscale state. In particular, we show that the new formulation provides
consistent statistics with the original system in the explicit leading moment
equations, as well as higher-order statistics in the stochastic equation. In
addition, the new formulation also enjoys a more tractable dynamical structure
to be adapted to the filtering framework.

###  2.1 The statistical and stochastic equations as macroscopic and
microscopic states

In order to identify the detailed multiscale interactions in the general
system (1.1), we decompose the random model state $u_{t}$ into the composition
of a statistical mean $\bar{u}_{t}$ and stochastic fluctuations
$u_{t}^{\prime}$ in a finite-dimensional representation under an orthonormal
basis $\left\\{\hat{v}_{k}\right\\}_{k=1}^{d}$ with
$\hat{v}_{k}\cdot\hat{v}_{l}=\delta_{kl}$

$u_{t}=\bar{u}_{t}+u_{t}^{\prime}=\sum_{k=1}^{d}\bar{u}_{k,t}\hat{v}_{k}+\sum_{k=1}^{d}Z_{k,t}\hat{v}_{k},\quad\mathrm{with}\;\bar{u}_{k,t}=\hat{v}_{k}\cdot\bar{u}_{t},\;Z_{k,t}=\hat{v}_{k}\cdot\left(u_{t}-\bar{u}_{t}\right).$
(2.1)

Above, $\bar{u}_{t}\in\mathbb{R}^{d}$ represents the dominant largest-scale
structure (for example, the zonal jets in geophysical turbulence or the
coherent radial flow in fusion plasmas); and
$Z_{t}=\left[Z_{1,t},\cdots,Z_{d,t}\right]^{\intercal}\in\mathbb{R}^{d}$ are
stochastic coefficients projected on each eigenmode $\hat{v}_{k}$, whose
randomness illustrates the uncertainty in each single scale of
$u_{t}^{\prime}$. In particular, we will show that the dynamics of the
stochastic modes $Z_{t}$ contain nonlinear interactions among a large number
of coupled multiscale fluctuations, which demonstrate the charactering feature
in turbulent systems demonstrating forward and backward energy cascades [21,
37].

####  2.1.1 Statistical equations for the macroscopic states

First, we define the leading-order mean and covariance according to the state
decomposition (2.1)

$\bar{u}_{t}=\mathbb{E}_{p_{t}}\left[u\right]\coloneqq\int_{\mathbb{R}^{d}}up_{t}\left(\mathrm{d}u\right),\quad
R_{kl,t}=\mathbb{E}_{p_{t}}\left[Z_{k,t}Z_{l,t}\right]\coloneqq\int_{\mathbb{R}^{d}}\hat{v}_{k}\cdot\left(u-\bar{u}_{t}\right)\left(u-\bar{u}_{t}\right)\cdot\hat{v}_{l}p_{t}\left(\mathrm{d}u\right),\;1\leq
k,l\leq d,$ (2.2)

where the expectations are with respect to the law $p_{t}$ from the PDF
solution of (1.2). Above, statistical mean $\bar{u}_{t}$ and covariance
$R_{t}$ can represent the _macroscopic_ physical quantities that are easiest
to achieve from direct measurements. The dynamical equations for the mean and
covariance can be derived by the following statistical equations for $1\leq
k,l\leq d$

$\displaystyle\frac{\mathrm{d}\bar{u}_{k,t}}{\mathrm{d}t}=$
$\displaystyle\>\hat{v}_{k}\cdot\left[\Lambda\bar{u}_{t}+B\left(\bar{u}_{t},\bar{u}_{t}\right)\right]+\sum_{m,n=1}^{d}\gamma_{kmn}\mathbb{E}_{p_{t}}\left[Z_{m,t}Z_{n,t}\right]+\hat{v}_{k}\cdot
F_{t},$ (2.3a) $\displaystyle\frac{\mathrm{d}R_{kl,t}}{\mathrm{d}t}=$
$\displaystyle\sum_{m=1}^{d}\left[L_{km}\left(\bar{u}_{t}\right)R_{t,ml}+R_{t,km}L_{lm}\left(\bar{u}_{t}\right)\right]+Q_{t,kl},$
(2.3b) $\displaystyle+$
$\displaystyle\sum_{m,n=1}^{d}\gamma_{kmn}\mathbb{E}_{p_{t}}\left[Z_{m,t}Z_{n,t}Z_{l,t}\right]+\gamma_{lmn}\mathbb{E}_{p_{t}}\left[Z_{m,t}Z_{n,t}Z_{k,t}\right].$

Above, we define the nonlinear coupling coefficients
$\gamma_{kmn}=\hat{v}_{k}\cdot B\left(\hat{v}_{m},\hat{v}_{n}\right)$, and the
white noise coefficient is defined as
$\Sigma_{t}=\left[\left(\hat{v}_{1}^{\intercal}\sigma_{t}\right)^{\intercal},\cdots,\left(\hat{v}_{d}^{\intercal}\sigma_{t}\right)^{\intercal}\right]^{\intercal}$
with $Q_{t}=\Sigma_{t}\Sigma_{t}^{\intercal}\in\mathbb{R}^{d\times d}$. The
mean-fluctuation coupling operator
$L\left(\bar{u}_{t}\right)\in\mathbb{R}^{d\times d}$ dependent on the
statistical mean state $\bar{u}_{t}$ is defined as

$L_{kl}\left(u\right)=\hat{v}_{k}\cdot\left[\Lambda\hat{v}_{l}+B\left(u,\hat{v}_{l}\right)+B\left(\hat{v}_{l},u\right)\right].$
(2.4)

Notice that the right hand side of (2.3b) involves the fluctuation modes
$Z_{t}$ defined from $u_{t}$ in (2.1), then the expectations on third moments
are taken w.r.t. the PDF $p_{t}$. Therefore, the resulting statistical
equations (2.3) are not closed and need to be combined with the FPE (1.2) to
achieve a closed formulation for the leading-order mean and covariance in the
nonlinear system.

####  2.1.2 Stochastic equations for the microscopic processes

Second, we introduce the SDE describing the time evolution of the multiscale
stochastic processes $Z_{t}$ as the _microscopic_ state consisting of the many
subscale fluctuations

$\mathrm{d}Z_{t}=L\left(\bar{u}_{t}\right)Z_{t}\mathrm{d}t+\varGamma\left(Z_{t}Z_{t}^{\intercal}-R_{t}\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}W_{t}.$
(2.5)

Above, $L\left(\bar{u}_{t}\right)$ is the same mean-fluctuation coupling
operator defined in (2.4) involving the statistical mean $\bar{u}_{t}$, and we
define the quadratic coupling operator $\varGamma:\mathbb{R}^{d\times
d}\rightarrow\mathbb{R}^{d}$ as a linear combination of the entries of the
input matrix $R\in\mathbb{R}^{d\times d}$ describing the nonlinear multiscale
coupling involving the covariance $R_{t}$

$\varGamma_{k}\left(R\right)=\sum_{m,n=1}^{d}\hat{v}_{k}\cdot
B\left(\hat{v}_{m},\hat{v}_{n}\right)R_{mn}.$ (2.6)

Similar to the statistical equations, the dynamics on the right hand side of
(2.5) is linked to the macroscopic quantities $\bar{u}_{t}$ and $R_{t}$, which
in tern requires additional information of $p_{t}$ from the original model
state $u_{t}$. This makes the stochastic equation also unclosed requiring
additional information from the PDF solution (1.2) from the original system.
The derivations of the above statistical and stochastic equations (2.3) and
(2.5) from the original general formulation (1.1) are shown in the proof of
Proposition 2.

###  2.2 A coupled stochastic-statistical closure model with explicit higher-
order feedbacks

Combining the ideas in the closely related stochastic equation (2.5) and the
statistical equations (2.3), we propose a _statistically consistent
stochastic-statistical closure model_ based on the following self-consistent
coupling of the microscopic stochastic processes $Z_{t}$ and the macroscopic
statistics $\bar{u},R_{t}$

$\displaystyle\mathrm{d}Z_{t}=$
$\displaystyle\>L\left(\bar{u}_{t}\right)Z_{t}\mathrm{d}t+\varGamma\left(Z_{t}Z_{t}^{\intercal}-R_{t}\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}W_{t},$
(2.7) $\displaystyle\frac{\mathrm{d}\bar{u}_{t}}{\mathrm{d}t}=$
$\displaystyle\>M\left(\bar{u}_{t}\right)+Q_{m}\left(\mathbb{E}\left[Z_{t}\otimes
Z_{t}\right]\right)+F_{t},$
$\displaystyle\frac{\mathrm{d}R_{t}}{\mathrm{d}t}=$
$\displaystyle\>L\left(\bar{u}_{t}\right)R_{t}+R_{t}L\left(\bar{u}_{t}\right)^{\intercal}+Q_{v}\left(\mathbb{E}\left[Z_{t}\otimes
Z_{t}\otimes Z_{t}\right]\right)$
$\displaystyle+\Sigma_{t}\Sigma_{t}^{\intercal}+\epsilon^{-1}\left(\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]-R_{t}\right).$

Above, the expectations are all w.r.t. the PDF $\rho_{t}$ of the stochastic
states $Z_{t}$. In the first moment equation for $\bar{u}_{t}$, with a bit
abuse of notation, we denote
$\bar{u}_{t}=\left[\bar{u}_{1,t},\cdots,\bar{u}_{d,t}\right]^{\intercal}\in\mathbb{R}^{d}$
with each component $\bar{u}_{k,t}=\bar{u}_{t}\cdot\hat{v}_{k}$,
$M=\left[M_{1},\cdots,M_{d}\right]^{\intercal}\in\mathbb{R}^{d}$ where
$M_{k}\left(\bar{u}_{t}\right)=\sum_{p,q}\hat{v}_{k}\cdot\left[\Lambda\hat{v}_{p}\bar{u}_{p,t}+B\left(\hat{v}_{p},\hat{v}_{q}\right)\bar{u}_{p,t}\bar{u}_{q,t}\right]$
for $1\leq k\leq d$ and $F_{t}=\left[\hat{v}_{1}\cdot
F_{t},\cdots,\hat{v}_{d}\cdot F_{t}\right]$; in the second moment equation for
$R_{t}\in\mathbb{R}^{d\times d}$, the operator
$L\left(\bar{u}_{t}\right)\in\mathbb{R}^{d\times d}$ indicates mean-
fluctuation interactions defined in (2.4); and in the stochastic equation for
$Z_{t}$, $\varGamma:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}^{d}$ is the
quadratic coupling operator defined in (2.6). The two higher-moment feedbacks
for the mean and covariance, $Q_{m},Q_{v}$, related to the second and third
moments of $Z_{t}$ respectively, are defined as

$\displaystyle Q_{m,k}\left(\mathbb{E}\left[Z_{t}\otimes Z_{t}\right]\right)=$
$\displaystyle\sum_{p,q=1}^{d}\gamma_{kpq}\mathbb{E}\left[Z_{p,t}Z_{q,t}\right],$
(2.8) $\displaystyle Q_{v,kl}\left(\mathbb{E}\left[Z_{t}\otimes Z_{t}\otimes
Z_{t}\right]\right)=$
$\displaystyle\sum_{p,q=1}^{d}\gamma_{kpq}\mathbb{E}\left[Z_{p,t}Z_{q,t}Z_{l,t}\right]+\gamma_{lpq}\mathbb{E}\left[Z_{p,t}Z_{q,t}Z_{k,t}\right],$

for $1\leq k,l\leq d$ with coupling coefficients
$\gamma_{kmn}=\hat{v}_{k}\cdot B\left(\hat{v}_{m},\hat{v}_{n}\right)$. Above,
$Q_{m}$ models the feedback in the mean equation due to the second moments
$\mathbb{E}\left[Z_{t}\otimes Z_{t}\right]$, and $Q_{v,kl}$ is the symmetric
feedback in the covariance equation due to all the third moments
$\mathbb{E}\left[Z_{t}\otimes Z_{t}\otimes Z_{t}\right]$. Notice that
$Q_{m},Q_{v}$ can be both viewed as linear operators w.r.t. $\rho_{t}$.
Besides, in the second moment equation an additional relaxation term with a
parameter $\epsilon>0$ is introduced. This term will not modify the model
dynamics with statistical consistency, but is playing a crucial role as a
‘reinforcement’ term in maintaining stable numerical performance (for example,
see [47]) especially with highly instability induced by the strong mean-
fluctuation coupling from $L\left(\bar{u}_{t}\right)$.

Different from the inherently coupled stochastic and statistical equations
(2.5) and (2.3) through the PDF $p_{t}$ of the original model state $u_{t}$,
the new closure model (2.7) provides a clean self-consistent formulation for
tractable theoretical analysis and direct numerical implementations. A new PDF
$\rho_{t}$ of the stochastic process $Z_{t}$ is introduced to close the
system. The statistical states $\bar{u}_{t},R_{t}$ are first treated as new
individual processes subject to higher-order moments w.r.t. $\rho_{t}$. Then,
the microscopic stochastic equation for $Z_{t}$ models the high-dimensional
multiscale process with explicit dependence on the macroscopic states
$\bar{u}_{t}$ and $R_{t}$. Thus, the whole system is closed by the law of the
stochastic state $Z_{t}$ itself. No additional information about the
intractable PDF $p_{t}$ will be needed for solving the FPE of the original
system.

In the rest part of this section, we built precise link between the new
closure model (2.7) and the coupled equations (2.5) and (2.3) from the
original system. First, the following lemma provides the self-consistency in
the leading moments of the stochastic modes $Z_{k}$ and statistical states
$\bar{u}_{t},R_{t}$ in (2.7).

######  Lemma 1.

With consistent initial conditions $\mathbb{E}\left[Z_{0}\right]=0$ and
$\mathbb{E}\left[Z_{0}Z_{0}^{\intercal}\right]=R_{0}$, the leading moments of
the stochastic modes $Z_{t}$ of the closure model (2.7) satisfy that for all
$t>0$

$\mathbb{E}\left[Z_{t}\right]=0,\quad\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]=R_{t},$
(2.9)

where the expectation is taken w.r.t. the PDF $\rho_{t}$ of $Z_{t}$, and
$R_{t}$ is the solution of the second moment equation in (2.7).

The proof can be found through the direct application of Itô’s formula and we
put the detailed proof in Appendix B. Lemma 1 demonstrates that the law of
$Z_{t}$ indeed recovers the same covariance of the fluctuation modes, while it
also contains more information in the higher-order statistics. Furthermore,
the closed coupled stochastic-statistical model generates the same statistical
solution as the original system (1.1). The following proposition describes the
statistical consistency between the coupled model (2.7) and the original
system (1.1).

######  Proposition 2.

Assume that the solution $u_{t}$ of the system (1.1) has a continuous PDF
$p_{t}$, and the solution $\left\\{\bar{u}_{t},R_{t};Z_{t}\right\\}$ of the
stochastic-statistical closure model (2.7) has the continuous PDF $\rho_{t}$
for $Z_{t}$ together with the deterministic solutions for $\bar{u}_{t}$ and
$R_{t}$. Then from the same initial conditions, the two models give the same
statistical solution, that is, for all $t>0$

$\mathbb{E}_{p_{t}}\left[u_{t}\right]=\bar{u}_{t},\quad\mathbb{E}_{p_{t}}\left[u_{t}^{\prime}u_{t}^{\prime\intercal}\right]=R_{t},$
(2.10)

where $u_{t}^{\prime}=u_{t}-\mathbb{E}_{p_{t}}\left[u_{t}\right]$.
Furthermore, for both models with any function $\varphi\in
C^{2}\left(\mathbb{R}^{d}\right)$ we have

$\mathbb{E}_{p_{t}}\left[\varphi\left(u_{t}^{\prime}\right)\right]=\mathbb{E}\left[\varphi\left(\sum_{k=1}^{d}Z_{k,t}\hat{v}_{k}\right)\right].$
(2.11)

The proof of Proposition 2 can be found in Appendix B through detailed
computation of each moments of (1.1) compared with that of (2.7). Notice that
the left hand sides of (2.10) and (2.11) consist of the statistics requiring
solving the PDF $p_{t}$ of the original system, while the right hand sides are
purely from the PDF $\rho_{t}$ of the closure model. This confirms the direct
link between the new stochastic-statistical closure model (2.7) to the
statistics in the original fully coupled multiscale system. From this detailed
formulation, quantification of non-Gaussian statistics in $Z_{t}$ relies on
the accurate estimation of the leading-order mean and covariance, which can be
assisted from the observation data. This leads to effective filtering methods
for improved probability forecast for $Z_{t}$ that will be discussed next in
Section 3.

######  Remark.

The second-order closure model (2.7) provides explicit statistical dynamics
for the mean and covariance. We can also propose a first-order closure model
involving only the statistical mean equation coupled with the McKean-Vlasov
SDE

$\displaystyle\mathrm{d}Z_{t}$
$\displaystyle=L\left(\bar{u}_{t}\right)Z_{t}\mathrm{d}t+\varGamma\left(Z_{t}Z_{t}^{\intercal}-\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}W_{t},$
(2.12) $\displaystyle\frac{\mathrm{d}\bar{u}_{t}}{\mathrm{d}t}$
$\displaystyle=M\left(\bar{u}_{t}\right)+\sum_{p,q}\mathbb{E}\left[Z_{p,t}Z_{q,t}\right]B\left(\hat{v}_{p},\hat{v}_{q}\right)+F_{t},$

However, the above equations may introduce an unbalanced approximation to the
mean and fluctuation interactions since the dynamics of the SDE for $Z_{t}$
will directly involve expectation w.r.t. its law $\rho_{t}$. Thus, this model
will be prone to numerical errors in practical applications. Still, (2.12) can
serve as an intermediate model for the analysis of filtering schemes discussed
in the next section.

###  2.3 The stochastic McKean-Vlasov equation as a multiscale interacting
system

From the closed stochastic-statistical formulation (2.7), the SDE for $Z_{t}$
is given by a stochastic McKean-Vlasov equation depending on its own
probability distribution $\rho_{t}$. In particular, the resulting McKean-
Vlasov SDE can be viewed naturally as the mean-field limit of the ensemble
approximation of $N$ individual trajectories

$\mathrm{d}Z_{t}^{\left(i\right)}=L\left(\bar{u}_{t}^{N}\right)Z_{t}\mathrm{d}t+\varGamma\left(Z_{t}^{\left(i\right)}Z_{t}^{\left(i\right)\intercal}-R_{t}^{N}\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}W_{t}^{\left(i\right)},\quad
i=1,\cdots,N,$ (2.13)

where $\left\\{W_{t}^{\left(i\right)}\right\\}_{i=1}^{N}$ are independent
white noise processes, and the initial samples
$\left\\{Z_{0}^{\left(i\right)}\right\\}_{i=1}^{N}$ are drawn from the same
initial measure $\rho_{0}$ of $Z_{t}$. Notice that the ensemble members
$Z_{t}^{\left(i\right)}$ as interacting particles are not evolving
independently with each other, but are coupled through feedbacks of the
leading-order statistics $\bar{u}_{t}^{N}$ and $R_{t}^{N}$ according to the
statistical equations

$\displaystyle\frac{\mathrm{d}\bar{u}_{t}^{N}}{\mathrm{d}t}=$
$\displaystyle\>M\left(\bar{u}_{t}^{N}\right)+Q_{m}\left(\mathbb{E}^{N}\left[\mathbf{Z}_{t}\otimes\mathbf{Z}_{t}\right]\right)+F_{t},$
(2.14) $\displaystyle\frac{\mathrm{d}R_{t}^{N}}{\mathrm{d}t}=$
$\displaystyle\>L\left(\bar{u}_{t}^{N}\right)R_{t}^{N}+R_{t}^{N}L^{\intercal}\left(\bar{u}_{t}^{N}\right)+Q_{v}\left(\mathbb{E}^{N}\left[\mathbf{Z}_{t}\otimes\mathbf{Z}_{t}\otimes\mathbf{Z}_{t}\right]\right)$
$\displaystyle+\Sigma_{t}\Sigma_{t}^{\intercal}+\epsilon^{-1}\left(\mathbb{E}^{N}\left[\mathbf{Z}_{t}\mathbf{Z}_{t}^{\intercal}\right]-R_{t}^{N}\right).$

Above, the expectations are computed through the empirical average of the
interacting particles
$\mathbf{Z}_{t}=\left\\{Z_{t}^{\left(i\right)}\right\\}_{i=1}^{N}$ from the
ensemble simulation

$\mathbb{E}^{N}\left[\varphi\left(\mathbf{Z}_{t}\right)\right]=\frac{1}{N}\sum_{i=1}^{N}\varphi\left(Z_{t}^{\left(i\right)}\right).$

The coupled ensemble approximation equations (2.13) and (2.14) have advantages
in practical applications. Unlike the general McKean-Vlasov SDEs [41], (2.13)
avoids the direct inclusion of the PDF of $Z_{t}$, which is very difficult to
approximate accurately from finite particles. Instead, the PDF only implicitly
enters the mean and covariance equations (2.14) through the computation of
higher moments. Effective computational algorithms with consistent statistics
then can be proposed (such as using the efficient random batch methods [47,
48]) for the straightforward ensemble model approximation. Besides in
practical computation, the relaxation term in $R_{t}^{N}$ provides additional
restoring forcing as a correction term to numerical errors with finite sample
approximation to reinforce stable dynamics and consistent statistics
especially in the case where internal instability is involved.

In particular, it is well-know [25] that the empirical measure converges
weakly to the true distribution, $\rho_{t}^{N}\rightarrow\rho_{t}$, as well as
the leading-order statistics in (2.13),
$\bar{u}^{N}\rightarrow\bar{u},R_{t}^{N}\rightarrow R_{t}$, as
$N\rightarrow\infty$ under relatively weak assumptions (see Proposition 12 in
Section 5). The dynamical equation for the continuous density function
$\rho_{t}$ of $Z_{t}$ is given by the corresponding equation

$\frac{\partial\rho_{t}}{\partial
t}=\mathcal{L}_{t}^{*}\left(\bar{u}_{t},R_{t}\right)\rho_{t}\coloneqq-\nabla_{z}\cdot\left[L\left(\bar{u}_{t}\right)z\rho_{t}\left(z\right)+\varGamma\left(zz^{\intercal}-R_{t}\right)\rho_{t}\left(z\right)\right]+\frac{1}{2}\nabla_{z}\cdot\left[\nabla_{z}\cdot\left(\Sigma_{t}\Sigma_{t}^{\intercal}\rho_{t}\left(z\right)\right)\right],$
(2.15)

where $\mathcal{L}_{t}^{*}$ is the adjoint of the generator $\mathcal{L}_{t}$
that is also dependent on the law of $Z_{t}$ shown in the statistics of the
mean $\bar{u}_{t}$ and covariance $R_{t}$. We postpone the detailed analysis
for the convergence of the ensemble approximation model in Section 5 together
with the filter approximation.

In general, the probability density $\rho_{t}$ will demonstrate non-Gaussian
features due to the nonlinear stochastic coupling effects. On the other hand,
though the equations for the stochastic state $Z_{t}$ contain nonlinear
structures, its PDF $\rho_{t}$ is subject to only linear dynamics (2.15) with
$\mathcal{L}_{t}^{*}\left(\bar{u}_{t},R_{t}\right)$, which is dependent on the
mean and covariance. These desirable features inspires the construction of
effective filtering methods in the next section to include leading-order
statistical observations to improve forecast of highly non-Gaussian
statistics.

##  3 Filtering models using observations in mean and covariance

In this section, we exploit the ideas in filtering to propose improved
approximation of probability distributions of the model state containing
highly non-Gaussian statistics. The optimal filter is developed by combining
the stochastic forecast model describing unobserved microscopic states and
leading-order statistics introduced as macroscopic observations. We start with
a precise description of the optimal filter equations satisfying a conditional
Gaussian process, then a new statistical filtering model are proposed
approximating the optimal filtering solution with equivalent statistics.

###  3.1 Filtering probability distributions using statistical observations

We first formulate the filtering problem for predicting probability
distributions based on observations from the leading-order statistics. From
the stochastic-statistical equations (2.7), we can reformulate the general
multiscale system (1.1) for $u_{t}$ as a composition of the macroscopic state
from the first two moments $\bar{u}_{t},R_{t}$ and the microscopic stochastic
processes $Z_{t}$. In practice, measurements are often available in the
macroscopic statistical states (such as the locally averaged mean state and
variances as the deviation from the mean). Therefore, it is natural to
incorporate the statistical observation data to improve the estimation of the
unobserved microscopic processes, especially to recover the unobserved higher-
order statistics (such as the deviation from the normal distribution
indicating the occurrence of high impact extreme events).

Using similar idea of the Fokker-Planck filter introduced in [1], we assume
that the signal process of the filtering problem is from the model state PDF
$\rho_{t}\in\mathcal{P}\left(\mathbb{R}^{d}\right)$ belonging to the space of
continuous probability density functions on $\mathbb{R}^{d}$ with all finite
moments. The observation process is generated by the finite-dimensional
statistical moments denoted as $y_{t}\in\mathbb{R}^{p}$. This leads to the
general _infinite-dimensional filtering system with statistical observations_

$\displaystyle\mathrm{d}\rho_{t}=$
$\displaystyle\>\mathcal{L}_{t}^{*}\left(y_{t}\right)\rho_{t}\mathrm{d}t,\qquad\qquad\qquad\quad\rho_{t=0}\sim\mu_{0},$
(3.1a) $\displaystyle\mathrm{d}y_{t}=$
$\displaystyle\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Gamma_{t}\mathrm{d}B_{t},\quad
y_{t=0}=y_{0},$ (3.1b)

where $\mu_{0}$ is a probability measure of the
$\mathcal{P}\left(\mathbb{R}^{d}\right)$-valued random field $\rho$. Given
$y_{t}$, we will have $\rho_{t}\in\mathcal{P}\left(\mathbb{R}^{d}\right)$ for
all $t>0$. Above, $\mathcal{L}_{t}\left(y_{t}\right)$ is the infinitesimal
generator of the corresponding SDE for $Z_{t}$ with the explicit form given by
(2.15); the general observation process $y_{t}$ satisfies the dynamical
equation subject to a linear observation operator
$\mathcal{H}:\mathcal{P}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{p}$
acting on the PDF $\rho_{t}$, as well as
$h_{t}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$. Detailed equations for
$y_{t}=\left(\bar{u}_{t},R_{t}\right)$ will be given next in (3.3) based on
the first two moments. Notice here $\rho_{t}$ becomes a random field (with a
precise definition given in (3.5)) due to the randomness in $y_{t}$ as well as
its initial uncertainty.

#### Statistical observations from leading-order moments containing errors

Let $p_{t}$ be the (unknown) PDF of the state $u_{t}$ in (1.1), that is, the
deterministic solution of the FPE (1.2). Then, we can assume that observations
are drawn from the mean or covariance of the state $u_{t}$ as

$\bar{u}_{k,t}=\mathbb{E}_{p_{t}}\left[u_{t}\cdot\hat{v}_{k}\right],\quad
R_{kl,t}=\mathbb{E}_{p_{t}}\left[\left(\hat{v}_{k}\cdot
u_{t}^{\prime}\right)\left(u_{t}^{\prime}\cdot\hat{v}_{l}\right)\right],$
(3.2)

projected to the observed large-scale modes $\hat{v}_{k},k\leq d^{\prime}$ in
(2.1). We refer it as the _full observation case_ with $d^{\prime}=d$ and
_partial observation case_ with $d^{\prime}<d$. For simplicity, we may always
consider the full observation case $d^{\prime}=d$ (that is,
$y_{t}=\left(\bar{u}_{t},R_{t}\right)\in\mathbb{R}^{p}$ with $p=d+d^{2}$) in
this paper without confusion. The dynamical equations for $\bar{u}_{t},R_{t}$
can be introduced according to the closure model (2.7). According to
Proposition 2, the statistical equations for $\bar{u}_{t},R_{t}$ in (2.7)
provide statistical solutions consistent with the law $p_{t}$ of the state
$u_{t}$.

On the other hand, model errors are usually introduced due to the finite
truncation of the originally infinite-dimensional continuous system as well as
measurement errors. Therefore, additional correction terms, using independent
white noises, $B_{m,t}$ and $B_{v,t}$, are added to the equations (2.3) for
observed statistics accounting for errors from the imperfect model
approximations. The detailed equations for the observations (3.1b) can be
rewritten according to (2.3) as the following SDEs for
$\bar{u}_{t}\in\mathbb{R}^{d},R_{t}\in\mathbb{R}^{d^{2}}$

$\displaystyle\mathrm{d}\bar{u}_{t}=$
$\displaystyle\left[\mathcal{H}_{m}\rho_{t}+h_{m,t}\left(\bar{u}_{t}\right)\right]\mathrm{d}t\;+\Gamma_{m}\mathrm{d}B_{m,t},$
(3.3) $\displaystyle\mathrm{d}R_{t}=$
$\displaystyle\left[\mathcal{H}_{v}\rho_{t}+h_{v,t}\left(\bar{u}_{t},R_{t}\right)\right]\mathrm{d}t+\Gamma_{v}\mathrm{d}B_{v,t}.$

where $h_{m,t}\left(\bar{u}\right)=M\left(\bar{u}\right)+F_{t}$ and
$h_{v,t}\left(\bar{u},R\right)=L\left(\bar{u}\right)R+RL\left(\bar{u}\right)^{\intercal}+\Sigma_{t}\Sigma_{t}^{\intercal}$
are deterministic functions, while the linear _observation operators_ ,
$\mathcal{H}_{m},\mathcal{H}_{v}$, are defined by the high-order statistical
feedback functions (2.8)

$\displaystyle\mathcal{H}_{m,k}\rho=\int_{\mathbb{R}^{d}}H_{k}^{m}\left(z\right)\rho\left(z\right)\mathrm{d}z=$
$\displaystyle\sum_{p,q=1}^{d}\gamma_{kpq}\int_{\mathbb{R}^{d}}z_{p}z_{q}\rho\left(z\right)\mathrm{d}z,$
(3.4)
$\displaystyle\mathcal{H}_{v,kl}\rho=\int_{\mathbb{R}^{d}}H_{kl}^{v}\left(z\right)\rho\left(z\right)\mathrm{d}z=$
$\displaystyle\sum_{p,q=1}^{d}\left[\gamma_{kpq}\int_{\mathbb{R}^{d}}z_{p}z_{q}z_{l}\rho\left(z\right)\mathrm{d}z+\gamma_{lpq}\int_{\mathbb{R}^{d}}z_{p}z_{q}z_{k}\rho\left(z\right)\mathrm{d}z\right].$

We assume $\rho\in\mathcal{P}\left(\mathbb{R}^{d}\right)$ the probability
density with finite moments in all orders, thus we have
$\mathcal{H}_{m,k}\rho<\infty$ and $\mathcal{H}_{v,kl}\rho<\infty$ for all
$k,l$. Then (3.3) fits into the general observation equation (3.1b) by setting
$y_{t}=\left(\bar{u}_{t},R_{t}\right)^{\intercal}\in\mathbb{R}^{p}$ with
$p=d+d^{2}$ in a column vector, and letting
$\mathcal{H}=\left(\mathcal{H}_{m},\mathcal{H}_{v}\right)^{\intercal}$,
$h_{t}=\left(h_{m,t},h_{v,t}\right)^{\intercal}$, and
$\Gamma_{t}=\mathrm{diag}\left(\Gamma_{m},\Gamma_{v}\right)$. With this
explicit setup of the filtering problem, we will consider the optimal
filtering solution for the probability density $\rho_{t}$ of $Z_{t}$ based on
the statistical observation data
$Y_{t}=\left\\{\left(\bar{u}_{s},R_{s}\right),s\leq t\right\\}$.

###  3.2 The optimal filter with conditional Gaussian structure

Let $\left(\Omega,\mathcal{F},\mathbb{P}\right)$ be the complete probability
space, and denote $\mathcal{P}\left(\mathbb{R}^{d}\right)$ as the space of
probability density functions with bounded all moments. We first define the
$\mathcal{P}\left(\mathbb{R}^{d}\right)$-valued stochastic process $\rho_{t}$
as

$\rho_{t}:\mathbb{R}^{d}\times\Omega\rightarrow\mathbb{R}^{+},\;\left(z,\omega\right)\mapsto\rho_{t}\left(z;\omega\right),\;\mathrm{with}\;\rho_{t}\left(\cdot;\omega\right)\in\mathcal{P}\left(\mathbb{R}^{d}\right),$
(3.5)

which is thereafter referred to as a _random field_. In contrast to the
standard filtering problem concerning the nonlinear SDE of the random model
states $Z_{t}$, for derivation purpose of the exact optimal equations, we lift
the problem into filtering the random field $\rho_{t}$ based on the
observation information in $y_{s},s\leq t$ as in (3.1). A stochastic model on
$\mathbb{R}^{d}$ (3.16) will be then proposed for practical implementations
next in Section 2. Let $\mathcal{G}_{t}=\sigma\left\\{\omega:y_{s},s\leq
t\right\\}$ be the $\sigma$-algebra generated by the observations up to time
$t$. We define the space as the collection of $\mathcal{G}_{t}$-measurable
square-integrable random fields

$\mathcal{V}_{t}\coloneqq
L^{2}\left(\Omega,\mathcal{G}_{t},\mathbb{P};\mathcal{P}\left(\mathbb{R}^{d}\right)\cap
L^{2}\left(\mathbb{R}^{d}\right)\right),$ (3.6)

satisfying
$\int\left\|\nu\left(\cdot;\omega\right)\right\|_{L^{2}\left(\mathbb{R}^{d}\right)}^{2}\mathrm{d}\mathbb{P}\left(\omega\right)<\infty$
and $\nu\left(\cdot;\omega\right)\in\mathcal{P}\left(\mathbb{R}^{d}\right)\cap
L^{2}\left(\mathbb{R}^{d}\right)$ for $\nu\in\mathcal{V}_{t}$. In this
infinite-dimensional filtering problem, we aim to find the optimal
approximation of $\rho_{t}$ in the space $\mathcal{V}_{t}$. The _optimal
filtering solution_ $\hat{\rho}_{t}$ is then introduced as the least-square
estimate with the minimum variance as

$\hat{\rho}_{t}\coloneqq\underset{\nu\in\mathcal{V}_{t}}{\arg\min}\>\mathbb{E}\left[\left\|\rho_{t}-\nu\right\|_{L^{2}\left(\mathbb{R}^{d}\right)}^{2}\right]=\mathsf{P}_{\mathcal{V}_{t}}\left[\rho_{t}\right],$
(3.7)

where the optimal solution $\hat{\rho}_{t}$ can be viewed as the unbiased
projection of $\rho_{t}$ onto the space $\mathcal{V}_{t}$. (3.7) indicates
that $\hat{\rho}_{t}$ gives the PDF estimation closest to the true
distribution $\rho_{t}$ in the mean square sense in agreement with the
observations in $\mathcal{G}_{t}$.

Accordingly, we define the _optimal filter distribution_
$\mu_{t}:\mathcal{P}\left(\mathbb{R}^{d}\right)\times\Omega\rightarrow\left[0,1\right]$
as the regular conditional measure of the stochastic process $\rho_{t}$ given
$\mathcal{G}_{t}$. That is, for any Borel set
$A\in\mathcal{B}\left(\mathcal{P}\left(\mathbb{R}^{d}\right)\right)$,
$\mu_{t}$ gives the conditional probability of $\rho_{t}$ given
$\mathcal{G}_{t}$

$\mu_{t}\left(A;\omega\right)\coloneqq\mathbb{P}\left(\rho_{t}\in
A\mid\mathcal{G}_{t}\right)\left(\omega\right),\quad\mathrm{a.s.}\>\omega\in\Omega.$
(3.8)

Notice that $\mu_{t}\left(A;\cdot\right)\in\mathcal{G}_{t}$ is still a
stochastic process. For any functional $F\in
C\left(\mathcal{P}\left(\mathbb{R}^{d}\right)\right)$ and $t>0$, we can
introduce the conditional expectation w.r.t. the measure $\mu_{t}$ given
$\mathcal{G}_{t}$ as

$\mathbb{E}\left[F\left(\rho_{t}\right)\mid\mathcal{G}_{t}\right]\coloneqq\int_{\mathcal{P}\left(\mathbb{R}^{d}\right)}F\left(\rho\right)\mu_{t}\left(\mathrm{d}\rho\right).$

Therefore, the optimal filter solution (3.7) is a random field

$\hat{\rho}_{t}=\mathbb{E}\left[\rho_{t}\mid\mathcal{G}_{t}\right]$ (3.9)

given by the conditional expectation of $\rho_{t}$ w.r.t. $\mu_{t}$. For any
linear operator
$\mathcal{M}:L^{2}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{p}$, we
have
$\mathcal{M}\hat{\rho}_{t}=\mathbb{E}\left[\mathcal{M}\rho_{t}\mid\mathcal{G}_{t}\right]$.
Furthermore, second moment of $\mathcal{M}\rho_{t}$ is given by

$\mathbb{E}\left[\left(\mathcal{M}\rho_{t}-\mathcal{M}\hat{\rho}_{t}\right)\left(\mathcal{M}\rho_{t}-\mathcal{M}\hat{\rho}_{t}\right)^{\intercal}\mid\mathcal{G}_{t}\right]=\mathcal{M}\hat{\mathcal{C}}_{t}\mathcal{M}^{*},$

where
$\hat{\mathcal{C}}_{t}\left(\omega\right):L^{2}\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)\rightarrow
L^{2}\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$ with
$\hat{\mathcal{C}}_{t}^{*}=\hat{\mathcal{C}}_{t}$ is the self-adjoint
covariance operator for any $f\in
L^{2}\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$

$\hat{\mathcal{C}}_{t}f=\mathbb{E}\left[\left(\rho-\hat{\rho}_{t}\right)\int_{\mathbb{R}^{d}}\left(\rho-\hat{\rho}_{t}\right)\left(z\right)f\left(z\right)\mathrm{d}z\mid\mathcal{G}_{t}\right].$
(3.10)

Notice again $\hat{\mathcal{C}}_{t}f\left(z;\cdot\right)$ is also a random
field conditional on $\mathcal{G}_{t}$. For clarification of notations, we
will call $\hat{\rho}_{t}$ and $\hat{\mathcal{C}}_{t}$ the optimal filter
solution for the mean and covariance, and $\mu_{t}$ the optimal filter
distribution in the rest part of the paper.

In particular, we can characterize the optimal filter solution
$\hat{\rho}_{t}$ as the best estimate in each order of moments. The following
result describes the accuracy of the filter approximations under any finite-
dimensional projections.

######  Proposition 3.

Let $\rho_{t}$ be the random field from the system (3.1) and
$\hat{\rho}_{t}=\mathbb{E}\left[\rho_{t}\mid\mathcal{G}_{t}\right]$ the
optimal filter solution given the observations in $\mathcal{G}_{t}$. For any
linear operator
$\mathcal{M}:\mathcal{P}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{p}$
defined by $\mathcal{M}\rho=\int M\left(z\right)\rho\left(\mathrm{d}z\right)$
with $M\in C\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$,
$\mathcal{M}\hat{\rho}_{t}=\mathbb{E}\left[\mathcal{M}\rho_{t}\mid\mathcal{G}_{t}\right]$
gives the best unbiased estimate of $\mathcal{M}\rho_{t}$ in the sense of
minimum mean square error, that is,

$\mathbb{E}\left[\left|\mathcal{M}\rho_{t}-\mathcal{M}\hat{\rho}_{t}\right|^{2}\right]=\min_{\nu\in\mathcal{V}_{t}}\mathbb{E}\left[\left|\mathcal{M}\rho_{t}-\mathcal{M}\nu\right|^{2}\right],\quad\mathrm{with}\quad\mathbb{E}\left[\mathcal{M}\hat{\rho}_{t}\right]=\mathbb{E}\left[\mathcal{M}\rho_{t}\right].$
(3.11)

###### Proof.

For any $\mathcal{G}_{t}$-measurable square-integrable stochastic process
$\nu$, we have from direct computation

$\displaystyle\mathbb{E}\left[\left|\mathcal{M}\rho_{t}-\mathcal{M}\nu\right|^{2}\right]-\mathbb{E}\left[\left|\mathcal{M}\rho_{t}-\mathcal{M}\hat{\rho}_{t}\right|^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\mathcal{M}\hat{\rho}_{t}-\mathcal{M}\nu\right)\cdot\left(2\mathcal{M}\rho_{t}-\mathcal{M}\nu-\mathcal{M}\hat{\rho}_{t}\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\mathbb{E}\left[\left(\mathcal{M}\hat{\rho}_{t}-\mathcal{M}\nu\right)\cdot\left(2\mathcal{M}\rho_{t}-\mathcal{M}\nu-\mathcal{M}\hat{\rho}_{t}\right)\mid\mathcal{G}_{t}\right]\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\left(\mathcal{M}\hat{\rho}_{t}-\mathcal{M}\nu\right)\cdot\left[\mathbb{E}\left[2\mathcal{M}\rho_{t}\mid\mathcal{G}_{t}\right]-\left(\mathcal{M}\nu+\mathcal{M}\hat{\rho}_{t}\right)\right]\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left|\mathcal{M}\hat{\rho}_{t}-\mathcal{M}\nu\right|^{2}\right]\geq
0.$

Above, the third equality uses the fact
$\mathcal{M}\hat{\rho}_{t}-\mathcal{M}\nu=\mathbb{E}\left[\mathcal{M}\rho_{t}\mid\mathcal{G}_{t}\right]-\mathcal{M}\nu$
is $\mathcal{G}_{t}$-measurable. Thus, we get $\mathcal{M}\hat{\rho}_{t}$
minimizes the mean square error. The consistency under the expectation in
$\mathcal{M}\hat{\rho}_{t}$ and $\mathcal{M}\rho_{t}$ can be directly implied
by definition. ∎

By taking the operator $\mathcal{M}$ as the expectation on
$M\left(Z_{t}\right)=Z_{t}^{m}$ with any integer $m$, $\mathcal{M}\rho_{t}$
and $\mathcal{M}\hat{\rho}_{t}$ give the $m$-th order moments of $Z_{t}$ under
the random field $\rho_{t}$ in (3.5) and optimal filter approximation
$\hat{\rho}_{t}$. A direct implication from (3.11) shows that we have the
unbiased statistics in all high-order moments
$\mathbb{E}\left[\mathcal{M}\hat{\rho}_{t}\right]=\mathbb{E}\left[\mathcal{M}\rho_{t}\right]$
with the minimum error
$\mathbb{E}\left[\left|\mathcal{M}\rho_{t}-\mathcal{M}\hat{\rho}_{t}\right|^{2}\right]$
from the optimal filter solution.

Importantly, the model equations (3.1) satisfy the desirable conditional
Gaussian process [36], that is, given the observations of
$Y_{t}=\left\\{y_{s}=\left(\bar{u}_{s},R_{s}\right),s\leq t\right\\}$ and
Gaussian initial state $\rho_{0}$, the random field $\rho_{t}$ follows a
Gaussian distribution at each time $t$. Let
$\rho_{t}=\rho_{t}\left(\cdot;\omega\right)$ be the (unknown) signal state
satisfying linear dynamics (3.1a), and $y_{t}=y_{t}\left(\omega\right)$ the
observed statistical process subject to linear observation operators (3.1b).
The optimal filter distribution $\mu_{t}$ (3.8) conditional on $Y_{t}$ then
becomes an infinite-dimensional Gaussian distribution,
$\mu_{t}\left(\cdot,\omega\right)=\mathcal{N}\left(\hat{\rho}_{t},\hat{\mathcal{C}}_{t}\right)\left(\omega\right)$,
where the mean $\hat{\rho}_{t}\left(\cdot;\omega\right)$ and covariance
$\hat{\mathcal{C}}_{t}\left(\omega\right)$ give the solution to (3.9) and
(3.10) respectively. Therefore, the equations for the mean and covariance are
given by the generalized version of Kalman-Bucy (KB) filter for the infinite-
dimensional conditional Gaussian processes

$\displaystyle\mathrm{d}\hat{\rho}_{t}=$
$\displaystyle\>\mathcal{L}_{t}^{*}\left(\bar{u}_{t},R_{t}\right)\hat{\rho}_{t}\mathrm{d}t+\mathcal{\hat{C}}_{t}\mathcal{H}_{m}^{*}\Gamma_{m}^{-2}\left\\{\mathrm{d}\bar{u}_{t}-\left[\mathcal{H}_{m}\hat{\rho}_{t}+h_{m,t}\left(\bar{u}_{t}\right)\right]\mathrm{d}t\right\\}$
(3.12)
$\displaystyle+\hat{\mathcal{C}}_{t}\mathcal{H}_{v}^{*}\Gamma_{v}^{-2}\left\\{\mathrm{d}R_{t}-\left[\mathcal{H}_{v}\hat{\rho}_{t}+h_{v,t}\left(\bar{u},R\right)\right]\mathrm{d}t\right\\},$
$\displaystyle\mathrm{d}\hat{\mathcal{C}}_{t}=$
$\displaystyle\left[\mathcal{L}_{t}^{*}\left(\bar{u}_{t},R_{t}\right)\mathcal{\hat{C}}_{t}+\mathcal{\hat{C}}_{t}\mathcal{L}_{t}\left(\bar{u}_{t},R_{t}\right)\right]\mathrm{d}t-\mathcal{\hat{C}}_{t}\left(\mathcal{H}_{m}^{*}\Gamma_{m}^{-2}\mathcal{H}_{m}+\mathcal{H}_{v}^{*}\Gamma_{v}^{-2}\mathcal{H}_{v}\right)\mathcal{\hat{C}}_{t}\mathrm{d}t.$

The equations of the conditional Gaussian processes and the uniqueness of the
solutions are developed in Chapter 12 of [36] for finite dimensional systems.
The results are then generalized to infinite-dimensional Hilbert space [10,
20] (see a summary of the main results in Appendix A). The system (3.12) gives
a closed set of coupled SPDEs (due to the randomness in $\bar{u}_{t},R_{t}$)
for the optimal filter solution $\hat{\rho}_{t}$ enabling more detailed
analysis and development of practical methods for computing the optimal
solution.

######  Remark.

The filtering problem using statistical observations is originally introduced
as the ensemble Fokker-Planck filter in [1]. We propose the filtering
equations (3.12) as a further generalization where nonlinearly coupled
conditional processes are involved. In particular, this new filtering model
directly fits into the coupled stochastic-statistical modeling framework
(2.7).

###  3.3 A surrogate filtering model for approximating the optimal filter
solution

The resulting optimal filtering problem from (3.1) requires to solve the
infinite-dimensional system (3.12) concerning the function
$\hat{\rho}_{t}\in\mathcal{P}\left(\mathbb{R}^{d}\right)$ and operator
$\hat{\mathcal{C}}_{t}$ on $\mathcal{P}\left(\mathbb{R}^{d}\right)$. It
usually becomes intractable in finding such infinite-dimensional solutions
from direct methods. In developing practical strategies to realize the optimal
filter solution, it is more useful to find a surrogate model for the
stochastic process $\tilde{Z}_{t}$, based on which effective ensemble-based
approaches can be built. Therefore, we aim to construct an approximating
filtering model from designing a new dynamical equation for $\tilde{Z}_{t}$,
whose PDF $\tilde{\rho}_{t}$ can effectively represent that of the optimal
filter solution $\hat{\rho}_{t}$.

####  3.3.1 Filtering updating cycle in a split two-step procedure

For a clear characterization of the filtering process, we follow the general
procedure in [7] to first describe the filtering process by concatenated
iterations of transporting maps on the corresponding probability distribution.
We propose a new stochastic process $\tilde{Z}_{t}$, whose law
$\tilde{\rho}_{t}\in\mathcal{V}_{t}$ is still a
$\mathcal{P}\left(\mathbb{R}^{d}\right)$-valued random field in
$\mathcal{G}_{t}$ dependent on the same statistical observation $y_{t}$ as in
the optimal filter (3.12). Thus, the filtering updating cycle during the time
interval $\left[t,t+\tau\right]$ can be characterized by the transport of the
probability density $\tilde{\rho}_{t}$ of $\tilde{Z}_{t}$.

The discrete time update of the approximation filter PDF $\tilde{\rho}_{t}$ is
carried out in a split two-step procedure. First, the _forecast step_ can be
viewed as the push-forward operator acting on the probability density function
at time instant $t$ with time step $\tau>0$

$\tilde{\rho}_{t}\rightarrow\tilde{\rho}_{t+\tau}^{-}\coloneqq\mathcal{F}_{t}^{\tau}\tilde{\rho}_{t}=e^{\int_{t}^{t+\tau}\mathcal{L}_{s}^{*}\left(y_{s}\right)\mathrm{d}s}\tilde{\rho}_{t},$
(3.13)

where $\mathcal{F}_{t}^{\tau}$ represents the forecast updating operator with
forward time step $\tau$, and $\mathcal{L}_{t}\left(y_{t}\right)$ is the same
generator as in (2.15). Second, the _analysis step_ updates the prior
distribution $\tilde{\rho}^{-}$ to the posterior distribution
$\tilde{\rho}^{+}$ by incorporating the observation data up to
$Y_{t+\tau}=\left\\{y_{s},s\leq t+\tau\right\\}$, that is

$\tilde{\rho}_{t+\tau}^{-}\rightarrow\tilde{\rho}_{t+\tau}^{+}\coloneqq\mathcal{A}_{t}^{\tau}\left(\tilde{\rho}_{t+\tau}^{-};Y_{t+\tau}\right),$
(3.14)

where $\mathcal{A}_{t}^{\tau}$ represents the analysis updating operator.
Therefore, the full filtering cycle from $t$ to $t+\tau$ can be summarized as
the composition of the forecast and analysis maps

$\tilde{\rho}_{t+\tau}=\mathcal{A}_{t}^{\tau}\left(\mathcal{F}_{t}^{\tau}\tilde{\rho}_{t};Y_{t+\tau}\right).$
(3.15)

Notice that $\mathcal{F}_{t}^{\tau}$ is a linear operator on
$\tilde{\rho}_{t}$, while $\mathcal{A}_{t}^{\tau}$ could contain nonlinear
actions due to the normalization of the probability distribution. Continuous
equation for $\partial_{t}\tilde{\rho}_{t}=\lim_{\tau\rightarrow
0}\frac{1}{\tau}\left(\tilde{\rho}_{t+\tau}-\tilde{\rho}_{t}\right)$ is then
achieved by letting the discrete time step $\tau\rightarrow 0$. Next, we first
propose the general structure of the new filtering model for
$\tilde{Z}_{t}\sim\tilde{\rho}_{t}$ as a combination of the above two-step
procedure, then detailed analysis can be done according to the design of the
forecast and analysis step operators $\mathcal{F}_{t}^{\tau}$ and
$\mathcal{A}_{t}^{\tau}$ accordingly.

####  3.3.2 Construction of equivalent statistical approximating filter

For simplicity of notations, we still use the general statistical observation
processes (3.3) for $y_{t}=\left(\bar{u}_{t},R_{t}\right)$ taking the compact
form

$\mathrm{d}y_{t}=\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Gamma_{t}\mathrm{d}B_{t},$

where the general observation operator (3.4), $\mathcal{H}\rho_{t}=\int
H\left(z\right)\rho_{t}\left(z\right)\mathrm{d}z$, is defined with the general
observation function $H\in C\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$ acting
on the density function $\rho_{t}$. Following the general construction in [8,
45], we seek the approximating filter model adopting the following McKean-
Vlasov representation with undetermined functionals $a_{t},K_{t}$

$\displaystyle\mathrm{d}\tilde{Z}_{t}=$
$\displaystyle\>L\left(\bar{u}_{t}\right)\tilde{Z}_{t}\mathrm{d}t+\varGamma\left(\tilde{Z}_{t}\tilde{Z}_{t}^{\intercal}-R_{t}\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}\tilde{W}_{t}$
(3.16) $\displaystyle+$ $\displaystyle
a_{t}\left(\tilde{Z}_{t};\tilde{\rho}_{t}\right)\mathrm{d}t+K_{t}\left(\tilde{Z}_{t};\tilde{\rho}_{t}\right)\left\\{\mathrm{d}y_{t}-\left[H\left(\tilde{Z}_{t}\right)+h_{t}\left(y_{t}\right)\right]\mathrm{d}t-\Gamma_{t}\mathrm{d}\tilde{B}_{t}\right\\}.$

We use the name ‘statistical filtering’ to refer the above new model
emphasizing our main goal of filtering statistical moments different from the
common filtering case. The first row of the above equation models the forecast
step of the filtering process, while the second row is the analysis step. The
forecast step accepts the same dynamical model of (3.1a) dependent on the mean
and covariance $\left(\bar{u}_{t},R_{t}\right)$. On the other hand, the
analysis step serves as an additional control correction over statistical
observations $y_{t}$. New functionals known as the drift
$a_{t}:\mathbb{R}^{d}\times\mathcal{P}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{d}$
and the control gain operator
$K_{t}:\mathbb{R}^{d}\times\mathcal{P}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{d\times
p}$ are introduced, resulting in an approximating filtering model about the
process $\tilde{Z}_{t}$. Most importantly, as will be shown next in Theorem 7
under proper condition, $a_{t}$ and $K_{t}$ are only implicitly dependent on
$\tilde{\rho}_{t}$ through the leading moments, without the need to compute
the (potentially highly non-Gaussian) density function $\tilde{\rho}_{t}$
explicitly.

In a more clear identification of the filtering updates involving several
levels of approximations, we take the split-step strategy to analyze the
coupled forecast step and analysis step of the filtering equation (3.16)
separately. In particular, the forecast step in the first row of (3.16) is
given by the exactly same form as the stochastic-statistical closure equations
(2.7) developed in Section 2. Thus in practice, the updating step with the
forecast operator can be implemented adopting the efficient uncertainty
prediction methods such as [38, 47]. Then, the remaining task is to propose
proper analysis step update in the second line of (3.16) concerning consistent
statistics with the optimal solution $\hat{\rho}_{t}$ in (3.12).

###  3.4 Statistical consistency in analysis step update of the approximating
filter

Now, we focus on updating posterior PDF $\tilde{\rho}_{t}$ in (3.14) of the
proposed approximating filter (3.16) based on the statistical observation
$y_{t}$ satisfying (3.3). Concentrating on the analysis step, the resulting
optimal filter equations (3.12) for the mean and covariance
$\left(\hat{\rho}_{t},\hat{\mathcal{C}}_{t}\right)$ become

$\displaystyle\mathrm{d}\hat{\rho}_{t}=$
$\displaystyle\>\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\Gamma_{t}^{-2}\left\\{\mathrm{d}y_{t}-\left[\mathcal{H}\hat{\rho}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(3.17) $\displaystyle\mathrm{d}\mathcal{\hat{C}}_{t}=$
$\displaystyle-\hat{\mathcal{C}}_{t}\mathcal{H}^{*}\Gamma_{t}^{-2}\mathcal{H}\mathcal{\hat{C}}_{t}\mathrm{d}t.$

Correspondingly, the approximating statistical filtering model for
$\tilde{Z}_{t}$ satisfies the second line of the SDE (3.16) as

$\displaystyle\mathrm{d}\tilde{Z}_{t}=$
$\displaystyle\>a_{t}\left(\tilde{Z}_{t}\right)\mathrm{d}t+K_{t}\left(\tilde{Z}_{t}\right)\left\\{\mathrm{d}y_{t}-\left[H\left(\tilde{Z}_{t}\right)+h_{t}\left(y_{t}\right)\right]\mathrm{d}t-\Gamma_{t}\mathrm{d}\tilde{B}_{t}\right\\}.$
(3.18)

Following the similar idea in the McKean-Vlasov representation of the
filtering equation [45, 57], we expect the PDF $\tilde{\rho}_{t}$ of
$\tilde{Z}_{t}$ to satisfy the following Kushner-Stratonovich-type equation
(with requirements on $a_{t},K_{t}$ given next in (3.21))

$\frac{\partial\tilde{\rho}_{t}}{\partial
t}=\left[H\left(z\right)-\mathcal{H}\tilde{\rho}_{t}\right]^{\intercal}\Gamma_{t}^{-2}\left[\frac{\mathrm{d}y_{t}}{\mathrm{d}t}-\mathcal{H}\tilde{\rho}_{t}-h_{t}\left(y_{t}\right)\right]\tilde{\rho}_{t}.$
(3.19)

Again, the goal here is to approximate the optimal filter mean
$\hat{\rho}_{t}$ in (3.17) by $\tilde{\rho}_{t}$ generated by the surrogate
SDE model (3.18) in the sense of consistent statistics.

Unfortunately, the approximation (3.19) and the optimal filtering equation
(3.17) will in general have different continuous solutions for
$\tilde{\rho}_{t}$ and $\hat{\rho}_{t}$ due to their distinctive dynamics. In
order to compare the statistical moments of the two distributions, we apply
the linear operator $\mathcal{H}$ to the optimal equations (3.17) as a finite-
dimensional projection on leading moments. The resulting equations optimal
mean and covariance equations become finite dimensional as

$\displaystyle\mathrm{d}\left(\mathcal{H}\hat{\rho}_{t}\right)=$
$\displaystyle\left(\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\right)\Gamma_{t}^{-2}\left\\{\mathrm{d}y_{t}-\left[\mathcal{H}\hat{\rho}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(3.20)
$\displaystyle\mathrm{d}\left(\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\right)=$
$\displaystyle-\left(\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\right)\Gamma_{t}^{-2}\left(\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\right)\mathrm{d}t.$

Above, remind that the observation operator
$\mathcal{H}:L^{2}\left(\mathbb{R}^{d}\right)\rightarrow\mathbb{R}^{p}$ and
its adjoint $\mathcal{H}^{*}:\mathbb{R}^{p}\rightarrow
L^{2}\left(\mathbb{R}^{d}\right)$ are defined based on the observation
function $H\in C_{b}^{2}\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$ as

$\mathcal{H}\rho=\int
H\left(z\right)\rho\left(z\right)\mathrm{d}z,\quad\left[\mathcal{H}^{*}u\right]\left(z\right)=u\cdot
H\left(z\right),$

and the covariance operator
$\hat{\mathcal{C}}_{t}:L^{2}\left(\mathbb{R}^{d}\right)\rightarrow
L^{2}\left(\mathbb{R}^{d}\right)$ is defined in (3.10). Therefore, (3.20)
gives the equations for the finite-dimensional quantities
$\mathcal{H}\hat{\rho}_{t}\in\mathbb{R}^{p}$ and
$\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\in\mathbb{R}^{p\times p}$ as
the first two moments of $H$ w.r.t. $\hat{\rho}_{t}$. The idea is to design
the analysis step operator $\mathcal{A}_{t}^{\tau}$ for the approximating
filter process $\tilde{Z}_{t}$ in (3.18), so that consistency in the first and
second-order moments can be achieved.

Denote the expectation,
$\tilde{\mathbb{E}}\left[\cdot\right]\coloneqq\mathbb{E}\left[\cdot\mid\mathcal{G}_{t}\right]=\mathbb{E}_{\tilde{\rho}_{t}}\left[\cdot\right]$,
w.r.t. the conditional density $\tilde{\rho}_{t}$ in (3.19) given the same
observation process in $\mathcal{G}_{t}$. We define
$\bar{H}_{t}=\tilde{\mathbb{E}}H\left(\tilde{Z}_{t}\right)=\int
H\left(z\right)\tilde{\rho}_{t}\left(z\right)\mathrm{d}z$ and
$C_{t}^{H}=\tilde{\mathbb{E}}\left[H_{t}\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right]\left[H_{t}\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right]^{\intercal}$
as the first and second-order moments of $H$ w.r.t. $\tilde{\rho}_{t}$ .
Assume that the drift $a_{t}$ and gain $K_{t}$ in the SDE approximation (3.18)
satisfy the following identities

$a_{t}=\nabla\cdot\left(K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right)-K_{t}\Gamma_{t}^{2}\nabla\cdot
K_{t}^{\intercal},\quad-\nabla\cdot\left(K_{t}^{\intercal}\tilde{\rho}_{t}\right)=\tilde{\rho}_{t}\Gamma_{t}^{-2}\left(H\left(z\right)-\tilde{\mathbb{E}}H\right),$
(3.21)

where the divergence on a matrix is defined columnwise as $\left(\nabla\cdot
A\right)_{i}=\sum_{j}\partial_{z_{j}}A_{ij}.$We first have the following
result concerning the evolution equations of $\bar{H}_{t}$ and $C_{t}^{H}$
given the realization $Y_{t}=\left\\{y_{s},s\leq t\right\\}$.

######  Lemma 4.

Given that $\Gamma_{t}\succ 0$ in (3.1) and the identities (3.21) are
satisfied, the evolution equations for the mean and covariance of the
observation function $H\left(\tilde{Z}_{t}\right)$ associated with the SDE
(3.18) are given by

$\displaystyle\mathrm{d}\bar{H}_{t}=$
$\displaystyle\>C_{t}^{H}\Gamma_{t}^{-2}\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(3.22) $\displaystyle\mathrm{d}C_{t}^{H}=$
$\displaystyle\>Q_{t}^{H}\Gamma_{t}^{-2}\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\}-C_{t}^{H}\Gamma_{t}^{-2}C_{t}^{H}\mathrm{d}t,$

where $Q_{t}^{H}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p\times p}$ is defined
as

$Q_{t}^{H}=\tilde{\mathbb{E}}\left[\left(H_{t}^{\prime}H_{t}^{\prime\intercal}\right)\otimes
H_{t}^{\prime\intercal}\right],$

containing third moments of
$H_{t}^{\prime}=H\left(\tilde{Z}_{t}\right)-\bar{H}_{t}$.

We put the detailed derivation of (3.22) in Appendix B. Notice that (3.22)
goes back to the Kalman-Bucy filter if we set linear observation
$H\left(\tilde{Z}_{t}\right)=\tilde{Z}_{t}$ satisfying a normal distribution
as in [1, 9]. Although here we are considering the more general nonlinear
dynamics and observation functions from (3.4).

Comparing (3.22) and (3.20) implies that the same statistical solution can be
reached in $\left(\bar{H}_{t},C_{t}^{H}\right)$ and
$\left(\mathcal{H}\hat{\rho}_{t},\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}\right)$
if we have $Q_{t}^{H}=0$. In order to achieve this, we further introduce the
projection operator on the space of probability distributions using the
Kullback-Leibler (KL) divergence [30] as an unbiased metric.

######  Definition 5.

Define the operator $\mathcal{S}_{H}$ making symmetric projection on the
probability density $\rho\in\mathcal{P}\left(\mathbb{R}^{d}\right)$

$\mathcal{S}_{H}\rho=\underset{\nu\in\mathcal{V}_{H}}{\arg\min}\>d_{\mathrm{KL}}\left(\nu\parallel\rho\right),$
(3.23)

where $d_{\mathrm{KL}}$ is the KL divergence between two probability measures.
The minimization is among the probability measures in the following set

$\mathcal{V}_{H}\left(\bar{H},C^{H}\right)=\left\\{\nu\in\mathcal{P}\left(\mathbb{R}^{d}\right):\mathbb{E}_{\nu}H=\bar{H},\;\mathbb{E}_{\nu}\left[H^{\prime}H^{\prime\intercal}\right]=C^{H},\;\mathrm{and}\;\mathbb{E}_{\nu}\left[H_{l}^{\prime}H_{m}^{\prime}H_{n}^{\prime}\right]=0\right\\},$

for all $l,m,n\leq d$ and
$H^{\prime}\left(Z\right)=H\left(Z\right)-\mathbb{E}_{\nu}\left[H\left(Z\right)\right]$.

In Definition 5, $\mathcal{S}_{H}$ acts as a symmetric approximation of
probability measures with vanishing third-order moments of the observation
function $H$, while maintains consistent first two leading moments of $H$. It
is clear that given $\bar{H},C^{H}$, the set $\mathcal{V}_{H}$ is closed with
respect to weak convergence of measures. From Proposition 2.1 of [46] and
[16], we have for any $\rho$ and weakly convergent sequence
$\left\\{\nu_{n}\right\\}$ to $\nu_{*}$

$\underset{n\rightarrow\infty}{\lim\inf}\>d_{\mathrm{KL}}\left(\nu_{n}\parallel\rho\right)\geq
d_{\mathrm{KL}}\left(\nu_{*}\parallel\rho\right).$

It follows immediately that there exists $\nu_{*}\in\mathcal{V}_{H}$ that
reaches the minimum. Therefore, we have the following lemma guaranteeing the
existence of the minimizer in the proposed projection (3.23).

######  Lemma 6.

Assume that there is one $\nu\in\mathcal{V}_{H}$ such that the KL-divergence
$d_{\mathrm{KL}}\left(\nu\parallel\rho\right)<\infty$ for given
$\rho\in\mathcal{P}\left(\mathbb{R}^{d}\right)$. Then a minimizer exists in
(3.23).

Even though Lemma 6 does not guarantee the uniqueness of the minimizer, we can
always find one minimizer given the same mean and covariance
$\bar{H}=\tilde{\mathbb{E}}\left[H\left(Z\right)\right],C^{H}=\tilde{\mathbb{E}}\left[H^{\prime}\left(Z\right)H^{\prime}\left(Z\right)^{\intercal}\right]$.
This provides the desirable target density function satisfying the required
symmetric statistics about $H\left(Z_{t}\right)$. Denote the push-forward
operator for the new SDE (3.18) with the structure functions (3.21) as
$\tilde{\rho}_{t+s}=\mathcal{Q}_{t}^{s}\left(\tilde{\rho}_{t};Y_{t+s}\right)$
for any $s\geq 0$. Under the above construction, we can finally propose the
forward operator in analysis step (3.14) as

$\mathcal{A}_{t}^{s}\left(\tilde{\rho}_{t}\right)\coloneqq\mathcal{Q}_{t}^{s}\left(\mathcal{S}_{H}\tilde{\rho}_{t};Y_{t+s}\right),$
(3.24)

where the projection $\mathcal{S}_{H}$ in (3.23) is a linear operator acting
on the random fields in space $\mathcal{V}_{t}$ in (3.6). In a similar fashion
as in the proof of Lemma 4, by applying the addition projection
$\mathcal{S}_{H}\tilde{\rho}_{t}$ in the expectations on the SDEs for
$H\left(\tilde{Z}_{t}\right)$ and
$H_{t}^{\prime}\left(\tilde{Z}_{t}\right)H_{t}^{\prime}\left(\tilde{Z}_{t}\right)^{\intercal}$
(see also (B10) and (B11) in Appendix B), $Q_{t}^{H}=0$ is automatically
guaranteed in (3.22) w.r.t. the new projected density
$\mathcal{S}_{H}\tilde{\rho}_{t}$. In addition, the first two moments
$\bar{H}_{t}$ and $C_{t}^{H}$ (3.22) will stay the same w.r.t.
$\mathcal{S}_{H}\tilde{\rho}_{t}$. The continuous filtering process by letting
$s\rightarrow 0$ will satisfy the equations (3.22) with $Q_{t}^{H}\equiv 0$.
Therefore, under the same initial condition and the uniqueness of the
solution, the same solution will be reached in both (3.17) and (3.22). This
leads to the main result of this section concerning the analysis step update
in the approximating filter solution.

######  Theorem 7.

Consider the analysis step update (3.24) of the statistical filtering model
(3.18). Assume that $a_{t},K_{t}$ in the statistical filtering SDE are
designed to satisfy (3.21) and the probability set $\mathcal{V}_{H}$ according
to $H$ defined in (3.23) is not empty. Under the same statistical observations
$y_{t},t\in\left[0,T\right]$ and the same initial conditions, the following
relations hold for $t\in\left[0,T\right]$

$\mathcal{H}\hat{\rho}_{t}=\tilde{\mathbb{E}}\left[H\left(\tilde{Z}_{t}\right)\right],\quad\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{H}^{*}=\tilde{\mathbb{E}}\left[H^{\prime}\left(\tilde{Z}_{t}\right)H^{\prime}\left(\tilde{Z}_{t}\right)^{\intercal}\right],$
(3.25)

where $\left(\hat{\rho}_{t},\hat{\mathcal{C}}_{t}\right)$ is the solution of
(3.17), and $\tilde{\rho}_{t}$ is given by the solution of (3.19) with
$H^{\prime}=H-\mathcal{H}\hat{\rho}_{t}$.

Theorem 7 validates the use of the statistical filtering model density
$\tilde{\rho}_{t}$ by solving (3.16) to approximate the optimal filter
$\hat{\rho}_{t}$ from (3.12). Though restricted only on the first two moments
of the observation function $H$, the resulting consistent statistics during
analysis step play a key role in accurate statistical forecast. Notice that
based on the statistical model in (2.7), accurate prediction of the important
leading statistics, $\bar{u}_{t},R_{t}$, is determined by key higher-order
feedbacks in the related functional $\mathcal{H}\rho_{t}$ (more specifically,
the terms $\mathcal{H}_{m}\rho_{t}$ and $\mathcal{H}_{v}\rho_{t}$). According
to Proposition 3, the optimal $\mathcal{H}\hat{\rho}_{t}$ gives the least mean
square estimate of the random variable $\mathcal{H}\rho_{t}$ given the
statistical observations. Thus, consistent approximating filter
$\tilde{\mathbb{E}}H\left(\tilde{Z}_{t}\right)$ for
$\mathcal{H}\hat{\rho}_{t}$ as well as its error estimate guarantees accurate
recovery of key model statistics. For example, applying the explicit forms of
the observation function (3.4), the quadratic observation operator
$\mathcal{H}_{m}$ in the mean equation gives

$\tilde{\mathbb{E}}H^{m}\left(\tilde{Z}_{t}\right)=\mathcal{H}_{m}\hat{\rho}_{t}\;\Leftrightarrow\;\sum_{p,q}\gamma_{kpq}\tilde{\mathbb{E}}\left[\tilde{Z}_{p,t}\tilde{Z}_{q,t}\right]=\sum_{p,q}\gamma_{kpq}\int
z_{p}z_{q}\hat{\rho}_{t}\left(z\right)\mathrm{d}z,$

which implies consistent statistical feedbacks from the statistical filtering
model $\tilde{\rho}_{t}$ and the optimal filter solution $\hat{\rho}_{t}$.
This demonstrates that the new approximating filter maintains the accuracy in
the statistical mean prediction $\bar{u}_{t}$. In addition, the covariance
operator characterizes the essential uncertainty in the optimal filter
estimate $\hat{\rho}_{t}$ as a random field, that is,

$\mathcal{H}_{m}\mathcal{\hat{C}}_{t}\mathcal{H}_{m}^{*}=\mathbb{E}_{\mu_{t}}\left[\left(\mathcal{H}_{m}\rho_{t}-\mathcal{H}_{m}\hat{\rho}_{t}\right)\left(\mathcal{H}_{m}\rho_{t}-\mathcal{H}_{m}\hat{\rho}_{t}\right)^{\intercal}\right].$

It is also linked to the approximation by
$C_{t}^{H^{m}}=\tilde{\mathbb{E}}\left[\left(H_{t}^{m}-\bar{H}_{t}^{m}\right)\left(H_{t}^{m}-\bar{H}_{t}^{m}\right)^{\intercal}\right]=\mathcal{H}_{m}\mathcal{\hat{C}}_{t}\mathcal{H}_{m}^{*}$,
demonstrating a consistent error estimate in the statistical filtering model.
Similar conclusion can be reached for the accurate prediction in the model
covariance prediction for $R_{t}$ based on the cubic observation operator
$\mathcal{H}_{v}$.

######  Remark.

Still, the statistical consistency in the analysis step does not guarantee the
consistency in the entire two-step updating procedure in (3.15). In
particular, the forecast models of (3.12) and (3.16) have the following
updating equations during the forecast step

$\begin{aligned}
\partial_{t}\hat{\rho}_{t}=&\>\mathcal{L}_{t}^{*}\hat{\rho}_{t}\\\
\partial_{t}\tilde{\rho}_{t}=&\>\mathcal{L}_{t}^{*}\tilde{\rho}_{t}\end{aligned}\;\Rightarrow\;\begin{aligned}
\partial_{t}\left(\mathcal{H}\hat{\rho}_{t}\right)=&\>\mathcal{H}\mathcal{L}_{t}^{*}\hat{\rho}_{t}=\int\left(\mathcal{L}_{t}H\right)\hat{\rho}_{t}\left(z\right)\mathrm{d}z\\\
\partial_{t}\left(\tilde{\mathbb{E}}H\right)=&\>\tilde{\mathbb{E}}\mathcal{L}_{t}H=\int\left(\mathcal{L}_{t}H\right)\tilde{\rho}_{t}\left(z\right)\mathrm{d}z\end{aligned}$

where the generator is defined with $y_{t}=\left(\bar{u}_{t},R_{t}\right)$

$\mathcal{L}_{t}\left(y_{t}\right)\varphi=\nabla_{z}\varphi\cdot\left\\{\left[H\left(z\right)-R_{t}\right]+L\left(\bar{u}_{t}\right)z\right\\}+\frac{1}{2}\Sigma_{t}^{2}:\nabla_{z}\nabla_{z}\varphi.$

The analysis step update only gives consistent first two moments of $H$, while
higher moments may be included in $\mathcal{L}_{t}H$. It is not guaranteed
that the forecast model can give consistent forecast in
$\mathcal{H}\hat{\rho}_{t}$ and $\tilde{\mathbb{E}}H$ as well as their
covariances. More work is still needed for the complete consistency analysis
combining the approximations in both the forecast and analysis step of the
filtering method.

##  4 Stability and convergence of the statistical filtering model

Here, we discuss the long-time performance of the full filtering problem for
the statistical filtering system (3.1) for finding the optimal filter PDF
$\hat{\rho}_{t}$ based on statistical observations
$\mathcal{G}_{t}=\sigma\left\\{y_{s},s\leq t\right\\}$. In Section 3, it shows
that the statistical filtering model (3.16) constitutes the approximate filter
solution $\tilde{\rho}_{t}$ with consistent mean and covariance, $\bar{H}_{t}$
and $C_{t}^{H}$, in the analysis step update. We show further here that the
full filter approximation of the observation function $\bar{H}_{t}$ will
approach the optimal filter $\mathcal{H}\hat{\rho}_{t}$ at the long-time limit
as $t\rightarrow\infty$. This guarantees the stable performance of the
proposed new filtering strategy.

###  4.1 Complete statistical filter equations based on the observation
operator

We consider the optimal filter solution based on the conditional Gaussian
model (3.1). The finite-dimensional statistical states
$\mathcal{H}\hat{\rho}_{t}\in\mathbb{R}^{p}$ and
$\hat{\mathcal{C}}_{t}^{\mathcal{H}}=\mathcal{H}\hat{\mathcal{C}}_{t}\mathcal{H}^{*}\in\mathbb{R}^{p\times
p}$ under the observation operator are solved by the Kalman-Bucy equations
(3.12) as

$\displaystyle\mathrm{d}\left(\mathcal{H}\hat{\rho}_{t}\right)=$
$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle\mathrm{d}t+\hat{\mathcal{C}}_{t}^{\mathcal{H}}\Gamma^{-2}\left\\{\mathrm{d}y_{t}-\left[\mathcal{H}\hat{\rho}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(4.1) $\displaystyle\mathrm{d}\hat{\mathcal{C}}_{t}^{\mathcal{H}}=$
$\displaystyle\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\mathcal{C}}_{t}H^{\intercal}\right\rangle+\left\langle\hat{\mathcal{C}}_{t}H,\mathcal{L}_{t}\left(y_{t}\right)H^{\intercal}\right\rangle\right]\mathrm{d}t-\hat{\mathcal{C}}_{t}^{\mathcal{H}}\Gamma^{-2}\hat{\mathcal{C}}_{t}^{\mathcal{H}}\mathrm{d}t.$

Above, $\left\langle F,G\right\rangle=\int
F\left(z\right)G\left(z\right)\mathrm{d}z$ denotes the componentwise
integration of the product of matrix-valued functions $F,G$. For simplicity,
we use constant observation noise, $\Gamma_{t}\equiv\Gamma$. In the first
equation for the mean, we rewrite the forecast step dynamics as

$\mathcal{H}\mathcal{L}_{t}^{*}\left(y_{t}\right)\hat{\rho}_{t}=\int
H\left(z\right)\mathcal{L}_{t}^{*}\left(y_{t}\right)\hat{\rho}_{t}\left(z\right)\mathrm{d}z=\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle.$

Similarly in the covariance equation, we rewrite using the definition of
$\hat{\mathcal{C}}_{t}$ in (3.10) under the conditional measure
$\mu_{t}\left(\cdot\right)=\mathbb{P}\left(\rho\in\cdot\mid\mathcal{G}_{t}\right)$

$\displaystyle\mathcal{H}\mathcal{L}_{t}^{*}\left(y_{t}\right)\mathcal{\hat{C}}_{t}\mathcal{H}^{*}$
$\displaystyle=\mathbb{E}_{\rho\sim\mu_{t}}\int
H\left(x\right)\mathcal{L}_{t}^{*}\left(\rho-\hat{\rho}_{t}\right)\left(x\right)\mathrm{d}x\int\left(\rho-\hat{\rho}_{t}\right)\left(z\right)H\left(z\right)^{\intercal}\mathrm{d}z=\left\langle\mathcal{L}_{t}H,\hat{\mathcal{C}}_{t}H^{\intercal}\right\rangle,$
$\displaystyle\mathcal{H}\mathcal{\hat{C}}_{t}\mathcal{L}_{t}\left(y_{t}\right)\mathcal{H}^{*}$
$\displaystyle=\mathbb{E}_{\rho\sim\mu_{t}}\int
H\left(z\right)\left(\rho-\hat{\rho}_{t}\right)\left(z\right)\mathrm{d}z\int\left(\rho-\hat{\rho}_{t}\right)\left(x\right)\mathcal{L}_{t}H\left(x\right)^{\intercal}\mathrm{d}x=\left\langle\hat{\mathcal{C}}_{t}H,\mathcal{L}_{t}H^{\intercal}\right\rangle,$

where
$\hat{\mathcal{C}}_{t}H\left(x\right)=\mathbb{E}_{\rho}\left[\left(\rho-\hat{\rho}_{t}\right)\left(x\right)\int\left(\rho-\hat{\rho}_{t}\right)\left(z\right)H\left(z\right)\mathrm{d}z\right]\in\mathbb{R}^{d}$
can be viewed as a unnormalized density.

Correspondingly, consider the approximating filter model (3.16) with
$\Sigma_{t}\equiv 0$. Let $\tilde{\rho}_{t}$ be the PDF after the symmetric
projection (3.23), that is, satisfying vanishing third-order moments
$\tilde{\mathbb{E}}\left[H_{n}^{\prime}H_{p}^{\prime}H_{q}^{\prime}\right]=0$
for all $n,p,q$. The moments
$\bar{H}_{t}=\tilde{\mathbb{E}}\left[H\left(\tilde{Z}_{t}\right)\right]$ and
$C_{t}^{H}=\tilde{\mathbb{E}}\left[H_{t}^{\prime}\left(\tilde{Z}_{t}\right)H_{t}^{\prime}\left(\tilde{Z}_{t}\right)^{\intercal}\right]$
according to the observation function $H$ satisfy the following equations (by
combining the forecast step update with generator $\mathcal{L}_{t}$ and the
analysis step dynamics (3.22))

$\displaystyle\mathrm{d}\bar{H}_{t}=$
$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle\mathrm{d}t+C_{t}^{H}\Gamma^{-2}\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(4.2) $\displaystyle\mathrm{d}C_{t}^{H}=$
$\displaystyle\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}H_{t}^{\prime\intercal}\right\rangle+\left\langle\tilde{\rho}_{t}H_{t}^{\prime},\mathcal{L}_{t}\left(y_{t}\right)H^{\intercal}\right\rangle\right]\mathrm{d}t-C_{t}^{H}\Gamma^{-2}C_{t}^{H}\mathrm{d}t,$

where $H_{t}^{\prime}\left(z\right)=H\left(z\right)-\bar{H}_{t}$.

First, we introduce the following assumptions for the approximating filter PDF
$\tilde{\rho}_{t}$ approximating the optimal filter solution $\hat{\rho}_{t}$
as random fields according to the same observation process
$Y=\left\\{y_{t},t\geq 0\right\\}$:

######  Assumption.

Assume that the optimal filter (3.12) and the approximating filter model
(3.16) have a continuous probability density functions $\hat{\rho}_{t}$ and
$\tilde{\rho}_{t}$ respectively that satisfy the following conditions:

  * •

Unique equilibrium solutions
$\left(\mathcal{H}\hat{\rho}_{\infty},\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)$
and $\left(\mathcal{H}\tilde{\rho}_{\infty},C_{\infty}^{H}\right)$ exist to
the systems (4.1) and (4.2) respectively.

  * •

There exist deterministic matrices
$L_{\infty}^{m},L_{\infty}^{v}\in\mathbb{R}^{p\times p}$. The generator
$\mathcal{L}_{t}$ (2.15) reaches the same statistical limit under $H$ w.r.t.
both $\hat{\rho}_{t}$ and $\tilde{\rho}_{t}$, that is,

$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle\rightarrow
L_{\infty}^{m}\mathcal{H}\hat{\rho}_{\infty},$
$\displaystyle\;\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\mathcal{C}}_{t}H^{\intercal}\right\rangle\rightarrow
L_{\infty}^{v}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}},$ (4.3)
$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle\rightarrow
L_{\infty}^{m}\mathcal{H}\tilde{\rho}_{\infty},$
$\displaystyle\;\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}H_{t}^{\prime\intercal}\right\rangle\rightarrow
L_{\infty}^{v}C_{\infty}^{H},$

a.s. as $t\rightarrow\infty$ given the observation data $y_{t}$.

  * •

Further, the real parts of eigenvalues of the limit matrices
$L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}$ and
$L_{\infty}^{v}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}$ are all
negative.

The first condition in assumption (4.3) guarantees that both filter equations
will finally converge to finite solutions. The second condition requires
consistent first and second-order moments under $H$ w.r.t. the optimal and
approximate model PDF, such as for the covariance as $t\rightarrow\infty$

$\displaystyle\int\mathcal{L}_{t}\left(y_{t}\right)H\left(z\right)\hat{\mathcal{C}}_{t}H\left(z\right)^{\intercal}\mathrm{d}z\rightarrow$
$\displaystyle L_{\infty}^{v}\int
H\left(z\right)\left(\hat{\mathcal{C}}_{\infty}H^{\intercal}\right)\left(z\right)\mathrm{d}z,$
$\displaystyle\int\mathcal{L}_{t}\left(y_{t}\right)H_{t}^{\prime}\left(z\right)H_{t}^{\prime}\left(z\right){}^{\intercal}\tilde{\rho}_{t}\left(z\right)\mathrm{d}z\rightarrow$
$\displaystyle L_{\infty}^{v}\int
H_{\infty}^{\prime}\left(z\right)H_{\infty}^{\prime}\left(z\right){}^{\intercal}\tilde{\rho}_{\infty}\left(z\right)\mathrm{d}z.$

In addition, we may further introduce the convergence rate, that is, there
exist constants $\bar{\lambda}>0$ and $K>0$ such that a.s.

$\displaystyle\left|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle-
L_{\infty}^{m}\mathcal{H}\hat{\rho}_{\infty}\right|\leq Ke^{-\bar{\lambda}t},$
$\displaystyle\quad\left\|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\mathcal{C}}_{t}H^{\intercal}\right\rangle-
L_{\infty}^{v}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|\leq
Ke^{-\bar{\lambda}t},$ (4.4)
$\displaystyle\left|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle-
L_{\infty}^{m}\mathcal{H}\tilde{\rho}_{\infty}\right|\leq
Ke^{-\bar{\lambda}t},$
$\displaystyle\quad\left\|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}H_{t}^{\prime\intercal}\right\rangle-
L_{\infty}^{v}C_{\infty}^{H}\right\|\leq Ke^{-\bar{\lambda}t},$

where $\left\|\cdot\right\|$ is the matrix norm. And the third condition asks
that the filter solutions will be stabilized at the long time limit. Next, we
ask the limit behaviour in the mean and covariance
$\left(\bar{H}_{t},C_{t}^{H}\right)$ from the approximation model (4.2) in
comparison with the optimal filter solution
$\left(\mathcal{H}\hat{\rho}_{t},\hat{\mathcal{C}}_{t}^{\mathcal{H}}\right)$
from (4.1).

###  4.2 Asymptotic stability of the equilibrium covariance matrix

With the assumption (4.3), we have consistent equilibrium covariance at
$t\rightarrow\infty$ in the two model solutions

$\displaystyle
L_{\infty}^{v}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}L_{\infty}^{v\intercal}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}$
$\displaystyle=0,$ $\displaystyle
L_{\infty}^{v}C_{\infty}^{H}+C_{\infty}^{H}L_{\infty}^{v\intercal}-C_{\infty}^{H}\Gamma^{-2}C_{\infty}^{H}$
$\displaystyle=0.$

Uniqueness of the solution directly implies that the final equilibrium
covariances satisfy $\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}=C_{\infty}^{H}$
with no randomness. Further, we have that the covariance $C_{t}^{H}$ in the
approximating filter will approach the optimal equilibrium covariance
$\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}$ as described in the following
result.

######  Lemma 8.

Suppose that Assumption Assumption is satisfied, $C_{t}^{H}$ is the covariance
solution to the statistical filtering model (4.2) and
$\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}$ is the unique equilibrium solutions
to the optimal model (4.1). Then there is

$\left\|C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|\rightarrow
0,$ (4.5)

a.s. as $t\rightarrow\infty$. Further, the convergence rate will be
exponential if (4.4) is also satisfied.

###### Proof.

Combining the covariance equation in (4.2) and the equilibrium equation of
(4.1), we have

$\displaystyle\mathrm{d}\left(C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)$
$\displaystyle=\left[L_{\infty}^{v}-\frac{1}{2}\left(C_{t}^{H}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)\Gamma^{-2}\right]\left(C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)\mathrm{d}t$
$\displaystyle+\left(C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)\left[L_{\infty}^{v\intercal}-\frac{1}{2}\Gamma^{-2}\left(C_{t}^{H}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)\right]\mathrm{d}t$
$\displaystyle+\left[\left\langle\mathcal{L}_{t}H,\tilde{\rho}_{t}H_{t}^{\prime}{}^{\intercal}\right\rangle-
L_{\infty}^{v}C_{t}^{H}\right]\mathrm{d}t+\left[\left\langle\tilde{\rho}_{t}H_{t}^{\prime},\mathcal{L}_{t}H^{\intercal}\right\rangle-
C_{t}^{H}L_{\infty}^{v\intercal}\right]\mathrm{d}t.$

For the last row of the above equation, using the uniqueness of the solution
we have $C_{t}^{H}\rightarrow C_{\infty}^{H}$ as $t\rightarrow\infty$. Then
denote

$F_{t}=\left\langle\mathcal{L}_{t}H,\tilde{\rho}_{t}H_{t}^{\prime}{}^{\intercal}\right\rangle-
L_{\infty}^{v}C_{t}^{H}.$

We get $\left|F_{t}\right|\rightarrow 0$ as $t\rightarrow\infty$ by using

$\left\|\left\langle\mathcal{L}_{t}H,\tilde{\rho}_{t}H_{t}^{\prime}{}^{\intercal}\right\rangle-
L_{\infty}^{v}C_{t}^{H}\right\|\leq\left\|\left\langle\mathcal{L}_{t}H,\tilde{\rho}_{t}H_{t}^{\prime}{}^{\intercal}\right\rangle-
L_{\infty}^{v}C_{\infty}^{H}\right\|+\left\|L_{\infty}^{v}\left(C_{\infty}^{H}-C_{t}^{H}\right)\right\|\rightarrow
0.$

Also by taking $\
\lambda_{v}=\min\left\\{\mathrm{Re}\lambda:\lambda\>\mathrm{is\>the\>eighenvalue\>of}\>-L_{\infty}^{v}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right\\}$,
we have $\lambda_{v}>0$ from Assumption Assumption. This implies
$\left\|e^{\int_{s}^{t}\left[L_{\infty}^{v}-\frac{1}{2}\left(C_{\tau}^{H}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right)\Gamma^{-2}\right]\mathrm{d}\tau}\right\|\leq
Ke^{-\lambda_{v}\left(t-s\right)}$ . Together with $F_{t}$ vanishing as
$t\rightarrow\infty$, we have for any $t>T$

$\left\|C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|\leq
K\left\|C_{T}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|e^{-2\lambda_{v}\left(t-T\right)}+2K\int_{T}^{t}e^{-2\lambda_{v}\left(t-s\right)}\left|F_{s}\right|\mathrm{d}s\leq
K_{1}e^{-2\lambda_{v}t}+K_{2}\sup_{s\geq T}\left|F_{s}\right|.$

Therefore, by first letting $t\rightarrow\infty$ then letting
$T\rightarrow\infty$, we reach a.s.
$\lim_{t\rightarrow\infty}\left\|C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|=0$.

Further, if we assume exponential convergence rate $\bar{\lambda}$ in the
covariances under the generator $\mathcal{L}_{t}$ as in (4.4), we can have
exponential convergence rate in both $\hat{\mathcal{C}}_{t}$ and $C_{t}^{H}$
as $t\rightarrow\infty$ a.s.

$\left\|\hat{\mathcal{C}}_{t}^{\mathcal{H}}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|\leq
K_{H}e^{-\min\left\\{\lambda_{v},\bar{\lambda}\right\\}t},\;\left\|C_{t}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|\leq
K_{H}e^{-\min\left\\{\lambda_{v},\bar{\lambda}\right\\}t}.$ (4.6)

∎

In addition, from the definition of the optimal filter solution (3.11), the
covariance is defined as

$\displaystyle\mathrm{tr}\mathbb{E}\hat{\mathcal{C}}_{t}^{\mathcal{H}}$
$\displaystyle=\mathbb{E}_{\mu_{t}}\mathrm{tr}\left[\left(\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right)\left(\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right)^{\intercal}\right]$
$\displaystyle=\mathbb{E}_{\mu_{t}}\mathrm{tr}\left[\left(\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right)^{\intercal}\left(\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right)\right]$
$\displaystyle=\mathbb{E}_{\mu_{t}}\left[\left|\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right|^{2}\right].$

Notice that $\mathcal{H}\hat{\rho}_{t}$ and
$\hat{\mathcal{C}}_{t}^{\mathcal{H}}$ are still random field in
$\mathcal{G}_{t}$. Under Assumption Assumption
$\hat{\mathcal{C}}_{t}^{\mathcal{H}}\rightarrow\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}$
a.s. with the observations in $\mathcal{G}_{t}$, we have as
$t\rightarrow\infty$

$\mathbb{E}\left[\left|\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right|^{2}\right]=\mathrm{tr}\mathbb{E}\left[\hat{\mathcal{C}}_{t}^{\mathcal{H}}\right]\rightarrow\mathrm{tr}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}.$
(4.7)

This confirms that the total uncertainty at equilibrium in the optimal filter
solution $\mathcal{H}\hat{\rho}_{t}$ is estimated by the total variance
$\mathrm{tr}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}$.

###  4.3 Convergence of the statistical state under the observation operator

Next, we consider the convergence of the statistical observation function
$\bar{H}_{t}$ from the approximating filter to the optimal filter solution
$\mathcal{H}\hat{\rho}_{t}$. We have the long-term stability in the
statistical solution in the following theorem.

######  Theorem 9.

Suppose that Assumption Assumption holds and the covariances goes to the same
deterministic limit

$\hat{\mathcal{C}}_{t}^{\mathcal{H}}\rightarrow\hat{\mathcal{C}}_{\infty}^{\mathcal{H}},\quad
C_{t}^{H}\rightarrow\hat{\mathcal{C}}_{\infty}^{\mathcal{H}},$

a.s. as $t\rightarrow\infty$. Then there is

$\mathbb{E}\left[\left|\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right|^{2}\right]\rightarrow
0,$ (4.8)

as $t\rightarrow\infty$. Furthermore, assume exponential convergence rate in
$\hat{\mathcal{C}}_{t}^{\mathcal{H}},C_{t}^{H}$ as in (4.6), and let
$\lambda=\min\left\\{\lambda_{m},\lambda_{v}\right\\}$ where $\lambda_{m}$ and
$\lambda_{v}$ are the minimums of the real parts of the eigenvalues of
$-L_{\infty}^{m}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}$ and
$-L_{\infty}^{v}+\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}$
respectively. There is also exponential convergence as

$\mathbb{E}\left[\left|\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right|^{2}\right]\leq
K_{H}e^{-\lambda_{1}t},$ (4.9)

with some $\lambda_{1}\leq\min\left\\{\lambda,\bar{\lambda}\right\\}$ and
$K_{H}$ a constant only dependent on the observation function $H$.

###### Proof.

By taking the difference of the mean equations in (4.1) and (4.2), we have

$\displaystyle\mathrm{d}\left(\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right)$
$\displaystyle=\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)\left(\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right)\mathrm{d}t+\left(\hat{\mathcal{C}}_{t}^{\mathcal{H}}-C_{t}^{H}\right)\Gamma^{-2}\left[\mathrm{d}y_{t}-h_{t}\left(y_{t}\right)\mathrm{d}t\right]$
$\displaystyle+\left[\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-\hat{\mathcal{C}}_{t}^{\mathcal{H}}\right)\Gamma^{-2}\mathcal{H}\hat{\rho}_{t}-\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-C_{t}^{H}\right)\Gamma^{-2}\bar{H}_{t}\right]\mathrm{d}t$
$\displaystyle+\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle-
L_{\infty}^{m}\mathcal{H}\hat{\rho}_{t}\right]\mathrm{d}t+\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle-
L_{\infty}^{m}\bar{H}_{t}\right]\mathrm{d}t.$

By applying Itô’s formula to the above equation, there is

$\displaystyle\mathrm{d}\left[e^{-t\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right)\right]$
$\displaystyle=e^{-t\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left\\{\left(\hat{\mathcal{C}}_{t}^{\mathcal{H}}-C_{t}^{H}\right)\Gamma^{-2}\left[\mathrm{d}y_{t}-h_{t}\left(y_{t}\right)\mathrm{d}t\right]+G_{t}\left(y_{t}\right)\mathrm{d}t\right\\}$
$\displaystyle+e^{-t\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left[\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-\hat{\mathcal{C}}_{t}^{\mathcal{H}}\right)\Gamma^{-2}\mathcal{H}\hat{\rho}_{t}-\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-C_{t}^{H}\right)\Gamma^{-2}\bar{H}_{t}\right]\mathrm{d}t.$

Above, we denote the residual term as

$G_{t}\left(y_{t}\right)=\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle-
L_{\infty}^{m}\mathcal{H}\hat{\rho}_{t}\right]+\left[\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle-
L_{\infty}^{m}\bar{H}_{t}\right].$

By similar computation in Lemma 8 according to assumption (4.3), there is

$\left|G_{t}\right|\leq\begin{aligned}
\left|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle-
L_{\infty}^{m}\mathcal{H}\hat{\rho}_{\infty}\right|&+\left|L_{\infty}^{m}\left(\mathcal{H}\hat{\rho}_{\infty}-\mathcal{H}\hat{\rho}_{t}\right)\right|\\\
+\left|\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle-
L_{\infty}^{m}\bar{H}_{\infty}\right|&+\left|L_{\infty}^{m}\left(\bar{H}_{\infty}-\bar{H}_{t}\right)\right|\end{aligned}\leq
Ke^{-\bar{\lambda}t}\rightarrow 0$

a.s. as $t\rightarrow\infty$ using the uniqueness of the solutions
$\mathcal{H}\hat{\rho}_{t}\rightarrow\mathcal{H}\hat{\rho}_{\infty}$ and
$\bar{H}_{t}\rightarrow\bar{H}_{\infty}$ (with exponential decay in the
stronger convergence case (4.6)). Equivalently, the above SDE can be written
as (ignoring initial condition by assuming
$\mathcal{H}\hat{\rho}_{0}=\bar{H}_{0}$)

$\displaystyle\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}$
$\displaystyle=\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left\\{\left(\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\left[\mathrm{d}y_{s}-h_{s}\left(y_{s}\right)\mathrm{d}s\right]+G_{s}\left(y_{s}\right)\mathrm{d}s\right\\}$
$\displaystyle+\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left[\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-\hat{\mathcal{C}}_{s}^{\mathcal{H}}\right)\Gamma^{-2}\mathcal{H}\hat{\rho}_{s}-\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\bar{H}_{s}\right]\mathrm{d}s.$

Therefore, by taking the expectation for each term on the right hand side of
the above identity

$\displaystyle\mathbb{E}\left[\left|\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right|^{2}\right]$
$\displaystyle\leq
5\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\mathrm{d}y_{s}\right|^{2}$
$\displaystyle+5\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}G_{s}\left(y_{s}\right)\mathrm{d}s\right|^{2}$
$\displaystyle+5\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}h_{s}\left(y_{s}\right)\mathrm{d}s\right|^{2}$
$\displaystyle+5\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-\hat{\mathcal{C}}_{s}^{\mathcal{H}}\right)\Gamma^{-2}\mathcal{H}\hat{\rho}_{s}\mathrm{d}s\right|^{2}$
$\displaystyle+5\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\bar{H}_{s}\mathrm{d}s\right|^{2}.$

The second line above follows the same argument as in Lemma 8 and using
Cauchy-Schwarz inequality

$\displaystyle\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}G_{s}\left(y_{s}\right)\mathrm{d}s\right|^{2}$
$\displaystyle\leq\int_{0}^{t}\left\|e^{\frac{1}{2}\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\right\|^{2}\mathrm{d}s\mathbb{E}\int_{0}^{t}\left|e^{\frac{1}{2}\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}G_{s}\left(y_{s}\right)\right|^{2}\mathrm{d}s.$
$\displaystyle\leq\lambda_{m}^{-1}\mathbb{E}\int_{0}^{t}e^{-\lambda_{m}\left(t-s\right)}\left|G_{s}\left(y_{s}\right)\right|^{2}\mathrm{d}s$
$\displaystyle\leq Ke^{-\lambda_{1}t}\rightarrow 0.$

Similar results can be achieved for line three to five following the
convergence (or exponential convergence) of the integrants. Finally, for the
first line using the observation equation in (3.1), that is,
$\mathrm{d}y_{t}=\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Gamma\mathrm{d}B_{t}$,
there is

$\displaystyle\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\mathrm{d}y_{s}\right|^{2}$
$\displaystyle\leq$
$\displaystyle\mathbb{E}\left|\int_{0}^{t}e^{\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left(\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right)\Gamma^{-2}\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}s\right|^{2}$
$\displaystyle+$
$\displaystyle\mathbb{E}\int_{0}^{t}e^{2\left(t-s\right)\left(L_{\infty}^{m}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\right)}\left\|\hat{\mathcal{C}}_{s}^{\mathcal{H}}-C_{s}^{H}\right\|^{2}\mathrm{d}s$
$\displaystyle\leq$ $\displaystyle
K_{1}\int_{0}^{t}e^{-\lambda_{m}\left(t-s\right)}\left(\mathbb{E}\left\|\hat{\mathcal{C}}_{s}^{\mathcal{H}}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|^{2}+\mathbb{E}\left\|C_{s}^{H}-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\right\|^{2}\right)\mathrm{d}s$
$\displaystyle\leq$ $\displaystyle Ke^{-\lambda_{1}t}\rightarrow 0.$

Above in the second to the last inequality, we use the uniform boundedness of
$\mathbb{E}\left|\mathcal{H}\rho_{t}\right|$ and
$\mathbb{E}\left|h_{t}\left(y_{t}\right)\right|$ as $t\rightarrow\infty$, and
the last line uses the convergence (or exponential convergence) of the
covariance (4.5) or (4.6). Combining all the above bounds, we finally get
(4.8) and (4.9). ∎

Together with the equilibrium estimate in (4.7) combining the result in
Theorem 9, we can also get the same error estimate compared with the target
field as $t\rightarrow\infty$

$\mathbb{E}\left[\left|\mathcal{H}\rho_{t}-\bar{H}_{t}\right|^{2}\right]\leq\mathbb{E}\left[\left|\mathcal{H}\rho_{t}-\mathcal{H}\hat{\rho}_{t}\right|^{2}\right]+\mathbb{E}\left[\left|\mathcal{H}\hat{\rho}_{t}-\bar{H}_{t}\right|^{2}\right]\rightarrow\mathrm{tr}\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}.$

Thus, we show that the statistical filtering model solution $\bar{H}_{t}$
converges to the optimal filter solution with the same mean square error, and
exponential convergence is reached if the forecast model has exponential
convergence to the equilibrium.

As a final comment, we can further relax Assumption Assumption as there exist
deterministic uniformly continuous functions, $L_{\infty}^{m}\left(y\right)$
and $L_{\infty}^{v}\left(y\right)$, so that for any $y_{t}\rightarrow
y_{\infty},t\rightarrow\infty$

$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\hat{\rho}_{t}\right\rangle\rightarrow
L_{\infty}^{m}\left(y_{\infty}\right)\mathcal{H}\hat{\rho}_{\infty},$
$\displaystyle\;\left\langle\mathcal{L}_{t}H,\hat{\mathcal{C}}_{t}H^{\intercal}\right\rangle\rightarrow
L_{\infty}^{v}\left(y_{\infty}\right)\hat{\mathcal{C}}_{\infty}^{\mathcal{H}},$
(4.10)
$\displaystyle\left\langle\mathcal{L}_{t}\left(y_{t}\right)H,\tilde{\rho}_{t}\right\rangle\rightarrow
L_{\infty}^{m}\left(y_{\infty}\right)\mathcal{H}\tilde{\rho}_{\infty},$
$\displaystyle\;\left\langle\mathcal{L}_{t}H,\tilde{\rho}_{t}H_{t}^{\prime}{}^{\intercal}\right\rangle\rightarrow
L_{\infty}^{v}\left(y_{\infty}\right)C_{\infty}^{H},$

a.s. as $t\rightarrow\infty$. And the limiting matrices are uniformly bounded
by negative-definite matrices $A_{m},A_{v}$

$L_{\infty}^{m}\left(y\right)-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\preceq
A_{m}\prec 0,\quad
L_{\infty}^{v}\left(y\right)-\hat{\mathcal{C}}_{\infty}^{\mathcal{H}}\Gamma^{-2}\preceq
A_{v}\prec 0.$ (4.11)

In addition, non-zero noise in (3.16) can be included satisfying
$\Sigma_{t}\rightarrow 0$ as $t\rightarrow\infty$. Then, the same result
applies for the convergence of the observation mean and covariance at long-
time limit as in Lemma 8 and Theorem 9.

##  5 Ensemble approximation of the statistical filtering model

Finally, we discuss the construction of practical ensemble methods for solving
the statistical filtering model (3.16) with explicit model parameters that is
easy to implement. We demonstrate the accuracy of the numerical approximation
at the limit of large number of samples and small integration time step.

###  5.1 Numerical methods for the approximating filter with discrete
observations

Assume that the observation data comes at discrete times $t_{n}=\delta n$ with
a constant observation frequency $\delta$. We can introduce the linear
interpolation of
$y_{n}=\left(\bar{u}_{t_{n}},R_{t_{n}}\right)\in\mathbb{R}^{p}$ as

$\frac{\mathrm{d}y_{t}^{\delta}}{\mathrm{d}t}=\frac{y_{n+1}-y_{n}}{t_{n+1}-t_{n}}=\frac{\Delta
y_{n+1}}{\delta},$ (5.1)

during time interval $t\in\left[t_{n},t_{n+1}\right]$. We propose ensemble
algorithms to approximate the filtering distribution
$\tilde{\rho}_{t_{n}}\sim\tilde{Z}_{t_{n}}$ conditional on the statistical
observations $Y_{t}^{\delta}=\left\\{y_{s}^{\delta},s\leq t\right\\}$ based on
the statistical filter equation (3.16). First, $N$ independent particles,
$\tilde{\mathbf{Z}}_{t}=\left\\{\tilde{Z}_{t}^{\left(i\right)}\right\\}_{i=1}^{N}$,
are drawn to sample the initial distribution of the stochastic state. Then,
the particles are evolved according to the following SDE with drift terms
$a_{t}^{m},a_{t}^{v}$ and control gains $K_{t}^{m},K_{t}^{v}$

$\displaystyle\mathrm{d}\tilde{Z}_{t}^{\left(i\right)}=$
$\displaystyle\;L\left(\bar{u}_{t}^{N}\right)\tilde{Z}_{t}^{\left(i\right)}\mathrm{d}t+\varGamma\left(\tilde{Z}_{t}^{\left(i\right)}\tilde{Z}_{t}^{\left(i\right)\intercal}-R_{t}^{N}\right)\mathrm{d}t+\Sigma_{t}\mathrm{d}\tilde{W}_{t}^{\left(i\right)}$
(5.2)
$\displaystyle+a_{t}^{m}\left(\tilde{Z}_{t}^{\left(i\right)}\right)\mathrm{d}t+K_{t}^{m}\left(\tilde{Z}_{t}^{\left(i\right)}\right)\left\\{\mathrm{d}\bar{u}_{t}^{\delta}-\left[H^{m}\left(\tilde{Z}_{t}^{\left(i\right)}\right)+h_{m,t}\left(\bar{u}_{t}^{\delta}\right)\right]\mathrm{d}t-\Gamma_{m,t}\mathrm{d}\tilde{B}_{m,t}^{\left(i\right)}\right\\}$
$\displaystyle+a_{t}^{v}\left(\tilde{Z}_{t}^{\left(i\right)}\right)\mathrm{d}t+K_{t}^{v}\left(\tilde{Z}_{t}^{\left(i\right)}\right)\left\\{\mathrm{d}R_{t}^{\delta}-\left[H^{v}\left(\tilde{Z}_{t}^{\left(i\right)}\right)+h_{v,t}\left(\bar{u}_{t}^{\delta},R_{t}^{\delta}\right)\right]\mathrm{d}t-\Gamma_{v,t}\mathrm{d}\tilde{B}_{v,t}^{\left(i\right)}\right\\},$

where the expressions for $h_{m},h_{v}$ are defined in (3.3), $H^{m},H^{v}$
are defined in (3.4), and
$\tilde{B}_{m,t}^{\left(i\right)},\tilde{B}_{v,t}^{\left(i\right)}$ are
independent white noises. Above, in the first line of (5.2) for the forecast
step of the filter, the first two moments
$\left(\bar{u}_{t}^{N},R_{t}^{N}\right)$ can be explicitly solved by the
statistical equations according to the stochastic-statistical model (2.7)

$\displaystyle\frac{\mathrm{d}\bar{u}_{t}^{N}}{\mathrm{d}t}=$
$\displaystyle\>M\left(\bar{u}_{t}^{N}\right)+F_{t}+\mathbb{E}^{N}\left[H^{m}\left(\tilde{\mathbf{Z}}_{t}\right)\right],$
(5.3) $\displaystyle\frac{\mathrm{d}R_{t}^{N}}{\mathrm{d}t}=$
$\displaystyle\>L\left(\bar{u}_{t}^{N}\right)R_{t}^{N}+R_{t}^{N}L\left(\bar{u}_{t}^{N}\right)^{\intercal}+\Sigma_{t}\Sigma_{t}^{\intercal}+\mathbb{E}^{N}\left[H^{v}\left(\tilde{\mathbf{Z}}_{t}\right)\right]+\epsilon^{-1}\left(\mathbb{E}^{N}\left[\tilde{\mathbf{Z}}_{t}\tilde{\mathbf{Z}}_{t}^{\intercal}\right]-R_{t}^{N}\right).$

The expectation is computed through the empirical average of the ensemble,
$\mathbb{E}^{N}f\left(\tilde{\mathbf{Z}}\right)=\frac{1}{N}\sum_{i=1}^{N}f\left(\tilde{Z}^{\left(i\right)}\right)$.
In this way, the particle simulation of (5.2) can be carried out easily for
each individual sample $\tilde{Z}_{t}^{\left(i\right)}$, and the dependence on
the distribution of the whole interacting particles is only introduced through
the empirical average in the statistical equations (5.3). We summarize the
ensemble filtering strategy in Algorithm 1.

1:Get the interpolated sequence of observation (5.1) of the mean and
covariance
$y^{\delta}_{t}=\left\\{\bar{u}_{t}^{\delta},R_{t}^{\delta}\right\\},t\in\left[0,T\right]$;
and introduce discrete time integration step $\tau$, and the initial
distribution $\rho_{0}$ of the model state.

2:At the initial time $t=0$, draw an ensemble of samples
$\left\\{\tilde{Z}_{0}^{\left(i\right)}\right\\}_{i=1}^{N}$ from the initial
distribution $\rho_{0}$, and compute initial mean and covariance
$\left\\{\bar{u}^{N}_{0},R^{N}_{0}\right\\}$ consistent with the statistics
w.r.t. $\rho_{0}$.

3:for $n=0$ while $n<\left\lfloor T/\tau\right\rfloor$, during the time
updating interval $t\in\left[t_{n},t_{n+1}\right]$ with $t_{n}=n\tau$ do

4: Compute the gain functions $K^{m}_{t_{n}}$ and $K^{v}_{t_{n}}$ using (5.6)
and the associated drift terms $a^{m}_{t_{n}}$ and $a^{v}_{t_{n}}$.

5: Update the samples
$\left\\{\tilde{Z}_{t_{n+1}}^{\left(i\right)}\right\\}_{i=1}^{N}$ using (5.2)
with the statistical states $\left\\{\bar{u}^{N}_{t},R^{N}_{t}\right\\}$ and
observation data $y^{\delta}_{t}$.

6: Update the statistical mean and covariance
$\left\\{\bar{u}^{N}_{t_{n+1}},R^{N}_{t_{n+1}}\right\\}$ by integrating (5.3)
to the next time step using the average of all samples.

7:end for

Algorithm 1 Ensemble statistical filter with observations in mean and
covariance

######  Remark.

Solving the equations (5.2) may still demand high computational cost for
resolving the multiple nonlinear coupling terms in high dimension $d\gg 1$.
One potential approach to address the computational challenge is to adopt the
efficient random batch approach [47, 48] developed for the coupled models
(2.7). A detailed investigation of the numerical methods will be performed in
the follow-up research.

###  5.2 Construction of explicit model operators in the analysis step

In the second and third lines of (5.2) for the analysis step update of
filtering, we still need to propose explicit expressions for functions
$a_{t}^{m},K_{t}^{m}$ and $a_{t}^{v},K_{t}^{v}$ according to the observations
of the mean and covariance respectively. According to Theorem 7, the gain
function $K_{t}$ needs to be solved from equation (3.21), that is,

$-\nabla\cdot\left(K_{t}^{\intercal}\tilde{\rho}_{t}\right)=\tilde{\rho}_{t}\Gamma_{t}^{-2}\left(H\left(z\right)-\bar{H}_{t}\right).$

Then the drift function $a_{t}$ can be directly computed from the solution of
$K_{t}$ as

$a_{t}=\nabla\cdot\left(K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right)-K_{t}\Gamma_{t}^{2}\nabla\cdot
K_{t}^{\intercal}.$

In general, it is still difficult to find solutions of the above equations. By
multiplying $H$ on both sides and integrating about $z$, the identity for
$K_{t}$ implies a necessary condition

$\tilde{\mathbb{E}}\left[K_{t}^{\intercal}\nabla
H\right]=\Gamma_{t}^{-2}C_{t}^{H},$ (5.4)

where
$C_{t}^{H}=\tilde{\mathbb{E}}\left[\left(H\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right)\left(H\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right)^{\intercal}\right]$
is the covariance of $H$. Therefore, we can first design proper gain functions
$K_{t}$ by solving (5.4) according to the specific structures of
$H^{m}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$H^{v}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d^{2}}$ required in our problem in
(3.4)

$\displaystyle H_{k}^{m}\left(z\right)=$
$\displaystyle\sum_{m,n}\gamma_{kmn}z_{m}z_{n}=z^{\intercal}A_{k}z,$ (5.5)
$\displaystyle H_{kl}^{v}\left(z\right)=$
$\displaystyle\sum_{m,n}\gamma_{kmn}z_{m}z_{n}z_{l}+\gamma_{lmn}z_{m}z_{n}z_{k}=\left(z^{\intercal}A_{k}z\right)z_{l}+z_{k}\left(z^{\intercal}A_{l}z\right),$

for all $1\leq k,l\leq d$ where we rewrite the quadratic and cubic functions
using the symmetric coefficient matrix,
$A_{k}^{\intercal}=A_{k}\in\mathbb{R}^{d\times d}$, using the assumed
structural symmetry in the coupling coefficient $\gamma_{kmn}$. The resulting
gain functions are then constructed with the following specific expressions.

######  Proposition 10.

The matrix-valued functions $K_{t}^{m}=\tilde{K}^{m}\Gamma_{m,t}^{-2}$ and
$K_{t}^{v}=\tilde{K}^{v}\Gamma_{v,t}^{-2}$ with
$\tilde{K}^{m}\left(z\right)\in\mathbb{R}^{d\times d}$ and
$\tilde{K}^{v}\left(z\right)\in\mathbb{R}^{d\times d^{2}}$ in the following
expressions

$\displaystyle\tilde{K}_{j,k}^{m}\left(z\right)=$
$\displaystyle\frac{1}{2}z_{j}\left[\left(z^{\intercal}A_{k}z\right)-\bar{H}_{k}^{m}\right],$
(5.6) $\displaystyle\tilde{K}_{j,kl}^{v}\left(z\right)=$
$\displaystyle\frac{1}{3}z_{j}\left(z^{\intercal}A_{k}z\right)z_{l}+\frac{1}{3}z_{j}\left(z^{\intercal}A_{l}z\right)z_{k}-\frac{1}{3}\bar{H}_{kl}^{v}.$

for $1\leq k,l\leq d$ and $1\leq j\leq d$ satisfy the equation (5.4) according
to the structures of the functions $H^{m}$ and $H^{v}$ in the form of (5.5)
respectively, and
$\bar{H}^{m}=\tilde{\mathbb{E}}\left[H^{m}\left(\tilde{Z}\right)\right]$ and
$\bar{H}^{v}=\tilde{\mathbb{E}}\left[H^{v}\left(\tilde{Z}\right)\right]$.

The proof of Proposition 10 is put in Appendix B. The average terms,
$\bar{H}^{m},\bar{H}^{v}$, are already computed in the statistical equations
(5.3) thus no additional computational cost is needed. On the other hand, it
is noticed that (5.6) can only give a necessary condition for the gain
operators and may not guarantee the original identity for $K_{t}$ in general.
However, in the proof of Theorem 7, it shows that (5.4) is the main relation
used to derive the consistent analysis statistics on the mean of
$H\left(\tilde{Z}_{t}\right)$. Therefore, (5.6) provides a desirable candidate
for practical implementations of the algorithm concerning the consistency in
the leading moments.

###  5.3 Ensemble approximation in the discrete filtering model

At last, we discuss the accuracy in the finite ensemble approximation. Let
$\tilde{\rho}_{t}^{N}\left(z\right)=\frac{1}{N}\sum_{i=1}^{N}\delta\left(z-\tilde{Z}_{t}^{\left(i\right)}\right)$
be the random field from the empirical PDF of the ensemble approximation with
finite ensemble size $N$, and $\tilde{\rho}_{t}$ is the random field from
(3.16) conditional on the observations $\mathcal{G}_{t}$. Additional
assumptions are needed for the structures of the model.

######  Assumption 11.

We assume that the functions $M:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$L:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ in the mean and covariance
equations (5.10) are Lipschitz continuous, that is, there is a constant
$\beta>0$ so that

$\left|M\left(u\right)-M\left(v\right)\right|\leq\beta\left|u-v\right|,\quad\left\|L\left(u\right)-L\left(v\right)\right\|\leq\beta\left|u-v\right|.$

In addition, the coefficients in
$\Gamma_{k}\left(R\right)=\sum_{m,n}\gamma_{kmn}R_{mn}$ are uniformly bounded,
that is, there exists a constant $C>0$, so that for all $k,m,n$

$\left|\gamma_{kmn}\right|\leq C.$

First, we direct use the conclusion from the McKean-Vlasov limit of the
$N$-particle system [43, 25]. Under the above assumptions of the model
coefficients, the ensemble representation (5.2) will converge to the
continuous system (3.16) as $N\rightarrow\infty$. In particular, assume that
the a unique solution exists for the McKean-Vlasov SDE, and the initial
ensemble $\left\\{\tilde{Z}_{0}^{\left(i\right)}\right\\}_{i=1}^{N}$ is drawn
from i.i.d. random samples. Given the observations $\mathcal{G}_{t}$, we have
for any $\varphi\in C_{b}^{2}\left(\mathbb{R}^{d}\right)$

$\left\langle\varphi,\tilde{\rho}_{t}^{N}\right\rangle=\frac{1}{N}\sum_{i=1}^{N}\varphi\left(\tilde{Z}_{t}^{\left(i\right)}\right)\rightarrow\mathbb{E}\left[\varphi\left(\tilde{Z}_{t}\right)\mid\mathcal{G}_{t}\right]=\left\langle\varphi,\tilde{\rho}_{t}\right\rangle,$
(5.7)

a.s. as $N\rightarrow\infty$. Furthermore, there is the error estimate for the
empirical estimate
$\mathbb{E}^{N}\varphi=\frac{1}{N}\sum\varphi\left(\tilde{Z}_{t}^{\left(i\right)}\right)$

$\mathbb{E}\left[\sup_{0\leq t\leq
T}\left|\left\langle\varphi,\tilde{\rho}_{t}^{N}\right\rangle-\left\langle\varphi,\tilde{\rho}_{t}\right\rangle\right|^{2}\right]\leq\frac{C_{T}}{N}\left\|\varphi\right\|_{\infty}^{2}.$
(5.8)

Detailed proofs on (5.7) and (5.8) can be found in such as Chapter 9 of [2].

Next, we consider the finite ensemble and discrete time estimation of the
leading-order mean and covariance from the filter model. With the ensemble
approximation of the density function $\tilde{\rho}_{t}^{N}$, the mean and
covariance estimates are computed from the equations (5.3)

$\displaystyle\frac{\mathrm{d}\bar{u}_{t}^{N,\tau}}{\mathrm{d}t}=$
$\displaystyle\>M\left(\bar{u}_{\sigma\left(t\right)}^{N,\tau}\right)+F_{\sigma\left(t\right)}+\frac{1}{N}\sum_{i=1}^{N}H^{m}\left(\tilde{Z}_{\sigma\left(t\right)}^{\left(i\right)}\right),$
(5.9) $\displaystyle\frac{\mathrm{d}R_{t}^{N,\tau}}{\mathrm{d}t}=$
$\displaystyle\>L\left(\bar{u}_{\sigma\left(t\right)}^{N,\tau}\right)R_{\sigma\left(t\right)}^{N,\tau}+R_{\sigma\left(t\right)}^{N,\tau}L\left(\bar{u}_{\sigma\left(t\right)}^{N,\tau}\right)^{\intercal}+\Sigma_{\sigma\left(t\right)}^{2}+\frac{1}{N}\sum_{i=1}^{N}H^{v}\left(\tilde{Z}_{\sigma\left(t\right)}^{\left(i\right)}\right),$

where a finite time integration scheme is also applied with
$\sigma\left(t\right)=n\tau$ for $t\in\left[t_{n},t_{n+1}\right]$. Above, we
neglect the last relaxation term with $\epsilon$ since it will automatically
vanish with the resulting consistency. For clarity, we show below again the
continuous forecast model from (2.7)

$\displaystyle\frac{\mathrm{d}\bar{u}_{t}}{\mathrm{d}t}=$
$\displaystyle\>M\left(\bar{u}_{t}\right)+F_{t}+\tilde{\mathbb{E}}H^{m}\left(\tilde{Z}_{t}\right),$
(5.10) $\displaystyle\frac{\mathrm{d}R_{t}}{\mathrm{d}t}=$
$\displaystyle\>L\left(\bar{u}_{t}\right)R_{t}+R_{t}L\left(\bar{u}_{t}\right)^{\intercal}+\Sigma_{t}^{2}+\tilde{\mathbb{E}}H^{v}\left(\tilde{Z}_{t}\right).$

Notice that $\bar{u}_{t}^{N},R_{t}^{N}$ and $\bar{u}_{t},R_{t}$ are stochastic
processes due to the random samples
$\left\\{\tilde{Z}_{t}^{\left(i\right)}\right\\}$ and $\tilde{Z}_{t}$ is still
dependent on the observations $\mathcal{G}_{t}$. We have the following
convergence result for the finite ensemble $N$ and finite time step $\tau$
approximation to the continuous model.

######  Proposition 12.

If Assumption 11 is satisfied, under the same initial condition the
statistical solution of the ensemble approximation model (5.9) with discrete
time step $\tau$ converges to that of the continuous model (5.10) with

$\displaystyle\mathbb{E}\left[\sup_{n\tau\leq
T}\left|\bar{u}_{t_{n}}^{N}-\bar{u}_{t_{n}}\right|^{2}\right]$
$\displaystyle\leq\left(C_{1,T}\tau+\frac{C_{2,T}}{N}\right)\left\|H^{m}\right\|_{\infty},$
(5.11) $\displaystyle\mathbb{E}\left[\sup_{n\tau\leq
T}\left|R_{t_{n}}^{N}-R_{t_{n}}\right|^{2}\right]$
$\displaystyle\leq\left(C_{1,T}^{\prime}\tau+\frac{C_{2,T}^{\prime}}{N}\right)\left(\left\|H^{m}\right\|_{\infty}+\left\|H^{v}\right\|_{\infty}\right).$

###### Proof.

First, consider the mean equations from the same initial state, we have

$\displaystyle\bar{u}_{t}^{N,\tau}-\bar{u}_{t}=$
$\displaystyle\int_{0}^{t}\left[M\left(\bar{u}_{\sigma\left(s\right)}^{N,\tau}\right)-M\left(\bar{u}_{s}\right)\right]\mathrm{d}s+\int_{0}^{t}\left[F_{\sigma\left(s\right)}-F_{s}\right]\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\left[\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{s}\right\rangle\right]\mathrm{d}s.$

Therefore, using Hölder’s inequality and Lipschitz condition for $M$ there is

$\displaystyle\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{t}^{N,\tau}-\bar{u}_{t}\right|^{2}$ $\displaystyle\leq
3T\beta\mathbb{E}\int_{0}^{T}\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}\right|^{2}\mathrm{d}s$
$\displaystyle+3T\mathbb{E}\int_{0}^{T}\left|\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{s}\right\rangle\right|^{2}\mathrm{d}s+C_{T}\tau$
$\displaystyle\leq C_{1}\int_{0}^{T}\mathbb{E}\sup_{s^{\prime}\leq
s}\left|\bar{u}_{s^{\prime}}^{N,\tau}-\bar{u}_{s^{\prime}}\right|^{2}\mathrm{d}s+C_{2}\tau\left\|H^{m}\right\|_{\infty}^{2}+\frac{C_{3}}{N}\left\|H_{m}\right\|_{\infty}^{2}.$
(5.12)

In the first term above on the right hand side, we follow the same procedure
to estimate the error in the corresponding continuous solution
$\bar{u}_{s}^{N,\tau}$ compared to the time discretization solution
$\bar{u}_{\sigma\left(s\right)}^{N,\tau}$

$\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}\right|^{2}\leq\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}^{N,\tau}\right|^{2}+\left|\bar{u}_{s}^{N,\tau}-\bar{u}_{s}\right|^{2}\leq
C\tau\left\|H^{m}\right\|_{\infty}^{2}+\left|\bar{u}_{s}^{N,\tau}-\bar{u}_{s}\right|^{2}.$

And in the second line for the term related to $H^{m}$, combining the discrete
time estimate with (5.8) gives

$\mathbb{E}\left[\sup_{t\leq T}\left|\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{t}\right\rangle\right|^{2}\right]\leq\left(C\tau^{2}+\frac{C_{T}}{N}\right)\left\|H^{m}\right\|_{\infty}^{2}.$

Therefore, applying Grönwall’s inequality to (5.12), we get the mean estimate
in (5.11).

Next, under a similar fashion, we can compute from the covariance equation and
using the Lipschitz condition for $L$

$\displaystyle\mathbb{E}\sup_{t\leq T}\left|R_{t}^{N,\tau}-R_{t}\right|^{2}$
$\displaystyle\leq C_{1}\beta\mathbb{E}\sup_{t\leq
T}\left\|\bar{u}_{t}\right\|_{\infty}^{2}\int_{0}^{T}\left|R_{\sigma\left(s\right)}^{N,\tau}-R_{s}\right|^{2}\mathrm{d}s$
$\displaystyle+C_{2}\beta\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{\sigma\left(t\right)}^{N,\tau}-\bar{u}_{t}\right|^{2}\sup_{t\leq
T}\left\|R_{\sigma\left(t\right)}\right\|^{2}$
$\displaystyle+C_{3}\mathbb{E}\int_{0}^{T}\left|\left\langle
H^{v},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{v},\tilde{\rho}_{s}\right\rangle\right|^{2}\mathrm{d}s+C_{T}\tau.$

Using the uniform boundedness of $\bar{u}_{t},R_{t}$ and (5.8) for $H^{v}$
together with the previous estimate for $\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{\sigma\left(t\right)}^{N,\tau}-\bar{u}_{t}\right|^{2}$, we
reach the final covariance estimate in (5.11). ∎

Theorem 12 shows that the discrete approximation scheme of the approximating
filter model can recover the key model statistics in mean and covariance. It
again demonstrates the central role in achieving accurate prediction of the
nonlinear observation functions $H^{m},H^{v}$ in the filtering method.

##  6 Summarizing discussions

We developed a systematic statistical filtering strategy that enables
effective ensemble approximation of non-Gaussian probability distributions of
multiscale turbulent states using observations in the leading-order moments of
the mean and covariance. The filtering model is based on a closed set of
coupled stochastic-statistical equations established for modeling general
turbulent dynamical systems involving nonlinear coupling. The non-Gaussian
features then can be characterized by a McKean-Vlasov SDE taking into account
both the stochastic forecast equation and corrections according to the
observed first two moments. Importantly, the proposed McKean-Vlasov SDE for
the finite-dimensional stochastic state does not require the explicit
computation of the infinite-dimensional probability distribution, but just
relies on the feedbacks from the leading moments that can be computed from the
associated statistical equations. This leads to straightforward numerical
algorithms using ensemble approximations. The stochastic–statistical
formulation offers a flexible approach for recovering essential model
statistics, making it applicable to a wide range of problems in uncertainty
quantification and data assimilation. In the immediate applications of this
research, the performance of the new filtering strategy will be tested on a
series of turbulent systems, starting from prototype models to realistic
applications in really high-dimensional systems. Also, currently,we are only
able to show the statistical consistency of the approximating filter in the
first two moments of the observation function under restricted conditions.
Further explorations exploiting specific model structures such as the
conservation properties will be used to provide a thorough understanding of
the approximation skill of the filter predictions.

## Acknowledge

The research of J.-G. L. is partially supported by the National Science
Foundation (NSF) under Grant No. DMS-2106988. The research of D. Q. is
partially supported by Office of Naval Research (ONR) Grant No.
N00014-24-1-2192.

##  Appendix A General background about filtering

Here, we summarize the useful results needed in the main text of this paper
following mostly [2, 36, 20].

###  A.1 Filtering equations for general stochastic systems

Let $\left(\Omega,\mathcal{F},\mathbb{P}\right)$ be the complete probability
space. The _signal process_ $u_{t}\in\mathbb{R}^{d}$ and the _observation
process_ $y_{t}\in\mathbb{R}^{p}$ are defined on the probability space
satisfying the following SDEs

$\displaystyle\mathrm{d}u_{t}$
$\displaystyle=F\left(u_{t}\right)\mathrm{d}t+\Sigma\mathrm{d}W_{t},\quad
u_{t=0}=u_{0},$ (A1a) $\displaystyle\mathrm{d}y_{t}$
$\displaystyle=H\left(u_{t}\right)\mathrm{d}t+\mathrm{d}B_{t},\quad
y_{t=0}=y_{0},$ (A1b)

where $F:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$H:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ are bounded and globally Lipschitz
continuous functions, and $W_{t}\in\mathbb{R}^{s},B_{t}\in\mathbb{R}^{p}$ are
independent standard Wiener processes with matrix coefficient
$\Sigma\in\mathbb{R}^{d\times s}$. The aim of the general filtering problem is
to determine the conditional probability distribution $\mu_{t}$ of the signal
process, $u_{t}$, given the accumulated observation process,
$Y_{t}=\left\\{y_{s},s\leq t\right\\}$.

Define the observation filtration $\mathcal{G}_{t}=\sigma\left\\{y_{s},s\leq
t\right\\}$. The random conditional distribution
$\mu_{t}:\mathbb{R}^{d}\times\Omega\rightarrow\left[0,1\right]$ is defined as
the $\mathcal{P}\left(\mathbb{R}^{d}\right)$-valued stochastic process which
is measurable w.r.t. $\mathcal{G}_{t}$, so that for any function $\varphi\in
C_{b}^{2}\left(\mathbb{R}^{d}\right)$ a.s.

$\mathbb{E}\left[\varphi\left(u_{t}\right)\mid\mathcal{G}_{t}\right]=\mu_{t}\left(\varphi\right)\coloneqq\int_{\mathbb{R}^{d}}\varphi\left(u\right)\mu_{t}\left(\mathrm{d}u\right).$

In particular, the optimal filter solution,
$\hat{u}_{t}=\mathbb{E}\left[u_{t}\mid\mathcal{G}_{t}\right]$, can be defined
based on $\mu_{t}$. It shows that $\hat{u}_{t}$ is the minimizer
$\mathbb{E}\left[\left|\hat{u}_{t}-u_{t}\right|^{2}\right]=\min_{v}\mathbb{E}\left[\left|v-u_{t}\right|^{2}\right]$
among all $v\in L^{2}\left(\Omega,\mathcal{G}_{t},\mathbb{P}\right)$ in the
set of $\mathcal{G}_{t}$-measurable square-integrable random variables for any
fixed $t$. The filtering equation for the conditional probability distribution
$\mu_{t}$ is verified to satisfy the Kushner-Stratonovich equation

$\mathrm{d}\mu_{t}\left(\varphi\right)=\mu_{t}\left(\mathcal{L}\varphi\right)\mathrm{d}t+\sigma\left(H,\varphi;\mu_{t}\right)\mathrm{d}\nu_{t}.$
(A2)

On the right hand side of the above equation, the first term is the drift due
to the infinitesimal generator
$\mathcal{L}=F\cdot\nabla+\frac{1}{2}\Sigma\Sigma^{\intercal}:\nabla\nabla$ of
the signal process (A1a); the second term represents the correction from the
observation process (A1b). The innovation process,
$\mathrm{d}\nu_{t}=\mathrm{d}y_{t}-\mu_{t}\left(H\right)\mathrm{d}t\in\mathbb{R}^{p}$,
is a $\mathcal{G}_{t}$-Brownian motion under the probability measure
$\mathbb{P}$, and $\sigma\left(H,\varphi;\mu_{t}\right)=\mu_{t}\left(\varphi
H^{\intercal}\right)-\mu_{t}\left(\varphi\right)\mu_{t}\left(H^{\intercal}\right)\in\mathbb{R}^{1\times
p}$ gives the coefficient with finite quadratic variation, where
$\mu_{t}\left(H\right)$ is the componentwise measure of the vector-valued
function $H$. In addition, the filtering equation (A2) is shown to have a
unique solution under proper conditions (Theorem 3.30 and 4.19 in [2] and
Theorem 7.7 in [56]) that is also stable (Theorem 2.7 in [6]). Therefore, this
guarantees that the solution of the Kushner-Stratonovich equation (A2)
uniquely characterizes the filter distribution $\mu_{t}$ as a
$\mathcal{P}\left(\mathbb{R}^{d}\right)$-valued stochastic process.

Next, assume that the conditional probability $\mu_{t}$ possesses a square
integrable density,
$\mu_{t}\left(\mathrm{d}x\right)=\rho_{t}\left(x\right)\mathrm{d}x$, with
respect to the Lebesgue measure. It can be shown under proper conditions
(Corollary 7.18 in [2] and [31]), the conditional probability solution
$\mu_{t}$ of (A1) has a probability density $\rho_{t}\in
W_{k}^{2}\left(\mathbb{R}^{d}\right)$. According to (A2) for the conditional
probability, $\rho_{t}$ can be found to be the unique solution to the
following SPDE

$\partial_{t}\rho_{t}=\mathcal{L}^{*}\rho_{t}\mathrm{d}t+\rho_{t}\left(H-\bar{H}_{t}\right)^{\intercal}\left(\mathrm{d}y_{t}-\bar{H}_{t}\mathrm{d}t\right),\quad\rho_{t=0}=\rho_{0},$
(A3)

where
$\bar{H}_{t}=\int_{\mathbb{R}^{d}}H\left(x\right)\rho_{t}\left(x\right)\mathrm{d}x$
and $\rho_{0}\in L^{2}\left(\mathbb{R}^{d}\right)$ is the absolute continuous