Hierarchically Coherent Multivariate Mixture Networks

Kin G. Olivares* David Luo* Cristian Challu Stefania La Vattiata

Max Mergenthaler Artur Dubrawski

Auton Lab Carnegie Mellon University Nixtla

###### Abstract

Large collections of time series data are often organized into hierarchies
with different levels of aggregation; examples include product and
geographical groupings. Probabilistic coherent forecasting is tasked to
produce forecasts consistent across levels of aggregation. In this study, we
propose to augment neural forecasting architectures with a coherent
multivariate mixture output. We optimize the networks with a composite
likelihood objective, allowing us to capture time series’ relationships while
maintaining high computational efficiency. Our approach demonstrates 13.2%
average accuracy improvements on most datasets compared to state-of-the-art
baselines. We conduct ablation studies of the framework components and provide
theoretical foundations for them. To assist related work, the code is
available at this http URL.

##  1 Introduction

Time series data is often organized into hierarchical structures spanning
different levels of aggregation. Independently forecasting any series from
such a hierarchy is unlikely to produce coherent results across the levels of
the hierarchy, that is, forecasts that satisfy the aggregation constraints of
the original data. In recent years, new hierarchical reconciliation has become
standard to achieve forecasts’ coherence (Babai et al., 2022), with
applications to electricity generation (Ben Taieb et al., 2021),
macroeconomics, and tourism (Eckert et al., 2021; Kourentzes and
Athanasopoulos, 2019).

(a) _Aggregation Constraints_

(b) _Error Accumulation_ (c) _Coherent Multivariate Mixture_

Figure 1: (a) Probabilistic coherent forecasting is a multivariate regression
problem with aggregation constraints. (b) Summing disaggregated level’s
forecasts can accumulate errors. (c) Augmenting neural forecasting
architectures with an accurate and coherent multivariate mixture is a
solution.

Previous neural hierarchical forecasting research is limited by restrictive
probabilistic assumptions, weak coherence enforcement, poor scalability, and
sometimes only modest accuracy improvements compared to statistical baselines.
To address these limitations, we introduce the _HIerarchically coherent
multivariate mixture NeTworks_ (HINT). Our contributions are:

  1. (i)

Coherent Multivariate Mixture for accurate and coherent forecasts, optimized
with composite likelihood. Our approach captures series’ relationships
maintaining high computational efficiency.

  2. (ii)

Probabilistic Reversible Instance Normalization to enhance forecast
distributions against stark variations in scales of the series. Our approach
normalizes inputs and reassembles probabilistic outputs via a global skip
connection.

  3. (iii)

State-of-the-art results on relevant benchmark data: Australian labour, SF Bay
Area traffic rates, Australian tourist visits, Wikipedia articles views.

This paper is structured as follows. Section 2 reviews literature and
introduces notation, Section 3 describes the methodology, Section 4 analyzes
our empirical results. Section 5 discusses results and future research
directions, while Section 6 concludes. Finally, in the Appendix, we study the
effects of the individual components of HINT and its theoretical foundations.

##  2 Related Work

Figure 2: A simple example of a four bottom-level series hierarchical
structure with geographic levels: the top level corresponds to the country,
the middle level to states, and the bottom level to regions.

###  2.1 Hierarchical Reconciliation

Classic hierarchical forecasting methods involve a two-stage process in which
a set of univariate statistical base forecasts are reconciled. There is a long
literature on these reconciliation strategies that include BottomUp (Orcutt et
al., 1968; Dunn et al., 1976), TopDown (Gross and Sohl, 1990; Fliedner, 1999),
and more recent optimal reconciliation strategies like Comb (Hyndman et al.,
2011), MinTrace (Wickramasuriya et al., 2019) and ERM (Ben Taieb and Koo,
2019). Probabilistic forecast reconciliation is at the forefront of
hierarchical forecasting research. Among the few methods capable of
probabilistic coherence, there is PERMBU that infuses multivariate
dependencies to bottom-level probabilities using copulas (Ben Taieb et al.,
2017), NORMALITY that reconciles a multivariate model under Gaussian
assumptions (Wickramasuriya, 2023), and BOOTSTRAP that generates reconciled
forecasts with bootstrap reconciliation (Panagiotelis et al., 2023).

###  2.2 Hierarchical Neural Forecasting

Neural network based methods have gained popularity in forecasting
applications, outperforming most alternatives. As surveys show, in recent
years, the academic community has greatly renovated interest in the topic
(Benidis et al., 2020). The literature has permeated into hierarchical
forecasting, with contributions such as SHARQ (Han et al., 2021), HIRED (Paria
et al., 2021), and PROFHIT (Kamarthi et al., 2022) approximate coherent
methods using variants of bottom-up aggregation regularization. Fully coherent
approaches include HierE2E (Rangapuram et al., 2021) a multivariate approach
that incorporates MinTrace-like reconciliation in the network’s optimization,
TDProb (Das et al., 2022) that learns TopDown proportions to probabilistically
reconcile univariate base models.

Despite recent progress in extending neural networks toward hierarchical
forecasting, existing solutions still face challenges: (i) their
implementations rely on restrictive probabilistic assumptions or are not
entirely coherent; (ii) computational complexity of multivariate approaches
scales poorly; (iii) the forecast accuracy improvements over statistical
baselines are still modest.

###  2.3 Mathematical Notation

A hierarchical time series (HTS) is a multivariate time series under
aggregation constraints. We denote the HTS by the vector
$\mathbf{y}_{[i],t}=[\,\mathbf{y}^{\intercal}_{[a],t}\,|\,\mathbf{y}^{\intercal}_{[b],t}\,]^{\intercal}\in\mathbb{R}^{N_{a}+N_{b}}$,
for time step $t$, where $[a],[b]$ denote respectively the aggregate and
bottom level indices. The total number of series in the hierarchy is
$|[i]|=(N_{a}+N_{b})=N_{i}$. We distinguish between the time indices $[t]$ and
forecast indices $\tau\in[t+1:t+h]$, and hierarchical, bottom and aggregate
indexes $\iota\in[i],\beta\in[b]$ , $\alpha\in[a]$.

At any time $t$, the constraints are
$\mathbf{y}_{[a],t}=\mathbf{A}_{[a][b]}\mathbf{y}_{[b],t}$ where
$\mathbf{A}_{[a][b]}$ denotes the relationship between the bottom-level series
to the upper-level series. We can write the HTS as

$\mathbf{y}_{[i],t}=\mathbf{S}_{[i][b]}\mathbf{y}_{[b],t}\quad\Leftrightarrow\quad\begin{bmatrix}\mathbf{y}_{[a],t}\\\
\mathbf{y}_{[b],t}\end{bmatrix}=\begin{bmatrix}\mathbf{A}_{[a][b]}\\\
\mathbf{I}_{[b][b]}\end{bmatrix}\mathbf{y}_{[b],t}$ (1)

where $\mathbf{S}_{[i][b]}$ and $\mathbf{I}_{[b][b]}$ are summing and identity
matrices. Figure 2 exemplifies $N_{i}=7$, $N_{b}=4$ and $N_{a}=3$:

$\begin{split}\mathbf{y}_{[a],t}&=\left[y_{\mathrm{Total},t},\;y_{\beta_{1},t}+y_{\beta_{2},t},\;y_{\beta_{3},t}+y_{\beta_{4},t}\right]^{\intercal},\\\
\mathbf{y}_{[b],t}&=\left[y_{\beta_{1},t},\;y_{\beta_{2},t},\;y_{\beta_{3},t},\;y_{\beta_{4},t}\right]^{\intercal},\end{split}$
(2)

where
$y_{\mathrm{Total},t}=y_{\beta_{1},t}+y_{\beta_{2},t}+y_{\beta_{3},t}+y_{\beta_{4},t}$.
The summing matrix associated to Figure 2 is:

$\mathbf{S}_{[i][b]}=\begin{bmatrix}\\\ \mathbf{A}_{\mathrm{[a][b]}}\\\ \\\
\\\ \mathbf{I}_{\mathrm{[b][b]}}\\\ \\\
\end{bmatrix}=\begin{bmatrix}\;1&1&1&1\\\ \;1&1&0&0\\\ \;0&0&1&1\\\
\hline\cr\;1&0&0&0\\\ \;0&1&0&0\\\ \;0&0&1&0\\\ \;0&0&0&1\\\ \end{bmatrix}.$
(3)

######  Definition 2.1.

(Probabilistic Coherence). Let
$(\Omega_{[b]},\mathcal{F}_{[b]},\mathbb{P}_{[b]})$ be a probabilistic
forecast space, with $\mathcal{F}_{[b]}$ a $\sigma$-algebra on
$\mathbb{R}^{N_{b}}$. Let $\mathbf{S}(\cdot):\Omega_{[b]}\mapsto\Omega_{[i]}$
be the constraints’ implied transformation. A coherent probabilistic forecast
space $(\Omega_{[i]},\mathcal{F}_{[i]},\mathbb{P}_{[i]})$ satisfies:

$\mathbb{P}_{[i]}\left(\mathbf{S}(\mathcal{B})\right)=\mathbb{P}_{[b]}\left(\mathcal{B}\right)\text{for
set }\mathcal{B}\in\mathcal{F}_{[b]}\text{, image
}\mathbf{S}(\mathcal{B})\in\mathcal{F}_{[i]}$ (4)

that is, it assigns a zero probability to any set without coherent forecasts
(Panagiotelis et al., 2023).

######  Definition 2.2.

(Hierarchical Reconciliation). For time $t$, horizon $h$, and forecast indexes
$\tau\in[t+1:t+h]$. Reconciliation for point forecasts
$\hat{\mathbf{y}}_{[i],\tau}$, is denoted by:

$\tilde{\mathbf{y}}_{[i],\tau}=\mathbf{S}_{[i][b]}\mathbf{P}_{[b][i]}\hat{\mathbf{y}}_{[i],\tau}=\mathbf{SP}(\hat{\mathbf{y}}_{[i],\tau})$
(5)

where $\mathbf{P}_{[b][i]}$ is defined by the reconciliation technique. And
$\mathbf{SP}(\cdot):\Omega_{[i]}\mapsto\Omega_{[b]}\mapsto\Omega_{[i]}$ is the
reconciliation’s transformation (Hyndman and Athanasopoulos, 2018).

##  3 Hierarchically Coherent Multivariate Mixture Network Methodology

Figure 3:  The HINT framework leverages a specialized multivariate mixture
probability output that achieves coherency via bootstrap sample reconciliation
$\mathbf{\tilde{y}}_{[i],\tau}=\mathbf{SP}(\mathbf{\hat{y}}_{[i],\tau})$.
Additionally HINT incorporates the TemporalNorm module to augment networks
with scale-invariance through a global skip connection that enables input
normalization and output scale recomposition.

The HINT framework estimates the following conditional probability under
coherency constraints:

$\mathbb{P}\left(\mathbf{y}_{[i][t+1:t+h]}\,|\,\boldsymbol{\theta}\,\right)$
(6)

where $\boldsymbol{\theta}$ depends on
$\\{\mathbf{x}^{h}_{[i][:t]},\;\mathbf{x}^{f}_{[i][:t+h]},\;\mathbf{x}^{s}_{[i]}\\}$
historic, future and static variables. Here we describe our proposed approach,
its main principles of operation are depicted in Figure 3. A HINT network
consists of a coherent probability output and a scale-robustified neural
forecast architecture.

###  3.1 Coherent Multivariate Mixture

HINT is a highly modular system that supports a wide range of probabilistic
outputs. We leverage HINT’s flexibility to accommodate a multivariate Gaussian
mixture model specialized in hierarchical forecasting. Its conditional
forecast distribution is described by:

$\begin{split}\hat{\mathbb{P}}\left(\mathbf{y}_{[i][t+1:t+h]}|\,\hat{\boldsymbol{\theta}}\right)\qquad\qquad\qquad\qquad\qquad\qquad\\\
=\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\prod_{(\iota,\tau)\in[i][t+1:t+h]}\mathcal{N}\left(y_{\iota,\tau}|\;\hat{\mu}_{\iota,\kappa,\tau}\;\hat{\sigma}_{\iota,\kappa,\tau}\right)\end{split}$
(7)

This multivariate mixture has advantageous theoretical properties, proven in
Appendix A. It can arbitrarily approximate univariate distributions and
describe the series’ correlations.

Optimization: HINT achieves high computational efficiency, because we optimize
it through composite-likelihood (Lindsay, 1988; Varin et al., 2011) of the
series in the SGD batches approximating the full joint distribution. Let
$\mathcal{B}=\\{[b_{i}]\\}$ be time-series SGD batches, and
$\boldsymbol{\omega}$ the neural network parameters, then HINT’s negative log
composite likelihood is:

$\begin{split}\mathcal{L}(\boldsymbol{\omega})=-\sum_{[b_{i}]\in\mathcal{B}}\mathrm{log}\Bigg{[}\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}(\boldsymbol{\omega})\prod_{(\iota,\tau)\in[b_{i}][t+1:t+h]}\qquad\qquad\\\
\left(\frac{1}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})\sqrt{2\pi}}\exp{\bigl{\\{}-\frac{1}{2}\left(\frac{y_{\iota,\tau}-\hat{\mu}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}\right)^{2}\bigl{\\}}}\right)\Bigg{]}\end{split}$
(8)

Bootstrap Reconciliation: We ensure the coherence via bootstrap reconciliation
(Panagiotelis et al., 2023). Figure 1.c shows it restoring the aggregation
constraints into the base samples. Let $\mathcal{H}$ be a coherent forecast
set, and $\mathbf{SP}^{-1}(\cdot)$ a reconciliation’s inverse image, Theorem
3.1 analytically derives the coherent forecast distribution:

$\tilde{\mathbb{P}}\left(\tilde{\mathbf{y}}_{[i],\tau}\in\mathcal{H}|\,\tilde{\boldsymbol{\theta}}\right)=\hat{\mathbb{P}}\left(\hat{\mathbf{y}}_{[i],\tau}\in\mathbf{SP}^{-1}(\mathcal{H})|\,\hat{\boldsymbol{\theta}}\right)$
(9)

Prior research has used properties of base forecast distributions and
analytically reconciled probabilities to achieve efficient inference times
(Olivares et al., 2022b).

The _bootstrap reconciliation_ technique offers a distribution-agnostic
reconciliation method. We leverage this flexibility to study HINT with
different distributions outputs in our ablation studies in Figure 7.

######  Theorem 3.1.

Consider a reconciliation from Definition 2.2 where the entire hierarchy’s
forecasts are combined into reconciled bottom-level forecasts using a
composition of linear transformations
$\mathbf{SP}(\cdot)=\mathbf{S_{[i][b]}P_{[b][i]}}(\cdot)$. The reconciled
probability for the entire hierarchical series is given by:

$\begin{split}\tilde{\mathbb{P}}_{[b]}\left(\mathbf{y}_{[b],\tau}\right)=|\mathbf{P}^{*}|\int\hat{\mathbb{P}}_{[i]}\left(\mathbf{P}_{\perp}\mathbf{y}_{[a],\tau}+\mathbf{P}^{-}\mathbf{y}_{[b],\tau}\right)\delta\mathbf{\mathbf{y}_{[a],\tau}}\\\
\tilde{\mathbb{P}}\left(\mathbf{y}_{[i],\tau}\right)=|\mathbf{S}^{*}|\tilde{\mathbb{P}}_{[b]}\left(\mathbf{S}^{-}\mathbf{y}_{[i],\tau}\right)\mathbbm{1}\\{\mathbf{y}_{[i],\tau}\in\mathcal{H}\\}\qquad\qquad\end{split}$
(10)

where $\tilde{\mathbb{P}}_{[b]}(\cdot)$ is the reconciled bottom forecast
distribution, $\mathbbm{1}(\tilde{\mathbf{y}}_{[i],\tau}\in\mathcal{H})$
indicates if realization belongs in the $N_{b}$-dimensional hierarchically
coherent subspace $\mathcal{H}$, $\mathbf{P}^{-}\in\mathbb{R}^{N_{i}\times
N_{b}},\mathbf{S}^{-}\in\mathbb{R}^{N_{b}\times N_{i}}$ are
$\mathbf{P}_{[b][i]},\mathbf{S}_{[i][b]}$ Moore-Penrose inverses and
$\mathbf{P}_{\perp}\in\mathbb{R}^{N_{a}\times
N_{i}},\mathbf{S}_{\perp}\in\mathbb{R}^{N_{i}\times N_{a}}$ their orthogonal
complements,
$\mathbf{P}^{*}=\left[\mathbf{P}_{\perp}\,|\,\mathbf{P}^{-}\right]$, and
$\mathbf{S}^{*}=\left[\mathbf{S}^{-}_{\perp}\,|\,\mathbf{S}\right]^{-1}$.
Proof available in Appendix A.

######  Theorem 3.2.

The multivariate mixture from Equation (7) models the time series’
relationships through the structure of its latent variables, its non-diagonal
covariance is the following $N_{k}-1$ rank matrix:

$\mathrm{Cov}(\mathbf{Y}_{[i],\tau})=\sum^{N_{k}}_{\kappa=1}\hat{\mathbf{w}}_{\kappa}(\hat{\boldsymbol{\mu}}_{[i],\kappa,\tau}-\bar{\boldsymbol{\mu}}_{[i],\tau})(\hat{\boldsymbol{\mu}}_{[i],\kappa,\tau}-\bar{\boldsymbol{\mu}}_{[i],\tau})^{\intercal}$
(11)

Proof in Appendix A.

Figure 5 shows the importance of the flexibility of the multivariate mixture.
When combined with the BottomUp reconciliation strategy, a univariate
independence simplification of Equation (7) translates into overestimating the
aggregate forecast’s uncertainty.

Estimating the series’ relationships can help sharpen the reconciled forecast
distribution. Unlike previous research on high-dimensional multivariate
techniques with covariance matrix low-rank approximations (Salinas et al.,
2019), our distribution does not store the covariance for greatly improved
memory complexity.

$\mathcal{O}(\mathrm{batch\;size})\text{ vs }\mathcal{O}(N^{2}_{i})$ (12)

Figure 4: Temporal normalization (left), layer normalization (center) and
batch normalization (right). The entries in green are in the normalizing
statistics.

###  3.2 Probabilistic Reversible Instance Normalization

Well-performing neural architectures, such as TFT (Lim et al., 2021), DeepAR
(Salinas et al., 2020), and NBEATSx (Oreshkin et al., 2020; Olivares et al.,
2022a), covertly incorporated scale-robustified cross-learning optimization.
Our TemporalNorm module standardizes the approach and makes it available to
other architectures. Kim et al. (2022) recently proposed a general temporal
normalization technique. We expanded on the reversible instance normalization
approach by increasing the available scaling methods and extending them to
probabilistic outputs.

(a) HINT Univariate Mixture

(b) HINT Multivariate Mixture

Figure 5:  HINT univariate/multivariate estimation comparison. Different rows of plots show forecasts for different levels of hierarchical aggregation. Light and dark blue represent 99% and 75% prediction intervals. Modeling the series’ correlations sharpens the reconciled forecast distributions. Table 1:  Empirical evaluation of probabilistically coherent forecasts. Scaled continuous ranked probability score (sCRPS) and relative mean squared error (relMSE), averaged over 10 random seeds, over all hierarchy series. The best results are bold-faced (lower values of metrics used are preferred). † The PROFHIT results differ from Kamarthi et al. 2022, as the only available implementation appears to be numerically unstable in its optimization. $\ddagger$ PatchTST/NBEATS use deterministic coherent normality intervals. * Best performing variant of TopDown (avg. proportions, proportions avg.), and MinTrace (ols, wls, shrinkage) reported. ** PERMBU only available for strict hierarchies.  |  |  HINT (Ours) | Neural Forecast |  ARIMA-BOOTSTRAP |  ARIMA-PERMBU**
---|---|---|---|---|---
| Dataset | Coherent | Base | HierE2E | PROFHIT† | PatchTST$\ddagger$ | NBEATS$\ddagger$ | BottomUp | MinTrace* | BottomUp | MinTrace*
sCRPS | Labour | .0067±.0001 | .0070±.0001 | .0171±.0003 | .2138±.0070 | 0.0094 | 0.0077 | .0078±.0001 | .0073±.0000 | .0077±.0001 | .0069±.0001
Traffic | .0589±.0004 | .0633±.0003 | .0426±.0008 | .1137±.0022 | 0.0872 | 0.0781 | .0736±.0024 | .0608±.0014 | .0849±.0009 | .0651±.0008
Tourism | .0536±.0004 | .0574±.0004 | .0761±.0007 | .1358±.0033 | 0.1093 | 0.0736 | .0682±.0018 | .0703±.0017 | .0649±.0016 | .0680±.0016
Tourism-L | .1176±.0002 | .1253±.0002 | .1424±.0019 | .2139±.0014 | 0.1277 | 0.1291 | .1375±.0013 | .1313±.0009 | - | -
Wiki2 | .2447±.0007 | .2395±.0006 | .2592±.0031 | .4009±.0028 | 0.2999 | 0.3621 | .2894±.0038 | .2808±.0035 | .3920±.0044 | .3821±.0049
relMSE | Labour | .5802±.0131 | .6333±.0141 | .8165±0.0353 | $6.774\times 10^{3}$ | 0.6132 | 0.4932 | 0.5382±.0000 | 0.3547±.0000 |  |
Traffic | .1212±.0051 | .1340±.0062 | .0328±0.0019 | 0.4536±.0224 | 0.1589 | 0.1489 | 0.1392±.0000 | 0.0744±.0000 |  |
Tourism | .0387±.0007 | .0574±.0008 | .1471±0.0046 | 0.9745±.0803 | 0.0671 | 0.0585 | 0.1002±.0000 | 0.1235±.0000 |  |
Tourism-L | .0577±.0009 | .0660±.0010 | .2449±0.0096 | 1.0401±.0296 | 0.0742 | 0.0692 | 0.3070±.0000 | 0.1375±.0000 |  |
Wiki2 | .1884±.0012 | .1966±0.0017 | .6598±0.0249 | 0.7901±.0384 | 0.2014 | 0.2137 | 1.0163±.0000 | 1.0068±.0000 |  |

##  4 Experimental Results

Hierarchical Forecast Datasets. We follow experimental protocols established
in previous research by Rangapuram et al. (2021). The benchmark datasets are:
Monthly Australian Labour (Australian Bureau of Statistics, 2019), SF Bay Area
daily Traffic (Dua and Graff, 2017), Quarterly Australian Tourism visits
(Tourism Australia, Canberra, 2005), Monthly Australian Tourism-L visits
(Tourism Australia, Canberra, 2019), and daily Wiki2 views (Anava et al.,
2018). Appendix B includes a detailed exploration of this data.

Baselines. In our main experiment, we compare with SoTA probabilistic coherent
methods. Neural forecasting baselines include (1) HierE2E (Rangapuram et al.,
2021), (2) PROFHIT (Kamarthi et al., 2022), (3) PatchTST (Nie et al., 2023),
(4) NBEATSx (Oreshkin et al., 2020; Olivares et al., 2022a), while statistical
baselines include variants of (5) BOOTSTRAP (Panagiotelis et al., 2023), and
(6) PERMBU probabilistic reconciliation (Ben Taieb et al., 2017) in
combination with BottomUp (Orcutt et al., 1968), and MinTrace (Wickramasuriya
et al., 2019) reconcilers. We use HierarchicalForecast library baselines
implementation (Olivares et al., 2022b).

Evaluation Metrics. To assess the forecast accuracy of our method, we compute
the scaled Continuous Ranked Probability Score (sCRPS; Matheson and Winkler
1976; Makridakis et al. 2022) and the Relative Mean Squared Error (relMSE;
Hyndman and Koehler 2006; Olivares et al. 2023). Note that probabilistic
coherence naturally implies mean hierarchical coherence.

$\begin{split}\mathrm{sCRPS}(\mathbb{P},\mathbf{y}_{[i],\tau})=\frac{2}{|[i\,]|}\sum_{i}\frac{\int^{1}_{0}\mathrm{QL}(\mathbb{P}_{i,\tau},y_{i,\tau})_{q}dq}{\sum_{i}|y_{i,\tau}|}\qquad\\\
\mathrm{relMSE}(\mathbf{y}_{[i]},\hat{\mathbf{y}}_{[i]},\mathbf{\check{y}}_{[i]})=\frac{\mathrm{MSE}(\mathbf{y}_{[i]},\mathbf{\hat{y}}_{[i]})}{\mathrm{MSE}(\mathbf{y}_{[i]},\mathbf{\check{y}}_{[i]})}\qquad\qquad\qquad\end{split}$
(13)

where $\mathrm{QL}(\hat{\mathbb{P}}_{i,\tau},y_{i,\tau})_{q}$ stands for the
q-quantile loss.

(a) _Example of a Coherent Multivariate Mixture_

(b) _Contributions of the Mixture Components_

Figure 6: (a) Modularity of HINT can accommodate various distribution outputs
and allows us to leverage hierarchically coherent multivariate mixture output
module. (b) Validation sCRPS results as a function of mixture components. The
sCRPS performance curves of HINT exhibit a bias-variance tradeoff shape, with
an optimal number of components (10).

