# Adaptive Time Series Forecasting with Markovian Variance Switching

Baptiste Abélès
Universidad Pompeu Fabra
<EMAIL_ADDRESS>&Joseph de Vilmarest
Viking Conseil
<EMAIL_ADDRESS>
&Olivier Wintemberger
Sorbonne Université
<EMAIL_ADDRESS>

###### Abstract

Adaptive time series forecasting is essential for prediction under regime
changes. Several classical methods assume linear Gaussian state space model
(LGSSM) with variances constant in time. However, there are many real-world
processes that cannot be captured by such models. We consider a state-space
model with Markov switching variances. Such dynamical systems are usually
intractable because of their computational complexity increasing exponentially
with time; Variational Bayes (VB) techniques have been applied to this
problem. In this paper, we propose a new way of estimating variances based on
online learning theory; we adapt expert aggregation methods to learn the
variances over time. We apply the proposed method to synthetic data and to the
problem of electricity load forecasting. We show that this method is robust to
misspecification and outperforms traditional expert aggregation.

_K_ eywords Time Series $\cdot$ Online Learning $\cdot$ Kalman Filter $\cdot$
Markov Chains

##  1 Introduction

Time series forecasting is a fundamental issue. A few examples are weather
forecasts for farmers, sales and inventory forecasting in retail, price and
cost forecasting for industrials. In this paper, a special motivation is
electricity load forecasting, a crucial task for grid operators as the
production must balance the consumption in real time [4, 15].

State-space models have been widely employed to model the temporal behavior of
data. In particular, linear Gaussian state-space models (LGSSMs) yield
tractable posterior distributions obtained by the Kalman filter (KF) [19].
Many situations are filled with non linearity both in dynamics and observation
equation. A vast literature has been devoted to deal with non-linear dynamical
systems, e.g. the extended Kalman filter (EKF) [10] which essentially consists
to linearize the model using Taylor first-order expansion, and the unscentend
Kalman filter (UKF) [30] which is radically different and consists to apply
the unscented transformation to both prediction and filtering step.

Although the inference in LGSSMs is known [19], the Kalman filter crucially
depends on the variances of the observation and state noises. The choice of
these variances may be seen as the parametrization of a second-order gradient
descent algorithm, once we remark the equivalence of the Kalman recursions and
gradient updates [26, 8]. The most widely-employed setting of these variances
is time-invariant [11, 10, 6], and the natural choice of these variances is
maximum-likelihood. However, constant variances mean a smooth evolution of the
environment, while in some cases there are changes of regime. To tackle that
problem, various methods have been introduced, often under the paradigm
adaptive Kalman filter [22, 16, 17].

Alternatively the switching Kalman filter (SKF) [24] considers several linear
dynamic models and try to combine them with a switching mechanism (such as a
Markov chain), selecting over time one among several regimes. This latter
approach is of a particular interest in many real-world problems that cannot
be approximated with a single regime, and where multiple behaviors are needed
instead. For instance, time series forecasting during a period covering stable
as well as unstable moments motivates switches between different regimes. Our
application to electricity load forecasting is one of these, as we include
very unstable periods such as lockdowns due to the coronavirus in 2020, during
which it was best to switch to a more reactive method [25].

The complexity of these hybrid dynamical systems, i.e. systems that combine
discrete and continuous valued latent variables, increases exponentially with
time. State-of-the art approaches propose to learn a switching linear
dynamical system using Variational Bayesian (VB) techniques [2]. However, VB
approaches have some limits under model misspecification [31].

In this paper, we propose to consider the SKF as an aggregation of experts,
where each expert corresponds to a different regime. We apply online learning
techniques to track the best expert as in [14, 5]. After introducing various
methods in Section 2, we define our algorithm in Section 3, and we provide a
regret bound in expectation. Then, in Section 4, we conduct various
experiments on synthetic data, as well as a real data set, forecasting the
French national electricity load. We study well-specified and misspecified
experiments and we show that this new framework is robust to misspecification
and outperforms Kalman aggregation as in [1].

##  2 Background

##### Notations

Given probability measures $p$ and $q$ with a common dominating measure $\mu$,
we denote the Kullback-Leibler divergence as $D(p||q):=\int p\log(p/q)d\mu$.
For any positive integer $t$ we denote $[t]:=\\{1,...,t\\}$ and
$\mathcal{F}_{t}:=\sigma(\boldsymbol{x}_{1},y_{1},...,\boldsymbol{x}_{t},y_{t},\boldsymbol{x}_{t+1})$
is the sigma algebra generated by the pair
$(\boldsymbol{x}_{i},y_{i})_{i=1}^{t}$ and $\boldsymbol{x}_{t+1}$. For some
$\mu,\sigma\in\mathbb{R}$, the notation $\mathcal{N}(.|\mu,\sigma)$ refers to
as the normal distribution with mean $\mu$ and variance $\sigma^{2}$ (we keep
the same notation if $\boldsymbol{\mu}$ is a vector and $\Sigma$ is the
covariance matrix associated). We denote $\Delta(\alpha,\beta,\gamma)$ as the
diagonal matrix in $\mathbb{R}^{3\times 3}$ with coefficients
$\alpha,\beta,\gamma$.

Let $y_{t}\in\mathbb{R}$ be the variable of interest,
$\boldsymbol{x}_{t}\in\mathbb{R}^{d}$ a vector of explainable variables and
$\boldsymbol{\theta}_{t}\in\mathbb{R}^{d}$ the vector of latent variables at
time $t$. The classical LGSSM writes :

$\displaystyle
y_{t}=\boldsymbol{\theta}_{t}^{T}\boldsymbol{x}_{t}+\epsilon_{t},\quad\text{where
}\epsilon_{t}\sim\mathcal{N}(0,\sigma_{t}^{2}),$ (1)
$\displaystyle\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}+\nu_{t},\quad\text{where
}\nu_{t}\sim\mathcal{N}(0,Q_{t}),$ (2)

Our goal is to estimate jointly the state and variances knowing the past that
we denote
$\boldsymbol{\hat{\theta}}_{t|t-1}=\mathbb{E}[\boldsymbol{\theta}_{t}|\mathcal{F}_{t-1}]$
and
$P_{t|t-1}=\mathbb{E}[(\boldsymbol{\theta}_{t}-\boldsymbol{\hat{\theta}}_{t|t-1})(\boldsymbol{\theta}_{t}-\boldsymbol{\hat{\theta}}_{t|t-1})^{T}|\mathcal{F}_{t-1}]$.
The well-known KF gives the exact updates if $\sigma_{t}^{2}$ and $Q_{t}$ are
known :

######  Theorem 1.

(Kalman Filter)

$\displaystyle
P_{t|t}=P_{t|t-1}-\frac{P_{t|t-1}\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{T}P_{t|t-1}}{\boldsymbol{x}_{t}^{T}P_{t|t-1}\boldsymbol{x}_{t}+\sigma_{t}^{2}},$
$\displaystyle\boldsymbol{\hat{\theta}}_{t|t}=\boldsymbol{\hat{\theta}}_{t|t-1}-\frac{P_{t|t}}{\sigma_{t}^{2}}\left(\boldsymbol{x}_{t}(\boldsymbol{\hat{\theta}}_{t|t-1}^{T}\boldsymbol{x}_{t}-y_{t})\right),$
$\displaystyle P_{t+1|t}=P_{t|t}+Q_{t},$
$\displaystyle\boldsymbol{\hat{\theta}}_{t+1|t}=\boldsymbol{\hat{\theta}}_{t|t}.$

Usual choice in the literature is to assume that the variances are time-
invariant [10, 6] and apply the Expectation-Maximisation (EM) algorithm. A
second paradigm which has frequently been called Adaptive Kalman Filter (AKF)
[22] consists in estimating variances over time and applying the Kalman
recursive updates using such adaptive variances.

A more recent approach consists in finding the complete posterior of both
state and variances [9]. The authors assume that $Q_{t},\sigma_{t}^{2}$ are
random variables and try to estimate the joint distribution of the state and
the variances in a dynamical way. They assume a model on the dynamics of the
variances; then, starting from a prior on
$(\boldsymbol{\theta}_{1},\sigma_{1}^{2},Q^{(1)})$, they propose to
recursively estimate the posterior distribution of
$(\boldsymbol{\theta}_{t},\sigma_{t}^{2},Q_{t})$. However there is no natural
parametric class of distributions on the joint distribution such that the
posterior remains in the class considered. To handle this issue, they apply
the Variational Bayesian (VB) approach [27]; it consists in approximating the
posterior with the best factorized form in the sense of the minimal Kullback-
Leibler divergence. We believe that the limits of such methods lie on too much
degree of freedom and we aim to constrain the model.

To reduce the degree of freedom of the problem, one can assume the existence
of a hidden Markov chain $(Z_{t})\in[K]$ parameterizing the dynamics of the
classical LGSSM. At each time $t$,
$\sigma_{t}^{2}=\sigma_{Z_{t}}^{2},Q_{t}=Q_{Z_{t}}$. This model writes :

$\displaystyle
y_{t}=\boldsymbol{\theta}_{t}^{T}\boldsymbol{x}_{t}+\epsilon_{t},\quad\text{where
}\epsilon_{t}\sim\mathcal{N}(0,\sigma_{Z_{t}}^{2}),$ (3)
$\displaystyle\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}+\nu_{t},\quad\text{where
}\nu_{t}\sim\mathcal{N}(0,Q_{Z_{t}}),$ (4)
$\displaystyle\mathbb{P}(Z_{t+1}=k)=\sum_{i=1}^{K}M_{ik}\mathbb{P}(Z_{t}=i),$
(5)

where $M_{ik}=\mathbb{P}(Z_{t+1}=k|Z_{t}=i)$ for $i,k\in[K]$ is an entry of
the transition matrix $\boldsymbol{M}\in\mathbb{R}^{K\times K}$.

However, the complexity of the posterior distribution of such dynamical
systems grows exponentially in time, making it intractable. Indeed the
posterior writes :

$\displaystyle p(\boldsymbol{\theta}_{1:t},Z_{1:t}|\mathcal{F}_{t})$
$\displaystyle\propto\pi(Z_{1})\mathcal{N}(\boldsymbol{\theta}_{1}|\boldsymbol{\hat{\theta}}_{1},P_{1})\prod_{s=1}^{t}\mathcal{N}(y_{s}|\boldsymbol{\theta}_{s}^{T}\boldsymbol{x}_{s},\sigma_{Z_{s}}^{2})$
$\displaystyle\prod_{s=1}^{t-1}\mathcal{N}(\boldsymbol{\theta}_{s+1}|\boldsymbol{\theta}_{s},Q_{Z_{s}})M_{Z_{s+1},Z_{s}},$

where $\pi(Z_{1})$ denotes the prior law on $Z_{1}$. Several approaches,
referred to as Generalized Pseudo Bayesian algorithm [24], approximate this
exponential mixture of Gaussians with a smaller mixture of Gaussians. Recent
works use VB techniques to approximate the posterior [2]. The main idea is to
look for a tractable approximation of the posterior. If we denote
$p_{\mathcal{F}_{t}}$ the posterior we look for
$(q_{Z}^{*},q_{\theta}^{*})\in\mathcal{D}_{t}$ such that :

$(q_{Z}^{*},q_{\theta}^{*})\in\arg\min_{(q_{Z},q_{\theta})\in\mathcal{D}_{t}}D(q_{Z}\times
q_{\theta}\ ||\ p_{\mathcal{F}_{t}}),$

where $\mathcal{D}_{t}$ is the parametric class considered. We provide the
derivation and results of this method that performs poorly in our experiments
in supplementary materials. An important limitation of this method is the lack
of robustness to misspecification [31].

To overcome this problem we propose a new approach to adaptive variance
estimation relying on expert aggregation. Aggregation of experts is a well-
established framework in machine learning [5, 29]. In particular, it yields
excellent results in the context of electricity load forecasting [1]. Our work
follows this classical setting, we consider a learner, $K$ experts and a loss
function $\ell:\boldsymbol{\mathcal{E}}\times\mathcal{Y}\to\mathbb{R}$ where
$\boldsymbol{\mathcal{E}}$ is the space of experts predictions and
$\mathcal{Y}$ is the space of observation. At each time $t$, each expert
$k\in[K]$ predicts $\boldsymbol{e}_{k,t}\in\boldsymbol{\mathcal{E}}$. The
learner picks a distribution $\boldsymbol{p}_{t}$ over the set of experts’
predictions, predicts
$\boldsymbol{e}_{t}=\sum_{k=1}^{K}\boldsymbol{p}_{k,t}\boldsymbol{e}_{k,t}\in\boldsymbol{\mathcal{E}}$
and observes his loss $\ell_{t}:=\ell(\boldsymbol{e}_{t},y_{t})$. The learner
has also access to the loss of each expert at time $t$, that we denote
$\ell_{k,t}:=\ell(\boldsymbol{e}_{k,t},y_{t})$. The goal of the learner is to
minimize his cumulative loss on sequence of length $T$. A first way to measure
his performance is to compare his cumulative loss with other expert during the
sequence.

