# Nested smoothing algorithms for inference and tracking of heterogeneous
multi-scale state-space systems

Sara Pérez-Vieites<EMAIL_ADDRESS>Harold Molina-Bulla
<EMAIL_ADDRESS>Joaquín Míguez<EMAIL_ADDRESS>Department of
Signal Theory & Communications, Universidad Carlos III de Madrid. Avenida de
la Universidad 30, 28911 Leganés, Madrid, Spain.

###### Abstract

Multi-scale problems, where variables of interest evolve in different time-
scales and live in different state-spaces. can be found in many fields of
science. Here, we introduce a new recursive methodology for Bayesian inference
that aims at estimating the static parameters and tracking the dynamic
variables of these kind of systems. Although the proposed approach works in
rather general multi-scale systems, for clarity we analyze the case of a
heterogeneous multi-scale model with 3 time-scales (static parameters, slow
dynamic state variables and fast dynamic state variables). The proposed
scheme, based on nested filtering methodology of [26], combines three
intertwined layers of filtering techniques that approximate recursively the
joint posterior probability distribution of the parameters and both sets of
dynamic state variables given a sequence of partial and noisy observations. We
explore the use of sequential Monte Carlo schemes in the first and second
layers while we use an unscented Kalman filter to obtain a Gaussian
approximation of the posterior probability distribution of the fast variables
in the third layer. Some numerical results are presented for a stochastic two-
scale Lorenz 96 model with unknown parameters.

###### keywords:

filtering; Kalman; Monte Carlo; Bayesian inference, parameter estimation.

##  1 Introduction

Multi-scale problems can be found naturally in many fields of science, such as
biology, chemistry, fluid dynamics and material science, where processes at
different time and/or spatial scales may be described by diverse laws [33, 34,
23]. Thanks to the improvements in computational power and the need for ever
more faithful models of real-world complex systems, the interest in multi-
scale modelling techniques has increased in recent years.

Multi-scale models are aimed at providing an efficient and more accurate
representation of a complex system by coupling (sub)models that address
different features/attributes at different levels of detail. However, the
modeling of such systems is far from straightforward, not only because they
are made up of dynamical systems of different characteristics, but also
because of often intricate cross dependencies among the relevant physical
processes at work [33]. As a consequence, the problem of tracking the
evolution of a multi-scale dynamical system involves the prediction and
estimation of several sets of variables that live in different state-spaces
and evolve in different time scales. Moreover, the tracking of the variables
of interest usually has to be performed from the observation of partial and
noisy observations. Efficient recursive algorithms for this task are badly
needed.

The simplest case of a multi-scale problem consists of a system with unknown
static parameters and dynamic state variables, since the parameters may be
considered as state variables that evolve at a greater time scale. Hence, it
is a multi-scale problem with only two time scales. In general, carrying out
both parameter estimation and state tracking at once implies several practical
and theoretical difficulties. Within this framework, a few well-principled
methods have been proposed in the last few years, though. Two examples of
schemes that yield theoretically-guaranteed solutions to this problem are
sequential Monte Carlo square (SMC2) [7] or particle Markov chain Monte Carlo
(PMCMC) [3]. They aim at computing the joint posterior probability
distribution of all the unknown variables and parameters of the system, which
provides all the information to obtain both point estimates and
quantifications of the estimation error. Unfortunately, both SMC2 and PMCMC
are batch techniques. In other words, every time a new observation arrives,
the whole sequence of observations may have to to be re-processed from scratch
in order to update the estimates, leading to a quadratic increase of the
computational effort over time. As an alternative, nested particle filters
(NPFs) [8] apply the same principles as SMC2 in a recursive way. An NPF
estimates both parameters and states using a scheme of two intertwined layers
of Monte Carlo methods, one inside the other, where the first layer estimates
the parameters while the second layer tracks the state dynamical variables.
This methodology is better suited for long sequences of observations, however,
the use of Monte Carlo in both layers of filters still makes its computational
cost prohibitive in high-dimensional problems. Nested hybrid filters (NHFs)
[26] are extensions of the NPF that introduce Gaussian filtering techniques in
the second layer of the algorithm, reducing the computational cost
considerably and making the methodology more appealing for online processing.

In the last few years, other algorithms with nested, or layered, structures
(in the vein of SMC2, NPF or NHF) have been proposed in ordert to address
inference in high-dimensional, complex models. The most recent examples are
the space-time particle filter (ST-PF) [5] and the nested sequential Monte
Carlo (NSMC) [22]. Both methods are intended to outperform classical
sequential Monte Carlo (SMC) in high dimensional systems. They rely on spatial
structures within the state space (a Markov random field in [22] and an auto-
regressive structure in [5]) and, therefore, they may be useful to tackle
multiple spatial scales.

One of the most typical examples of multi-scale systems is the two-scale
Lorenz 96 model [18, 6, 4]. This is a simplified weather model which includes
different spatial and temporal scales. Specifically, it is a heterogeneous
multi-scale model, i.e. a model where the macro-scale level description needs
to be completed by the data extracted from a micro-scale level [30]. The
coupling of both scales can be handled in two different ways:

  * •

applying parameterization [4, 23, 32] or

  * •

using a macro-micro solver [34, 35].

The former is a method used to replace the contribution of the micro-scale
level of the model by a simplified process that depends only on the slow state
variables of the macro scale. Once the model is simplified to a single scale,
many other algorithms can be used to estimate the evolution of the macro-scale
process, that is usually the scale of interest. However, the latter method
aims at estimating both scales and avoids any simplification of the model.
Inference in the two-scale Lorenz 96 model has been addressed using algorithms
such as particles filters [36, 21, 16], Gaussian filters [28, 29, 24, 13] or a
combination of different methods [14, 27].

In this paper, we propose a generalization of the nested hybrid filter (NHF)
methodology aimed at performing recursive Bayesian inference for a class of
heterogeneous multi-scale state-space models [1] that can be numerically
approximated with a micro-macro solver.We analyse the case of a Lorenz 96
system that displays three time scales (static parameters, slow dynamic state
variables at the macro-scale and fast dynamic state variables at the micro-
scale), but the methodology works in the same way for more general examples
(namely, systems with $n$ scales either in time or space).

The new scheme can be described as a three-layer nested smoother that
approximates, in a recursive manner, the posterior probability distributions
of the parameters and the two sets of state variables given the sequence of
available observations. Specifically, we approximate the posterior probability
distribution of the parameters in a first layer of computation, the posterior
probability distribution of the slow state variables in a second layer, and
the posterior probability distribution of the fast state variables in a third
layer. The computations on the second layer are conditional on the candidate
parameter values generated on the first layer, while the calculations on the
third layer are conditional on the candidates drawn at the first and second
layers. The inference techniques used in each layer can vary, leading to
different computational costs and degrees of accuracy. As examples that
illustrate the methodology, we propose two methods. The first one uses SMC
algorithms in the first and second layers, intertwined with an unscented
Kalman filter (UKF) [25] in the third layer. Similarly, the second method uses
a sequential Monte Carlo (SMC) algorithm in the first layer, but incorporates
the use of ensemble Kalman filters (EnKFs) [12] and a extended Kalman filters
(EKFs) in the second and third layers of the scheme, respectively.

The rest of the paper is organized as follows. We state the problem to be
addressed in Section 2. In Section 3, we describe the optimal smoother for
multi-scale systems with static parameters and two sets of dynamic state
variables. Two specific methods derived from the general methodology are shown
in Section 4. Numerical results for the stochastic two-scale Lorenz 96 model
are shown in Section 5 and conclusions are drawn in Section 6.

##  2 Problem Statement

###  2.1 State space models

In this chapter we place our attention on state space models that result from
the analysis of physical systems that display (intertwined) dynamical features
at different time scales. To be specific, let us consider the class of
multidimensional stochastic differential equations (SDEs) that can be written
as

$\displaystyle
d\boldsymbol{x}=f_{\boldsymbol{x}}(\boldsymbol{x},\boldsymbol{\theta})d\tau+g_{\boldsymbol{x}}(\boldsymbol{z},\boldsymbol{\theta})d\tau+\boldsymbol{Q}_{x}d\boldsymbol{v},$
(1) $\displaystyle
d\boldsymbol{z}=f_{\boldsymbol{z}}(\boldsymbol{x},\boldsymbol{\theta})d\tau+g_{\boldsymbol{z}}(\boldsymbol{z},\boldsymbol{\theta})d\tau+\boldsymbol{Q}_{z}d\boldsymbol{w},$
(2)

where

  * •

$\tau$ denotes continuous time,

  * •

$\boldsymbol{x}(\tau)\in\mathbb{R}^{d_{x}}$ and
$\boldsymbol{z}(\tau)\in\mathbb{R}^{d_{z}}$ are the slow and fast states of
the system, respectively,

  * •

$f_{\boldsymbol{x}}\colon\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{x}}$,
$g_{\boldsymbol{x}}\colon\mathbb{R}^{d_{z}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{x}}$,
$f_{\boldsymbol{z}}\colon\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{z}}$
and
$g_{\boldsymbol{z}}\colon\mathbb{R}^{d_{z}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{z}}$
are (possibly nonlinear) transition functions parameterized by a fixed vector
of unknown parameters, $\boldsymbol{\theta}\in\mathbb{R}^{d_{\theta}}$,

  * •

$\boldsymbol{Q}_{x}$ and $\boldsymbol{Q}_{z}$ are known scaling matrices that
control the intensity and covariance of the stochastic perturbations,

  * •

and $\boldsymbol{v}(\tau)$ and $\boldsymbol{w}(\tau)$ are vectors of
independent standard Wiener processes with dimension $d_{x}$ and $d_{z}$,
respectively.

Equations (1)–(2) do not have closed form solutions for general nonlinear
functions $f_{\boldsymbol{x}}$, $f_{\boldsymbol{z}}$, $g_{\boldsymbol{x}}$ and
$g_{\boldsymbol{z}}$ and they have to be discretized for their numerical
integration. In order to handle the slow and fast time scales, we apply a
macro-micro solver [31, 35] that runs an Euler-Maruyama scheme for each set of
state variables, albeit with different integration steps. To be specific, we
use $\Delta_{z}$ as the integration step of $\boldsymbol{z}$ while
$\Delta_{x}\gg\Delta_{z}$ is the integration step of $\boldsymbol{x}$. Then,
we can simulate $\boldsymbol{x}$ and $\boldsymbol{z}$ using the pair of
difference equations

$\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{x}_{t-1}+\Delta_{x}(f_{\boldsymbol{x}}(\boldsymbol{x}_{t-1},\boldsymbol{\theta})+g_{\boldsymbol{x}}(\bar{\boldsymbol{z}}_{t},\boldsymbol{\theta}))+\sqrt{\Delta_{x}}\boldsymbol{Q}_{x}\boldsymbol{v}_{t},$
(3) $\displaystyle\boldsymbol{z}_{n}$
$\displaystyle=\boldsymbol{z}_{n-1}+\Delta_{z}(f_{\boldsymbol{z}}(\boldsymbol{x}_{\lfloor\frac{n-1}{h}\rfloor},\boldsymbol{\theta})+g_{\boldsymbol{z}}({\boldsymbol{z}}_{n-1},\boldsymbol{\theta}))+\sqrt{\Delta_{z}}\boldsymbol{Q}_{z}\boldsymbol{w}_{n},$
(4)

where $\boldsymbol{x}_{t}\approx\boldsymbol{x}(t\Delta_{x})$ and
$\boldsymbol{z}_{n}\approx\boldsymbol{z}(n\Delta_{z})$ are the state signals,
$t\in\mathbb{N}$ denotes discrete time in the time scale of the slow
variables, $n\in\mathbb{N}$ denotes discrete time in the fast time scale,
$h=\frac{\Delta_{x}}{\Delta_{z}}\in\mathbb{Z}^{+}$ is the number of fast steps
(in the scale of $\boldsymbol{z}$) per slow step (in the scale of
$\boldsymbol{x}$), $\boldsymbol{v}_{t}$ and $\boldsymbol{w}_{n}$ are Gaussian
random variables (r.v.s) of zero mean and covariance matrices
$\boldsymbol{I}_{d_{x}}$ and $\boldsymbol{I}_{d_{z}}$ respectively, and
$\bar{\boldsymbol{z}}_{t}$ is an average of the fast signal computed as

$\bar{\boldsymbol{z}}_{t}=\frac{1}{h}\sum_{i=h(t-1)+1}^{ht}\boldsymbol{z}_{i}.$
(5)

We assume that the available observations may be directly related to both sets
of state variables $\boldsymbol{x}_{t}$ and $\boldsymbol{z}_{n}$, but only in
the (slow) time scale of $\boldsymbol{x}$. To be specific, the $t$-th
observation is a $d_{y}$-dimensional r.v.,
$\boldsymbol{y}_{t}\in\mathbb{R}^{d_{y}}$, which we model as

$\boldsymbol{y}_{t}=l(\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})+\boldsymbol{r}_{t},$
(6)

where
$l\colon\mathbb{R}^{d_{z}}\times\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{y}}$
is a transformation that maps the states into the observation space, and
$\boldsymbol{r}_{t}$ is a zero-mean observational-noise vector with covariance
matrix $\boldsymbol{R}$.