HINT configurations. We showcase HINT’s modularity by augmenting three well-
established neural forecast architectures: NHITS (Oreshkin et al., 2020;
Challu et al., 2023), TFT (Lim et al., 2021), and TCN (Bai et al., 2018). We
optimally select HINT’s underlying neural architecture using temporal cross-
validation; Appendix D describes analyzed configurations and hyperparameter
selection process. The experiments involve two alternative scenarios: the
first is a base HINT with no reconciliation with no guaranteed coherent
forecasts, the second is HINT with end-to-end reconciliation.

###  4.1 Empirical Results

Main Results: In Table 1, HINT demonstrates superior performance across all
datasets, except for Traffic where it ranks as the second-best, behind only
HierE2E. We execute all methods eight times and present its average and
standard deviation of the sCRPS and relMSE scores. Overall HINT improves sCRPS
by 13.2% on average on all datasets except Traffic. The relMSE results are
highly correlated, although relMSE is more susceptible to outliers. We extend
these results in Appendix E where we explore the performance of the methods
across levels of hierarchical aggregation.

Ablation Studies: We conducted ablation studies on different variants of our
model. Appendix C.1, examines the impact of the TemporalNorm module. We
observe that scale-decoupled optimization can improve accuracy of neural
forecasting by an order of magnitude. Figure 8 also shows how robust scaling
significantly improves over vanilla REVIN (Kim et al., 2022).

Figure 7: Ablation study of forecast distribution’s accuracy in the Tourism-L
dataset. Our multivariate mixture model translates into clear accuracy gains.

Appendix C.2 ablation study guides our selection for the number of mixture
components, as shown in Figure 6, there is a clear advantage in the usage of a
flexible multivariate mixture, as our method is both capable of improving over
restrictive distribution assumptions and improve the forecast accuracy by
capturing the relationships among series in the hierarchy.

Figure 8: Ablation study of scaling variants in the TemporalNorm module and
its effects on forecasts’ accuracy in Tourism-L dataset. In the presence of
noise, robust scaling improves over REVIN baseline.

Our ablation studies place the origins of the observed improvements in the
HINT framework: (i) the mixture distribution improves performance upon simpler
and more restrictive probabilistic output assumptions. (ii) scale-robustified
optimization is a crucial enabler of cross-learning under stark scale
variations; (iii) advanced reconciliation techniques can overcome potential
systemic biases of base distributions; (iv) the latest neural forecast
architectural innovations tend to improve the forecast accuracy.

##  5 Future Research and Discussion

The Hierarchically Coherent Multivariate Mixture Network (HINT) framework aims
to learn the joint distribution for the entire hierarchy, while guaranteeing
its coherence via bootstrap reconciliation. The proposed framework is highly
modular and enables many possible extensions, some of which we briefly discuss
below.

Concerning hierarchical reconciliation, while most hierarchical forecast
methods rely on post hoc techniques, the potential to learn the reconciliation
method using a machine learning framework has received relatively little
attention. However, initial findings along this line of investigation show
very promising results (Ben Taieb and Koo, 2019; Das et al., 2023).

In our experiments, the only case where HINT did not outperform the HierE2E
(Rangapuram et al., 2021) baseline was the Traffic dataset. HierE2E
outperforms all the other methods on those data due to including all the
series as predictor variables in a VAR-like multivariate forecasting system
under clear Granger causal relationships. Our results using univariate inputs
challenge the effectiveness of existing multivariate input hierarchical
forecasting approaches; the design of algorithms capable of maintaining good
performance in the absence of Granger causalities while leveraging them, if
present, is an exciting line of future research.

The demand for large-scale hierarchical forecasting systems is increasing
across multiple important application domains. However, existing methods
typically rely on access to the entire hierarchical structure for
reconciliation. Innovative variations of the TopDown and BottomUp
reconciliation approaches could be employed to develop a parallelizable
approach that could leverage distributed computing frameworks.

##  6 Conclusion

We introduced HINT, a new framework that augments proven neural forecasting
architectures to generate probabilistically coherent hierarchical forecasts.
Two highly modular novelties are foundations of our approach: (i) a highly
efficient multivariate mixture optimized with composite likelihood and
transformed via bootstrap sample reconciliation, and (ii) incorporation of
robust feature extraction and recomposition of probabilistic output scales
within the network’s architecture.

We empirically demonstrated that the flexibility of our multivariate mixture
output can yield significant accuracy improvements (13.2% on average on
representative collection of benchmark data) compared to state-of-the-art
baselines. These findings highlight the effectiveness of the presented
approach in producing accurate and coherent forecasts and pave the way for
further advancements in hierarchical forecasting.

HINT is open-sourced, and practitioners can easily incorporate this framework
into their work by accessing the code under this link http URL.

## References

  * Anava et al. (2018) Oren Anava, Vitaly Kuznetsov, and (Google Inc. Sponsorship).  Web traffic time series forecasting, forecast future traffic to wikipedia pages.  Kaggle Competition, 2018.  URL https://www.kaggle.com/c/web-traffic-time-series-forecasting/.
  * Australian Bureau of Statistics (2019) Australian Bureau of Statistics.  Labour force, australia.  Accessed Online, 2019.  URL https://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/6202.0Dec%202019?OpenDocument.
  * Babai et al. (2022) M. Zied Babai, John E. Boylan, and Bahman Rostami-Tabar.  Demand forecasting in supply chains: a review of aggregation and hierarchical approaches.  _International Journal of Production Research_ , 60(1):324–348, 2022.  doi: 10.1080/00207543.2021.2005268.  URL https://doi.org/10.1080/00207543.2021.2005268.
  * Bai et al. (2018) Shaojie Bai, J. Zico Kolter, and Vladlen Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  _Computing Research Repository_ , abs/1803.01271, 2018.  URL http://arxiv.org/abs/1803.01271.
  * Ben Taieb and Koo (2019) Souhaib Ben Taieb and Bonsoo Koo.  Regularized regression for hierarchical forecasting without unbiasedness conditions.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’19, page 1337–1347, New York, NY, USA, 2019. Association for Computing Machinery.  ISBN 9781450362016.  doi: 10.1145/3292500.3330976.  URL https://doi.org/10.1145/3292500.3330976.
  * Ben Taieb et al. (2017) Souhaib Ben Taieb, James W. Taylor, and Rob J. Hyndman.  Coherent probabilistic forecasts for hierarchical time series.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 3348–3357. PMLR, 06–11 Aug 2017.  URL http://proceedings.mlr.press/v70/taieb17a.html.
  * Ben Taieb et al. (2021) Souhaib Ben Taieb, James W. Taylor, and Rob J. Hyndman.  Hierarchical probabilistic forecasting of electricity demand with smart meter data.  _Journal of the American Statistical Association_ , 116(533):27–43, 2021.  doi: 10.1080/01621459.2020.1736081.  URL https://doi.org/10.1080/01621459.2020.1736081.
  * Benidis et al. (2020) Konstantinos Benidis, Syama Sundar Rangapuram, Valentin Flunkert, Bernie Wang, Danielle Maddix, Caner Turkmen, Jan Gasthaus, Michael Bohlke-Schneider, David Salinas, Lorenzo Stella, Laurent Callot, and Tim Januschowski.  Neural forecasting: Introduction and literature overview.  _Computing Research Repository_ , 2020.
  * Bergstra et al. (2011) James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.  Algorithms for hyper-parameter optimization.  In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 24, pages 2546–2554. Curran Associates, Inc., 2011.  URL https://proceedings.neurips.cc/paper/2011/file/86e8f7ab32cfd12577bc2619bc635690-Paper.pdf.
  * Challu et al. (2023) Cristian Challu, Kin G. Olivares, Boris N. Oreshkin, Federico Garza, Max Mergenthaler, and Artur Dubrawski.  NHITS: Neural Hierarchical Interpolation for Time Series forecasting.  In _The Association for the Advancement of Artificial Intelligence Conference 2023 (AAAI 2023)_ , 2023.  URL https://arxiv.org/abs/2201.12886.
  * Das et al. (2022) Abhimanyu Das, Weihao Kong, Biswajit Paria, and Rajat Sen.  A deep top-down approach to hierarchically coherent probabilistic forecasting, 2022.  URL https://arxiv.org/abs/2204.10414.
  * Das et al. (2023) Abhimanyu Das, Weihao Kong, Biswajit Paria, and Rajat Sen.  Dirichlet proportions model for hierarchically coherent probabilistic forecasting.  In _The 39th Conference on Uncertainty in Artificial Intelligence_ , 2023.  URL https://openreview.net/forum?id=O11qR0Qrcc.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Dunn et al. (1976) D. M. Dunn, W. H. Williams, and T. L. Dechaine.  Aggregate versus subaggregate models in local area forecasting.  _Journal of the American Statistical Association_ , 71(353):68–71, 1976.
  * Eckert et al. (2021) Florian Eckert, Rob J. Hyndman, and Anastasios Panagiotelis.  Forecasting swiss exports using bayesian forecast reconciliation.  _European Journal of Operational Research_ , 291(2):693–710, 2021.  ISSN 0377-2217.  doi: https://doi.org/10.1016/j.ejor.2020.09.046.  URL https://www.sciencedirect.com/science/article/pii/S037722172030850X.
  * Fliedner (1999) Gene Fliedner.  An investigation of aggregate variable time series forecast strategies with specific subaggregate time series statistical correlation.  _Computers and Operations Research_ , 26(10–11):1133–1149, September 1999.  ISSN 0305-0548.  doi: 10.1016/S0305-0548(99)00017-9.  URL https://doi.org/10.1016/S0305-0548(99)00017-9.
  * Garza et al. (2022) Federico Garza, Max Mergenthaler Canseco, Cristian Challú, and Kin G. Olivares.  StatsForecast: Lightning fast forecasting with statistical and econometric models.  PyCon Salt Lake City, Utah, US 2022, 2022.  URL https://github.com/Nixtla/statsforecast.
  * Gross and Sohl (1990) Charles W. Gross and Jeffrey E. Sohl.  Disaggregation methods to expedite product line forecasting.  _Journal of Forecasting_ , 9(3):233–254, 1990.  doi: 10.1002/for.3980090304.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/for.3980090304.
  * Han et al. (2021) Xing Han, Sambarta Dasgupta, and Joydeep Ghosh.  Simultaneously reconciled quantile forecasting of hierarchically related time series.  In Arindam Banerjee and Kenji Fukumizu, editors, _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pages 190–198. PMLR, 13–15 Apr 2021.  URL http://proceedings.mlr.press/v130/han21a.html.
  * Hyndman and Athanasopoulos (2018) Rob J Hyndman and George Athanasopoulos.  _Forecasting: Principles and Practice_.  OTexts, Melbourne, Australia, 2018.  available at https://otexts.com/fpp2/.
  * Hyndman and Khandakar (2008) Rob J. Hyndman and Yeasmin Khandakar.  Automatic time series forecasting: The forecast package for r.  _Journal of Statistical Software, Articles_ , 27(3):1–22, 2008.  ISSN 1548-7660.  doi: 10.18637/jss.v027.i03.  URL https://www.jstatsoft.org/v027/i03.
  * Hyndman and Koehler (2006) Rob J. Hyndman and Anne B. Koehler.  Another look at measures of forecast accuracy.  _International Journal of Forecasting_ , 22(4):679 – 688, 2006.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2006.03.001.  URL http://www.sciencedirect.com/science/article/pii/S0169207006000239.
  * Hyndman et al. (2011) Rob J. Hyndman, Roman A. Ahmed, George Athanasopoulos, and Han Lin Shang.  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 http://www.sciencedirect.com/science/article/pii/S0167947311000971.
  * Kamarthi et al. (2022) Harshavardhan Kamarthi, Lingkai Kong, Alexander Rodriguez, Chao Zhang, and B. Prakash.  PROFHIT: Probabilistic robust forecasting for hierarchical time-series.  _Computing Research Repository_ , 06 2022.  URL https://arxiv.org/abs/2206.07940.
  * Kim et al. (2022) Taesung Kim, Jinhee Kim, Yunwon Tae, Cheonbok Park, Jang-Ho Choi, and Jaegul Choo.  Reversible instance normalization for accurate time-series forecasting against distribution shift.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=cGDAkQo1C0p.
  * Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba.  ADAM: A method for stochastic optimization, 2014.  URL http://arxiv.org/abs/1412.6980.  cite arxiv:1412.6980Comment: Published as a conference paper at the 3rd International Conference for Learning Representations (ICLR), San Diego, 2015.
  * Kourentzes and Athanasopoulos (2019) Nikolaos Kourentzes and George Athanasopoulos.  Cross-temporal coherent forecasts for australian tourism.  _Annals of Tourism Research_ , 75:393–409, 2019.  ISSN 0160-7383.  doi: https://doi.org/10.1016/j.annals.2019.02.001.  URL https://www.sciencedirect.com/science/article/pii/S0160738319300167.
  * Lim et al. (2021) Bryan Lim, Sercan O. Arık, Nicolas Loeff, and Tomas Pfister.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  _International Journal of Forecasting_ , 37(4):1748–1764, 2021.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2021.03.012.  URL https://www.sciencedirect.com/science/article/pii/S0169207021000637.
  * Lindsay (1988) B. G. Lindsay.  Composite likelihood methods.  _Contemporary Mathematics_ , 80:221–239, 1988.
  * Makridakis et al. (2022) Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos, Zhi Chen, Anil Gaba, Ilia Tsetlin, and Robert L. Winkler.  The m5 uncertainty competition: Results, findings and conclusions.  _International Journal of Forecasting_ , 38(4):1365–1385, 2022.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2021.10.009.  URL https://www.sciencedirect.com/science/article/pii/S0169207021001722.  Special Issue: M5 competition.
  * Matheson and Winkler (1976) James E. Matheson and Robert L. Winkler.  Scoring rules for continuous probability distributions.  _Management Science_ , 22(10):1087–1096, 1976.  ISSN 00251909, 15265501.  URL http://www.jstor.org/stable/2629907.
  * Nguyen and McLachlan (2018) Hien D. Nguyen and Geoffrey J. McLachlan.  On approximations via convolution-defined mixture models.  _Computing Research Repository_ , 2018.  URL https://arxiv.org/abs/1611.03974.
  * Nie et al. (2023) Yuqi Nie, Nam H. Nguyen, Phanwadee Sinthong, and Jayant Kalagnanam.  A time series is worth 64 words: Long-term forecasting with transformers, 2023.
  * Olivares et al. (2022a) Kin G. Olivares, Cristian Challu, Grzegorz Marcjasz, Rafał Weron, and Artur Dubrawski.  Neural basis expansion analysis with exogenous variables: Forecasting electricity prices with NBEATSx.  _International Journal of Forecasting_ , 2022a.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2022.03.001.  URL https://www.sciencedirect.com/science/article/pii/S0169207022000413.
  * Olivares et al. (2022b) Kin G. Olivares, Federico Garza, David Luo, Cristian Challú, Max Mergenthaler, Souhaib Ben Taieb, Shanika L. Wickramasuriya, and Artur Dubrawski.  HierarchicalForecast: A reference framework for hierarchical forecasting.  _Journal of Machine Learning Research, submitted_ , abs/2207.03517, 2022b.  URL https://arxiv.org/abs/2207.03517.
  * Olivares et al. (2023) Kin G. Olivares, Nganba Meetei, Ruijun Ma, Rohan Reddy, Mengfei Cao, and Lee Dicker.  Probabilistic hierarchical forecasting with deep poisson mixtures.  _International Journal of Forecasting, accepted_ , Preprint version available at arXiv:2110.13179, 2023.  URL https://arxiv.org/abs/2110.13179.
  * Orcutt et al. (1968) Guy H. Orcutt, Harold W. Watts, and John B. Edwards.  Data aggregation and information loss.  _The American Economic Review_ , 58(4):773–787, 1968.  ISSN 00028282.  URL http://www.jstor.org/stable/1815532.
  * Oreshkin et al. (2020) Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-BEATS: neural basis expansion analysis for interpretable time series forecasting.  In _8th International Conference on Learning Representations, ICLR 2020_ , 2020.  URL https://openreview.net/forum?id=r1ecqn4YwB.
  * Panagiotelis et al. (2023) Anastasios Panagiotelis, Puwasala Gamakumara, George Athanasopoulos, and Rob J. Hyndman.  Probabilistic forecast reconciliation: Properties, evaluation and score optimisation.  _European Journal of Operational Research_ , 306(2):693–706, 2023.  ISSN 0377-2217.  doi: https://doi.org/10.1016/j.ejor.2022.07.040.  URL https://www.sciencedirect.com/science/article/pii/S0377221722006087.
  * Paria et al. (2021) Biswajit Paria, Rajat Sen, Amr Ahmed, and Abhimanyu Das.  Hierarchically Regularized Deep Forecasting.  In _Submitted to Proceedings of the 39th International Conference on Machine Learning_. PMLR. Working Paper version available at arXiv:2106.07630, 2021.
  * Paszke et al. (2019) Paszke et al.  Pytorch: An imperative style, high-performance Deep Learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Rangapuram et al. (2021) Syama Sundar Rangapuram, Lucien D. Werner, Konstantinos Benidis, Pedro Mercado, Jan Gasthaus, and Tim Januschowski.  End-to-end learning of coherent probabilistic forecasts for hierarchical time series.  In Maria Florina Balcan and Marina Meila, editors, _Proceedings of the 38th International Conference on Machine Learning_ , Proceedings of Machine Learning Research. PMLR, 06–11 Aug 2021.
  * Salinas et al. (2019) David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, Jan Gasthaus, and Roberto Medico.  High-dimensional multivariate forecasting with low-rank gaussian copula processes.  In _NeurIPS 2019_ , 2019.  URL https://www.amazon.science/publications/high-dimensional-multivariate-forecasting-with-low-rank-gaussian-copula-processes.
  * Salinas et al. (2020) David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 36(3):1181–1191, 2020.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2019.07.001.  URL https://www.sciencedirect.com/science/article/pii/S0169207019301888.
  * Titterington et al. (1985) D Michael Titterington, Smith Afm, Adrian FM Smith, UE Makov, et al.  _Statistical analysis of finite mixture distributions_ , volume 198.  John Wiley & Sons Incorporated, 1985.
  * Tourism Australia, Canberra (2005) Tourism Australia, Canberra.  Tourism Research Australia (2005), Travel by Australians.  https://www.kaggle.com/luisblanche/quarterly-tourism-in-australia/, 2005.
  * Tourism Australia, Canberra (2019) Tourism Australia, Canberra.  Detailed tourism Australia (2005), Travel by Australians, Sep 2019.  Accessed at https://robjhyndman.com/publications/hierarchical-tourism/.
  * Varin et al. (2011) Cristiano Varin, Nancy Reid, and David Firth.  An overview of composite likelihood methods.  _Statistica Sinica_ , 21(1):5–42, 2011.  ISSN 10170405, 19968507.  URL http://www.jstor.org/stable/24309261.
  * Wickramasuriya (2023) Shanika L. Wickramasuriya.  Probabilistic forecast reconciliation under the Gaussian framework.  _Accepted at Journal of Business and Economic Statistics_ , 2023.
  * Wickramasuriya et al. (2019) Shanika L. Wickramasuriya, George Athanasopoulos, and Rob J. Hyndman.  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://robjhyndman.com/publications/mint/.
  * Yao et al. (2007) Yuan Yao, Lorenzo Rosasco, and Caponneto Andrea.  On early stopping in gradient descent learning.  _Constructive Approximation_ , 26(2):289–315, 2007.  URL https://doi.org/10.1007/s00365-006-0663-2.

