Forecast–assimilation (FA) process with unstable dynamics]
Asymptotic behavior of the forecast–assimilation process with unstable
Dan Crisan
Department of Mathematics, Imperial College London, 180 Queen's Gate, London SW7 2AZ, UK
Michael Ghil
Geosciences Department and Laboratoire de Météorologie Dynamique (CNRS and IPSL), École Normale Supérieure and PSL University, Paris, France
Atmospheric & Oceanic Sciences Department, University of California at Los Angeles, Los Angeles, CA, USA

Extensive numerical evidence shows that the assimilation of observations has
a stabilizing effect on unstable dynamics, in numerical weather prediction
and elsewhere. In this paper, we apply mathematically rigorous
methods to showing why this is so. Our stabilization results do not assume
a full set of observations and we provide examples where it suffices
to observe the model's unstable degrees of freedom.

This paper is dedicated to the memory of Anna Trevisan
and to her contributions to data assimilation

Department of Mathematics, Imperial College London, London, UK.

Geosciences Department and Laboratoire de Météorologie Dynamique (CNRS and IPSL),
Ecole Normale Supérieure and PSL University, F-75231 Paris Cedex 05, France
Department of Atmospheric and Oceanic Sciences, University of
California,Los Angeles, CA 90095-1565, USA

V. Bjerknes first described weather prediction as an initial-value problem
in 1904 <cit.>. As J. von Neumann and associates started using
computers to implement this idea immediately after World War II, it quickly
became apparent that the requisite initial data available were incomplete
<cit.>. The appearance of weather
satellites in the 1960s led further on to the concept of time-continuous
assimilation of remote-sensing data <cit.>.
Nowadays, data assimilation (DA) is being applied across all the areas of
the climate sciences and much beyond <cit.>.
Three crucial issues are still insufficiently well understood: (i) standard
proofs for the convergence of the DA process rely on the stability of the
model dynamics, even in the linear case, while atmospheric and oceanic
dynamics are famously unstable and chaotic <cit.>; (ii) data availability over time
appears to successfully compensate for insufficient instantaneous coverage
in space <cit.>; and, last
but not least, (iii) it appears that observations of a model's unstable
manifold are sufficient for the convergence of the time-continuos
forecast–assimilation cycle <cit.>.
The present paper uses concepts and methods from
the stochastic calculus <cit.>, random dynamical systems <cit.>
and nonlinear filtering <cit.>
to achieve several significant steps in
clarifying all three of these issues.

§ INTRODUCTION AND MOTIVATION

§.§ The forecast–assimilation cycle in meteorology

A key metaproblem of data assimilation (DA) in atmospheric, oceanic and
climate problems is to show that sequential filters of various degrees of
sophistication are stable, and that they converge to solutions with suitable
properties that lie sufficiently close to the observations, such as they
are, in the case of the unstable dynamics that characterizes these
problems. Moreover, a full solution to this metaproblem should allow one to
compare, with relative ease, the efficiency and accuracy of several filters.

Heuristically, the motivation for this metaproblem being soluble is the
success of practical DA methods in numerical weather prediction (NWP), and
in related oceanographic and climate problems, in keeping track of a
system's observed state <cit.>. We outline
herein some simple ideas of why sequential filters do have a chance of being
stable and convergent, even in the presence of dynamically unstable modes.

As the importance of DA methodology in the climate sciences and an
increasing number of other areas, all the way to finance, has been growing rapidly,
DA has attracted more and more attention in relevant areas of the mathematical
sciences <cit.>. As we shall argue
further below, this increase of attention has not covered as yet the issues of
instability of the basic dynamics that one wishes to track nor that of partial
observations. The main point of this paper is to take substantial steps in
addressing these two issues with the same degree of rigor as that used so
far in addressing the considerably simpler situation of stable dynamics
and complete measurements.

To lead up to the full complexity of the setting involved, consider first, for simplicity,
a scalar, linear DA problem in continuous time, using “sloppy” notation.
The model equation is

\begin{equation}
\dot x = Ax+u(t),  \label{model}
\end{equation}

and the observation equation is

\begin{equation}
\dot y = Hx+v(t),  \label{obs}
\end{equation}

where $(u,v)$ are noises that are white in time $t$ and have “nice” densities
— i.e., centered and with finite variance, e.g., Gaussian — while the
dot notation $(x, y)^\cdot$ stands for the time derivative. We keep caps $%
(A,H)$ for the model and observation operators because the scalar case is
supposed to be just shorthand for the vector–matrix case. The fully
nonlinear and high-dimensional problems arising in the actual applications
involve also much more complex noise processes, of course; see,
for instance, the work of C. Nicolis and coworkers
in ref. <cit.> and further references therein.

The forecast–assimilation (FA) process for $\hat{x}$, the best linear
unbiased estimate of $x$, obeys<cit.>

\begin{equation}  \label{FA}
\dot{\hat{x}} = A\hat{x} + K(y - H\hat{x}).
\end{equation}

Here $K$ is a weight matrix, which equals the Kalman-Bucy optimum <cit.> in the Gaussian-noise case, and $y - H\hat{x}$
is the innovation vector that equals the difference between the actually
observed value of $x$ and the one expected by forecasting this value based
on past observations.

The intuitive
motivation to hope for convergence of this FA cycle to the true evolution —
or for its synchronization with the observations<cit.> — is
simply re-writing equation (<ref>) as

\begin{equation}
\dot{\hat{x}}=(A-KH)\hat{x}+Ky.  \label{sync}
\end{equation}

This equation exhibits the new, and hopefully stabler, dynamics $(A-KH)$ of
the FA process vs. the original, pure-evolution dynamics $A$. It also suggests
using random dynamical system (RDS) theory<cit.> for the FA problem, given the time-dependent forcing
by $Ky$, in which the observations are subject to random errors.
An exhaustive presentation of RDS theory is given in the Ludwig
Arnold monograph <cit.>, which makes, however, somewhat
difficult reading for the non-specialist. More accessible presentations
for DA practitioners can be found in refs. <cit.>.

One has to show that the nonlinear, multidimensional — and possibly even
infinite-dimensional FA evolution, as in the generalization of Eq.
(<ref>) to partial differential equations (PDEs) — is stable <cit.> even in the presence of dynamic
For the way that it might still suffice, in the presence of dynamic
instabilities, to have dim$\{y\}\le $ dim$\{x\}$, see numerical results
for simplified atmospheric and oceanic models in ref. <cit.>
and in Section 2 of the Ghil (1997)<cit.> review paper, for instance.
Frank and Zhuk <cit.> did obtain
a mathematically rigorous result along these lines for a deterministic
system of nonlinear ordinary differential equations.

We restrict ourselves here to the finite-dimensional case: in the
operational practice of numerical weather prediction (NWP), the PDEs
governing atmospheric and oceanic flows are discretized in physical space
— using finite differences, finite elements or spectral and
pseudo-spectral methods <cit.>. These days,
the number of resulting finite-difference equations in time is very large
— up to order of $10^8$–$10^9$ — but still finite and will stay so for
the foreseeable future.

The nonlinearity of atmospheric and ocean dynamics <cit.> compels us to deal not just with the mean
and variance of the estimated state $\hat x$, as in the classical
Kalman-Bucy filter <cit.>, but with the
entire probability distribution function (pdf) of the state $x$, conditioned on the observations $z$.
This pdf may be multimodal or, more generally, non-Gaussian <cit.>
and include the presence of long tails due to extreme events <cit.>.

To fully describe the pdf of the observed state $x(t)$ given the data $z(t)$, one needs more than this pdf's mean and variance. Indeed, absent the linearity and Gaussianity assumptions, the system of equations satisfied by the mean and the variance of the pdf of the observed state is no longer closed and we need more that just these two quantities to describe the FA process and compute its evolution. This point of view is, by now, widely shared by the operational NWP community <cit.>.

In this paper, we assume that the model $x$ is a, possibly unstable, stochastic process
explicitly defined in section <ref>. We also assume that
the model is observed only partially and that the observations arrive
continuously in time. The latter assumption is consistent with the already
mentioned continuous flow of observations in the satellite era <cit.>.

In fact, shortly after the advent of meteorogical satellites in the late
1960s, Charney et al. <cit.> formulated the conjecture that
a knowledge of the continuous time history of the atmospheric temperature
field will allow one to determine the other state variables, in particular
the winds. M. Ghil and coauthors provided analytical arguments for the
correctness of this “Charney conjecture” in two-dimensional
(2-D) geophysical fluid dynamics (GFD) models <cit.> and
documented its usefulness numerically with time-continous DA of actual
remotely sensed temperatures in a fairly realistic NWP-type,
three-dimensional (3-D) atmospheric model <cit.>. E. Titi and
coauthors provided rigorous proofs in a purely deterministic setting for
both 2-D and 3-D models of GFD interest <cit.>.

Figure <ref> here illustrates the FA process's evolution in the
presence of observations that are partial in both their nature — i.e.,
temperature vs. winds — and their spatial coverage — ocean vs. land
in this simple example. The figure represents DA results using the Kalman-Bucy
(hereafter KF) for a linear, mid-latitude shallow-water model in one space
dimension with a simplified geometry <cit.>. In this geometry,
there are two data-rich regions of equal size that stand for the North
American and Eurasian land masses, alternating with two data-poor regions of
the same size that stand for the North Atlantic and North Pacific.

In the particular numerical experiment selected here for illustration purposes,
all three model variables — namely the geopotential height $\phi = gh$ of the
free surface, where $h$ is the actual height and $g$ the acceleration of
gravity, along with the cartesian velocity components $(u, v)$ — were
available at the so-called synoptic times of noon and midnight GMT over the
land areas, while no data at all were available over the ocean areas. In
this case, all three curves for $u (t), v(t)$ and $\phi (t)$ — as well as
for the energy $E = u^2 + v^2 + \phi^2/\Phi$, where $\Phi$ is the
equilibrium value of $\phi$ about which the equations are linearized —
have exactly the same behavior.

Over the data-rich land, the error drops sharply at the first observing time,
12 h after the start, and it grows parabolically between each observing time
and the next one, due to the KF's Gaussian and white-in-time model noise.
The parabolic growth of the root-mean-square (RMS) error in this case
is due to the additive model noise, as in the case of scalar Brownian motion
<cit.>, in which it is the variance of the process that is
proportional to time $t$; compare with the dashed red line in Fig. <ref>
below. Over land, though, the evolution of the RMS error
asymptotes very quickly, in 1–2 model days, to one in which the FA error
level oscillates around that of the observational noise, which equals
roughly 0.9 in nondimensional units.

Typical results of a forecast–assimilation (FA) cycle.
(a) Expected root-mean-square (RMS) error over land;
(b) expected RMS error over the ocean; and
(c) expected RMS error over the entire model domain.
This figure is based on Fig. 6 in Ghil et al. <cit.>,
with permission, and E. Bach provided the version herein.

Over the data-poor oceans, the RMS error still decreases, due to the
advection of information from the land areas by the mean westerly winds, $U
= - \partial \Phi/\partial y$. But this error decrease is slower than over
land, the asymptotically periodic behavior is only reached after 4–5
model days, and the mean values of RMS errors stay above the observational
noise. Finally, the RMS error behavior over the entire area is essentially a
weighted average of the results over land and over ocean.

This type of behavior was modified by the advent of time-continuous
satellite data over the oceans, as shown, for instance, by the work of Halem
et al.<cit.>, as part of the DA studies associated with the
Global Atmospheric Research Experiment (GARP). While the Ghil et al. <cit.> model was one-dimensional, linear, stable and had only a
rather small spatial resolution, that of Halem et al.<cit.>
and many others were fully 3-D, nonlinear, unstable, and had rather high
resolution by the standards prevailing at that point in time.

Figure 5 in Halem et al.<cit.> (not shown here) clearly
indicates the improvement in 6-hr forecasts over the Western U.S. from
initial states that do use the time-continuous satellite data over the North
Pacific vs. those that use only the conventional data available over land.
The mechanism that advects information by the westerly winds — or misinformation,
from ocean to land, in the case of the conventional observing network — is
clearly still working, in spite of the unstable and nonlinear dynamics of
the forecast model and of the ad-hoc sequential-estimation method in ref.
<cit.>. The latter was a successive-correction method
<cit.> in the Halem et al.<cit.> 3-D
model, while it was a KF in the Ghil et al. <cit.> 1-D model.

An idealized version of the error components that play a role in the
FA process are illustrated in Fig. <ref> and discussed further
in Appendix I. The major point
of the figure and of the associated appendix is that high-end, operational
NWP models have various sources of instability, and that NWP would be
impossible if DA did not stabilize the FA dynamics and thus achieve
real-time low-error forecasts. The main purpose of the present paper is to
justify rigorously the numerically observed fact that this is, indeed, the

Schematic diagram of forecast error growth in a numerical weather
prediction (NWP) model; see Appendix I for explanations. The perfect-model
(dashed blue) and real-model (solid blue) curves here should be compared
with the blue curve in the two panels of Fig. <ref>
there, for $S = 0$ and $S \neq 0$, respectively. The straight red line
(dash-dotted) labeled $Qt$ here represents the linear growth of forecast
error variance due to additive white noise, as is the case in Fig. <ref>(a),
between updates and over land,
while the Lorenz <cit.> and Dalcher and Kalnay <cit.> models in the appendix only assume a constant
deterministic model error $S$. E. Bach kindly provided this figure.

§.§ A rough sketch of the mathematical formulation

In the present paper's context, the “best estimate” of the model $x$ is
its conditional distribution with respect to the observational data $z(t)$
available up to the current time $t $. We will denote this probability
distribution by $\pi (t) \equiv \pi_t$ and the distribution of the model $x$
in the absence of any observational data by $p(t) \equiv p_t$. The pdf $p$
is called the prior distribution and $\pi$ the posterior distribution, where
we dropped for simplicity the dependence on time $t$. The two distributions
$\pi$ and $p$ can be viewed as dynamical systems that both evolve in the
infinite-dimensional space of probability measures $\mathcal{P}$ over
the model state space $\mathbb{R}^{d}$. Put in simple words, our results
refer to the ideal, truly optimal filter and not to any specific approximations
thereof, such as the extended Kalman filter (EKF <cit.>) or
ensemble Kalman filter (EnKF <cit.>) or any of the variational methods <cit.>.

For the comfort of the interested reader coming from the NWP community
or from other areas where DA is being used — or its use is being contemplated
— we are summarizing in Table  I  the correspondence of
key symbols and terms in this paper vs. their counterparts in NWP use.

Notation    Data Assimilation Language      Stochastic Filtering Language

$p_t$    Probability distribution of pure forecast     Prior distribution of the signal

$\pi_t$     Probability distribution of the analysis      Conditional distribution of the signal

$ \lambda_t^f$     Deterministic forecast model $M$ <cit.>      Push-forward operator

Correspondence of symbols and terms in this paper vs. their counterparts in NWP

Depending on the model for $x \in \mathbb{R}^{d}$, the prior system may be
unstable: Starting the prior system $p$ from two different initial conditions, the $p_0$
that results in $p_t$ and $\mu \neq p_0$, and using the
same forecast operator for both
produces two probability measures, $p_t$ and $p^\mu_{t}$, that will diverge from
each other in time, in the sense that the Wasserstein distance $W_2$ between the two
tends to infinity. The distance $W_2$ between two probability measures
$\mu$ and $\nu$ is given by

\begin{equation}\label{eq:W2}
{W_{2}(\mu ,\nu ):=\left( \inf_{\gamma \in \Gamma (\mu ,\nu )}\int_{{{%
\mathbb{R}}^{d}}\times {{\mathbb{R}}^{d}}}|x-y|^{2}\,\mathrm{d}\gamma
(x,y)\right) ^{1/2},}
\end{equation}

where ${\Gamma (\mu ,\nu )}$ denotes the collection of all measures on ${{%
\mathbb{R}}^{d}}\times {{\mathbb{R}}^{d}}$ with marginals $\mu $ and $\nu $ on the first and second factor, respectively.

This distance is often called in applications the
“earth mover's distance,” following the original motivation of G. Monge<cit.>.
Ghil <cit.> originally proposed the idea of using the Wasserstein distance
in the context of the climate sciences as a way to generalize the traditional concept
of equilibrium climate sensitivity <cit.> in the presence of
a time-dependent forcing, such as seasonal or anthropogenic forcing.
Further details are given in Appendix II.A and considerably more information on the definitions and methods used herein can be found in the Panaretos and Zemel monograph <cit.>.

As stated already several times, we concentrate on the prior process's $p_t$ being unstable,
i.e., starting it from the initial condition $\mu \neq p_o$ will lead to divergence of the
trajectory $p^\mu_{t}$ from $p_t$ in Wasserstein distance $W_2$. We show in
Sec. <ref> that, to the contrary, starting the posterior process $\pi^\mu_{t}$
from the initial condition $\mu \neq p_o$ and evolving it with the same FA operator
as that used to evolve $\pi_{t}$ will generate a probability measure $\pi^\mu_{t}$ that
will keep the  $W_2$ distance between $\pi^\mu_{t}$ and $\pi _{t}$ bounded in
expectation. Moreover, in the linear case, we show that  the $W_2$
distance between $\pi^\mu_{t}$ and $\pi _{t}$ actually tends to zero.

