captionUnknown document class

Ben-Gurion University of the Negev

The Faculty of Natural Sciences

The Department of Computer Science

Hierarchical Time Series Forecasting with Bayesian Modeling

Thesis submitted in partial fulfillment of the requirements

for the Master of Sciences degree

Gal Elgavish

Under the supervision of Prof. Eyal Shimony and Dr. David Tolpin

August 2023

Hierarchical Time Series Forecasting with Bayesian Modeling

Gal Elgavish

Master of Sciences Thesis

Ben-Gurion University of the Negev

2023

Abstract

We encounter time series data in many domains such as finance, physics,
business, and weather. One of the main tasks of time series analysis, one that
helps to take informed decisions under uncertainty, is forecasting. Time
series are often hierarchically structured, e.g., a company sales might be
broken down into different regions, and each region into different stores. In
some cases the number of series in the hierarchy is too big to fit in a single
model to produce forecasts in relevant time, and a decentralized approach is
beneficial.

One way to do this is to train independent forecasting models for each series
and for some summary statistics series implied by the hierarchy (e.g. the sum
of all series) and to pass those models to a reconciliation algorithm to
improve those forecasts by sharing information between the series.

In this work we focus on the reconciliation step, and propose a method to do
so from a Bayesian perspective — Bayesian forecast reconciliation. We also
define the common case of linear Gaussian reconciliation, where the forecasts
are Gaussian and the hierarchy has linear structure, and show that we can
compute reconciliation in closed form. We evaluate these methods on synthetic
and real data sets, and compare them to other work in this field.

## Acknowledgements

I would like to thank Dr. David Tolpin who not only supervised me in this
research with great patience and attention, but also taught me how to learn
new stuff and how to approach problems.

I will also like to thank Prof. Eyal Shimony who made this research possible
by giving great advices, comments and challenged me to improve my work.

###### Contents

  1. 1 Introduction
  2. 2 Background
    1. 2.1 Time Series Forecasting
      1. 2.1.1 ARIMA
      2. 2.1.2 State-Space Models
      3. 2.1.3 Gaussian Process Regression
    2. 2.2 Bayesian Graphical Models
      1. 2.2.1 Bayesian Machine Learning
      2. 2.2.2 Probabilistic Graphical Models
      3. 2.2.3 Soft Evidence
    3. 2.3 Other Related Work
      1. 2.3.1 Forecast Reconciliation
  3. 3 Problem Statement
  4. 4 Theoretical Analysis
    1. 4.1 Generative Model
      1. 4.1.1 Linear Gaussian Reconciliation
  5. 5 Empirical Evaluation
    1. 5.1 Metrics
    2. 5.2 Synthetic Data Set
      1. 5.2.1 Data — Experiment A
      2. 5.2.2 Results - Experiment A
      3. 5.2.3 Data — Experiment B
      4. 5.2.4 Results - Experiment B
    3. 5.3 Real Data Set
      1. 5.3.1 Data
      2. 5.3.2 Results
  6. 6 Conclusion

###### List of Figures

  1. 2.1 A network representing the graphical model described in the example in 2.2.3. Undirected snake edge between $e$ and $C$ denotes soft evidence about $C$ upon observing $e$.
  2. 3.1 Example of a system with four subsystems belong to two groups.
  3. 4.1 A forecasting model. The latent state $s$ is observed through observations $y$. A forecasted parameter of interest $x$ depends on $s$.
  4. 4.2 Bayesian network for forecast reconciliation. The undirected snake edge between $u$ and $y$ denotes soft evidence about $u$ upon observing $y$.
  5. (a) full network
  6. (b) relevant subnetwork
  7. 5.1 The synthetic data set for Experiment 5.2.1. The four systems $\\{x^{(n)}\\}_{n=1}^{4}$ are connected to the same parent $u$.
  8. 5.2 Example of a data set and its base forecasts for $\sigma_{\epsilon}=\\{0.1,0.2,0.5\\}$. Black x’s denote training data (observations $\boldsymbol{y}^{(n)}$), green x’s denote test data ($\boldsymbol{x}^{(n)}$), blue lines and shaded blue areas denote the mean and 95% confidence interval of the base forecasts, respectively.
  9. 5.3 Example of a data set and its ”shifted” base forecasts for $\sigma_{\epsilon}=\\{0.1,0.2,0.5\\}$. Black x’s denote training data (observations $\boldsymbol{y}^{(n)}$), green x’s denote test data ($\boldsymbol{x}^{(n)}$), blue lines and shaded blue areas denote the mean and 95% confidence interval of the base forecasts, respectively.
  10. 5.4 The first level of the hierarchy of the dataset. Top graph: the number of visitors by month to each of the seven states of Australia. Bottom graph: The sum of the series in the top graph.

###### List of Tables

  1. 5.1 Results for the experiment on Data Set 5.2.1.
  2. 5.2 Results for the experiment on Data Set 5.2.3.
  3. 5.3 Results for the experiment on Data Set of Section 5.3.1.

##  Chapter 1 Introduction

###

Time series data appears in many domains, such as: finance, physics, business,
weather, and basically every field that has temporal measurements. In many of
these fields decisions must be taken today based on what will happen in the
future, which makes Time Series Forecasting an important task in helping
stakeholders to make informed decisions under uncertainty.

Different time series forecasting problems might have different
characteristics:

  * •

Dimensionality of the time series are either univariate or multivariate.

  * •

Time stamps may be regularly or irregularly sampled.

  * •

The forecast horizon may be either short term (nowcasting [1]) or long term.

  * •

Time series might have seasonality, trend, auto-correlation as well as other
features.

  * •

Forecasts might be point-forecasts or probabilistic forecasts.

Methods for time series forecasting treat the series as a stochastic process.
Classical methods of forecasting include auto-regressive integrated moving
average (ARIMA) models [2] and Kalman filtering [3], while more modern methods
include Gaussian processes [4, 5] and deep neural networks [6, 7]. We will
elaborate on some techniques in Chapter 2.

###

There are two types of multivariate time series: heterogeneous and
homogeneous. An example of a heterogeneous time series is, for example,
medical monitoring — blood pressure, heart rate, saturation etc. The time
series are related but they are of different natures, and the number of
components is small. Examples of a homogeneous time series are stock prices,
exchange rates, precipitation levels by area and so on. All of the components
have the same ‘nature’, and the number of components can be very large. While
for relatively low-dimensional heterogeneous (and homogeneous) time series,
methods like Gaussian processes and deep neural networks can be used, for very
high dimensional homogeneous time series we often cannot predict all of the
time series together, and for reasons of efficiency and resource distribution
we want to predict each component separately and then ‘consolidate’ the
predictions. One approach to that is hierarchical time series forecasting.
Some examples of hierarchical time series and their unique characteristics
are:

###### Demand Forecasting

Trying to predict customer demand in order to optimize supply chains. The
forecasts help answering business questions such as: production planning,
inventory management, entering a new market.

Characteristics: series comprised of non-negative integers (number of units
for each product); forecast horizon might be both for short and long term;
time series might have seasonality; series might be irregularly sampled, i.e.,
not have historical data for the entire period, for example, in a grocery
store not all fruits are being sold throughout the year.

The hierarchy of these forecasts can be defined by different possible
groupings: (a) the different regions the products are sold; (b) the different
customer segments the products are sold to; (c) the different categories each
product is associated with.

###### Stock Prices

Trying to predict stock prices in order to decide on a trading strategy that
will optimize returns.

Characteristics: series are continuous (stock price); forecast horizon might
be both for short and long term; time series usually have seasonality and
trend; prices can be affected by psychology, i.e. forecasts themselves might
change investors’ behavior.

The hierarchy of these forecasts can be defined by the different indices each
stock is associated with (S&P500, FTSE100, etc.). A stock index (see [8]) is a
financial instrument that measures a stock market by grouping multiple stock
prices into one number. The aggregation of the different prices into one
number may be different between indices, e.g., it can be a simple average or
it can be weighted average on the different market cap (total value) of each
stock.

###### Electricity Demand

In order to be efficient in managing a power grid, it is vital to be able to
forecast electricity demand. The increasing adoption of smart meters provides
more data points at the customer level that can increase forecasts accuracy.

Characteristics: series are continuous; forecast horizon might be both for
short and long term; time series usually have daily, weekly, and yearly
seasonality.

The hierarchy of these forecasts can be defined by geographical regions and
customer segments (private, business, etc.). The hierarchy usually has many,
could be even millions, of series at the bottom level which might affect
computational complexity.

###

When trying to forecast a hierarchy of many time series simultaneously, and to
update forecasts in real time with every new data point that arrives, an issue
of computational complexity arises. In this case, a single forecasting model
on a single machine may take too long to train and produce forecasts, and one
might choose to use a technique such as Federated Learning [9] which takes a
decentralized approach, i.e., distribute the training to multiple machines,
each train on different series independently and then pass the trained models
to the central machine to produce forecasts. Producing forecasts for each
series independently is of course not the best solution since the information,
imposed by the hierarchy and the relation between the series, is not taken
into account.

We can see that, for example, in the electricity demand use case, where we
have both spatial and semantic hierarchical relation between time series
(users in the grid), sharing information between time series should improve
individual forecasts. Training that many series (might be millions) in a
single-model-single-machine is not tractable in relevant time, and a
distributed approach is needed. In the case of distributed approach, an
algorithm for sharing information between the forecasts would be beneficial.

The problem we will tackle in this thesis is the following: we are given
independent forecasts for each time series in a hierarchy that were (probably)
produced by distributed machines; we need to develop a method to improve those
forecasts (in a post-processing step) by sharing information between the
series using knowledge imposed by the hierarchy.