##  Appendix A Coherent Multivariate Mixture Properties

###  A.1 Base Multivariate Mixture Probability

|
---|---

(a) Gaussian Mixture (b) Poisson Mixture

Figure 1: HINT’s multivariate joint distribution has advantageous properties
that make it uniquely suited for hierarchical forecasting. It is highly
flexible, capable of efficiently modeling series’ relationships, and under
minimal restrictions, guarantees probabilistic coherence.

In Section 3, we highlight that HINT boasts a flexible and modular framework
that can handle various probabilistic outputs. We use this flexibility to
enhance the networks with a multivariate mixture density. Specifically, the
mixture model describes the joint probability of the hierarchical time series
$\mathbf{Y}_{[i][t+1:t+h]}$ as follows:

$\mathbb{P}\left(Y_{[i][t+1:t+h]}=\mathbf{y}_{[i][t+1:t+h]}|\,\hat{\boldsymbol{\theta}}\right)=\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\prod_{(\iota,\tau)\in[i][t+1:t+h]}\mathcal{N}\left(y_{\iota,\tau}|\;\hat{\mu}_{\iota,\kappa,\tau}\;\hat{\sigma}_{\iota,\kappa,\tau}\right)$

where the mixture describes individual series through the location and
variance parameters. For simplicity we denote the combined parameters
$\hat{\boldsymbol{\theta}}_{[i][k][t+1:t+h]}=[\hat{\boldsymbol{\mu}}_{[i][k][t+1:t+h]}\;|\;\hat{\boldsymbol{\sigma}}_{[i][k][t+1:t+h]}]$.

Under reasonable assumptions for the underlying probability, the mixture
distribution offers arbitrary approximation guarantees (Titterington et al.,
1985; Nguyen and McLachlan, 2018). We can control its flexibility by adjusting
the number of components $|[k]|=N_{k}$. Furthermore, the mixture is not
limited to Gaussian components; we can extend it to include discrete
variables. Figure 1 presents an example of its marginal probabilities.

Conditional Independence: A key consequence of the multivariate mixture
probability in Equation (A.1), is the assumption that the modeled series
$\mathbf{y}_{[i][t+1:t+h]}$ are conditionally independent given the parameters
$\hat{\boldsymbol{\mu}}_{[i][k][t+1:t+h]}$ and
$\hat{\boldsymbol{\sigma}}_{[i][k][t+1:t+h]}$. That is for any series and
horizons $(\iota,\tau)\neq(\iota^{\prime},\tau^{\prime})\text{,
}(\iota,\tau),(\iota^{\prime},\tau^{\prime})\in[b][t+1:t+h]$ and
$\kappa\in[k]$:

$\mathbb{P}(Y_{\iota,\tau},Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota,\kappa,\tau},\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})=\mathbb{P}(Y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau})\hat{\mathbb{P}}(Y_{\iota^{\prime},\tau^{\prime}},|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})$
(14)

Computational Efficiency: To handle large-scale data scenarios, we explicitly
avoid using a multivariate covariance matrix, which has an
$\mathcal{O}(N_{i}^{2})$ complexity. Instead, we rely on the mixture latent
variables $\kappa$ and its associated weights
$\hat{\mathbf{w}}_{[k]}\in[0,1]^{N_{k}},\;\hat{\mathbf{w}}_{[k]}\geq 0$ and
$\sum^{N_{k}}_{\kappa=1}\hat{w}_{\kappa}=1$ to model the series correlations.
We show the relationship between the mixture components and the covariance in
Appendix A.3. Instead of relying on the Markov assumption, we have adopted a
joint multi-step forecasting approach that can significantly enhance the
computational efficiency of our algorithm. By making predictions in a single
forward pass, we can avoid the need for recurrent computations.

In the following subsections, we delve into the properties of the multivariate
mixture, such as the analytic version of its implied marginal probability, the
relationship between its covariance and the number of mixture components, the
bootstrap sample reconciled probability, and the optimization strategies that
make it well-suited for large-scale applications. The proofs are inspired on
previous work by (Olivares et al., 2023), generalized and extended.

###  A.2 Marginal Distributions

We define the joint distribution of all hierarchical time series in Equation
(A.1). By integrating the joint probability on the remaining series and time
indices, we can obtain the marginal distribution for a single future horizon
$\tau\in[t+1:t+h]$ and series $\iota\in[i]$. We express the resulting marginal
distribution as follows:

$\mathbb{P}(Y_{\iota,\tau}=y_{\iota,\tau}|\,\hat{\boldsymbol{\theta}})=\sum_{\kappa=1}^{N_{k}}w_{\kappa}\mathcal{N}\left(y_{\iota,\tau}|\;\hat{\mu}_{\iota,\kappa,\tau}\;\hat{\sigma}_{\iota,\kappa,\tau}\right)$
(15)

###### Proof.

$\begin{split}\mathbb{P}(Y_{\iota,\tau}=y_{\iota,\tau}|\,\hat{\boldsymbol{\theta}})&=\int^{+\infty}_{-\infty}\ldots\int^{+\infty}_{-\infty}\mathbb{P}\left(Y_{[i][t+1:t+h]}=\mathbf{y}_{[i][t+1:t+h]}|\,\hat{\boldsymbol{\theta}}\right)\delta
y_{\iota^{\prime},\tau^{\prime}}\setminus\delta y_{\iota,\tau}\\\
&=\int^{+\infty}_{-\infty}\ldots\int^{+\infty}_{-\infty}\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\;\mathbb{P}(y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau})\times\prod_{(\iota^{\prime},\tau^{\prime})\in[i][t+1:t+h]\setminus(\iota,\tau)}\quad\sum_{y_{\iota^{\prime},\tau^{\prime}}}\mathbb{P}(y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})\delta
y_{\iota^{\prime},\tau^{\prime}}\\\
&=\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\;\mathbb{P}(y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau})\times\int^{+\infty}_{-\infty}\ldots\int^{+\infty}_{-\infty}\prod_{(\iota^{\prime},\tau^{\prime})\in[i][t+1:t+h]\setminus(\iota,\tau)}\quad\sum_{y_{\iota^{\prime},\tau^{\prime}}}\mathbb{P}(y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})\delta
y_{\iota^{\prime},\tau^{\prime}}\\\
&=\sum_{\kappa=1}^{N_{k}}w_{\kappa}\;\mathbb{P}(y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau})=\sum_{\kappa=1}^{N_{k}}w_{\kappa}\mathcal{N}\left(y_{\iota,\tau}|\;\hat{\mu}_{\iota,\kappa,\tau}\;\hat{\sigma}_{\iota,\kappa,\tau}\right)\end{split}$