######  Definition 1.

We define the regret with respect to the expert $k$ for a sequence of length
$T$ as :

$R_{T}(k)=\sum_{t=1}^{T}\ell_{t}-\ell_{k,t}$

###  2.1 Exponential-Weight-Average (EWA) algorithm for Kalman Filter

To better understand our framework, we present the famous EWA to show how a
learner can optimize his strategy to control his regret.

We make a small remark here that is we will essentially focus on the variance
of the state Equation 4. Indeed Kalman Filter gives exactly the same update
for the pair $(\sigma^{2},Q)$ and $(1,Q^{*})$ with
$Q^{*}=\frac{Q}{\sigma^{2}}$ (refer to [3]).

Let’s consider the notations of the state space model and assume each expert
is a KF derived with $Q_{t}^{*}=Q^{*(k)}$ which starting from some
$\boldsymbol{\hat{\theta}}_{0|0}$ yields the prediction
$\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)}$ of the state at time $t$, for
$t\in[T]$. In this case, let’s identify the space of experts as the
observation space , i.e $\boldsymbol{\mathcal{E}}=\mathcal{Y}$ and each expert
predicts $\boldsymbol{e}_{k,t}=\hat{y}_{k,t}$ where
$\hat{y}_{k,t}=\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)T}\boldsymbol{x}_{t}$. We
consider the quadratic loss $\ell_{k,t}=(y_{t}-\hat{y}_{k,t})^{2}$.

Algorithm 1 EWA Algorithm

Set $\boldsymbol{p}_{1}=K^{-1}\mathbf{1}_{K}$

for $t\in[T]$ do

select $\boldsymbol{p}_{t}$, forecast
$\hat{y}_{t}=\sum_{k=1}^{K}\boldsymbol{p}_{t}(k)\hat{y}_{k,t}$

derive loss $\ell_{t}$ and $\ell_{k,t}$ for $k\in[K]$

for $k\in[K]$,update:
$\boldsymbol{p}_{t+1}(k)=\frac{\boldsymbol{p}_{t}(k)e^{-\eta\ell_{k,t}}}{\sum_{j=1}^{K}\boldsymbol{p}_{t}(j)e^{-\eta\ell_{j,t}}}$

end for

A first important result is given by the following proposition.

######  Proposition 1.

Assume that the loss function $\ell$ is convex and takes value in $[0,1]$:

$R_{T}^{*}\leq\frac{\log(K)}{\eta}+\frac{\eta T}{8}\,.$

Refer to [5], Theorem 2.2 for a proof.

Algorithm 1 belongs to the so-called class of Static Expert algorithm which
tend to converge to the best expert. However in our case since we assume
switches; we look for aggregation algorithms able to switch from one expert to
the other.

###  2.2 Fixed-Share (FS) algorithm

For this reason we should consider a more refined measure of performance.

######  Definition 2.

Let $k^{T}=(k_{1},...,k_{T})$ be a sequence of experts we can define the
regret with respect to this sequence as:

$R_{T}(k^{T})=\sum_{t=1}^{T}\ell_{t}-\ell_{k_{t},t}$

and

$R_{T}^{*}=\sup_{k^{T}}R_{T}(k^{T})$

A first and naïve approach would be to consider each sequence as a new expert
and try to recover the best expert on the set of every possible sequence of
expert using EWA algorithm starting from a prior distribution $p_{1}(k^{T})$
over this set. However in practice we will face a problem of intractability
since we need to store and update $O(K^{T})$ weights to compare to
$R_{T}(k^{T})$. The problem of controlling the regret with respect to
sequences of experts, known as tracking the best expert, was introduced by
Herbster and Warmuth [14], who proposed the simple FS algorithm with good
regret guarantees. A key fact highlighted in [28] (section 3.3) is that FS,
can be interpreted as the EWA algorithm on sequences of experts under a markov
prior.

Below we present a generalization of FS named as Markov-Hedge (MH) algorithm
which starting from markov chain prior
$\boldsymbol{p}_{1}(k_{1},...,k_{T})=\pi(k_{1})M_{k_{2}|k_{1}}...M_{k_{T}|k_{T-1}}$
collapses to efficient algorithm with bounded regret [23]. For $t\in[T]$,
$M_{k_{t}|k_{t-1}}$ correspond to the entry of the transition matrix
$\boldsymbol{M}\in\mathbb{R}^{K\times K}$ of the markov chain.

Algorithm 2 Markov-Hedge Algorithm

Require:
$\boldsymbol{M}=(M_{k_{t}|k_{t-1}}),\boldsymbol{p}_{1}=K^{-1}\mathbf{1}_{K}$

for $t\in[T]$ do

select $\boldsymbol{p}_{t}$, forecast
$\hat{y}_{t}=\sum_{k=1}^{K}p_{k,t}\hat{y}_{k,t}$

derive loss $\ell_{t}$ and $\ell_{k,t}$ for $k\in[K]$

for $k\in[K]$,update:
$\boldsymbol{p}_{t+1}(k)=\frac{\boldsymbol{p}_{t}(k)e^{-\eta\ell_{k,t}}}{\sum_{j=1}^{K}\boldsymbol{p}_{t}(j)e^{-\eta\ell_{j,t}}}$

update: $\boldsymbol{p}_{t+1}=M\boldsymbol{p}_{t+1}$ i.e for $k\in[K]$:
$\boldsymbol{p}_{t+1}(k)=\sum_{j=1}^{K}M_{k|j}\boldsymbol{p}_{t+1}(j)$

end for

We have the following proposition :

######  Proposition 2.

Algorithm Markov-Hedge, with initial distribution
$\boldsymbol{p}_{1}(k_{1},...,k_{T})=\pi(k_{1})M_{k_{2}|k_{1}}...M_{k_{T}|k_{T-1}}$
and transition matrix $\boldsymbol{M}\in\mathbb{R}^{K\times K}$, guarantees
the following regret bound. For $T\geq 1$ and any sequence
$k^{T}=(k_{1},...,k_{T})$ :

$R_{T}(k^{T})\leq\frac{1}{\eta}\log\frac{1}{\pi(k_{1})}+\frac{1}{\eta}\sum_{t=2}^{T}\log\frac{1}{M_{k_{t}|k_{t-1}}}+\frac{\eta
T}{8}\,.$

Refer to [5] Theorem 5.1 for a proof

##  3 Methodology

Following the remark made in Subsection 2.1 we will first consider
$\sigma_{t}^{2}=1$ for $t\in[T]$.

###  3.1 Well-specified setting

As we saw in Section 2 what is usually done [1] is to perform an aggregation
over the space of observations. This means that we consider our $K$ experts as
$K$ Kalman Filters derived with constant variances
$Q_{t}=Q^{(k)},\sigma_{t}^{2}=1$. The expert $k$ predicts $\hat{y}_{k,t}$ for
all $t\in[T]$ and updates his weights using Algorithm 2. The main contribution
of this paper is to apply the aggregation framework over the space of
variances, as displayed by diagrams in Figure 1.

Figure 1: Aggregation of Kalman (left) versus Kalman with aggregation of
variances (right).

The algorithm that we design is a Kalman filter that competes the best
covariance matrix at each time $t$ using expert aggregation on
$\boldsymbol{\mathcal{Q}}=\\{Q^{(1)},...,Q^{(K)}\\}$. We consider the space of
experts $\boldsymbol{\mathcal{E}}=\boldsymbol{\mathcal{Q}}$ so that at time
$t$, each expert $k$ computes
$P_{t|t-1}^{(k)},\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)}$ using equations of
Theorem 1 and predicts:

$\hat{y}_{k,t}=\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)T}\boldsymbol{x}_{t}.$

The learner observes the loss $\ell_{k,t}=(y_{t}-\hat{y}_{k,t})^{2}$ of each
expert and updates his weights as in Algorithm 2. Then he selects $Q^{(k)}$
with probability $\boldsymbol{p}_{t}(k)$ to compute
$P_{t+1|t},\boldsymbol{\theta}_{t+1|t}$ given his previous predictions
$P_{t|t-1},\boldsymbol{\theta}_{t|t-1}$ and predicts:

$\hat{y}_{t+1}=\boldsymbol{\hat{\theta}}_{t+1|t}^{T}\boldsymbol{x}_{t+1}.$

This yields the randomized version of our algorithm, Kalman Filter Markov-
Hedge (KFMH), detailed in Algorithm 3.

Algorithm 3 Randomized KFMH

Require:
$\boldsymbol{p}_{1}=K^{-1}\mathbf{1}_{K},\boldsymbol{M}=(M_{i,j})_{i,j=1}^{K}$

for $t\in[T]$ do

Reveal $y_{t}$ and ${\boldsymbol{x}}_{t-1}$

for $k\in[K]$ do

compute $P_{t|t-1}^{(k)}$, $\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)}$, and
$\hat{y}_{k,t}=\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)T}\boldsymbol{x}_{t}$

derive loss $\ell_{k,t}$

end for

update:
$\boldsymbol{p}_{t+1}(k)=\frac{\boldsymbol{p}_{t}(k)e^{-\eta\ell_{k,t}}}{\sum_{j=1}^{K}\boldsymbol{p}_{t}(j)e^{-\eta\ell_{j,t}}}$

update: $\boldsymbol{p}_{t+1}=M\boldsymbol{p}_{t+1}$

Select $Q_{t}=Q^{(k)}$ with proba $p_{t+1}(k)$

Compute $P_{t+1|t},\boldsymbol{\hat{\theta}}_{t+1|t}$

Forecast
$\hat{y}_{t+1}=\boldsymbol{\hat{\theta}}_{t+1|t}^{T}\boldsymbol{x}_{t+1}$

end for

The main parameters of this algorithm are $\eta$ and $M$. To reduce the
complexity of the problem we will assume in our experience that $M$ only
depends of a parameter $\alpha$ such that :

$M_{\alpha}=\begin{pmatrix}1-\alpha&\frac{\alpha}{K-1}&\dots&\frac{\alpha}{K-1}\\\
\frac{\alpha}{K-1}&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\frac{\alpha}{K-1}\\\
\frac{\alpha}{K-1}&\dots&\frac{\alpha}{K-1}&1-\alpha\end{pmatrix}$

We can then calibrate $\eta$ and $\alpha$ by grid-search based on the best
cumulative loss.

######  Proposition 3.

Algorithm KFMH randomized, with transition matrix $M_{\alpha}$ and fined
tuning for $\eta$ and $\alpha$, satisfies:

$\mathbb{E}[R_{T}^{*}]\leq{\cal O}(\sqrt{ST\log T})\,,\qquad T\geq 1,$

where $S$ is the number of regime switching of the SKF competitor.

###### Proof.

Without loss of generality we consider that the random selection at step $t$
following the distribution ${\bf p}_{t}$ is independent of ${\cal F}_{t}$ but
included in ${\cal F}_{t+1}$. The tower property yields the identity

$\mathbb{E}\Big{[}\sum_{t=1}^{T}\ell_{t}\Big{]}=\mathbb{E}\Big{[}\sum_{t=1}^{T}\mathbb{E}[\ell_{t}\mid{\cal
F}_{t}]\Big{]}=\mathbb{E}\Big{[}\sum_{t=1}^{T}\sum_{k=1}^{K}p_{t}(k)\ell_{t,k}\Big{]}\,.$

Rewriting the regret bound of Proposition 2 for the linearized loss and for
the specific transition matrix $M_{\alpha}$, we obtain

$\mathbb{E}[R_{T}]\leq\dfrac{S+1}{\eta}\log\dfrac{K}{\alpha}+\dfrac{T-S}{\eta}\log\Big{(}\dfrac{1}{1-\alpha}\Big{)}+\dfrac{\eta
T}{8}\,.$

The desired result follows after an optimisation of the regret bound in
$\alpha$ and $\eta$. ∎

###  3.2 Misspecified setting

In well-specified cases, the SKF competitor is optimal and has small
cumulative quadratic losses, see e.g. [1]. However, in real-world applications
one needs to be more robust to specification. We introduce a new setting where
the prediction is provided by a Kalman filter acting on a state equation
driven by the covariance matrix

$\hat{Q}_{t}=\sum_{k=1}^{K}\boldsymbol{p}_{t}(k)Q^{(k)}\,.$