This problem is an active field of research and, very often, this post-
processing method of improving forecasts is referred to as reconciliation. In
many works (such as [10]) the goal of this reconciliation is to assure
coherency of the forecasts, i.e., the forecasts should obey the hierarchical
structure. For example, in the demand forecasting use case, the forecasts for
all stores in a region should add up to the forecast of the aggregated time
series of the region. As stated above, the problem we are trying to solve is
related but different. We are trying to improve each of the individual
forecasts as opposed to assuring coherency.

###

The next Chapter 2 reviews the relevant background needed for the scope of
this work. Chapter 3 formally defines the problem of hierarchical time series
forecasting. Chapter 4 presents the new theory we developed to deal with the
problem. Chapter 5 presents empirical results from experiments we conducted.

##  Chapter 2 Background

###  2.1 Time Series Forecasting

Time series is any data that is indexed by time, it may be univariate or
multivariate. The main tasks of time series analysis are smoothing, filtering
and forecasting. Smoothing is the task of approximating a state in time, based
on all other past and future states. Filtering is the task of “predicting the
present”, as in estimating the current not-noisy state, based on all previous
and current state. Forecasting is predicting the future, based on the past.

Forecasting can be divided into two types: point forecasting and probabilistic
or uncertainty forecasting. Point forecasts result in a single value that
might be univariate or multivariate, while probabilistic forecasts result in a
probability distribution that quantifies uncertainty.

An overview of methods for time series forecasting is provided below.

####  2.1.1 ARIMA

Autoregressive integrated moving average (ARIMA) is a method for analyzing
time series. It consists of three different models: autoregressive (AR),
differencing (integration) (I), and moving average (MA). ARIMA model is
parameterized by the order term of each of its “sub” models, i.e.,
ARIMA($p,d,q$) is the equivalent of combining: AR($p$), I($d$) and MA($q$).
[2] introduced a method for modeling time series with ARIMA models. Details of
these models appear below:

#####  2.1.1.1 AR

The autoregressive model of order $p$ (AR($p$)) is defined as

$y_{t}=\mu+\sum_{i=1}^{p}\phi_{i}y_{t-i}+\epsilon_{t}$

where $\mu$ is a constant term, $\\{\phi_{i}\\}_{1}^{p}$ are the
autoregressive parameters and $\epsilon_{t}$ is a white noise term for state
$t$.

#####  2.1.1.2 I

Differencing is a transformation applied to a non-stationary time series in
order to make it stationary in the mean sense, i.e.:

$mean_{X}(t)=mean_{X}(t+\tau)\quad\forall\tau\in\mathbb{R}$

For example, an I($1$) differencing is

$y^{\prime}_{t}=y_{t}-y_{t-1}$

and an I($2$) differencing is

$y^{\prime}_{t}=y_{t}-2y_{t-1}+y_{t-2}$

.

#####  2.1.1.3 MA

The moving average model of order $q$ (MA($q$)) is defined as

$y_{t}=\mu+\sum_{i=1}^{q}\theta_{i}\epsilon_{t-i}+\epsilon_{t}$

where $\mu$ is the mean of the series, $\\{\theta_{i}\\}_{1}^{q}$ are the
moving average parameters and $\epsilon_{t}$ is a white noise term for state
$t$.

#####  2.1.1.4 Full Model

The full ARIMA($p,d,q$) is given by combining all models together, so, for
example, an ARIMA($1,1,2$) is given by:

$y_{t}=\mu+y_{t-1}+\phi_{1}(y_{t-1}-y_{t-2})+\theta_{1}\epsilon_{t-1}+\theta_{2}\epsilon_{t-2}+\epsilon_{t}$

Parameters $\\{\phi_{i}\\}_{1}^{p},\\{\theta_{i}\\}_{1}^{q},\mu$ are estimated
by maximum likelihood estimation (MLE). Parameters $p,d,q$ are usually
optimized by comparing an information criterion such as Akaike information
criterion (AIC) or Bayesian information criterion (BIC).

#####

ARIMA models try to deal with non-stationary time series by the differencing
step, but some trends cannot be removed with differencing. Another approach to
deal with non-stationary time series is to use state-space models (Section
2.1.2).

####  2.1.2 State-Space Models

State-space models generally attempt to describe a system by assuming that
there is a latent state with time-varying dynamic relationship. The
characteristics of the system is that:

  * •

It has a Markov property, i.e., the state of the system in time $t$ is enough
to infer the state in future time $t+1$.

  * •

We cannot observe the latent state itself other than through a noisy
observation.

This leads to two main equations describing a state-space model — the
transition equation which describes how the latent state $\alpha_{t}$ evolves
through time:

$\alpha_{t+1}=f_{t}(\alpha_{t})$

and the observation equation which relates between observed data $y_{t}$ and
latent state:

$y_{t}=g_{t}(\alpha_{t})+\epsilon_{t}\quad\quad\epsilon_{t}\sim\mathcal{D}_{t}$

As with ARIMA (2.1.1) models, parameters of state space models can be
estimated with MLE, but can also be estimated with Bayesian inference (see
Section 2.2.1).

####  2.1.3 Gaussian Process Regression

Gaussian processes (GPs) are models suitable for both regression and
classification tasks. Since time series forecasting is a regression task,
i.e., we are predicting values based on time-stamped input, GPs are a great
tool for time series forecasting. Since this is the subject of this thesis, we
continue the background and intuition on regression problems. One can get more
information on GPs for classification (and for regression too) in [4].

In regression problems we are given a set of (possibly noisy) observations:

$\mathcal{D}=\\{\boldsymbol{x}_{t}\in\mathbb{R}^{D},y_{t}\\}_{t=1}^{T}$

or in matrix form:

$X\in\mathbb{R}^{T\times D},\boldsymbol{y}\in\mathbb{R}^{T}$

where we assume the observation model is:

$y_{t}=f(\boldsymbol{x}_{t})+\epsilon_{t},\quad\epsilon_{t}\sim\mathcal{N}(0,\sigma^{2})$

$\epsilon_{t}$ is random noise i.i.d. normally distributed with variance
$\sigma^{2}$. We wish to infer the function at new positions $X_{*}$ which we
denote as $\boldsymbol{f}_{*}\in\mathbb{R}^{T}$, or in symbols:

$\boldsymbol{f}_{*}\mid\boldsymbol{y};X,X_{*}$

Remarks:

  * •

$y_{t}$ is scalar but of course $\boldsymbol{y}_{t}$ can be multivariate, in
which case we should use what is called Multi-task Gaussian process (MTGP).
For more information about MTGP see [11].

  * •

$\epsilon$ is Gaussian distributed according to the assumption of a GP. It is
possible to have a different noise model, not examined in this work.

  * •

we define $\epsilon_{t}$ as homoscedastic noise but it is possible to define a
model with heteroscedastic noise, i.e., $\epsilon_{t}$ is parameterized by
$\sigma_{t}$ that can change over time. This is not crucial for the
understanding of this work and will make derivations more complicated so we
stick with homoscedastic noise.

Returning to GPs, a GP is a collection of random variables, any finite number
of which have a joint Gaussian distribution. A GP is fully defined by its mean
function $m(\boldsymbol{x})$ and covariance function (kernel function)
$K(\boldsymbol{x},\boldsymbol{x}^{\prime})$:

$\displaystyle m(\boldsymbol{x})$
$\displaystyle=\mathbb{E}[f(\boldsymbol{x})]$ $\displaystyle
K(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=\mathbb{E}[(f(\boldsymbol{x})-m(\boldsymbol{x}))(f(\boldsymbol{x}^{\prime})-m(\boldsymbol{x}^{\prime}))]$

and is written:

$f(\boldsymbol{x})\sim\mathcal{GP}(m(\boldsymbol{x}),K(\boldsymbol{x},\boldsymbol{x}^{\prime}))$

i.e., the function $f$ is GP distributed with mean function $m$ and covariance
function $K$. With this explanation a GP can also be thought of as a
distribution over functions. Since we are dealing with normal distributions it
must hold that the matrix $K$ which elements are given by:
$K_{i,j}=K(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$ is a symmetric and positive
semi-definite matrix, and thus, the kernel
$K(\boldsymbol{x},\boldsymbol{x}^{\prime})$ must be a symmetric and positive
semi-definite kernel.

Returning to the regression problem described here, the full GP model relevant
to the problem is defined as:

$\begin{bmatrix}\boldsymbol{y}\\\
\boldsymbol{f}_{*}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(X)\\\
m(X_{*})\end{bmatrix},\begin{bmatrix}K(X,X)+\sigma^{2}I&K(X,X_{*})\\\
K(X_{*},X)&K(X_{*},X_{*})\end{bmatrix}\right)$

i.e. $\boldsymbol{y}$ and $\boldsymbol{f}_{*}$ are jointly Gaussian
distributed. This is of course after we defined functions $m$ and $K$, which
will be examined in Section 2.1.3.1. We can then use the rules of conditioning
in a multivariate normal distribution to get:

$\displaystyle\boldsymbol{f}_{*}\mid\boldsymbol{y};X,X_{*}\sim\mathcal{N}($
$\displaystyle
m(X_{*})+K(X_{*},X)[K(X,X)+\sigma^{2}I]^{-1}(\boldsymbol{y}-m(X)),$
$\displaystyle K(X_{*},X_{*})-K(X_{*},X)[K(X,X)+\sigma^{2}I]^{-1}K(X,X_{*}))$

Inference in GPs is very simple and requires only linear algebra operations.
The downside is that the operation involves inverting a $T\times T$ matrix,
where $T$ is the number of observations, which results in $\mathcal{O}(T^{3})$
complexity.

The kernel determines almost all the generalization properties of a GP model,
and that is why many times the mean function is chosen to be a constant, i.e.
$m(\boldsymbol{x})=\boldsymbol{c}$.

We now expand on different kernel functions and how to choose them.

#####  2.1.3.1 Kernels

A kernel function that is used as a covariance function can be chosen upon the
specific problem, but it must be positive semi-definite and symmetric to be
valid. I will elaborate on some popular kernels, how to combine kernels in
general and on their parameters optimization, some of it is taken from [5].

######  2.1.3.1.1 Standard Kernels

Some popular kernels and possible derivations are:

###### Radial Basis Function

Radial basis function (RBF) is the default or “go-to” kernel in most
implementations, and can be defined as

$K_{RBF}(x,x^{\prime})=\sigma^{2}\exp\left({-\frac{(x-x^{\prime})^{2}}{2\ell^{2}}}\right)$

where $\ell$, the length scale, controls the length of the “wiggles” of the
function, and $\sigma^{2}$ is a scale factor of the distance from the mean.

###### Periodic

A periodic kernel is good for functions that repeat themselves exactly. A
periodic kernel can be defined as:

$K_{Per}(x,x^{\prime})=\sigma^{2}\exp\left({-\frac{2\sin{(\pi|x-x^{\prime}|p)}^{2}}{\ell^{2}}}\right)$

where $p$ determines the periodicity and $\ell,\sigma^{2}$ are the same as in
RBF kernels.

###### Linear

The linear kernel is non stationary as opposed to RBF and periodic kernels,
i.e., its values are not dependant of the distance between the inputs but on
the inputs themselves. A linear kernel can be defined as:

$K_{Lin}(x,x^{\prime})=\sigma^{2}xx^{\prime}$

here, again, $\sigma$ is the output scale.

######  2.1.3.1.2 Combining Kernels

Since positive definiteness is closed under addition and multiplication, to
create new kernels, one can combine kernels through addition or
multiplication. Generally speaking, one can think of addition of kernels as an
‘or’ operator and multiplying kernels as an ‘and’ operator. For example to
create a periodic function with increasing variance in time, one can multiply
linear and periodic kernels to create a new kernel, i.e.,
$K(x,x^{\prime})=K_{Per}(x,x^{\prime})K_{Lin}(x,x^{\prime})$ as defined above.

######  2.1.3.1.3 Optimization of Hyperparameters

Choosing the right hyperparameters, like $\ell$ and $\sigma$ for the RBF
kernel, may be difficult, and in reality practitioners will either optimize
those hyperparameters over the set of observations $\mathcal{D}$ with, e.g.,
maximum likelihood estimation (MLE), or conduct Bayesian inference.

For example MLE in our regression problem would be:

$\boldsymbol{\theta}^{MLE}=\underset{\boldsymbol{\theta}}{\operatorname{\operatorname*{arg\,max}}}\;\mathcal{N}(\boldsymbol{y}|X,\boldsymbol{\theta})$

where $\boldsymbol{\theta}$ are the hyperparameters of the kernel function
chosen, e.g., if the kernel is an RBF kernel and $m(X)=0$, then
$\boldsymbol{\theta}=\\{\sigma,\ell\\}$.

###  2.2 Bayesian Graphical Models

####  2.2.1 Bayesian Machine Learning

Let us recall the Bayes’ rule:

$P(A|B)=\frac{P(B|A)P(A)}{P(B)}$ (2.1)

In Bayesian modeling we usually use the following notations:

  * •

$y$ \- the data or observations.

  * •

$\theta$ \- the latent (unobserved) parameters.

  * •

$\tilde{y}$ \- the predictions of the model.

So, when replacing $(A,B)$ with $(\theta,y)$ in (2.1) respectively, we get:

$P(\theta|y)=\frac{P(y|\theta)P(\theta)}{P(y)}$ (2.2)

These probabilities are referred to as:

$\text{posterior}=\frac{\text{likelihood}\times\text{prior}}{\text{evidence}}$

Note also that the definition of conditional probabilities allows us to write
(2.2) as:

$P(\theta|y)=\frac{P(y|\theta)P(\theta)}{P(y)}=\frac{P(\theta,y)}{P(y)}\propto
P(\theta,y)$ (2.3)

note that $P(y)$ is constant with respect to $\theta$ and only affects the
posterior as a normalizing constant.

Bayesian modeling for machine learning is done in three main steps:

  1. 1.

Specifying the generative model, which is the joint probability of all
observable and unobservable quantities in a problem, i.e., $p(\theta,y)$. Can
be thought of as the story of how the data came to be.

  2. 2.

Conditioning on observed data for two main missions: 1) to infer the
parameters or getting the posterior distribution, i.e., $p(\theta|y)$, also
called Bayesian updating; 2) to get the predictive distribution
$p(\tilde{y}|y)$.

  3. 3.

Evaluating the fit of the model: How well does the model fit the data? Are the
conclusions reasonable? How sensitive are the results to the modeling
assumptions?

One can refer to [12, 13] to read more on Bayesian data analysis and machine
learning.

####  2.2.2 Probabilistic Graphical Models

A probabilistic graphical model (PGM) is a model that defines a set of random
variables in a graph and their joint probability distribution. There are two
types of PGMs: 1) Directed Acyclic Graphs (DAGs) or Bayesian Networks (BNs);
2) Markov Random Fields (MRFs) which are undirected graphs. We will expand on
BNs since they are what we will use in this thesis. More on MRFs (and BNs) can
be found in [14].

#####  2.2.2.1 Bayesian Networks

A Bayesian network is a DAG whose nodes represent some random variables, e.g.,
$(X_{1},...,X_{N})$. The random variables represented in the BN can either be
discrete or continuous. For each node $X_{i}$ we define the conditional
probability distribution as $P(X_{i}|Pa(X_{i}))$, where $Pa(X_{i})$ denotes
the direct parents of $X_{i}$ (and $P(X_{i})$ if $X_{i}$ has no parents). With
that we can express the joint distribution of the variables by the chain rule
of BNs:

$P(X_{1},...,X_{N})=\prod_{i=1}^{N}P(X_{i}|Pa(X_{i}))$

i.e., the assumption is that each $X_{i}$ is independent of all other
nondescendents given its parents.

In the discrete case we would represent the conditional probability
distribution of a specific node in the form of a conditional probability
table. For example, if we have two Bernoulli random variables $(X_{1},X_{2})$,
with the following graph:

$X_{1}$$X_{2}$

we would need to build these conditional probability tables:

$X_{1}=0$ | $X_{1}=1$
---|---
$p$ | $1-p$
| $X_{2}=0$ | $X_{2}=1$
---|---|---
$X_{1}=0$ | $q_{0}$ | $1-q_{0}$
$X_{1}=1$ | $q_{1}$ | $1-q_{1}$

Where $(p,q_{0},q_{1})$ are some probabilities.

In the continuous case we would represent the conditional probability
distribution of a specific node in the form of a conditional probability
density function. For example, if we have two continuous random variables
$(\theta,y)$, with the following graph:

$\theta$$y$

we would need to define $p(\theta)$ and $p(y|\theta)$, and this defines the
generative model of the problem, i.e.:

$p(\theta,y)=p(\theta)p(y|\theta)$

####  2.2.3 Soft Evidence

Evidence is a term used for the occurrence of an event, e.g., $X=x$, i.e., the
event where random variable $X$ takes the value $x$. This is called hard
evidence. Sometimes we are not completely certain that the event has occurred,
but we believe that the event has happened with some probability, i.e., we
have moved from the prior belief of $p(X)$ to a new belief $p^{\prime}(X)$.
This is called soft evidence and is a form of uncertain evidence.

The idea of how to incorporate this new knowledge is discussed in [15] and
explained in relation to graphical models in [16].

We will use the example in [15, 16] to explain this: consider an agent that
inspects the color $C$ of a piece of cloth by candle light. That cloth has
three possible colors we will denote as $C_{1},C_{2},C_{3}$. Before the
inspection the agent had prior probability $p(C)$ of the discrete distribution
of the color, and after the inspection, he changed his beliefs to
$p^{\prime}(C)$. Of course he could not be sure of the color he saw since it
was seen with the light of the candle.

Notice that a usual relation in a graphical model (Section 2.2.2) will be that
the color $C$ is a latent parameter that is being observed in the evidence
$e$, but with soft evidence we have the other way around where the evidence
$e$ changed our belief in the latent parameter $C$, i.e., we have distribution
$C|e$, rather than a $e|C$, in our generative model.

Back to our example, let us imagine that we have a new parameter $S$ that
denotes the chance that this cloth will be sold and that $S$ depends only on
the cloth’s color, i.e., $S|C$. A picture of this network is presented in
Figure 2.1. We would like to update our belief in $p(S)$ to $p^{\prime}(S)$,
i.e., the new belief after the soft evidence $e$ — that is done with Jeffrey’s
Rule of updating:

$p^{\prime}(S)=\sum_{i}p(S|C_{i})p^{\prime}(C_{i})$

In that way we can incorporate the knowledge acquired by the evidence on our
parameter $S$ which is not directly connected to the evidence.

$C$$e$$S$ Figure 2.1: A network representing the graphical model described in
the example in 2.2.3. Undirected snake edge between $e$ and $C$ denotes soft
evidence about $C$ upon observing $e$.