Mathematically, we consider the forecast operator $\lambda _{t}^{f }$ and
the FA operator $\lambda _{t}$ associated with the FA process, respectively,
which are defined by

\begin{align}
& \lambda _{t}^{f }: \mathcal{P}(\mathbb{R}^{d})\rightarrow \mathcal{P}(
\mathbb{R}^{d}),\ p_{t}=\lambda _{t}^{f }p_{0},  \label{eq:fcst} \\
& \lambda _{t}: \Omega \times \mathcal{P}(\mathbb{R}^{d})\rightarrow
\mathcal{P}(\mathbb{R}^{d}),~\pi _{t}\left( \omega \right) =\lambda
_{t}\left( \omega \right) \pi _{0},  \label{eq:FA}
\end{align}

and we show that, under certain conditions,

\begin{equation}
\sup_{t\in \lbrack 0,\infty )}\mathbb{E}\left[ W_2\left( \lambda
_{t}\mu,\lambda _{t}\pi_{0}\right) \right] <\infty ,  \label{eq:expl}
\end{equation}

whilst possibly $\lim_{t\rightarrow \infty }W_2( \lambda _{t}^f\mu,\lambda _{t}^fp_{0}) =\infty$.
Moreover, in the linear case, we will show that $\lim_{t\rightarrow \infty
}W_2\left( \lambda _{t}\mu,\lambda _{t}\pi_{0}\right) =0$.

Note that, in Eqs. (<ref>) above, both the push-forward operators $\lambda _{t}^{f }$ and $\lambda _{t}$ are stochastic processes and
act on the full probability measures $p$ and $\pi$, respectively. In the present  framework, we do
not limit ourselves just to the mean and variance of the state $x(t)$, as is
the case for the linear KF <cit.>.  The distance $W_2\left( \lambda
_{t}\mu,\lambda _{t}\pi_{0}\right)$ between the RDS $\lambda
_{t}\mu$ starting from $\mu$ and the RDS $\lambda
_{t}\pi_0$ starting from $\pi_0$  is of course random: It depends on the realization of the observation process,
and hence the expectation in Eq. (<ref>) is taken with respect to the pdf of  this process.

We will prove rigorously, in the precise sense described above, that the incorporation
of observational data into the FA process does indeed have a stabilizing
effect on unstable dynamics, as shown abundantly by NWP practice and
suggested, in particular, by the work of Anna Trevisan and her collaborators
on assimilation in the unstable manifold (AUS)<cit.>. The latter work was one source of
inspiration for the present paper.

Another source was RDS theory, according to the arguments presented by one
of the authors (MG) at the “Symposium Honoring the Legacy of Anna
Trevisan” held in Bologna, Italy in October 2017. A more specific source of
inspiration for the rigorous mathematics herein was the paper of L. Arnold
<cit.> on stabilization by noise, which was presented to MG by
Franco Flandoli during the trimester on “The Mathematics of Climate and the
Environment”, held at the Institut Henri Poincaré
in Paris in Fall 2019.

It turned out, however, rather quickly that the
latter paper's arguments could not be applied directly to the DA problem at
hand, since the deterministic component of the process under study there is
linear, and we did not see how to extend Arnold's arguments <cit.> to fully nonlinear processes.
Moreover, the latter arguments are only valid for finite-dimensional
systems, while we are dealing here with RDSs that evolve in the infinite-dimensional space of probability
measures $\mathcal{P}(\mathbb{R}^d)$.

The layout of the paper is dictated by the intent to bring the two
communities — of DA practitioners, on the one hand <cit.> —
and of the rapidly increasing numbers of applied mathematicians and
physicists interested in DA on the other <cit.> — closer together. Hence, after this fairly long
introduction, we describe in Sec. <ref> the precise
mathematical framework that is used herein. The main rigorous results are
outlined in Sec. <ref> and further details on definitions and proofs
appear in Appendix II. Conclusions and some thoughts on further work appear
in Sec. <ref>.

§ MATHEMATICAL FRAMEWORK

§.§ Prior results

The two stabilization theorems presented here, Theorems <ref> and <ref>, are related to
and do use in their proofs certain arguments from the existing
stability results in the nonlinear filtering literature, e.g., refs. <cit.>. Such results have not, by and large, matched the objectives of the applied DA community when studying the asymptotic behavior of the FA process.

To be more precise, the applied DA community is interested in results for a forecast cycle that is unstable — as is the
case in meteorology and oceanography — and for which applying DA has the mysterious but
salutary effect of stabilizing the FA process. In addition, the DA community has to rely,
typically, on results where only a small subset of the forecast cycle's degrees of freedom can be
observed. Without being comprehensive, of course, we give here a classification of the conditions under which the existing results on filter stability for nonlinear forecast processes hold.

* The forecast process is assumed to be ergodic or to have good mixing properties, e.g. Atar<cit.>, Atar and Zeitouni<cit.>, Budhiraja<cit.>, Chigansky and Liptser<cit.>,   Chigansky, Liptser and Van Handel<cit.>, Del Moral and Miclo<cit.>,  and Del Moral, Doucet and Singh<cit.>. Please note that ergodicity in these papers  is meant in the sense of convergence of the pdf of the push-forward process $\lambda^f$ as defined  in Eq. (<ref>), not in the pathwise sense of RDSs <cit.>. For example, the stochastic Lorenz
model<cit.> is ergodic in the pathwise sense of RDSs, but does not have a stable push-forward process $\lambda^f$. Note that ergodicity of the pdf of the push-forward process $\lambda^f$ does imply its stability. For further results on ergodicity of the push-forward process, see, for instance, Bakry, Cattiaux and Guillin <cit.>.

* The forecast cycle is fully observed, e.g., Picard<cit.>,
Crisan and Heine<cit.>, Stannat<cit.>,
Chigansky, Liptser and Van Handel<cit.>, and Van Handel<cit.>.

* The forecast process has a drift term which is the gradient of a convex function, or a perturbation thereof,  and  its noise term is strictly elliptic, e.g., Stannat<cit.>. This is a very restrictive class of forecast processes, which  are expected to be  stable; see Bakry & all<cit.>. Specific drifts belonging to this class are used in the classical Metropolis-adjusted Langevin algorithm known to converge faster to the diffusion's invariant measure; e.g., Roberts and Tweedie<cit.>.

To try to answer the question raised by the applied community, we no longer insist on proving
that the FA process is (exponentially) stable in the sense advocated by the theoretical
community. We relax the definition and only require that the FA process initialised from the
wrong distribution does not diverge too strongly from the correctly initialized FA process, even
when the forecast process does so; this is why the ergodicity assumption for the forecast process is
not useful.

We substantially strengthen, however, the stabilization result, in the sense that we
want to control the mean and the second moment of the FA process. Again this is needed for
practical reasons. The practitioners want to know that the pointwise estimate of their algorithm
of choice does not diverge from the theoretical mean of the FA process. But they also want to
know that their error bars are not too different from the theoretical ones.

Furthermore, the results presented herein have the advantage that the signal $X(t)$ is not required to be fully observable. In particular, the dimension $n$ of the observation $Y(t)$ and the dimension $d$ of the signal $X(t)$, in the notation of Sec. <ref>, do not need to coincide; in operational NWP and many other applications, $d \ll n$ <cit.>.

Heuristically we need to be able to observe all the "unstable" directions, as suggested by A. Trevisan and coworkers
<cit.>, who gave several fairly realistic examples of this idea working quite well. The connection with our results is provided by some simple illustrative examples in  Appendix II.D.

In the linear case, the stability of the FA process is better understood;
see, for instance Ocone and Pardoux<cit.>,  Picard<cit.>, Van Handel<cit.>, and
The stability of the associated Riccati equation has been studied in Bishop and Del Moral<cit.>. The stabilization result obtained in Theorem <ref> below for the linear signal is weaker than many of the  existing results discussed above. For example, no (exponential) rates of convergence are deduced herein. This situation clearly leaves considerable room for proving stronger results, given the hypotheses of Theorem <ref>.

§.§ The present setting

Having presented in the previous subsection a quick review of previous mathematical results
on the stability of the FA process, we proceed now by introducing the setting of our two theorems
in the next section.
Let $X=(X^{i})_{i=1}^{d}$ be the solution of the following stochastic
differential equation driven by a $p$-dimensional Brownian motion process $%

\begin{equation}
X_{t}=X_{0}+\int_{0}^{t}f\left( X_{s}\right) \mathrm{d}s+\int_{0}^{t}\sigma
\left( X_{s}\right) \mathrm{d}V_{s}.  \label{signal}
\end{equation}

Here, we assume that $f=(f^{i})_{i=1}^{d}:\mathbb{R}%
^{d}\rightarrow \mathbb{R}^{d}$ and $\sigma =(\sigma ^{ij})_{i=1,\ldots
,d,j=1,\ldots ,p}:\mathbb{R}^{d}\rightarrow \mathbb{R}^{d\times p}$ are
globally Lipschitz. This will ensure that the equation (<ref>) has a unique solution.

As stated in Sec. <ref>, we are interested in tracking the evolution of
the full pdf of the prior and posterior processes, namely $p_t$ and $\pi_t$, respectively.
To do so, we recall that the process $X$ is a diffusion process with infinitesimal generator given by

\begin{equation*}
A\varphi =\sum_{i,j}\frac{1}{2}a_{ij}\partial _{i}\partial _{j}\varphi
+\sum_{i}f_{i}\partial _{i}\varphi ,
\end{equation*}

where $a_{ij}=\sum_{k}\sigma _{ik}\sigma _{jk}$. The prior distribution of $X_{t}$
is also called its law in the context of filtering <cit.>.

For an arbitrary measurable function $\varphi :\mathbb{R}^{d}\rightarrow \mathbb{R}$
that is integrable with respect to the law of $X_{t}$, one has

\begin{equation*}
p_{t}(\varphi )=\mathbb{E}\left[ \varphi \left( X_{t}\right) \right] .
\end{equation*}

By restricting $\varphi $ further to lie in a suitably chosen space of functions denoted by $%
\mathcal{D}(A)$,

\begin{equation}
p_{t}(\varphi )=p_{0}(\varphi )+\int_{0}^{t}p_{s}(A\varphi )\mathrm{d}s.
\label{priordistributionequation}
\end{equation}

Let $\lambda _{t}^{f}:\mathcal{P}(\mathbb{R}^{d})\rightarrow \mathcal{P}(%
\mathbb{R}^{d})$ be the push-forward operator associated with equation
(<ref>). In particular, for arbitrary $\mu \in
\mathcal{P}(\mathbb{R}^{d})$, the map $t\rightarrow \lambda _{t}^{f}\mu
=:p_{t}^{\mu }$ is the solution of (<ref>) with
the initial condition $\mu $.

Let $Y$ be an $n$-dimensional process, where $n$ is the number of
observations that are taken to be one-dimensional and have measurement
noises modelled by independent Brownian motions:

\begin{equation}
s+W_{t}^{i},~~~i=1,...,n,  \label{observation}
\end{equation}

and $h_{i}\left( x\right) $ are the corresponding observation operators,
which generalize the observation matrix from the linear framework; see

Appendix II.C for details.

Let $\pi =\{\pi _{t},t\geq 0\}$ be probability measure–valued process that
gives us, at time $t\geq 0$, the conditional distribution of the signal $%
X_{t}$ given the observations accumulated up to time $t$, $\{Y_{s},\ s\in
\lbrack 0,t]\}$. It is this process that we call the FA process in our rigorous mathematical
setting. In other words, $\pi _{t}$ satisfies

\begin{equation*}
\pi _{t}(\varphi )=\mathbb{E}\left[ \varphi \left( X_{t}\right) |Y_{s},\
s\in \lbrack 0,t]\right] ,
\end{equation*}

where $\varphi :\mathbb{R}^{d}\rightarrow \mathbb{R}$ is an arbitrary
measurable function that is integrable with respect to the law of $\pi _{t}$. The
FA process satisfies the following stochastic partial
differential equation, formulated here in the integral form

\begin{eqnarray}
\pi _{t}(\varphi ) &=& \ \pi _{0}(\varphi )+\int_{0}^{t}\pi _{s}(A\varphi )\,%
\mathrm{d}s+\int_{0}^{t}\left( \pi _{s}(\varphi h^{\top })-\pi _{s}(h^{\top
})\pi _{s}(\varphi )\right) (\mathrm{d}Y_{s}-\pi _{s}(h)\,\mathrm{d}s), \label{ks0} \\
&=& \ \pi _{0}(\varphi )+\int_{0}^{t}\pi _{s}(A\varphi )\,\mathrm{d}%
s+\int_{0}^{t}\left( \pi _{s}(\varphi h^{\top })-\pi _{s}(h^{\top })\pi
_{s}(\varphi )\right) \mathrm{d}I_{s}, \label{ks1}
\end{eqnarray}

for any test function $\varphi \in \mathcal{D}(A)$. Here, $I$ is
the innovation process, defined as

\begin{equation*}
I_{t}^{i}=Y_{t}^{i}-\int_{0}^{t}\pi _{s}\left( h_{i}\right) \mathrm{d}%
s,~~~i=1,...,n,~~~t\geq 0,
\end{equation*}

which is the rigorous analog of the innovation vector in Eq. (<ref>)
of Sec. <ref>. The innovation process is a Brownian motion, see e.g., Chapter 3 in Bain&Crisan<cit.>. In particular it is a martingale and stochastic integrals with respect to the  innovation process are easier to manipulate. In particular upper bounds for the stochastic integrals with respect to martingales are easier to obtain than those  for  stochastic integrals with respect to general  semi-martingales. The observation process is a semi-martingale (it is a Brownian motion plus a drift term) so harder to handle.

More on stochastic partial differential equations and their difficulties can be found, for instance, in Liu and Röckner<cit.>, or Rozovsky and Lototsky<cit.>.

In order to analyze the asymptotic behavior of the FA process, we can recast
the solution of  (<ref>) as an RDS. More
precisely, there exists a measurable map $\lambda :[0,\infty )\ \times
\mathcal{P}(\mathbb{R}^{d})\times \Omega\rightarrow \mathcal{P}(\mathbb{R}%
^{d})$, $(t,\mu ,\omega )\mapsto \lambda(t,\omega )\mu$ such that $\lambda
(0,\omega )=I$, namely the identity map on $\mathcal{P}(\mathbb{R}^{d})$, and

\begin{equation}
\lambda (t+s,\omega )=\lambda \left( t,\vartheta _{s}\omega \right) \circ
\lambda (s,\omega )  \label{cocycle}
\end{equation}

for all $t,s\in \lbrack 0,\infty )$ and for all $\omega \in \Omega $. In (
<ref>) the symbol $\circ $ means map composition. $A$ family of maps
$\lambda (t,\omega )$ satisfying (<ref>) is called a cocycle, and (
<ref>) is the cocycle property. The map $\left\{ \vartheta
_{t}:\Omega \rightarrow \Omega \right\} ,t\in \lbrack 0,\infty )$ is a
family of measure-preserving transformations of a probability space $(\Omega
,\mathcal{F},P)$ termed the shift operators;  see, for
instance, Section 2.5 in Karatzas and Shreve <cit.> for
further details on the shift operators.

Using this map, the solution of (<ref>) can be expressed as

\begin{equation}\label{eq:rds}
\pi _{t}(\omega )=\lambda (t,\omega )\pi ^{0}.
\end{equation}

Moreover for an arbitrary $\mu \in \mathcal{P}(\mathbb{R}^{d})$, the process
$\pi^\mu =\{\pi^\mu _{t},t\geq 0\}$, defined as

\begin{equation*}
\pi^{\mu}_t(\omega):=\lambda(t, \omega) \mu
\end{equation*}

is the solution of the SPDE (<ref>) with initial condition $%
\mu$. Finally the map $\lambda(t, \omega)$ is a continuous map when we endow
$\mathcal{P}(\mathbb{R}^{d})$, (or, rather, the set of probability measures
with  second moment) with the topology induced by the Wasserstein

Since the FA process is infinite-dimensional, as explained in Sec. <ref>,
its RDS characterization is not immediate. RDS theory is well developed
for finite-dimensional processes <cit.>. A
subclass of these questions is settled in a fairly satisfactory manner by
the theory of stochastic flows; see, for instance, Ikeda and Watanabe <cit.> and Kunita <cit.>. A stochastic flow needs jointly
continuous dependence of the solutions of the stochastic differential equation
under consideration on time and on the initial state, except for a set of measure zero.
This often does not hold for infinite dimensions. Some infinite-dimensional systems do generate a
stochastic flow, others do not. For further definitions of possibly
infinite-dimensional RDSs, as well as for other related
results, we refer to Crauel and Flandoli<cit.> and Flandoli <cit.>. The RDS characterization of the FA process is discussed in
Appendix II.E.

In this paper, we show that, despite the possible divergence of the prior
distributions, that is,

\begin{equation*}
\lim_{t\rightarrow \infty }W_{2}(p^{\mu}_{t},p_{t})=\infty ,
\end{equation*}

the FA process has a stabilizing effect, in the sense
that it keeps the distance $W_{2}(\pi^{\mu}_{t},\pi _{t})$ uniformly bounded in
expectation. Moreover, in the linear case, it makes the distance $%
W_{2}(\pi^{\mu}_{t},\pi _{t})$ vanish asymptotically. The main results of the
paper are Theorems <ref> and <ref> below.