###  2.2 Model inference

We aim at performing joint Bayesian estimation of the parameters,
$\boldsymbol{\theta}$, and all states, $\boldsymbol{x}$ and $\boldsymbol{z}$,
for the state-space model described by Eqs. (3)-(4) and (6). Typically, the
three vectors of unknowns are tightly coupled. The estimation of the fixed
parameters is necessary to track both sets of state variables and, at the same
time, tracking the slow state variables is needed for predicting the time
evolution of the fast states and vice versa.

While in many practical applications one is typically interested in filtering,
i.e., the computation of the posterior probability density function (pdf) of
$\boldsymbol{\theta}$, $\boldsymbol{x}_{t}$ and $\boldsymbol{z}_{n}$ (for
$n=ht$) given the data sequence
$\boldsymbol{y}_{1:t}=\\{\boldsymbol{y}_{1},\boldsymbol{y}_{2},\ldots,\boldsymbol{y}_{t}\\}$,
we find more convenient to tackle the smoothing pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t},\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$.
Using the chain rule, we can factorize the latter density as

$p(\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t},\boldsymbol{\theta}|\boldsymbol{y}_{1:t})=p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t}),$
(7)

where we identify the three key conditional distributions that we seek to
compute (or approximate). Each one of these pdfs can be handled in a different
layer of computation. Hence, we aim at designing a nested inference algorithm
(in the vein of [26]) with three layers. In the first layer we compute
$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$, in the second one we obtain
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$, and in the
third layer we tackle
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$.

Hereafter we describe the methodology for the optimal (yet impractical)
calculation of the posterior pdf in Eq. (7) as well as two approximate
numerical solutions that admit feasible computational implementations.

###  2.3 Notation

We use lower-case, regular-face letters, e.g., $x$, to denote scalar
quantities and bold-face, e.g., $\boldsymbol{x}$, to denote vectors. Matrices
are represented by bold-face upper-case letters, e.g., $\boldsymbol{X}$. We
make no notational difference between between deterministic quantities and
random variables (r.v.’s).

If $\boldsymbol{x}$ is a $d$-dimensional random vector on $\mathbb{R}^{d}$, we
denote its pdf with respect to the Lebesgue measure as $p(\boldsymbol{x})$.
This is an argument-wise notation, i.e., if $\boldsymbol{x}$ and
$\boldsymbol{y}$ are random vectors, we denote their pdf’s as
$p(\boldsymbol{x})$ and $p(\boldsymbol{y})$ even if they are possibly
different functions (when $\boldsymbol{x}$ and $\boldsymbol{y}$ obey different
probability laws). Similary, $p(\boldsymbol{x},\boldsymbol{y})$ denotes the
joint pdf of $\boldsymbol{x}$ and $\boldsymbol{y}$, while
$p(\boldsymbol{x},\boldsymbol{y})$ is the conditional pdf of $\boldsymbol{x}$
given $\boldsymbol{y}$. This notation, which has been broadly used in the
field if particle filtering [17, 10, 9, 11], is simple yet sufficient to
describe the methodologies in this paper.

We also resort to a more specific notation for Gaussian pdf’s. If
$\boldsymbol{x}$ is a $d$-dimensional Gaussian random vector with mean
$\bar{\boldsymbol{x}}$ and covariance matrix $\boldsymbol{C}>0$ then we can
explicitly write the pdf $p(\boldsymbol{x})$ as

${\mathcal{N}}(\boldsymbol{x}|\bar{\boldsymbol{x}},\boldsymbol{C})=\frac{1}{\left(2\pi\right)^{\frac{d}{2}}|\boldsymbol{C}|^{\frac{1}{2}}}\exp\left\\{-\frac{1}{2}\left(\boldsymbol{x}-\bar{\boldsymbol{x}}\right)^{\top}\boldsymbol{C}^{-1}\left(\boldsymbol{x}-\bar{\boldsymbol{x}}\right)\right\\},$

where the superscript ⊤ indicates transposition and $|\boldsymbol{C}|$ denotes
the determinant of matrix $\boldsymbol{C}$.

##  3 Optimal nested smoother

We introduce the optimal nested smoothing algorithm, consisting of three
layers, that computes each of the pdfs in Eq. (7). The scheme is summarized in
Fig. 1. As a result, we obtain the posterior smoothing density
$p(\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t},\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$
which, in turn, can be used to compute the optimal filtering pdf,
$p(\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$,
by marginalization if necessary. When the exact computations demanded by this
algorithm are not feasible (for general nonlinear and/or non-Gaussian
dynamical systems) it serves as a template for approximate numerical schemes,
as shown in Section 4.

###  3.1 First layer: static parameters

The aim of this layer is to compute the posterior pdf of the parameters,
$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$, recursively. We assume that the
a priori density $p(\boldsymbol{\theta})$ is known.

At time $t$, assume that $p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t-1})$ has
been calculated. When a new observation, $\boldsymbol{y}_{t}$, is obtained, we
need to compute the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ in order to
obtain the posterior pdf of $\boldsymbol{\theta}$ at time $t$ as

$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})\propto
p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t-1}).$
(8)

However, the parameter likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ cannot be
computed directly. Instead, we decompose it as

$p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})=\int
p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{x}_{0:t},$
(9)

where
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
and $p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ are
the likelihood and the predictive pdf of the state sequence
$\boldsymbol{x}_{0:t}$, respectively, conditional on the previous observations
and the parameters. These pdfs are computed in the second layer of the
algorithm.

Optimal filter$t-1$$t$$t+1$$\boldsymbol{y}_{t-1}$$\boldsymbol{y}_{t}$ 1st
layer Initialization $p(\boldsymbol{\theta})$, $p(\boldsymbol{x}_{0})$,
$p(\boldsymbol{z}_{0})$ $p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t-1})$
Likelihood $p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Posterior pdf $p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$ 2nd layer
$p(\boldsymbol{x}_{0:t-1}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Predictive pdf
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Posterior pdf
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$ 3rd layer
$p(\boldsymbol{z}_{h(t-2)+1:h(t-1)}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Predictive pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
Likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$
Posterior pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$
ab Figure 1: Schematic depiction of the optimal smoother. Each column
represents a layer of computation and the dependencies among pdfs are
indicated by arrows. The dashed arrows are used to show relations among
different layers while the solid arrows represent dependencies in the same
layer. Arrows a and b indicate that some intermediate computations are needed
to relate both pdfs.

###  3.2 Second layer: slow states

Computations in this layer are conditional on the parameter vector
$\boldsymbol{\theta}$. We seek to compute the smoothing posterior
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$ as well as
the predictive density
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ and the
likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$,
which are needed in the first layer –see Eq. (9). We assume that the prior
density $p(\boldsymbol{x}_{0})$ is known and the posterior pdf of the slow
states at time $t-1$ (conditional on $\boldsymbol{\theta}$),
$p(\boldsymbol{x}_{0:t-1}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$, is
available at time $t$.

We first seek the predictive density of $\boldsymbol{x}_{0:t}$, namely,

$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})=p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})p(\boldsymbol{x}_{0:t-1}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta}),$
(10)

which is obtained recursively from the posterior at time $t-1$,
$p(\boldsymbol{x}_{0:t-1}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$, but
requires the evaluation of the marginal density
$p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$.
The latter is not directly available. It has to be computed as an integral
with respect to (w.r.t.) the fast state variables, in particular

$\displaystyle
p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})=$
$\displaystyle\int
p(\boldsymbol{x}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})\times$
$\displaystyle\times
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-1)+1:ht}.$
(11)

The two densities in the integrand of Eq. (11), which involve the fast state
variables $\boldsymbol{z}_{h(t-1)},\ldots,\boldsymbol{z}_{ht}$, are calculated
in the third layer. Recall that $h$ is the number of discrete-time steps of
the fast states per each single time step of the slow variables (i.e., the
$\boldsymbol{z}_{n}$’s are $h$ time faster than the $\boldsymbol{x}_{t}$’s).

As for the likelihood, when $\boldsymbol{y}_{t}$ becomes available we update
the posterior density of $\boldsymbol{x}_{0:t}$ (conditional on
$\boldsymbol{\theta}$) as

$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})\propto
p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta}).$
(12)

In the equation above, the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
can be computed as an integral w.r.t. the fast state variables, specifically,

$\displaystyle
p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})=$
$\displaystyle\int
p(\boldsymbol{y}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})\times$
$\displaystyle\times\quad
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-1)+1:ht}$
$\displaystyle=$ $\displaystyle\int
p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-1)+1:ht}.$
(13)

The likelihood function
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$
can be obtained directly from the state-space model described by Eqs. (3)-(6),
while the conditional pdf of the subsequence $\boldsymbol{z}_{h(t-1)+1:ht}$
can be further decomposed as

$\displaystyle
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
$\displaystyle=$
$\displaystyle\frac{p(\boldsymbol{x}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{t-1},\boldsymbol{\theta})}{p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})}$
(14) $\displaystyle\times
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta}).$

Both the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
of Eq. (13) and the predictive density
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
of Eq. (14) are explicitly computed in the third layer.

###  3.3 Third layer: fast states

Computations on this layer are conditional on the parameter vector
$\boldsymbol{\theta}$ and the sequence of slow states $\boldsymbol{x}_{0:t}$.
In particular, we seek to compute the conditional posterior pdfs of
$\boldsymbol{z}_{h(t-1)+1:ht}$, including the predictive densities,

$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})\quad\text{and}\quad
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta}),$

as well as the updated density
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$.
We also evaluate the plain likelihood function
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$.
We assume that the prior pdf of the fast states, $p(\boldsymbol{z})$, is known
and the posterior up to time $t-1$,
$p(\boldsymbol{z}_{h(t-2)+1:h(t-1)}|\boldsymbol{x}_{t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$,
is available to enable recursive computations.

The first predictive pdf is computed recursively from the posterior up to time
$t-1$ as the integral

$\displaystyle
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
$\displaystyle=$ $\displaystyle\int
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{z}_{h(t-2)+1:h(t-1)},\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})\times$
$\displaystyle\times
p(\boldsymbol{z}_{h(t-2)+1:h(t-1)}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-2)+1:h(t-1)}$
$\displaystyle=$ $\displaystyle\int
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{z}_{h(t-1)},\boldsymbol{x}_{t-1},\boldsymbol{\theta})\times$
$\displaystyle
p(\boldsymbol{z}_{h(t-2)+1:h(t-1)}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-2)+1:h(t-1)},$
(15)

where the transition pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{z}_{h(t-1)},\boldsymbol{x}_{t-1},\boldsymbol{\theta})$
is obtained immediately by iterating Eq. (3) in the state-space model $h$
times and
$p(\boldsymbol{z}_{h(t-2)+1:h(t-1)}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})d\boldsymbol{z}_{h(t-2)+1:h(t-1)}$
is the posterior pdf of the fast states in the previous time step. Besides,
the second predictive density,
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$,
is obtained by substituting the first predictive pdf of Eq. (15) into Eq.
(14)111Note that, in Eq. (14), the density
$p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
is the normalization constant for the conditional pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$,
while
$p(\boldsymbol{x}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{t-1},\boldsymbol{\theta})$
results from the iteration of Eq. (3)..

Finally, when the observation $\boldsymbol{y}_{t}$ becomes available, we
compute the plain likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$
(from Eq. (6) in the state-space model) and then update the conditional
posterior pdf of the fast state variables, namely,

$\displaystyle
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})=$
$\displaystyle\frac{p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})p(\boldsymbol{x}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{t-1},\boldsymbol{\theta})}{p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})}\times$
$\displaystyle\times
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$

or, simply,

$\displaystyle
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})\propto$
$\displaystyle
p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})p(\boldsymbol{x}_{t}|\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{t-1},\boldsymbol{\theta})\times$
$\displaystyle\times
p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
(16)

if we skip the normalization constant that is typically not needed explicitly
for numerical implementations.

###  3.4 Outline of the optimal nested smoother

The optimal nested smoother uses each layer of computation to track a subset
of r.v.s that evolve over their own time scale, by computing the corresponding
predictive and updated pdfs (when observations are collected), as well as the
necessary likelihoods. To be specific:

  * •

The third layer tracks the fast state variables, $\boldsymbol{z}_{n}$, and
computes the predictive pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
of Eq. (15) and the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$.
They are used to track the conditional posterior density
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$
of Eq. (16).

  * •

The second layer takes the pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
and the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{z}_{ht},\boldsymbol{x}_{t},\boldsymbol{\theta})$
in order to compute the predictive pdf
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ in Eq.
(10) and the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
in Eq. (13). These are used to track the posterior pdf of the slow state,
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$, of Eq.
(12).

  * •

The first layer takes the pdfs
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ and
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
to track the posterior pdf of the parameters,
$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$, of Eq. (8).