The interpretation is to transfer the confidence on each experts computed by
$\boldsymbol{p}_{t}$ directly on the state Equation 4. The gain is to design
an estimation of the state variance and a forecaster that is a KF and not an
aggregation. One also has to design each expert $k$ in a different way than
pure independent KF. Indeed, while in Algorithm 3 each expert was computing
his prediction independently from the learner; now the expert starts from a
prior knowledge $P_{t-2|t-2},\boldsymbol{\hat{\theta}}_{t-1|t-2}$ of learner’s
previous predictions and derives new KF updates as follows:

$\displaystyle P_{t-1|t-2}^{(k)}=P_{t-2|t-2}+Q^{(k)},$ $\displaystyle
P_{t-1|t-1}^{(k)}=P_{t-1|t-2}^{(k)}-\frac{P_{t-1|t-2}^{(k)}\boldsymbol{x}_{t-1}\boldsymbol{x}_{t-1}^{T}P_{t-1|t-2}^{(k)}}{\boldsymbol{x}_{t-1}^{T}P_{t-1|t-2}^{(k)}\boldsymbol{x}_{t-1}+1},$
$\displaystyle\boldsymbol{\hat{\theta}}_{t-1|t-1}^{(k)}=\boldsymbol{\hat{\theta}}_{t-1|t-2}-P_{t-1|t-1}^{(k)}\boldsymbol{x}_{t}(\boldsymbol{\hat{\theta}}_{t-1|t-2}^{T}\boldsymbol{x}_{t}-y_{t}),$
$\displaystyle\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)}=\boldsymbol{\hat{\theta}}_{t-1|t-1}^{(k)},$

to predict:

$\hat{y}_{k,t}=\boldsymbol{\hat{\theta}}_{t|t-1}^{T(k)}\boldsymbol{x}_{t}.$

Then the learner observes the loss $\ell_{k,t}$ of each expert, updates the
weights and computes $\boldsymbol{\hat{\theta}}_{t+1|t},P_{t+1|t}$ applying KF
updates of Theorem 1 with $\hat{Q}_{t+1}$.

###  3.3 Choice of loss and optimization over $\sigma$

In Section 2 we considered the square loss. However, standard methods
optimizing time-invariant variances in state-space models optimize the
likelihood as the natural objective function [10, 6]. This leads us to
consider the log-likelihood as loss function in our algorithm. Using the chain
rule, the log-likelihood writes:

$\displaystyle\log p(y_{1:T},\boldsymbol{x}_{1:T}|Q_{1:T-1},\sigma_{1:T}^{2})$
$\displaystyle\qquad=\sum_{t=1}^{T}\log
p(y_{t}|\boldsymbol{x}_{t},y_{1:t-1},\boldsymbol{x}_{1:t-1},Q_{1:(t-1)},\sigma_{1:t}^{2})+c$
$\displaystyle\qquad=\sum_{t=1}^{T}\ell_{t}(\boldsymbol{\hat{\theta}}_{t|t-1},P_{t|t-1},\sigma_{t}^{2},y_{t})+c\,.$

where $c$ is a constant term and the loss function writes:

$\displaystyle\ell_{t}(\boldsymbol{\hat{\theta}}_{t|t-1},P_{t|t-1},\sigma_{t}^{2},y_{t})$
$\displaystyle\quad=\log\mathcal{N}(y_{t}|\boldsymbol{\hat{\theta}}_{t|t-1}^{T}\boldsymbol{x}_{t},\sigma_{t}^{2}+\boldsymbol{x}_{t}^{T}P_{t|t-1}\boldsymbol{x}_{t})\,.$

The dependence of this loss to $Q_{t-1}$ is hidden in the estimated state
distribution represented by $\boldsymbol{\hat{\theta}}_{t\mid t-1},P_{t\mid
t-1}$. The peculiarity of the log-likelihood is that it strongly depends of
$\sigma_{t}^{2}$. This means that the assumption $\sigma_{t}^{2}=1$ is too
restrictive. In our experiments we observed that it did not perform well to
apply the expert framework to both variances. We believe that it’s due to the
fact that $Q_{t-1}$ and $\sigma_{t}^{2}$ are strongly related. Furthermore,
estimating observation noise variance is an easier task compared to the
estimation of the covariance matrix of a latent process. Therefore we directly
apply a gradient descent algorithm, namely ADAM [20]. Equation $(*)$ in
Algorithm 4 and function $f(.)$ refers to this ADAM gradient descent which
involves 4 hyper-parameters (namely $\alpha_{1},\beta_{1},\beta_{2},\epsilon)$
and the gradient of the loss over $\sigma$ which writes :

$\frac{\partial\ell_{t}}{\partial\sigma}=-\frac{\sigma(y-\boldsymbol{\hat{\theta}}_{t|t-1}^{T}\boldsymbol{x}_{t})^{2}}{\sigma^{2}+\boldsymbol{x}_{t}^{T}P_{t|t-1}\boldsymbol{x}_{t}}+\frac{\sigma}{\sigma^{2}+\boldsymbol{x}_{t}^{T}P_{t|t-1}\boldsymbol{x}_{t}}.$

Note that by decoupling the estimation of both variances, we don’t need a
priori on the possible values of $\sigma_{t}^{2}$ and their relation to the
possible values of $Q_{t}$.

###  3.4 Sliding Window

We propose another improvement of Algorithm 4 motivated by [17]. The main
intuition behind this is to better discriminate between the different experts.
We note that $P_{t-1\mid t-1}$ depends on $Q_{1},Q_{2},\ldots Q_{t-2}$. Our
initial algorithm proposes that expert $k$ depends on $Q_{1},Q_{2},\ldots
Q_{t-3},Q^{(k)}$, that is, we replace only $Q_{t-2}$. Instead, we propose to
replace the last $\tau$ covariance matrices by $Q^{(k)}$ for each expert $k$.
This replacement of a window $\tau$ is also motivated by the fact we assume
$\alpha\ll 1$ in our experiments (typically $\alpha=0.01$), that is, the
covariance matrix is not switching too frequently from one mode to another.
Considering this sliding window helps capturing the dynamics by discriminating
better between the experts.

Formally, instead of starting from
$\boldsymbol{\hat{\theta}}_{t-2|t-2},P_{t-2|t-2}$, we start from
$\boldsymbol{\hat{\theta}}_{t-\tau-1|t-\tau-1},P_{t-\tau-1|t-\tau-1}$. For
each expert $k$, we derive Kalman updates using $Q_{s}=Q^{(k)}$ for $s\geq
t-\tau-1$, yielding $\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)},P_{t\mid
t-1}^{(k)}$. Consequently, this sliding window is integrated in Algorithm 4 as
the introduction of a new loss function
$\ell^{(\tau)}_{k,t}=\mathcal{N}(y_{t}|\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)T}\boldsymbol{x}_{t},\sigma_{t}^{2}+\boldsymbol{x}_{t}^{T}P_{t|t-1}^{(k)}\boldsymbol{x}_{t})$
where the dependency in $\tau$ refers to the way we create
$\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)},P_{t|t-1}^{(k)}.$

Algorithm 4 KFMH

Require:
$\boldsymbol{p}_{1}=K^{-1}\mathbf{1}_{K},\boldsymbol{M}=(M_{i,j})_{i,j=1}^{K},\sigma_{0}$

for $t\in[T]$ do

Reveal $y_{t}$ and ${\boldsymbol{x}}_{t-1}$.

for $k\in[K]$ do

compute $P_{t|t-1}^{(k)}$, $\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)}$, and
$\hat{y}_{k,t}=\boldsymbol{\hat{\theta}}_{t|t-1}^{(k)T}\boldsymbol{x}_{t}$

derive loss $\ell_{k,t}^{(\tau)}$

end for

for $k\in[K]$ update:
$\boldsymbol{p}_{t+1}(k)=\frac{\boldsymbol{p}_{t}(k)e^{-\eta\ell_{k,t}}}{\sum_{j=1}^{K}\boldsymbol{p}_{t}(j)e^{-\eta\ell_{j,t}}}$

update: $\boldsymbol{p}_{t+1}=M\boldsymbol{p}_{t+1}$

Select $Q_{t}=\sum_{k=1}^{K}p_{t}(k)Q^{(k)}$

Compute $P_{t+1|t},\boldsymbol{\hat{\theta}}_{t+1|t}$

Forecast
$\hat{y}_{t+1}=\boldsymbol{\hat{\theta}}_{t+1|t}^{T}\boldsymbol{x}_{t+1}$

Observe loss
$\ell_{t}(\boldsymbol{\hat{\theta}}_{t|t-1},P_{t|t-1},\sigma_{t}^{2},y_{t})$

Update
$\sigma_{t}=\sigma_{t-1}+f\left(\frac{\partial\ell_{t}(\boldsymbol{\hat{\theta}}_{t|t-1},P_{t|t-1},y_{t})}{\partial\sigma}\right)$
(*)

end for

##  4 Experiments

We compare the different methods introduced on various data sets. We begin
with synthetic data; first we generate the data following the state-space
model introduced in Section 2, we call that well-specified and it is the most
favourable setting; then, we generate data not following the state-space
assumption anymore, we call that experiment misspecified and it is an
intermediate step before the application to real data. We conclude by a real
data set motivating our setting: electricity load forecasting. We show that
our method outperforms the standard Kalman aggregation as long as the experts
are well-known even in the misspecified context.

###  4.1 Experiments on Synthetic Data

In the following we consider $K=2$.

####  4.1.1 Well-Specified (WS) Data with $\sigma_{t}^{2}=\sigma^{2}$

##### Data Generation

We consider the model given by Equations $(\ref{SKM::var::eq_obs})$,(4) and
(5) and generate data as follow :

$\displaystyle\mathbb{P}(Z_{t+1}=Z_{t})$ $\displaystyle=1-\alpha,$
$\displaystyle\boldsymbol{\theta}_{t+1}-\boldsymbol{\theta}_{t}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0}_{\mathbb{R}^{3}},Q_{Z_{t}}),$
$\displaystyle y_{t}-\boldsymbol{\theta}_{t}^{T}\boldsymbol{x}_{t}$
$\displaystyle\sim\mathcal{N}(0,1),$

with $Q_{Z_{t}}\in\\{Q^{(1)},Q^{(2)}\\}$ where:

$\displaystyle Q^{(1)}$ $\displaystyle=\Delta(10^{-4},10^{-1},10^{-2}),$
$\displaystyle Q^{(2)}$ $\displaystyle=\Delta(10^{-1},10^{-4},10^{-2}).$

We distinguish three ways of generating $\boldsymbol{x}_{t}$: i.i.d Gaussian
design:
$\boldsymbol{x}_{t}\sim\mathcal{N}(\boldsymbol{0}_{\mathbb{R}^{d}},\boldsymbol{1}_{\mathbb{R}^{d}}),$
i.i.d uniform design:
$\boldsymbol{x}_{t}\sim\mathcal{U}(\boldsymbol{0}_{\mathbb{R}^{d}},\boldsymbol{1}_{\mathbb{R}^{d}})$
and non-i.i.d design: $\boldsymbol{x}_{1}$ is generated as before. Then for
$j\in[\\![1,3]\\!]$ and $t\geq 2$, we consider
$\boldsymbol{a}_{t,j}=\boldsymbol{x}_{t-1,j}+\epsilon_{t,j}$ where
$\epsilon_{t,j}\sim\mathcal{N}(0_{\mathbb{R}^{3}},10^{-3})$ and we use
$\boldsymbol{x}_{t,j}=\boldsymbol{a}_{t,j}$ if $0\leq\boldsymbol{a}_{t,j}\leq
1$, otherwise
$\boldsymbol{x}_{t,j}=\lceil\boldsymbol{a}_{t,j}\rceil-\boldsymbol{a}_{t,j}$.
We generate sequence of length $T=1000$ with a transition parameters
$\alpha=10^{-2}$.

Algorithm Parameterization
For both Algorithm 2 and 4 we set $\eta=1$, and take latence $\tau=5$ for KFMH
and initial $\sigma_{0}=0.8$. For ADAM optimizer $f(.)$ we take the step-size
$\alpha_{1}=0.01$ and the others standard hyper-parameters as in [20]
$(\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=10^{-8}).$

####  4.1.2 Misspecified (MS) Data

In this subsection we forget the Markov representation and consider that we
play against a sequence of expert which is generated by some specific
function.

##### Data generation : Sinusoidal

As discussed before, we assume here that $\sigma_{t}^{2}$ and $Q_{t}$ are
independent. We generate $Q_{t}$ and $\sigma_{t}^{2}$ in the following way :