§ MAIN RESULTS

We introduce now the Wasserstein metric on $\mathcal{P}^{2}(\mathbb{R}^{d})$
defined by (<ref>) on
the set of all probability measures on the collection of Borel sets $\mathcal{B}(\mathbb{R}^{d}) $
that have a finite second moment. Recall that

the set ${\Gamma (\mu ,\nu )}$ in (<ref>) denotes the collection of all measures on ${{%
\mathbb{R}}^{d}}\times {{\mathbb{R}}^{d}}$ with marginals $\mu $ and $\nu $
on the first and second factor, respectively; it
is called the set of all couplings of the measures $\mu $ and $\nu $.

The Wasserstein metric is equivalently defined by

\begin{equation}  \label{equiv}
{W_{2}(\mu ,\nu )=\left( \inf \mathbb{E}{\big [}|X-Y|^{2}{\big ]}\right)
\end{equation}

where $\mathbb{E}{[Z]}$ denotes the expected value of a random
variable or vector $Z$ and the infimum is taken over all joint distributions of
the random variables $X$ and $Y$ with marginals $\mu $ and $\nu $

The main results of the paper are:

For nonlinear coefficients $f$, $h$ and $\sigma$ and measures $\pi_0$ and
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ that satisfy the conditions stated in  Appendix II, Sec. B,
there exists a bound $R=R\left( \pi_0,\mu\right)$ such that

\begin{equation} \label{eq:Thm1}
\sup_{t\ge 0}\mathbb E[W_{2}(\pi _{t}^{\mu },\pi_{t})]\leq R.
\end{equation}

The complete proof of this theorem is given in Appendix II, Sec. B. For the benefit of the curious but hasty reader, we provide here a brief sketch of the argument. First, we give a bound on the difference $%
\hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}$ between the mean of the FA process initialized from $\mu $ and the original
FA process, initialized from $\pi _{0}$. This is done in two steps: For arbitrary $\delta>0$,
we deduce that there exists a constant $c_{\delta }$ independent of $k$ such that

\begin{equation*}
\sup_{t\in \lbrack k\delta ,\left( k+1\right) \delta ]}\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \leq
c_{\delta }\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}%
_{k\delta }\right\vert \right] . %\label{cd}
\end{equation*}

This inequality appears as Eq. (<ref>) in the appendix.

Next, we show that there exists $R_{\delta }$ such that $\sup_{k\ge 0}\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}%
_{k\delta }\right\vert \right] \leq R_{\delta }.  $ These two inequalities give us a uniform bound, over $t\in [0,\infty)$, of the difference $\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right]$. Finally the uniform bound, for all positive times, of $\mathbb{E}[W_{2}(\pi _{t}^{\mu },\pi_{t})] $ comes by means of Lemma <ref> in the appendix
from the bound on the difference $\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right]$ and that of the  covariance matrices of the
FA process initialized from $\mu $ and, respectively, the original
FA process, initialized from $\pi _{0}$.

For linear coefficients $f$, $h$ and $\sigma$ and measures $\pi_0$ and
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{d})$ that satisfy the conditions stated in Appendix II, Sec. C,
we have the much stronger result that

\begin{equation} \label{eq:Thm2}
\lim_{t\rightarrow \infty }
W_{2}(\pi _{t}^{\mu },\pi_{t})=0.
\end{equation}

The complete proof of Theorems <ref> is
given in Appendix II, C. Again, we provide here a brief sketch of the argument.
We define $\pi _{t}^{0}:=N\left( \hat{x}_{t}^{0},P_{t}^{0}\right)$
to be a suitably chosen probability measure process that serves as a "reference."
Then we show that $\pi ^{\mu }$ gets asymptotically close to the
reference process $\pi ^{0}$, regardless of the initial condition, and
since this holds true for $\mu =\pi _{0}$ too, we immediately deduce that $\lim_{t\rightarrow \infty }W_{2}(\pi _{t}^{\mu },\pi _{t})=0$.

The reference process $\pi ^{0}$ is convenient to work with: its centered version  $N\left( 0,P_{t}^{0}\right) $ converges weakly, as well as in Wasserstein distance, to $\pi ^{\infty }=N\left( 0,P_{\infty }\right) $.
Using the equivalent definition of the Wasserstein distance
(<ref>), we deduce that $\lim_{t\rightarrow \infty}W_{2}(\pi _{t}^{\mu },\pi _{t}^{0})=0$,
if and only if the following three properties hold true:

* $\lim_{t\rightarrow \infty }| \hat{\pi}_{t}^{\mu }-\hat{x}%
_{t}^{0}| =0;$

* $\lim_{t\rightarrow \infty }\left\vert P_{\pi _{t}^{\mu
}}-P_{t}^{0}\right\vert =0;$ and

* $\lim_{t\rightarrow \infty }\left\vert \pi _{t}^{\mu }\left( \varphi
_{t}\right) -\pi _{t}^{0}\left( \varphi _{t}\right) \right\vert =0$ for any
bounded uniformly continuous function $\varphi $, where $\varphi
_{t}$ is the same function shifted by the mean $\hat{x}_{t}^{0},$ that is $%
\varphi _{t}\left( x\right) :=\varphi _{t}\left( x+\hat{x}_{t}^{0}\right) ,$
$x\in \mathbb{R}^{d}$.

These three properties are then shown to hold, thus completing the proof.

In the linear case of Theorem <ref>, it is fairly easy to verify that the Kalman-Bucy filter <cit.> satisfies the assumptions of the theorem. For the nonlinear case of Theorem <ref>, it is the subject of future research to find numerical criteria that guarantee the required assumptions.   This is a challenging  problem as there are many suboptimal approximations of an optimal filter and the verification of the corresponding hypotheses may prove more difficult.

Stability properties of suboptimal filters — e.g., the extended Kalman-Bucy filters (EKFs) <cit.>, as opposed to the truly optimal filter studied herein  — have also been studied under the assumption of uniformly stable and fully observable signals. The stability constraint for EKFs has been removed in <cit.>; see also <cit.> for a study of the stability of the mean-squared filtering error.

§ CONCLUSIONS AND FURTHER WORK

§.§ Summary

The main results of this work are given by Theorems 1 and 2 in Sec. <ref>. Essentially,

* For nonlinear dynamics or observations — including unstable dynamics of the prior process and subject to certain  technical but plausible assumptions — the supremum of the  expectation of the  Wasserstein distance $W_2$ between the true posterior solution and a solution of the FA process with the wrong initial conditions remains bounded at all future times; and

* For linear dynamics and observations — including unstable dynamics of the prior process and subject to certain technical but plausible assumptions — the Wasserstein distance $W_2$ between the two posterior distributions tends to zero.

§.§ Discussion

These results, to the best of our knowledge, are the first to address the stability of the posterior FA process given an unstable prior process.  In the linear case, convergence in $W_2$ of the posterior processes starting from correct and incorrect initial data, $p_0$ and $\mu$, has been demonstrated (Theorem <ref>) and the applicability to partial observations has been illustrated in Appendix II.D. In particular, observing the unstable components of the prior process, as originally proposed by A. Trevisan and her collaborators <cit.>, seems to be an excellent idea. The results in the nonlinear case only guarantee $W_2$-boundedness of the difference between the  two posterior processes starting from distinct initial probability measures. A considerable amount of practical DA work also indicates that the FA process can track the correct solution <cit.>, in particular when using observations from the unstable subspace <cit.>.

The first results herein towards a more realistic mathematical treatment of the unstable-dynamics case open the door to a whole slew of additional results, both theoretical and practical.

Deterministic and stochastic EnKF. Operational DA in NWP relies these days more and more on the EnKF <cit.>. But, in practice, most operational EnKF algorithms  randomize only over the NWP model's initial states and not over observations, too, as done herein.

The operational NWP literature on DA distinguishes, in fact, between the deterministic EnKF, which only takes into account random errors in the observations via the covariance matrix $\mathbf R$ of observational errors, and the stochastic EnKF, which explicitly simulates random errors in the observations. An excellent review of the EnKF for atmospheric DA in general appears in Houtekamer and Zhang <cit.>, with particular attention to this issue in its Section 2b. Lawson and Hansen <cit.> give interesting examples of the two versions of EnKF being applied to relatively simple examples of atmospheric and oceanic flows in one and two spatial dimensions, and Hoteit et al. <cit.> discuss some of the problems that might arise in the stochastic EnKF by introducing these observational random errors into the FA process. P.J. van Leeuwen <cit.> has proposed recently a self-consistent way of applying the stochastic EnKF.

Given the novel convergence results obtained herein in the presence of a random observation process and some of the renascent interest in the NWP literature, it might be worthwhile revisiting the usefulness of the stochastic EnKF. In particular, retaining random perturbations in the observations might obviate the need for artificial inflation of the ensemble's rapidly lost variance for the deterministic EnKF.

Multiple models and model error. In practice, in NWP and elsewhere, prediction can be served by more than one model. The models can differ by their spatial resolution, by the physical processes taken into consideration and by the numerical discretization of the PDEs governing them. Multi-model DA is discussed in some detail by Bach and Ghil <cit.>, including the issue of model error growth in this situation. It would be of considerable interest to extend the rigorous results herein to such a broader setting.

Typically, given the fact that DA is more expensive than straight forecasting <cit.>, it is natural to use lower-resolution models for the FA process than for the forecasting. In the set-up of Sec. <ref> herein, doing so corresponds to distinct prior processes $p_t$ and posterior processes $\pi_t$ and appropriate consideration of such issues would be quite worthwhile.

Parameter estimation and the synchronization point of view. In Sec. <ref>, we have mentioned already the view of the FA process as the synchronization of the forecast model with the observed process <cit.>. This point of view has been used as a unifying principle between DA and supermodelling, namely the use of ensembles of models that do not only serve for a posteriori averaging of their results but learn from each other in the process of a simulation or prediction run <cit.>. This learning is clearly related to the estimation of imperfectly known model parameters <cit.>.

In this broader perspective, one could try to demonstrate, given suitable hypotheses and observations, (i) the convergence of a single model's parameter estimation process; and (ii) the convergence of a supermodel to the observed process.

Practical examples. We presented in Appendix II.D,  a simple linear model and two observation schemes to illustrate the fact that stabilization by the FA process does not require observing all of a model's degrees of freedom and that observing just the unstable ones suffices. In future work, we aim to apply these results to realistic models, for example models that numerically approximate the PDEs of geophysical fluid dynamics.

Particle filters. The results presented in this paper are theoretical in nature. In practice, the FA process cannot be computed exactly: numerical approximations are required to estimate the posterior distribution of the signal given the data. Among these numerical approximations, particle filters have the crucial property of being theoretically justified in the sense that the numerical error can be controlled by the computational effort. Moreover, they are asymptotically consistent, i.e., as the number of data points used increases, the sequence of estimates converges in probability

A particle filter is a sequential Monte Carlo method in which the posterior distributiuon is approximated using a set of particles, yielding a random measures of the form
$\displaystyle\sum_{\ell} \mathrm{w}_{t}^{\ell}\delta({x_{t}^{\ell}})$,
where $\delta$ is the Dirac delta function, $\mathrm{w}_{t}^1, \mathrm{w}_{t}^2, \ldots $ are the weights of the  particles and $x_{t}^1, x_{t}^2, \ldots$ are their corresponding positions,<cit.>
centered around the state vector $x_{t}$. The approximations evolve in time, by following the time evolution induced by the prior model, and are corrected by the observations to keep them close to the evolution of the FA process <cit.>.

Particle filters have been very successful in many applications, including engineering, economics and finance; see, for instance, Doucet et al. <cit.> and the references therein. In recent years, applications of particle filtering to DA problems for planetary flows have flourished. For in-depth reviews of the most recent efforts in this direction, see refs. <cit.>. Such applications require enhancements of the classical particle filters in order to eable them to tackle the so-called curse of dimensionality, e.g., by relying on optimal transport ideas <cit.>, tempering <cit.>, localization <cit.>, model reduction <cit.>,  jittering<cit.>, nudging <cit.>, and judicious proposal densities <cit.>. Some of these approaches have been tested in operational NWP weather prediction systems <cit.>.  The suitability of particle  filters for high-dimensional problems has been studied in <cit.> and tested in <cit.>. For example, in <cit.>, the method is used for the stochastic incompressible two-dimensional Euler model with forcing and damping, while in <cit.> it is tested for a two-dimensional quasi-geostrophic model.

The theoretical results herein pave the way for the stability analysis of particle
filters under the same assumptions. More precisely, one can attempt to show that particle
filters have numerical errors that can be controlled uniformly in time. Again, we will be guided in  pursuing such results by existing ones in nonlinear filtering, stochastic analysis and applied probability <cit.>.    The bound will be in expectation, as in Theorem <ref> of this paper. Coupled with Remark <ref> in Appendix II, Sec. B, such results will offer theoretical validation to applying particle filters for long-run DA problems.

§ APPENDIX I. FORECAST ERROR GROWTH IN NWP

In order to better understand the nature and role of forecast error growth
in the FA process, we consider here three different models for error growth in
NWP, namely those of C. E. Leith <cit.>, E. N. Lorenz <cit.> and of A. Dalcher and E. Kalnay <cit.>.

The forecast error model of Leith<cit.> is

\begin{equation}  \label{eq:Leith}
\dot V = \alpha V + S,
\end{equation}

where $V$ is the mean-square error, $S$ is the systematic model error, and $t
$ is the lead time. The $V$ in this appendix should not be confused with the
Brownian motion process in Eq. (<ref>) of Sec. <ref>.
The forecast error growth in the Leith model is given by

\begin{equation}
V(t) = \left(V_0 + \frac{S}{\alpha}\right) e^{\alpha t} - \frac{S}{\alpha},
\end{equation}

where $V_0 = V(0)$ is the initial error. Note that short-time forecast
errors grow exponentially and that the systematic model error acts to
increase the coefficient of this growth.

Leith's forecast error model can only apply for short-time error growth,
since it does not saturate. Note that under certain statistical assumptions,
the mean-square error saturation value of a single forecast will be 2$C$,
where $C$ is the climatological variance. In a real NWP model with $N$
variables, the scalar $C$ will be equal to the mean trace of the
climatological covariance matrix $\mathbf{C}$, where $\mathbf{C}$ has
dimension $N \times N$. For ensemble forecasts, the saturation value becomes
$(1 + 1/m) C$, where $m$ is the ensemble size, see Leith1974<cit.>

Lorenz's model of forecast error growth <cit.> is

\begin{equation}  \label{eq:Lorenz}
\dot E = aE(E_\infty - E),
\end{equation}

where $E$ is the root-mean-square error and $E_\infty$ is its saturation
value. To compare this model directly to models based on mean-square error,
like Eq. (<ref>), we can change variables to $V = E^2$, and get

\begin{equation}  \label{eq:Lorenz_V}
\dot V = 2a V^{1/2}_{\infty} V \left(1 - (V/V_{\infty})^{1/2} \right).
\end{equation}

Lorenz's model includes a nonlinear saturation term, but does not
incorporate systematic model error $S$.

The error model proposed by Dalcher and Kalnay <cit.>
(henceforth DK) combines the key features of the Leith <cit.> and
Lorenz <cit.> models,

\begin{equation}  \label{eq:DK}
\dot V = (\alpha V + S)(1 - V/V_\infty);
\end{equation}

it thus includes both saturation $V_\infty$ and systematic model error $S$.
For short-time error growth, we can take $V_\infty\to \infty$, recovering
Leith's model. For $S = 0$, the model is similar to that of Lorenz, but with
$V_{\infty}$ having unit power in the saturation term, rather than $1/2$.