Finally, the three conditional posterior pdfs are needed to compute the joint
smoothing density
$p(\boldsymbol{z}_{h(t-1)+1:ht},\boldsymbol{x}_{0:t},\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$
in Eq. (7).

Figure 1 is a schematic representation of the optimal smoother, which displays
each layer in a different column. Most of the pdfs that need to be computed
are included in this scheme, showing the dependencies among them with arrows.
These relations are direct except for the arrows labeled a and b, which
require one or more intermediate computations. In the case of arrow a, the
predictive pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
is used to compute
$p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
in Eq. (11), that is necessary to calculate the predictive density
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$ of the
second layer in Eq. (10). As for the arrow labeled b, the predictive density
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t-1},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
is used to compute the pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
in Eq. (14), that is used, in turn, to obtain the likelihood
$p(\boldsymbol{y}_{t}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t-1},\boldsymbol{\theta})$
in Eq. (13).

##  4 Approximate smoothing algorithms

The optimal algorithm described in Section 3 cannot be implemented exactly for
most practical models. Instead, one needs to devise suitable approximations
that can be implemented numerically in an efficient way. One possible approach
is a full-blown SMC implementation that extends the nested particle filter of
[8]. However, such a scheme with three layers of computation results in a
prohibitive computational cost. Instead, we introduce herein two different
algorithms that combine SMC and Gaussian approximations at the different
layers. The resulting algorithms can be implemented numerically in a more
efficient manner and are suitable for parallelization, which leads to very
fast runtimes.

The first method involves using SMC schemes both at the first and second
layer, together with a bank of unscented Kalman filters (UKFs) [15, 20] to
approximate (as Gaussians) the conditional densities to be computed at the
third layer. This implementation has great potential for parallelization, but
it is computationally costly nevertheless. Hence, we also introduce a second,
less demanding scheme that utilizes the same SMC scheme at the first layer but
employs EnKFs [12] at the second layer and simple EKFs [2] to approximate the
densities needed in the second and third layer, respectively. A numerical
study of performance is carried out in Section 5 for both implementations.

###  4.1 First scheme

We introduce a numerical approximate smoother where the probability measures
$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})d\boldsymbol{\theta}$ and
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})d\boldsymbol{x}_{0:t}$
are approximated using SMC while we replace the conditional smoothing pdf
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$
by a sequence of Gaussian approximations of the densities
$p(\boldsymbol{z}_{n}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$,
for $n=h(t-1)+1,\ldots,ht$, computed using a bank of UKFs.

#### First layer

Algorithm 1 describes the first layer of the nested smoother, which aims at
the approximation of the posterior distribution of the parameters. It receives
as inputs the prior pdfs of the parameters, $p(\boldsymbol{\theta})$, and the
two subsets of state variables, $p(\boldsymbol{x}_{0})$ and
$p(\boldsymbol{z}_{0})$. In the initialization step, they are used to generate
starting Monte Carlo particles (for the SMC schemes) and sigma-points (for the
UKFs) needed at each layer. Specifically, we generate $N$ parameter samples
$\\{\boldsymbol{\theta}_{0}^{i}\\}_{1\leq i\leq N}$, $J$ slow state particles
per each parameter sample, $\\{\boldsymbol{x}_{0}^{i,j}\\}_{1\leq j\leq J}$,
and $L$ sigma-points of the fast state per each slow state sample,
$\\{\boldsymbol{z}_{0}^{i,j,l}\\}_{0\leq l\leq L-1}$, to obtain a set of the
form
$\\{\boldsymbol{\theta}_{0}^{i},\\{\boldsymbol{x}_{0}^{i,j},\\{\boldsymbol{z}_{0}^{i,j,l}\\}_{0\leq
l\leq L-1}\\}_{1\leq j\leq J}\\}_{1\leq i\leq N}$. All particles are
independent at time $t=n=0$, provided the priors are independent.

Additionally, a Markov kernel
$\kappa_{N}^{\boldsymbol{\theta}^{\prime}}(d\boldsymbol{\theta})$ is needed
for the jittering of parameter samples [8], i.e., to draw a new set of
particles, $\\{\bar{\boldsymbol{\theta}}_{t}^{i}\\}_{1\leq i\leq N}$, at each
discrete-time step. This is needed to preserve the diversity of the particles,
otherwise after a few resampling steps the parameter particles would be
reduced to just a few distinct values and the filter would collapse.

At every time step $t$ (in the slow time scale), we compute the approximate
likelihood for each particle $\bar{\boldsymbol{\theta}}_{t}^{i}$, namely

$\hat{p}^{J,L}(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})\approx
p(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i}),$

in order to obtain the non-normalized weights $\\{\tilde{v}_{t}^{i}\\}_{1\leq
i\leq N}$. The superscripts $J$ and $L$ indicate the dependence of the
approximation on the number of particles generated for the second layer ($J$)
and the number of sigma-points employed by the UKFs in the third layer ($L$).
The states
$\\{\boldsymbol{x}_{t-1}^{i,j},\\{\boldsymbol{z}_{h(t-1)}^{i,j,l}\\}_{0\leq
l\leq L-1}\\}_{1\leq j\leq J}$ are propagated to time $t$ in the nested layers
of filtering in step 2b. Finally, we normalize the weights in order to
resample not only the parameter particles $\bar{\boldsymbol{\theta}}_{t}^{i}$,
but also their associated sets of state variables.

###### Algorithm 1

SMC approximation of $p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$ in the
first method

Inputs

  * -

Prior distributions $p(\boldsymbol{\theta})$, $p(\boldsymbol{x}_{0})$ and
$p(\boldsymbol{z}_{0})$.

  * -

A Markov kernel
$\kappa_{N}^{\boldsymbol{\theta}^{\prime}}(d\boldsymbol{\theta})$ which, given
$\boldsymbol{\theta}^{\prime}$, generates jittered parameters
$\boldsymbol{\theta}\in\mathbb{R}^{d_{\theta}}$.

Initialization: this is a joint initialization for all three layers.

  * -

Draw $N$ i.i.d. sample $\boldsymbol{\theta}_{0}^{i}$, $i=1,\ldots,N$ from the
prior distribution $p(\boldsymbol{\theta})$.

  * -

Draw $J$ i.i.d. samples $\boldsymbol{x}_{0}^{i,j}$, $i=1,\ldots,N$,
$j=1,\ldots,J$, from the prior distribution $p(\boldsymbol{x}_{0})$.

  * -

Compute $L=2d_{z}+1$ sigma-points, $\boldsymbol{z}_{0}^{i,j,l}$, with their
respective weights, $\lambda_{0}^{i,j,l}$, $i=1,\ldots,N$, $j=1,\ldots,J$,
$l=0,\ldots,L-1$, from the prior distribution
$p(\boldsymbol{z}_{0}|\hat{\boldsymbol{z}}_{0},\boldsymbol{C}_{0}(\boldsymbol{z}))$
as

$\displaystyle\boldsymbol{z}_{0}^{i,j,0}$
$\displaystyle=\hat{\boldsymbol{z}}_{0},\quad$
$\displaystyle\lambda_{0}^{i,j,0}$ $\displaystyle=\frac{1}{1+d_{z}},$
$\displaystyle\boldsymbol{z}_{0}^{i,j,l}$
$\displaystyle=\hat{\boldsymbol{z}}_{0}+\boldsymbol{S}_{l},\quad$
$\displaystyle\lambda_{0}^{i,j,l}$
$\displaystyle=\frac{1-\lambda_{0}^{i,j,0}}{2d_{z}},$
$\displaystyle\boldsymbol{z}_{0}^{i,j,l+d_{z}}$
$\displaystyle=\hat{\boldsymbol{z}}_{0}-\boldsymbol{S}_{l},\quad$
$\displaystyle\lambda_{0}^{i,j,l+d_{z}}$
$\displaystyle=\frac{1-\lambda_{0}^{i,j,0}}{2d_{z}},$

for $l=1,\ldots,d_{z}$, where $\boldsymbol{S}_{l}$ is the $l$-th row or column
of the matrix square root of
$\frac{d_{z}}{1-\lambda_{0}^{i,j,0}}\boldsymbol{C}_{0}(\boldsymbol{z})$.

Procedure For $t\geq 0$:

  1. 1.

Draw $N$ i.i.d samples $\bar{\boldsymbol{\theta}}_{t}^{i}$ from
$\kappa_{N}^{\boldsymbol{\theta}_{t-1}^{i}}(d\boldsymbol{\theta})$.

  2. 2.

For $i=1,\ldots,N$:

    1. (a)

Compute

$\tilde{v}_{t}^{i}=\hat{p}^{J,L}(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i}),$
(17)

where the approximate likelihood is evaluated at layer 2.

    2. (b)

Obtain new particles
$\\{\boldsymbol{x}_{t}^{i,j},\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq l\leq
L-1}\\}_{1\leq j\leq J}$ at time $t$ (from layers 2 and 3).

    3. (c)

Normalize the weights

$\quad v_{t}^{i}=\frac{\tilde{v}_{t}^{i}}{\sum_{i=1}^{N}\tilde{v}_{t}^{i}}.$
(18)

  3. 3.

Resample: set for each $m=1,\ldots,N$ and with probability $v_{t}^{i}$

$\\{\boldsymbol{\theta}_{t}^{m},\\{\boldsymbol{x}_{t}^{(m,j)},\\{{\boldsymbol{z}}_{ht}^{m,j,l}\\}_{0\leq
l\leq L-1}\\}_{1\leq j\leq
J}\\}=\\{\bar{\boldsymbol{\theta}}_{t}^{i},\\{\boldsymbol{x}_{t}^{(i,j)},\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq
l\leq L-1}\\}_{1\leq j\leq J}\\}.$ (19)

Outputs:
$\\{\boldsymbol{\theta}_{t}^{i},\\{\boldsymbol{x}_{t}^{(i,j)},\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}\\}_{1\leq
j\leq J}\\}_{1\leq i\leq N}$.

#### Second layer

Algorithm 2 describes the implementation of a bank of conditional SMC schemes
in the second layer of the multi-scale nested smoother, one for each particle
$\bar{\boldsymbol{\theta}}_{t}^{i}$, $i=1,\ldots,N$. In this second layer we
approximate the posterior distribution with density
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})$.
The procedure is similar to the one in Algorithm 1, starting with the
generation of the particles $\bar{\boldsymbol{x}}_{t}^{i,j}$, $j=1,\ldots,J$,
the computation of the approximate likelihood

$\hat{p}^{L}(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})\approx
p(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$

and the non-normalized weights $\\{\tilde{u}_{t}^{i,j}\\}_{1\leq j\leq J}$ in
step 1a. By averaging the latter weights we can obtain $\tilde{v}_{t}^{i}$ for
its use in the first layer222The average
$\tilde{v}_{t}^{i}=\frac{1}{J}\sum_{j=1}^{J}\tilde{u}_{t}^{i,j}$ is an
approximation of the integral of Eq. (9).. After propagating the fast state
variables in the third layer (as described below), one can resample the set
$\\{\boldsymbol{x}_{0:t}^{i,j},\\{\boldsymbol{z}_{ht}^{i,j,l}\\}_{0\leq l\leq
L-1}\\}_{1\leq j\leq J}$ using the normalized weights
$\\{u_{t}^{i,j}\\}_{j=1}^{J}$ obtained in step 25.

#### Third layer

Algorithm 3 outlines the implementation of a bank of UKFs [15] conditional on
each parameter sample $\bar{\boldsymbol{\theta}}_{t}^{i}$ and the set of slow
states $\\{\bar{\boldsymbol{x}}_{t}^{i,j}\\}\cup\boldsymbol{x}_{0:t-1}^{i,j}$.
If we follow the template of the optimal smoother, then we should seek an
approximation of the density
$p(\boldsymbol{z}_{{h(t-1)+1:ht}}|\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$.
However, performing this calculation with a UKF-like scheme implies that the
dimension of the filter should be $d_{z}\times h$, in order to include the
whole subsequence of states $\boldsymbol{z}_{{h(t-1)+1:ht}}$. Such approach
would demand $2d_{z}h+1$ sigma-points for each conditional UKF algorithm, and
the computation of $NJL$ covariance matrices with dimension $2d_{z}h\times
2d_{z}h$ each, which is impractical even for moderate $d_{z}$ and $h$. For
simplicity, in order to avoid operations with large matrices, we choose to
compute Gaussian approximations of the marginal predictive densities
$p(\boldsymbol{z}_{q}|\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i}),$
for $q=h(t-1)+1,\ldots,ht$, and then use these marginals to estimate the
average of the fast states $\bar{\boldsymbol{z}}_{t}$ which is necessary in
the micro-macro solver of Eq. (3). The complete procedure is outlined in
Algorithm 3, with further details below.

Step 1 of Algorithm 3 generates new sigma-points in the space of fast states,
$\tilde{\boldsymbol{z}}_{q}^{i,j,l}$, for $q=h(t-1)+1,\ldots,ht$ and
$l=0,\ldots,L-1$, conditional on the parameters
$\bar{\boldsymbol{\theta}}_{t}^{i}$ and slow variables
$\boldsymbol{x}_{t-1}^{i,j}$. At each time $q$, we compute a predictive mean
and a covariance matrix as

$\displaystyle\check{\boldsymbol{z}}_{q}^{i,j}$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{q-1}^{i,j,l}\tilde{\boldsymbol{z}}_{q}^{i,j,l}\quad\text{and}$
(20) $\displaystyle\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{q-1}^{i,j,l}(\tilde{\boldsymbol{z}}_{q}^{i,j,l}-\check{\boldsymbol{z}}_{q}^{i,j})(\tilde{\boldsymbol{z}}_{q}^{i,j,l}-\check{\boldsymbol{z}}_{q}^{i,j})^{\top}+\Delta_{z}\boldsymbol{Q}_{z},$
(21)

where the $\lambda_{q-1}^{i,j,l}$’s are the weights333These weights are
deterministic and can be computed a priori in different ways. See [20] for a
survey of methods. of the sigma points $\breve{\boldsymbol{z}}_{q-1}^{i,j,l}$
and $\Delta_{z}\boldsymbol{Q}_{z}$ is the covariance matrix of the noise in
Eq. (4). The mean in Eq. (20) and the covariance in Eq. (21) yield the
approximation

$\mathcal{N}(\boldsymbol{z}_{q}|\check{\boldsymbol{z}}_{q}^{i,j},\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z}))\approx
p(\boldsymbol{z}_{q}|\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$
(22)

and we compute a new weighted set of sigma-points
$\\{\breve{\boldsymbol{z}}_{q}^{i,j,l},\lambda_{q}^{i,j,l}\\}$ to represent
the Gaussian density in Eq. (22).

###### Algorithm 2

SMC approximation of
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$

Inputs

  * -

Known parameter $\bar{\boldsymbol{\theta}}_{t}^{i}$ and known initial states,
$\boldsymbol{x}_{t-1}^{i,j}$ and $\boldsymbol{z}_{h(t-1)}^{i,j,l}$, for
$j=1,\ldots,J$ and $l=0,\ldots,L-1$.

Procedure For $t\geq 0$:

  1. 1.

For $j=1,\ldots,J$:

    1. (a)

Compute

$\displaystyle\tilde{u}_{t}^{i,j}$ $\displaystyle=$
$\displaystyle\hat{p}^{L}(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i}),$
(23) $\displaystyle\tilde{v}_{t}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{J}\sum_{j=1}^{J}\tilde{u}_{t}^{i,j},$ (24)