Another way to incorporate uncertain evidence in a graphical model is with
virtual evidence. The basic idea is to add a virtual node $V$ in the graph
that always takes a certain value (with probability 1), and is connected to
the node of the parameter it convey evidence to (in our example it is node
$C$). Then the likelihood $p(V=1|e)$ needs to be defined and added that to the
joint distribution of the graph. This is not the method we use in this work,
so we will not add more on this subject, but one can refer to [16] for more
information.

###  2.3 Other Related Work

####  2.3.1 Forecast Reconciliation

Reconciliation is an active field of research, and in other related works it
is done in order to ensure coherency of the different forecasts of the
hierarchy, i.e. forecast for each level of the hierarchy should be equal to
the aggregation of its descendants.

As stated in Chapter 1, we would take a different approach of reconciling the
individual forecasts by sharing information between the series in order to
improve each individual forecast, rather then assuring coherency. This
approach is well described in Chapter 3. Nonetheless, we will note main works
in this field. Works can be divided to two main approaches: Linear
Reconciliation and Probabilistic Reconciliation.

#####  2.3.1.1 Linear Reconciliation

[17] introduced a least squares method, that was later improved by MinT [10],
which assumes the base forecasts are unbiased, while trying to produce the
minimum variance unbiased revised forecasts. A different approach, ERM [18],
does not assume unbiased forecasts, and tries to minimize the mean squared
revised forecast errors. Other linear reconciliation works include [19, 20].

#####  2.3.1.2 Probabilistic Reconciliation

[21] propose using copulas and solving a LASSO problem to get a hierarchical
coherent probabilistic forecasts. [22] propose a two level Bayesian approach,
firstly learning the distribution of bottom level forecasts through estimation
of the covariance matrix by historical errors, and then obtain aggregate
consistent point forecasts, by defining a problem-specific loss function. [23]
propose a fully Bayesian approach of reconciling the forecasts in order to
make them coherent, and we will evaluate our method against theirs in Chapter
5.

##  Chapter 3 Problem Statement

Let us define forecast reconciliation formally. In this work, we limit the
discussion to reconciliation of scalar real-valued forecasts. Reconciliation
of multiple multidimensional forecasts follows the same lines but involves a
much more elaborate theoretical analysis, with limited contribution to
understanding of the principles of Bayesian reconciliation.

######  Definition 1 (Forecast Reconciliation).

Consider a system consisting of $N$ subsystems. Assume that the parameter of
interest of the $i$-th subsystem is $x_{i}\in\mathbb{R}$, so that the vector
of parameters of interest of the whole system is
$\boldsymbol{x}\in\mathbb{R}^{N}$. Forecasts at a new location or a certain
time in the future are obtained

  * •

for each $x_{i}$

  * •

for $u=f(\boldsymbol{x})$ for some $f:\mathbb{R}^{N}\to\mathbb{R}^{M}$, $M\leq
N$

in the form of distributions $\Theta_{i}$ with density $\theta_{i}$ for each
$x_{i}$ and $\mathrm{H}$ with density $\eta$ for $u$. Then, forecast
reconciliation consists of inferring $\boldsymbol{x}|\Theta,\mathrm{H}$.

Function $f$ can be interpreted as computing summary statistics of some form
on $\boldsymbol{x}$. A popular form of $f$ is $A\boldsymbol{x}$ where $A$ is a
matrix of shape $M\times N$ with all elements equal to either 1 or 0, such
that each element of $u$ is the sum of a subset of components of
$\boldsymbol{x}$. In this case, $A$ apparently reflects the grouping of
subsystems such that all subsystems corresponding to entries of 1 in a row of
$A$ belong to the same group. Groups may overlap. However, other forms of $f$,
both linear and general, are conceivable; for example, $f$ may compute the
mean, the product, or the sum of squares of $\boldsymbol{x}$. It is worth
noting that $x_{i}$ are forecast independently, but $u$ must be forecast as a
whole; and that is why we set $\boldsymbol{x}$, but not $u$, in bold.
Intuitively, $u$ is a small set of summary statistics of $\boldsymbol{x}$ that
conveniently admits simultaneous forecasting.

$x_{1}$$x_{2}$$x_{3}$$x_{4}$$u_{1}$$u_{2}$ Figure 3.1: Example of a system
with four subsystems belong to two groups.

Consider the following example of the system appearing in Figure 3.1. In this
system we have $N=4$ subsystems in $M=2$ groups. In this example $f$ can be
any function such that $f:\mathbb{R}^{4}\to\mathbb{R}^{2}$. Continuing the
example, we define $f(\boldsymbol{x})=A\boldsymbol{x}$, and if the summary
statistics represent a summation, then:

$\boldsymbol{x}=\begin{bmatrix}x_{1}&x_{2}&x_{3}&x_{4}\end{bmatrix}^{T}$

and

$A=\begin{bmatrix}1&1&0&0\\\ 0&0&1&1\end{bmatrix}$

As stated in Definition 1, we are also given forecasts for each subsystem and
for $u$, and in our example those forecasts could be normally distributed:

  * •

$\boldsymbol{x}_{i}\sim\Theta_{i}\equiv\mathcal{N}(\mu_{i},\sigma_{i}^{2})\quad\forall
i\in\\{1,2,3,4\\}$

  * •

$u\sim\mathrm{H}\equiv\mathcal{MVN}(\boldsymbol{\mu}_{u},\Sigma_{u})$

and we would like to infer $\boldsymbol{x}|\Theta,\mathrm{H}$.

Notice that although $M=2$, we have a single (multivariate normally
distributed) forecast for $u$. Connecting this example to the introduction in
Chapter 1, we are assuming that although it is intractable to forecast
$\boldsymbol{x}$ in a single model when $N$ is very large (of course in this
toy example it is not the case), it is tractable to forecast $M$ when $M<<N$
in a single model, hence we have a single forecast distribution for $u$.

##  Chapter 4 Theoretical Analysis

###  4.1 Generative Model

$s$$y$$x$

Figure 4.1: A forecasting model. The latent state $s$ is observed through
observations $y$. A forecasted parameter of interest $x$ depends on $s$.

We approach the problem of forecast reconciliation using the methodology of
Bayesian generative modelling (Section 2.2). First, we describe the full
model. Then, we reduce the full network to a subnetwork relevant for
reconciliation of forecasts of $x_{i}$. Before turning to forecast
reconciliation, let us describe in general terms the structure of a Bayesian
network for forecasting.

######  Definition 2.

A Bayesian forecasting model is a generative model (Figure 4.1):

$\displaystyle s$ $\displaystyle\sim D_{s}$ (4.1) $\displaystyle y|s$
$\displaystyle\sim D_{y}(s)$ $\displaystyle x|s$ $\displaystyle\sim D_{x}(s)$

where

  * •

$s$ is a latent variable representing the state of the system;

  * •

$y$ is the observation through which $s$ is observed, i.e., historical values;

  * •

$x$ is the parameter of interest to be forecasted, i.e., future states;

  * •

$D_{s}$, $D_{y}(s)$, and $D_{x}(s)$ are distributions.

A forecast is the marginal distribution of $x$ given $y$.

There are $N+1$ forecasting models involved in the reconciliation problem
(Definition 1):

  * •

one model for forecasting $u$;

  * •

$N$ models for forecasting each of $x_{i}$, $i\in 1...N$.

We assume that all forecasting models have access to the same observation $y$.
Instead of the prior and conditional distributions for each of the models, we
are given the forecasts — the posterior marginal distributions $\Theta_{i}$ of
$x_{i}$ and $\mathrm{H}$ of $u$.

Based on the forecast of $u$, we aim to update the forecasts of $x_{i}$. To
formalize this, we imagine yet another forecasting model, for the whole
system, for which neither the prior or the conditional distributions, nor the
marginal posterior of $\boldsymbol{x}$ given $y$ are given, but which is
probabilistically related to the forecasting models for $x_{i}$ and $u$. The
full model, combining the individual forecast models for $x_{i}$ and $u$ as
well as the imagined forecasting model for $\boldsymbol{x}$, is shown in
Figure 4.2(a).

$x$$u$$\hat{x}_{i}$$\hat{u}$$s$$s_{i}$$s_{u}$$y$$N$

(a) full network

$x$$u$$\hat{x}_{i}$$\hat{u}$$s$$s_{i}$$s_{u}$$y$$N$

(b) relevant subnetwork

Figure 4.2: Bayesian network for forecast reconciliation. The undirected snake
edge between $u$ and $y$ denotes soft evidence about $u$ upon observing $y$.

In the full model,

  * •

$s$ and $\boldsymbol{x}$ are the latent state and the forecasted parameter of
interest of the imagined forecasting model for the whole system;

  * •

$s_{i}$ and $\hat{x}_{i}$, $i\in 1...N$, are the latent state and the
forecasted parameter of interest of the $i$-th subsystem;

  * •

$s_{u}$ and $\hat{u}$ are the latent state and the forecasted parameter of the
summary $u=f(\boldsymbol{x})$.

Observation $y$ is shared by all forecasting models. Dashed edges denote
probabilistic relationships between $\boldsymbol{x}$ and $\hat{x}_{i}$ and
$u=f(\boldsymbol{x})$ and $\hat{u}$. However, Definition 1 does not specify
these relationships. We assume that

  * •

the forecast $\hat{u}$ is accurate, that is $p(u|y)\equiv p(\hat{u}|y)$;

  * •

although $\boldsymbol{x}|y$ cannot be constructed from $\hat{x}_{i}|y$ because
interdependencies between $x_{i}$ are unknown, the product distribution of
$\hat{x}_{i}$ is suitable as a prior for $\boldsymbol{x}$.