To compare the three models graphically, we set $\alpha = 1$, $V_0 = 1$, and
$V_\infty = 100$. To match the Lorenz model's short-term error growth to
that of the other two models, we set $a = \alpha/\left(2

Comparison between the model error growth $V(t)$ in the Leith
<cit.>, Lorenz <cit.> and DK
<cit.> models. The systematic model error is (a)
$S=0$, dashed curves; and (b) $S = 6$, solid curves. See legend for color

Figure <ref> shows a comparison between the three error
models, both with a perfect model for which $S=0$ and with an imperfect
model with $S=6$. The Leith curve is the same in both cases, since it does
not account for $S$. For the perfect model case, all three curves experience
similar short-term exponential error growth. However, the Leith and DK
models grow faster in the imperfect case than in the perfect case. In both
cases, the Leith model diverges from the other two curves in the medium
range, due to its lack of saturation. The Lorenz and DK curves both saturate
to $V_\infty$, although DK saturates more quickly due to the difference in
the functional forms of the saturation terms in the two models.

DK actually used their model's three parameters, $(\alpha, S, V_\infty)$, to
match their error growth curve to the archived real-time performance for the
years 1980-1981 of the European Centre for Medium-range Weather Forecasts
(ECMWF) model. Their Figure 9 shows separate, near-perfect fits out to 10
days, for boreal winter, when the weather is more active, and boreal summer,
when it is less so. Stroe and Royer <cit.> subsequently
generalized the DK model, by introducing the power $V^p$, with $p \neq 2$,
in the saturation term, and took the limit $p \to \infty$ to obtain

\begin{equation}  \label{eq:Royer}
\dot V = -aV \log(V/V_\infty).
\end{equation}

These authors found that Eq. (<ref>) gave better fits for
extended-range, 45-day experimental weather forecasts than either Eq.
(<ref>) or (<ref>).

Simmons et al. <cit.> also obtained rather good fits to
operational NWP model performance with Lorenz's quadratic error growth
model, but with a smaller error growth exponent that in the DK paper <cit.>. Trevisan et al. <cit.>, though,
showed — by using an intermediate, quasigeostrophic two-layer model on the
$\beta$-plane <cit.> — that only very small initial
errors in such a model obey Lorenz's quadratic error growth model and that
the error growth curve in general depends significantly on the magnitude of
the initial errors.

Savijärvi <cit.> combined features of the Lorenz and DK
models in the study of the (then) U.S. National Meteorological Center's
(NMC's) Medium-Range Forecast (MRF) Model's 0–10-day forecasts for
1988–1993. Growth parameters, as well as model and analysis errors for this
data set, were estimated using the quadratic error growth assumption.
Savijärvi showed that both the MRF model error and analysis error nearly
halved during the six years under study but, at the same time, the growth
parameters nearly doubled, since smaller errors grow faster.

Model error growth is thus a complex topic with much more to be said about;
see, for instance, the line of inquiry developed by C. Nicolis and coworkers, which
includes transient bimodality of the error's pdf <cit.>.
The topic's quick review in this appendix suffices, though, to
show the presence of error growth–generating instabilities in high-end,
operational NWP models. The ground covered here in Secs. <ref>
sec:framework and <ref>sec:wd and in Appendix II shows that DA can
overcome these instabilities,  in theory as well as in practice.

§ APPENDIX II. RIGOROUS DEFINITIONS AND PROOFS

§.§ A. The Wasserstein topology

In this appendix, we present the reason for showing the stabilizing effect of the FA process on unstable dynamics with respect to the Wasserstein topology — i.e., the topology generated by
the Wasserstein distance — and not with respect to the more popular weak
topology. Before doing so, we provide a little more information on the
Wasserstein distance for the benefit of NWP practitioners who might not yet
be familiar with it.

Gaspard Monge, an artillery officer in Napoleon's armies, as well as one of the founders of France's Ecole Polytechnique, introduced it as early as the 1780s <cit.>, and Leonid V. Kantorovich<cit.> used it during World War II in optimizing the transport of resources within the Soviet Union. The contemporary developments of this distance and of its applications are largely due to R. L. Dobrushin <cit.>, who coined the name  Wasserstein <cit.> distance for it, and to C. Villani <cit.>.

In the climate sciences, Ghil <cit.> illustrated the use of the Wasserstein distance for measuring the parameter sensitivity of simple models with time-dependent forcing, thus providing a link between nonautonomous dynamical systems theory <cit.> and optimal transport <cit.>. Robin et al. <cit.> then used this distance to compute the difference between the snapshot attractors of the Lorenz <cit.> model for different time-dependent forcings, while Vissio et al. <cit.> used it to help intercomparing climate models and evaluating their performance against given benchmarks in the Coupled Model Intercomparison Project that is part of the  Intergovernmental Panel on Climate Change process.

Returning to the main purpose of this appendix, let us introduce first some notation. For a measure $%
\mu\in \mathcal{P}^2({\mathbb{R}^d})$, we use the following notation:

* means: $\hat{\mu}=(\hat{\mu}^{i})_{i\in 1,\dots d}$ is the mean vector of $\mu$, i.e.,

\begin{equation} \label{eq:means}
\hat{\mu^{i}}=\int_{\mathbb{R}^d}x_{i}\mu(\mathrm{d }x), \quad {i\in
1,\dots d}; \quad|\hat{\mu}| = \left(\sum_{i=1}^d(\hat{\mu^{i}})^2\right)^{1/2};
\end{equation}

* second moments: $\mu^2$ is the sum of the second moments of $\mu$,  i.e.,

\begin{equation} \label{eq:moments}
\mu^{2}=\sum_{i=1}^d\int_{\mathbb{R}^d}x_{i}^2\mu(\mathrm{d }x);
\end{equation}

* covariance matrix: $P_\mu=(P_\mu^{ij})_{i,j\in 1,\dots d}$ is the covariance matrix of $%
\mu$, i.e.,

\begin{align}
& P_\mu^{ij} = \int_{\mathbb{R}^d}(x_i-\hat{\mu^{i}})(x_j-\hat{\mu^{j}})\mu(%
\mathrm{d }x), \ \ {i,j\in 1,\dots d,} \\
& |P_\mu| = \left(\sum_{i,j=1}^d(P_\mu^{ij})^2 \right)^{1/2}.
\end{align}

Next, we recall the definition of the weak topology on the space of probability
measures $\mathcal{P}(\mathbb{R}^{d})$ and note that,
of course, $\mathcal{P}^2(\mathbb{R}^{d})\subset \mathcal{P}(\mathbb{R}^{d})$.
We can thus consider also weak convergence of probability measures belonging to
the smaller space $\mathcal{P}^2(\mathbb{R}^{d})$:

A sequence of probability measures $\left(
\mu _{n}\right) _{n}\in \mathcal{P}(\mathbb{R}^{d})$, converges weakly
to $\mu \in \mathcal{P}(\mathbb{R}^{d})$ if and only if $\left( \mu
_{n}\left( \varphi \right) \right) _{n}$ converges to $\mu \left( \varphi
\right) $ as $n\rightarrow \infty $ for all $\varphi \in C_{b}(\mathbb{R}%
^{d})$. The weak topology on the space $\mathcal{P}(\mathbb{R}^{d})$ is
defined to be the weakest topology such that for all $f\in C_{b}(\mathbb{R}%
^{d})$, the function $\mu \mapsto \mu \left( f\right) $ is continuous.

The weak convergence of $\left( \mu _{n}\right) _{n}$ to $\mu $ is denoted $%
\mu _{n}\Rightarrow \mu $. A set of probability measures $\mathcal{A}\subset
\mathcal{P}(\mathbb{R}^{d})$ is relatively compact in the weak topology if and only if
for all $\varepsilon >0$ there exists $K_{\varepsilon }$ such that $\mu
\left( K_{\varepsilon }\right) \geq 1-\varepsilon $ for all $\mu
\in \mathcal{A}$. If $\mathcal{A}\subset \mathcal{P}^2(\mathbb{R}^{d})$, one
can show that the set $\mathcal{A}$ will be relatively compact in the weak
topology if the means and the covariance matrices of the probability
measures in the set $\mathcal{A}$ are uniformly bounded.

The set $\mathcal{A}$ can, however, still be relatively compact, even if the means of the
probability measures in it do not remain bounded. For
example, if we choose $\mu _{n}=\left( 1-\frac{1}{n}\right) \delta _{0}+%
\frac{1}{n}\delta _{n^{2}}, $ then the sequence $\mu _{n}$ is relatively
compact — in fact, $\mu _{n}\Rightarrow \delta_0 $ — but the
corresponding sequence of means $\hat \mu _{n} $ is not bounded, since $\mu _{n}\left( \varphi
\right) =n.$ On the other hand, if we choose $\mu _{n}=\left( 1-\frac{1}{n}\right)
\delta _{0}+\frac{1}{n}\delta _{n}, $ then the sequence $\mu _{n}$ is
relatively compact, the means $\hat\mu\equiv1$ form a
trivially bounded sequence, $\mu _{n}\left( \varphi \right) =1$, but the
second moments are not, as $\mu _{n}^2=n$. This state of affairs is not
satisfactory for our purposes.

The Wasserstein topology, though, adds the convergence of the first and second moments to
the weak convergence of the measures. To be precise, we have $%
\lim_{n\rightarrow \infty }{W_{2}({\mu }_{n},\mu )=0}$ for $\mu_{n},\mu \in
\mathcal{P}^2({\mathbb{R}^d})$, if and only if ${\mu }_{n} $ converges to $%
\mu $ in the weak topology and the first and second moments converge as well,
using the notation of Eqs. (<ref>)–(<ref>) above.
Moreover, a set of probability measures $\mathcal{%
A}\in \mathcal{P}^2({\mathbb{R}^d})$ is relatively compact in the topology
given by the Wasserstein distance if and only if

\begin{equation}
\lim_{R\mapsto\infty}\sup_{\mu\in\mathcal{A}}\int_{|x|>R}|x|^2\mu(\mathrm{d }%
\end{equation}

Finally, we have the following lemma which follows immediately from the alternative
definition (<ref>) of the Wasserstein distance $W_2$:
There exists a constant $C=C(d)$ such that, for any  $\mu ,\nu\in \mathcal{P}^2({\mathbb{R}^d})$,

\begin{eqnarray*}
(W_{2}(\mu ,\nu ))^2&\le& C (\mu^2+ \nu^2 ),\\
W_{2}(\mu ,\nu )&\le& C(|P_\mu|^{\frac{1}{2}}+|P_\nu|^{\frac{1}{2}}
\end{eqnarray*}

§.§ B. Assumptions and proof of Theorem <ref>

To start, we formulate here the set of assumptions on the coefficients of
the signal and observation equations (<ref>)+(<ref>)
under which Theorem <ref> holds:

* We assume that the coefficients $f$ and $h$ can be decomposed into a
linear part and a bounded nonlinear part.

In other words, we will assume that

\begin{equation}
f=F{\mathcal{I}} +\tilde{f},\ \ h=H{\mathcal{I}} +\tilde{h},  \label{Fhsigma}
\end{equation}

* ${\mathcal{I}} :\mathbb{R}^{d}\rightarrow \mathbb{R}^{d}$ is the
identity function defined as ${\mathcal{I}} (x)=x$ for any $x\in \mathbb{R}%

* $F\in \mathbb{R}^{d\times d}$, $H\in \mathbb{R}^{d\times n}$ are
given matrices; and

* $\tilde{f}:\mathbb{R}^{d}\rightarrow \mathbb{R}^{d}$ , $\tilde{h}:%
\mathbb{R}^{d}\rightarrow \mathbb{R}^{n}$ are bounded measurable functions
that incorporate the nonlinear parts of the coefficients of the system

* The covariance matrices of the processes $\pi ^{\mu }$ and $%
\pi $, respectively, are uniformly bounded in expectation:

\begin{equation}
\sup_{t\geq 0}\mathbb{E}[|P_{t}^{\mu }|^{8}]=C^{\mu }<\infty ,\ \
\sup_{t\geq 0}\mathbb{E}[|P_{t}^{\pi _{0}}|^{8}]=C^{\pi _{0}}<\infty .
\label{controlpt}
\end{equation}

* The matrix-valued process $Q_{s}^{\mu }:=F-P_{s}^{\mu }H^{\top
}H-\lambda _{s}^{\mu }H$ is exponentially stable in expectation, where $%
\lambda ^{\mu }$ is the matrix-valued process defined as $\lambda _{s}^{\mu
} := \pi _{s}^{\mu }({\mathcal{I}} -\hat{\pi}^{\mu })\tilde{h}^{\top }$ for
$s\geq 0$. In other words, if $\psi _{s:t}$ is the solution of the linear
matrix ordinary differential equation

\begin{equation*}
\mathrm{d}_{t}\psi _{s:t} =Q_{t}^{\mu }\psi _{s:t},\ \ \psi _{s:s}=I,
\end{equation*}

where $I$ is the identity matrix, then there exist some constant $c>0$ such

\begin{equation}
\mathbb{E}\left[ \left\vert \psi _{s:t}\right\vert ^{2}\right] \leqq
e^{-c\left( t-s\right) }.  \label{exp}
\end{equation}

Moreover, we assume that there exists $\delta >0$ such that

\begin{equation}  \label{psi-}
\int_{k\delta }^{\left( k+1\right) \delta }\mathbb{E}\left[ \left\vert \psi
_{k\delta :s}^{-1}\right\vert ^{4}\right] \mathrm{d}s\leq C^{inv}
\end{equation}

where $C^{inv}$ is a constant independent of $k$.

The proof of Theorem <ref> requires the bound of the difference $%
\hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}$ between the mean of the
FA process initialized from $\mu $ and the original
FA process, initialized from $\pi _{0}$. We deduce from
that, for $t\in \lbrack k\delta ,\left( k+1\right)
\delta ]$,

\begin{eqnarray*}
\hat{\pi}_{tk\delta }^{\mu } &=&\hat{\pi}_{k\delta }^{\mu }+\int_{k\delta
}^{t}\pi _{s}^{\mu }(f)\,\mathrm{d}s+\int_{k\delta }^{t}\varkappa _{s}^{\mu
}(\mathrm{d}Y_{s}-\pi _{s}^{\mu }(h)\,\mathrm{d}s), \\
\hat{\pi}_{t}^{\mu } &=&\hat{\pi}_{k\delta }^{\mu }+\int_{k\delta }^{t}\pi
_{s}^{\mu }(f)\,\mathrm{d}s+\int_{k\delta }^{t}\varkappa _{s}^{\mu }\left(
\mathrm{d}I_{s}+\left( \pi _{s}(h)-\pi _{s}^{\mu }(h)\right) \mathrm{d}%
s\right), \\
\hat{\pi}_{t} &=&\hat{\pi}_{k\delta }+\int_{k\delta }^{t}\pi _{s}(f)\,%
\mathrm{d}s+\int_{k\delta }^{t}\varkappa _{s}\mathrm{d}I_{s},
\end{eqnarray*}

\begin{eqnarray*}
\varkappa _{s}^{\mu } &:&=\pi _{s}^{\mu }({\mathcal{I}} h^{\top })-\pi
_{s}^{\mu }({\mathcal{I}} )\pi _{s}^{\mu }(h^{\top })=P_{s}^{\mu }H^{\top
}+\lambda _{s}^{\mu }, \\
\varkappa _{s} &:&=\pi _{s}({\mathcal{I}} h^{\top })-\pi _{s}({\mathcal{I}}
)\pi _{s}(h^{\top })=P_{s}^{\pi _{0}}H^{\top }+\lambda _{s},
\end{eqnarray*}

$P_{s}^{\mu }$, $P_{s}^{\pi _{0}}$ are the covariance matrices of $\pi
_s^{\mu }$ and of $\pi $, respectively, and $\lambda _{s}^{\mu }:=\pi _{s}^{\mu
}(({\mathcal{I}} -\hat{\pi}^{\mu })\tilde{h}^{\top })$ and $%
\lambda _{s}:=\pi _{s}(({\mathcal{I}} -\hat{\pi})\tilde{h}^{\top })$,  respectively. It
follows that

\begin{eqnarray}
\hat{\pi}_{t}^{\mu }-\hat{\pi}_{t} &=&\left( \hat{\pi}_{k\delta }^{\mu }-%
\hat{\pi}_{k\delta }\right) +\int_{k\delta }^{t}\left( \pi _{s}^{\mu }-\pi
_{s}\right) (f-\varkappa _{s}^{\mu }h)\,\mathrm{d}s+\int_{k\delta
}^{t}\left( \varkappa _{s}^{\mu }-\varkappa _{s}\right) \mathrm{d}I_{s}
\notag \\
&=&\left( \hat{\pi}_{k\delta }^{\mu }-\hat{\pi}_{k\delta }\right)
+\int_{k\delta }^{t}Q_{s}^{\mu }\left( \hat{\pi}_{s}^{\mu }-\hat{\pi}%
_{s}\right) \mathrm{d}s+\int_{k\delta }^{t}\left( P_{s}^{\mu }-P_{s}\right)
H^{\top }\mathrm{d}I_{s}+z_{k\delta :t}^{\mu },  \label{deltapi}
\end{eqnarray}

where $z_{k\delta :t}^{\mu }$ is a process that contains the nonlinearities
in the evolutions $\hat{\pi}^{\mu }$ and $\hat{\pi}$

\begin{equation*}
z_{k\delta :t}^{\mu }=\int_{k\delta }^{t}\left( \pi _{s}^{\mu }-\pi
_{s}\right) (\tilde{f}-\varkappa _{s}^{\mu }\tilde{h})\mathrm{d}%
s+\int_{k\delta }^{t}\left( \lambda _{s}^{\mu }-\lambda _{s}\right) \mathrm{d%
}I_{s},~~~t\in \lbrack k\delta ,\left( k+1\right) \delta ].
\end{equation*}

Replacing the observation process $Y_{t}$ by the innovation process $I_{t}$
in the evolution equations for $\hat{\pi}_{t}^{\mu }$ and $\hat{\pi}%
_{t}$, respectively, is important: Unlike $Y_t$, $I_{t}$ is a standard Brownian motion,
which enables us to use classical stochastic calculus properties to bound
the moments of the stochastic integrals appearing in Eq. (<ref>).

Using an argument based on the Grönwall inequality, one deduces that there exists a
constant $c_{\delta }$ independent of $k$ such that

\begin{equation}
\sup_{t\in \lbrack k\delta ,\left( k+1\right) \delta ]}\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \leq
c_{\delta }\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}%
_{k\delta }\right\vert \right] . \label{cd}
\end{equation}

We will show that there exists $R_{\delta }$ such that

\begin{equation}
\sup_{k\ge 0}\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}%
_{k\delta }\right\vert \right] \leq R_{\delta }.  \label{rd}
\end{equation}

From (<ref>) and (<ref>), one can then deduce that