$Q_{t}=\cos^{2}\left(\frac{3\pi t}{T}\right)Q^{(1)}+\sin^{2}\left(\frac{3\pi
t}{T}\right)Q^{(2)},$ (6) $\sigma_{t}=\sigma^{(1)}\cos^{2}\left(\frac{5\pi
t}{T}\right)+\sigma^{(2)}\sin^{2}\left(\frac{5\pi t}{T}\right),$ (7)

where $Q^{(1)},Q^{(2)}$ are as in Subsection 4.1.1 and
$\sigma^{(1)}=1,\sigma^{(2)}=2$.

###  4.2 Real Data

Lastly we present our results on Real Data. We apply our methodology to the
case study of [25]. We forecast the half-hourly French national electricity
consumption published by Réseau de Transport d’Électricité
111https://opendata.reseaux-energies.fr/, as well as meteorological data from
MétéoFrance222https://donneespubliques.meteofrance.fr/. We consider the
adaptation of the following generalized additive model, built independently
for the 48 half-hours of the day:

$\displaystyle
y_{t}=f(Toy_{t})+\sum\limits_{i=1}^{7}\alpha_{i}\mathbf{1}_{DayType_{t}=i}$
(8) $\displaystyle\quad+f(TempS95_{t})+f(TempS99_{t})+f(t,Temp_{t})$
$\displaystyle\quad+s(CPWind_{t})+f(t)+\beta Lag1_{t}+\gamma
Lag7_{t}+\varepsilon_{t}$

where at each time $t$, $y_{t}$ is the consumption, $DayType_{t}$ is the day
of the week between $1$ and $7$, $Temp_{t}$ is the temperature, $TempS95_{t}$
and $TempS99_{t}$ are exponential smoothings of the temperature of exponential
parameters $0.95$ and $0.99$ respectively, $CPWind_{t}$ is the cooling power
of wind (c.f. [21], Section 2.2), $Lag1_{t}$ and $Lag7_{t}$ are the delayed
consumption (one day ago and seven days ago, respectively).

We follow the framework of [25]: we consider the adaptation of a linear
combination of the frozen effects. Formally, Equation (8) yields a GAM of 9
effects, and we define $\boldsymbol{x}_{t}\in\mathbb{R}^{9}$ the vector
composed of these linear and non-linear effects. Our state-space paradigm is
integrated naturally as

$\displaystyle
y_{t}=\boldsymbol{\theta}_{t}^{T}\boldsymbol{x}_{t}+\epsilon_{t},\qquad\epsilon_{t}\sim\mathcal{N}(0,\sigma_{t}^{2})$
$\displaystyle\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}+\nu_{t},\qquad\nu_{t}\sim\mathcal{N}(0,Q_{t})$

Our data ranges from January 2012 to December 2022. It has been extensively
studied that the coronavirus crisis has implied a brutal drop in the
electricity consumption in many countries, and France in particular [18, 12].
State-space models have provided a very efficient framework to tackle the poor
performances of traditional forecasting methods during this period [25, 7]. We
propose to apply a switching strategy, in order to win in stable as well as
unstable periods. Intuitively, defining a slow covariance matrix $Q^{(slow)}$
and a fast covariance matrix $Q^{(fast)}$, a switching Kalman-based algorithm
is able to use $Q^{(slow)}$ in slow mode and $Q^{(fast)}$ in fast mode,
yielding improvements in both worlds.

To define $Q^{(slow)}$ and $Q^{(fast)}$, we rely on an optimization of
constant variances introduced in [25], Section 2.A.2. However, we optimize the
variances on different time periods. Optimizing on pre-covid data (2012-2019)
yields the slow mode, and optimizing on covid times (the year 2020) yields the
fast mode. We apply the same optimization method to define $\sigma^{2(slow)}$
and $\sigma^{2(fast)}$ associated. We compute our test performances on the
period ranging from January 2021 to December 2022.

###  4.3 Results

Figure 2: Evolution of the weights $\boldsymbol{p}_{t}$ for WS (left) and MS
Data (right).

Figure 3: Evolution of $\sigma_{t}$ starting from 2 different initial points
for WS (left) and MS Data (right).