where the new particle $\bar{\boldsymbol{x}}_{t}^{i,j}$ is generated, and the
approximate likelihood is evaluated, at layer 3.

    2. (b)

Obtain new particles $\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq l\leq L-1}$
at time $t$, from layer 3.

    3. (c)

Normalize the weights

$u_{t}^{i,j}=\frac{\tilde{u}_{t}^{i,j}}{\sum_{j=1}^{J}\tilde{u}_{t}^{i,j}}.$
(25)

  2. 2.

Resample: set

$\\{\boldsymbol{x}_{t}^{i,m},\\{{\boldsymbol{z}}_{ht}^{i,m,l}\\}_{0\leq l\leq
L-1}\\}=\\{\bar{\boldsymbol{x}}_{t}^{i,j},\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq
l\leq L-1}\\}$ (26)

with probability $u_{t}^{i,j}$ for each $m=1,\ldots,J$.

Outputs:
$\\{\boldsymbol{x}_{t}^{i,j},\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq l\leq
L-1}\\}_{1\leq j\leq J}$ and $\tilde{v}_{t}^{i}$.

###### Algorithm 3

UKF approximation of
$p(\boldsymbol{z}_{h(t-1)+1:ht}|\boldsymbol{x}_{0:t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$

Inputs

  * -

Integration steps $\Delta_{x}$, $\Delta_{z}$ and time scale ratio
$h=\frac{\Delta_{x}}{\Delta_{z}}\in\mathbb{Z}^{+}$.

  * -

Known parameter vector $\bar{\boldsymbol{\theta}}_{t}^{i}$ and initial slow
state $\boldsymbol{x}_{t-1}^{i,j}$. Weighted sigma-points for the fast state
at time $h(t-1)$, denoted
$\\{\boldsymbol{z}_{h(t-1)}^{i,j,l},\lambda_{h(t-1)}^{i,j,l}\\}_{l=0}^{L-1}$.

Procedure For $t>0$:

  1. 1.

Set
$\breve{\boldsymbol{z}}_{h(t-1)}^{i,j,l}={\boldsymbol{z}}_{h(t-1)}^{i,j,l}$.
For $l=0,\ldots,L-1$ and for $q=h(t-1)+1,...,ht$:

    1. (a)

Integrate with step $\Delta_{z}$

$\tilde{\boldsymbol{z}}_{q}^{i,j,l}=\breve{\boldsymbol{z}}_{q-1}^{i,j,l}+\Delta_{z}(f_{\boldsymbol{z}}(\breve{\boldsymbol{z}}_{q-1}^{i,j,l},\bar{\boldsymbol{\theta}}_{t}^{i})+g_{\boldsymbol{z}}(\boldsymbol{x}_{t-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})),$
(27)

and compute the predictive mean, $\check{\boldsymbol{z}}_{q}^{i,j}$, and the
predictive covariance matrix,
$\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})$, using Eqs. (20) and (21).

    2. (b)

Approximate the predictive pdf of the fast states as Gaussian density,

$p(\boldsymbol{z}_{q}|\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})\approx\mathcal{N}(\boldsymbol{z}_{q}|\check{\boldsymbol{z}}_{q}^{i,j},\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})).$

Represent this Gaussian distribution by a set of weighted sigma-points denoted
$\\{\breve{\boldsymbol{z}}_{q}^{i,j,l},\lambda_{q}^{i,j,l}\\}_{l=0}^{L-1}$.

  2. 2.

In the space of the slow state variables:

    1. (a)

For $l=0,\ldots,L-1$, project the sigma-points
$\breve{\boldsymbol{z}}_{h(t-1)+1:ht}^{i,j,l}$ to obtain sigma-points in the
space of the slow states,

$\tilde{\boldsymbol{x}}_{t}^{i,j,l}=\boldsymbol{x}_{t-1}^{i,j}+\Delta_{x}(f_{\boldsymbol{x}}(\boldsymbol{x}_{t-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})+g_{\boldsymbol{x}}(\bar{\boldsymbol{z}}_{t}^{i,j,l},\bar{\boldsymbol{\theta}}_{t}^{i})),$
(28)

where
$\bar{\boldsymbol{z}}_{t}^{i,j,l}=\frac{1}{h}\sum_{q=h(t-1)+1}^{ht}{\breve{\boldsymbol{z}}_{q}^{i,j,l}}$.
Then, compute a mean vector $\check{\boldsymbol{x}}_{t}^{i,j}$ and a
covariance matrix $\check{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{x})$ using
Eqs. (34) and (35).

    2. (b)

Sample
$\bar{\boldsymbol{x}}^{i,j}_{t}\sim\mathcal{N}(\boldsymbol{x}_{t}|\check{\boldsymbol{x}}_{t}^{i,j},\check{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{x}))$.

  3. 3.

Once we collect a new observation $\boldsymbol{y}_{t}$,

    1. (a)

For $l=0,\ldots,L-1$, project the sigma-points
$\breve{\boldsymbol{z}}_{ht}^{i,j,l}$ and the new sample
$\bar{\boldsymbol{x}}^{i,j}_{t}$ into the observation space,

$\tilde{\boldsymbol{y}}_{t}^{i,j,l}=l(\breve{\boldsymbol{z}}_{ht}^{i,j,l},\bar{\boldsymbol{x}}^{i,j}_{t},\bar{\boldsymbol{\theta}}_{t}^{i}),$
(29)

then compute the mean vector $\hat{\boldsymbol{y}}_{ht}^{i,j}$ and the
covariance matrix ${\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{y})$ using Eqs. (36)
and (37).

    2. (b)

Compute
$\tilde{w}_{t}^{i,j,l}=p(\boldsymbol{y}_{t}|\breve{\boldsymbol{z}}_{ht}^{i,j,l},\bar{\boldsymbol{x}}_{t}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})p(\bar{\boldsymbol{x}}^{i,j}_{t}|\breve{\boldsymbol{z}}_{h(t-1)+1:ht}^{i,j,l},{\boldsymbol{x}}_{t-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})$
and the weights for the second layer

$\tilde{u}_{t}^{i,j}=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}\tilde{w}_{t}^{i,j,l}.$
(30)

  4. 4.

Update the mean and the covariance matrix of the fast variables

$\displaystyle\boldsymbol{K}_{t}(\boldsymbol{z})$
$\displaystyle={\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z},\boldsymbol{y})\big{(}{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{y})\big{)}^{-1},$
(31) $\displaystyle\hat{\boldsymbol{z}}_{ht}^{i,j}$
$\displaystyle=\check{\boldsymbol{z}}_{ht}^{i,j}+\boldsymbol{K}_{t}(\boldsymbol{z})\big{(}\boldsymbol{y}_{t}-\hat{\boldsymbol{y}}_{t}^{i,j}\big{)}\quad\text{and}$
(32) $\displaystyle\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})$
$\displaystyle=\check{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})+\boldsymbol{K}_{t}(\boldsymbol{z})\boldsymbol{C}_{t}^{i,j}(\boldsymbol{y})\big{(}\boldsymbol{K}_{t}(\boldsymbol{z})\big{)}^{\top},$
(33)

where ${\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z},\boldsymbol{y})$ is the
cross-covariance matrix computed in Eq. (38).

  5. 5.

From the new Gaussian pdf
$\mathcal{N}(\boldsymbol{z}_{ht}|\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z}))$,
generate $L=2d_{z}+1$ sigma-points and weights
$\\{\boldsymbol{z}_{ht}^{i,j,l},\lambda_{ht}^{i,j,l}\\}_{0\leq l\leq L-1}$.

Outputs: $\\{{\boldsymbol{z}}_{ht}^{i,j,l}\\}_{0\leq l\leq L-1}$,
$\bar{\boldsymbol{x}}_{t}^{i,j}$ and $\tilde{u}_{t}^{i,j}$.

In step 2 of Algorithm 3 we use the sigma-points at time $q=ht$ to generate
new particles for the slow states at time $t$. Specifically, we project the
$\breve{\boldsymbol{z}}_{ht}^{i,j,l}$’s through the state equation of the slow
state variables to obtain sigma-points in the space of the slow variables,
denoted $\tilde{\boldsymbol{x}}_{t}^{i,j,l}$. From these sigma-points, we
obtain a mean vector and a covariance matrix, respectively,

$\displaystyle\check{\boldsymbol{x}}_{t}^{i,j}$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}\tilde{\boldsymbol{x}}_{t}^{i,j,l}\quad\text{and}$
(34) $\displaystyle\check{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{x})$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}(\tilde{\boldsymbol{x}}_{t}^{i,j,l}-\check{\boldsymbol{x}}_{t}^{i,j})(\tilde{\boldsymbol{x}}_{t}^{i,j,l}-\check{\boldsymbol{x}}_{t}^{i,j})^{\top}+\Delta_{x}\boldsymbol{Q}_{x},$
(35)

where $\Delta_{x}\boldsymbol{Q}_{x}$ is the covariance matrix of the noise in
Eq. (3). Eqs. (34) and (35) yield a Gaussian approximation of the predictive
pdf of the slow states, namely,

$p(\boldsymbol{x}_{t}|\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{0:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})\approx{\mathcal{N}}(\boldsymbol{x}|\check{\boldsymbol{x}}_{t}^{i,j},\check{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{x})).$

We generate the new particle $\bar{\boldsymbol{x}}_{t}^{i,j}$ from this
Gaussian density.

In step 3 of the algorithm we propagate the sigma-points
$\breve{\boldsymbol{z}}_{ht}^{i,j,l}$ and the particle
$\bar{\boldsymbol{x}}_{t}^{i,j}$ through the observation function $l(\cdot)$
to obtain projected sigma-points (on the observation space)
$\\{\tilde{\boldsymbol{y}}_{t}^{i,j,l}\\}_{0\leq l\leq L-1}$. We use these
projected sigma-points to obtain a predictive mean and covariance matrix for
the observation $\boldsymbol{y}_{t}$, namely,

$\displaystyle\hat{\boldsymbol{y}}_{t}^{i,j}$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}\tilde{\boldsymbol{y}}_{t}^{i,j,l}\quad\text{and}$
(36) $\displaystyle{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{y})$
$\displaystyle=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}(\tilde{\boldsymbol{y}}_{t}^{i,j,l}-\hat{\boldsymbol{y}}_{t}^{i,j})(\tilde{\boldsymbol{y}}_{t}^{i,j,l}-\hat{\boldsymbol{y}}_{t}^{i,j})^{\top}+\boldsymbol{R},$
(37)

where $\boldsymbol{R}$ is the covariance matrix of the noise in the
observation equation. At this step we also compute the non-normalized
importance weight $\tilde{u}_{t}^{i,j}$ which is output to layer 2.

In step 4, we compute the Kalman gain using the observation covariance matrix
of Eq. (37) and the cross-covariance matrix