\begin{eqnarray*}
\sup_{t\geq 0}\mathbb{E}\left[ \left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}%
_{t}\right\vert \right] &\leq& \sup_{k \ge 0}\sup_{t\in \lbrack k\delta
,\left( k+1\right) \delta ]}\mathbb{E}\left[ \left\vert \hat{\pi}_{t}^{\mu }-%
\hat{\pi}_{t}\right\vert \right] \\
&\leq& c_{\delta }\sup_{k \ge 0}\mathbb{E}%
\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}_{k\delta }\right\vert %
\right] \leq c_{\delta }R_{\delta }<\infty .
\end{eqnarray*}

Finally, from Lemma <ref> and (<ref>), it follows that

\begin{eqnarray*}
\sup_{t\geq 0}\mathbb{E}[W_{2}(\pi _{t}^{\mu },\pi _{t})] &\leq &\sup_{t\geq
0}\mathbb{E}\left[ \sqrt{|P_{t}^{\mu }|}\right] +\sup_{t\geq 0}\mathbb{E}%
\left [\sqrt{|P_{t}^{\pi _{0}}|}\right]+\sup_{t\geq 0}\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \\
&\leq &\left( \sup_{t\geq 0}\mathbb{E}\left[ |P_{t}^{\mu }|^{8}\right]
\right) ^{\frac{1}{16}}+\left( \sup_{t\geq 0}\mathbb{E}[|P_{t}^{\pi
_{0}}|^{8}]\right) ^{\frac{1}{16}}+\sup_{t\geq 0}\mathbb{E}\left[ \left\vert
\hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \\
&\leq &\left( C^{\mu }\right) ^{\frac{1}{16}}+\left( C^{\pi _{0}}\right) ^{%
\frac{1}{16}}+c_{\delta }R_{\delta }<\infty,
\end{eqnarray*}

which yields our claim.

To complete the proof it remains to show the validity of (<ref>). From (<ref>),
one obtains a form of Duhamel's principle for the difference between the values
of the mean $\hat{\pi}$ of the FA process at the steps $\left( k+1\right) \delta$ and
$k \delta$:

\begin{eqnarray}
\hat{\pi}_{\left( k+1\right) \delta }^{\mu }-\hat{\pi}_{k\delta }^{\mu }
&=&\psi _{k\delta ,\left( k+1\right) \delta }\left( \hat{\pi}_{k\delta
}^{\mu }-\hat{x}_{k\delta }\right) +\psi _{k\delta ,\left( k+1\right) \delta
}\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta
,,s}^{-1}(P_{s}^{\mu }-P_{s})H^{\top }\mathrm{d}I_{s}  \notag \\
&&+\psi _{k\delta ,\left( k+1\right) \delta }\int_{k\delta }^{\left(
k+1\right) \delta }\psi _{k\delta ,s}^{-1}\mathrm{d}z_{s}^{\mu }.
\label{deltapi2}
\end{eqnarray}

We analyze next each of the three terms on the right-hand side of (<ref>). For the first term, we use (<ref>) to derive the inequality

\begin{equation}
\mathbb{E}\left[ \left\vert \psi _{k\delta ,\left( k+1\right) \delta }\left(
\hat{\pi}_{k\delta }^{\mu }-\hat{x}_{k\delta }\right) \right\vert \right]
\leq e^{-c\delta }\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-%
\hat{x}_{k\delta }\right\vert \right].  \label{ld1}
\end{equation}

For the second term, we use the so-called Itô's integral isometry
property, cf. Karatzas and Shreve<cit.>, to obtain

\begin{eqnarray}
&&\hspace{-2.5cm}\mathbb{E}\left[ \left\vert \psi _{k\delta ,\left(
k+1\right) \delta }\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta
,,s}^{-1}(P_{s}^{\mu }-P_{s})H^{\top }\mathrm{d}I_{s}\right\vert \right] ^{2}
\notag \\
&\leq &\mathbb{E}\left[ \left\vert \psi _{k\delta ,\left( k+1\right) \delta
}\right\vert ^{2}\right] \mathbb{E}\left[ \left\vert \int_{k\delta }^{\left(
k+1\right) \delta }\psi _{k\delta ,,s}^{-1}(P_{s}^{\mu }-P_{s})H^{\top }%
\mathrm{d}I_{s}\right\vert ^{2}\right]  \notag \\
&\leq &e^{-2c\delta }\mathbb{E}\left[ \int_{k\delta }^{\left( k+1\right)
\delta }\left\vert \psi _{k\delta ,,s}^{-1}(P_{s}^{\mu }-P_{s})H^{\top
}\right\vert ^{2}\mathrm{d}s\right]  \notag \\
&\leq &e^{-2c\delta }\mathbb{E}\left[ \int_{k\delta }^{\left( k+1\right)
\delta }\left\vert \psi _{k\delta ,,s}^{-1}\right\vert ^{4}\mathrm{d}s\right]
+e^{-2c\delta }\mathbb{E}\left[ \int_{k\delta }^{\left( k+1\right) \delta
}\left\vert (P_{s}^{\mu }-P_{s})H^{\top }\right\vert ^{4}\mathrm{d}s\right]
\notag \\
&\leq &C^{I}e^{-2c\delta },  \label{ld2}
\end{eqnarray}

where $C^{I}=C^{inv}+\delta \left\vert H\right\vert ^{4}\sqrt{C^{\mu }}%
+\delta \left\vert H\right\vert ^{4}\sqrt{C^{\pi _{0}}}$.

Next, (<ref>) implies that there exists positive constants $%
C^{\varkappa ^{\mu }},C^{\varkappa },C^{\lambda ^{\mu }},C^{\lambda }$ such

\begin{equation}
\sup_{t\geq 0}\mathbb{E}[|\varkappa _{t}^{\mu }|^{8}]=C^{\varkappa ^{\mu
}}<\infty ,\ \ \sup_{t\geq 0}\mathbb{E}[|\varkappa _{t}|^{8}]=C^{\varkappa
}<\infty ,~~~\sup_{t\geq 0}\mathbb{E}[|\lambda _{t}^{\mu }|^{8}]=C^{\lambda
^{\mu }}<\infty ,~~~\sup_{t\geq 0}\mathbb{E}[|\lambda _{s}|^{8}]=C^{\lambda
}<\infty .  \label{othercontrol}
\end{equation}

For the third term, one gets

\begin{eqnarray}
\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta ,s}^{-1}\mathrm{d}%
z_{s}^{\mu }&=&\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta
,s}^{-1}\left( \pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa
_{s}^{\mu }\tilde{h})\mathrm{d}s  \notag \\
&&+\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta ,s}^{-1}\left(
\pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa _{s}^{\mu }\tilde{h}%
)\left( \lambda _{s}^{\mu }-\lambda _{s}\right) \mathrm{d}I_{s},  \label{ld3}
\end{eqnarray}

and one proceeds to bound separately the two terms in (<ref>). For the first term in (
<ref>), we get

\begin{eqnarray}
&&\hspace{-2cm}\mathbb{E}\left[ \left\vert \psi _{k\delta ,\left( k+1\right)
\delta }\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta
,s}^{-1}\left( \pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa
_{s}^{\mu }\tilde{h})\mathrm{d}s\right\vert \right]  \notag \\
&\leq &e^{-c\delta }\mathbb{E}\left[ \int_{k\delta }^{\left( k+1\right)
\delta }\left\vert \psi _{k\delta ,s}^{-1}\left( \pi _{s}^{\mu }-\pi
_{s}\right) (\tilde{f}-\varkappa _{s}^{\mu }\tilde{h})\right\vert \mathrm{d}s%
\right]  \notag \\
&\leq &e^{-c\delta }\left( \mathbb{E}\left[ \int_{k\delta }^{\left(
k+1\right) \delta }\left\vert \psi _{k\delta ,s}^{-1}\right\vert ^{2}\mathrm{%
d}s\right] +\mathbb{E}\left[ \int_{k\delta }^{\left( k+1\right) \delta
}\left\vert \left( \pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa
_{s}^{\mu }\tilde{h})\right\vert ^{2}\mathrm{d}s\right] \right)  \notag \\
&\leq &e^{-c\delta }\left( C^{inv}+\delta +\left\vert \tilde{f}\right\vert
^{2}\delta +\left\vert \tilde{h}\right\vert ^{2}\left( C^{\varkappa ^{\mu
}}\right) ^{\frac{1}{4}}\delta \right) .  \label{ld4}
\end{eqnarray}

For the second term in (<ref>), we get that

\begin{eqnarray}
&&\hspace{-1cm}\mathbb{E}\left[ \left\vert \psi _{k\delta ,\left( k+1\right)
\delta }\int_{k\delta }^{\left( k+1\right) \delta }\psi _{k\delta
,s}^{-1}\left( \pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa
_{s}^{\mu }\tilde{h})\left( \lambda _{s}^{\mu }-\lambda _{s}\right) \mathrm{d%
}I_{s}\right\vert \right]  \notag \\
&\leq &e^{-c\delta }\mathbb{E}\left[ \left\vert \int_{k\delta }^{\left(
k+1\right) \delta }\psi _{k\delta ,s}^{-1}\left( \pi _{s}^{\mu }-\pi
_{s}\right) (\tilde{f}-\varkappa _{s}^{\mu }\tilde{h})\left( \lambda
_{s}^{\mu }-\lambda _{s}\right) \mathrm{d}I_{s}\right\vert ^{2}\right]
\notag \\
&\leq &e^{-c\delta }\int_{k\delta }^{\left( k+1\right) \delta }\mathbb{E}%
\left[ \left\vert \psi _{k\delta ,s}^{-1}\left( \pi _{s}^{\mu }-\pi
_{s}\right) (\tilde{f}-\varkappa _{s}^{\mu }\tilde{h})\left( \lambda
_{s}^{\mu }-\lambda _{s}\right) \right\vert ^{2}]\mathrm{d}s\right]  \notag
\\
&\leq &e^{-c\delta }\int_{k\delta }^{\left( k+1\right) \delta }\mathbb{E}%
\left[ \left\vert \psi _{k\delta ,s}^{-1}\right\vert ^{4}]\mathrm{d}s\right]
+e^{-c\delta }\int_{k\delta }^{\left( k+1\right) \delta }\mathbb{E}\left[
\left\vert \left( \pi _{s}^{\mu }-\pi _{s}\right) (\tilde{f}-\varkappa
_{s}^{\mu }\tilde{h})\left( \lambda _{s}^{\mu }-\lambda _{s}\right)
\right\vert ^{4}]\mathrm{d}s\right]  \notag \\
&\leq &e^{-c\delta }\left( C^{inv}+\delta \left\vert \tilde{f}\right\vert
^{4}\left( \sqrt{C^{\lambda ^{\mu }}}+\sqrt{C^{\lambda }}\right) +\delta
\left\vert \tilde{h}\right\vert ^{4}\left( C^{\varkappa ^{\mu
}}+C^{\varkappa }\right) \right) .  \label{ld5}
\end{eqnarray}

From (<ref>), (<ref>) and (<ref>), it follows that there exists a

\begin{equation*}
C^{z^{\mu }}=C^{z^{\mu }}\left( \delta ,C^{inv},\tilde{f},C^{\lambda ^{\mu
}},C^{\lambda },\tilde{h},C^{\varkappa ^{\mu }},C^{\varkappa }\right)
\end{equation*}

independent of $k$ such that

\begin{equation}
\mathbb{E}\left[ \left\vert \int_{k\delta }^{\left( k+1\right) \delta }\psi
_{k\delta ,s}^{-1}\mathrm{d}z_{s}^{\mu }\right\vert \right] \leq C^{z^{\mu
}}e^{-c\delta } .  \label{ld6}
\end{equation}

Finally, from (<ref>), (<ref>), (<ref>) and (<ref>) we
deduce that

\begin{equation}
\mathbb{E}\left[ \left\vert \hat{\pi}_{\left( k+1\right) \delta }^{\mu }-%
\hat{x}_{\left( k+1\right) \delta }\right\vert \right] \leq e^{-c\delta
}\left( \mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}%
_{k\delta }\right\vert \right] +C^{I}+C^{z^{\mu }}\right) .  \label{ld7}
\end{equation}

Choose now $R=R\left( \delta \right) =\max \left( \left\vert \hat{\mu}-\hat{%
\pi}_{0}\right\vert ,\frac{C^{I}+C^{z^{\mu }}}{c\delta }\right) $ and use
induction to prove (<ref>). From the definition of $R$, we deduce that $%
\left\vert \hat{\mu}-\hat{\pi}_{0}\right\vert \leq R$. Next assume that $%
\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-\hat{x}_{k\delta
}\right\vert \right] \leq R$. From this and (<ref>), one can obtain that

\begin{equation*}
\mathbb{E}\left[ \left\vert \hat{\pi}_{\left( k+1\right) \delta }^{\mu }-%
\hat{x}_{\left( k+1\right) \delta }\right\vert \right] \leq e^{-c\delta
}\left( R+C^{I}+C^{z^{\mu }}\right) \leq e^{-c\delta }R\left( 1+\frac{%
C^{I}+C^{z^{\mu }}}{R}\right) \leq e^{-c\delta }R\left( 1+c\delta \right)
\leq R.
\end{equation*}

It follows that $\mathbb{E}\left[ \left\vert \hat{\pi}_{k\delta }^{\mu }-%
\hat{x}_{k\delta }\right\vert \right] \leq R$ holds true for any $k \ge 0$
and so does (<ref>). The proof of Theorem <ref> is now

The proof of the Theorem <ref> relies on the application of Duhamel's principle, as in (<ref>). The Duhamel principle can be applied on an arbitrary interval $[0,t]$ to deduce that