Based on these assumptions, we reduce the full network in Figure 4.2(a) to a
subnetwork in Figure 4.2(b) required for reconciliation, dropping parts of the
full network with unknown prior or conditional probabilities. In particular,
since neither the prior distributions of state variables $s$, $s_{i}$,
$s_{u}$, nor the conditional distributions of observations or forecasts given
the state variables are known, we omit them in the subnetwork in Figure
4.2(b). The snake-shaped undirected link between $y$ and $u$ denotes that $y$
provides soft evidence 2.2.3 about $u$, that is, the posterior distribution
$u|y$, rather than the conditional distribution $y|u$, is given.

######  Definition 3.

A Bayesian forecast reconciliation model is a Bayesian generative model
(Figure 4.2(b)):

$\displaystyle\boldsymbol{x}$
$\displaystyle\sim\Theta_{1}\times\Theta_{2}\times...\times\Theta_{N}$ (4.2)
$\displaystyle u|\boldsymbol{x}$ $\displaystyle\sim\delta(f(\boldsymbol{x}))$
$\displaystyle u|y$ $\displaystyle\sim\mathrm{H}$

That is,

  * •

$\boldsymbol{x}$ is the latent variable to be inferred, with the prior
constructed as the product of marginal forecasts of its components;

  * •

$u$ is a nuisance latent variable with Dirac delta density
$p(u|\boldsymbol{x})=\delta(u-f(\boldsymbol{x}))$;

  * •

observation $y$ provides soft evidence about $u$.

The objective of inference on (4.2) is the posterior distribution
$\boldsymbol{x}|y$. The computation of $p(\boldsymbol{x}|y)$ is formalized by
the following proposition:

######  Proposition 4.1.

Under the assumptions of model (4.2),

$p(\boldsymbol{x}|y)=p(\boldsymbol{x})\frac{p(u=f(\boldsymbol{x})|y)}{p(u=f(\boldsymbol{x}))},$
(4.3)

where

$\displaystyle p(\boldsymbol{x})$
$\displaystyle\triangleq\prod_{i=1}^{N}\theta_{i}(x_{i})$ (4.4) $\displaystyle
p(u|y)$ $\displaystyle\triangleq\eta(u)$ (4.5) $\displaystyle p(u)$
$\displaystyle\triangleq\int_{X}\delta(u-f(\boldsymbol{x}))p(\boldsymbol{x})d\boldsymbol{x}$
(4.6)

######  Proof.

$p(\boldsymbol{x}|y)=\int_{U}p(\boldsymbol{x}|u)p(u|y)du$ (4.7)

$p(u|y)=\eta(u)$ is known, thus we only need to express $p(\boldsymbol{x}|u)$
in known terms. By the Bayes rule,

$p(\boldsymbol{x}|u)=p(\boldsymbol{x})\frac{p(u|\boldsymbol{x})}{p(u)}=p(\boldsymbol{x})\frac{\delta(u-f(\boldsymbol{x}))}{p(u)}$
(4.8)

Substituting (4.8) into (4.7), we obtain

$\displaystyle p(\boldsymbol{x}|y)$
$\displaystyle=\int_{U}p(\boldsymbol{x})\frac{\delta(u-f(\boldsymbol{x}))}{p(u)}p(u|y)du$
(4.9)
$\displaystyle=p(\boldsymbol{x})\frac{p(u=f(\boldsymbol{x})|y)}{p(u=f(\boldsymbol{x}))}$

$\square$

####  4.1.1 Linear Gaussian Reconciliation

In applications, such as those mentioned in Chapter 3, reconciliation is often
performed with $f(\boldsymbol{x})=A\boldsymbol{x}$ and Gaussian forecast
distributions. In this case, $p(\boldsymbol{x}|y)$ can be obtained in the
closed form of the density of a multivariate Gaussian distribution, as
formalized by the following proposition:

######  Proposition 4.2.

Let $f(\boldsymbol{x})=A\boldsymbol{x}$. Let the forecast distributions
$\Theta$ and $\mathrm{H}$ be the Gaussian distributions, parameterized by
$(\boldsymbol{\mu}_{\theta},\Sigma_{\theta})$,
$(\boldsymbol{\mu}_{\eta},\Sigma_{\eta})$, i.e.,

  * •

$p(\boldsymbol{x})=p_{\mathcal{N}}(\boldsymbol{x}|\boldsymbol{\mu}_{\theta},\Sigma_{\theta})$,

  * •

$p(f(\boldsymbol{x})|y)=p_{\mathcal{N}}(f(\boldsymbol{x})|\boldsymbol{\mu}_{\eta},\Sigma_{\eta})$.

Then,

$p(\boldsymbol{x}|y)\propto
p_{\mathcal{N}}(\boldsymbol{x}|\boldsymbol{\mu}_{LG},\Sigma_{LG})$ (4.10)

where

$\displaystyle\Sigma_{LG}$
$\displaystyle=(\Sigma_{\theta}^{-1}+A^{T}[\Sigma_{\eta}^{-1}-(A\Sigma_{\theta}A^{T})^{-1}]A)^{-1}$
(4.11) $\displaystyle\boldsymbol{\mu}_{LG}$
$\displaystyle=\boldsymbol{\mu}_{\theta}+\Sigma_{LG}A^{T}\Sigma_{\eta}^{-1}(\boldsymbol{\mu}_{\eta}-A\boldsymbol{\mu}_{\theta})$
(4.12)

######  Proof.

$f(\boldsymbol{x})=A\boldsymbol{x}$ is an affine transformation on the
multivariate Gaussian distributed parameter $\boldsymbol{x}$, which implies:

$p(f(\boldsymbol{x})|\boldsymbol{\theta})=p_{\mathcal{N}}(f(\boldsymbol{x})|A\boldsymbol{\mu}_{\theta},A\Sigma_{\theta}A^{T})$
(4.13)

Using (4.3) with (4.13) and the distributions in Proposition 4.2:

$\displaystyle p(\boldsymbol{x}|y)$
$\displaystyle=p_{\mathcal{N}}(\boldsymbol{x}|\boldsymbol{\mu}_{\theta},\Sigma_{\theta})\frac{p_{\mathcal{N}}(f(\boldsymbol{x})|\boldsymbol{\mu}_{\eta},\Sigma_{\eta})}{p_{\mathcal{N}}(f(\boldsymbol{x})|A\boldsymbol{\mu}_{\theta},A\Sigma_{\theta}A^{T})}$
(4.14)
$\displaystyle\propto\exp\left[(\boldsymbol{x}-\boldsymbol{\mu}_{\theta})^{T}\Sigma_{\theta}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{\theta})\right]\frac{\exp\left[(A\boldsymbol{x}-\boldsymbol{\mu}_{\eta})^{T}\Sigma_{\eta}^{-1}(A\boldsymbol{x}-\boldsymbol{\mu}_{\eta})\right]}{\exp\left[(A\boldsymbol{x}-A\boldsymbol{\mu}_{\theta})^{T}(A\Sigma_{\theta}A^{T})^{-1}(A\boldsymbol{x}-A\boldsymbol{\mu}_{\theta})\right]}$
(4.15) $\displaystyle=\exp(g(\boldsymbol{x}))$ (4.16)

Here $g(\boldsymbol{x})$ represents the grouping of all the terms in the
exponents of (4.15).
Expanding the terms in the exponent of (4.15):

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}\Sigma_{\theta}^{-1}\boldsymbol{x}+\boldsymbol{x}^{T}A^{T}\Sigma_{\eta}^{-1}A\boldsymbol{x}-\boldsymbol{x}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{x}$
(4.17)
$\displaystyle-2\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}\boldsymbol{x}-2\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A\boldsymbol{x}+2\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{x}$
$\displaystyle+\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}\boldsymbol{\mu}_{\theta}+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}\boldsymbol{\mu}_{\eta}-\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{\mu}_{\theta}$

Let us denote the terms in the third line of 4.17 by $C$ defined as follows:

$C=\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}\boldsymbol{\mu}_{\theta}+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}\boldsymbol{\mu}_{\eta}-\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{\mu}_{\theta}$

Here $C$ does not depend on $\boldsymbol{x}$ and only affects the
normalization constant of the density of Equation (4.10). Replacing $C$ in
4.17:

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}\Sigma_{\theta}^{-1}\boldsymbol{x}+\boldsymbol{x}^{T}A^{T}\Sigma_{\eta}^{-1}A\boldsymbol{x}-\boldsymbol{x}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{x}$
(4.18)
$\displaystyle-2\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}\boldsymbol{x}-2\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A\boldsymbol{x}+2\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A\boldsymbol{x}$
$\displaystyle+C$

Grouping the terms in (4.18):

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}(\Sigma_{\theta}^{-1}+A^{T}\Sigma_{\eta}^{-1}A-A^{T}(A\Sigma_{\theta}A^{T})^{-1}A)\boldsymbol{x}$
(4.19)
$\displaystyle-2(\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A-\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A)\boldsymbol{x}$
$\displaystyle+C$

Adding and subtracting the term
$(\boldsymbol{\mu}_{\theta}^{T}A^{T}\Sigma_{\eta}^{-1}A)$ from the parentheses
in the second line of (4.19):

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}(\Sigma_{\theta}^{-1}+A^{T}\Sigma_{\eta}^{-1}A-A^{T}(A\Sigma_{\theta}A^{T})^{-1}A)\boldsymbol{x}$
(4.20)
$\displaystyle-2(\boldsymbol{\mu}_{\theta}^{T}\Sigma_{\theta}^{-1}+\boldsymbol{\mu}_{\theta}^{T}A^{T}\Sigma_{\eta}^{-1}A-\boldsymbol{\mu}_{\theta}^{T}A^{T}(A\Sigma_{\theta}A^{T})^{-1}A+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A-\boldsymbol{\mu}_{\theta}^{T}A^{T}\Sigma_{\eta}^{-1}A)\boldsymbol{x}$
$\displaystyle+C$