${\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z},\boldsymbol{y})=\sum_{l=0}^{L-1}\lambda_{ht}^{i,j,l}(\breve{\boldsymbol{z}}_{ht}^{i,j,l}-\check{\boldsymbol{z}}_{ht}^{i,j})(\tilde{\boldsymbol{y}}_{t}^{i,j,l}-\hat{\boldsymbol{y}}_{t}^{i,j})^{\top}.$
(38)

Then we update the mean, $\hat{\boldsymbol{z}}_{ht}^{i,j}$, and covariance
matrix, $\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})$, of the fast state
variables to obtain the approximation

$p(\boldsymbol{z}_{ht}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})\approx{\mathcal{N}}(\boldsymbol{z}_{ht}|\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})).$
(39)

Finally, in step 5 we generate new weighted sigma-points to characterize the
Gaussian pdf in Eq. (39).

###  4.2 Second scheme

The method in Section 4.1 may still have a prohibitive computational cost, as
it generates a total of $N\times J\times L$ particles in the joint space of
the parameters, the slow states and the fast states (if we count sigma-points
as deterministic particles). In this section we describe a computationally-
lighter procedure that replaces the SMC procedure in layer 2 by an EnKF [12]
and the UKF in layer 3 by a simpler EKF [2]. The complete procedure is
described below.

#### First layer

We describe the use of a SMC algorithm for the first layer of the nested
algorithm in order to approximate
$p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$. This is the same as in
Algorithm 1, except that we need to take into account the initializations
needed for layers 2 and 3. Algorithm 4 receives as inputs the prior
distributions of the parameters, $p(\boldsymbol{\theta})$, and both subsets of
state variables, $p(\boldsymbol{x}_{0})$ and $p(\boldsymbol{z}_{0})$. They are
used to generate

  * •

the initial particles from $p(\boldsymbol{\theta})$ for the SMC scheme,
denoted $\\{\boldsymbol{\theta}_{0}^{i}\\}_{i=1}^{N}$,

  * •

the samples from $p(\boldsymbol{x}_{0})$ utilized to build the ensembles
$\boldsymbol{X}_{0}^{i}$, $i=1,\ldots,N$, for the EnKFs in the second layer,

  * •

and the mean and covariance matrix of $\boldsymbol{z}$ for the EKFs in the
third layer, denoted $\boldsymbol{z}_{0}^{i,j}=\mathbb{E}[\boldsymbol{z}_{0}]$
and $\boldsymbol{C}_{0}^{i,j}(\boldsymbol{z})=\text{Cov}(\boldsymbol{z}_{0})$,
respectively (note that they are the same for all $i$ and $j$).

The rest of the procedure is the same as in Algorithm 1.

#### Second layer

In Algorithm 5 we employ an EnKF to obtain ensemble approximations of
$p(\boldsymbol{x}_{t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$
and
$p(\boldsymbol{x}_{t}|\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})$.
The ensembles are denoted $\bar{\boldsymbol{X}}_{t}^{i}$ and
$\boldsymbol{X}_{t}^{i}$, respectively, and they are used to approximate the
computations that involved the joint pdfs
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$
and
$p(\boldsymbol{x}_{0:t}|\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})$
in the optimal smoother. Note that all calculations are conditional on the
$i$-th parameter particle, $\bar{\boldsymbol{\theta}}_{t}^{i}$.

The scheme is similar to Algorithm 2. At step 1a, we retrieve the new samples
$\bar{\boldsymbol{x}}_{t}^{i,j}$ and the approximate likelihood
$\tilde{u}_{t}^{i,j}$ from layer 3, and compute the non-normalized importance
weight $\tilde{v}_{t}^{i}$ which is output to layer 1.

At steps 2 and 3 we generate the predictive ensemble for the slow states,
$\bar{\boldsymbol{X}}^{i}_{t}$, and the observations,
$\boldsymbol{Y}_{t}^{i}$, respectively. These ensembles are then used, when
the new observation $\boldsymbol{y}_{t}$ is collected, to compute an updated
ensemble $\boldsymbol{X}_{t}^{i}$ which yields non-weighted particle
approximation of the distribution with pdf
$p(\boldsymbol{x}_{t}|\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})$.
The update step of the EnKF can be implemented in different ways. Here we
follow the scheme in [19] which avoids the direct computation of the inverse
of the covariance observation matrix ($d_{y}\times d_{y}$), being better
suited for high-dimensional systems.

#### Third layer

In Algorithm 6 we describe how to use an EKF to obtain Gaussian approximations
$\mathcal{N}(\boldsymbol{z}_{q}|\check{\boldsymbol{z}}_{q},\check{\boldsymbol{C}}_{q}(\boldsymbol{z}))\approx
p(\boldsymbol{z}_{q}|\boldsymbol{x}_{t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$,
$q=h(t-1)+1,\ldots,ht$, and the updated pdf
$\mathcal{N}(\boldsymbol{z}_{ht}|\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z}))\approx
p(\boldsymbol{z}_{ht}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{t-1}^{i,j},\boldsymbol{y}_{1:t},\bar{\boldsymbol{\theta}}_{t}^{i})$.
We also generate the new slow states at time $t$, denoted
$\bar{\boldsymbol{x}}_{t}^{i,j}$, and the likelihood estimates
$\tilde{u}_{t}^{i,j}\approx
p(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$.
Note that all computations in this layer are conditional on
$\boldsymbol{x}_{t-1}^{i,j}$ and $\bar{\boldsymbol{\theta}}_{t}^{i}$.

In step 1, the algorithm propagates the mean,
$\hat{\boldsymbol{z}}_{h(t-1)}^{i,j}$, and the covariance matrix,
$\hat{\boldsymbol{C}}_{h(t-1)}^{i,j}(\boldsymbol{z})$, for
$q=h(t-1)+1,\ldots,ht$, conditional on the parameters
$\bar{\boldsymbol{\theta}}_{t}^{i}$ and slow variables
$\boldsymbol{x}_{t-1}^{i,j}$. At each time step $q$, we obtain a predictive
mean, $\check{\boldsymbol{z}}_{q}^{i,j}$, and the predictive covariance
matrix, $\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})$. The average of the
predictive means,
$\bar{\boldsymbol{z}}_{t}^{i,j}=\frac{1}{h}\sum_{q=h(t-1)+1}^{ht}\check{\boldsymbol{z}}_{q}^{i,j}$,
is then used to propagate the slow state $\boldsymbol{x}_{t-1}^{i,j}$ and
generate the new sample $\bar{\boldsymbol{x}}_{t}^{i,j}$ from the Gaussian
approximation

${\mathcal{N}}(\boldsymbol{x}_{t}|\check{\boldsymbol{x}}_{t}^{i,j},\Delta_{x}\boldsymbol{Q}_{x})\approx
p(\boldsymbol{x}_{t}|\boldsymbol{x}_{t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$

at step 2. Note that the covariance of the $\boldsymbol{z}_{q}^{i,j}$’s is
neglected for simplicity in this computation and we use just the covariance
matrix $\Delta_{x}\boldsymbol{Q}_{x}$ of the slow state in Eq. (3).

In step 3 we project the predictive mean $\check{\boldsymbol{z}}_{ht}^{i,j}$
and the sample $\bar{\boldsymbol{x}}_{t}^{i,j}$ into the observation space to
obtain the predictive observation $\hat{\boldsymbol{y}}_{t}^{i,j}$. When the
actual observation $\boldsymbol{y}_{t}$ is available we also estimate the
likelihood
$p(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$
as

$\tilde{u}_{t}^{i,j}=p(\boldsymbol{y}_{t}|\check{\boldsymbol{z}}_{ht}^{i,j},\bar{\boldsymbol{x}}_{t}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})p(\boldsymbol{x}_{t}|\check{\boldsymbol{z}}_{h(t-1)+1:ht}^{i,j},{\boldsymbol{x}}_{t-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i}).$

Note that we use the predictive means
$\check{\boldsymbol{z}}_{h(t-1):ht}^{i,j}$ for simplicity, instead of actually
integrating w.r.t. the random states $\boldsymbol{z}_{h(t-1):ht}$.

Finally, in step 4, we compute the Kalman gain using the predictive covariance
matrix of Eq. (50) and the Jacobian matrix
$\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})$. Then, we update the mean,
$\hat{\boldsymbol{z}}_{ht}^{i,j}$, and covariance matrix,
$\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})$, of the fast variables for
the next time step.

###### Algorithm 4

SMC approximation of $p(\boldsymbol{\theta}|\boldsymbol{y}_{1:t})$ in the
second method Inputs

  * -

Prior distributions $p(\boldsymbol{\theta})$, $p(\boldsymbol{x}_{0})$ and
$p(\boldsymbol{z}_{0})$.

  * -

A Markov kernel
$\kappa_{N}^{\boldsymbol{\theta}^{\prime}}(d\boldsymbol{\theta})$ which, given
$\boldsymbol{\theta}^{\prime}$, generates jittered parameters
$\boldsymbol{\theta}\in\mathbb{R}^{d_{\theta}}$.

Initialization: this is a joint initialization for all three layers.

  * -

Draw $N$ i.i.d. sample $\boldsymbol{\theta}_{0}^{i}$, $i=1,\ldots,N$ from the
prior distribution $p(\boldsymbol{\theta})$.

  * -

Draw $NJ$ i.i.d. samples $\boldsymbol{x}_{0}^{i,j}$, $i=1,\ldots,N$,
$j=1,\ldots,J$, from the prior distribution $p(\boldsymbol{x}_{0})$, and build
the ensembles $\boldsymbol{X}_{0}^{i}$, $i=1,\ldots,N$, as

$\boldsymbol{X}_{0}^{i}=[\boldsymbol{x}_{0}^{i,1},\ldots,\boldsymbol{x}_{0}^{i,J}].$
(40)

  * -

Set $\hat{\boldsymbol{z}}_{0}^{i,j}=\bar{\boldsymbol{z}}_{0}$ and
$\hat{\boldsymbol{C}}_{0}^{i,j}(\boldsymbol{z})=\boldsymbol{C}_{0}$, for
$i=1,\ldots,N$, $j=1,\ldots,J$, where $\bar{\boldsymbol{z}}_{0}$ and
$\boldsymbol{C}_{0}$ are the prior mean and prior covariance of
$\boldsymbol{z}_{0}$, respectively, obtained from the prior density
$p(\boldsymbol{z}_{0})$.

Procedure For $t\geq 0$:

  1. 1.

Draw $N$ i.i.d samples $\bar{\boldsymbol{\theta}}_{t}^{i}$ from
$\kappa_{N}^{\boldsymbol{\theta}_{t-1}^{i}}(d\boldsymbol{\theta})$.

  2. 2.

For $i=1,\ldots,N$:

    1. (a)

Retrieve

$\tilde{v}_{t}^{i}=\hat{p}^{J}(\boldsymbol{y}_{t}|\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i}),$
(41)

where the approximate likelihood is evaluated at layer 2.

    2. (b)

Obtain new particles
$\\{\hat{\boldsymbol{X}}_{t}^{i},\\{\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z})\\}_{1\leq
j\leq J}\\}$ at time $t$ (from layers 2 and 3).

    3. (c)

Normalize the weights

$\quad v_{t}^{i}=\frac{\tilde{v}_{t}^{i}}{\sum_{i=1}^{N}\tilde{v}_{t}^{i}}.$
(42)

  3. 3.

Resample: set for each $m=1,\ldots,N$

$\displaystyle\\{\boldsymbol{\theta}_{t}^{m},$
$\displaystyle\boldsymbol{X}_{t}^{m},\\{\hat{\boldsymbol{z}}_{ht}^{m,j},\hat{\boldsymbol{C}}_{t}^{m,j}(\boldsymbol{z})\\}_{j=1}^{J}\\}=\\{\bar{\boldsymbol{\theta}}_{t}^{i},\hat{\boldsymbol{X}}_{t}^{i},\\{\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z})\\}_{j=1}^{J}\\}$
(43)

with probability $v_{t}^{i}$.

Outputs:
$\\{\boldsymbol{\theta}_{t}^{i},\boldsymbol{X}_{t}^{i},\\{\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z})\\}_{j=1}^{J}\\}_{i=1}^{N}$.

###### Algorithm 5

EnKF approximation of
$p(\boldsymbol{x}_{t}|\boldsymbol{y}_{1:t},\boldsymbol{\theta})$

Inputs

  * -

Parameter vector $\bar{\boldsymbol{\theta}}_{t}^{i}$; ensemble of slow states,
$\boldsymbol{X}_{t-1}^{i}=[\boldsymbol{x}_{t-1}^{i,1},\ldots,\boldsymbol{x}_{t-1}^{i,J}]$;
mean vector $\hat{\boldsymbol{z}}_{h(t-1)}^{i,j}$ and the covariance matrix
$\hat{\boldsymbol{C}}^{i,j}_{t-1}(\boldsymbol{z})$, for $j=1,\ldots,J$.

Procedure For $t\geq 0$:

  1. 1.

For $j=1,\ldots,J$:

    1. (a)

Retrieve the new sample $\bar{\boldsymbol{x}}_{t}^{i,j}$ and the likelihood
estimate $\tilde{u}_{t}^{i,j}$ from layer 3 and compute the non-normalized
importance weight