\begin{equation}
\hat{\pi}_{t }^{\mu }-\hat{\pi}_{k\delta }^{\mu }
=\psi _{0 ,t }\left( \hat{\pi}_{0
}^{\mu }-\hat{x}_{0 }\right) +\psi _{0 ,t }\int_{0 }^{t}\psi _{0,s}^{-1}(P_{s}^{\mu }-P_{s})H^{\top }\mathrm{d}I_{s}
+\psi _{0,t }\int_{0 }^{t}\psi _{0,s}^{-1}\mathrm{d}z_{s}^{\mu }.
\label{deltapi2'}
\end{equation}

For a deterministic $\psi$, we can rewrite (<ref>) as

\begin{equation}
\hat{\pi}_{t }^{\mu }-\hat{\pi}_{k\delta }^{\mu }
=\psi _{0 ,t }\left( \hat{\pi}_{0
}^{\mu }-\hat{x}_{0 }\right) +\int_{0 }^{t}\psi _{s,t}(P_{s}^{\mu }-P_{s})H^{\top }\mathrm{d}I_{s}
+\int_{0 }^{t}\psi _{s,t}\mathrm{d}z_{s}^{\mu }.
\label{deltapi2''}
\end{equation}

In this case, the bound of the stochastic terms in (<ref>) would follow directly from (<ref>) or, more precisely, from the deterministic exponential decay of $\psi_{s,t}$. However, in our case, $\psi$ is not deterministic. Moreover the process $\psi_{s,t}$ is not adapted with respect to the filtration generated by the Brownian motion $I$.  To be more precise, it does not depend on $\{I_r, r\in [0,s]\}$ but on $\{I_r, r\in [0,t]\}$.      Hence, the stochastic integrals in (<ref>) do not make sense as standard Itô integrals. One can interpret them using a more general definition of Skorohod integration, but the control of such resulting integrals is no longer immediately obvious. A bound using Malliavin calculus may be possible;  see the monograph of Nualart<cit.> for details of the methodology.

In another approach, one could attempt to keep out of the stochastic integrals the part that is
not adapted, i.e. $\psi_{0,t}$, and only use (<ref>) but not (<ref>). To do so, one would need to bound the exponential blow-up of the stochastic integrals in (<ref>) as $t\rightarrow \infty$
over the entire positive half line $[0,\infty)$. Instead, we limit ourselves to apply Duhamel's priciple on intervals of the form $[k\delta, (k+1)\delta]$, where $\delta$ is small enough to be able to rely on the bound given by (<ref>).

Finally, note that, using again an argument based on the Grönwall inequality, one obtains
that, for any $T>0$, there exists a constant $c_{T }$ such that the expected difference between the means of the FA process started from different initial distributions $\mu$ and $\pi_0$ satisties

\begin{equation}
\sup_{t\in [0,T]}\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \leq
c_{T }. \label{cd'}
\end{equation}

Thus, we only need to check that there exists a time horizon $T$ for which

\begin{equation}
\sup_{t\in [T,\infty)}\mathbb{E}\left[
\left\vert \hat{\pi}_{t}^{\mu }-\hat{\pi}_{t}\right\vert \right] \leq
c_{T } \label{cd''}
\end{equation}

The bound of the eighth moments of the covariance matrices $P_{t}^{\mu }$ and $P_{t}^{\pi _{0}}$,  respectively, in (<ref>), as well as of the fourth moments of $\psi_{k\delta :s}^{-1}$, suffices for the result above, but it is perhaps not necessary. We leave it for future research to find the optimal bounds.

In the linear and Gaussian case, $\tilde f$ and $\tilde h$ are constant and $\lambda^\mu_s=0$. Therefore

\[
z^{\mu}=\int_{0}^{t}(\tilde f-P_{s}^\mu H^{\top }\tilde h)\mathrm{d}s;
\]

in other words, $z^\mu$ is a process with bounded variation. Moreover, $P_{s}^\mu$ is deterministic. In this case, all the technical difficulties described above vanish and Theorem <ref> holds true with a much simpler proof. In this case, though, a much stronger result holds, namely Theorem <ref>, which is proved in the next section of this appendix.

Theorem <ref> gives us only an upper bound on the expected value of the distance between $\pi_{t}^\mu$ and  $\pi_t$. From it, we can derive the following alternative control: Let us denote by $T(t,R)$ the average time spent by the process $\pi^{\mu}_t $ outside the ball $B(\pi_t,R)$, that is
\[
T(t,R)=\frac{1}{t}\int_0^tI_{\{W_{2}(\pi _{t}^{\mu },\pi_{t}^{0})>R\}}\mathrm{d}s.
\]
Then, for any $\varepsilon>0$, there exists a constant
$R=R({\mu},\varepsilon)$ independent of $t$ such that

\begin{equation*}
\sup_{t\ge 0}\mathbb E\left [T(t,R)\right]<\varepsilon.
\end{equation*}

We note that we cannot expect $T(t,R)$ to decrease to zero as $t$ tends to $\infty$, since $\hat\pi^{\mu}_t-\hat{x}_t $  can be viewed as a  $d$-dimensional Ornstein-Uhlenbeck process with random coefficients and perturbed  by the random residue process  $z$.

The concept of stabilization described in Theorem <ref> is much weaker than the one described in many classical stabilization results, e.g. such as those covered by Budhiraja <cit.>, Picard<cit.> or Van Handel<cit.>. On the other hand, it has the advantage that the signal $X(t)$ is not required to be fully observable. In particular, the dimension $n$ of the observation $Y(t)$ and the dimension $d$ of the signal $X(t)$, in the notation of Sec. <ref>, do not need to coincide. Heuristically we need to be able to observe all the "unstable" directions. See Section D in this Appendix
for some simple illustrative examples in this direction. This feature of our results is very advantageous in practical applications as, in many situations, the dimension of the observation space is considerably smaller than that of the state space, $n \ll d$  <cit.>.

One can deduce a criterion for the bound (<ref>) to hold. The derivation of such a criterion can be carried out in terms of the centralized third moments of $\pi_t^{\mu}$, along the lines of the arguments in Section 6.2 of Bain and Crisan <cit.>, and it is the subject of future research.

§.§ C. Assumptions and proof of Theorem <ref>

The assumptions and the arguments in this section  are based on the framework of
Ocone and Pardoux<cit.> and, in part, on their results. In particular, we impose the
following set of assumptions on the coefficients of the signal and observation equations
(<ref>)+(<ref>) under which Theorem <ref> holds:

* We assume that the coefficients $f$ and $h$ are linear.

In other words, we will assume that

\begin{equation}
f=F{\mathcal{I}}+\tilde{f},\ \ h=H{\mathcal{I}}+\tilde{h},  \label{Fhsigmal}
\end{equation}

* ${\mathcal{I}}:\mathbb{R}^{d}\rightarrow \mathbb{R}^{d}$ is the
identity function defined as ${\mathcal{I}}(x)=x$ for any $x\in \mathbb{R}%

* $F\in \mathbb{R}^{d\times d}$, $H\in \mathbb{R}^{d\times n}$ are
given matrices;

* $\tilde{f}$ and $\tilde{h}$ are $d$-dimensional and $n$
-dimensional vectors, respectively; and

* the function $\sigma $ is a constant $d\times d$-matrix.

* We assume that the measure $\mu $ has finite second
moments — and thus it it belongs to the Wasserstein space, as discussed in
Appendix I — and that it is absolutely continuous
with respect to $\pi _{0}$. We denote by $\theta _{\mu }$ the density of $%
\mu $ with respect to $\pi _{0}$ which is  integrable
with respect to $\pi _{0}.$
Using a standard probabilistic result (see, for
example, Problem 3.20 in Karatzas & Shreve <cit.>), it follows that there
exists a random variable $\Upsilon _{\mu }\geq 0$ such that $%
\lim_{t\rightarrow \infty }\mathbb{E}\left[ \theta _{\mu }\left(
X_{0}\right) |\mathcal{Y}_{t}\right] =\Upsilon _{\mu }$. Moreover both $%
\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}\right] $
and $\Upsilon _{\mu }$  have a pathwise representation in terms of the observation path $Y_{t}.$

We assume that, for the path for which we do the analysis, $\lim_{t\rightarrow \infty }%
\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}\right]
=\Upsilon _{\mu }>0.$ Similarly, we will be assuming that both $\lim_{t\rightarrow
\infty }\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) X_{0}|\mathcal{Y}%
_{t}\right] $ and $\lim_{t\rightarrow \infty }\mathbb{E}\left[ \theta _{\mu
}\left( X_{0}\right) X_{0}X_{0}^*|\mathcal{Y}_{t}\right] $ exist,
based on the same standard probabilistic result.

* We define $\pi _{t}^{0}:=N\left( \hat{x}_{t}^{0},P_{t}^{0}\right),$
where $\hat{x}$ satisfies the stochastic differential equation

\begin{multline}
\mathrm{d}\hat{x}_{t}^{0} = \left( F\hat{x}_{t}^{0}+\tilde{f}\right) \,\mathrm{%
d}t+P_{t}H^{\top }(\mathrm{d}Y_{t}-(H\hat{x}_{t}^{0}+\tilde{h})\,\mathrm{d}%
t) \\
= (Q_{t}\hat{x}_{t}^{0}+\tilde{f}-P_{t}H^{\top }\tilde{h})\,\mathrm{d}%
t+P_{t}H_{t}^{\top }\mathrm{d}Y_{t},\quad \hat{x}_{0}=0
\label{eq:kalman:conditionalmean}
\end{multline}

and $P^{0}$ satisfies the deterministic matrix Riccati equation

\begin{multline}
\frac{\mathrm{d}P_{t}^{0}}{\mathrm{d}t} = \sigma _{t}\sigma _{t}^{\top
}+FP_{t}^{0}+P_{t}^{0}F^{\top }-P_{t}^{0}H^{\top}HP_{t}^{0}\\
= FQ_{t}+Q_{t}F^{\top }+\sigma _{t}\sigma _{t}^{\top
}+P_{t}H^{\top }HP_{t},  \label{eq:kalman:conditionalcovariance}
\end{multline}

with $P_{0}=0\in \mathbb{R}^{d\times d}$ being the matrix with null entries.
We assume, furthermore, that there exists a unique solution $P_{\infty }\geq 0$ to
the algebraic Riccati equation

\begin{equation}
\sigma \sigma ^{\top }+FP+PF^{\top }-PH^{\top }HP = 0,  \label{assym}
\end{equation}

and that $Q_{t}:=F-P_{t}^{0}H^{\top }H$ is asymptotically stable.
Using this unique non-negative solution, we
define $Q_{\infty }:=F-P_{\infty }H^{\top }H$.

Following from Remark 2.1 and Lemma 2.2 in Ocone and Pardoux<cit.>, which
in turn uses the classical result of Theorem 3.7 in Kwakernaak and Sivan <cit.>, we make the following remark:

If $(F,H)$ is detectable and  $(F, \sigma_{0})$ is stabilizable, then (<ref>) has a unique non-negative definite solution $P_{\infty}$ and one has both $P_{t} \rightarrow P_{\infty}$ and $Q_{t} \rightarrow Q_{\infty}$  for any initial condition $P_{0}^{} \geq 0$, with $Q_t$ asymptotically stable. See Chapter 1 in Kwakernaak and Sivan <cit.> for the definition of detectability and stabilizability. In particular, the eigenvalues of $Q_\infty$ have negative real parts and  for any $0<a<\inf \left\{-\operatorname{Re} \lambda ; {\rm over} \,\, \lambda\right.$ being an eigenvalue of $\left. Q_\infty \right\}$, there is a constant $K_{a}$ such that

\left\|P_{t}^{}-P_{\infty}\right\| \leq K_{a} e^{-at}

The last fact can be proved by observing that

\[
\frac{d}{d t}\left(P_{t}^{}-P_{\infty}\right)=\left[F-1 / 2\left(P_{t}^{}+P_{\infty}\right) H^{\top} H\right]\left(P_{t}^{}-P_{\infty}\right)
+\left(P_{t}^{}-P_{\infty}\right)\left[F-1 / 2\left(P_{t}^{R}+P_{\infty}\right) H^{\top} H\right]^{\top}
\]

and carrying out an analysis similar to that in Theorem 2.3 of Ocone and Pardoux<cit.>.

We will show below that $\pi ^{\mu }$ gets asymptotically close to the
“reference process” $\pi ^{0}$, regardless of the initial condition, and
since this holds true for $\mu =\pi _{0}$ too, we immediately deduce that

\begin{equation*}
\lim_{t\rightarrow \infty }W_{2}(\pi _{t}^{\mu },\pi _{t})=0.
\end{equation*}

Note that $\pi $ coincides with $\pi ^{0}$ if the initial condition $\pi
_{0} $ is the Dirac delta distribution at $0$, that is, $\pi _{0}=\delta _{0}$. The
reference process $\pi ^{0}$ is convenient to work with: From Remark
<ref> we deduce that the centered version of $\pi ^{0,-\hat{x}%
_{t}^{0}}:=N\left( 0,P_{t}^{0}\right) $ converges weakly, as well as in
Wasserstein distance, to $\pi ^{\infty }=N\left( 0,P_{\infty }\right) $. Using
the equivalent definition of the Wasserstein distance
(<ref>), it follows that

\begin{equation*}
\lim_{t\rightarrow \infty }W_{2}(\pi _{t}^{\mu },\pi
_{t}^{0})=0\Longleftrightarrow \lim_{t\rightarrow \infty }W_{2}(\pi
_{t}^{\mu ,-\hat{x}_{t}^{0}},\pi _{t}^{0,-\hat{x}_{t}^{0}})=0%
\Longleftrightarrow \lim_{t\rightarrow \infty }W_{2}(\pi _{t}^{\mu ,-\hat{x}%
_{t}^{0}},\pi ^{\infty })=0.
\end{equation*}

To prove the last limit, it suffices to show that $\pi _{t}^{\mu ,-\hat{x}%
_{t}^{0}}$ converges to $\pi ^{\infty }$ in the weak topology and the first
and second moments of $\pi _{t}^{\mu ,-\hat{x}_{t}^{0}}$ converge to the
first and second moments of $\pi ^{\infty }$, respectively. Equivalently, it
suffices to show that:

* $\pi _{t}^{\mu ,-\hat{x}_{t}^{0}}$ converges to $\pi^{\infty }$ in the weak topology;

* $\lim_{t\rightarrow \infty}| \hat{\pi}_{t}^{\mu }-\hat{x}_{t}^{0}| =0$ ,
i.e., the distance betweeen the mean of $\pi _{t}^{\mu }$ and that of $\pi _{t}^{0}$
tends to $0$; and

* $\lim_{t\rightarrow \infty }| P_{\pi_{t}^{\mu }}-P_{t}^{0}| =0$, i.e., the distance
betweeen the covariance matrix of $\pi _{t}^{\mu }$ and that of $\pi _{t}^{0}$ tends to $0$.

Since the set of bounded uniformly continuous functions are convergence determining, to justify that
$\pi _{t}^{\mu ,-\hat{x}_{t}^{0}}$ converges to $\pi ^{\infty }$ in the weak topology it suffices to
show that
$\lim_{t\rightarrow \infty }|\pi _{t}^{\mu ,-%
\hat{x}_{t}^{0}}\left( \varphi \right) -\pi ^{\infty }\left( \varphi \right)
| =0$, for any $\varphi $ that is a bounded and uniformly continuous function.
To recap, we have that $\lim_{t\rightarrow \infty}W_{2}(\pi _{t}^{\mu },\pi _{t}^{0})=0$,
if and only if the following three
properties hold true:

* $\lim_{t\rightarrow \infty }| \hat{\pi}_{t}^{\mu }-\hat{x}%
_{t}^{0}| =0;$

* $\lim_{t\rightarrow \infty }\left\vert P_{\pi _{t}^{\mu
}}-P_{t}^{0}\right\vert =0;$ and

* $\lim_{t\rightarrow \infty }\left\vert \pi _{t}^{\mu }\left( \varphi
_{t}\right) -\pi _{t}^{0}\left( \varphi _{t}\right) \right\vert =0$ for any
bounded uniformly continuous function $\varphi $, where $\varphi
_{t}$ is the same function shifted by the mean $\hat{x}_{t}^{0},$ that is $%
\varphi _{t}\left( x\right) :=\varphi _{t}\left( x+\hat{x}_{t}^{0}\right) ,$
$x\in \mathbb{R}^{d}$.

Note that the reference process $\pi ^{0}$ is not unique: we can replace
it by any other process $\tilde{\pi}^{0}$ with the property that $%
\lim_{t\rightarrow \infty }W_{2}(\tilde{\pi}_{t}^{0},\pi _{t}^{0})=0$.

The proof of Theorem <ref> is largely based on the arguments in the
proof of Theorem 2.6 in Ocone and Pardoux<cit.>.
Using the linearity of the equation (<ref>), we can
decompose the signal as follows:

\begin{eqnarray*}
X_{t} &=&e^{Ft}X_{0}+\tilde{X}_{t}, \\
\tilde{X}_{t} &:&=\int_{0}^{t}(F\tilde{X}_{s}+f)\mathrm{d}%
s+\int_{0}^{t}\sigma \mathrm{d}V_{s}=\int_{0}^{t}e^{F(t-s)}(\tilde{f}%
\mathrm{d}s+\sigma \mathrm{d}V_{s}), \\
Y_{t} &=&\int_{0}^{t}(He^{Fs}X_{0}+H\tilde{X}_{s}+\tilde{h})\mathrm{d}%
\end{eqnarray*}

and we introduce the new measure $\overline{\mathbb{P}}$ defined by

\begin{equation*}
\left.\frac{\mathrm{d}\overline{\mathbb{P}}}{\mathrm{d}\mathbb{P}}\right |_{\mathcal{F}_t}=\exp \left[ \int_{0}^{t}-\left%
\langle He^{Fs}X_{0},\mathrm{d}W_{s}\right\rangle -\frac{1}{2}\int_{0}^{t}\left\vert
He^{Fs}X_{0}\right\vert ^{2}\mathrm{d}s\right],\ \ \ t\ge 0.
\end{equation*}

By Girsanov's theorem, under this measure $\overline{\mathbb{P}},$ the process

\begin{equation*}
\bar{W}_{t}:=\int_{0}^{t}He^{Bs}X_{0}\mathrm{d}s+W_{t},\quad t\leq T
\end{equation*}

is a Brownian motion, and $X_{0}$ is independent of $\left( V_{t},\bar{W}%
_{t}\right) _{t\leq T}$.  Following the proof of Proposition 3.13  in Bain and Crisan<cit.>, the law of $X_{0}$ remains unchanged under $\bar {\mathbb P}$. Let

\begin{equation*}
L_{t}:= \left.\frac{\mathrm{d}\mathbb{P}}{\mathrm{d}\overline{\mathbb{P}}}\right |_{\mathcal{F}_t}=\exp \left[ \int_{0}^{t}\left\langle He^{Bs}X_{0},\mathrm{d}\bar{W}%
_{s}\right\rangle -\frac{1}{2}\int_{0}^{t}\left\vert He^{Bs}X_{0}\right\vert
^{2}\mathrm{d}s\right] =\exp \left( \left\langle X_{0},W^{\circ }\right\rangle -\frac{%
1}{2}\left\langle M_{t}X_{0},X_{0}\right\rangle \right) ,
\end{equation*}