∎

By removing all other time series and forecast horizons from the joint
probability product, the conditional independence expressed in Equation A.1
can efficiently generate forecast distributions for individual variables.

###  A.3 Efficient Covariance Matrix Low-Rank Approximation

Due to the computational challenges of estimating high-dimensional covariance
matrices, existing multivariate methods are limited in their ability to handle
a large number of series. To overcome this challenge, our method utilizes a
low-rank covariance structure implied by the latent variables of the mixture
probability, thereby avoiding the need to compute the covariance matrix
explicitly. By doing so, we significantly reduce the number of parameters and
enable the modeling of time-varying correlations across millions of time
series.

Let a multivariate random variable
$\mathbf{Y}_{[i][t+1:t+h]}\in\mathbb{R}^{N{i}\times h}$ distribution be
described by the mixture from Equation (A.1), the non-diagonal terms of its
implied covariance is the following $N_{k}-1$ rank matrix:

$\mathrm{Cov}(\mathbf{Y}_{[i],\tau})=\sum^{N_{k}}_{\kappa=1}\hat{\mathbf{w}}_{\kappa}(\hat{\boldsymbol{\mu}}_{[i],\kappa,\tau}-\bar{\boldsymbol{\mu}}_{[i],\tau})(\hat{\boldsymbol{\mu}}_{[i],\kappa,\tau}-\bar{\boldsymbol{\mu}}_{[i],\tau})^{\intercal}$
(16)

###### Proof.

We will start by showing that for a pair of series the covariance function is
given by:

$\mathrm{Cov}(Y_{\iota,\tau},Y_{\iota^{\prime},\tau^{\prime}})=\overline{\sigma}_{\iota,\tau}\mathbbm{1}(\iota=\iota^{\prime})\mathbbm{1}(\tau=\tau^{\prime})+\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\left(\hat{\mu}_{\iota,\kappa,\tau}-\overline{\mu}_{\iota,\tau}\right)\left(\hat{\mu}_{\iota^{\prime},\kappa,\tau^{\prime}}-\overline{\mu}_{\iota^{\prime},\tau^{\prime}}\right)$
(17)

By the law of total covariance:

$\mathrm{Cov}(Y_{\iota,\tau},Y_{\iota^{\prime},\tau^{\prime}})=\mathbb{E}\left[\mathrm{Cov}(Y_{\iota,\tau},Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota,\kappa,\tau},\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})\right]+\mathrm{Cov}\left(\mathrm{E}\left[Y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau}\right],\mathbb{E}\left[Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}}\right]\right)$

Using the conditional independence from Equation (14). We can rewrite
conditional covariance expectation:

$\begin{split}\mathbb{E}\left[\mathrm{Cov}(Y_{\iota,\tau},Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota,\kappa,\tau},\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}})\right]&=\mathbb{E}\left[\mathrm{Var}(Y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau})\right]\mathbbm{1}(\iota=\iota^{\prime})\mathbbm{1}(\tau=\tau^{\prime})\\\
&=\mathbb{E}\left[\hat{\sigma}_{\iota,\kappa,\tau}\right]\mathbbm{1}(\iota=\iota^{\prime})\mathbbm{1}(\tau=\tau^{\prime})\\\
&=\overline{\sigma}_{\iota,\tau}\mathbbm{1}(\iota=\iota^{\prime})\mathbbm{1}(\tau=\tau^{\prime})\end{split}$

where
$\overline{\sigma}_{\iota,\tau}=\mathbb{E}\left[\hat{\sigma}_{\iota,\kappa,\tau}\right]=\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\hat{\sigma}_{\iota,\kappa,\tau}$.

In the second term, because the conditional distributions are Normal we have

$\mathrm{E}\left[Y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau}\right]=\hat{\mu}_{\iota,\kappa,\tau}\quad\text{
and
}\quad\mathrm{E}\left[Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}}\right]=\hat{\mu}_{\iota^{\prime},\kappa,\tau^{\prime}}$

Which implies

$\mathrm{Cov}\left(\mathbb{E}\left[Y_{\iota,\tau}|\hat{\theta}_{\iota,\kappa,\tau}\right],\mathbb{E}\left[Y_{\iota^{\prime},\tau^{\prime}}|\hat{\theta}_{\iota^{\prime},\kappa,\tau^{\prime}}\right]\right)=\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}\left(\hat{\mu}_{\iota,\kappa,\tau}-\bar{\mu}_{\iota,\tau}\right)\left(\hat{\mu}_{\iota^{\prime},\kappa,\tau^{\prime}}-\bar{\mu}_{\iota^{\prime},\tau}\right)$

Combining the two partial results we recover the pair-wise covariance formula
in Equation (17), which can be easily extended to the multivariate case from
Equation (16). The rank of the matrix can be infered by observing that
Equation (16) is the sum of $N_{k}$ vectors centered around their means. ∎

###  A.4 Bootstrap Reconciled Probabilities

Let
$(\Omega_{[i]},\mathcal{F}_{[i]},\hat{\mathbb{P}}(\cdot\,|\,\hat{\boldsymbol{\theta}}))$
be a probabilistic forecast space, with $\mathcal{F}_{[i]}$ a $\sigma$-algebra
on $\mathbb{R}^{N_{i}}$. Let a hierarchical reconciliation transformation be
denoted by
$\mathbf{SP}(\cdot):\Omega_{[i]}\mapsto\Omega_{[b]}\mapsto\Omega_{[i]}$.
Consider $\hat{\mathbf{y}}^{s}_{[i],\tau},\;s=1,\dots,S$ samples drawn from an
unconstrained base probability
$\hat{\mathbb{P}}(\cdot\,|\,\hat{\boldsymbol{\theta}})$, and the transformed
samples
$\tilde{\mathbf{y}}^{s}_{[i],\tau}=\mathbf{SP}\left(\hat{\mathbf{y}}^{s}_{[i],\tau}\right)$.

The probability distribution of the reconciled samples is given by:

$\tilde{\mathbb{P}}\left(\tilde{\mathbf{y}}_{[i],\tau}\in\mathcal{H}|\,\tilde{\boldsymbol{\theta}}\right)=\hat{\mathbb{P}}\left(\hat{\mathbf{y}}_{[i],\tau}\in\mathbf{SP}^{-1}(\mathcal{H})|\,\hat{\boldsymbol{\theta}}\right)$
(18)

with $\mathcal{H}$ be a coherent forecast measurable set, and
$\mathbf{SP}^{-1}(\cdot)$ the reconciliation’s inverse image.

###### Proof.

This proof makes only minor modifications to the arguments presented in
Panagiotelis et al. (2023).

$\begin{split}\tilde{\mathbb{P}}\left(\tilde{\mathbf{y}}_{[i],\tau}\in\mathcal{H}|\,\tilde{\boldsymbol{\theta}}\right)&=\lim_{S\to\infty}\frac{1}{S}\sum^{S}_{s=1}\mathbbm{1}\\{\tilde{\mathbf{y}}_{[i],\tau}\in\mathcal{H}\\}\\\
&=\lim_{S\to\infty}\frac{1}{S}\sum^{S}_{s=1}\mathbbm{1}\\{\hat{\mathbf{y}}_{[i],\tau}\in\mathbf{SP}^{-1}(\mathcal{H})\\}\\\
&=\hat{\mathbb{P}}\left(\hat{\mathbf{y}}_{[i],\tau}\in\mathbf{SP}^{-1}(\mathcal{H})|\,\hat{\boldsymbol{\theta}}\right)\end{split}$
(19)

The first and final equalities follow from the weak law of large numbers, as
by definition the indicator functions are independent, identically distributed
with a finite mean and variance. The second equality follows from the
definition of the inverse image
$\mathbf{SP}^{-1}(\mathcal{H})=\\{\hat{y}_{[i],\tau}\,|\,\mathbf{SP}(\hat{y}_{[i],\tau})\in\mathcal{H}\\}$.
∎

A general analytic reconciled probability is derived in Appendix A.5. It is
worth noting that the bootstrap reconciliation induces a tradeoff between
reduced inference speed and the requirement for knowledge of the reconciled
parameters $\tilde{\boldsymbol{\theta}}$.

###  A.5 Analytical Reconciled Probabilities

We use the bootstrap sample reconciliation technique (Panagiotelis et al.,
2023) to ensure the probabilistic coherence of HINT. This technique restores
aggregation constraints to base samples, regardless of their distribution. It
enhances HINT’s modularity by ensuring its probabilistic coherence on a wide
range of base probabilities, including non-parametric ones, without requiring
any modifications to the original algorithm. We show how a reconciled
probability can be recovered analytically through change of variables and
marginalization. For simplicity the proofs refer to $\mathbf{y}_{[b]}$ as
$\mathbf{b}$, and $\mathbf{y}_{[a]}$ as $\mathbf{a}$.

Lemma. Consider the classic reconciliation approach where the entire
hierarchy’s forecasts are combined into reconciled bottom-level forecasts
using a composition of linear transformations
$\mathbf{SP}(\cdot)=\mathbf{S_{[i][b]}P_{[b][i]}}(\cdot)$. The reconciled
probability for the new bottom-level series is given by:

$\tilde{\mathbb{P}}_{[b]}\left(\mathbf{b}\right)=|\mathbf{P}^{*}|\int\hat{\mathbb{P}}_{[i]}\left(\mathbf{P}_{\perp}\mathbf{a}+\mathbf{P}^{-}\mathbf{b}\right)\delta\mathbf{a}$
(20)

where $\hat{\mathbb{P}}\left(\cdot\right)$ is the unconstrained base forecast
distribution, $\mathbf{P}_{\perp}\in\mathbb{R}^{N_{i}\times N_{a}}$,
$\mathbf{P}^{-}\in\mathbb{R}^{N_{i}\times N_{b}}$ are the orthogonal
complement of $\mathbf{P}_{[i][b]}$ and its Moore-Penrose inverse; matrix
$\mathbf{P}^{*}=\left[\mathbf{P}_{\perp}\,|\,\mathbf{P}^{-}\right]$, and
bottom level $\mathbf{b}$ and aggregate level $\mathbf{a}$ vectors are
obtained through the following variable change:

$\hat{\mathbf{y}}_{[i]}=\mathbf{P}^{*}\begin{bmatrix}\mathbf{a}\\\
\mathbf{b}\end{bmatrix}$ (21)

###### Proof.

Using the multivariate change of variables theorem, and properties of the
Jacobian of a linear mapping:

$\tilde{\mathbb{P}}(\mathbf{a},\mathbf{b})=\left|\det\left[\left.{\frac{d\mathbf{P}^{*}(\mathbf{z})}{d\mathbf{z}}}\right|_{\mathbf{z}=(\mathbf{a},\mathbf{b})}\right]\right|\mathbb{P}{\Bigl{(}}\mathbf{P}^{*}\begin{bmatrix}\mathbf{a}\mathbf{b}\end{bmatrix}{\Bigr{)}}=|\mathbf{P}^{*}|\mathbb{P}_{[i]}\left(\mathbf{P}_{\perp}\mathbf{a}+\mathbf{P}^{-}\mathbf{b}\right)\qquad\qquad\qquad\qquad$
(22)

Marginalizing $\mathbf{a}$ we obtain the reconciled probability for the bottom
level series. ∎

Theorem 3.1: Consider the classic reconciliation approach where the entire
hierarchy’s forecasts are combined into reconciled bottom-level forecasts
using a composition of linear transformations
$\mathbf{SP}(\cdot)=\mathbf{S_{[i][b]}P_{[b][i]}}(\cdot)$. The reconciled
probability for the entire hierarchical series is given by:

$\tilde{\mathbb{P}}\left(\mathbf{y}_{[i]}\right)=|\mathbf{S}^{*}|\tilde{\mathbb{P}}_{[b]}\left(\mathbf{S}^{-}\mathbf{y}_{[i]}\right)\mathbbm{1}\\{\mathbf{y}_{[i]}\in\mathcal{H}\\}$
(23)