$\displaystyle\tilde{v}_{t}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{J}\sum_{j=1}^{J}\tilde{u}_{t}^{i,j},$ (44)

    2. (b)

Obtain the new mean $\hat{\boldsymbol{z}}_{ht}^{i,j}$ and covariance matrix
$\hat{\boldsymbol{C}}^{i,j}_{t}(\boldsymbol{z})$ at time $t$, from layer 3.

  2. 2.

Compute the predictive mean $\check{\boldsymbol{x}}_{t}^{i}$ and construct the
predictive ensemble $\bar{\boldsymbol{X}}_{t}^{i}$ as

$\check{\boldsymbol{x}}_{t}^{i}=\frac{1}{J}\sum_{j=1}^{J}\bar{\boldsymbol{x}}_{t}^{i,j}\quad\text{and}\quad\bar{\boldsymbol{X}}_{t}^{i}=[\bar{\boldsymbol{x}}_{t}^{i,1},\ldots,\bar{\boldsymbol{x}}_{t}^{i,J}].$
(45)

  3. 3.

Obtain predictive observations $\tilde{\boldsymbol{y}}_{t}^{i,j}$ from layer
3, then compute the mean $\hat{\boldsymbol{y}}_{t}^{i}$ and the ensemble
$\boldsymbol{Y}_{t}^{i}$ as

$\displaystyle\hat{\boldsymbol{y}}_{t}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{J}\sum_{j=1}^{J}\tilde{\boldsymbol{y}}_{t}^{i,j}\quad\text{and}\quad{\boldsymbol{Y}}_{t}^{i}=[\tilde{\boldsymbol{y}}_{t}^{i,1},\ldots,\tilde{\boldsymbol{y}}_{t}^{i,J}],\quad\text{for
}j=1,\ldots,J.$ (46)

  4. 4.

Update the ensemble of slow variables

$\displaystyle{\boldsymbol{C}}_{t}^{i}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=$
$\displaystyle\frac{1}{J}\tilde{\boldsymbol{X}}_{t}^{i}(\tilde{\boldsymbol{Y}}_{t}^{i})^{\top},$
$\displaystyle\big{(}{\boldsymbol{C}}_{t}^{i}(\boldsymbol{y})\big{)}^{-1}$
$\displaystyle=$
$\displaystyle\boldsymbol{R}^{-1}-\boldsymbol{R}^{-1}\frac{1}{J}\tilde{\boldsymbol{Y}}_{t}^{i}\bigg{(}\boldsymbol{I}_{J}+(\tilde{\boldsymbol{Y}}_{t}^{i})^{\top}\boldsymbol{R}^{-1}\frac{1}{J}\tilde{\boldsymbol{Y}}_{t}^{i}\bigg{)}^{-1}(\tilde{\boldsymbol{Y}}_{t}^{i})^{\top}\boldsymbol{R}^{-1},$
$\displaystyle\boldsymbol{K}_{t}(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle{\boldsymbol{C}}_{t}^{i}(\boldsymbol{x},\boldsymbol{y})\big{(}{\boldsymbol{C}}_{t}^{i}(\boldsymbol{y})\big{)}^{-1},$
(47) $\displaystyle\hat{\boldsymbol{X}}_{t}^{i}$ $\displaystyle=$
$\displaystyle\bar{\boldsymbol{X}}_{t}^{i}+\boldsymbol{K}_{t}(\boldsymbol{x}){\big{(}\boldsymbol{y}_{t}\mathds{1}_{d_{y}\times
J}+\boldsymbol{T}_{t}^{i}-{\boldsymbol{Y}}_{t}^{i}\big{)}},$ (48)

where
$\tilde{\boldsymbol{X}}_{t}^{i}=\bar{\boldsymbol{X}}_{t}^{i}-\check{\boldsymbol{x}}_{t}^{i}\mathds{1}_{d_{x}\times
J}$ and
$\tilde{\boldsymbol{Y}}_{t}^{i}={\boldsymbol{Y}}_{t}^{i}-\hat{\boldsymbol{y}}_{t}^{i}\mathds{1}_{d_{y}\times
J}$, ${\boldsymbol{C}}_{t}^{i}(\boldsymbol{x},\boldsymbol{y})$ is the cross
covariance matrix, ${\boldsymbol{C}}_{t}^{i}(\boldsymbol{y})$ is the
covariance matrix of the observation, $\boldsymbol{R}$ is the covariance
matrix of the noise in the observation equation, $\mathds{1}_{a\times b}$ is a
$a\times b$ matrix of ones and
$\boldsymbol{T}_{t}^{i}=\left[\boldsymbol{r}_{t}^{1},\ldots,\boldsymbol{r}_{t}^{J}\right]$
with
$\boldsymbol{r}_{t}^{j}\sim\mathcal{N}(\boldsymbol{r}_{t}|\boldsymbol{0},\boldsymbol{R})$
is a matrix of Gaussian perturbations.

Outputs: $\hat{\boldsymbol{X}}_{t}^{i}$,
$\\{\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{t}^{i,j}(\boldsymbol{z})\\}_{j=1}^{J}$
and $\tilde{v}_{t}^{i}$.

###### Algorithm 6

EKF approximation of
$p(\boldsymbol{z}_{ht}|\boldsymbol{x}_{t},\boldsymbol{y}_{1:t},\boldsymbol{\theta})$

Inputs

  * -

Integration steps $\Delta_{x}$, $\Delta_{z}$ and time scale ratio
$h=\frac{\Delta_{x}}{\Delta_{z}}\in\mathbb{Z}^{+}$.

  * -

Parameter vector $\bar{\boldsymbol{\theta}}_{t}^{i}$; slow states
${\boldsymbol{x}}_{t-1}^{i,j}$ (i.e. the $j$-th column of
$\boldsymbol{X}_{t-1}^{i}$) and mean fast state
$\hat{\boldsymbol{z}}_{h(t-1)}^{i,j}$ with covariance matrix
$\hat{\boldsymbol{C}}_{t-1}^{i,j}(\boldsymbol{z})$.

Procedure For $t\geq 0$:

  1. 1.

Set $\check{\boldsymbol{z}}_{h(t-1)}^{i,j}={\boldsymbol{z}}_{h(t-1)}^{i,j}$
and
$\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})={\boldsymbol{C}}_{h(t-1)}^{i,j}(\boldsymbol{z})$.
For $q=h(t-1)+1,...,ht$, compute

$\displaystyle\check{\boldsymbol{z}}_{q}^{i,j}$ $\displaystyle=$
$\displaystyle\check{\boldsymbol{z}}_{q-1}^{i,j}+\Delta_{z}(f_{\boldsymbol{z}}(\check{\boldsymbol{z}}_{q-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})+g_{\boldsymbol{z}}({\boldsymbol{x}_{t-1}^{i,j}},\bar{\boldsymbol{\theta}}_{t}^{i})),$
(49) $\displaystyle\check{\boldsymbol{C}}_{q}^{i,j}(\boldsymbol{z})$
$\displaystyle=$
$\displaystyle\boldsymbol{J}_{q}^{i,j}(\boldsymbol{z})\check{\boldsymbol{C}}_{q-1}^{i,j}(\boldsymbol{z})\big{(}\boldsymbol{J}_{q}^{i,j}(\boldsymbol{z})\big{)}^{\top}+\Delta_{z}\boldsymbol{Q}_{z},$
(50)

where $\boldsymbol{J}_{q}^{i,j}(\boldsymbol{z})$ is the Jacobian matrix of the
transition function of $\boldsymbol{z}$ and $\Delta_{z}\boldsymbol{Q}_{z}$ is
the covariance matrix of the noise in Eq. (4).

  2. 2.

In the space of the slow state variables:

    1. (a)

Project the predictive means $\check{\boldsymbol{z}}_{h(t-1)+1:ht}^{i,j}$ into
the space of slow variables,

$\check{\boldsymbol{x}}_{t}^{i,j}={\boldsymbol{x}_{t-1}^{i,j}}+\Delta_{x}(f_{\boldsymbol{x}}({\boldsymbol{x}_{t-1}^{i,j}},\bar{\boldsymbol{\theta}}_{t}^{i})+g_{\boldsymbol{x}}(\bar{\boldsymbol{z}}_{t}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})),$
(51)

where
$\bar{\boldsymbol{z}}_{t}^{i,j}=\frac{1}{h}\sum_{q=h(t-1)+1}^{ht}{\check{\boldsymbol{z}}_{q}^{i,j}}$.

    2. (b)

Sample
$\bar{\boldsymbol{x}}^{i,j}_{t}\sim\mathcal{N}(\boldsymbol{x}_{t}|\check{\boldsymbol{x}}_{t}^{i,j},\Delta_{x}\boldsymbol{Q}_{x})$,
where $\Delta_{x}\boldsymbol{Q}_{x}$ is the covariance matrix of the noise in
Eq. (3).

  3. 3.

When a new observation $\boldsymbol{y}_{t}$ is collected:

    1. (a)

Project the predictive mean $\check{\boldsymbol{z}}_{ht}^{i,j}$ and the slow
state $\bar{\boldsymbol{x}}_{t}^{i,j}$ into the observation space

$\tilde{\boldsymbol{y}}_{t}^{i,j}=l(\check{\boldsymbol{z}}_{ht}^{i,j},\bar{\boldsymbol{x}}^{i,j}_{t},\bar{\boldsymbol{\theta}}_{t}^{i}).$
(52)

    2. (b)

Compute
$\tilde{u}_{t}^{i,j}=p(\boldsymbol{y}_{t}|\check{\boldsymbol{z}}_{ht}^{i,j},\bar{\boldsymbol{x}}_{t}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})p(\boldsymbol{x}_{t}|\check{\boldsymbol{z}}_{h(t-1)+1:ht}^{i,j},{\boldsymbol{x}}_{t-1}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})$.
This is an estimate of the likelihood
$p(\boldsymbol{y}_{t}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{1:t-1},\bar{\boldsymbol{\theta}}_{t}^{i})$.

  4. 4.

Update the mean and the covariance matrix of the fast variables

$\displaystyle\boldsymbol{K}_{t}(\boldsymbol{z})$ $\displaystyle=$
$\displaystyle\check{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})^{\top}\bigg{(}\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})\check{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})^{\top}+\boldsymbol{R}\bigg{)}^{-1},$
(53) $\displaystyle\hat{\boldsymbol{z}}_{ht}^{i,j}$ $\displaystyle=$
$\displaystyle\check{\boldsymbol{z}}_{ht}^{i,j}+{\boldsymbol{K}_{t}(\boldsymbol{z})\big{(}\boldsymbol{y}_{t}-\tilde{\boldsymbol{y}}_{t}^{i,j}\big{)}}\quad\text{and}$
(54) $\displaystyle\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})$
$\displaystyle=$
$\displaystyle\bigg{(}\boldsymbol{I}_{d_{z}}-\boldsymbol{K}_{t}(\boldsymbol{z})\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})\bigg{)}\check{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z}),$
(55)

where $\boldsymbol{H}_{t}^{i,j}(\boldsymbol{z})$ is the Jacobian matrix of
function
$l(\cdot,\bar{\boldsymbol{x}}_{t}^{i,j},\bar{\boldsymbol{\theta}}_{t}^{i})$
w.r.t. $\check{\boldsymbol{z}}_{ht}^{i,j}$ and $\boldsymbol{R}$ is the
covariance matrix of the noise in the observation equation. We obtain the
approximation
$\mathcal{N}(\boldsymbol{z}_{ht}|\hat{\boldsymbol{z}}_{ht}^{i,j},\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z}))\approx
p(\boldsymbol{z}_{ht}|\bar{\boldsymbol{x}}_{t}^{i,j},\boldsymbol{x}_{0:t-1}^{i,j},\boldsymbol{y}_{0:t},\bar{\boldsymbol{\theta}}_{t}^{i})$.

Outputs:
$\hat{\boldsymbol{z}}_{ht}^{i,j}$,$\hat{\boldsymbol{C}}_{ht}^{i,j}(\boldsymbol{z})$,
$\bar{\boldsymbol{x}}_{t}^{i,j}$ and $\tilde{u}_{t}^{i,j}$.

##  5 Example

###  5.1 Stochastic two-scale Lorenz 96 model

In order to illustrate the application of the methods described in Section 4,
we consider a stochastic version of the two-scale Lorenz 96 model [4], which
depends on a set of fixed parameters, a set of fast variables and a set of
slow variables. The slow variables are represented by a $d_{x}$-dimensional
vector, $\boldsymbol{x}$, while the fast variables, $\boldsymbol{z}$, are
$d_{z}$-dimensional. Let us assume there are $R$ fast variables per slow
variable, therefore $d_{z}=Rd_{x}$. The system is described, in continuous-
time $\tau$, by the SDEs