\begin{equation*}
\left( W_{t}^{o},M_{t}\right) =\left( e^{F^*t}\int_{0}^{t}H^*\mathrm{d}\bar{W}%
\end{equation*}

Let $%
\psi :$ $\mathbb{R}^{d}\times \mathbb{R}^{d}\rightarrow \mathbb{R}$ be a function that is
integrable with respect to the joint law of $\left( X_{0},\tilde{X}%
_{t}\right)$. Then

\begin{equation*}
\mathbb{E}\left[ \psi \left( X_{0},\tilde{X}_{t}\right) |\mathcal{Y}_{t}%
\right] =\frac{\bar{\mathbb{E}}\left[ \psi \left( X_{0},\tilde{X}_{t}\right)
L_{t}|\mathcal{Y}_{t}\right] }{\bar{\mathbb{E}}\left[ L_{t}|\mathcal{Y}_{t}%
\right] }.
\end{equation*}

From the above formula for $\tilde{X}_{t}$ and the fact that $X_{0}$ and $(%
\bar{W},Y)$ are $\overline{\mathbb{P}}$-independent, we deduce that

\begin{equation}
\mathbb{E}\left[ \psi \left( X_{0},\tilde{X}_{t}\right) |\mathcal{Y}_{t}%
\right] =\frac{\int_{\mathbb{R}^{d}}e^{-1/2\left( M_{t}x,x\right) }\bar{%
\mathbb{E}}\left[ \psi \left( x,\tilde{X}_{t}\right) \exp \left(
\left\langle x,W^{\circ }\right\rangle \right) |\mathcal{Y}_{t}\right] \pi
_{0}(\mathrm{d}x)}{\int_{\mathbb{R}^{d}}e^{-1/2\left( M_{t}x,x\right) }\bar{\mathbb{E}%
}\left[ \exp \left( \left\langle x,W^{\circ }\right\rangle \right) |\mathcal{%
Y}_{t}\right] \pi _{0}(\mathrm{d}x)}.  \label{pair}
\end{equation}

In particular, for any function $\varphi :\mathbb{R}^{d}\rightarrow \mathbb{R%
}$ such that the random variable $\varphi \left( X_{t}\right) =\varphi
\left( e^{Ft}X_{0}+\tilde{X}_{t}\right) $ is integrable, we deduce from (<ref>) that

\begin{eqnarray*}
\mathbb{E}\left[ \varphi \left( X_{t}\right) |\mathcal{Y}_{t}\right] &=&%
\frac{\int_{\mathbb{R}^{d}}e^{-1/2\left( M_{t}x,x\right) }\bar{\mathbb{E}}%
\left[ \varphi _{x}^{1}\left( \tilde{X}_{t},W^{\circ }\right) |\mathcal{Y}%
_{t}\right] \pi _{0}(\mathrm{d}x)}{\int_{\mathbb{R}^{d}}e^{-1/2\left( M_{t}x,x\right)
}\bar{\mathbb{E}}\left[ \varphi _{x}^{2}\left( \tilde{X}_{t},W^{\circ
}\right) |\mathcal{Y}_{t}\right] \pi _{0}(\mathrm{d}x)}; \\
\varphi _{x}^{1}\left( a,b\right) &=&\varphi \left( e^{Ft}x+a\right) \exp
\left\langle x,b\right\rangle ,~~a\in \mathbb{R}^{d}, \, b\in \mathbb{R}^{n}; \,\, {\rm and} \\
\varphi _{x}^{2}\left( a,b\right) &=&\exp \left\langle x,b\right\rangle,
~~~~a\in \mathbb{R}^{d}, \, b\in \mathbb{R}^{n}.
\end{eqnarray*}

The pair $\left( \tilde{X},W^{\circ }\right) $ satisfies a linear system of
stochastic differential equations driven by the $\left( d+m\right) $-Brownian motion $\left( V,W^{\circ
}\right) $ with initial condition $\left( 0,0\right) $. Moreover, under $%
\overline{\mathbb{P}},$ the process $Y$ satisfies

\begin{equation}\label{newobs}
\end{equation}

It follows that we can express $\bar{\mathbb{E}}\left[ \varphi
_{x}^{i}\left( \tilde{X}_{t},W_{t}^{\circ }\right) |\mathcal{Y}_{t}\right]
,i=1,2,$ as integrals with respect to a Gaussian distribution $\eta $
with mean vector  $\left( \hat{x}_{t}^{0},\hat{x}%
_{t}^{W}\right)^*$ and covariance matrix

\begin{equation*}
\begin{array}{cc}
P_{t}^{0} & S_{t} \\
S_{t} & R_{t}%
\end{array}%
\right)
\end{equation*}

that satisfy the equations satisfied by a Kalman-Bucy filter with the signal equations identical to those for the pair
$\left( \tilde{X}_{t},W_{t}^{\circ }\right)$ and the observation equation identical to (<ref>).

More precisely,   $\hat{x}_{t}^{0}$ is the solution of equation (<ref>), $P_{t}^{0}$ is the solution of the equation (<ref>), $\hat{x}_{t}^{W}$ is the solution of the

\begin{equation*}
\mathrm{d}\hat{x}_{t}^{W}=\left( e^{Ft}+S_{t}\right) ^{\ast }H^{\ast }\left(
\mathrm{d}Y_{t}-\left( H\hat{x}_{t}^{0}+\tilde{h}\right) \mathrm{d}t\right) ,\quad \hat{x}%
\end{equation*}

and the pair $\left( S_{t},Q_{t}\right) $
solve the equation s

\begin{eqnarray*}
\frac{\mathrm{d}S_{t}}{\mathrm{d}t} &=&FS_{t}-P_{t}^{0}H^{\ast }H\left( e^{Ft}+S_{t}\right)
,\quad S_{0}=0, \\
\frac{\mathrm{d}R_{t}}{\mathrm{d}t} &=&-e^{F^{\ast }t}H^{\ast }HS_{t}-S_{t}^{\ast }H^{\ast
}He^{Ft}-S_{t}^{\ast }H^{\ast }HS_{t},\quad R_{0}=0;
\end{eqnarray*}

see formulae (27–33) in Ocone and Pardoux<cit.> for details.
One can check that

\begin{equation*}
\frac{\mathrm{d}S^\circ_t}{\mathrm{d}t} =\left( F-P_{t}^{0}H^{\ast }H\right)
\end{equation*}

where $S^\circ_t:=e^{Ft}+S_{t}$ and, since $Q_{t}$ is asymptotically stable, we conclude that

\begin{equation}
\lim_{t\rightarrow \infty }\left\vert \left\vert S^\circ_t\right\vert
\right\vert =0.  \label{eq:sumas}
\end{equation}

In fact, similar to Remark <ref>, the convergence in (<ref>) has an exponential decay rate to 0.

We  deduce that

\begin{eqnarray}
\bar{\mathbb{E}}\left[ \varphi _{x}^{2}\left( \tilde{X}_{t},W^{\circ
}\right) |\mathcal{Y}_{t}\right] &=&\bar{\mathbb{E}}\left[ \exp \left\langle
x,W^{\circ }\right\rangle |\mathcal{Y}_{t}\right] = e^{\left\langle x,\hat{x}%
_{t}^{W}\right\rangle +1/2\left( Q_{t}x,x\right) } ; \label{help} \\
\bar{\mathbb{E}}\left[ \varphi _{x}^{1}\left( \tilde{X}_{t},W^{\circ
}\right) |\mathcal{Y}_{t}\right] &=&\bar{\mathbb{E}}\left[ \varphi \left(
e^{Ft}x+\tilde{X}_{t}\right) \exp \left\langle x,W^{\circ }\right\rangle |%
\mathcal{Y}_{t}\right]  \notag \\
&=&e^{\left\langle x,\hat{x}_{t}^{W}\right\rangle +1/2\left( Q_{t}x,x\right)
}\int \varphi \left( S^\circ_tx+a\right) \pi
_{t}^{0}\left( \mathrm{d}a\right) .  \notag
\end{eqnarray}

It follows that

\begin{equation*}
\pi _{t}\left( \varphi \right) =\mathbb{E}\left[ \varphi \left( X_{t}\right)
/\mathcal{Y}_{t}\right] =\frac{\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right)
\int \varphi \left( S^\circ_tx+a\right) \pi _{t}^{0}\left(
\mathrm{d}a\right) \pi _{0}(\mathrm{d}x)}{\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \pi
\end{equation*}

where $\Xi _{t}\left( x\right) :=e^{-1/2\left( M_{t}x,x\right) +\left\langle
x,\hat{x}_{t}^{W}\right\rangle +1/2\left( Q_{t}x,x\right) }.$ Moreover, by
choosing a function $\psi :$ $\mathbb{R}^{d}\times \mathbb{R}^{d}\rightarrow
\mathbb{R}$ in (<ref>) that is independent in the second component, $%
\psi \left( x,y\right) =\varphi \left( x\right) $, $\left( x,y\right) \in
\mathbb{R}^{d}\times \mathbb{R}^{d}$, we obtain from (<ref>) that

\begin{equation*}
\mathbb{E}\left[ \varphi \left( X_{0}\right) |\mathcal{Y}_{t}\right] =\frac{%
\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \varphi \left( x\right) \pi
_{0}(\mathrm{d}x)}{\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \pi _{0}(\mathrm{d}x)}=\frac{%
\pi _{0}\left( \Xi _{t}\varphi \right) }{\pi _{0}\left( \Xi _{t}\right) }.
\end{equation*}

We identify the FA process $\lambda :[0,\infty )\ \times
\mathcal{P}(\mathbb{R}^{d})\times \Omega \rightarrow \mathcal{P}(\mathbb{R}%
^{d})$ by the formula

\begin{eqnarray}
\pi _{t}^{\mu }\left( \varphi \right) &=&\lambda (t,\omega )\mu \left(
\varphi \right)  \notag \\
&:&=\frac{\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \int \varphi \left(
S^\circ_t x+a\right) \pi _{t}^{0}\left( \mathrm{d}a\right) \mu (\mathrm{d}x)}{%
\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \mu (\mathrm{d}x)}  \notag \\
&=&\frac{1}{C^{\pi _{t}^{\mu }}}\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right)
\int \varphi \left( S^\circ_tx+a\right) \pi _{t}^{0}\left(
\mathrm{d}a\right) \theta _{\mu }\left( x\right) \pi _{0}(\mathrm{d}x),  \label{FAlin}
\end{eqnarray}

where $C^{\pi _{t}^{\mu }}$ is the normalization constant

\begin{equation*}
C^{\pi _{t}^{\mu }}:=\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \theta
_{\mu }\left( x\right) \pi _{0}(\mathrm{d}x)=\pi _{0}\left( \Xi _{t}\theta _{\mu
}\right) =\mathbb{E}\left[ \theta
_{\mu } \left( X_{0}\right) |\mathcal{Y}_{t}%
\right] \pi _{0}\left( \Xi _{t}\right) .
\end{equation*}

To obtain the representation (<ref>), we used the fact that $\mu $ is
absolutely continuous with respect to $\pi _{0}$ and that $\theta _{\mu }$ is
the density of $\mu $ with respect to $\pi _{0}$.

*Convergence of the first moments. Observe that

\begin{eqnarray}
\hat{\pi}_{t}^{\mu } &=&\frac{1}{C^{\pi _{t}^{\mu }}}\int_{\mathbb{R}%
^{d}}\Xi \left( x\right) \int (S^\circ_tx+a)\pi
_{t}^{0}\left( \mathrm{d}a\right) \theta _{\mu }\left( x\right) \pi _{0}(\mathrm{d}x)  \notag
\\
&=&\hat{x}_{t}^{0}+\frac{S^\circ_t \int_{\mathbb{R}%
^{d}}x\Xi _{t}\left( x\right) \theta _{\mu }\left( x\right) \pi _{0}(\mathrm{d}x)}{%
\int_{\mathbb{R}^{d}}\Xi _{t}\left( x\right) \theta _{\mu }\left( x\right)
\pi _{0}(\mathrm{d}x)}=\hat{x}_{t}^{0}+\frac{S^\circ_t \left(
\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) X_{0}|\mathcal{Y}_{t}%
\right] \pi _{0}\left( \Xi _{t}\right) \right) }{\left( \mathbb{E}\left[
\theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}\right] \pi _{0}\left( \Xi
_{t}\right) \right) }.  \label{meanhat}
\end{eqnarray}

It follows that

\begin{equation*}
\left\vert \hat{\pi}_{t}^{\mu }-\hat{x}_{t}^{0}\right\vert \leq \frac{\left\vert
\left\vert S^\circ_t\right\vert \right\vert \left\vert \left\vert \mathbb{%
E}\left[ \theta _{\mu }\left( X_{0}\right) X_{0}|\mathcal{Y}_{t}\right]
\right\vert \right\vert }{\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right)
|\mathcal{Y}_{t}\right] }.
\end{equation*}

Since the processes $t\rightarrow $ $\mathbb{E}\left[ \theta _{\mu }\left(
X_{0}\right) X_{0}|\mathcal{Y}_{t}\right] $ and $t\rightarrow $ $\mathbb{E}%
\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}\right] $ converge
and the second limit is positive, we have, by using (<ref>),

\begin{equation}\label{com}
\lim_{t\rightarrow \infty }\left\vert \hat{\pi}_{t}^{\mu }-\hat{x}%
_{t}^{0}\right\vert =\lim_{t\rightarrow \infty }\left\vert \left\vert
S^\circ_t\right\vert \right\vert \frac{\lim_{t\rightarrow \infty
}\left\vert \left\vert \mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right)
X_{0}|\mathcal{Y}_{t}\right] \right\vert \right\vert }{\lim_{t\rightarrow
\infty }\left\vert \left\vert \mathbb{E}\left[ \theta _{\mu }\left(
X_{0}\right)|\mathcal{Y}_{t}\right] \right\vert \right\vert }=0.
\end{equation}

Again, similar to Remark <ref>, the convergence in (<ref>) is exponentially fast. Since $\lim_{t\rightarrow \infty }e^{-\bar{\epsilon}t}\left\vert \pi
_{t}^{\mu }\right\vert =\lim_{t\rightarrow \infty }e^{-\bar{\epsilon}%
t}\left\vert \hat{x}_{t}^{0}\right\vert =0$, for any $\bar{\epsilon}>0$, it
follows that, for any $i,j=1,...,d$,

\begin{equation*}
\lim_{t\rightarrow \infty }\left\vert \left( \pi _{t}^{\mu }\right)
^{i}\left( \pi _{t}^{\mu }\right) ^{j}-\left( \hat{x}_{t}^{0}\right)
^{i}\left( \hat{x}_{t}^{0}\right) ^{j}\right\vert =0.
\end{equation*}

The last limit is used in the proof of the convergence of the covariance matrix below.

*Convergence of the covariance matrix. Choose

\begin{equation*}
\varphi ^{ij}:\mathbb{R}^{d}\rightarrow \mathbb{R},~~~\varphi ^{ij}\left(
x^{1},..,x^{d}\right) :=x^{i}x^{j}-\left( \pi _{t}^{\mu }\right) ^{i}\left(
\pi _{t}^{\mu }\right) ^{j}
\end{equation*}

We have that

\begin{eqnarray*}
( P_{\pi _{t}^{\mu }}) ^{ij}-\left( P_{t}^{0}\right) ^{ij} &=&\pi
_{t}^{\mu }\left( \varphi ^{ij}\right) -\left( P_{t}^{0}\right) ^{ij} \\
&=&\frac{1}{C^{\pi _{t}^{\mu }}}\int_{\mathbb{R}^{d}}\Xi \left( x\right)
\int \varphi ^{ij}\left( \left( S^\circ_t\right) x+a\right) \pi
_{t}^{0}\left( da\right) \theta _{\mu }\left( x\right) \pi _{0}(dx)-\left(
P_{t}^{0}\right) ^{ij} \\
&=&\sum_{k,l}\frac{\left( S^\circ_t\right) ^{il}\left(
S^\circ_t\right) ^{jk}\left( \mathbb{E}_{\mu }\left[ \theta _{\mu }\left(
X_{0}\right) X_{0}^{l}X_{0}^{k}|\mathcal{Y}_{t}\right] \right) }{\mathbb{E}%
\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}\right] } \\
&&+\sum_{l}\frac{\left( S^\circ_t\right) ^{il}\left( \mathbb{E}_{\mu }%
\left[ \theta _{\mu }\left( X_{0}\right) X_{0}^{l}|\mathcal{Y}_{t}\right]
\right) }{\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}%
\right] }\left( \hat{x}_{t}^{0}\right) ^{j} \\
&&+\sum_{l}\frac{\left( S^\circ_t\right) ^{jl}\left( \mathbb{E}_{\mu }%
\left[ \theta _{\mu }\left( X_{0}\right) X_{0}^{l}|\mathcal{Y}_{t}\right]
\right) }{\mathbb{E}\left[ \theta _{\mu }\left( X_{0}\right) |\mathcal{Y}_{t}%
\right] }\left( \hat{x}_{t}^{0}\right) ^{i} \\
&&+\left( \hat{x}_{t}^{0}\right) ^{i}\left( \hat{x}_{t}^{0}\right)
^{j}-\left( \pi _{t}^{\mu }\right) ^{i}\left( \pi _{t}^{\mu }\right) ^{j},
\end{eqnarray*}

which gives the required convergence to 0.

*Convergence for bounded uniformly continuous test functions. We
have that

\begin{equation*}
\left\vert \pi _{t}^{\mu }\left( \varphi _{t}\right) -\pi _{t}^{0}\left(
\varphi _{t}\right) \right\vert \leq \frac{1}{C^{\pi _{t}^{\mu }}}\int_{%
\mathbb{R}^{d}}\Xi \left( x\right) \int \left\vert \varphi _{t}\left( \left(
S^\circ_t\right) x+a\right) -\varphi _{t}\left( a\right) \right\vert \pi
_{t}^{0}\left( da\right) \theta _{\mu }\left( x\right) \pi _{0}(dx).
\end{equation*}

Decomposing next the integral in the numerator into the sum of the
integral over the region $\left\vert ( S^\circ_t) x\right\vert
<\delta $ and the integral over the region $\left\vert \left(
S^\circ_t\right) x\right\vert \geq \delta $ yields