where $\tilde{\mathbb{P}}_{[b]}(\cdot)$ is the reconciled bottom forecast
distribution, $\mathbbm{1}(\tilde{\mathbf{y}}_{[i]}\in\mathcal{H})$ indicates
if realization belongs in the $N_{b}$-dimensional hierarchically coherent
subspace $\mathcal{H}$, $\mathbf{S}^{-}\in\mathbb{R}^{N_{b}\times N_{i}}$ is
$\mathbf{S}_{[i][b]}$ Moore-Penrose inverse and
$\mathbf{S}_{\perp}\in\mathbb{R}^{N_{i}\times N_{a}}$ its orthogonal
complement.

###### Proof.

This proof follows closely that provided in Panagiotelis et al. (2023). Given
bottom level forecast distribution from the Lemma, one can create a degenerate
distribution for the entire hierarchy by adding additional dimensions
$\mathbf{u}\in\mathbb{R}^{N_{a}}$.

$\tilde{\mathbb{P}}_{[i]}(\mathbf{u},\mathbf{b})=\tilde{\mathbb{P}}_{[b]}\left(\mathbf{b}\right)\mathbbm{1}\\{\mathbf{u}=\mathbf{0}\\}$
(24)

Let $\mathbf{S}=\mathbf{S}_{[i][b]}$ and $\mathbf{S}^{\intercal}_{\perp}$ its
orthogonal complement, and $\mathbf{S}^{-}$ and $\mathbf{S}^{-}_{\perp}$ the
respective Moore-Penrose inverses. Using the following change of variables
$\mathbf{b}=\mathbf{S}^{-}\mathbf{y}_{[i]}$ and
$\mathbf{u}=\mathbf{S}^{\intercal}_{\perp}\mathbf{y}_{[i]}$ we obtain

$\mathbf{y}_{[i]}=\left[\mathbf{S}^{-}_{\perp}\,|\,\mathbf{S}\right]\begin{bmatrix}\mathbf{u}\\\
\mathbf{b}\end{bmatrix}\quad\iff\quad\mathbf{S}^{*}\mathbf{y}_{[i]}=\begin{bmatrix}\mathbf{S}^{\intercal}_{\perp}\\\
\mathbf{S}^{-}\end{bmatrix}\mathbf{y}_{[i]}=\begin{bmatrix}\mathbf{u}\\\
\mathbf{b}\end{bmatrix}$ (25)

$\tilde{\mathbb{P}}\left(\mathbf{y}_{[i]}\right)=|\mathbf{S}^{*}|\tilde{\mathbb{P}}_{[b]}\left(\mathbf{S}^{\intercal}_{\perp}\mathbf{0}+\mathbf{S}^{-}\mathbf{y}_{[i]}\right)\mathbbm{1}\\{\mathbf{S}^{\intercal}_{\perp}\mathbf{y}_{[i]}=\mathbf{0}\\}=|\mathbf{S}^{*}|\tilde{\mathbb{P}}_{[b]}\left(\mathbf{S}^{-}\mathbf{y}_{[i]}\right)\mathbbm{1}\\{\mathbf{y}_{[i]}\in\mathcal{H}\\}$
(26)

By definition of the orthogonal complement if
$\mathbf{S}^{\intercal}_{\perp}\mathbf{y}_{[i]}=\mathbf{0}$, that means that
$\mathbf{y}_{[i]}\in\mathrm{span}(\mathbf{S})$, that matches the definition of
the hierarchically coherent subspace $\mathcal{H}$.

∎

As mentioned earlier, the analytical version of reconciled probabilities can
provide highly efficient inference times, depending on the properties of the
reconciliation and the base forecast distributions. For instance, recent
studies have utilized Gaussian distributions (Panagiotelis et al., 2023;
Wickramasuriya, 2023), and Poisson Mixtures (Olivares et al., 2022b).

###  A.6 HINT Parameter Estimation

Maximum Likelihood Estimation

To estimate the HINT’s parameters, we can use maximum likelihood estimation
for the multivariate probability in Equation (7). Let $\boldsymbol{\omega}$
denote HINT’s parameters that condition the probabilistic output layer
parameters. Then, we express the negative log-likelihood function as follows:

$\mathcal{L}(\boldsymbol{\omega})=-\mathrm{log}\Bigg{[}\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}(\boldsymbol{\omega})\prod_{(\iota,\tau)\in[i][t+1:t+h]}\left(\frac{1}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})\sqrt{2\pi}}\exp{\bigl{\\{}-\frac{1}{2}\left(\frac{y_{\iota,\kappa,\tau}-\hat{\mu}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}\right)^{2}\bigl{\\}}}\right)\Bigg{]}$
(27)

Although standard maximum likelihood estimation can model relationships
between multiple time series across the forecast horizon, a scalability
challenge arises when the number of series and forecast horizon increase
significantly. Since its computation requires access to the entire
multivariate series, MLE can become computationally intensive and time-
consuming.

Maximum Composite Likelihood Estimation

Composite likelihood provides a computationally efficient alternative to
maximum likelihood estimation for optimizing the parameters of HINT. Unlike
MLE, which computes the whole multivariate likelihood, composite likelihood
decomposes the hierarchical variable high-dimensional space support into sub-
spaces and optimizes the weighted product of the subspaces’ marginal
likelihood. When defining the sub-spaces in composite likelihood, the
probabilistic model is restricted to learning relationships within each sub-
space while assuming independence across non-overlapping sub-spaces. These
sub-spaces can be defined based on the user’s application needs. For instance,
they can be guided by the geographic proximity of the time series data. In
order to simplify the HINT algorithm, we randomly assign each series to the
sub-spaces defined by the stochastic gradient batches. Let
$\mathcal{B}=\\{[b_{i}]\\}$ be time-series SGD batches, then HINTś negative
log composite likelihood is:

$\mathcal{L}(\boldsymbol{\omega})=-\sum_{[b_{i}]\in\mathcal{B}}\mathrm{log}\left[\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}(\boldsymbol{\omega})\prod_{(\iota,\tau)\in[b_{i}][t+1:t+h]}\left(\frac{1}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})\sqrt{2\pi}}\exp{\bigl{\\{}-\frac{1}{2}\left(\frac{y_{\iota,\kappa,\tau}-\hat{\mu}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}\right)^{2}\bigl{\\}}}\right)\right]$
(28)

Composite likelihood’s independent sub-spaces can be defined as each series
forecast, leading to the univariate estimation approach.

$\mathcal{L}_{univ}(\boldsymbol{\omega})=-\sum_{\iota\in[i]}\mathrm{log}\left[\sum_{\kappa=1}^{N_{k}}\hat{w}_{\kappa}(\boldsymbol{\omega})\prod_{\tau\in[t+1:t+h]}\left(\frac{1}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})\sqrt{2\pi}}\exp{\bigl{\\{}-\frac{1}{2}\left(\frac{y_{\iota,\kappa,\tau}-\hat{\mu}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}{\hat{\sigma}_{\iota,\kappa,\tau}(\boldsymbol{\omega})}\right)^{2}\bigl{\\}}}\right)\right]$
(29)

Figure 5 compares univariate and composite likelihood estimation.

##  Appendix B Hierarchical Dataset’s Exploration

In this Appendix we complement the description of the benchmark datasets from
Section 4.

Table A1: Summary, of experiment hierarchical datasets. Dataset | Total | Aggregate | Bottom | Frequency. | h | Levels
---|---|---|---|---|---|---
Labour | 57 | 25 | 32 | Month | 8 | 4
Traffic | 207 | 7 | 200 | Daily | 7 | 3
Tourism | 89 | 33 | 56 | Quarterly | 4 | 4
Tourism-L | 555 | 175 | $76\,/\,304$ | Month | 12 | 4/5
Wiki2 | 199 | 49 | 150 | Daily | 7 | 5

Labour reports monthly Australian employment from February 1978 to December
2020. It contains a structure built by the labour categories (Australian
Bureau of Statistics, 2019). Traffic measures the occupancy of 963 traffic
lanes in the Bay Area, the data is grouped into a year of daily observations
and organized into a 207 hierarchical structure (Dua and Graff, 2017). Tourism
consists of 89 Australian location quarterly visits series; it covers from
1998 to 2006. Several studies have used this dataset in the past (Tourism
Australia, Canberra, 2005). Tourism-L summarizes an Australian visitor survey
managed by the Tourism Research Australia, the dataset contains 555 monthly
series from 1998 to 2016, and it is organized into geographic and purpose of
travel (Tourism Australia, Canberra, 2019). Wiki2 contains the daily views of
145,000 Wikipedia articles from July 2015 to December 2016. The dataset is
filtered and processed into 150 bottom series and 49 aggregate series (Anava
et al., 2018). Figure 3 shows each dataset’s most aggregated series along with
its training methodology partition. Figure 2 shows each dataset’s hierarchical
aggregation constraint matrices.

(a) Labour

(b) Traffic

(c) TourismS

(d) TourismL

(e) Wiki2

Figure 2: Dataset’s hierarchical constraints. (a) Labour groups 32 occupation
series by gender and geography. (b) Traffic groups 200 highways’ occupancy
series into quarters, halves and total. (c) Tourism groups 56 quarterly
Australian tourist visits by geographic levels. (d) Tourism-L groups its 555
monthly Australian regional visit series, into a combination travel purpose,
zones, states and country geographical aggregations. (e) Wiki2 groups 150
daily visits to Wikipedia articles by language and article categorical
taxonomy. Figure 3: Datasets’ partition into train, validation, and test sets
used in our experiments. All use the last horizon window as defined in Table
A1 (marked by the second dotted line), and the previous window preceding the
test set as validation (between the first and second dotted lines). Validation
provides the signal for hyperparameter optimization.

##  Appendix C Ablation Studies

In this Appendix, we perform ablation studies on the validation set of five
hierarchical datasets Labour, Traffic, Tourism, Tourism-L, and Wiki2. For
these experiments, we change minimally the HINT settings defined in Table A4,
comparing normalization techniques, varying the number of Mixture components,
and exploring different hierarchical reconciliation strategies to understand
their contribution to the performance of the method.

###  C.1 Scaled Decoupled Optimization

Section 3.2, introduced HINT’s scale decouple optimization strategy with the
TemporalNorm transform. Here we study the effects of different temporal
normalization strategies on the forecast accuracy performance of the model,
measured with the overall sCRPS. For simplicity consider a network with only
temporal input $\mathbf{x}^{h}_{[i][:t][c]}$, with $[i]$ batch, $[:t]$ time,
and $[c]$ feature channel indexes, we consider transformations of the general
form:

$\begin{split}\check{\mathbf{x}}^{h}_{[i][:t+h][c]}&=\mathrm{TemporalNorm}(\mathbf{x}^{h}_{[i][:t+h][c]})=\frac{\mathbf{x}^{h}_{[i][:t][c]}-\mathbf{a}}{\mathbf{b}}\\\
\hat{\theta}(\check{\mathbf{x}}^{h}_{[i][:t+h][c]})&=\mathrm{TemporalNorm}^{-1}(\boldsymbol{\omega}(\check{\mathbf{x}}_{[i][:t+h][c]}))=\mathbf{b}\boldsymbol{\omega}_{[i][t+h]}+\mathbf{a}\end{split}$
(30)

where $\mathbf{a},\mathbf{b}\in\mathbb{R}^{(N_{a}+N_{b})\times N_{c}}$ is the
shift and the scale of the historic inputs. In this experiment we augment
NHITS, Autoformer, LSTM, MLP, and TFT with three different TemporalNorm
normalization schemes, along with the recently proposed Reversible Instance
Normalization (RevIN) technique:

$\displaystyle\begin{aligned}
\mathrm{minmax}&:\frac{(\mathbf{x}^{h}_{[i][:t][c]}-\mathrm{min}(\mathbf{x}^{h}_{[i][:t][c]})_{[i][c]})}{(\mathrm{max}(\mathbf{x}^{h}_{[i][:t][c]})_{[i][c]}-\mathrm{min}(\mathbf{x}^{h}_{[i][:t][c]})_{[i][c]})}&\quad\quad\mathrm{standard}&:\frac{(\mathbf{x}^{h}_{[i][:t][c]}-\bar{\mathbf{X}}_{[i][c]})}{\hat{\sigma}_{[i][c]}}\\\
&&\\\
\mathrm{robust}&:\frac{(\mathbf{x}^{h}_{[i][:t][c]}-\mathrm{median}(\mathbf{x}^{h}_{[i][:t][c]})}{\mathrm{mad}(\mathbf{x}^{h}_{[i][:t][c]})_{[i][c]}}&\quad\quad\mathrm{revin}&:\lambda_{[i][c]}*\frac{(\mathbf{x}^{h}_{[i][:t][c]}-\bar{\mathbf{X}}_{[i][c]})}{\hat{\sigma}_{[i][c]}}+\beta_{[i][c]}\\\
\end{aligned}$

Figure 4 compares the four different normalization techniques across different
architectures for different levels of noise, on the Tourism-L dataset. It is
evident that the performance improvements vary depending on the type of
temporal normalization employed. In cases where signals are smooth, the type
of normalization used makes little difference. However, as the signal becomes
noisier and more volatile, the robust normalization strategy consistently
outperforms all other normalization techniques, with the minmax strategy
lagging behind. This underscores the advantages of employing median shift and
mad scale to enhance the model’s robustness against large outliers in the
data. As a result, we have selected robust Temporal Normalization as the
default option for all subsequent experiments.

|  |
---|---|---
|  |

(a) Autoformer (b) LSTM

(c) MLP

(d) NHITS (e) TFT

Figure 4: Validation _scaled Continuous Ranked Probability Score_ (sCRPS)
curves on Tourism-L for five different architectures and varying noise levels.
For each noise level _p_ , we transformed a fraction _p_ of the training
dataset by a random scaling factor. We see the robust normalization (green
curve) is the lowest for all five architectures, showing how it consistently
outperforms the other normalization techniques at larger noise levels.

###  C.2 Mixture Size Exploration

As mentioned in Section 3.1 and proven in Appendix A, our multivariate mixture
probability model is capable of capturing the relationships among the
hierarchical series and its expressivity directly determined by the number of
components in the Gaussian mixture. Here we conduct a study to compare the
accuracy effects of different mixture sizes. We vary the number of components
monotonically and follow the model’s performance. An NHITS model configuration
is automatically selected as defined from Table A4.

During our Labour, Traffic, and Tourism-L experiments, we observed a bias-
variance trade-off relationship between the number of mixture components and
the validation sCRPS. If we set the number of components to 1, the sCRPS score
is the highest, but as the number of components increases, the sCRPS improves.
However, if there are too many components, the CRPS score worsens again after
a certain point. Our theory is that if the mixture has too few components, the
probability does not have enough parameters to accurately depict the data’s
correlations; this leads to a high bias and reflects in a poor sCRPS. On the
other hand, if the mixture has too many components, the model becomes too
complicated and quickly overfits the training data, resulting in high variance
and poor performance on the validation data. The results of this ablation
experiment show a clear benefit of using a flexible multivariate mixture
distribution, contrasting it to the simpler approach of using a single
component (Gaussian/Poisson regression). We explain these improvements from
the flexibility of the mixture approach that operates as a Kernel density
estimation and can arbitrarily approximate a wide variety of target
distributions.

Figure 5:  Validation _scaled Continuous Ranked Probability Score_ (sCRPS) for
Labour, Traffic, and Tourism-L. We show performance curves for NHITS as a
function of the number of Mixture Components. We observe a bias-variance
tradeoff, where initially, the sCRPS decreases as the number of components
increases and reaches an optimal value at K=10 components. From there, after
that, we see the sCRPS worsening, thus giving us the classic U-shaped tradeoff
pattern.

###  C.3 Architectures and Reconciliation Ablation Study

As mentioned in Section 3.1 and proven in Appendix A, the flexibility of the
multivariate mixture probability model is compatible with most probabilistic
reconciliation techniques. In this ablation study, we compare the accuracy
effect of different reconciliation strategies. The reconciliation strategies
considered are BottomUp (Orcutt et al., 1968; Dunn et al., 1976), TopDown
(Gross and Sohl, 1990; Fliedner, 1999), and MinTrace (Hyndman et al., 2011;
Wickramasuriya et al., 2019) variants. We describe them in detail below.

Consider the base forecasts
$\hat{\mathbf{y}}_{[i],\tau}\in\mathbb{R}^{N_{a}+N_{b}}$, a reconciliation
process uses a matrix
$\mathbf{P}_{[b][i]}\in\mathbb{R}^{N_{b}\times(N_{a}+N_{b})}$ that collapses
the original base forecasts into bottom-level forecast that are later
aggregated for the upper levels of the hierarchy into the reconciled forecasts
$\tilde{\mathbf{y}}_{[i],\tau}$. Here we use the convenient representation of
the reconciliation strategies introduced in Section 2.

$\tilde{\mathbf{y}}_{[i],\tau}=\mathbf{S}_{[i][b]}\mathbf{P}_{[b][i]}\hat{\mathbf{y}}_{[i],\tau}=\mathbf{SP}\left(\hat{\mathbf{y}}_{[i],\tau}\right)$

Bottom-Up. The most basic hierarchical reconciliation consists of simply
aggregating the bottom-level base forecasts $\hat{\mathbf{y}}_{[b],\tau}$. By
construction, it satisfies the hierarchical aggregation constraints.

$\mathbf{P}_{[b][i]}=[\mathbf{0}_{\mathrm{[b][a]}}\;|\;\mathbf{I}_{\mathrm{[b][b]}}]$
(31)

Top-Down. The TopDown strategy distributes an aggregate level forecast into
the bottom-level forecasts using proportions $\mathbf{p}_{[b]}$. Proportions
can be historical values, or they can be forecasted.

$\mathbf{P}_{[b][i]}=[\mathbf{p}_{\mathrm{[b]}}\;|\;\mathbf{0}_{\mathrm{[b][a,b\;-1]}}]$
(32)

MinTrace. Newer reconciliation strategies use all the information available
throughout the hierarchy optimally. In particular, the MinTrace reconciliation
is proven to be the optimizer of a mean squares error objective that
transforms base predictions into hierarchically coherent predictions under an
unbiasedness assumption. Its reconciliation matrix is given by:

$\mathbf{P}_{[b][i]}=\left(\mathbf{S}^{\intercal}\hat{\boldsymbol{\Sigma}}_{\tau}\mathbf{S}\right)^{-1}\mathbf{S}^{\intercal}\hat{\boldsymbol{\Sigma}}^{-1}_{\tau}$
(33)

We summarize the ablation study results for the different reconciliation
strategies in Table A2. We report the overall sCRPS across five datasets for
the NHITS, TCN, and TFT architectures. We obtain the probabilistic predictions
using bootstrap (Panagiotelis et al., 2023). We observe clear advantages from
adopting novel reconciliation techniques such as MinTrace as it improves
accuracy over BottomUp by 20 to 30 percent margins across well-established
neural forecast architectures. We find that post processing reconciliation is
capable of improving complex end-to-end approaches that integrate the
hierarchical constraints into the training procedure Rangapuram et al. 2021;
Kamarthi et al. 2022; Han et al. 2021. Based on the results of these ablation
studies, we conducted the main experiments of this work with the MinTrace and
the BottomUp reconciliation techniques.

Table A2:  Empirical evaluation of probabilistic coherent forecasts. Mean
scaled continuous ranked probability score (sCRPS), averaged over 10 random
seeds, at each aggregation level. The best result is highlighted (lower
measurements are preferred).
* The TopDown reconciliation is only available for strictly hierarchical datasets.  | Dataset | Base |  MinTrace-ols |  MinTrace-wls |  TopDown-ap* |  TopDown-pa* | BottomUp
---|---|---|---|---|---|---|---
NHITS | Labour | 0.0082 | 0.0082±0.0001 | 0.0084±0.0001 | 0.0092±0.0001 | 0.0091±0.0000 | 0.0094±0.0001
Traffic | 0.0629 | 0.0635±0.0011 | 0.0643±0.0010 | 0.0651±0.0010 | 0.0650±0.0013 | 0.0660±0.0008
Tourism | 0.0791 | 0.0806±0.0011 | 0.0771±0.0014 | 0.0920±0.0010 | 0.0913±0.0014 | 0.0756±0.0012
Tourism-L | 0.1274 | 0.1281±0.0004 | 0.1261±0.0006 | - | - | 0.1351±0.0005
Wiki2 | 1.4531 | 1.3165±0.0302 | 1.8399±0.0904 | 0.5165±0.0159 | 0.5178±0.0088 | 3.3351±0.1690
TCN | Labour | 0.0213 | 0.0243±0.0004 | 0.0202±0.0003 | 0.0237±0.0003 | 0.0237±0.0004 | 0.0187±0.0004
Traffic | 0.0566 | 0.0569±0.0007 | 0.0577±0.0007 | 0.0605±0.0010 | 0.0605±0.0009 | 0.0623±0.0008
Tourism | 0.0664 | 0.0660±0.0009 | 0.0641±0.0011 | 0.0803±0.0012 | 0.0807±0.0015 | 0.0678±0.0018
Tourism-L | 0.1632 | 0.1638±0.0010 | 0.1640±0.0008 | - | - | 0.1677±0.0007
Wiki2 | 2.7345 | 2.0942±0.0539 | 2.9305±0.1247 | 1.6188±0.0374 | 1.6232±0.0476 | 4.2341±0.1338
TFT | Labour | 0.0073 | 0.0071±0.0001 | 0.0074±0.0000 | 0.0084±0.0001 | 0.0084±0.0001 | 0.0087±0.0001
Traffic | 0.0632 | 0.0641±0.0007 | 0.0638±0.0008 | 0.0650±0.0015 | 0.0646±0.0011 | 0.0658±0.0009
Tourism | 0.0944 | 0.1015±0.0010 | 0.0895±0.0007 | 0.0834±0.0008 | 0.0832±0.0006 | 0.0922±0.0013
Tourism-L | 0.1360 | 0.1364±0.0007 | 0.1362±0.0008 | - | - | 0.1442±0.0010
Wiki2 | 0.2609 | 0.2608+0.0010 | 0.2641±0.0025 | 0.2560±0.0016 | 0.2560±0.0013 | 0.2755±0.0041

##  Appendix D Software and Training Methodology

###  D.1 Hyperaparameters and Training Methodology

Table A3: HINT fixed hyperparameters. Hyperparameter | Fixed Values
---|---
Architecture | NHITS | TFT | TCN
Activation | ReLU | ReLU | ReLU
Encoder units | 256 | 256 | 256
Encoder layers* | 4 | 3 | 4
Encoder type | MLP | LSTM | Conv1D
Train Objective | Comp.Lik. | Comp.Lik. | Comp.Lik.
Table A4: HINT optimized hyperparameters. Hyperparameter | Considered Values
---|---
Initial learning rate. | {1e-3,5e-4,1e-4}
Number of learning rate decays. | {None, 3}
Training steps. | {.5e3, 1e3, 1.5e3, 2e3, 2.5e3, 3e3}
Input size multiplier (L=m*H). | $m\in\\{2,3,4\\}$
Reconciliation strategy. | {BottomUp,
| MinTraceOLS, MinTraceWLS }

Training HINT and the benchmark models involves dividing the data into
training, validation, early stopping, and test sets, as shown in Figure 3. The
training set consists of the observations before the last two horizon windows;
validation is the window between the train and test sets, with test being the
last window. The model’s performance on the validation set guides the
exploration of the hyperparameter space (HYPEROPT, Bergstra et al. 2011).
During the recalibration phase, we retrain the models to incorporate new
information before being tested.

We followed a standard two-stage approach for hyperparameter selection. In the
first stage, based on validation ablation studies from Appendix C.3, we fixed
the architecture and the probability distribution to be estimated; Table A3
describes the hyperparameters. Then conditional on a good performing
architecture, in the second stage, we optimized the training procedure of the
architecture, optimally exploring the space defined in Table A4 with HYPEROPT.
This approach allowed us to explore the hyperparameter space while keeping it
computationally tractable. It also demonstrated the HINT’s robustness, broad
applicability, and potential to achieve high accuracy with only slight
adjustments.

We train HINT to maximize the composite likelihood from Equation (28) using
the ADAM (Kingma and Ba, 2014) stochastic gradient algorithm. An early
stopping strategy (Yao et al., 2007) is employed to halt training if there is
no improvement in the overall sCRPS metric on the validation set.

###  D.2 Software Implementation

All statistical baselines use StatsForecast’s AutoARIMA (Hyndman and
Khandakar, 2008; Garza et al., 2022) and HierarchicalForecast’s reconciliation
methods implementations (Olivares et al., 2022b). We created a unified Python
re implementation of various widely-used hierarchical forecasting techniques
that we make publicly available in the HierarchicalForecast library (Olivares
et al., 2022b). The shared implementation allows us to standardize the
comparison of the methods, controlling for experimental details and
guaranteeing the quality of statistical baselines. The code is publicly
available in a dedicated repository to support reproducibility and related
research.

Regarding the hierarchical neural forecast baselines, HierE2E (Rangapuram et
al., 2021) is available in the GluonTS library, while PROFHIT (Kamarthi et
al., 2022) is available in a PROFHIT dedicated repository. As mentioned
earlier the only available implementation for PROFHIT suffers from significant
numerical instability in its optimization. We use the optimal configurations
reported in HierE2E and PROFHIT repositories.

The HINT framework is implemented in PyTorch (Paszke et al., 2019) and can be
run on both CPUs and GPUs. We have made the HINT source code available, along
with all the experiments in the following HINT dedicated repository.

##  Appendix E Extended Main Results

Table 1 summarizes probabilistic and point forecast accuracy. We provide
measurements of the best results between NHITS, TFT and a simple TCN for
different hierarchy levels, following the training and hyperparameter
selection methodologies in Appendix D. The top row of each panel reports the
overall sCRPS or relMSE. Moving from aggregate to disaggregate levels
increases forecast errors. HINT outperforms the second-best alternative
overall sCRPS by an average of 8.14% across datasets, and 13.2% without
Traffic. We observed improvements of 8.2% on Labour, 17.4% on Tourism, 10.9%
on Tourism-L, and 15.5% on Wiki2. However, for Traffic, HierE2E outperforms
HINT by -11.3% due to the clear Granger causalities in the dataset, which
merit the HierE2E’s VAR approach. The sCRPS and relMSE results are highly
correlated, although relMSE is less robust to outliers.

Table A5:  Empirical evaluation of probabilistic coherent forecasts. Mean scaled continuous ranked probability score (sCRPS) and mean relative squared error (relMSE), averaged over 10 random seeds, at each aggregation level. The best result is highlighted. † The PROFHIT results differ from Kamarthi et al. (2022), as the only available implementation suffers from significant numerical instability in its optimization. * Best performing variant of TopDown (avg. proportions, proportions avg.), and MinTrace (ols, wls, shrinkage) reported. ** The PERMBU/TopDown only available for strictly hierarchical datasets.  |  |  | HINT | OTHER | BOOTSTRAP | PERMBU**
---|---|---|---|---|---|---
| Dataset | Level | (Ours) | HierE2E | PROFHIT† | BottomUp | TopDown* | MinTrace* | BottomUp | TopDown* | MinTrace*
sCRPS |  Labour | Overall | .0067±.0000 | .0171±.0003 | .2138±.0007 | .0078±.0001 | .0668±.0000 | .0073±.0000 | .0077±.0001 | .0623±.0001 | .0069±.0001
Country | .0017±.0000 | .0052±.0003 | .2097±.0038 | .0021±.0001 | .0012±.0001 | .0013±.0001 | .0026±.0001 | .0014±.0001 | .0016±.0001
Region | .0047±.0001 | .0181±.0003 | .2150±.0028 | .0058±.0001 | .0458±.0000 | .0045±.0001 | .0060±.0001 | .0420±.0001 | .0045±.0001
Region/Gender | .0075±.0000 | .0188±.0003 | .2161±.0012 | .0088±.0001 | .0798±.0001 | .0087±.0001 | .0084±.0001 | .0739±.0001 | .0078±.0001
Reg/Gndr/Empl | .0128±.0000 | .0262±.0004 | .2142±.0027 | .0145±.0001 | .1403±.0000 | .0148±.0001 | .0137±.0001 | .1320±.0002 | .0137±.0001
Traffic | Overall | .0589±.0004 | .0426±.0008 | .1137±.0022 | .0736±.0024 | .0741±.0012 | .0608±.0014 | .0849±.0009 | .0708±.0008 | .0651±.0008
Level1 | .0340±.0008 | .0276±.0011 | .0899±.0059 | .0468±.0031 | .0301±.0020 | .0299±.0020 | .0651±.0012 | .0373±.0011 | .0367±.0011
Level2 | .0347±.0006 | .0287±.0009 | .0879±.0034 | .0483±.0030 | .0329±.0017 | .0323±.0017 | .0622±.0013 | .0367±.0009 | .0357±.0008
Level3 | .0392±.0005 | .0297±.0009 | .0926±.0032 | .0530±.0025 | .0360±.0013 | .0385±.0014 | .0614±.0010 | .0383±.0009 | .0405±.0010
Level4 | .1275±.0002 | .0845±.0003 | .1842±.0014 | .1463±.0017 | .1975±.0017 | .1424±.0015 | .1507±.0004 | .1709±.0010 | .1473±.0004
Tourism | Overall | .0536±.0004 | .0761±.0007 | .1358±.0033 | .0682±.0018 | .1040±.0014 | .0703±.0017 | .0649±.0016 | .0898±.0012 | .0680±.0016
Country | .0147±.0004 | .0400±.0009 | .0941±.0151 | .0290±.0028 | .0333±.0025 | .0335±.0026 | .0267±.0023 | .0329±.0021 | .0333±.0025
Purpose | .0360±.0004 | .0609±.0012 | .1300±.0069 | .0490±.0027 | .0782±.0017 | .0507±.0023 | .0450±.0017 | .0697±.0021 | .0497±.0018
State/Purpose | .0709±.0006 | .0914±.0008 | .1323±.0076 | .0828±.0016 | .1399±.0010 | .0845±.0016 | .0793±.0014 | .1176±.0013 | .0806±.0014
Region/Purpose | .0929±.0006 | .1122±.0007 | .1867±.0031 | .1118±.0012 | .1646±.0010 | .1124±.0013 | .1087±.0017 | .1390±.0014 | .1085±.0016
Tourism-L | Overall | .1176±.0002 | .1424±.0019 | .2139±.0014 | .1375±.0013 | - | .1313±.0009 | - | - | -
Country | .0325±.0006 | .0698±.0029 | .1353±.0090 | .0622±.0026 | - | .0471±.0018 | - | - | -
State | .0606±.0006 | .0936±.0019 | .1610±.0020 | .0820±.0019 | - | .0723±.0011 | - | - | -
Zone | .1025±.0004 | .1260±.0017 | .1893±.0034 | .1207±.0010 | - | .1143±.0007 | - | - | -
Region | .1457±.0003 | .1653±.0016 | .2277±.0022 | .1646±.0007 | - | .1591±.0006 | - | - | -
Purpose | .0706±.0006 | .0996±.0028 | .1845±.0071 | .0788±.0018 | - | .0723±.0014 | - | - | -
State/Purpose | .1088±.0003 | .1317±.0021 | .2160±.0031 | .1268±.0017 | - | .1243±.0014 | - | - | -
Zone/Purpose | .1772±.0003 | .1926±.0015 | .2679±.0019 | .1949±.0010 | - | .1919±.0008 | - | - | -
Region/Purpose | .2426±.0005 | .2606±.0017 | .3296±.0010 | .2698±.0008 | - | .2694±.0006 | - | - | -
Wiki2 | Overall | .2447±.0007 | .2592±.0031 | .4009±.0028 | .2894±.0038 | .3231±.0037 | .2808±.0035 | .3920±.0044 | .4269±.0036 | .3821±.0049
World | .1247±.0016 | .1007±.0046 | .1244±.0085 | .1796±.0069 | .1777±.0084 | .1793±.0067 | .1777±.0125 | .1945±.0109 | .1801±.0123
Country | .1805±.0011 | .1963±.0037 | .2775±.0141 | .2392±.0047 | .2437±.0058 | .2232±.0043 | .2778±.0073 | .3036±.0029 | .2684±.0066
Access | .2546±.0010 | .2784±.0038 | .4405±.0034 | .2966±.0032 | .3379±.0026 | .2781±.0028 | .4196±.0059 | .4621±.0071 | .4006±.0059
Agent | .2699±.0010 | .2900±.0043 | .4526±.0084 | .3036±.0033 | .3427±.0026 | .2855±.0029 | .4255±.0058 | .4669±.0070 | .4073±.0060
Topic | .3938±.0016 | .4307±.0039 | .7094±.0109 | .4282±.0038 | .5134±.0037 | .4379±.0029 | .6595±.0060 | .7074±.0049 | .6540±.0058
relMSE |  Labour | Overall | .5802±.0131 | .5667±.0265 | $6.774\times 10^{3}$ | .5382±.0000 | 16.8204±.0000 | .3547±.0000 |  |  |
Country | .1032±.0130 | .1536±.0284 | $9.424\times 10^{2}$ | .2362±.0000 | .0542±.0000 | .0729±.0000 |  |  |
Region | .6502±.0276 | 1.1486±.0413 | $6.635\times 10^{2}$ | .8281±.0000 | 14.6118±.0000 | .3740±.0000 |  |  |
Region/Gender | 1.6870±.0241 | 1.1206±.0395 | $4.113\times 10^{2}$ | .9021±.0000 | 35.6038±.0000 | .7519±.0000 |  |  |
Reg/Gndr/Empl | 1.8103±.0131 | 1.4491±.0361 | $1.664\times 10^{2}$ | .8069±.0000 | 5.6047±.0000 | .8041±.0000 |  |  |
Traffic | Overall | .1212±.0051 | .0340±.0051 | .4536±.0224 | .1394±.0000 | .0614±.0000 | .0744±.0000 |  |  |
Level1 | .1004±.0057 | .0253±.0057 | .4591±.0413 | .1296±.0000 | .0491±.0000 | .0634±.0000 |  |  |
Level2 | .1156±.0051 | .0302±.0050 | .4291±.0336 | .1342±.0000 | .0625±.0000 | .0690±.0000 |  |  |
Level3 | .1602±.0041 | .0529±.0038 | .4612±.0240 | .1582±.0000 | .0730±.0000 | .0958±.0000 |  |  |
Level4 | .8893±.0042 | .4206±.0048 | .7709±.0081 | .6457±.0000 | .6525±.0000 | .6194±.0000 |  |  |
Tourism | Overall | .0387±.0007 | .1471±.0046 | .9745±.0803 | .1002±.0000 | .1919±.0000 | .1235±.0000 |  |  |
Country | .0200±.0013 | .1821±.0094 | 1.2240±.1474 | .0841±.0000 | .1328±.0000 | .1233±.0000 |  |  |
Purpose | .0342±.0009 | .1038±.0040 | .8208±.0487 | .0778±.0000 | .1669±.0000 | .0957±.0000 |  |  |
State/Purpose | .0750±.0020 | .1550±.0032 | .8511±.0489 | .1563±.0000 | .3482±.0000 | .1620±.0000 |  |  |
Region/Purpose | .0928±.0015 | .1772±.0027 | .7227±.0942 | .2000±.0000 | .3628±.0000 | .2007±.0000 |  |  |
Tourism-L | Overall | .0577±.0009 | .2449±.0096 | 1.0401±.0296 | .3070±.0000 | - | .1375±.0000 |  |  |
Country | .0336±.0009 | .2918±.0194 | 1.3473±.1430 | .4399±.0000 | - | .1268±.0000 |  |  |
State | .0598±.0012 | .2850±.0107 | 1.0854±.0358 | .3504±.0000 | - | .1564±.0000 |  |  |
Zone | .1263±.0019 | .3620±.0087 | 1.0821±.0397 | .3950±.0000 | - | .2664±.0000 |  |  |
Region | .1777±.0024 | .3594±.0055 | .9447±.0300 | .3996±.0000 | - | .3211±.0000 |  |  |
Purpose | .0478±.0006 | .1581±.0086 | .8257±.0758 | .1624±.0000 | - | .0759±.0000 |  |  |
State/Purpose | .0752±.0009 | .1884±.0055 | .7839±.0491 | .1860±.0000 | - | .1332±.0000 |  |  |
Zone/Purpose | .1638±.0017 | .2872±.0049 | .8178±.0187 | .2932±.0000 | - | .2550±.0000 |  |  |
Region/Purpose | .2296±.0016 | .3360±.0047 | .8013±.0215 | .3661±.0000 | - | .3464±.0000 |  |  |
Wiki2 | Overall | .1884±.0012 | .6598±.0249 | .7901±.0384 | 1.0163±.0000 | 1.4482±.0000 | 1.0068±.0000 |  |  |
World | .1955±.0037 | .2738±.0301 | .3274±.0798 | .9245±.0000 | 1.6135±.0000 | .9883±.0000 |  |  |
Country | .1648±.0015 | .7427±.0430 | .9150±.0886 | 1.0204±.0000 | 1.3529±.0000 | 1.0252±.0000 |  |  |
Access | .1921±.0012 | .9575±.0331 | 1.1538±.0429 | 1.1267±.0000 | 1.4159±.0000 | 1.0267±.0000 |  |  |
Agent | .1972±.0010 | .9384±.0381 | 1.1312±.0695 | 1.1008±.0000 | 1.3562±.0000 | 1.0047±.0000 |  |  |
Topic | .1923±.0010 | 1.0305±.0145 | 1.1881±.0492 | 1.0603±.0000 | 1.2282±.0000 | 1.0182±.0000 |  |  |