$\displaystyle{dx_{j}}$
$\displaystyle=\bigg{[}-x_{j-1}(x_{j-2}-x_{j+1})-x_{j}+F-\frac{HC}{B}\sum_{l=(j-1)R}^{Rj-1}z_{l}\bigg{]}d\tau+\sigma_{x}dv_{j},$
(56) $\displaystyle{dz_{l}}$
$\displaystyle=\bigg{[}-CBz_{l+1}(z_{l+2}-z_{l-1})-Cz_{l}+\frac{CF}{B}+\frac{HC}{B}x_{{\lfloor(l-1)/R\rfloor}}\bigg{]}d\tau+\sigma_{z}dw_{l},$
(57)

where $j=0,\ldots,d_{x}-1$, $l=0,\ldots,d_{z}-1$;
$\boldsymbol{v}=(v_{0},\ldots,v_{d_{x}-1})^{\top}$ and
$\boldsymbol{w}=(w_{0},\ldots,w_{d_{z}-1})^{\top}$ are, respectively,
$d_{x}$\- and $d_{z}$-dimensional vectors of independent standard Wiener
processes; $\sigma_{x}>0$ and $\sigma_{z}>0$ are known scale parameters and
$\boldsymbol{\theta}=(F,H,C,B)^{\top}\in\mathbb{R}$ are static model
parameters. Using a micro-macro solver [31, 35] that runs an Euler-Maruyama
scheme at each time-scale to integrate Eqs. (56)–(57), the discrete-time state
equation can be written as

$\displaystyle x_{t+1,j}$ $\displaystyle=$ $\displaystyle
x_{t,j}+\Delta_{x}(f_{\boldsymbol{x},j}(\boldsymbol{x}_{t},\boldsymbol{\theta})+g_{\boldsymbol{x},j}(\bar{\boldsymbol{z}}_{t+1},\boldsymbol{\theta}))+\sqrt{\Delta_{x}}\sigma_{x}v_{t+1,j},$
(58) $\displaystyle z_{n+1,l}$ $\displaystyle=$ $\displaystyle
z_{n,l}+\Delta_{z}(f_{\boldsymbol{z},l}(\boldsymbol{x}_{\lfloor\frac{n}{h}\rfloor},\boldsymbol{\theta})+g_{\boldsymbol{z},l}(\boldsymbol{z}_{n},\boldsymbol{\theta}))+\sqrt{\Delta_{z}}\sigma_{z}w_{n+1,l},$
(59)

where

$\boldsymbol{x}_{t}=(x_{t,0},\ldots,x_{t,d_{x}-1})^{\top}\quad\text{and}\quad{\boldsymbol{z}_{n}=(z_{n,0},\ldots,z_{n,d_{z}-1})^{\top}}$

are the discrete-time slow and fast variables, respectively;
$\bar{\boldsymbol{z}}_{t}$ is the time-average

$\bar{\boldsymbol{z}}_{t}={\frac{1}{h}\sum_{n=h(t-1)+1}^{ht}\boldsymbol{z}_{n}}$

and we denote
$\bar{\boldsymbol{z}}_{t}=(\bar{z}_{t,0},\ldots,\bar{z}_{t,d_{z}-1})^{\top}$;
the terms $v_{t,j}$ and $w_{n,l}$ are independent Gaussian variables with
identical ${\mathcal{N}}(\cdot|0,1)$ pdf for all $t$, $j$, $n$ and $l$, and
the functions

$\displaystyle f_{\boldsymbol{x},j}$ $\displaystyle:$
$\displaystyle\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{x}},$
$\displaystyle g_{\boldsymbol{x},j}$ $\displaystyle:$
$\displaystyle\mathbb{R}^{d_{z}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{x}},$
$\displaystyle f_{\boldsymbol{z},l}$ $\displaystyle:$
$\displaystyle\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{z}}\quad\text{and}$
$\displaystyle g_{\boldsymbol{z},l}$ $\displaystyle:$
$\displaystyle\mathbb{R}^{d_{z}}\times\mathbb{R}^{d_{\theta}}\to\mathbb{R}^{d_{z}}$

can be expressed as

$\displaystyle f_{\boldsymbol{x},j}(\boldsymbol{x}_{t},\boldsymbol{\theta})$
$\displaystyle=-x_{t,j-1}(x_{t,j-2}-x_{t,j+1})-x_{t,j}+F,$ $\displaystyle
g_{\boldsymbol{x},j}(\bar{\boldsymbol{z}}_{t},\boldsymbol{\theta})$
$\displaystyle=-\frac{HC}{B}\sum_{l=(j-1)R}^{Rj-1}\bar{z}_{t,l},$
$\displaystyle f_{\boldsymbol{z},l}(\boldsymbol{x}_{t},\boldsymbol{\theta})$
$\displaystyle=\frac{HC}{B}x_{t,{\lfloor(l-1)/R\rfloor}}\quad\text{and}$
$\displaystyle g_{\boldsymbol{z},l}(\boldsymbol{z}_{n},\boldsymbol{\theta})$
$\displaystyle=-CBz_{n,l+1}(z_{n,l+2}-z_{n,l-1})-Cz_{n,l}+\frac{CF}{B}.$

We assume that the observations are linear and Gaussian, namely,

$\boldsymbol{y}_{t}=\boldsymbol{A}_{t}\begin{bmatrix}\boldsymbol{x}_{t}\\\
\boldsymbol{z}_{ht}\end{bmatrix}+\boldsymbol{r}_{t},$ (60)

where $\boldsymbol{A}_{t}$ is a known $d_{y}\times(d_{x}+d_{z})$ matrix and
$\boldsymbol{r}_{t}$ is a $d_{y}$-dimensional Gaussian random vector with
known covariance matrix

$\boldsymbol{R}=\begin{bmatrix}\sigma_{y,x}^{2}\boldsymbol{I}_{d_{x}}&\boldsymbol{0}\\\
\boldsymbol{0}&\sigma_{y,z}^{2}\boldsymbol{I}_{d_{z}}\end{bmatrix},$ (61)

and $\sigma_{y,x}^{2},\sigma_{y,z}^{2}>0$ are known variances.

###  5.2 Numerical results

We have run simulations for the two-scale Lorenz 96 model of Section 5.1, with
dimensions $d_{x}=10$ and $d_{z}=50$. The time steps for the Euler-Maruyama
integrators are $\Delta_{x}=10^{-3}$ and $\Delta_{z}=10^{-4}$ continuous-time
units. We set the fixed parameters as $F=8$, $H=0.75$, $C=10$ and $B=15$. In
order to obtain the initial states $\boldsymbol{x}_{0}$ and
$\boldsymbol{z}_{0}$, we simulate a deterministic version of Eqs. (58)–(59)
($\sigma_{x}=\sigma_{z}=0$) for $20$ continuous-time units. We set the initial
states as the values of variables $\boldsymbol{x}$ and $\boldsymbol{z}$ at the
last time step of this simulation. This initialization is used in all
simulations of this computer experiment in order to generate both “ground
truth” sequences of $\boldsymbol{x}_{t}$ and $\boldsymbol{z}_{n}$ and the
associated sequences of observations $\boldsymbol{y}_{t}$. We set the matrix
$\boldsymbol{A}_{t}=\boldsymbol{I}_{d_{y}}$, for $d_{y}=d_{x}+d_{z}$.

In the experiments, we compare the performance of both methods proposed (the
first one of Section 4.1 and the second one of Section 4.2). We experiment
with different number of samples $N$ and $J$ in the first and second layers of
the former methods. Additionally, for the first method (SMC-SMC-UKF) we run
the multi-scale hybrid filter with $L=2d_{z}+1=101$ sigma-points for the UKF
in the third layer. We need to estimate $\boldsymbol{\theta}=[F,C,H,B]^{\top}$
(hence, $d_{\theta}=4$). The prior for the unknown parameters is uniform,
namely $p(\boldsymbol{\theta})=\mathcal{U}([2,20]^{2})$, while the priors used
in the filtering algorithm for both unknown state variables are Gaussian,
namely
$p(\boldsymbol{x}_{0})=\mathcal{N}(\boldsymbol{x}_{0},0.1\boldsymbol{I}_{d_{x}})$
and
$p(\boldsymbol{z}_{0})=\mathcal{N}(\boldsymbol{z}_{0},10\boldsymbol{I}_{d_{z}})$.
The noise scaling factors, $\sigma_{x}=\frac{1}{2}$,
$\sigma_{z}=\frac{1}{16}$, $\sigma_{y,x}=10^{-1}$ and $\sigma_{y,z}=10^{-3}$,
are known. The jittering kernel is
$\kappa_{N}^{\boldsymbol{\theta}^{\prime}}(d\boldsymbol{\theta})=\mathcal{N}(\boldsymbol{\theta}|\boldsymbol{\theta}^{\prime},\tilde{\sigma}^{2}\boldsymbol{I}_{d_{\theta}})$,
where $\tilde{\sigma}^{2}=\frac{0.05}{\sqrt{N^{3}}}$ is selected following
[8].

We assess the accuracy of the algorithms in terms of the normalized mean
square error (NMSE) of the estimators of the parameters, the slow state
variables and the fast state variables. In the plots, we show the NMSEs
computed at time $t$,

$\displaystyle\text{NMSE}_{\boldsymbol{\theta},t}$
$\displaystyle=\frac{\|\boldsymbol{\theta}_{t}-\hat{\boldsymbol{\theta}}_{t}\|^{2}}{\|\boldsymbol{\theta}_{t}\|^{2}},$
(62) $\displaystyle\text{NMSE}_{\boldsymbol{x},t}$
$\displaystyle=\frac{\|\boldsymbol{x}_{t}-\hat{\boldsymbol{x}}_{t}\|^{2}}{\|\boldsymbol{x}_{t}\|^{2}},$
(63) $\displaystyle\text{NMSE}_{\boldsymbol{z},t}$
$\displaystyle=\frac{\|\boldsymbol{z}_{ht}-\hat{\boldsymbol{z}}_{ht}\|^{2}}{\|\boldsymbol{z}_{ht}\|^{2}},$
(64)

averaged over 50 independent simulation runs of 20 continuous-time units each,
where the estimators take the form

$\displaystyle\hat{\boldsymbol{\theta}}_{t}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}v_{t}^{i}{\boldsymbol{\theta}}^{i},$ (65)
$\displaystyle\hat{\boldsymbol{x}}_{t}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{J}v_{t}^{i}u_{t}^{i,j}{\boldsymbol{x}}^{i,j}_{t}\quad\text{and}$
(66) $\displaystyle\hat{\boldsymbol{z}}_{ht}$ $\displaystyle=$
$\displaystyle{\sum_{i=1}^{N}\sum_{j=1}^{J}\sum_{l=0}^{L}v_{t}^{i}u_{t}^{i,j}\lambda_{ht}^{i,j,l}{\boldsymbol{z}}^{i,j,l}_{ht},}$
(67)

for the first method. For the second method the estimators of the state
variables are

$\displaystyle\hat{\boldsymbol{x}}_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{J}\sum_{i=1}^{N}\sum_{j=1}^{J}v_{t}^{i}{\boldsymbol{x}}^{i,j}_{t}\quad\text{and}$
(68) $\displaystyle\hat{\boldsymbol{z}}_{ht}$ $\displaystyle=$
$\displaystyle\frac{1}{J}\sum_{i=1}^{N}\sum_{j=1}^{J}v_{t}^{i}{\boldsymbol{z}}^{i,j}_{ht},$
(69)

where $\boldsymbol{x}_{t}^{i,j}$ is the $j$-th member of the ensemble
$\boldsymbol{X}_{t}^{i}$ in Algorithm 5.

Figure 2 shows the performance of the proposed methods for different values of
$J$ (number of samples in the second layer) and $N=20$. This is evaluated in
terms of averaged NMSE θ, NMSE x and NMSE z together with the running time in
hours. The first method (SMC-SMC-UKF) shows an improvement in the accuracy as
the number of samples $J$ increases, although this improvement is only
significant for the slow state (Fig. 2(b)). The second method (SMC-EnKF-EKF)
remains stable with $J$. The second method outperforms the first one in
accuracy of the parameter estimation (Fig. 2(a)) as well as the slow state
estimation (Fig. 2(b)). However, the first method obtains a better NMSE z.
Additionally, the second method runs faster since the computational cost is
considerably lower.

2050100$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$JNMSEθ (a) NMSE θ for different
number of particles/ensembles ($J$) in the second layer of the filter.

2050100$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$JNMSEx (b) NMSE x for different
number of particles/ensembles ($J$) in the second layer of the filter.

2050100$10^{-2}$$10^{-1}$JNMSEz (c) NMSE z for different number of
particles/ensembles ($J$) in the second layer of the filter.

205010002468101214Jrunning time (hours) (d) Running time for different number
of particles/ensembles ($J$) in the second layer of the filter.

Figure 2: Averaged NMSE θ (2(a)), NMSE x (2(b)), NMSE z (2(c)) and average
running time (2(d)) of SMC-SMC-UKF (in blue) and SMC-EnKF-EKF (in red),
averaged over 50 simulation runs. The number of particles of the first layer
(SMC) is set to $N=20$.