\begin{eqnarray*}
\left\vert \pi _{t}^{\mu }\left( \varphi _{t}\right) -\pi _{t}^{0}\left(
\varphi _{t}\right) \right\vert  &\leq &\sup_{\left\vert y-y^{\prime
}\right\vert <\delta }\left\vert \varphi (y)-\varphi \left( y^{\prime
}\right) \right\vert +2\Vert \varphi \Vert _{\infty }\frac{\left( \mathbb{E}%
_{\mu }\left[ \theta _{\mu }\left( X_{0}\right) \mathbf{1}_{\left\{
\left\vert \left( S^\circ_t\right) X_{0}\right\vert \geq \delta \right\}
}|\mathcal{Y}_{t}\right] \right) }{\mathbb{E}\left[ \theta _{\mu }\left(
X_{0}\right) |\mathcal{Y}_{t}\right] } \\
&\leq &\sup_{\left\vert y-y^{\prime }\right\vert <\delta }\left\vert \varphi
(y)-\varphi \left( y^{\prime }\right) \right\vert +2\frac{\Vert \varphi
\Vert _{\infty }}{\delta ^{2}}\frac{\left( S^\circ_t\right) \mathbb{E}%
_{\mu }\left[ \theta _{\mu }\left( X_{0}\right) \left\vert X_{0}\right\vert
^{2}|\mathcal{Y}_{t}\right] }{\mathbb{E}\left[ \theta _{\mu }\left(
X_{0}\right) |\mathcal{Y}_{t}\right] }.
\end{eqnarray*}

As above, it follows that

\begin{equation*}
\limsup_{t\rightarrow \infty }\left\vert \pi _{t}^{\mu }\left( \varphi
_{t}\right) -\pi _{t}^{0}\left( \varphi _{t}\right) \right\vert \leq
\sup_{\left\vert y-y^{\prime }\right\vert <\delta }\left\vert \varphi
(y)-\varphi \left( y^{\prime }\right) \right\vert
\end{equation*}

and then using the uniform continuity of $\varphi $, we deduce the result and
thus complete the proof.

§.§ D. A simple example

In this subsection, we provide an illustrative example of unstable dynamics and partial observations that still yield convergence of the FA process; see also refs. <cit.>. We choose a two-dimensional linear signal $\left( X^{1},X^{2}\right) $ with $%
\lambda _{1}<0$ and $\lambda _{2}>0$:

\begin{equation*}
X_{t}^{i} =x_{0}^{i}+\int_{0}^{t}\lambda _{i}X_{s}^{i}\mathrm{d}%
s+\int_{0}^{t}\sigma ^{i}\mathrm{d}V_{s}^{i} ,\ \ i=1,2,
\end{equation*}

so that the law of $\left( X^{1},X^{2}\right) $ is not stable. More precisely,
if we choose two systems $\left( X^{1},X^{2}\right) $, $\left( \tilde{X}^{1},%
\tilde{X}^{2}\right) $ starting from $\left( x_{0}^{1},x_{0}^{2}\right) $ and $%
\left( \tilde{x}_{0}^{1},\tilde{x}_{0}^{2}\right) $, respectively, then
their corresponding expected values drift away from each other. In

\begin{equation*}
\lim_{t\rightarrow \infty }|E\left[ X_{t}^{2}\right] -E\left[ X_{t}^{2}%
\right] |=\lim_{t\rightarrow \infty }e^{\lambda _{2}t}\left\vert x_{0}^{2}-%
\tilde{x}_{0}^{2}\right\vert =\infty .
\end{equation*}

As a result, the Wasserstein distance $W_2$ between $p_t$ and $\tilde p_t$ tends
to $\infty$.

In fact, one can easily show that the first component is stable, whilst
the second one is unstable. However, we can stabilize the system by observing
only the second component, namely the unstable one. For example, we can
choose a one-dimensional observation process of the form

\begin{equation*}
\end{equation*}

which will guarantee that

\begin{equation*}
\lim_{t\rightarrow \infty }d_{W}\left( \pi _{t},\tilde{\pi}_{t}\right) = 0.
\end{equation*}

One can justify this either by checking that the assumptions of Theorem <ref> are satisfied or by explicit calculation. Using the second approach, one notices that the probability
measures $\pi _{t}$ and $\tilde{\pi}_{t}$ are both Gaussian and
shows that (i) the distance between the coresponding means converges to $0$ and
(ii) that the covariance matrices of $\pi _{t}$ and $\tilde{\pi}_{t}$
coincide. The latter two matrices are given by diag$\left( q_{t}^{11},q_{t}^{22}\right) $. They
are diagonal and

\begin{eqnarray*}
\frac{\mathrm{d}}{\mathrm{d}t}\left(
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right) &=&\left(
\begin{array}{cc}
1 & 0 \\
0 & 1%
\end{array}%
\right) +\left(
\begin{array}{cc}
\lambda _{1} & 0 \\
0 & \lambda _{2}%
\end{array}%
\right) \left(
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right) +\left(
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right) \left(
\begin{array}{cc}
\lambda _{1} & 0 \\
0 & \lambda _{2}%
\end{array}%
\right) \\
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right) \left(
\begin{array}{cc}
0 & 0 \\
0 & h^{2}%
\end{array}%
\right) \left(
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right), \\
\frac{\mathrm{d}}{\mathrm{d}t}\left(
\begin{array}{cc}
q_{t}^{11} & 0 \\
0 & q_{t}^{22}%
\end{array}%
\right) &=&\left(
\begin{array}{cc}
1+\lambda _{1}p_{t}^{22} & 0 \\
0 & 1+2\lambda _{2}p_{t}^{11}-h^{2}\left( p_{t}^{11}\right) ^{2}%
\end{array}%
\right),
\end{eqnarray*}

which implies that

\begin{equation*}
\lim_{t\rightarrow \infty }\left(
\begin{array}{c}
p_{t}^{11} \\
\end{array}%
\right) = \left(
\begin{array}{c}
- {\displaystyle \frac{1}{2\lambda _{1}}} \\
{\displaystyle \frac{ \sqrt{\left( \lambda _{2}\right) ^{2}+h^{2}}-\lambda _{2}}{h^{2}}}%
\end{array}%
\right)
\end{equation*}

and the limit is valid independent of the initial condition. This result
together with the Gaussianity property of the processes $\pi _{t}$ and $%
\tilde{\pi}_{t}$ implies the convergence in Wasserstein distance.

In the reverse situation, in which we would observe the stable component and
not the unstable one, the posterior distribution would not be stable.

In the above example, one can still apply existing results to the
one-dimensional unstable system and, after coupling it to the unobserved
component, obtain the stability of the pair. This works because the
observation process depends only on the observed component. In multidimensional,
nonlinear applications with very large dimensions $d$ and $n$, though,
it would be quite difficult to do this as it may not be possible to identify
the unstable components a priori.

The result presented in this paper does not assume that one needs to identify
a system's unstable components and observe those. For example, if
the observation process has the form

\begin{equation*}
dY_{t}=h\left( X_{s}^{1}+X_{s}^{2}\right) \mathrm{d}t+\mathrm{d}W_{t},
\end{equation*}

the FA process will still stabilize the system. More precisely, it is still the case
that both $\pi _{t}$ and $\tilde{\pi}%
_{t}$ are Gaussian and that

\begin{equation*}
\lim_{t\rightarrow \infty }d_{W}\left( \pi _{t},\tilde{\pi}_{t}\right) = 0.
\end{equation*}

§.§ E. The forecast–assimilation (FA) process as a random dynamical
system (RDS)

In this subsection, we give a brief justification of the fact that the solution
of (<ref>) can be recast as an RDS evolving in the space of
probability measures $\mathcal{P}(\mathbb{R}^{d})$. This result is
related to similar  approaches in the nonlinear filtering literature; see,
for instance, Budhiraja<cit.>, Kunita<cit.>, Picard<cit.>,
Ocone and Pardoux<cit.>, etc. The cornerstone of the argument is that the solution of
equation (<ref>) can be expressed as

\begin{equation}
\pi _{t}^{\mu }=\lambda (t,\omega )\mu =\frac{1}{C^{\pi _{0}}}\rho
_{0:t}(\mu ),  \label{kallstrie}
\end{equation}

where $\rho _{t_{0},t_{1}}(\mu )$ is a two-parameter measure-valued process
defined as

\begin{equation}
\rho _{t_{0},t_{1}}(\mu )\left( \varphi \right) =\int_{\mathbb{R}^{d}}%
\mathbb{\tilde{E}}_{t_{0},x}[\tilde{Z}_{t_{0},t_{1}}\varphi (\tilde{X}%
_{t_{0},t_{1}})\mid \mathcal{Y}]\mu \left( dx\right) ,  \label{ztilde}
\end{equation}

where $\varphi $ is an arbitrary bounded Borel measurable map and $C^{\mu }$ is
the normalization constant $C^{\mu }=\rho _{t_{0},t_{1}}(\mu )(1)$.
In (<ref>), we have the following definitions:

* the process $\tilde{Z}_{t_{0},t_{1}}=\{\tilde{Z}_{t_{0},t_{1}},\ t\geq
0\}$ is defined by

\begin{equation}
\tilde{Z}_{t_{0},t_{1}}=\exp \left( \sum_{i=1}^{m}\int_{t_{0}}^{t_{1}}h^{i}(%
\tilde{X}_{t_{0},s})\,\mathrm{d}Y_{s}^{i}-\frac{1}{2}\sum_{i=1}^{m}%
\int_{t_{0}}^{t_{1}}h^{i}(\tilde{X}_{t_{0},s})^{2}\,\mathrm{d}s\right) ;
\label{ztilde2}
\end{equation}

* the process $\tilde{X}_{t_{0},(\cdot) }=\{\tilde{X}_{t_{0},t_{1}},\
t_{0},t_{1}\geq 0\}$ is a stochastic process independent of $Y$ satisfying
the signal equation (<ref>) on $[t_{0},\infty )$; and

* $\mathbb{\tilde{E}}_{t_{0},x}$ is the expectation with respect to a
probability measure $\mathbb{P}_{t_{0},x}$ under which $Y$ is a Brownian
motion independent of $\tilde{X}$ and $\tilde{X}%
_{t_{0},t_{0}}\equiv x$.

The independence of $\tilde{X}$ and $Y$ under $\mathbb{P}_{t_{0},x}$ in
formula (<ref>) enables us to show that $\rho _{t_{0},t_{1}}$ is an
RDS, which will immediately imply that $\pi _{t_{0},t_{1}}
$ is one, too. To justify this we introduce $\Theta ^{y_{(\cdot) }}(t_{0},t_{1})
$ to be the following two-parameter family of random variables

\begin{equation}
\Theta ^{y_{(\cdot) }}(t_{0},t_{1})\triangleq \exp \!\left( h(\tilde{X}%
_{t_{0},t_{1}})^{\top }y_{t_{1}}-h(\tilde{X}_{t_{0},t_{0}})^{\top
}y_{t_{0}}+I_{t_{0},t_{1}}^{y_{(\cdot) }}-\frac{1}{2}\sum_{i=1}^{m}%
\int_{t_{0}}^{t_{1}}h^{i}(\tilde{X}_{t_{0},s})^{2}\,\mathrm{d}s\right) ,
\label{eq:robust:babar}
\end{equation}

where $I_{t_{0},t_{1}}^{y_{(\cdot) }}$, is a version of the stochastic
integral $\int_{t_{0}}^{t_{1}}y_{s}^{\top }\,\mathrm{d}h(\tilde{X}_{t_{0},s})
$ and $y_{(\cdot) }$ is a continuous path, $y_{(\cdot) }\in C_{\mathbb{R}%
^{m}}[0,\infty ).$ The argument of the exponent in the definition of $\Theta
^{y_{(\cdot) }}(t_{0},t_{1})$ is recognizable as a formal integration by
parts of the argument of the exponential in (<ref>).

Let $\rho _{t_{0},t_{1}}^{y_{(\cdot) }}\left( \mu \right) $ and $\pi
_{t_{0},t_{1}}^{y_{(\cdot) }}(\mu )$ be the following two-parameter measure
valued processes,

\begin{align}
\rho _{t_{0},t_{1}}^{y_{(\cdot) }}\left( \mu \right) (\varphi ) &= \int_{%
\mathbb{R}^{d}}\mathbb{\tilde{E}}_{t_{0},x}[\varphi (\tilde{X}%
_{t_{0},t_{1}})\Theta ^{y_{(\cdot) }}(t_{0},t_{1})\mid \mathcal{Y}]\mu \left(
dx\right),   \label{eq:robust:fdefn} \\
\pi _{t_{0},t_{1}}^{y_{(\cdot) }}(\mu )(\varphi ) &= \frac{\rho
_{t_{0},t_{1}}^{y_{(\cdot) }}\left( \mu \right) (\varphi )}{\rho
_{t_{0},t_{1}}^{y_{(\cdot) }}\left( \mu \right) (1)}.  \label{eq:robust:pidefn}
\end{align}

Then $\rho _{t_{0},t_{1}}^{Y_{(\cdot) }}(\mu )$ and $\pi
_{t_{0},t_{1}}^{Y_{(\cdot) }}(\mu )$ are versions of $\rho _{t_{0},t_{1}}(\mu )
$ and $\pi _{t_{0},t_{1}}(\mu )$.

Since $\rho _{t_{0},t_{1}}^{y_{(\cdot) }},$ $\pi _{t_{0},t_{1}}^{y_{(\cdot) }}$ can
be recast as time-inhomogenous dynamical systems, we can use them as
a basis for defining $\rho _{t_{0},t_{1}}(\mu )$ and $\pi
_{t_{0},t_{1}}(\mu )$. As a result, it is indeed the case that $\rho
_{t_{0},t_{1}}(\mu )$ and $\pi _{t_{0},t_{1}}(\mu )$ can be viewed as RDSs.
Moreover, one can show that $\rho _{t_{0},\cdot }(\mu
)=\{\rho _{t_{0},t_{1}}(\mu ),\ t_{1}\geq t_{0}\}$ satisfies the evolution

\begin{equation}
\rho _{t_{0},t_{1}}(\mu )(\varphi )=\ \mu (\varphi
)+\int_{t_{0}}^{t_{1}}\rho _{t_{0},s}(\mu )(A\varphi )\,\mathrm{d}%
s+\int_{t_{0}}^{t_{1}}\rho _{t_{0},s}(\mu )(\varphi h^{\top })\mathrm{d}%
Y_{s},  \label{zakai}
\end{equation}

for any $\varphi \in \mathcal{D}(A)$ and, in particular, that

\begin{equation}
\rho _{t_{0},t_{1}}(\mu )(1)=1+\int_{t_{0}}^{t_{1}}\rho _{t_{0},s}(\mu
)(h^{\top })\mathrm{d}Y_{s}=1+\int_{t_{0}}^{t_{1}}\rho _{t_{0},s}(\mu
)\left( 1\right) \pi _{t_{0},s}(\mu )(h^{\top })\mathrm{d}Y_{s}.
\label{zakaimass}
\end{equation}

From (<ref>) and (<ref>), one deduces that the ratio ${%
\rho _{t_{0},t_{1}}(\mu )(\varphi )/\rho _{t_{0},t_{1}}(\mu )(1)}$
satisfies (<ref>) with initial condition $\pi
_{t_{0},t_{0}}(\mu )=\mu $ and, therefore, by the uniqueness of the solution
of (<ref>), we obtain that the solution of (<ref>) can indeed be recast as an RDS.

Finally, notice that $\pi _{0,t}(\pi _{0})=\pi_{t_{0},t}(\pi _{t_{0}})$ is indeed the FA process
we considerd throughout this paper. For
the particular case of $\mu =\pi _{t_{0}}$, formula (<ref>)
is known as the Kallianpur-Striebel formula <cit.> and it is deduced directly
from the definition of the conditional expectation. For arbitrary $\mu $, (
<ref>) serves as definition for the RDS, which
is then shown to be the solution of the evolution equation (<ref>) starting from $\mu $ at time $t_{0}$.

The map $\pi _{t_{0},t_{1}}^{y_{(\cdot) }}:\mathcal{P}(\mathbb{R}%
^{d})\rightarrow \mathcal{P}(\mathbb{R}^{d})$ is a continuous map when we
endow $\mathcal{P}(\mathbb{R}^{d})$ or, rather, the set of probability
measures with finite second moment), with the topology induced by the
Wasserstein metric.

It is a pleasure to thank Eviatar Bach for several useful discussions, for providing the final version of Figs. <ref>fig:Ghil1981 and <ref>fig:FA_Cycle and a draft of Appendix I , as well as for comments on the near-final manuscript. Alberto Carrassi and Pierre del Moral also read the near-final manuscript and made constructive suggestions.
M.G. acknowledges, however belatedly, the Nelder Fellowship of the Imperial College's Mathematics Department, <http://www.imperial.ac.uk/mathematics/research/opportunities/nelder-visiting-fellowships/fellows/professor-michael-ghil/> that put the two authors in closer contact in Spring 2014. Both authors are pleased to acknowledge the Institut Henri Poincaré's trimester on “The Mathematics of Climate and the Environment” in Fall 2019,
<http://www.ihp.fr/en/CEB/T3-2019>,  which supported the real start of the collaboration leading to this paper. It is particularly gratifying to thank Franco Flandoli for an extended conversation during this trimester that set the two authors on the right path for the use of the tools in Crauel & Flandoli<cit.>. The comments of two reviewers — one from the DA community, the other from the mathematical one — have further improved the paper. D.C. was  partially supported by EU project STUOD - DLV-856408. The present paper is TiPES contribution #  xy; this project has received funding from the EU Horizon 2020 research and innovation program under grant agreement No. 820970, and it helps support the work of M.G. Work on this paper has also been supported by the EIT Climate-KIC (grant no.
190733); EIT Climate-KIC is supported by the European Institute of Innovation & Technology (EIT), a body of the European Union.