Grouping terms:

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}(\Sigma_{\theta}^{-1}+A^{T}\Sigma_{\eta}^{-1}A-A^{T}(A\Sigma_{\theta}A^{T})^{-1}A)\boldsymbol{x}$
(4.21)
$\displaystyle-2(\boldsymbol{\mu}_{\theta}^{T}[\Sigma_{\theta}^{-1}+A^{T}\Sigma_{\eta}^{-1}A-A^{T}(A\Sigma_{\theta}A^{T})^{-1}A]+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A-\boldsymbol{\mu}_{\theta}^{T}A^{T}\Sigma_{\eta}^{-1}A)\boldsymbol{x}$
$\displaystyle+C$

Introducing the term $\Sigma_{LG}$ as defined in Equation (4.11) and replacing
in (4.21):

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}$ (4.22)
$\displaystyle-2(\boldsymbol{\mu}_{\theta}^{T}\Sigma_{LG}^{-1}+\boldsymbol{\mu}_{\eta}^{T}\Sigma_{\eta}^{-1}A\Sigma_{LG}\Sigma_{LG}^{-1}-\boldsymbol{\mu}_{\theta}^{T}A^{T}\Sigma_{\eta}^{-1}A\Sigma_{LG}\Sigma_{LG}^{-1})\boldsymbol{x}$
$\displaystyle+C$

Grouping terms:

$\displaystyle g(\boldsymbol{x})=$
$\displaystyle\boldsymbol{x}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}$ (4.23)
$\displaystyle-2(\boldsymbol{\mu}_{\theta}+\Sigma_{LG}A^{T}\Sigma_{\eta}^{-1}\boldsymbol{\mu}_{\eta}-\Sigma_{LG}A^{T}\Sigma_{\eta}^{-1}A\boldsymbol{\mu}_{\theta})^{T}\Sigma_{LG}^{-1}\boldsymbol{x}$
$\displaystyle+C$

Introducing the term $\boldsymbol{\mu}_{LG}$ as defined in Equation (4.12) and
substituting it in (4.23):

$g(\boldsymbol{x})=\boldsymbol{x}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}-2\boldsymbol{\mu}_{LG}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}+C$
(4.24)

Plugging Equation (4.24) in Equation (4.15):

$\displaystyle p(\boldsymbol{x}|y)$
$\displaystyle\propto\exp\Big{(}\boldsymbol{x}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}-2\boldsymbol{\mu}_{LG}^{T}\Sigma_{LG}^{-1}\boldsymbol{x}+C\Big{)}$
(4.25)
$\displaystyle\propto\exp\Big{(}(\boldsymbol{x}-\boldsymbol{\mu}_{LG})^{T}\Sigma_{LG}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{LG})\Big{)}$
(4.26) $\displaystyle\propto
p_{\mathcal{N}}(\boldsymbol{x}|\boldsymbol{\mu}_{LG},\Sigma_{LG})$ (4.27)

Which completes the proof. $\square$

##  Chapter 5 Empirical Evaluation

In this chapter we evaluate the Linear Gaussian model described in Section
4.1.1.

In each evaluation (experiment) we compare two reconciliation methods — the
first method is the one presented in this thesis, Section 4.1.1, which we
denote by ET Reconciler, the second is the LG method from [23] which we denote
by Bayes Paper LG Reconciler. For each reconciliation method tested we produce
a reconciled probability $\tilde{\mathcal{D}}_{t}$ for each time step $t$.

We perform those experiments on two datasets — synthetic and real data set,
where these evaluations are described in Section 5.2 and 5.3 respectively.
Generally, we fit a GP to each series to obtain the base forecast, then we
apply the reconciliation method to try to improve these forecasts, and compute
the metric described in Section 5.1 before and after the reconciliation. On
the synthetic data set we performed two evaluations (or experiments) that
differ in the way the base forecasts are generated — this will be described in
Section 5.2.

###  5.1 Metrics

For each reconciliation method we compute the Negative Log Predictive Density
(NLPD):

$\text{NLPD}=\frac{1}{T}\sum_{t}-\log{p_{\tilde{\mathcal{D}}_{t}}(\boldsymbol{x}_{t})}$
(5.1)

i.e., the average, over the time steps, of the negative log of the true values
$\boldsymbol{x}$ under the reconciled probability density function
$\tilde{\mathcal{D}}$. For each data set we perform $N=50$ simulations and
compute the NLPD. We report the mean, the standard error of the mean (denoted
as SEM) and the standard deviation (denoted as SD) of the NLPDs of these $N$
simulations for each reconciliation method (and for the base forecasts as a
benchmark).

###  5.2 Synthetic Data Set

####  5.2.1 Data — Experiment A

$u$$x^{(2)}$$x^{(1)}$$x^{(3)}$$x^{(4)}$

Figure 5.1: The synthetic data set for Experiment 5.2.1. The four systems
$\\{x^{(n)}\\}_{n=1}^{4}$ are connected to the same parent $u$.

The first experiment includes a synthetic dataset consisting of a hierarchy of
4 series all connected to same parent, as shown in Figure 5.1. Each series is
a sine wave, corrupted by noise. The observations are the first period of the
sine wave and the task is to predict the second cycle.

Observations for each series are denoted by $y^{(i)}_{t}$, and the future
(hidden) states by $x^{(i)}_{t}$ for the bottom series and $u_{t}$ for the sum
of the series. A bold symbol, $\boldsymbol{x}^{(n)}\in\mathbb{R}^{T}$, is a
vector of all time steps, and $\boldsymbol{x}_{t}\in\mathbb{R}^{N}$, is a
vector of all series in time step $t$.

Data was created:

$\displaystyle y_{t}^{(n)}$ $\displaystyle=\sin(t)+\epsilon_{t}^{(n)},\quad
t\in[0,2\pi),\quad\epsilon_{t}^{(n)}\sim\mathcal{N}(0,\sigma_{\epsilon})\quad\forall
n=\\{1,2,3,4\\}$ (5.2) $\displaystyle x_{t}^{(n)}$
$\displaystyle=\sin(t)+\epsilon_{t}^{(n)},\quad
t\in[2\pi,4\pi],\quad\epsilon_{t}^{(n)}\sim\mathcal{N}(0,\sigma_{\epsilon})\quad\forall
n=\\{1,2,3,4\\}$ (5.3)

In the different experiments on this data set we experimented with the
following noise levels: $\sigma_{\epsilon}\in\\{0,0.1,0.2,0.5,1,2,5\\}$.

The sum of series is given by:

$u_{t}=f(\boldsymbol{x}_{t})=\sum_{i=1}^{4}x^{(i)}_{t}$ (5.4)

Base forecasts were created by fitting a GP (2.1.3) to the observations
$\boldsymbol{y}^{(n)}$ and are parameterized by $\\{\Theta^{(n)}\\}_{n=1}^{4}$
for each series $\\{\boldsymbol{x}^{(n)}\\}_{n=1}^{4}$ and $\mathrm{H}$ for
$\boldsymbol{u}$:

$\hat{\boldsymbol{x}}^{(n)}\equiv\boldsymbol{x}^{(n)}|\Theta^{(n)}\sim\mathcal{GP}^{(n)}\quad\forall
n=\\{1,2,3,4\\},\quad\boldsymbol{u}|\mathrm{H}\sim\mathcal{GP}$ (5.5)

The GP’s are set with a scaled periodic kernel (see 2.1.3.1.1). The period
length hyperparameter $p$ is assumed to be known and fixed at $2\pi$, which is
a common case in many applications where the seasonality is known.

Reconciliation for each time step is performed independently, so the
reconciler actually gets the marginalised prediction for time step $t$:

$\hat{x}_{t}^{(n)}\equiv
x_{t}^{(n)}|\Theta_{t}^{(n)}\sim\mathcal{N}_{t}^{(n)}\quad\forall
n=\\{1,2,3,4\\},\quad u_{t}|\mathrm{H}_{t}\sim\mathcal{N}_{t}$

A visualization of the data set and the base forecasts can be seen in Figure
5.2.

Figure 5.2: Example of a data set and its base forecasts for
$\sigma_{\epsilon}=\\{0.1,0.2,0.5\\}$. Black x’s denote training data
(observations $\boldsymbol{y}^{(n)}$), green x’s denote test data
($\boldsymbol{x}^{(n)}$), blue lines and shaded blue areas denote the mean and
95% confidence interval of the base forecasts, respectively.

####  5.2.2 Results - Experiment A

|  | Base Forecast | Bayes Paper LG Reconciler | ET Reconciler
---|---|---|---|---
$\sigma_{\epsilon}=0$ | mean | -15.73 | -15.85 | -15.04
SEM | 0.00 | 0.00 | 0.00
SD | 0.00 | 0.00 | 0.00
$\sigma_{\epsilon}=0.1$ | mean | 2.13 | 2.88 | 0.26
SEM | 0.52 | 0.59 | 0.35
SD | 3.67 | 4.18 | 2.48
$\sigma_{\epsilon}=0.2$ | mean | 3.77 | 4.37 | 2.43
SEM | 0.47 | 0.53 | 0.33
SD | 3.35 | 3.78 | 2.31
$\sigma_{\epsilon}=0.5$ | mean | 6.61 | 7.17 | 5.73
SEM | 0.37 | 0.42 | 0.27
SD | 2.59 | 2.95 | 1.91
$\sigma_{\epsilon}=1$ | mean | 44.38 | 48.55 | 39.92
SEM | 21.82 | 24.00 | 19.32
SD | 154.31 | 169.71 | 136.61
$\sigma_{\epsilon}=2$ | mean | 139.82 | 163.48 | 134.89
SEM | 109.48 | 127.55 | 106.00
SD | 774.13 | 901.95 | 749.51
$\sigma_{\epsilon}=5$ | mean | 29.84 | 34.14 | 29.74
SEM | 9.64 | 11.68 | 9.43
SD | 68.16 | 82.60 | 66.71
Table 5.1: Results for the experiment on Data Set 5.2.1.