The results of our experiments on synthetic data are gathered in table 1 (see
Supplementary materials). We generate $n_{iter}=1000$ sequences of length $T$
and compare our algorithm KFMH to differents methods using the Mean Square
Error (MSE). First we consider the KF-adaptive which is the optimal Kalman
filter knowing $Q_{t},\sigma_{t}^{2}$; this is the oracle or best possible
forecast. Then we observe the performance of experts $1$ and $2$ which are
Kalman filters with constant covariance matrix $Q_{t}=Q^{(1)}$ or
$Q_{t}=Q^{(2)}$; we denote them by KF-$Q^{(1)}$ and KF-$Q^{(2)}$. We consider
another Kalman filter with constant covariance matrix, KF-$Q^{(mean)}$, which
applies $Q_{t}=\frac{Q^{(1)}+Q^{(2)}}{2}$. This is a better containder for our
algorithm since in both WS and MS scenarios we generate data such that we
approximately spend the same time in both regimes. Finally we present the
results of KF-$agg$ the method that executes Algorithm 2. Except for our
algorithm, every method has access to the true $\sigma_{t}^{2}$. We see that
our algorithm outperforms the basic aggregation over the space of observation
in every case. We observe that in the MS setting, KF-$Q^{(mean)}$ with
knowledge of $\sigma_{t}^{2}$ beats our algorithm when generating i.i.d
uniform and non-i.i.d $\boldsymbol{x}_{t}$, although we are very close.
Overall we observe that $Q^{(mean)}$ yields good performances on synthetic
data; this was to expect because the covariance matrix of the latent process
is $Q^{(1)}$ half of the time and $Q^{(2)}$ the other half. The results on
real data are presented in table 2 (see Supplementary materials). For our
evaluation, we consider the first ten instants of the day over the $2021-2022$
period. Since we did not generate the data, KF-adaptive exists no longer and
we only compare KFMH to the Kalman filters of constant variances
(KF-$Q^{(slow)}$, KF-$Q^{(fast}$, KF-$Q^{(mean)}$), as well as the standard
expert aggregation. This time we don’t have access to the knowledge of
$\sigma_{t}^{2}$ so each constant expert assumes $\sigma_{t}^{2}=\sigma^{2}$
where $\sigma^{2}$ is respectively $\sigma^{2(slow)}$, $\sigma^{2(fast)}$ and
$\sigma^{2(mean)}=\frac{\sigma^{2(slow)}+\sigma^{2(fast)}}{2}$. We find that
we always outperform every constant expert, and beat the standard aggregation
in average.

Figure LABEL:fig:expe shows the evolution of the weights $\boldsymbol{p}_{t}$
obtained by Algorithm 4 (in blue) and Algorithm 2 (in green) in well-specified
as well as misspecified settings. We also display the evolution of
$\sigma_{t}$ estimated with Algorithm 4 starting from two different initial
points; respectively $\sigma_{0}=0.5$ (blue) and $\sigma_{0}=2$ (green) with
respect to the true $\sigma_{t}$ (black) both in WS case ($\sigma_{t}^{2}=1)$
and MS context where $\sigma_{t}$ is generated with Equation 7. The figures
suggest the consistency of the method and its robustness towards
misspecification.

Below we provide the results of the experiments described in Section 4. We see
that our algorithm is always better than a standard aggregation when
performing on synthetic datasets, and always very close to the optimal.

Table 1: Comparison of MSE for each method on synthetic datas in both WS and MS case. Method | KF-adpative | KF-$Q^{(1)}$ | KF-$Q^{(2)}$ | KF-$Q^{(mean)}$ | KF-$agg$ | KFMH
---|---|---|---|---|---|---
i.i.d gauss WS | $1.802$ | $5.091$ | $5.523$ | $1.990$ | $2.519$ | 1.916
i.i.d unif WS | $1.796$ | $4.021$ | $4.393$ | $1.931$ | $2.301$ | 1.891
non-i.i.d WS | $1.251$ | $1.471$ | $1.415$ | $1.267$ | $1.301$ | 1.266
i.i.d gauss MS | $1.908$ | $5.192$ | $5.526$ | $1.988$ | $2.799$ | 1.986
i.i.d unif MS | $1.370$ | $2.367$ | $2.343$ | 1.393 | $1.653$ | $1.402$
non-i.i.d WS | $2.758$ | $3.040$ | $2.964$ | 2.768 | $2.823$ | $2.777$
Table 2: Comparison of RMSE over 2021-2022 period for each method on real data for 10 instants of the day. Method | KF-$Q^{(slow)}$ | KF-$Q^{(fast)}$ | KF-$Q^{(mean)}$ | KF-$agg$ | KFMH
---|---|---|---|---|---
Instant 1 | $1005.572$ | $805.620$ | $796.485$ | $848.549$ | 788.284
Instant 2 | $888.951$ | $808.081$ | $797.431$ | $790.475$ | 786.788
Instant 3 | $896.954$ | $838.322$ | $827.544$ | 810.124 | $812.846$
Instant 4 | $820.090$ | $876.641$ | $827.933$ | 812.816 | $813.703$
Instant 5 | $842.068$ | $848.340$ | $844.974$ | $849.140$ | 838.871
Instant 6 | $842.665$ | $846.738$ | $845.461$ | $842.816$ | 840.530
Instant 7 | $854.612$ | $908.600$ | $858.957$ | $895.452$ | 849.253
Instant 8 | $853.396$ | $926.782$ | $860.292$ | 843.292 | $847.000$
Instant 9 | $863.809$ | $935.049$ | $871.084$ | 855.107 | $859.937$
Instant 10 | $958.609$ | $922.263$ | $912.158$ | $937.550$ | 901.185

##  5 Conclusion

In this paper, we have studied LGSSMs under unknown variances. We relied on
online learning to estimate the variances. We considered the paradigm of SKF
and used aggregation of experts to infer the variance of the process noise; to
approximate the observation noise variance we applied an online gradient
descent algorithm. Our experiments show good performances not only in the
favourable well-specified setting but also under misspecification on both
synthetic and real data, illustrating the robustness of our method.

This paper opens different leads for future research. First, we specified on
KF algorithms but we could apply our framework to other forecasting
algorithms, involving the expert aggregation confidence weights at the right
level, not necessarily on the predictions. Second, the choice of the diversity
among the experts is crucial in aggregation [13], and it is challenging to
design different KF algorithms in an optimal way.

## References

  * [1] Eric Adjakossa, Yannig Goude, and Olivier Wintenberger.  Kalman recursions aggregated online.  Statistical Papers, pages 1–36, 2023.
  * [2] Xavier Alameda-Pineda, Vincent Drouard, and Radu Patrice Horaud.  Variational inference and learning of piecewise linear dynamical systems.  IEEE Transactions on Neural Networks and Learning Systems, 33(8):3753–3764, 2021.
  * [3] Sileye Ba, Xavier Alameda-Pineda, Alessio Xompero, and Radu Horaud.  An on-line variational bayesian model for multi-person tracking from cluttered scenes.  Computer Vision and Image Understanding, 153:64–76, 2016.
  * [4] D Bunn and E Dillon Farmer.  Comparative models for electrical load forecasting.  1985\.
  * [5] Nicolo Cesa-Bianchi and Gábor Lugosi.  Prediction, learning, and games.  Cambridge university press, 2006.
  * [6] Peter J Brockwell Richard A Davis.  Introduction to time series and forecasting.  springer publication, 2016.
  * [7] Joseph De Vilmarest and Yannig Goude.  State-space models for online post-covid electricity load forecasting competition.  IEEE Open Access Journal of Power and Energy, 9:192–201, 2022.
  * [8] Joseph De Vilmarest and Olivier Wintenberger.  Stochastic online optimization using kalman recursion.  The Journal of Machine Learning Research, 22(1):10173–10227, 2021.
  * [9] Joseph de Vilmarest and Olivier Wintenberger.  Viking: Variational bayesian variance tracking.  arXiv preprint arXiv:2104.10777, 2021.
  * [10] James Durbin and Siem Jan Koopman.  Time series analysis by state space methods, volume 38.  OUP Oxford, 2012.
  * [11] Ludwig Fahrmeir, Gerhard Tutz, Wolfgang Hennevogl, and Eliane Salem.  Multivariate statistical modelling based on generalized linear models, volume 425.  Springer, 1994.
  * [12] Mostafa Farrokhabadi, Jethro Browell, Yi Wang, Stephen Makonin, Wencong Su, and Hamidreza Zareipour.  Day-ahead electricity demand forecasting competition: Post-covid paradigm.  IEEE Open Access Journal of Power and Energy, 9:185–191, 2022.
  * [13] Pierre Gaillard and Yannig Goude.  Forecasting electricity consumption by aggregating experts; how to design a good set of experts.  In Modeling and stochastic learning for forecasting in high dimensions, pages 95–115. Springer, 2015.
  * [14] Mark Herbster and Manfred K Warmuth.  Tracking the best expert.  Machine learning, 32:151–178, 1998.
  * [15] Tao Hong, Pierre Pinson, and Shu Fan.  Global energy forecasting competition 2012, 2014.
  * [16] Yulong Huang, Yonggang Zhang, Zhemin Wu, Ning Li, and Jonathon Chambers.  A novel adaptive kalman filter with inaccurate process and measurement noise covariance matrices.  IEEE transactions on Automatic Control, 63(2):594–601, 2017.
  * [17] Yulong Huang, Fengchi Zhu, Guangle Jia, and Yonggang Zhang.  A slide window variational adaptive kalman filter.  IEEE Transactions on Circuits and Systems II: Express Briefs, 67(12):3552–3556, 2020.
  * [18] IEA.  Year-on-year change in weekly electricity demand, weather corrected, in selected countries.  https://www.iea.org/data-and-statistics/charts/year-on-year-change-in-weekly-electricity-demand-weather-corrected-in-selected-countries-january-december-2020, 2020.
  * [19] Rudolph Emil Kalman.  A new approach to linear filtering and prediction problems.  1960\.
  * [20] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [21] Nicole Ludwig, Siddharth Arora, and James W Taylor.  Probabilistic load forecasting using post-processed weather ensemble predictions.  Journal of the Operational Research Society, 74(3):1008–1020, 2023.
  * [22] Raman Mehra.  On the identification of variances and adaptive kalman filtering.  IEEE Transactions on automatic control, 15(2):175–184, 1970.
  * [23] Jaouad Mourtada and Odalric-Ambrym Maillard.  Efficient tracking of a growing number of experts.  In International Conference on Algorithmic Learning Theory, pages 517–539. PMLR, 2017.
  * [24] Kevin P Murphy.  Switching kalman filters.  1998\.
  * [25] David Obst, Joseph De Vilmarest, and Yannig Goude.  Adaptive methods for short-term electricity load forecasting during covid-19 lockdown in france.  IEEE transactions on power systems, 36(5):4754–4763, 2021.
  * [26] Yann Ollivier.  Online natural gradient as a kalman filter.  2018\.
  * [27] Václav Šmídl and Anthony Quinn.  The variational Bayes method in signal processing.  Springer Science & Business Media, 2006.
  * [28] Vladimir Vovk.  Derandomizing stochastic prediction strategies.  In Proceedings of the tenth annual conference on Computational learning theory, pages 32–44, 1997.
  * [29] Vladimir G Vovk.  A game of prediction with expert advice.  In Proceedings of the eighth annual conference on Computational learning theory, pages 51–60, 1995.
  * [30] Eric A Wan and Rudolph Van Der Merwe.  The unscented kalman filter for nonlinear estimation.  In Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium (Cat. No. 00EX373), pages 153–158. Ieee, 2000.
  * [31] Yixin Wang and David Blei.  Variational bayes under model misspecification.  Advances in Neural Information Processing Systems, 32, 2019.

##  Appendix A VB Derivation

In this section we provide the details of the derivation of the VB approach,
mentioned in Section 2, that we found ineffective.

#####

As mentioned VB approach looks for the pair
$(q_{\theta}^{*},q_{Z}^{*})=(q_{\theta}^{*}(\boldsymbol{\theta}_{1:t}),q_{Z}^{*}(Z_{1:t}))$
that is the closest to the complete (i.e taken over the whole sequence) joint
posterior
$p_{\mathcal{F}_{t}}=p(\boldsymbol{\theta}_{1:t},Z_{1:t}|\mathcal{F}_{t})$ :

$(q_{\theta}^{*},q_{Z}^{*})\in\text{argmin}_{(q_{\theta},q_{Z})\in\mathcal{D}}\
KL(q_{\theta}\times q_{Z}\ ||\ p_{\mathcal{F}_{t}})$

One can show that :

$\ln
q_{\theta}^{*}(\boldsymbol{\theta}_{1:t})=\mathbb{E}_{q_{Z}(Z_{1:t})}\left[\ln
p_{\mathcal{F}_{t}}\right]+c_{1}\qquad\qquad\qquad\text{(E-$\theta$ step)}$
(9) $\ln
q_{Z}^{*}(Z_{1:t})=\mathbb{E}_{q_{\theta}(\boldsymbol{\theta}_{1:t})}\left[\ln
p_{\mathcal{F}_{t}}\right]+c_{2}\qquad\qquad\qquad\text{(E-$Z$ step)}$ (10)

The aim is then to iteratively approximate the marginal $q_{\theta}^{*}$ and
$q_{Z}^{*}$, assuming the other is known and fixed. We provide further the
closed form of what we call respectively E-$\theta$ step and E-$Z$ step.

###  A.1 E-$\theta$ step

######  Proposition 4.

(E-$\theta$ step)
For some $t>1$, for any $q_{Z}(Z_{1:t})$ if we define for $s\in[\\![1,t]\\!]$
:

$\displaystyle\overline{\sigma}_{s}^{-1}$
$\displaystyle=\sqrt{\overset{K}{\underset{k=1}{\sum}}q_{Z}(Z_{s}=k)\frac{1}{\sigma_{k}^{2}}}$
$\displaystyle\overline{Q}_{s}^{-1}$
$\displaystyle=\overset{K}{\underset{k=1}{\sum}}q_{Z}(Z_{s}=k)Q_{k}^{-1}$

then the distribution $q_{\theta}^{*}$ which minimizes $KL(q_{\theta}\times
q_{Z}||\ p_{\mathcal{F}_{t}})$ for $q_{Z}$ fixed writes:

$\ln
q_{\theta}^{*}(\boldsymbol{\theta}_{1:t})=\sum_{s=1}^{t}\Bigg{(}\ln\mathcal{N}(y_{s}|\boldsymbol{\theta}_{s}^{T}\boldsymbol{x}_{s},\overline{\sigma}_{s}^{2})+\ln\mathcal{N}(\boldsymbol{\theta}_{s}|\boldsymbol{\theta}_{s-1},\overline{Q}_{s-1})\Bigg{)}+c$
(11)

where we denote $\boldsymbol{\theta}_{0}=\boldsymbol{\hat{\theta}}_{1}$ and
$\overline{Q}_{0}=P_{1}$. .

###### Proof.

For $s\in[1:t]$ we have:

$\displaystyle\mathbb{E}_{q_{Z}}\left[\ln\mathcal{N}(y_{s}|\theta^{T}_{s}x_{s},\sigma^{2}_{Z_{s}})\right]$
$\displaystyle=-\frac{1}{2}\sum_{1\leq k_{1},\ldots,k_{t}\leq
K}q_{Z}(Z_{1:t}=k_{1},\ldots,k_{t})\left(\frac{\boldsymbol{\theta}_{s}^{T}\boldsymbol{x}_{s}\boldsymbol{x}_{s}^{T}\boldsymbol{\theta}_{s}}{\sigma_{k_{s}}^{2}}-2y_{s}\frac{\theta^{T}_{s}\boldsymbol{x}_{s}}{\sigma_{k_{s}}^{2}}\right)+c$
$\displaystyle=-\frac{1}{2}\sum_{k_{s}=1}^{K}q_{Z}(Z_{s}=k_{s})\left(\frac{\boldsymbol{\theta}^{T}_{s}\boldsymbol{x}_{s}x_{s}^{T}\boldsymbol{\theta}_{s}}{\sigma_{k_{s}}^{2}}-\frac{2y_{s}\boldsymbol{\theta}^{T}_{s}\boldsymbol{x}_{s}}{\sigma_{k_{s}}^{2}}\right)+c$
$\displaystyle\mathbb{E}_{q_{Z}}\left[\ln\mathcal{N}(\boldsymbol{\theta}_{s}|\boldsymbol{\theta}_{s-1},Q_{Z_{s}})\right]$
$\displaystyle=-\frac{1}{2}\sum_{k_{s}=1}^{K}q_{Z}(Z_{s-1}=k_{s-1})\left(\boldsymbol{\theta}_{s}^{T}Q_{k_{s-1}}^{-1}\boldsymbol{\theta}_{s}-2\boldsymbol{\theta}_{s}^{T}Q_{k_{s-1}}^{-1}\boldsymbol{\theta}_{s-1}+\boldsymbol{\theta}_{s-1}^{T}Q_{k_{s-1}}^{-1}\boldsymbol{\theta}_{s-1}\right)+c$

Thus by denoting
$\bar{\sigma}^{-1}_{s}=\sqrt{\sum_{k=1}^{K}q_{Z}(Z_{s}=k)\frac{1}{\sigma_{k}^{2}}}$
and $\bar{Q}_{s}^{-1}=\sum_{k=1}^{K}q_{Z}(Z_{s}=k)Q^{-1}_{k}$, we therefore
have:

$\ln
q_{\theta}^{*}(\boldsymbol{\theta}_{1:t})=\sum_{s=1}^{t}\ln\mathcal{N}(y_{s}|\boldsymbol{\theta}^{T}_{s}\boldsymbol{x}_{s},\bar{\sigma}_{s}^{2})+\ln\mathcal{N}(\boldsymbol{\theta}_{s}|\boldsymbol{\theta}_{s-1},\bar{Q}_{s-1})+c$

∎

This last proposition shows that we are exactly in the Kalman setting where
the covariance matrix at each time step is the inverse of the average of the
inverse covariance matrix over the law of $Z$. Indeed if we define the latent
space model as :

$\text{State
:}\qquad\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}+\boldsymbol{\nu}_{t}\
,\qquad\boldsymbol{\nu}_{t}\sim\mathcal{N}(\boldsymbol{0},\overline{Q}_{t})$
(12) $\text{Space : }\qquad
y_{t}=\boldsymbol{\theta}_{t}^{T}\boldsymbol{x}_{t}+\epsilon_{t}\
,\qquad\epsilon_{t}\sim\mathcal{N}(0,\overline{\sigma}_{t})$ (13)

The log-likelihoof of the model writes as in (11).

######  Proposition 5.

For some $t>1$ and $s\in[\\![1:t]\\!]$, if we denote the
$q_{\theta}^{F}(\boldsymbol{\theta}_{s})$ the filtering approximation
of$p(\boldsymbol{\theta}_{s}|\boldsymbol{x}_{1:s},y_{1:s})$ then we have:

$q_{\theta}^{F}(\boldsymbol{\theta}_{s})=\mathcal{N}(\boldsymbol{\theta}_{s}|\boldsymbol{\hat{\theta}}_{s|s},P_{s|s})$
(14)

with the following updates :

$\displaystyle P_{s|s-1}=P_{s-1|s-1}+\overline{Q}_{s-1}$
$\displaystyle\boldsymbol{\hat{\theta}}_{s|s-1}=\boldsymbol{\hat{\theta}}_{s-1|s-1}$
$\displaystyle
P_{s|s}=P_{s|s-1}-\frac{P_{s|s-1}\boldsymbol{x}_{s}\boldsymbol{x}_{s}^{T}P_{s|s-1}}{\boldsymbol{x}_{s}^{T}P_{s|s-1}\boldsymbol{x}_{s}+\overline{\sigma}_{s}^{2}}$
$\displaystyle\boldsymbol{\hat{\theta}}_{s|s}=\boldsymbol{\hat{\theta}}_{s|s-1}-\frac{P_{s|s}}{\overline{\sigma}_{s}^{2}}\left(\boldsymbol{x}_{s}(\boldsymbol{\hat{\theta}}_{s|s-1}^{T}\boldsymbol{x}_{s}-y_{s})\right)$

The proof is immediate considering Theorem 1 and the previous remark.

###  A.2 E-$Z$ step

######  Proposition 6.

(E-$Z$ step)
For some $t>1$, for any $q_{\theta}(\theta_{1:t})$, the distribution
$q_{Z}^{*}$ which minimizes $KL(q_{\theta}\times q_{Z}||\
p_{\mathcal{F}_{t}})$ for $q_{\theta}$ fixed writes:

$\displaystyle q_{Z}^{*}(Z_{1:t})$
$\displaystyle\propto\pi(Z_{1})\mathcal{N}(\boldsymbol{\mu}_{1}|\boldsymbol{\hat{\theta}}_{1},P_{1})\mathcal{N}(y_{1}|\boldsymbol{\mu}_{1}^{T}\boldsymbol{x}_{1},\sigma_{Z_{1}}^{2})\exp\left(-\frac{1}{2}\text{Tr}\left(\left[P_{1}^{-1}+\frac{\boldsymbol{x}_{1}\boldsymbol{x}_{1}^{T}}{\sigma_{Z_{1}}^{2}}\right]V_{1}\right)\right)$
$\displaystyle\prod_{s=2}^{t}\Bigg{[}\mathcal{N}(y_{s}|\boldsymbol{\mu}_{s}^{T}\boldsymbol{x}_{s},\sigma_{Z_{s}}^{2})\mathcal{N}(\boldsymbol{\mu}_{s}|\boldsymbol{\mu}_{s-1},Q_{Z_{s-1}})\exp\left(-\frac{1}{2}\text{Tr}\left[Q_{Z_{s-1}}^{-1}(V_{s}+V_{s-1}-2W_{s-1})\right]\right)M_{Z_{s}|Z_{s-1}}\Bigg{]}$