Figure 3 compares the performance of the proposed methods and the EnKF for
different values of $N$ (number of samples in the first layer) and $J=50$.
This is shown with the averaged NMSE θ, NMSE x and NMSE z together with the
running time in hours. Similar to the previous figure, the first method (SMC-
SMC-UKF) shows a slight improvement in the accuracy of the slow state
estimation as the number of samples $J$ increases (Fig. 3(b)). The second
method (SMC-EnKF-EKF) remains stable with $N$. The second method outperforms
the first one in accuracy of the parameter estimation (Fig. 3(a)) and the slow
state estimation (Fig. 3(b)), but not for the fast state estimation (Fig.
3(c)). Again, the second method runs faster since the computational cost is
considerably lower.

2050100$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$NNMSEθ (a) NMSE θ for different
number of particles/ensembles ($N$) in the first layer of the filter.

2050100$10^{-4}$$10^{-3}$$10^{-2}$NNMSEx (b) NMSE x for different number of
particles/ensembles ($N$) in the first layer of the filter.

2050100$10^{-2}$$10^{-1}$NNMSEz (c) NMSE z for different number of
particles/ensembles ($N$) in the first layer of the filter.

20501000102030Nrunning time (hours) (d) Running time for different number of
particles/ensembles ($N$) in the first layer of the filter.

Figure 3: Averaged NMSE θ (2(a)), NMSE x (2(b)), NMSE z (2(c)) and average
running time (2(d)) of SMC-SMC-UKF (in blue) and SMC-EnKF-EKF (in red),
averaged over 50 simulation runs. The number of particles/ensembles of the
second layer (SMC for the first method and EnKF for the second method) is set
to $J=50$.

Finally, we show results for a computer experiment in which we have used the
SMC-EnKF-EKF method to estimate the parameters $F$, $C$, $B$ and $H$ and track
the state variables of the two-scale Lorenz 96 system with dimension
$d_{x}=10$ and $d_{z}=50$. The number of particles used to approximate the
sequence of parameter posterior distributions is $N=50$ and the number of
samples in the ensembles of the second layer is $J=50$.

Figure 4 shows the true state trajectories, together with their estimates, for
the first slow state variable ($x_{1}$) and the first fast state variable
($z_{1}$) of the two-scale Lorenz 96 model. We note that although the accuracy
of the estimation of the fast variable is similar throughout the whole
simulation run (over 20 continuous-time units), we only show the last $2$
continuous-time units of the simulation.

$0$$5$$10$$15$$20$$-5$$0$$5$$10$$\tau$${x}_{1}$ (a) Time sequence of $x_{1}$

$18$$18.5$$19$$19.5$$20$$-1$$-0.5$$0$$0.5$$1$$\tau$${z}_{1}$ (b) Time sequence
of $x_{2}$

Figure 4: Sequences of state values (black line) and estimates (dashed red
line) in $x_{1}$ (plot 4(a)) and $z_{1}$ (plot 4(b)) over time.

In Fig. 5 we observe the estimated posterior pdfs of the fixed parameters $F$,
$C$, $B$ and $H$, together with the ground truth values. Figure 5(a) displays
the approximate posterior pdf of the parameter $F$ (red dashed line) together
with the true value $F=8$ (vertical black line), Fig. 5(b) displays the
approximate posterior pdf of the parameter $C$ (blue dashed line) together
with the true value $C=10$ (vertical black line), Fig. 5(c) displays the
approximate posterior pdf of the parameter $B$ (green dashed line) together
with the true value $B=15$ (vertical black line) and Fig. 5(d) displays the
approximate posterior pdf of the parameter $H$ (magenta dashed line) together
with the true value $H=0.75$ (vertical black line). We observe that for all
the pdfs, nearly all probability mass is allocated close to the true values,
except for the parameter $B$ (Fig. 5(c)). In this case, the pdf is slightly
shifted w.r.t. the true value.

$7$$8$$9$$0$$0.5$$1$$F$posterior density (a) Posterior pdf of $F$

91011$0$$0.5$$1$$C$posterior density (b) Posterior pdf of $C$

13141516$0$$0.5$$1$$B$posterior density (c) Posterior pdf of $B$

012$0$$0.5$$1$$H$posterior density (d) Posterior pdf of $H$

Figure 5: Posterior density of the parameters (dashed lines) at time
$\tau=20$. The true values are indicated by a black vertical line.

##  6 Conclusions

We have introduced a further generalization of the NHF methodology of [26]
that, using long sequences of observations collected over time, estimates the
static parameters and the stochastic dynamical variables of a class of
heterogeneous multi-scale state-space models [1]. This scheme combines three
layers of filters, one inside the other. It approximates recursively the
posterior probability distributions of the parameters and the two sets of
state variables given the sequence of available observations. In a first layer
of computation we approximate the posterior probability distribution of the
parameters, in a second layer we approximate the posterior probability
distribution of the slow state variables, and the posterior probability
distribution of the fast state variables is approximated in a third layer. The
inference techniques used in each layer can vary, leading to different
computational costs and degrees of accuracy. To be specific, we describe two
possible algorithms that derive from this scheme, combining Monte Carlo
methods and Gaussian filters at different layers. The first method involves
using sequential Monte Carlo (SMC) methods in both first and second layers,
together with a bank of unscented Kalman filter (UKFs) in the third layer
(i.e., the SMC-SMC-UKF). The second method employs a SMC in the first layer,
ensemble Kalman filters (EnKFs) at the second layer and introduces the use of
a bank of extended Kalman filters (EKFs) in the third layer (i.e., the SMC-
EnKF-EKF). We have presented numerical results for a two-scale stochastic
Lorenz 96 model with synthetic data and we have evaluated the performance of
the algorithm in terms of the normalized mean square errors (NMSEs) for the
parameters and the dynamic (slow and fast) state variables. The proposed
implementations (both of them) obtain good results in terms of accuracy,
having a considerably reduction in running time (i.e., the computational cost)
with the second method. Further research is still needed, studying the
stability of the multi-layer structure when the sequence of observations are
rare and/or few. Moreover, we can compare the proposed algorithms with other
methods, using the Lorenz 96 system but also other models.

## Acknowledgements

This research was partially supported by the Office of Naval Research (award
no. N00014-19-1-2226), Agencia Estatal de Investigación of Spain (ref.
RTI2018-099655-B-I00 CLARA) and Comunidad Autónoma de Madrid (ref.
Y2018/TCS-4705 PRACTICO).

## References

  * [1] Assyr Abdulle, E Weinan, Björn Engquist, and Eric Vanden-Eijnden.  The heterogeneous multiscale method.  Acta Numerica, 21:1–87, 2012.
  * [2] B. D. O. Anderson and J. B. Moore.  Optimal Filtering.  Englewood Cliffs, 1979.
  * [3] C. Andrieu, A. Doucet, and R. Holenstein.  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society B, 72:269–342, 2010.
  * [4] H. M. Arnold.  Stochastic parametrisation and model uncertainty.  PhD thesis, University of Oxford, 2013.
  * [5] Alexandros Beskos, Dan Crisan, Ajay Jasra, Kengo Kamatani, and Yan Zhou.  A stable particle filter for a class of high-dimensional state-space models.  Advances in Applied Probability, 49(1):24–48, 2017.
  * [6] Mallory Carlu, Francesco Ginelli, Valerio Lucarini, and Antonio Politi.  Lyapunov analysis of multiscale dynamics: The slow manifold of the two-scale lorenz’96 model.  arXiv preprint arXiv:1809.05065, 2018.
  * [7] Nicolas Chopin, Pierre E Jacob, and Omiros Papaspiliopoulos.  SMC2: an efficient algorithm for sequential analysis of state space models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(3):397–426, 2013.
  * [8] Dan Crisan, Joaquin Miguez, et al.  Nested particle filters for online parameter estimation in discrete-time state-space Markov models.  Bernoulli, 24(4A):3039–3086, 2018.
  * [9] P. M. Djurić, J. H. Kotecha, J. Zhang, Y. Huang, T. Ghirmai, M. F. Bugallo, and J. Míguez.  Particle filtering.  IEEE Signal Processing Magazine, 20(5):19–38, September 2003.
  * [10] A. Doucet, S. Godsill, and C. Andrieu.  On sequential Monte Carlo Sampling methods for Bayesian filtering.  Statistics and Computing, 10(3):197–208, 2000.
  * [11] Arnaud Doucet and Adam M Johansen.  A tutorial on particle filtering and smoothing: Fifteen years later.  Handbook of nonlinear filtering, 12(656-704):3, 2009.
  * [12] G. Evensen.  The ensemble Kalman filter: Theoretical formulation and practical implementation.  Ocean dynamics, 53(4):343–367, 2003.
  * [13] Ian Grooms and Yoonsang Lee.  A framework for variational data assimilation with superparameterization.  Nonlinear Processes in Geophysics, 22(5):601–611, 2015.
  * [14] Ian Grooms and Gregor Robinson.  A hybrid particle-ensemble Kalman filter for problems with medium nonlinearity.  Plos one, 16(3):e0248266, 2021.
  * [15] S. J. Julier and J. Uhlmann.  Unscented filtering and nonlinear estimation.  Proceedings of the IEEE, 92(2):401–422, March 2004.
  * [16] Nishanth Lingala, N Sri Namachchivaya, Nicolas Perkowski, and Hoong C Yeong.  Particle filtering in high-dimensional chaotic systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4):047509, 2012.
  * [17] J. S. Liu and R. Chen.  Sequential Monte Carlo methods for dynamic systems.  Journal of the American Statistical Association, 93(443):1032–1044, September 1998.
  * [18] Edward N Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on predictability, volume 1, 1996.
  * [19] Jan Mandel.  Efficient implementation of the ensemble Kalman filter.  University of Colorado at Denver and Health Sciences Center, Center for Computational Mathematics, 2006.
  * [20] Henrique MT Menegaz, João Y Ishihara, Geovany A Borges, and Alessandro N Vargas.  A systematization of the unscented Kalman filter theory.  IEEE Transactions on automatic control, 60(10):2583–2598, 2015\.
  * [21] J. Miguez, D. Crisan, and I. P. Mari no.  Particle filtering for Bayesian parameter estimation in a high dimensional state space model.  In Proceedings of the 23rd European Signal Processing Conference (EUSIPCO), September 2015.
  * [22] Christian A Naesseth, Fredrik Lindsten, and Thomas B Schön.  High-dimensional filtering using nested sequential Monte Carlo.  IEEE Transactions on Signal Processing, 67(16):4177–4188, 2019\.
  * [23] Grigoris Pavliotis and Andrew Stuart.  Multiscale methods: averaging and homogenization.  Springer Science & Business Media, 2008.
  * [24] Manuel Pulido, Pierre Tandeo, Marc Bocquet, Alberto Carrassi, and Magdalena Lucini.  Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods.  Tellus A: Dynamic Meteorology and Oceanography, 70(1):1–17, 2018\.
  * [25] J. Uhlmann S. J. Julier and H. F. Durrant-Whyte.  A new method for the non linear transformation of means and covariances in filters and estimators.  IEEE Transactions Automatic Control, 3:477–482, March 2000.
  * [26] J. Míguez S. Pérez-Vieites, I. P. Mariño.  Probabilistic scheme for joint parameter estimation and state prediction in complex dynamical systems.  Physical Review E, 98(6), 063305, 2017.
  * [27] Naratip Santitissadeekorn and Christopher Jones.  Two-stage filtering for joint state-parameter estimation.  Monthly Weather Review, 143(6):2028–2042, 2015.
  * [28] Zheqi Shen, Youmin Tang, Xiaojing Li, Yanqiu Gao, and Junde Li.  On the localization in strongly coupled ensemble data assimilation using a two-scale lorenz model.  Nonlinear Processes in Geophysics Discussions, pages 1–24, 2018\.
  * [29] Tadashi Tsuyuki.  Data assimilation in a two-scale model with Kalman filters.  CAS/JSC WGNE Research Activities in Atmospheric and Oceanic Modelling, pages 27–28, 2012.
  * [30] Eric Vanden-Eijnden.  Heterogeneous multiscale methods: a review.  Communications in Computational Physics2 (3), pages 367–450, 2007\.
  * [31] Eric Vanden-Eijnden et al.  Fast communications: Numerical techniques for multi-scale dynamical systems with stochastic effects.  Communications in Mathematical Sciences, 1(2):385–391, 2003.
  * [32] Gabriele Vissio and Valerio Lucarini.  A proof of concept for scale-adaptive parametrizations: the case of the lorenz’96 model.  Quarterly Journal of the Royal Meteorological Society, 144(710):63–75, 2018.
  * [33] E Weinan.  Principles of multiscale modeling.  Cambridge University Press, 2011.
  * [34] E Weinan and Bjorn Engquist.  Multiscale modeling and computation.  Notices of the AMS, 50(9):1062–1070, 2003.
  * [35] E Weinan, Di Liu, and Eric Vanden-Eijnden.  Analysis of multiscale methods for stochastic differential equations.  Communications on Pure and Applied Mathematics, 58(11):1544–1585, 2005.
  * [36] Hoong C Yeong, Ryne T Beeson, N Sri Namachchivaya, and Nicolas Perkowski.  Particle filters with nudging in multiscale chaotic systems: with application to the lorenz’96 atmospheric model.  Journal of Nonlinear Science, pages 1–34, 2020.