In Table 5.1 are the results for Experiment A.

The first takeaway is that the algorithm proposed in 4.2 works, i.e., the ET
Reconciler improves the NLPD score compared to the base forecasts (except for
the case where $\sigma_{\epsilon}=0$ which we will explain in a bit).

Another takeaway from looking at the results is that we can notice that except
for when the data was generated with $\sigma_{\epsilon}=0$, ET Reconciler
always outperformed Bayes Paper LG Reconciler. Here is some intuition why.

The base forecasts were generated using GPs (see 5.5), and the kernels of a GP
have an additive property, i.e., the sum of the four periodic kernels of the
bottom time series is also a periodic kernel. So, we expect the mean of the
base forecasts for the upper time series to be close to the sum of the means
of the base forecasts for the bottom time series, and in the notations of
4.1.1, we expect $A\boldsymbol{\mu}_{\theta}$ and $\boldsymbol{\mu}_{\eta}$ to
be close. In both reconciliation methods (ET Reconciler and Bayes Paper LG
Reconciler) there is a mechanism (see $\mu_{LG}$ in Proposition 4.12 and
Section 3.1 in [23]) to correct the means of the base forecasts if there is a
difference between the mean of the upper series and the means of the bottom
series after being passed through the mapping $f(\boldsymbol{x})$, in this
data set it is a summation. Since we used GPs to generate the forecasts there
is not a big difference between the means and both reconciliation methods
report a similar mean.

So the difference between the reconcilers must come from the covariance.

What happens is that ET Reconciler adds to the covariance of the base
predictions the difference between the covariance of the prediction for the
upper series and the transformed (after the linear mapping) covariance of the
prediction for the bottom series (see $\Sigma_{LG}$ in Proposition 4.11),
wheres in Bayes Paper LG Reconciler they always reduce the covariance of the
reconciled predictions (see Section 3.1 in [23]). This property of the
reconcilers makes ET Reconciler enlarge its total variance, whereas Bayes
Paper LG Reconciler reduces its total variance, and as the data set is getting
more ”noisy”, i.e., $\sigma_{\epsilon}$ gets larger (up to a certain level)
the reconciler with the larger total variance gets a better NLPD score.

As for the experiment where the data set has no noise ($\sigma_{\epsilon}=0$),
what happens is that the base forecasts are quite accurate to begin with and
their mean is accurate, so in order to make the NLPD score better one needs to
reduce the variance, and that is why this is the only case that Bayes Paper LG
Reconciler gets a better score. This is of course a very extreme case where
the base forecasts are accurate and there is no noise in the data set. Also,
the magnitude of improvement of the NLPD is very small compared to other
settings of $\sigma_{\epsilon}$.

####  5.2.3 Data — Experiment B

This setting is the exact same setting as in 5.2.1, with the only difference
that the base forecasts themselves, for the bottom level only, were ”shifted”
for each series by some random noise. So instead of Equation 5.5 we now have:

$\hat{\boldsymbol{x}}^{(n)}_{B}=\hat{\boldsymbol{x}}^{(n)}_{A}+\zeta^{(n)},\quad\zeta^{(n)}\sim\mathcal{N}(0,1)\quad\forall
n=\\{1,2,3,4\\},\quad\boldsymbol{u}|\mathrm{H}\sim\mathcal{GP}$ (5.6)

where $\hat{\boldsymbol{x}}^{(n)}_{A},\hat{\boldsymbol{x}}^{(n)}_{B}$ are the
base forecasts in Experiments A and B respectively. Reconciliation is done the
same as in 5.2.1.

A visualization of the data set and the ”shifted” base forecasts can be seen
in Figure 5.3.

Figure 5.3: Example of a data set and its ”shifted” base forecasts for
$\sigma_{\epsilon}=\\{0.1,0.2,0.5\\}$. Black x’s denote training data
(observations $\boldsymbol{y}^{(n)}$), green x’s denote test data
($\boldsymbol{x}^{(n)}$), blue lines and shaded blue areas denote the mean and
95% confidence interval of the base forecasts, respectively.

####  5.2.4 Results - Experiment B

|  | Base Forecast | Bayes Paper LG Reconciler | ET Reconciler
---|---|---|---|---
$\sigma_{\epsilon}=0$ | mean | 15947.89 | 15191.35 | 12166.44
SEM | 1727.36 | 1675.84 | 1549.79
SD | 12214.29 | 11849.97 | 10958.69
$\sigma_{\epsilon}=0.1$ | mean | 367.63 | 348.50 | 265.19
SEM | 46.69 | 44.55 | 38.77
SD | 330.15 | 315.02 | 274.13
$\sigma_{\epsilon}=0.2$ | mean | 109.28 | 103.90 | 79.01
SEM | 10.64 | 9.83 | 7.75
SD | 75.21 | 69.51 | 54.83
$\sigma_{\epsilon}=0.5$ | mean | 20.82 | 20.59 | 16.75
SEM | 1.86 | 1.80 | 1.60
SD | 13.15 | 12.70 | 11.30
$\sigma_{\epsilon}=1$ | mean | 52.74 | 58.21 | 43.30
SEM | 23.05 | 27.23 | 18.19
SD | 163.00 | 192.53 | 128.64
$\sigma_{\epsilon}=2$ | mean | 164.91 | 177.92 | 153.76
SEM | 102.91 | 112.67 | 98.58
SD | 727.67 | 796.71 | 697.07
$\sigma_{\epsilon}=5$ | mean | 107.73 | 119.53 | 106.10
SEM | 41.73 | 47.53 | 42.00
SD | 295.09 | 336.06 | 296.97
Table 5.2: Results for the experiment on Data Set 5.2.3.

In Table 5.2 are the results for Experiment B.

Here we can notice that again for all settings of $\sigma_{\epsilon}$ ET
Reconciler outperforms Bayes Paper LG Reconciler, even when
$\sigma_{\epsilon}=0$. This happens here because for all settings the
reconcilers must change the mean of the base forecasts, as well as the
covariance, so the fact that for $\sigma_{\epsilon}=0$ there is no noise does
not affect the improvement of the NLPD compared to the base forecasts. As in
the previous experimental results (Table 5.1), the NLPD improves up to a
certain level of noise in the data. Also notice that here even for a very
large noise ($\sigma_{\epsilon}=5$), ET Reconciler improves the NLPD whereas
Bayes Paper LG Reconciler does not.

###  5.3 Real Data Set

####  5.3.1 Data

Figure 5.4: The first level of the hierarchy of the dataset. Top graph: the
number of visitors by month to each of the seven states of Australia. Bottom
graph: The sum of the series in the top graph.

The real data set we tested on is a data set consisting of the number of
tourists who entered Australia on a monthly basis, by their reason of entrance
(Holiday, Visitor, Business, Other), between 1998 to 2016. On the bottom level
the data set contains 304 time series referring to Australia’s 76 regions and
4 reasons of entrance. The data is then aggregated into 27 macro zones and
then into 7 states on the upper level. We chose to test the 7 series of the
total tourists entering each state. A visualization of the data set is shown
in Figure 5.4.

Base forecasts were created using the same methodology of the synthetic data
set (see 5.2.1), with the only difference that the period length
hyperparameter $p$ was set to $12$ and that we added an RBF kernel (see
2.1.3.1.1).

Reconciliation here is also performed independently as in 5.2.1, and got the
same marginalised prediction for time step $t$. The main change here regarding
train and test data, is that we are splitting the data set into training and
test at a different point in each simulation we perform, i.e., we perform 50
simulations, each with 60 time steps (5 years) as training data and 24 time
steps (2 years) as test data, starting from the first period.

####  5.3.2 Results

| mean | SEM | SD
---|---|---|---
Base Forecast | 60.56 | 0.41 | 2.88
Bayes Paper LG Reconciler | 61.19 | 0.41 | 2.87
ET Reconciler | 59.35 | 0.42 | 2.99
Table 5.3: Results for the experiment on Data Set of Section 5.3.1.

Table 5.3 depicts the results for the experiment on the real data set. Here we
can see that our proposed method works for the real data set as well, i.e.,
the NLPD of the reconciled forecasts decreases with respect to the base
forecast’s NLPD.

In this dataset, the observations for the bottom series are quite noisy, while
for the upper series they are less noisy (see Figure 5.4). The base forecasts
were generated with a GP with periodic and RBF kernels, which generated quite
good base forecasts in the mean sense, i.e., they were quite accurate in their
mean but not in their variance. Due to the additive property of the GPs, the
mean of the sum of the base forecasts for the bottom level was close to the
mean of the upper level series, but the base forecasts for the bottom level
were too “certain”, i.e., had relatively small variance, that had to be taken
care of. Bayes Paper LG Reconciler, as mentioned in Section 5.2.2, always
decreased the total variance of the reconciled forecasts, while our ET
Reconciler did not. So, ET Reconciler made the reconciled forecasts less
”certain”, while Bayes Paper LG Reconciler made the base forecasts more
”certain”, i.e., decreasing and increasing total variance respectively. Our
metric is a probabilistic metric, i.e., it takes into account the density of
the prediction, and that property along with the fact that the base forecasts
were too “certain” but quite accurate in the mean sense, made ET Reconciler
improve the NLPD and Bayes Paper LG Reconciler worsen the NLPD, compared to
that of the base forecasts.