where, for $s\in[\\![1,t]\\!]$ :
$\boldsymbol{\mu}_{s}=\mathbb{E}_{q_{\theta}(\boldsymbol{\theta}_{1:t})}[\boldsymbol{\theta}_{s}],V_{s}=\mathbb{E}_{q_{\theta}(\boldsymbol{\theta}_{1:t})}[\boldsymbol{\theta}_{s}\boldsymbol{\theta}_{s}^{T}]-\boldsymbol{\mu}_{s}\boldsymbol{\mu}_{s}^{T}$
and
$W_{s}=\mathbb{E}_{q_{\theta}(\boldsymbol{\theta}_{1:t})}[\boldsymbol{\theta}_{s}\boldsymbol{\theta}_{s+1}^{T}]-\boldsymbol{\mu}_{s}\boldsymbol{\mu}_{s+1}^{T}$
which can be rewritten as:

$q_{Z}(Z_{1:t})=\rho_{1}(Z_{1})\prod_{s=2}^{t}\rho_{s}(Z_{s})\gamma_{Z_{s}|Z_{s-1}}$
(15)

where for $s>1$:

$\displaystyle\rho_{1}(Z_{1})$
$\displaystyle\propto\pi(Z_{1})\mathcal{N}(\boldsymbol{\mu}_{1}|\boldsymbol{\hat{\theta}}_{1},P_{1})\mathcal{N}(y_{1}|\boldsymbol{\mu}_{1}^{T}\boldsymbol{x}_{1},\sigma_{Z_{1}}^{2})\exp\left(-\frac{1}{2}\text{Tr}\left(\left[P_{1}^{-1}+\frac{\boldsymbol{x}_{1}\boldsymbol{x}_{1}^{T}}{\sigma_{Z_{1}}^{2}}\right]V_{1}\right)\right)$
$\displaystyle\rho_{s}(Z_{s})$
$\displaystyle\propto\mathcal{N}(y_{s}|\boldsymbol{\mu}_{s}^{T}\boldsymbol{x}_{s},\sigma_{Z_{s}}^{2})\exp\left(-\frac{1}{2}\text{Tr}\left(\frac{\boldsymbol{x}_{s}\boldsymbol{x}_{s}^{T}}{\sigma_{Z_{s}}^{2}}V_{s}\right)\right)$
$\displaystyle\gamma_{Z_{s}|Z_{s-1}}$
$\displaystyle\propto\mathcal{N}(\boldsymbol{\mu}_{s}|\boldsymbol{\mu}_{s-1},Q_{Z_{s-1}})\exp\left(-\frac{1}{2}\text{Tr}\left[Q_{Z_{s-1}}^{-1}(V_{s}+V_{s-1}-2W_{s-1})\right]\right)\tau_{Z_{s}|Z_{s-1}}$

This last expression (15) is equivalent to a Hidden Markov Model (HMM). We
define the filtration
$\mathcal{H}_{s}=\sigma(x_{1:s},y_{1:s},\mu_{1:s},V_{1:s},W_{1:s-1})$ and
$p^{\prime}$ such that $q_{Z}(Z_{1:t})=p^{\prime}(Z_{1:t}|H_{t})$. Therefore,
by computing the standard forward-backward algorithm for HMMs, one can easily
obtain the forward $q_{Z}^{F}(Z_{s})=p^{\prime}(Z_{s}|\mathcal{H}_{s})\approx
p(Z_{s}|\mathcal{F}_{s})$ and backward
$q_{Z}^{B}(Z_{s})=p^{\prime}(Z_{s}|\mathcal{H}_{s+1:t})\approx
p(Z_{s}|\boldsymbol{x}_{s+1:t},y_{s+1:t})$ and obtain the posterior
approximation $q_{Z}(Z_{s})\approx p(Z_{s}|\mathcal{F}_{t})$ :

$q_{Z}(Z_{s})\propto q_{Z}^{F}(Z_{s})q_{Z}^{B}(Z_{s})$ (16)

######  Proposition 7.

The forward-backward recursions write for some $t>1$ and $s\in[\\![1,t]\\!]$ :

$q_{Z}^{F}(Z_{s})\propto\rho_{s}(Z_{s})\odot(q_{Z}^{F}(Z_{s-1})\gamma_{Z_{s}|Z_{s-1}})$
(17)
$q_{Z}^{B}(Z_{s})\propto\gamma_{Z_{s+1}|Z_{s}}\left(q_{Z}^{B}(Z_{s+1})\odot\rho_{s+1}(Z_{s+1})\right)$
(18)

###### Proof.

Starting from the joint distribution we have :

$\displaystyle\alpha_{s}(Z_{s})$
$\displaystyle=p(Z_{s},x_{1:s},y_{1:s},\mu_{1:s},V_{1:s},W_{1:s-1})$
$\displaystyle=\sum_{k=1}^{K}p(Z_{s},Z_{s-1}=k,y_{1:s},x_{1:s},\mu_{1:s},V_{1:s},W_{1:s-1})\qquad\qquad\text{(marginalizing
over $Z_{s-1}$)}$
$\displaystyle=\sum_{k=1}^{K}p(y_{s}|Z_{s},Z_{s-1}=k,y_{1:s-1},x_{1:s},\mu_{1:s},V_{1:s},W_{1:s-1})p(Z_{s},Z_{s-1}=k,y_{1:s-1},x_{1:s},\mu_{1:s},V_{1:s},W_{1:s-1})$
$\displaystyle=\sum_{k=1}^{K}p(y_{s}|Z_{s},Z_{s-1}=k,x_{s},\mu_{s},V_{s},W_{s-1})p(Z_{s}|Z_{s-1}=k)$
$\displaystyle
p(x_{s}|Z_{s-1}=k,x_{1:s-1},y_{1:s-1},\mu_{1:s-1},V_{1:s-1},W_{1:s-2})p(Z_{s-1}=k,x_{1:s-1},y_{1:s-1},\mu_{1:s-1},V_{1:s-1},W_{1:s-2})$
$\displaystyle\propto\sum_{k=1}^{K}\rho_{s}(Z_{s})\gamma_{Z_{s}|Z_{s-1}=k}\alpha_{s-1}(Z_{s-1}=k)$
$\displaystyle\propto\rho_{s}(Z_{s})\odot\left(\alpha_{s-1}(Z_{s-1})\gamma_{Z_{s}|Z_{s-1}}\right)$
$\displaystyle\Longrightarrow q_{Z}^{F}(Z_{s})$
$\displaystyle=p(Z_{s}|x_{1:s},y_{1:s})$
$\displaystyle\propto\rho_{s}(Z_{s})\odot(q_{Z}^{F}(Z_{s-1})\gamma_{Z_{s}|Z_{s-1}})$
$\displaystyle q_{Z}^{B}(Z_{s})$ $\displaystyle=p(y_{s+1:t},x_{s+1:t}|Z_{s})$
$\displaystyle=\sum_{k=1}^{K}p(y_{s+1:t},x_{s+1:t},Z_{s+1}=k|Z_{s})\qquad\qquad\text{(marginalizing
over $Z_{s+1}$)}$
$\displaystyle=\sum_{k=1}^{K}p(y_{s+1}|x_{s+1:t},y_{s+2:t},Z_{s+1}=k,Z_{s})p(x_{s+1:t},y_{s+2:t},Z_{s+1}=k,Z_{s})$
$\displaystyle=\sum_{k=1}^{K}p(y_{s+1}|x_{s+1},Z_{s+1}=k,Z_{s})p(Z_{s+1}=k|Z_{s},x_{s+1:t},y_{s+2:t})p(x_{s+1:t},y_{s+2:t},Z_{s+1}=k)$
$\displaystyle\propto\sum_{k=1}^{K}\rho_{s+1}(Z_{s+1}=k)q_{Z}^{B}(Z_{s+1}=k)\gamma_{Z_{s+1}=k|Z_{s}}\qquad\qquad\text{(\\`{a}
$p(x_{s+1}|.)$ pr\\`{e}s)}$
$\displaystyle\propto\gamma_{Z_{s+1}|Z_{s}}(q_{Z}^{B}(Z_{s+1})\odot\rho_{s+1}(Z_{s+1}))$

∎

###  A.3 VB Algorithm

Finally the VB algorithm can writes as follow. The parameter $n_{iter}$
corresponds to the number of iteration we need to converge and usually is of
order $5$.

Algorithm 5 VB

Require:
$q_{Z}^{(0)},(\boldsymbol{\mu_{s}^{(0)}})_{s=1}^{t},(V_{s}^{(0)})_{s=1}^{t},(W_{s}^{(0)})_{s=1}^{t-1}(\boldsymbol{x}_{s})_{s=1}^{t},(y_{s})_{s=1}^{t}$

for $k\in[\\![1,n_{iter}]\\!]$ do:

For $s\in[\\![1,t]\\!]$, compute $\bar{\sigma}_{s},\bar{Q}_{s-1}$.

For $s\in[\\![1,t]\\!]$, derive $q_{\theta}^{F}(\boldsymbol{\theta}_{s})$ with
Proposition 5.$\triangleright$ E-$\theta$ step

For $s\in[\\![1,t]\\!]$ compute $\boldsymbol{\mu}_{s},V_{s},W_{s}$.

For $s\in[\\![1,t]\\!]$, derive $q_{Z}(\boldsymbol{\theta}_{s})\propto
q_{Z}^{F}(Z_{s})q_{Z}^{B}(Z_{s})$ with Proposition 6. $\triangleright$ E-$Z$
step

end for