##  Chapter 6 Conclusion

Time series forecasting is a key task in the ability to make informed
decisions under uncertainty. There are many cases where the number of series
we need to predict is too large to fit in a single model. In that case, a
distributed approach combined with a way to utilize the knowledge of the
hierarchical structure of these series, allows to get better predictions in a
reasonable runtime. The process is to learn base forecasts for each series
independently and for some summary statistics series based on the hierarchy
structure, and pass those base forecasts through a reconciliation algorithm to
improve each forecast. In this work we tackled the reconciliation problem,
with the assumption that the base forecasts are given to us.

Reconciliation is an active field of research, but it is worth noting that
works in this field, including the ones cited in Section 2.3.1, are trying to
solve a somewhat different problem — assuring coherency of the base forecasts,
i.e., that the summary statistics series forecasts are aligned with both the
forecasts for each series and the hierarchical structure. As stated before, we
tried to solve a somewhat different problem of reconciliation in a distributed
setting.

In this work we presented a method for Bayesian forecast reconciliation of
hierarchical time series. We defined the reconciliation problem as a
generative model, and derived the posterior distribution of the time series of
interest given the historical observations that is given to us through the
individual forecasts of the series and the forecast for the summary
statistics. We highlighted the assumptions in our case about the accuracy of
the base forecasts for the summary statistics, and about the prior for the
base forecasts for the series of interest. We finally addressed the specific
case of linear Gaussian reconciliation — a case where the base forecasts are
Gaussian distributed and the summary statistics have a linear mapping. We
showed, by deriving the formulas, that reconciliation for this case can be
achieved in closed form of a multivariate Gaussian distribution.

A limitation of the proposed method is that we assume that although it is
intractable to forecast the series in a single model where the number of
series is very large it is possible to forecast their summary statistics in a
single model. This is an assumption that we have to make for the validity of
the proposed method.

We evaluated our linear Gaussian reconciliation method on both a synthetic
data set and a real data set. We saw that our method improved the base
forecasts on the metric described in Section 5.1. We then compared our method
to the one proposed in [23] on both data sets by the same metric. The results
showed that our method outperformed theirs on both the synthetic data set and
the real data set. This difference might come from the fact that the data was
noisy and the adjustments that had to be made where in the variance of the
base forecasts and not on the mean of the forecasts, and the method proposed
in [23] does not do well in that case.

Our experiments were done on data sets that have a small number of series
while our proposed method was set to work on a setting with a very large
number of series, so, future work should definitely include testing the
methods proposed in this work on a large scale data set. Another direction is
reconciliation through time — we tested our method on reconciliation of each
time step independently, and future work could try to test this as well.

## References

  * [1] (economics), Nowcasting.  Nowcasting (economics) - wikipedia, the free encyclopedia, 2023.  URL https://en.wikipedia.org/wiki/Nowcasting(economics).
  * Box and Jenkins [1970] Box, George and Jenkins, Gwilym.  _Time series analysis: forecasting and control_.  San Francisco: Holden-Day, 1970.
  * Kalman [1960] Kalman, R. E.  A New Approach to Linear Filtering and Prediction Problems.  _Journal of Basic Engineering_ , 82(1):35–45, 03 1960.  ISSN 0021-9223.  doi: 10.1115/1.3662552.  URL https://doi.org/10.1115/1.3662552.
  * Rasmussen and Williams [2005] Rasmussen, Carl Edward and Williams, Christopher K. I.  _Gaussian Processes for Machine Learning_.  The MIT Press, 11 2005.  ISBN 9780262256834.  doi: 10.7551/mitpress/3206.001.0001.  URL https://doi.org/10.7551/mitpress/3206.001.0001.
  * Duvenaud [2014] Duvenaud, David.  _Automatic Model Construction with Gaussian Processes_.  PhD thesis, Computational and Biological Learning Laboratory, University of Cambridge, 2014.
  * Flunkert et al. [2017] Flunkert, Valentin, Salinas, David, and Gasthaus, Jan.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _CoRR_ , abs/1704.04110, 2017.  URL http://arxiv.org/abs/1704.04110.
  * Salinas et al. [2019] Salinas, David, Bohlke-Schneider, Michael, Callot, Laurent, Medico, Roberto, and Gasthaus, Jan.  High-dimensional multivariate forecasting with low-rank gaussian copula processes.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019\.  URL https://proceedings.neurips.cc/paper/2019/file/0b105cf1504c4e241fcc6d519ea962fb-Paper.pdf.
  * market index [2023] market index, Stock.  Stock market index - wikipedia, the free encyclopedia, 2023.  URL https://en.wikipedia.org/wiki/Stock-market-index.
  * Li et al. [2020] Li, Tian, Sahu, Anit Kumar, Talwalkar, Ameet, and Smith, Virginia.  Federated learning: Challenges, methods, and future directions.  _IEEE Signal Processing Magazine_ , 37(3):50–60, 2020.  doi: 10.1109/MSP.2020.2975749.
  * Wickramasuriya et al. [2019] Wickramasuriya, Shanika L., Athanasopoulos, George, and Hyndman, Rob J.  Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization.  _Journal of the American Statistical Association_ , 114(526):804–819, 2019.  doi: 10.1080/01621459.2018.1448825.  URL https://doi.org/10.1080/01621459.2018.1448825.
  * Bonilla et al. [2007] Bonilla, Edwin V, Chai, Kian, and Williams, Christopher.  Multi-task gaussian process prediction.  In Platt, J., Koller, D., Singer, Y., and Roweis, S., editors, _Advances in Neural Information Processing Systems_ , volume 20. Curran Associates, Inc., 2007.  URL https://proceedings.neurips.cc/paper/2007/file/66368270ffd51418ec58bd793f2d9b1b-Paper.pdf.
  * Gelman et al. [2013] Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., and D.B., Rubin.  _Bayesian Data Analysis (3rd ed.)_.  Chapman and Hall/CRC, 2013.  URL https://doi.org/10.1201/b16018.
  * McElreath [2020] McElreath, R.  _Statistical Rethinking: A Bayesian Course with Examples in R and Stan_.  A Chapman & Hall book. CRC Press, 2020.  ISBN 9780367139919.
  * Koller and Friedman [2009] Koller, Daphne and Friedman, Nir.  _Probabilistic Graphical Models: Principles and Techniques_.  MIT press, 2009.
  * Jeffrey [1990] Jeffrey, Richard C.  _The logic of decision_.  University of Chicago press, 1990.
  * Pearl [2009] Pearl, Judea.  _Probabilistic reasoning in intelligent systems : networks of plausible inference_.  Morgan Kaufmann, San Francisco, Calif., 2009.  ISBN 9781558604797 1558604790.
  * Hyndman et al. [2011] Hyndman, Rob J., Ahmed, Roman A., Athanasopoulos, George, and Shang, Han Lin.  Optimal combination forecasts for hierarchical time series.  _Computational Statistics & Data Analysis_, 55(9):2579–2589, 2011.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2011.03.006.  URL https://www.sciencedirect.com/science/article/pii/S0167947311000971.
  * Ben Taieb and Koo [2019] Ben Taieb, Souhaib and Koo, Bonsoo.  Regularized regression for hierarchical forecasting without unbiasedness conditions.  In _KDD 2019_ , pages 1337–1347, United States of America, 2019. Association for Computing Machinery (ACM).  doi: https://dl.acm.org/doi/proceedings/10.1145/3292500.  URL https://www.kdd.org/kdd2019/.  ACM International Conference on Knowledge Discovery and Data Mining 2019, KDD 2019 ; Conference date: 04-08-2019 Through 08-08-2019.
  * Hyndman et al. [2016] Hyndman, Rob J., Lee, Alan J., and Wang, Earo.  Fast computation of reconciled forecasts for hierarchical and grouped time series.  _Computational Statistics & Data Analysis_, 97:16–32, 2016.  ISSN 0167-9473.  doi: https://doi.org/10.1016/j.csda.2015.11.007.  URL https://www.sciencedirect.com/science/article/pii/S016794731500290X.
  * Panagiotelis et al. [2021] Panagiotelis, Anastasios, Athanasopoulos, George, Gamakumara, Puwasala, and Hyndman, Rob J.  Forecast reconciliation: A geometric view with new insights on bias correction.  _International Journal of Forecasting_ , 37(1):343–359, 2021.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2020.06.004.  URL https://www.sciencedirect.com/science/article/pii/S0169207020300911.
  * Taieb et al. [2017] Taieb, Souhaib Ben, Taylor, James W., and Hyndman, Rob J.  Coherent probabilistic forecasts for hierarchical time series.  In Precup, Doina and Teh, Yee Whye, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 3348–3357. PMLR, 2017\.  URL https://proceedings.mlr.press/v70/taieb17a.html.
  * Novak et al. [2017] Novak, Julie, McGarvie, Scott, and Garcia, Beatriz Etchegaray.  A bayesian model for forecasting hierarchically structured time series, 2017.
  * Corani et al. [2021] Corani, Giorgio, Azzimonti, Dario, Augusto, João P. S. C., and Zaffalon, Marco.  Probabilistic reconciliation of hierarchical forecast via bayes’ rule.  In Hutter, Frank, Kersting, Kristian, Lijffijt, Jefrey, and Valera, Isabel, editors, _Machine Learning and Knowledge Discovery in Databases_ , pages 211–226, Cham, 2021. Springer International Publishing.  ISBN 978-3-030-67664-3.