* Ishwaran et al. [2008] Ishwaran, H., Kogalur, U. B., Blackstone, E. H., & Lauer, M. S. (2008).  Random survival forests.  The Annals of Applied Statistics, 2, 841 -- 860. URL: https://doi.org/10.1214/08-AOAS169. doi:10.1214/08-AOAS169.
  * Jakubowski et al. [2022a] Jakubowski, J., Stanisz, P., Bobek, S., & Nalepa, G. (2022a).  Roll wear prediction in strip cold rolling with physics-informed autoencoder and counterfactual explanations.  In J. Zhexue Huang, Y. Pan, B. Hammer, M. Khurram Khan, X. Xie, L. Cui, & Y. He (Eds.), 2022 IEEE 9th International Conference on Data Science and Advanced Analytics (DSAA).  [s.l.]: IEEE.
  * Jakubowski et al. [2021] Jakubowski, J., Stanisz, P., Bobek, S., & Nalepa, G. J. (2021).  Explainable anomaly detection for hot-rolling industrial process.  In 2021 IEEE 8th International Conference on Data Science and Advanced Analytics (DSAA).  IEEE.  URL: https://doi.org/10.1109%2Fdsaa53316.2021.9564228. doi:10.1109/dsaa53316.2021.9564228.
  * Jakubowski et al. [2022b] Jakubowski, J., Stanisz, P., Bobek, S., & Nalepa, G. J. (2022b).  Anomaly detection in asset degradation process using variational autoencoder and explanations.  Sensors, 22. URL: https://www.mdpi.com/1424-8220/22/1/291. doi:10.3390/s22010291.
  * Jakubowski et al. [2022c] Jakubowski, J., Stanisz, P., Bobek, S., & Nalepa, G. J. (2022c).  Performance of explainable AI methods in asset failure prediction.  In Computational Science – ICCS 2022 (pp. 472--485).  Springer International Publishing.  URL: https://doi.org/10.1007%2F978-3-031-08760-8_40. doi:10.1007/978-3-031-08760-8_40.
  * Jamshidi et al. [2023] Jamshidi, S., Nowaczyk, S., Fanaee-T, H., & Rahat, M. (2023).  A systematic approach for tracking the evolution of xai as a field of research.  In I. Koprinska, & et al. (Eds.), Machine Learning and Principles and Practice of Knowledge Discovery in Databases (pp. 461--476).  Springer Nature Switzerland.
  * Jiang & Tieu [2004] Jiang, Z., & Tieu, A. (2004).  A 3-d finite element method analysis of cold rolling of thin strip with friction variation.  Tribology International, 37, 185--191. URL: https://www.sciencedirect.com/science/article/pii/S0301679X03000495. doi:https://doi.org/10.1016/S0301-679X(03)00049-5.  Austrib 2002.
  * Kaplan & Meier [1958] Kaplan, E. L., & Meier, P. (1958).  Nonparametric estimation from incomplete observations.  Journal of the American Statistical Association, 53, 457--481. URL: http://www.jstor.org/stable/2281868.
  * Kaufman & Rousseeuw [1987] Kaufman, L., & Rousseeuw, P. J. (1987).  Clustering by means of medoids.
  * Kaufmann & Rousseeuw [1987] Kaufmann, L., & Rousseeuw, P. (1987).  Clustering by Means of Medoids.  Data Analysis based on the L1-Norm and Related Methods, (pp. 405--416).
  * Kaur et al. [2020] Kaur, H., Nori, H., Jenkins, S., Caruana, R., Wallach, H., & Wortman Vaughan, J. (2020).  Interpreting Interpretability: Understanding Data Scientists’ Use of Interpretability Tools for Machine Learning.  In Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems (pp. 1--14).  New York, NY, USA: Association for Computing Machinery.  URL: https://doi.org/10.1145/3313831.3376219.
  * Kharal [2020] Kharal, A. (2020).  Explainable artificial intelligence based fault diagnosis and insight harvesting for steel plates manufacturing.  CoRR, abs/2008.04448. URL: https://arxiv.org/abs/2008.04448. arXiv:2008.04448.
  * Kim et al. [2016a] Kim, B., Khanna, R., & Koyejo, O. (2016a).  Examples Are Not Enough, Learn to Criticize! Criticism for Interpretability.  In Proceedings of the 30th International Conference on Neural Information Processing Systems NIPS’16 (pp. 2288--2296).  Red Hook, NY, USA: Curran Associates Inc.
  * Kim et al. [2016b] Kim, B., Khanna, R., & Koyejo, O. O. (2016b).  Examples are not enough, learn to criticize! criticism for interpretability.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, & R. Garnett (Eds.), Advances in Neural Information Processing Systems.  Curran Associates, Inc. volume 29.  URL: https://proceedings.neurips.cc/paper/2016/file/5680522b8e2bb01943234bce7bf84534-Paper.pdf.
  * Kim et al. [2016c] Kim, B., Khanna, R., & Koyejo, O. O. (2016c).  Examples are not enough, learn to criticize! criticism for interpretability.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, & R. Garnett (Eds.), Advances in Neural Information Processing Systems.  Curran Associates, Inc. volume 29.  URL: https://proceedings.neurips.cc/paper/2016/file/5680522b8e2bb01943234bce7bf84534-Paper.pdf.
  * Kleyner et al. [2017] Kleyner, A., Vasan, A., & Pecht, M. (2017).  A new application for failure prognostics--reduction of automotive electronics reliability test duration.  In Annual Conference of the Prognostics and Health Management Society.
  * Koh & Liang [2017] Koh, P. W., & Liang, P. (2017).  Understanding Black-box Predictions via Influence Functions.  URL: https://arxiv.org/abs/1703.04730. doi:10.48550/ARXIV.1703.04730.
  * Le Son et al. [2013] Le Son, K., Fouladirad, M., Barros, A., Levrat, E., & Iung, B. (2013).  Remaining useful life estimation based on stochastic deterioration models: A comparative study.  Reliability Engineering & System Safety, 112, 165--175.
  * Lei et al. [2020] Lei, Y., Yang, B., Jiang, X., Jia, F., Li, N., & Nandi, A. K. (2020).  Applications of machine learning to machine fault diagnosis: A review and roadmap.  Mechanical Systems and Signal Processing, 138, 106587.
  * Lei & Zuo [2009] Lei, Y., & Zuo, M. J. (2009).  Gear crack level identification based on weighted k nearest neighbor classification algorithm.  Mechanical Systems and Signal Processing, 23, 1535--1547.
  * Lenard [2014] Lenard, J. G. (2014).  9 - tribology.  In J. G. Lenard (Ed.), Primer on Flat Rolling (Second Edition) (pp. 193--266).  Oxford: Elsevier. (Second edition ed.).  URL: https://www.sciencedirect.com/science/article/pii/B9780080994185000093. doi:https://doi.org/10.1016/B978-0-08-099418-5.00009-3.
  * Li & Jung [2022] Li, G., & Jung, J. J. (2022).  Deep learning for anomaly detection in multivariate time series: Approaches, applications, and challenges.  Information Fusion, .
  * Li et al. [2014] Li, X.-Y., Liu, L., Kang, R., Xu, D., Sun, F., & Lee, J. (2014).  Experiments for phm: Needs, developments and challenges.  In Safety and Reliability: Methodology and Applications (pp. 605--612).  CRC Press.
  * Li et al. [2022] Li, Z., Zhu, Y., & van Leeuwen, M. (2022).  A survey on explainable anomaly detection.  arXiv preprint arXiv:2210.06959, .
  * Liao et al. [2021] Liao, Q. V., Pribić, M., Han, J., Miller, S., & Sow, D. (2021).  Question-Driven Design Process for Explainable AI User Experiences.  arXiv:2104.03483.
  * Linardatos et al. [2021] Linardatos, P., Papastefanopoulos, V., & Kotsiantis, S. B. (2021).  Explainable AI: A Review of Machine Learning Interpretability Methods.  Entropy, 23.
  * Liu et al. [2013] Liu, K., Gebraeel, N. Z., & Shi, J. (2013).  A data-level fusion model for developing composite health indices for degradation modeling and prognostic analysis.  IEEE Transactions on Automation Science and Engineering, 10, 652--664.
  * Liu et al. [2018] Liu, R., Yang, B., Zio, E., & Chen, X. (2018).  Artificial intelligence for fault diagnosis of rotating machinery: A review.  Mechanical Systems and Signal Processing, 108, 33--47.
  * Lundberg & Lee [2017] Lundberg, S. M., & Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, & R. Garnett (Eds.), Advances in Neural Information Processing Systems 30 (pp. 4765--4774).  Curran Associates, Inc.
  * Lundgren [2020] Lundgren, A. (2020).  Data-driven engine fault classification and severity estimation using residuals and data.
  * Madhikermi et al. [2019] Madhikermi, M., Malhi, A. K., & Främling, K. (2019).  Explainable artificial intelligence based heat recycler fault detection in air handling unit.  In International Workshop on Explainable, Transparent Autonomous Agents and Multi-Agent Systems (pp. 110--125).  Springer.
  * Mahajan et al. [2020] Mahajan, D., Tan, C., & Sharma, A. (2020).  Preserving causal constraints in counterfactual explanations for machine learning classifiers.  arXiv:1912.03277.
  * Markus et al. [2021] Markus, A. F., Kors, J. A., & Rijnbeek, P. R. (2021).  The role of explainability in creating trustworthy artificial intelligence for health care: A comprehensive survey of the terminology, design choices, and evaluation strategies.  Journal of Biomedical Informatics, 113, 103655. URL: http://dx.doi.org/10.1016/j.jbi.2020.103655. doi:10.1016/j.jbi.2020.103655.
  * Mashhadi et al. [2020] Mashhadi, P. S., Nowaczyk, S., & Pashami, S. (2020).  Stacked ensemble of recurrent neural networks for predicting turbocharger remaining useful life.  Applied Sciences, 10. doi:10.3390/app10010069.
  * Matyas et al. [2017] Matyas, K., Nemeth, T., Kovacs, K., & Glawar, R. (2017).  A procedural approach for realizing prescriptive maintenance planning in manufacturing industries.  CIRP Annals, 66, 461--464. URL: https://www.sciencedirect.com/science/article/pii/S0007850617300070. doi:10.1016/j.cirp.2017.04.007.
  * May & Sawodny [2021] May, M. P., & Sawodny, O. (2021).  Multi-class sensor fault isolation for vehicle dynamics with decision trees.  In IECON 2021–47th Annual Conference of the IEEE Industrial Electronics Society (pp. 1--6).  IEEE.
  * Meissner et al. [2021a] Meissner, R., Meyer, H., & Wicke, K. (2021a).  Concept and Economic Evaluation of Prescriptive Maintenance Strategies for an Automated Condition Monitoring System.  International Journal of Prognostics and Health Management, 12. URL: https://papers.phmsociety.org/index.php/ijphm/article/view/2911. doi:10.36001/ijphm.2021.v12i3.2911.  Number: 3.
  * Meissner et al. [2021b] Meissner, R., Rahn, A., & Wicke, K. (2021b).  Developing prescriptive maintenance strategies in the aviation industry based on a discrete-event simulation framework for post-prognostics decision making.  Reliability Engineering & System Safety, 214, 107812. URL: https://www.sciencedirect.com/science/article/pii/S0951832021003331. doi:10.1016/j.ress.2021.107812.
  * Mey & Neufeld [2022] Mey, O., & Neufeld, D. (2022).  Explainable ai algorithms for vibration data-based fault detection: Use case-adadpted methods and critical evaluation.  arXiv preprint arXiv:2207.10732, .
  * Miguéis et al. [2022] Miguéis, V. L., Borges, J. L. et al. (2022).  Automatic root cause analysis in manufacturing: an overview & conceptualization.  Journal of Intelligent Manufacturing, (pp. 1--18).
  * Miller et al. [2021] Miller, A. C., Foti, N. J., & Fox, E. B. (2021).  Breiman’s two cultures: You don’t have to choose sides.  arXiv:2104.12219.
  * Miller [2019] Miller, T. (2019).  Explanation in artificial intelligence: Insights from the social sciences.  Artificial Intelligence, 267, 1--38. URL: https://www.sciencedirect.com/science/article/pii/S0004370218305988. doi:https://doi.org/10.1016/j.artint.2018.07.007.
  * Molnar [2019] Molnar, C. (2019).  Interpretable Machine Learning.
  * Molnar [2020] Molnar, C. (2020).  Interpretable machine learning.  Lulu.com.
  * Moraffah et al. [2020] Moraffah, R., Karami, M., Guo, R., Raglin, A. J., & Liu, H. (2020).  Causal interpretability for machine learning - problems, methods and evaluation.  ACM SIGKDD Explorations Newsletter, 22, 18 -- 33.
  * Mothilal et al. [2020] Mothilal, R. K., Sharma, A., & Tan, C. (2020).  Explaining machine learning classifiers through diverse counterfactual explanations.  In Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency (pp. 607--617).
  * Mueller et al. [2019] Mueller, S. T., Hoffman, R. R., Clancey, W., Emrey, A., & Klein, G. (2019).  Explanation in Human-AI Systems: A Literature Meta-Review, Synopsis of Key Ideas and Publications, and Bibliography for Explainable AI.  arXiv:1902.01876.
  * Munir et al. [2018] Munir, M., Siddiqui, S. A., Dengel, A., & Ahmed, S. (2018).  Deepant: A deep learning approach for unsupervised anomaly detection in time series.  Ieee Access, 7, 1991--2005.
  * Myrtakis et al. [2021] Myrtakis, N., Christophides, V., & Simon, E. (2021).  A comparative evaluation of anomaly explanation algorithms.  In 24th International Conference on Extending Database Technology (EDBT’2021).
  * Nectoux et al. [2012] Nectoux, P., Gouriveau, R., Medjaher, K., Ramasso, E., Chebel-Morello, B., Zerhouni, N., & Varnier, C. (2012).  Pronostia: An experimental platform for bearings accelerated degradation tests.  In IEEE International Conference on Prognostics and Health Management, PHM’12. (pp. 1--8).  IEEE Catalog Number: CPF12PHM-CDR.
  * Nemeth et al. [2018] Nemeth, T., Ansari, F., Sihn, W., Haslhofer, B., & Schindler, A. (2018).  PriMa-X: A reference model for realizing prescriptive maintenance and assessing its maturity enhanced by machine learning.  Procedia CIRP, 72, 1039--1044. URL: https://www.sciencedirect.com/science/article/pii/S2212827118304566. doi:10.1016/j.procir.2018.03.280.
  * Nori et al. [2019] Nori, H., Jenkins, S., Koch, P., & Caruana, R. (2019).  InterpretML: A Unified Framework for Machine Learning Interpretability.  URL: https://arxiv.org/abs/1909.09223. doi:10.48550/ARXIV.1909.09223.
  * Nowaczyk et al. [2020] Nowaczyk, S., Rögnvaldsson, T., Fan, Y., & Calikus, E. (2020).  Towards autonomous knowledge creation from big data in smart cities.  In J. C. Augusto (Ed.), Handbook of Smart Cities (pp. 1--35).  Cham: Springer International Publishing.  doi:10.1007/978-3-030-15145-4_38-1.
  * Pang et al. [2021] Pang, G., Shen, C., Cao, L., & Hengel, A. V. D. (2021).  Deep learning for anomaly detection: A review.  ACM Computing Surveys (CSUR), 54, 1--38.
  * Pang et al. [2020] Pang, Y., He, Q., Jiang, G., & Xie, P. (2020).  Spatio-temporal fusion neural network for multi-class fault diagnosis of wind turbines based on scada data.  Renewable Energy, 161, 510--524.
  * Pawelczyk et al. [2021] Pawelczyk, M., Bielawski, S., van den Heuvel, J., Richter, T., & Kasneci, G. (2021).  CARLA: A Python Library to Benchmark Algorithmic Recourse and Counterfactual Explanation Algorithms.  URL: https://arxiv.org/abs/2108.00783. doi:10.48550/ARXIV.2108.00783.
  * Pearl & Mackenzie [2018] Pearl, J., & Mackenzie, D. (2018).  The Book of Why: The New Science of Cause and Effect.  (1st ed.).  USA: Basic Books, Inc.
  * Piltan & Kim [2021] Piltan, F., & Kim, J.-M. (2021).  Crack size identification for bearings using an adaptive digital twin.  Sensors, 21, 5009\.
  * Prytz et al. [2015] Prytz, R., Nowaczyk, S., Rögnvaldsson, T., & Byttner, S. (2015).  Predicting the need for vehicle compressor repairs using maintenance records and logged vehicle data.  Engineering Applications of Artificial Intelligence, 41, 139--150.
  * Qiu et al. [2019] Qiu, G., Gu, Y., & Cai, Q. (2019).  A deep convolutional neural networks model for intelligent fault diagnosis of a gearbox under different operational conditions.  Measurement, 145, 94--107.
  * Rajaoarisoa & Sayed-Mouchaweh [2017] Rajaoarisoa, L., & Sayed-Mouchaweh, M. (2017).  Adaptive online fault diagnosis of manufacturing systems based on devs formalism.  IFAC-PapersOnLine, 50, 6825--6830. URL: https://www.sciencedirect.com/science/article/pii/S2405896317317081. doi:https://doi.org/10.1016/j.ifacol.2017.08.1202.  20th IFAC World Congress.
  * Randriarison et al. [2023] Randriarison, J., Rajaoarisoa, L., & Sayed-Mouchaweh, M. (2023).  Faults explanation based on a machine learning model for predictive maintenance purposes.  In Proceedings of the 7th edition in the series of the International Conference on Control, Automation and Diagnosis, (p. 1).
  * Ranganath et al. [2016] Ranganath, R., Perotte, A., Elhadad, N., & Blei, D. (2016).  Deep survival analysis.  In F. Doshi-Velez, J. Fackler, D. Kale, B. Wallace, & J. Wiens (Eds.), Proceedings of the 1st Machine Learning for Healthcare Conference (pp. 101--114).  Northeastern University, Boston, MA, USA: PMLR volume 56 of Proceedings of Machine Learning Research.  URL: https://proceedings.mlr.press/v56/Ranganath16.html.
  * Remadna et al. [2020] Remadna, I., Terrissa, S. L., Zemouri, R., Ayad, S., & Zerhouni, N. (2020).  Leveraging the power of the combination of cnn and bi-directional lstm networks for aircraft engine rul estimation.  In 2020 Prognostics and Health Management Conference (PHM-Besançon) (pp. 116--121).  IEEE.
  * Ribeiro et al. [2016] Ribeiro, M. T., Singh, S., & Guestrin, C. (2016).  "why should i trust you?": Explaining the predictions of any classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining KDD ’16 (p. 1135–1144).  New York, NY, USA: Association for Computing Machinery.  doi:10.1145/2939672.2939778.
  * Ribeiro et al. [2023] Ribeiro, R. P., Mastelini, S. M., Davari, N., Aminian, E., Veloso, B., & Gama, J. (2023).  Online anomaly explanation: A case study on predictive maintenance.  In Machine Learning and Principles and Practice of Knowledge Discovery in Databases: International Workshops of ECML PKDD 2022, Grenoble, France, September 19–23, 2022, Proceedings, Part II (pp. 383--399).  Springer.
  * Rigamonti et al. [2016] Rigamonti, M. M., Baraldi, P., Zio, E. et al. (2016).  Echo state network for the remaining useful life prediction of a turbofan engine.  In annual conference of the prognostics and health management society 2015 (pp. 255--270).
  * Rögnvaldsson et al. [2018] Rögnvaldsson, T., Nowaczyk, S., Byttner, S., Prytz, R., & Svensson, M. (2018).  Self-monitoring for maintenance of vehicle fleets.  Data mining and knowledge discovery, 32, 344--384.
  * Rudin [2019] Rudin, C. (2019).  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nature Machine Intelligence, 1, 206--215. doi:10.1038/s42256-019-0048-x.
  * Sayed-Mouchaweh & Rajaoarisoa [2022] Sayed-Mouchaweh, M., & Rajaoarisoa, L. (2022).  Explainable decision support tool for iot predictive maintenance within the context of industry 4.0.  In 2022 21st IEEE International Conference on Machine Learning and Applications (ICMLA) (pp. 1492--1497).  doi:10.1109/ICMLA55696.2022.00234.
  * Schröder et al. [2022] Schröder, L., Dimitrov, N. K., Verelst, D. R., & Sørensen, J. A. (2022).  Using transfer learning to build physics-informed machine learning models for improved wind farm monitoring.  Energies, 15. URL: https://www.mdpi.com/1996-1073/15/2/558. doi:10.3390/en15020558.
  * Sejr & Schneider-Kamp [2021] Sejr, J. H., & Schneider-Kamp, A. (2021).  Explainable outlier detection: What, for whom and why?  Machine Learning with Applications, 6, 100172.
  * Semenova et al. [2021] Semenova, L., Rudin, C., & Parr, R. (2021).  A study in Rashomon curves and volumes: A new perspective on generalization and model simplicity in machine learning.  arXiv:1908.01755.
  * Senoner et al. [2022] Senoner, J., Netland, T., & Feuerriegel, S. (2022).  Using explainable artificial intelligence to improve process quality: Evidence from semiconductor manufacturing.  Management Science, 68, 5704--5723. URL: https://doi.org/10.1287/mnsc.2021.4190. doi:10.1287/mnsc.2021.4190.
  * Serdio et al. [2015] Serdio, F., Lughofer, E., Pichler, K., Pichler, M., Buchegger, T., & Efendic, H. (2015).  Fuzzy fault isolation using gradient information and quality criteria from system identification models.  Information Sciences, 316, 18--39.
  * Serradilla et al. [2021] Serradilla, O., Zugasti, E., Ramirez de Okariz, J., Rodriguez, J., & Zurutuza, U. (2021).  Adaptable and explainable predictive maintenance: Semi-supervised deep learning for anomaly detection and diagnosis in press machine data.  Applied Sciences, 11, 7376.
  * Si et al. [2011] Si, X.-S., Wang, W., Hu, C.-H., & Zhou, D.-H. (2011).  Remaining useful life estimation--a review on the statistical data driven approaches.  European journal of operational research, 213, 1--14.
  * Siddiqui et al. [2019] Siddiqui, M. A., Fern, A., Dietterich, T. G., & Wong, W.-K. (2019).  Sequential feature explanations for anomaly detection.  ACM Transactions on Knowledge Discovery from Data (TKDD), 13, 1--22.
  * Srinivasan et al. [2021] Srinivasan, S., Arjunan, P., Jin, B., Sangiovanni-Vincentelli, A. L., Sultan, Z., & Poolla, K. (2021).  Explainable ai for chiller fault-detection systems: gaining human trust.  Computer, 54, 60--68.
  * Steenwinckel et al. [2021] Steenwinckel, B., De Paepe, D., Vanden Hautte, S., Heyvaert, P., Bentefrit, M., Moens, P., Dimou, A., Van Den Bossche, B., De Turck, F., Van Hoecke, S., & Ongenae, F. (2021).  Flags: A methodology for adaptive anomaly detection and root cause analysis on sensor data streams by fusing expert knowledge with machine learning.  Future Generation Computer Systems, 116, 30--48. URL: https://www.sciencedirect.com/science/article/pii/S0167739X20329927. doi:https://doi.org/10.1016/j.future.2020.10.015.
  * Szelążek et al. [2020] Szelążek, M., Bobek, S., Gonzalez-Pardo, A., & Nalepa, G. J. (2020).  Towards the modeling of the hot rolling industrial process. preliminary results.  In C. Analide, P. Novais, D. Camacho, & H. Yin (Eds.), Intelligent Data Engineering and Automated Learning – IDEAL 2020 (pp. 385--396).  Cham: Springer International Publishing.
  * Tercan & Meisen [2022] Tercan, H., & Meisen, T. (2022).  Machine learning and deep learning based predictive quality in manufacturing: a systematic review.  Journal of Intelligent Manufacturing, 33, 1879--1905. URL: https://doi.org/10.1007/s10845-022-01963-8. doi:10.1007/s10845-022-01963-8.
  * Trindade et al. [2021] Trindade, J., Vinagre, J., Fernandes, K., Paiva, N., & Jorge, A. (2021).  Partially Monotonic Learning for Neural Networks.  In P. H. Abreu, P. P. Rodrigues, A. Fernández, & J. Gama (Eds.), Advances in Intelligent Data Analysis XIX (pp. 12--23).  Cham: Springer International Publishing.
  * Veloso et al. [2023] Veloso, B., Ribeiro, R., Gama, J., & Pereira, P. (2023).  The metropt dataset for predictive maintenance.  Scientific data, to appear, 1--10.
  * Venkatasubramanian et al. [2003] Venkatasubramanian, V., Rengaswamy, R., Yin, K., & Kavuri, S. N. (2003).  A review of process fault detection and diagnosis: Part i: Quantitative model-based methods.  Computers & chemical engineering, 27, 293--311.
  * Vilone & Longo [2021] Vilone, G., & Longo, L. (2021).  Notions of explainability and evaluation approaches for explainable artificial intelligence.  Information Fusion, 76, 89--106. URL: https://www.sciencedirect.com/science/article/pii/S1566253521001093. doi:https://doi.org/10.1016/j.inffus.2021.05.009.
  * Wachter et al. [2017] Wachter, S., Mittelstadt, B., & Russell, C. (2017).  Counterfactual explanations without opening the black box: Automated decisions and the gdpr, .  URL: https://arxiv.org/abs/1711.00399. doi:10.48550/ARXIV.1711.00399.
  * Wang et al. [2019] Wang, G., Ledwoch, A., Hasani, R. M., Grosu, R., & Brintrup, A. (2019).  A generative neural network model for the quality prediction of work in progress products.  Applied Soft Computing, 85, 105683. URL: https://www.sciencedirect.com/science/article/pii/S1568494619304648. doi:https://doi.org/10.1016/j.asoc.2019.105683.
  * Wang [2010] Wang, T. (2010).  Trajectory similarity based prediction for remaining useful life estimation.  Ph.D. thesis University of Cincinnati.
  * Wang et al. [2008] Wang, T., Yu, J., Siegel, D., & Lee, J. (2008).  A similarity-based prognostics approach for remaining useful life estimation of engineered systems.  In Prognostics and Health Management, 2008. PHM 2008. International Conference on (pp. 1--6).  IEEE.
  * Willard et al. [2022] Willard, J., Jia, X., Xu, S., Steinbach, M., & Kumar, V. (2022).  Integrating scientific knowledge with machine learning for engineering and environmental systems.  ACM Computing Surveys, . doi:10.1145/3514228.
  * Zhang et al. [2022] Zhang, T., Chen, J., Li, F., Zhang, K., Lv, H., He, S., & Xu, E. (2022).  Intelligent fault diagnosis of machines with small & imbalanced data: A state-of-the-art review and possible extensions.  ISA transactions, 119, 152--171.
  * Zhang et al. [2021] Zhang, X., Fujiwara, T., Chandrasegaran, S., Brundage, M. P., Sexton, T., Dima, A., & Ma, K.-L. (2021).  A visual analytics approach for the diagnosis of heterogeneous and multidimensional machine maintenance data.  In 2021 IEEE 14th Pacific Visualization Symposium (PacificVis) (pp. 196--205).  doi:10.1109/PacificVis52677.2021.00033.
  * Zhao et al. [2017] Zhao, Z., Liang, B., Wang, X., & Lu, W. (2017).  Remaining useful life prediction of aircraft engine based on degradation pattern learning.  Reliability Engineering & System Safety, 164, 74--83.
  * Zheng et al. [2017] Zheng, S., Ristovski, K., Farahat, A., & Gupta, C. (2017).  Long short-term memory network for remaining useful life estimation.  In 2017 IEEE International Conference on Prognostics and Health Management (ICPHM) (pp. 88--95).  IEEE.