$\clubsuit$ CLOVER $\clubsuit$: Probabilistic Forecasting with Coherent Learning Objective Reparameterization
Authors: Kin G. Olivares
Geoffrey Négiar
Ruijun Ma
Oinam Nganba Meetei
Mengfei Cao
Michael W. Mahoney
Venue: Accepted by TMLR ()

Abstract: Obtaining accurate probabilistic forecasts is an operational challenge in many applications, such as energy management, climate forecasting, supply chain planning, and resource allocation.
Many of these applications present a natural hierarchical structure over the forecasted quantities; and forecasting systems that adhere to this hierarchical structure are said to be coherent. 
Furthermore, operational planning benefits from the accuracy at all levels of the aggregation hierarchy. However, building accurate and coherent forecasting systems is challenging: classic multivariate time series tools and neural network methods are still being adapted for this purpose. In this paper, we augment an MQForecaster neural network architecture with a modified multivariate Gaussian factor model that achieves coherence by construction. The factor model samples can be differentiated with respect to the model parameters, allowing optimization on arbitrary differentiable learning objectives that align with the forecasting system's goals, including quantile loss and the scaled Continuous Ranked Probability Score (CRPS). We call our method the Coherent Learning Objective Reparametrization Neural Network (CLOVER). In comparison to state-of-the-art coherent forecasting methods, 
CLOVER achieves significant improvements in scaled CRPS forecast accuracy, with average gains of 15%, as measured on six publicly-available datasets.

Published in Transactions on Machine Learning Research (12/2024)
CLOVER +: Probabilistic Forecasting with Coherent Learning Objective Reparameterization
Kin G. Olivares*
kigutic Samazon.com
Amazon Forecasting Science, New York
Geoffrey Négiar
The Forecasting Compony
Ruijun Ma² Amazon Forecasting Science, New York
rujjumma@omazzo.ccom
O. Nganba Meetei Peloton, New York
nganbo Caannbol@mmll..om com
Mengfei Cao Amazon Porecasting Science, New York
mffcao Cawaazo..com
Michael W. Mahoney
Amazon Forecasting Science, New York
Reviewed on OpenReview: https: openreview. net/ forum? id= FOXE
Abstract
Obtaining accurate probabilistic forecasts is an operational challenge in many applications, such as energy management, climate forecasting, supply chain planning, and resource allocation. Many of these applications present natural hierarchical structure over the forecasted quantitics; and forecasting systems that adhere to this hierarchical structure said to be coherent. Furthermore, operational planning benefits from the accuracy at all levels of the aggregation hierarchy. However, building accurate and coherent forecasting systems is challenging: classic multivariate time scries tools and neural network methods are still being adapted for this purpose. In this paper, we augment an MQForecaster neural network architecture with a modified multivariate Gaussian factor model that achieves coherence by construction. The factor model samples can be differentiated with respect to the model parameters, allowing optimization arbitrary differentiable learning objectives that align with the forecasting system's 's goals, including quantile loss and the scaled Continuous Ranked Probability Score (CRPS). We call our method the Coherent Learning Objective Reparamelrization Neural Network (CLOVER). In comparison to state-of-the-arl coherent forecasting methods, CLOVER achieves significant improvements scaled CRPS forecast accuracy, with average gains of 15%, as measured on six publicyy-vailable datascts.
Equal contribution.
Work completed Amazon Forecasting Science, where Gealfrey interned as Berkeley Ph.D. candidate.

---

Published in Transactions on Machine Learning Research (12/2024)
Coherent
Multivariate Series Inputs
Multivariate Output Distribution
Arbitrary Learning Objective
Method
End-to- Bad-to-End -nddEnd
PERMBU (Ben Taieb et al.
2017
x
x
x
X
Bootstrap Panagiotelis et al.
2023
XXX
X
X
X
Normality Wickramasuriya
2023
X
X
X
DPMN OOlvarres el al., 2023)
X
X
HierE2E [Rangapuram et al.
2021
/
X
CLOVER (ours)
/
Table 1: Coherent forecast methods' desirable properties.
1 Introduction
Obtaining accurate forecasts is important step for planning in complex and uncertain environments, with applications ranging from energy to supply chain management (Chen al. 2022
Wolff et al. 2024). from
transportation to climate prediction (Hong ct al., 2014 Gneiting & Katzfuss 2014 Makridakis et al. 2022
Going beyond point forecasts such as means and medians, probabilistic forecasting provides a key tool for predicting uncertain future events. This involves, for example, forecasting that there is a 90% chance of rain on a certain day or that there is a 99% chance that people will want to buy fewer than 100 items a certain store in a given week. Providing more detailed predictions of this form permits finer uncertainty quantification. This, in turn, allows planners prepare for different scenarios and allocate resources according to their anticipated likelihood and cost structure. This can lead to better resource allocation, improved decision-making and less waste.
In many forecasting applications, there exist natural hierarchies over the quantities one wants to predict, such as energy consumption at various temporal granularities (from monthly to weekly), different geographic levels (from building- to city-level to state-level), or retail demand for specific items (in a hicrarchical product taxonomy). Typically, most or all levels of the hierarchy are important: the bottom levels are key for operational short-term planning, while higher levels of aggregation provide insight into longer-term broader trends. Moreover, probabilistic forecasts are often desired to be coherent (or consistent) to ensure efficient decision-making at all levels (Hong et al. 2014 Jeon et al.) 2019). Coherence is achieved when the forecast distribution assigns zero probability to forecasts that do not satisfy the hierarchy's constraints (Ben Taieb et al. 2017; Panagiotelis et al.) 2023 Olivares et al. 2023) (see Definition 2.1). Designing an accurate model, capable of leveraging information from all hierarchical levels while enforcing coherence is a well-known and challenging task (Hyndman et al. 2011).
The hierarchical forecasting literature has been dominated by two-stage reconciliation approaches, where univariate methods are first. fitted and later reconciled towards coherence. For many years, most research has
focused on mean reconciliation (Hyndman et al. 2011 2024 Vitullo 2011
Hyndman et al. 2016 Dangerfield
& Morris
1992
Wickramasuriya et al. 2019
Mishchenko et al., 2019). More recent statistical methods
consider coherent probabilistic forecasts through variants of the bootstrap reconciliation technique (Ben
Taieb et al.] 2017 Panagiotelis al.) 2023) or the clever use of the properties of the Gaussian forecast.
distributions Wickramasuriya, 2023). A detailed hierarchical forecast review is provided by AAthaassoooulos
et al. (2024). Large-scale applications of hierarchical forecasting require practitioners to simplify the two-stage reconciliation process by favoring end-to-end approaches that simultaneously fit all levels of the hierarchy, while still achieving coherence. The end-to-end approach refers to training model constrained to achieve coherence by directly optimizing for accuracy. End-to-end methods offer advantages such as reduced complexity, improved computational efficiency, and improved adaptability by streamlining the entire forecasting pipeline into a single, unified model. More importantly, end-to-cnd models generally achieve better accuracy compared to two-stage models that are first trained independently for optimized accuracy and then made coherent. through various reconciliation approaches (Rangapuram et al. 2021 Olivares et al. 2023
2

---

Published in Transactions on Machine Learning Research (12/2024)
To our knowledge, only three methods produce coherent probabilistic forecasts and allow models to be trained
end-to-end: Rangapuram et al. (202T 9. Olivares et al. (2023), and Das et al. [023]] (There is also parallel
research on hierarchical forecasts with relaxed constraints (Han et al. 2021 Paria et al.) 2021 Kamarthi
ct al. 2022 2024 Umagami ct al. 2023 but we limit our analysis of this line of work since our focus is
on strictly coherent forecasting met hods). In particular, Olivares et al. 2023 considers a finite mixture of
Poisson distributions that captures correlations implicitly through latent variables and does not leverage cross-time series information. Donti et al. (2021) achieves coherence by constructing probabilistic forecasts for aggregated time series through the completion of equality. On the other hand, Rangapuram et al. 22021 leverages multivariate time series information but achieves coherence through differentiable projection layer, which could degrade forecast accuracy: it does not directly model correlations between the multivariate outputs, but rather couples them through its projection layer. Dedicated effort is still necessary to capture these hierarchical relationships to improve forecast accuracy, Such probabilistic methods can benefit from the ability optimize for arbitrary loss functions by differentiating samples, as demonstrated by Rangapuran et al. [2021) The capacily to optimize any loss computed from the forecast. samples can help align the forecasting system's goals with the neural network's learning objective.
Summarized, an idcal hicrarchical forccasting method should satisfy several desiderata: 1) be coherent end- to-end; 2) model a joint multivariate probability distribution, capturing the intricate relationships between series within the hierarchy; 3) leverage cross-time series information and accurately reflect these relationships; and 4) generate differentiable samples, to enable the method to optimize for arbitrary learning objectives that align with the forecasting system's goals. In this paper, we present method that satisfies all these ideal properties; see Table [ for summary.
We introduce the Coherent Learning Objective Reeparamederizationn Newral Network (CLOVER), method for producing probabilistic coherent forecasts that satisfies all the desired properties stated above. Our main contributions are the following.
We design a multivariate Gaussian factor model that achieves forecast coherence exactly by con- struction, representing the joint forecast distribution over the bottom-level series. This factor model is generic and can augment most neural forecasting models with minimal modifications, with the condition that all the hierarchy's series available during training and inference.
In contrast to other joint distribution models, our factor model samples' differentiablltyy-nnna enabled by the reparameterization trick offers greater versatility, allowing optimization of a wider range of learning objectives beyond negative log likelihood. We leverage this versatility by training the model using both the Continuous Ranked Probability Score (CRPS) and the Energy Score (Matheson & Winkler 1976). For hicrarchical forccasting tasks, aligning the training objective with the CRPS evaluation metric leads to significant improvements in accuracy.
3. We combine the multivariate factor model with an augmented MQCNN architecture (Wen al. 2017 Olivares et al. 2022c) into the Coherent Learning Objective Reparamcterization Newral Network (CLOVER). Our augmented architecture accepts vector autoregressive (VAR) inputs through multi layer perceptron module.
CLOVER achieves state-of-thheart performance in six publicly available benchmark datasets, demon- strating substantial improvement in both probabilistic and mean forecast accuracy. On average, CLOVER putperforms the second-best. method by 15 percent in probabilistic forecasting accuracy and delivers a 30 percent increase in mean forecast accuracy.
3

---

Published in Transactions on Machine Learning Research (12/2024)
TToall
1
1
0
0
[a]]b]]
0
1 0 0
1
Y1,
2z
Y3,r + Yau
0
0 0
0
1
0
V1.
Y4.
1
(a) Graph representation
(b) Matrix representation
Figure 1: A simple time series hierarchical structure with Na aggregates over No bottom time scries. Figure la shows the disaggregated bottom variables with blue background. Figure 1b (right) shows the corresponding hierarchical aggregation constraints matrix with horizontal lines to separate levels of the hierarchy. We decompose our evaluation throughout the levels.
2 Hierarchical Forecast Task.
Notation. We introduce the hierarchical forecast task following Olivares et (0023) We denote a hierarchical multivariate time series vector by Ypi, |iale y1b,a NNN,+N,, where [a] [6]. [a], and [6] denote the set of full, aggregate, and bottom indices of the time series, respectively. There are | = Na No time series in total, with |[a]] = Na aggregates from the [[]]] N, bottom time series, at the finest level of granularity. We use as time index. In our notation, we keep track of the shape of lensors using square brackets in subscripts. Since each aggregated time series is a linear transformation of the multivariate bottom series, we write the hierarchical aggregation constraint as
[+] [ab
Y(i, = MMM
= Y(a)..
(1)
The aggregation matrix A Aalb represents the collection linear transformations for deriving the aggregates and sums the bottom series to the aggregate levels. The hierarchical aggregation constraints matrix Skipp is obtained by stacking A and the N, N, identity matrix Iыjы-
For simple example, consider N, = 4 bottom-series, so [b] = 2,,2,,,44) and 3Tooaal = Σ: Vi.-- Figure shows an example of such hierarchical structure, where the multivariate hicrarchical time series defined as:
(ale [Vrotal, V1.4 V2,1 Va, va.al".
YM,. [V1,e, y2,1 Va,t, va.el
(2)
Probabilistic Forecast Task. Consider historical temporal features x (h)
[b]]] known future information
(flt+1++ and static data (e) forecast creation date and forecast horizons in + 1 No)- A classic forecasting task aims to estimate the following conditional probability
(f)
Pi Y 1g,tto
( 1l(t-1+++,,)
X for 17 = 1,...,,.. N
[i
The hierarchical forecasting task augments the forecast probability in Eqn. with coherence constraints
in Eqn. 2.1 [Bnn Taieb et al.] 2020 Panagiotelis et al.
2023 Olivares el al. 2022c), by restricting the
probabilistic forecast space to assign zero probability to non-coherent forecasts. Definition 2.1 formalizes the intuition, stating that the distribution of given aggregate random variable exactly the distribution defined the aggregates of the bottom-series distributions through the summation matrix Skj.--
The conditional independence assumption across in Eqn. maintains the computational tractability of the forecasting task, and an assumption used forecasting methods.
4

---

Published in Transactions on Machine Learning Research (12/2024)
Definition 2.1. Let (Rpp). Fi)) (p)) be probabilistic space (on the bottom-level series Y(j). with sample space Rbb. event space Flb), and Pb] a forecast probability. Let ((gjy :21b + S4g be the linear transformation implied by the aggregation constraints matrix. A coherent forecast space (S2 Fij P(t)) satisfies
P (Swr (B).
(4)
for any set B E Fig and image San F()--
Hierarchical Forecast Scoring Rule. In this work and most of the hierarchical forecast literature, the performance of probabilistic forecasts primarily evaluated by the Continuous Ranked Probability Score
(CRPS), e.g. Ben Taich al. (2017); Rangapuram al. (2021) : Olivares et al. (2023); Panagiotelis al.
(2023): : Das et al. (2023); Wickramasuriya (2023). The CRPS between a target y and distributional forecast. Y defined as
1
CRPS (v. Y) = Ey [| - vll 2--2BBx By, -
where Y' is distributed as Y, but is independent of it (Matheson & Winkler
1976
Gneiting & Raftery 2007
CRPS is commonly used because it is a strictly proper scoring rule (Gnciting & Raftery 2007 and is agnostic to model and distributional assumptions. The metric also gives a summary performance on all quantile forecasts (Laio & Tamea 2007), as
CRPS(y. Y) Y)= 2 / QU,
(6)
where F (q) is the q-quantile of variable Y. The q-quantile loss is defined as QL, (4) = q(v F (a))+ + (1 (F0 (g) - y)
(7)
3 Methodology
In this section, we describe main method, the Coherent Learning Objective Reparameterization Newral Network (CLOVER). It consists of a multivariate probabilistic model, underlying neural network structure, and a coherent end-to-end model estimation procedure.
3.1 Coherent Probabilistic Model
Our predicted probabilistic forecasts at all hierarchical levels are jointly represented by Gaussian factor model. Our neural network maps the known information (past, static and known future) to the location, scale, and shared factor parameters, and the forecasted factor model parameters are designed to model correlations between the bottom-level series, while conditioning on all known information. Our factor mode combined with the coherent aggregation in Equ. [0 directly estimates the multivariate probability of Nottom-leve series (h)
(plttea,, conditioning on historical, nhown-future, and static covariates x [b]la])
x (f)
X0] i.e.,
(s)
[b]1++1+++N))
P [[(++1N
(n) [p][ej] x (p(+ x x = IP (S (M (iR4L4NN)
Mooiin a (IAA) ([MMM))))
(8)
At a given forecast creation date 6, the model uses the location Anpe) E N. scale d NoxN and shared factor F PonRMMJ N.N.N parameters, along with samples from standard normal variables
Zbii. ~ (0 Ip(()) and Epp),n ~ (( (uje)) to compose the following multivariate variables for cach
horizon:
1 = (6))),). + Diag( G 1p.. + 1,
1 L...,,
99)
Farly work factor forecast models augmenting neural networks done by Wang al. 2019 does ensure coherence.
5

---

Published in Transactions on Machine Learning Research (12/2024)
Decoder
Splret
S
Denoiar
Decodier
Deroder
Sosss Posne
Diagi Monn Pontiie PAl
r-2
--1
Len Mll
Loes TAlp42
b pj
Dists
degs pir
ThmpConollrion Cre
MILPi
MLLP
RD 161213
L 181122
40
wif 1M
(R)) (Njel
wesk [e]
VD TRE
Encoder
The Coherent Learning Objective Reyuanumelerization Neural Network is Sequence-lo-Sequeee with Context network that uses dilated temporal convolutions as the primary encoder and multilayer perceptron based decoders for the creation of the multi-step forecast. CLOVER coherently aggregates the samples of the factor model f1,,.,, = (aln & j... We mark in red the standard normal samples that are parameter-free, the reparameterization trick allows to apply backpropagation through the factor model outputs. CLOVER extends upon the univariate MQCNN, through the cross series multi layer perceptron.
Figure 2:
After sampling from the multivariate factor, we coherently aggregate the clippod outputs of the network,
M.n. = Sp1 ($pM))=))
(10)
The shared factors enable the factor model to capture the relationships across the disaggregated series, and the covariance structure of the disaggregated series follows:
Cov (8p)...) = Diag( (I.o.)) + PiM.. [blk),,
(11)
Definition 3.1. Any bottom-level multivariate distribution can be transformed into a coherent distribution through coherent aggregation 3 Given a sample 8t6 ~ Proj. coherent PN distribution can be constructed with the following sample Iransformation
[ (KM (8()).
(12)
In other words, it is enough to aggregate the bottom-level forecasts in bottom-up manner. We include in Appendix A proof of the approach's 's coherence properly and the covariance structure.
3.2 Neural Network Architecture
As mentioned in Section [ the factor model can augment most neural forecasting models if all series in the hierarchy available during training and inference, as this is sufficient condition to obtain the factor model parameters. Appendix shows augmented subset of NeuralForecast models (Olivares et al. 2022b
Coherent aggregation can thought of special of bootstrap reconciliation Panagiotells that only relies bottom-level forecast. distribution.

---

Published in Transactions on Machine Learning Research (12/2024)
while the main paper focuses on the MQCNN-based model (Wen et al. 2017 Olivares et al, 2022c We select MQCNN for its outstanding performance in multi-step forecasting problems. Our architecture based on MQCNN has a main encoder which consists of a stack of dilated temporal convolutions, and it is applied to historical information for all series. In addition, it uses a global multi layer perceptron (MLP) to encode the static and future information. The encoder at time Equ. 13 is applied to each disaggregated series:
(h)
TempConvolution
( (1) (8) (01:aa))
f
h jb]
(s)
MLP1 (xx(²
(13)
(f) b).
= MLP₂
(x(0)
We use residual cross series MLP to capture vector autoregressive relationships in the hierarchy with minimal modifications the architecture D
(h) M
MrossScriesMLP
(b (h)
(14)
(1.+))
CLOVER uses a two-stage MLP decoder: the first decoder summarizes information in the horizon-agnostie (ag) context (b), and the horizon-specific context Cbb[M]e and the second stage decoder transforms the contexts
[sp)
into the factor model parameters (MiMM (bllk), ((lliii)))
Equ. 15 describes the operations:
(h)
(+)
(
[6),,
(ag)
MLP3 (Mm),)
and
[sp) [b[ih],
= MLP4 (h(pa))
(15)
[b],c
(poli (ill), ((M)))))) MLPs
(sp)
(ag) (b)a
( (o11+11++
(b)
For the last step, the network composes the factor model samples using Equ.
g and aggregates them in Equ. 10
3.3 Learning Objective Reparameterization
Let 0 be a model that resides in the class of models θ defined by the model architecture. Here, 0 can be thought as mapping from the model feature space to the parameter set for all target horizons, we have
(n) (Moj))) pilN) F (IRIM)) 0xx [M(a]] (0|++1+N
(16)
Let (0) be the random variable parameterized by O. In some problems, multi-step coherent forecasts for multiple items needed (c.g., in retail business, coherent regional demand forecasts are required for each product). Let u be the index of such an item within an index set {1,..., N₂) of interest, and let (0) be the coherent forecast random variable for the target Vu,isnn- As we explain in Appendix [ using the reparamcterization trick (Kingma & Welling 2013), within the class of paramcters θ defined by the ncural network architecture, we optimize for either CRPS
min OEE
Σ CRPS (yu.sem Pu)m)(())= Σ
E [luneeto Yusmall L Puiin Pusnel
(17)
M.F.m.A
u.i.m.
the Energy Score
min Σ ES (Pa,)) (.I,H((())) ece
Σ [1ba..10 - Y 1u 2 (
B [1
V
(18)
We train the model using stochastic gradient descent (Adam Kingma & Ba 2014 with carly stopping Yao
et al. 2007 Appendix C contains details of the network's optimization and hyperparameter selection.
"Temporal exogenous data only aggregates the target signal, other features (e.g., calendar) maintained without aggregation. 5 Approach also inspired by regional forecast Amazon's Supply Chain Optimisation Technologies (Savorgnan 2024
7

---

Published in Transactions on Machine Learning Research (12/2024)
3.4 Discussion
Here, we discuss differences between CLOVER (our method) with two coherent end-to-end probabilistic forecasting basclines in HierE2E (Rangapuram al. 2021 and DPMN (Olivares al. 2023
The HierE2E method of Rangapuram el al. (2021) is "loo general." It consists of augmented DeepVAR neural network model (Flunkert et al. 2017) that produces marginal forecasts for all hierarchical series. HierE2E claims to be more general than hicrarchical forecasting, sinco it is designed to enforce any convex constraint sat- isfied by the forecasts; due to the constraining operation in the method, it has to revise the optimized forecasts. It does not leverage the specifics of the hierarchical constraints, which are more structured than a general convex constraint. HierE2E produces samples from independent Gaussian distributions for each time-series in the hierarchy; since the samples are not guaranteed to be hierarchically coherent, HierE2E couples samples by projecting them onto the space of coherent probabilistic forecasts. Both the sampling operation (Kingma. & Welling 2013) and the projection are differentiable, allowing the method to be trained end-to-end; HierE2E restricts to optimize likelihood or CRPS. Both HierE2E and CLOVER can allow different distribution choices, since Gaussians can be replaced by any distribution that can be sampled in a differentiable way, ie.. almost any continuous distribution (Ruiz et al. 2016 Figurnov et al., 2018 Jankowiak & Obermeyer 2018). In Rangapu- ram et al. 2021), the projection operator ensures coherence, and correlations between bottom-levels are learned only by optimizing the neural network. In contrast, CLOVER produces forecasts for bottom-leve series only, while relying common factors to encode correlations. This removes the need to forecast all levels simulta- ncously, therefore, reducing computational requirements if we are only interested in a subset of the aggregates.
On the other hand, the DPMN baseline (Olivares et al. 2023) is "too restrictive," " in particular, as a Poisson mixture can be prone to distribution misspecification problems. It is known that when a probability model is misspecified, optimizing log likelihood equivalent to minimizing Kullback-Leiber (KL) divergence with respect to the true probabilistic distribution, KL divergence measures change in probability space, while optimizing CRPS is equivalent to minimizing the Cramer-von Mises criterion (Gneiting & Raftery 2007);
which quantifies the distance with respect to the probability model in the sample space. CLOVER's learning objective for the probabilistic model is resilient to distributional misspecification (Bellemare et al. 2017 Moreover, CLOVER can be optimized to accommodate other evaluation metrics of interest.
Finally, DPMN estimates the covariance among time series, but does not take advantage of the multivariate input when encoding historical time series. Similarly to other ARIMA-based baselines, on specific hierarchical benchmark datasets such Traffic, DPMN produces suboptimal bottom-series forecasts. We improve the encoder for historical time scries by adding CrossSeriesMMP after the Temporal convolution encoder, which bridges the accuracy gap between HierE2E and our MQCNN-based approach.
4 Empirical Evaluation
In this section, we present our main empirical results. First, we describe the empirical setup. Second, we evaluate CLOVER and compare it with state-of-the-art hierarchical forecast models. Third, we present the results of the ablation study that further analyze the source of improvements in variants of CLOVER.
4.1 Setting
Datasets. We analyze six qualitatively different public datasots: Labour, Traffic, Tourism-S, Tourism-L, Wiki, and Favorita, each requiring significant modeling flexibility due to their varied properties. The Favorita dataset is the largest dataset evaluated which includes count and real-valued regional sales data, with over 340,000 series. The Tourism-S and Tourism-L datasets, report quarterly and monthly visitor numbers to Australian regions, respectively, and they are grouped by region and travel purposc. The Traffic dataset contains daily highway occupancy rates from the San Francisco Bay Area, leaturing highly correlated series with strong Granger causalities. The smallest Labour dataset tracks monthly Australian employment by status, gender, and goography, the series in this dataset are highly cointegratod. Lastly, the Wiki datasct summarizes daily views of online articles by country, topic, and access type. We provide more dataset details in Appendix D
8

---

Published in Transactions on Machine Learning Research (12/2024)
Dataset
Items (N₂)
Bottom (N6) Levels Aggregated (Na + N6)
Time range
Frequency Horizon (Na)
Labour
32
57
2/1978- 22/2019
Monthly
12
Traffic
200
207
1/2008-3/ 2009
Daily
Touriam-S
56
89
1998-2006
Quarterly
Tourism-L
304
4/5
555
1998- 2016
Monthly
Wiki
150
5
199
1/2016-12/
Daily
7
Favorita
4036
54
93
1/2013 8/2017
Daily
Table 2: Summary of publicly-availabl data used in empirical evaluation.
Evaluation metrics. Our main evaluation metric is the mean scaled CRPS from Equ. defined as the score described in Equ. 19 divided by the sum of all target values. Let (19) be vector of length N. + No consisting of binary indicators for hicrarchical level 9, where for each j [t]. 1/9) = 1 if aggregated series j is included in the hierarchical level g, and 0 otherwise. Then sCRPS for the hierarchical level g is defined as
Eist
CRPS( Vine)) (69)
sCRPS (y114)144))) YIDD ((0))=
(19)
N.+ Di
|| -1.0)
Baseline Models. We compare method with the following coherent probabilistic methods: (1) DPMN-
GroupBU OOivares et al.] 2023). (2) HierE2E (Rangapuram et al. 2021 (3) ARIMA-PERMBU-Mint (Ben)
Taieb et al. 2017), (4) ARIMA-Booten MLIMA-Bootttrap-BU (Panagiotelis et al. 2023), and (5) an ARIMA. In addition, in
Appendix G we compare our method with the following coherent mean methods: (1) DPMN-GroupBU, (2)
ARIMA-ERM Ben Taieb & Koo 2019), (3) ARIMA-Min! (Wickramasuriya et al. 2019), (4) ARIMA-BU, (5) an
ARIMA, and (6) Scasonal Naive. We use the implementation of statistical methods available in the StatsForecast.
and HlierarchicalForeass libraries OOlivares et al. 2022c Garza al. 2022
4.2 Forecasting Results
As mentioned above, we compare the proposed model the DPMN (Olivares 2023), the HierE2E (Ranga puram et al. 2021), and two ARIMA-based reconciliation methods (Wiekramasuriya et al. 2019 Panagiotelis
Table 3 The reconciliation results for ARIMA are generated using HiierarchiaalForecastt Olivares et al. (0022c), with each confidence interval calculated based on 10 independent runs. The results for HierE2E are generated based on three independent runs using hyperparameterss tuned by Olivares et al. (2022c). All metrics for DPMN are quoted from Olivares et al. (2023 with identical experimental settings on all datascts.
et al. 2023). Following previous work, we report the sCRPS at all levels of the defined hierarchies;
As the Overall row Table 3 shows that CLOVER improves sCRPS compared to the best alternative five of six datasets, with gains of 27.67% on Favorita, 10.24% on Tourism-S, 4.16% on Tourism-L, 4.5% on Wiki and 54.40% in Traffic. There is degradation of 8% in the reconciliation that performs the best on Labour, the smallest data set. For five out of six datasets, our model achieves the best or second-best sCRPS accuracy across all the levels of the hierarchy; it is important to consider that aggregate levels are much smaller in sample size for which we prefer the bottom-level measurements as an indicator of the methods' accuracy.
In Traffic, our model achieves remarkably better results than the DPMN and HierE2E baselines, which we explain by the ability model VAR relationships accurately because of the smoothing the time serics features before using them inputs for other series. For the smallest dataset, Australian Labour, all deep learning models, including CLOVER, showed a decline in accuracy compared to statistical baselines. We hypothesize that this performance drop is due to the limited data available, as 57 monthly series may not be sufficient. to effectively train complex models like deep learning architectures. Another explanation for the Labour degradation is the presence of strong trends and cycles, which can lead to artificially high correlations between the series, resulting spurious relationships that may confuse the VAR modules.
We complement the main results with mean forecast evaluation section in Appendix G As we see in Table the sCPSS accuracy gains mostly mirrored in the relative squared error metric (relSF). We qualitatively
9

---

Published in Transactions on Machine Learning Research (12/2024)
Table 3: Empiriaal evaluation of probabilistic coherent forecasts. Mean scaled CRPS (sCRPS) averaged over 5 runs, at each aggregation level, the best result. is highlighted (lower values are preferred). We report 95% confidence intervals, the methods without standard deviation have deterministic solutions.
resulte daffer from Ranzaparam
(2023) here the sCRPS quantile hyparparamettrr au
space fnnr korizon=4, instead Hiectl''s original
percent instead of
Bangapuram (2021). We
Peumu Mtat on Tousten- 4. to unavatlable because tox implemencation cannor natvely be applied to dasaantt with multiple heararchine, this is becaisee the bottom up ogga epatitito steutegy cannol ideusify unique way to obtain the upper level distribution..
DATA LEVEL
CLOVER (crps)
(energy)
CLOVER
DPONNGGrouPB
HierE2e
PERMBU-KinT
Bootstrap-BU
(not coherent)
ARIMA
Overall
0.0074++0.0000
0.0075=0.0000
0.0102++0.00
0.0171+.0003
0.0067=0.000
0.0076++0.000
0.0070
Laboar
1 (geo.)
0.00324:.0011
0.0035=0.0000
0.0033+:0 0011
2.0052:+.0003
0.0013=000001
0.0017#:0 0001
0.0017
[goo.)
0.0054++0.0000
0.0054+0.0004
0.0088++0.000
0.0181+00003
0.0045++0.000
0.0053++0.000
0.0044
3 (geo.)
9.0075:+0.0000
0.0074=0.000 0000
0.0115=#0 0009
0.0188:t.0000
0.0076=00.000
0.0086::0.000
0.0076
(geo.)
0.0137+0.0004
0.013820.0002
0.01732:0
0.0262:
0.0133=000000
0.0148:10.000
0.0138
Overall
0.0171+#00036
0.0228=0.008 00005
0. 0.0907::0 0024
0,0375+0.000
0058
0.0677=0.0001
0.0736: 0024
0.0751
Trafric
1 (geo.)
0.0026:+0.0012
0.006120.0077
0.0397 00011
0.018320.0000 0090
0.033120.0085
0.046810.000 0031
0.0376
2 (geo.)
0.0029:10.0014
0.006420.0065
0.0537.10 0024
0.01831:0.000
0.0341 0.0081
0.0483:10. 0.330
0.0412
3 (geo.)
0.0044:100.022
0.0075+0.0005
0.0538++0.000 +00022
0.0209+000071
0.0417+0.0061
0.0530+0.0022
0.0549
4 (geo.)
0.0587+:0.0106
0.0709=0.0104
0.2155: .0222
0.097440.002
1.22100.0027
0.1463:+0.000 0.01
0.1665
Tour 1an-8
Overall
0.0631#00.001
0.066920.0022
0.0710::00011 0122
0.0761#0.000
0007
0. 0.0812=0.0010
0,0703: 0017
0.0759
1 (geo.)
0.0268:10.0033
0.0304=0.0048
0.0329:10 0074
0.0400:10.000
0.0375::.0011
0.0335:10. 0026
0.0354
[gro.)
0.0484:10.0019
0.0507
0.0502+:0.000 +00.058
0.0609+0.0012
0.0679-
0.0507 0-0023
0.0709
3 (geo.)
0.0784:+0.0022
0.08 0.0817=0.0022
0.09144t0 1:1.4+00.000 0171
0.09 0000
0. 0.0959=#0.00 0012
0.0845: 10
0.0848
(geo.)
0.0980=+0.0017
0.1050=0.0034
0.1212::: 0238
0.1122:0.0000 0.007
0.123120.0012
0.1121
0.1125
Overall
0.1197++0.0033
126 0045
0.1249---0.0 .0200
0.1472++0.00
0029
0.1375+ 30
0.1416
1 (geo.)
0.0292++0.0112
0.029110.0008
0.0131#00 00-10
0.08-12::0
0.05
0,0622::0. 0026
0.0263
Touriax-i
2 (geo.)
0.0593:10.000
0.0598=0.0085
0.0637:10 0032
0.1012
0,009
0.0820:100. 0019
0.0904
3 (geo.)
0.10111+000000
0.1077=0.000
0.10841:: 0033
0.1317
0020
0.1207= 0010
0.1389
4 (geo.)
0.15:1::0.0046
0.1620=0.000 0060
0.1554:: 0025
0.1705:000000 00002
0.1646::0. 0007
0.1878
[prp.]
0.0594110.0076
0.0600=0.0041
0.0700.10 00038
0.099510.0000
0.0788110.00
0.0770
6 (prp.)
0.1100++0.0049
0.1140+0.0034
0.1070=+:0.00
1.1336+00000
0.1268+00.00
0.1
7 (prp.)
0.1824+0.0024
0.1934=0.000 0045
0.1887=00 0032
0.1955+0.000 0005
0.1949-+0. 0010
0.2022
8 (prp.)
0.259120.0050
0.2902=000054
0.2629::00000
0,2615520.0016
0,2698:10.00 0004
0,2831
Overall
0.2476=+000076
0.3184 0262
0.3158+t0..m 02400
0. 0.25922+:000 00001
000006
0.2816=+0. 0036
0.3907
1 (geo.)
0.0823:10.016
0.1536 0666
0.1709
0354
0. 0.100720.000
06116
0. 0.1886220011
0.1630::0 0..00000.006. 0065
0.1981
VIX1
2 (geo.)
0.1765++0.0154
0.2381
0.2299-
0241
0.1933++0000 0-0077
0.2091+0.0077
0.2192++0. 0043
0.2566
3 (geo.)
0.27-33=20.0108
0.3213=0.0311
0,3311
0230
0.278/1=0.000
0.10:19:0.0002
0.2923::0. 0032
0,4100
4 (geo.)
0.2814::00100
0.3326=0.0333
0.3370-10 0223
00900:10.003
0.4236620.0062
0.2999110 0032
0.4182
5 (prp.)
0.0227
0.5466=0.1400
50088100.0667
0.430710..003
0.0067
0.4339:10.000
0.6708
Favorica
Overall
0.2908:10.000
0.3279
0.40201:0.00
0.5298100.000
0.4670
0.4110:10. D.885
0.4373
1 (geo.)
0.1841110.0053
0.2209-00.0445
0.2760-10.0144
0.47141:00103
0.2692 0.0076
0.2900:10.0065
0.3112
2 (geo.)
0.27544++0.002
0.31 1
0.3865++0.8 0.0207
0.5182+0.0100
1.3824+0.0000
0.3877+ 0002
0.4183
3 [geo.)
0.29-15+00000
0.3313=0.0000 0053
0. 0.1068#0..o 0.20
0,5291=00 0129
0. 0.6838
0.44900 0098
0.4146
4 (geo.)
0.4092++0.0022
0.4484 000068
0. 1.5387+000255
0.6012+0.0131
0.5532 + 0.011 16
0.5749+0.00000
0.5749
10

---

Published in Transactions on Machine Learning Research (12/2024)
0.09
0.175
Validation SCRPS
0.05
Validation SCRPS
0.150
0.06
0.100
0.05
0.050
0.03
0.025
CLOVER
MOCNN
LSTM
NHITS
KREATS
TFP
PCRAGA
CRPS
NLL
Neural Forecast Architecture
Learning Objective
Figure 3: Ablation studies on the Bay Area Traffic dataset: a) In highly correlated hierarchies, VAR inputs enabled by CLOVER significantly improve over the univariate MQCNN's accuracy. b) The factor model CRPS learning objective remonst.ates clear advantages over classic negalive log-likelihood. Full ablation studies described in Appendix F
show the CLOVER forecast distributions for a hierarchical structure in the Favorita dataset in Appendix L We also include in Appendix E a cestaabiitt-poobablli plot comparing the similarity of the empirical and forecast distributions. As can be seen in Figure 7 the forecast distributions are qualitatively well calibrated.
4.3 Ablation Studies
To analyze CLOVER's source of improvements, we performed two ablation studies on the Traffic dataset, where we investigated first the effects of the learning objective and then the effects of VAR inputs enabled by the CrossSeriesMLP. Here we report a summary and refer to the details Appendix F
In the first study, we explore the effects of cross-scrics MLP that mimics the vector antoregressive model. We compare the CLOVER architecture with and without the cross-series multi layer perceptron (Cros)SeriesMLP) introduced in Equ. 14 Such module enables the network to share information from the series in the hierarchy with minimal modifications to the architecture. Table and Figure 3 show that CrossSeries....PP improved Traffic forecast accuracy by 66% compared to variants without it. We attribute the effectiveness of the VAR approach to the presence of Granger-causal relationships in traffic intersections. The VA VAR-augmented CLOVER improves upon well-established univariate architectures, including LSTM (Sak et al.) 2014), NBEATS (Oreshkin
ct al. 2020
Olivares al., 2022a), NHITS (Challu et al.) 2023 [Lim et al.) 2021 ), and FCGAGA Oreshkin
et al.
2021
spalio-tempora specialized archilecture.
In the second ablation study learning objectives, we compared CRPS-based optimization, described in Equ. 17 with the negative log-likclihood estimation for the Gaussian factor model introduced in Section 3.1
and other peelihood-estiimatee distributions. Table 0 and Figure 3 show that the ERPSSoptimized factor model improves forecast accuracy by nearly 60% when compared to the log-likelihood optimized model.
5 Conclusion
In this work, we present novel multivariate factor forecasting model, integrated with the MQCNN neural network architecture, resulting in the Coherent Learning Objective Regarameteriaation Neural Network (CLOVER). Using CRPS as the learning objective, we achieve significant improvements in forecast accuracy over traditional negative likelikeli objectives.
Our experiments on six benchmark datasets Favorita grocery demand, Australian quarterly and monthly tourism, Australian labour, Wikipedia article visits, and Bay Area Traffic -raffic-showed consistent improvements in CRPS accuracy, averaging over percent. We observed Australian labour dataset. Ablation studies further confirmed the importance the CRPS learning objective and the enhancements made to support multivariate time series inputs, which were key to the model's success.
percent performance drop on the smaller
11

---

Published in Transactions on Machine Learning Research (12/2024)
Our model's success datasets with strong Granger causality is due to its vector autoregressive capabilities. However, in the presence of cointegrated series, vector autoregresive approaches can be ineffective. Developing multivariate time series methods that are robust to cointegration is promising area for future research.
This study uses the reparameterization trick to enable alternative learning objectives such as CRPS and the energy score for neural forecasting tasks and the use of coherent factor models to capture correlations among hierarchical series structures. While we focus on parameterizing the multivariate predictive distribution as a Gaussian multivariate factor model, the framework's fiexibility can also accommodate other distributions that support sample differentiabilityy This is of special interest for outlier quantiles that cannot be well approximated by Gaussian variables. Exciting future rescarch directions include extending the reparameterization trick to handle discrete distributions, which could further enhance the accuracy of forecast distributions built on this framework. Another promising direction for future research is to extend the usage of the reparameterization trick from the learning objectives, into the hierarchical aggregation structure itself, provided the aggregation structure is done through differentiable transformations.
Finally, we observed performance challenges of neural forecasting models when working with limited dataset sizes. Traditional tools like data augmentation and pretraining need to be adapted to handle multivariate time series and associated covariates. Extending these techniques to the hicrarchical forecasting domain is a promising research direction.
Acknowledgements
This work was supported by the Amazon Supply Chain Optimization (SCOT) Forecasting Team. We thank the TMLR reviewers and the Action Editor for their helpful feedbacks and suggestions. We express our gratitude to Stefania La Vattiata for her contribution to the illustrations, which effectively captured the essence of our methods. Special thanks to David Luo for his work on hierarchical forecast baselines. Thanks to Boris Oreshkin for his pointers to latest developments in graph neural networks, and recommendations to enhance our ablation studies. We also appreciate Riccardo Savorgnan and the regional team for enriching discussions on alternatives for adapting the MQCNN architecture to multivariate time series. The authors thank Utkarsh for his suggestion on optimizing the closed-forn CRPS for truncated normal variables. We thank Youxin Zhang for his ideas to improve the computational complexity of the Energy Score and the CRPS estimators. Thanks to Lee Dicker and Medha Agarwal for their input on the energy score part of the paper and qualitative evaluation of the forecast calibration.
References
Oren Anava, Vitaly Kuznetsov, and (Google Inc. Sponsorship). Web traffic time series forecasting, fore- cast future traffic to wikipedia pages. Kaggle Competition, 2018. URL
c/
George Athanasopoulos, Rob J. Hyndman, Nikolaos Kourentzes, and Anastasios Panagiotelis. Forecast. reconeiliation: A review. International Journal of Forecasting, 40(2):430 456, 2024. ISSN 0169-2070. doi: https:://oo.orr//D0 11..://///1111111./11111 http:///ooooooo/.//..///o///oota att.ccle/pp//500662777000000 wwwffoecaat..223.1....o. oewp:///e/////////oeoooeeeeeeeee wtwp://iionndreee..ooo//ooeeeeae
URL https://wwwt
Australian Bureau of Statistics. Labour Force, Australia. Accessed Online, 2019. URL http:: btttps://www.ar abs
gov.au/AUse
Marc G Bellemare, Ivo Danihelka, Will Dabney, Shakir Mohamed, Balaji Lakshninarayanan, Stephan Hoyer, and Rémi Munos. The cramer distance as solution to biased wasserstein gradients. arXiv preprint arXiv: 10743, 2017.
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, pp. 1337-1347, New York, NY, USA, 2019. Association for Computing Machinery. TSBN 9781450362 doi: 10. 00.1145//292500333300 http://///.//////////...././...0
12

---

Published in Transactions on Machine Learning Research (12/2024)
Souhaib Ben Taieb, James W. Taylor, and Rob J. Hyndman. Coherent probabilistic forecasts for hierarchical time series. In Doina Precup and Yee Whye Teh (eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 3348 3357. PMLR, 06-11 Aug 2017. URL http:://proceedingg wwwpeiiiinnmmm......mmmmmmmmmm www.m.rrr.re... http./pppeee//Ottaaalllaaaaaaaa
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:27 43, 2020. URL https://apl. seappp.ssmanccccoola... wrwaaar..rr/oooe 0rgg//orrrssD:42221477 432147722
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 Asso- ciation for the Advancement of Artificial Intelligence Conference 2023 (AAAT 2023), 2023. URL https://aaxiss
hts//bba/22000 12886
Kevin C. Chen, Lee Dicker, Carson Eisenach, and Dhruv Madeka. MQTransformer: Multi-Horizon Forecasts with Context Dependent Attention and Optimal Bregman Volatility, 8 2022. URL https://kcddmmiles
2etrss_22__paa________ 9623.pdf
Byron J. Dangerfield and John S. Morris. Top-down or bottom-up: Aggregate versus disaggregate extrapola-
tions. International Journal of Forecusting, 8(2):233-22) 241, 1992. URL htttp::/www.scieneeeiieeeee wwwenceeiieet.ccooo com/
science/ http:////rrrieee/////////////00/
A. Das, W. Kong, B. Paria, and R. Sen. Dirichlet proportions model for hicrarchically coherent probabilistic forecasting. In Robin J. Evans and Πya Shpitser (eds.), Proceedings of the Thirly-Ninth Conference on Uncertainty in Artificial Intelligence, volume 216 of Proceedings of Machine Leorning Research, pp. 518 528. PMLR, 31 Jul 04 Aug 2023. URL htttps:/pproceedings mlr ppess/v216/das2533 html
Priya L Donti, David Rolnick, and J Zico Kolter, De3: A learning method for optimization with hard constraints. arXiv preprint arXiv:2104. 2021.
Corporación Favorita, inversion, Julia Elliott, and Mark McDonald. Corporación favorita grocery sales forecasting, 2017. URL
Michael Figurnov, Shakir Mohamed, and Andriy Mnih. Implicit reparameteerization gradicnts. In Neural Information Processing Systems, 2018.
Valentin Flunkert, David Salinas, and Jan Gasthaus. DeepAR: Probabilistic forecasting with autoregressive recurrent networks. ArXiu, abs/1704.0411 170..04110, 2017.
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
Tilmann Gneiting and Matthias Katzfuss. Probabilistic forecasting. Annual Review of Statistics and Its
Application, 1(1):125-155) 1511 2014. URL sstps://ddoi.org/111 http://tt....rt.....1./////eeee/ 00t.://///mmmerrsttaaaaae1111110
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359-377 2007.
Xing Han, Sambarta Dasgupta, and Joydeep Ghosh. Simultaneously reconciled quantile forecasting of hierarchically related time series. In International Conference on Artificial Intelligence and Statistics, 2021.
Tao Hong, Pierre Pinson, and Shu Fan. Global energy forecasting competition 2012, 2014.
Rob Hyndman and Anne B. Kochler. Another look at measures of forccast accuracy. International Journal of Forecasting, 62(4)6679 688, 2006. ISSN 0169-2070. doi: http://www....aoo/////oooooooott URL
13

---

Published in Transactions on Machine Learning Research (12/2024)
Rob J. Hyndman, Roman A. Ahmed, George Athanasoponlos, and Han Lin Shang. Optimal combination forecasts for hierarchical time series. Compulational Stalistics & Data Analysis, 55(9):2511 (1)::2)7--288) 2011.
URL
Rob J. Hyndman, Alan J. Lee, and Earo Wang. Fast computation of reconciled forecasts for hierarchical and grouped time scries. Comput. Stat. Data Anol., 97(C):16-32, may 2016. URL https://dot.org::
1016/j. 2015. 11.007
Rob J Hyndman, George Athanasopoulos, Azul Garza, Cristian Challu, Max Mergenthaler, and Kin G. Olivares. Forecusting: Principles and Practice, the Pythonic Way. OTexts, Melbourne, Australia, 2024. available
Martin Jankowiak and Fritz Obermeyer. Pathwise derivatives beyond the reparameterization trick. ArXiv, 2018.
Jooyoung Jeon, Anastasios Panagiotelis, and Fotios Petropoulos. Probabilistic forecast reconciliation with applications to wind power and electric load. European Journal of Operational Research, (73(2):364-77) 2019.
Harshavardhan Kamarthi, Lingkai Kong, Alexander Rodriguez, Chao Zhang, and B Aditya Prakash. Profhit: Probabilistic robust forecasting for hierarchical time-series. orXiv preprint orxi::220600090 2022.
Harshavardhan Kamarthi, Aditya B. Sasanur, Xinjie Tong, Xingyu Zhou, James Peters, Joe Czyzyk, and B. Aditya Prakash. Large scale hierarchical industrial demand time-series forecasting incorporating sparsity. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD '24 pp. 5230 5239, New York, NY, USA, 2024. Association for Computing Machinery. ISBN 9798400704901. doi: 10. 36775528..66163... URL
3671632
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. CoRR, abs/ has/1412.6980, 2014.
Diederik P. Kingma and Max Welling, Auto-encoding variational bayes. CoRR, abs/1312.6114, 2013.
Francesco Laio and Stefania Tamea. Verification tools for probabilistic forecasts of continuous hydrological variables. Hydrology and Earth System Sciences, 1144::126711277, 2007.
Bryan Lim, Sercan Ö. Arik, Nicolas Loeff, and Tomas Pfister. Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting, 37(4):1748 170111778866444,,,--111111- 2021. ISSN 0169-2070. doi: https:///oi.orr//10 1ttp:///1..1111....fcoaaaaaaaa. 1..10.0..0.....2201............1 URL http:://www.sciemceeeee...o www.ssceeneerrret.coomm com/ acciene//rracce///1/////000000o
Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos. M5 accuracy competition: Results, findings, and conclusions. International Journal of Forecusting, 38(4):1346 (344):1346-1344 1364, 2022. URL https:
Special Issue: M5 competition.
James E. Matheson and Robert L. Winkler. Scoring rules for contimuuu probability distributions. Management Science, 22:1087-1096, 1976.
Konstantin Mishchenko, Mallory Montgomery, and Federico Vaggi. A self-supervised approach to hicrarchical forecasting with applications to groupwise synthetic controls. ArXiv, htss/1966.15886, 2019.
Kin G. Olivares, Cristian Challu, Grzegorz Marejasz, Rafal Weron, and Artur Dubrawski. Neural basis expansion analysis with exogenous variables: Forecasting electricity prices with NBEATSx. International URL
Journal of Forecasting, 2022a. ISSN 2070, doi: http:://cooorr///0000/////oooeee http://hhh...aa.o.o..oooooooottp www.:e/ceeeneeeeeeccooooooo
Kin G. Olivares, Cristian Challú, Foderico Garza, Max Mergenthaler Canscco, and Artur Dubrawski. NouralForecast: User friendly state-of-thh-art neural forecasting models. PyCon Salt Lake City, Utah, US 2022, 2022b. URL. htt.:p/p//gghtho..aa///saaaaaaaa
14

---

Published in Transactions on Machine Learning Research (12/2024)
Kin G. Olivares, Federico Garza, David Luo, Cristian Challú, Max Mergenthaler, Souhaib Ben Taieb, Shanika L. Wickramasuriya, and Artur Dubrawski. HierarchicalForccnt A reference framework for hierarchical forecasting in python. Work in progress paper, submitted to Journal of Machine Learning Research., htts/2207.0033 .03517, 2022c. URL
03517
Kin G. Olivares, O. Nganba Meetei, Ruijun Ma, Rohan Reddy, Mengfei Cao, and Lee Dicker. Probabilistic hierarchical forecasting with deep poisson mixtures. International Journal of Forecasting, 2023. URL
Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-BEATS: neural basis ex-
pansion analysis for interpretable time series forecasting. Tn 8th International Conference on Learning Representations, ICLR 2020, 2020. URL https://opeerreieer net:/foruntddrrreemmdl
Boris N. Oreshkin, Arezou Amini, Lucy Coyle, and Mark J. Coates. FC-GAGA: Fully connected gated graph architecture for spatio-temporal forecasting. In The Association for the Advancement of Artificial Intelligence Conference 2021 (AAAI 2021), 2021. URL https::/arxii... org/ httgg/abs/2000 15531
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-7006 2023.
Biswajit Paria, Rajat Sen, Amr Ahmed, and Abhimanyu Das. Hierarchically Regularized Deep Forecasting, In Submitted to Proceedings of the 39th International Conference Machine Learning. PMLR. Working Paper version available at arXiv:2106.0748 00.6000, 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjic Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In Advances in Neural Information Processing Systems 32, pp. 8024-8035. Curran Associates, Inc., 2019. URL. http::/papeers wwwpaapprrse.uuiss...... coouurps.ccppapaa wttpiiigpppooooeeeeeeeeeeeee -eannnen
9015
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 International Conference on Machine Learning, pp. 8832 8843. PMLR, 2021.
Francisco J. R. Ruiz, Michalis K. Titsias, and David M. Blei. The generalized reparameterization gradient. In NIPS, 2016.
Hasim Sak, Andrew W. Senior, and Françoise Beaufays. Long short-term memory based recurrent neural net- work architectures for large vocabulary speech recognition. Computing Research Repository, abs/ 0ts//4402.122 2014. URL http:///axiv..rg/// 1t////112..111 1128
Riccardo Savorgnan, Chao Guo, Quenneville-Belair Vincent, Rahul Gopalsamy, and Kenny Shirley. Cross region attention for regional forecasting. Demand Forocasting Team, Supply Chain Optimization Tech- nologies (SCOT). Amazon Machine Learning Conference, 2024. URL RRLLhttss::/// httpp::///waaaa......mm amazon. wwtpeeccc./mmeiiiiieeeee Forecasting/Delphi
Thordis L. Thorarinsdottir and Tilmann Gneiting. Probabilistic Forecasts of Wind Speed: Ensemble Model Output Statistics by using Heteroscedastic Censored Regression. Journal of the Royal Statistical Society Series A: Statistics in Society, 173())3371-388, 11 2009. ISSN 0964-1998. doi: 101111/11467-993.22000
00616.x. URL https:://di.orrgeeeee
2009 00616.xx
Tourism Australia, Canberra. Tourism Research Australia (2005), Travel by Australians. https://www
http://wwwoo/.aaaaaaaaaaaaraarrr itw.riinnnnnaaataaaaaaalaaaaaa/a 2005.
Ryo Umagami, Yu Ono, Yusuke Mukuta, and Tatsuya Harada. Hiperformer: Hierarchically permutation-
equivariant transformer for time series forecasting, 2023. URL https://arxiis http:/axxxvvoorr//ba 2atps//5300.00033.33
15

---

Published in Transactions on Machine Learning Research (12/2024)
Steven R. Vitullo. Disaggregating time series data for energy consumption by aggregate and individual customer. Deportment of Electricnl and Computer Engineering, Ph. D. Dissertation., 2011.
Yuyang Wang, Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski. Deep factors for forecasting. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 6607 6617. PMLR, 09 15 Jun 2019. URL httpp::/pproceediggg mlr wwt.eperrs//////ss//e//e htmll
Ruofeng Wen, Kari Torkkola, Balakrishnan Narayanaswamy, and Dhruv Madcka. Multi-horizon Quantile Recurrent Forecaster. In 31st Conference Neural Informalion Processing Systems NIPS 2017, Time Series Workshop, 2017. URL https:://aaxw...og///aaaaao wttp://wwww.aaab//....... 11053,
Shanika L. Wickramasuriya. Probabilistic forecast reconciliation under the Gaussian framework. Accepted at Journal of Business and Economic Statistics, 2023.
Shanika L. Wickramasuriya, George Athanasopoulos, and Rob J. Hyndman. Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. Journal of the Americon Statistical Association, (14((56):8044888))) 819, 2019.
Malcolm Wolff, Kin G. Olivares, Boris Oreshkin, Sunny Ruan, Silan Yang, Abhinav Katoch, Shankar Ramasubramanian, Youxin Zhang, Michael W. Mahoney, Dmitry Efimov, and Vincent Quenneville- Bélair. SPADE : Split peak attention decomposition. In Thirty-Eighhta Annual Conference on Neural Information Processing Systems NewrIPS 2024, volume Time Series in the Age of Large Models Workshop, Vancouver, Canada, 2024. NeurIPS 2024. URL https: http://aaxxw....oaabeeess/
05852
Yuan Yao, Lorenzo Rosasco, and Caponneto Andrea. On early stopping in gradient descent learning. Constructive Approrimation, 26(2):289-33) 211, 2007. URL
Michaël Zamo and Philippe Naveau. Estimation of the Continuous Ranked Probability Score with Limited Information and Applications to Ensemble Weather Forecasts. Mathematical Geosciences, 50(2):209 234, February 2018. doi: 10.1077/111044-017:9700-- URL https: //hal, science/hal-02297642
16

---

Published in Transactions on Machine Learning Research (12/2024)
A Multivariate Factor Model Coherence and Covariance
A.1 Coherent Aggregation Properties
The coherence of CLOVER is a special case of the bootstrap sample reconciliation technique Panagiotelis et al. 2023 as explored by Olivares et al. (2023).
Lemma A.1. Let (Rt Fiu; P(j)) be a probabilistic forecust space with Fib) a a-algebra in Mpp If a forecast distribution P13] assigns a zero probability to sets that do not contain coherent forecasts, it defines a coherent probabilistic forecast space (R. P(j)) with [ip [10p))
Pla] (y $ A Aal (B) [3) = 0
F (1) = P (B)) AB E FM..
(20)
Proof. We note the following:
Pg (B)) P|a1 (
A [allb]]
A lall (B) n
IR...
B (1))
(B) -P
R
F[a] to]b (B) B3) P(b) (B) (1
P(a] (y[a)
A (o B)) x PM (B) PM (B).
The first equality is the image of set B Sp: corresponding to the constraints matrix transformation, the second equality defines the spanned space as subspace intersection of the aggregate series and the bottom series, the third equality uses the conditional probability multiplication rule, and the final equality uses the zero probability assumption.
By construction the samples of our model У1 = assumptions of the lemma and proving the coherence of our approach.
San (fre)) and $1a = Ajali (910)) satisfying the
A.2 Covariance Structure
Here we prove the covariance structure of our factor model introduced in Section 3.1 Lemma A.2. Let our factor model be defined by
Spi.s P.pl).m. Diag( (0))) 2pp. Fl]]A) 17 = 1, ,,,,,,,,,,,,,,,,, Nn.
(21)
with independent factors Zibi.n
NO (paiee, and E(H),n NN(Op]) Iupa),, its covariance satisfies
Cov (8 ())))a))) Diag( (b.o.)) + 10 F [i
(22)
Proof. First, we observe that
Cov ( Mp,,))) Cov (Diag( (DDiag((()),m) (b) ( (b). Diag ( (2)))))))))))
2Cov
(Diag(O M. m. (I
(23)
Cov
F (M (111),.,
By bilincarity of covariance and independence of the sampled factors, it follows that
Cov (0p)).as pmm)) Diag (). Cov (z(0).m. 2(1)..s) Diag ( ( Fi]]MM,, Cov (qpp)mer (M))))
BT
We conclude that
Cov (19...::16))m. Diag(00)),. F Mo FT (o
(24)
17

---

Published in Transactions on Machine Learning Research (12/2024)
B Code Script for Sampling
sample(self, distr_args, sindov_size, nun_sampees-None):
one
**Parameter***
"dietr_arge Forecast Distribution argusents.
'vindou_ aize int=1, for reconciliation reshapes in sample method, 'num_ samples : int=500, number of samples for the empirical quantiles.
**Returne**
"samples's : tensor, shape (B,8, 'auu samples '].
"quantiles'; tensor, empirical quantiles defined by 'levels'.
factor, loading, stds distr_args callapsed_batch, H,
size()
(collapsed, batch, (B, -> Ns,N,
factor _ooddigg factor factor_loadee www.looddd....leee.e reshape
(-1, self __series, vindou_size, н, self.n_ f_ctoos))) contiguous factor_loading torch.einsee einsum ("iv. www.bbuhf--aae -bbhht"' self.8P, factor (
means means teans.reshin reshape gelf self.n_sete series, vindow_size H, 1) contiguous
stds = stds stds.reehapp reshape soriie,,, sindov_aize, H, 1). contiguous ()
# Loading factors for covariance (iag(atda) F -> (8PF) (SPFrt)
hidden_facton factor Normal(
boc=torch.zeron
(factor_loading shape [0), vindov_size, H, self n_factors),
scale=1.0)
sample_factore factors hidden_facton rsample shapee (self nun_ samp__mmmee sample_factors sample_factors per.aacctorr...u.eeeee...ee pernute(
(1, 3, 0)
) contiguous () (n_items, sindow_size, nunnsampples
sample wwwded_ddeeeee.... factors
einsum ("bvahf bwhfn buhfn->bvuhn", factor_loadiae
sample_factora
8 sample lwwdedddddddann means Baadl__roded_facttrrr
# Sample Normal
normal Normal (loc=torch. zeros_ like loaded_menns)) scale-1.0)
samples normal rsample() () samples
atds samples)
samples torch. einsum i"iv ("iv, samples reshape (collapsed_ batch, H, self num_samples). contiguous ()
bww.n b bbishhn,, self SP, sanples)
# Compute quantiles and mean guantiles_device self www.uuuiiee.ttoeeeeeeee..eeeeeee device)
quants torch quantile (input=samples, qrquantiles _eaaatiies_ddeieee___e_ quants quants pernute ((1, 0)) # [Q [, sample_mess torch. (eamples, dim=-1, keepdim-True) samples, sample_mean, quants
device, dim=-1)
Figure 4: PyTorcc function for sampling from Gaussian Factor model. Note that the factor samples are shared across all bottom-level distributions. The samples are differentiable with regard to the function inputs. We can easily adapt this function to sample from other distributions.
18

---

Published in Transactions on Machine Learning Research (12/2024)
Table 4: Coherent Learning Objective Reparameterization Neural Network (CLOVER) hyperparameters. We use small or large model configuration depending on the datasets' size.
PARAMRTER
Considered Values
TRAPPIC
TouRLSN-S
WikI
Lasoua
ToURISH-L
FAVORI7A
Activation Function.
ReLU
ReLU
ReLU
ReLU
ReLU
ReLU
Static Encoder Dimension.
10
Temporal Convolution Channel Size.
10
10
10
20
30
30
Future Encoder Dimension.
20
20
Horizon Specifle Decoder Dimensions.
5
5
5
Horizon Aguostic Decoder Dimensions.
20
10
20
10
20
Factor Model Components.
10
10
20
20
10
Cross Series MLP Hidden Size.
200
SGD Batch Size.
1
SGD Effective Batch Sixe.
207
199
57
555
744
SGD Max steps.
2e3
1e3
2e3
2e3
263
80e3
Early Stop Patience steps.
Learning Rate.
5e-3
5e-4
5e-4
1e-3
5c-4
5e-4
C Training Methodology and Hyperparameters
Here, we complement and exlend the description of our method in Section
3
To avoid information leakage, we perform ablation studies in the validation set preceding the test set, where we explored variants of the probabilistic method, as well as its optimization. We report these ablation studies in the Appendix F] For each dataset, given the prediction horizon h, the test set is composed of the last. h time steps. The validation set is composed of the time-steps precoding the test set time range. The training set is composed of all dates previous to the validation time-range. When reporting the final accuracy results of our model in the test set, we used the settings that perform the best in the validation set.
We minimally tune the architecture and its parameters varying only its size and the convolution kernel filters to match the seasonalities present in each dataset. For the data set Favorita, we use dilations of [1, 2, 4, 8, 16, 32) to match the weekly and monthly seasonalities. For the Tourism-L and Labour datasets we use dilations of [1, 2, 6, 12) to match the monthly and yearly seasonalities. For the Traffic dataset we use dilations [1, 7, 14, 28] as multiples of to match the weekly scasonalitics. For the Tourism-S dataset we use dilat ions of [1, 2, as multiples of to match the quarterly seasonalities. For the Wiki dataset we use minimal dilations [1.
The selection of the number of factors mostly follows the memory constraints of the GPU, the effective batch size implied by our probabilistic model grows rapidly as function of the multivariate series. In the Favorita dataset, more factors are likely to continue to improve accuracy but with the tradeoff of the computational speed. Similarly, the Cross-series MLP hidden size is selected following the GPU memory constraints.
We share a learning rate of 5c-4 constant across the three datasots, which shows that the method reasonably robust across different forecasting tasks. During the optimization of the networks we use adaplive moment stochastic gradient descent (Kingma & Ba 2014 with early stopping (Yao et al. 2007) guided by the sCRPS signal measured in the validation set. We use a learning rate scheduler that decimates the learning rate four times during optimization (SGD Maxsteps / 4), to ensure the convergence of the optimization.
The CLOVER model is implemented using PyTorch (Paszke et al. 2019), with the NeuralForecast library framework OOivaares al. 2022b). We run all experiments using single NVIDIA V100 GPU.
As mentioned earlier in Section 4.1 statistical methods available in StatsForeeast. and HierarchicalForeeast
libraries (Olivares et al., 2022c; Garza et al., 2022). In particular, use available code the hierarchical baselines repository and Inerarchical datasets repository As described in Table 2 we deviate slightly from the experimental setting in Rangapuram et al. (2021) h to ensure the reproducibility of the results of this paper. For Table g we rerun the HierE2E (Rangapuram et al. 2021 baseline on Labour, Traffic, Wiki and
Tourism-S using their best reported hyperparameters.
19

---

Published in Transactions on Machine Learning Research (12/2024)
10
15
20
25
30
0
10
20
30
40
50
0
100
200
300
10
50
100
150
10
100
20
50
200
30
100
30
40
300
400
40
200
60
70
500
80
(a) Labour
(b) Traffic
(c) Tourism-S
(d) Tourism-L L
0
10
20
30
40
50
0
20
4D
60
80 100
120
20
40
150
(e) Wiki
(f) Favorita
Figure 5: Hierarchical constraints of the empirical evaluation datasets. (a) Labour reports 57 series number of employees by full-time status, gender and geographic levels. (b) Traffic organizes the occupancy series of 200 highways into quarters, halves, and totals. (c) Tourism-S categorizes its 89 regional visit series based on travel purpose, zones, states, and country-level aggregations and urbanization within regions. (d) Tourism-L categorizes its 555 regional visit series based on travel purpose, zones, states, and country-level geographical
aggregations. (e) Wiki groups 150 daily visits to Wikipedia articles by language and article categorical
taxonomy. (f) Favorita classifies its grocery sales by store, city, state, and country levels.
D Dataset Details
Labour:
The Labour dataset (Australian Burean of Statistics 2019
tracks monthly Australian employment
from February 1978 to December 2019, reporting total employees by part-time/full-timme status, gender, and geography. It includes N- 57 series in total, with N₂- = 25 aggregate series and N₂- 32 bottom-level series. We use 8 months from May 2019 to December 2019 as test, and the rest of the data as training and validation.
Traffic: The Traffic dataset (en Taieb & Koo 2019) contains daily (aggregated from hourly rates) freeway occupancy rates for 200 car lanes in the San Francisco Bay Arca, aggregated from Jamuary 2008 to March 2009. The data is grouped into three levels: four groups of 50 lanes, two groups of 100 lanes, and one overall group of 200 lanes, following the random grouping in Rangapuram et al. (202I]; Olivares et al. (023). Consistent with prior work (Ben Taieb & Koo 2019 Rangapuram et al. 2021 Olivares et al. 2023), we split the dataset into 120 training, 120 validation, and 126 test samples, reporting accuracy for the last test date. We use geographic node dummies, weekend indicators, and proximity to Saturday for exogenous variables.
20

---

Published in Transactions on Machine Learning Research (12/2024)
TourismLange
45
40
Tourist Visits
35
30-
20
201 1-12-01
2013-12-01
2015-12-01
Period (January 01, 2011 December 01, 2016]
Frequency [MS]
Figure 6: Tourism-L dataset partition into train, validation, and test sets used in our experiments. All datasets use the last horizon window as defined in Table 2 (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 selection and the ablation studies.
Tourism-S: The Tourism-S dataset (Tourism Australia, Canberra 2005) records quarterly visits to Australia from 1996 to 2006. We use data from 2005 for validation, 2006 for testing, and the remaining quarters for training. Each series contains 28 quarterly observations and is structured with aggregated visit data by country, purpose of travel, state, regions, and urbanization level within regions. The most disaggregated level includes 56 regions by urbanization, while the aggregated data consists of 33 series. We use state dummies, and for the future exogenous variables, we use quarterly dummics and scasonal naive 4 and 8 anchors.
The Tourism-L dataset ( (Wickramasuriya et al. 2019) represents visits to Australia, at a monthly frequency, between January 1998 and December 2016. We use 2015 for validation, and 2016 for lesting, and all previous years for training. The dataset contains 228 monthly observations. For each month, we have the number of visits each of Australia's 78 regions, which are aggregated to the zone, state, and national level, and for cach of four purposes of travel. These two dimensions of aggregation total N 304 leaf entities (a region-purpose pair), with total of M = 555 series in the hierarchy.
Tourism-L:
We pre-process the data to include static features. We use purpose of travel as well as state dummies. For the historical information, we use month dummics, and for the future exogenous variables, we use month dummies and seasonal naive anchor forecast that helps greatly to account for the series seasonality.
Favorita: The Favorita dataset (Favorita et al. 2017) contains grocery sales of the Ecuadorian Corporación Favorita in N 54 stores. We perform geographical aggregation of the sales the store, city, stale, and national levels, following OOlivares et al. 2023 This yields total of M 94 aggregates. Concerning features, past unit sales and number of transactions as historical data. In the Favorita dataset, include item perishability static information, geographic state dummy variables, and for the historic exogenous features and future exogenous features, we use promotions and day of the week. During the model's optimization we consider a balanced dataset of items and stores, for 217,944 bottom level series (4,036 items * 54 stores), along with aggregate levels for a total of 371,312 time series. The dataset is at the daily level and starts from 2013-01-01 and ends on 2017-08-15, comprehending 1688 days. We keep 34 days (1654 to 1988 days) hold-out test and 34 days (1620 to 1654 days) as validation.
Wiki: The Wikipedia dataset (Anava al. 2018 contains daily views of 145,000 online articles from July 2015 to December 2016. The dataset is processed into N₂ = 150 bottom series and N, 49 aggregate series based on country, access, agent and article categories, following Ben Taieb & Koo 2019 Rangapuram et al. (201 processing. The last week December 2016 used as test, and the remaining data as training and validation. We use day of the week dummies to capture seasonalities, country, access, agent dummy static features, and a scasonal naive anchor.
21

---

Published in Transactions on Machine Learning Research (12/2024)
E
Forecast Distributions Visualization
1.0
1.0
1.0
Empirical CDF
Emplrical CDF
0.8
Empirical
0.6
0.66
0.2
0.2
0.2
0.00
0.25
0.50
0.76
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.56
0.75
1.00
Theoretical CDF
Theoretical CDP
Theoretical CDF
(a) Labour
(b) Traffic
(c) Tourism-S
1.0
1.0
Empirical CDF
0.0
Emplrica) CDF
0.8
Empirical
0.4
0.0
0.6
0.25
Thearntical CDF
Theoretical CDF
Theomtical CDP
(d) Tourism-L
(e) Wiki
(f) Favorita
Figure 7: Crobabilliy-Probbbbii plot comparing the cumulative probabilities of the empirical distribution and the CLOVER' forecast distribution. Points near the 45-degree line indicate similarity between the distributions. The plot uses single run from the models reported in the main results Table 3
(a) Labour
(b) Traffic
Tourism-S
125
300
160
16
(d) Tourism-L
(e) Wiki
(I) Favorita
Figure Estimated correlation matrices produced by CLOVER. The color scale indicates the strength of correlation between pairs of series the hicrarchy. Green represents negative correlations, red indicates positive correlations, and black highlights the most highly correlated series. The estimated covariances are sparse, with many correlations close to zero. Aggregated levels, seen in the top-left corner, show stronger positive correlations. The plot uses single run from the models in Table 3
22

---

Published in Transactions on Machine Learning Research (12/2024)
300
Ecuador
200
100
city [Puyo] state [Pastaza]
10
15
20
25
30
10
10
20
30
10
10
20
store [29]
25 Predictions E (t 1...., t N₁
10 15
20
30
Figure CLOVER forecast distributions on the Favorita dataset. We show the forccasted demand for a grocery item on a store of the Puyo City, in the State of Pastaza and the whole country demand in the top row. Forecast distributions show the 90% forecast intervals in light blue, and the forecasted median dark blue. The clipped Normal distribution achieves non-negative predictions and a point mass zero.
23

---

Published in Transactions on Machine Learning Research (12/2024)
Table 5: Ablation study on the Traffic dataset, empirical evaluation of probabilistic coherent forecasts. Mean scaled continuous ranked probability score (sCRPS) averaged over 5 runs, at each aggregation level, the best result is highlightod (lower measurements are preferred). We report 95% confidence intervals. Comparison of different neural forecasting architectures augmented with the Factor Model and trained with CRPS.
Level
CLUVER
MQCKN
TFT
NHITS
NBEATS
LSTH
FCGAGA
Overall
0.0242++0.0035
0.0613-0.0257
0.0698++0.010
0.0644+0.0095
0.0688++0.0156
0.0646+0.0066
0.0592+00.020
Total
0.0035=+00000 0015
0.0432=0.0296
0.0311+0.0
0.0339=0.0131
0.0423++.0117
0.0377=0.0138
0.0241+0.01 01--000.103
Halves
0.0018::0000 0030
0.0137=0.02$0
0.010240.0100
0.0376=0.0122
0.0113::0 0159
10402=000005
0.0259+0 0099
Quarters
0.0041110.002
0.0440=0.0281
0.044110.0111
0.0420=0.0084
0.047810.0132
0.0402=0.0084
0.02851:0 0083
Lanes
0.0905-10.0084
0.1145
0.1656:10.0199
0.1442:0.0013
0.1409:t0.0182
0.1404#0.0071
0.1583:40.099
Table 6: Ablation study on the Traffic dataset. empirical evaluation of probabilistic coherent forecasts.
Mean scaled continuous ranked probubility score (sCRPS) averaged over 5 runs, at each aggregation level, the best result is highlighted (lower measurements are preferred). We report 95% confidence intervals. Comparison of MQCNN- based architecture trained with different learning objectives.
The Narmal and Student non coherent forecast distributions, contrast to the Factor Model and the Poisson Mixture.
Level
CLOVER Energy Score
Other Distributions
CRPS
NLL
PoissonMixture
Student
Normal
Overall
0.0259200.0 0060
0. 0.0269=0.0011
10879::001116
0,0827=0, 1 1408
0.0600 03667
0.0734=0.022 0.253
Total
0.00231100000
0.0031=0.0030
0.0623:10.1335
0.0667=0.1817
0.0280:10
0.0556::.0052
Halves
0.0028 00028++000000 +-000118
0.0035-0.0022
0.0640+0.1316
0.0763+0.2088
0.0296 +0006-00000000
0.0510+0.0016
Quarters
0.0043++000015
0.0074=000019
0.0644++0.131
0.0742200.
0.0301+:0 0294
0.0450=0.0333
Lanes
0.0912==0.0195
0.0882=0.0159
0.1608=000611
0.1136=0.0056
0.1523=0.0752
0.1422=000051
F Ablation Studies Details
To analyze the sources of improvements in our model, we conducted ablation studies on variants of the KLOVER/MOCIN/DPMI (Wen et al.) 2017 Olivares et al., (022c). We used a simplified setup on the Traffic datasot, focusing on the same forecasting task as in the main experiment. We evaluated sCRPS from Equ. 19 in the validation set across five randomly initialized neural networks. The experiments Table 4 hyperparameters, and vary a single characteristic of the network, and measure its validation effects.
In our first ablation study, we explore the effects of the impact of including vector autoregrossive relationships of the hierarchy through the MrossSeriesMLP module described in Eqn. 14 In the experiment, we train a CLOVER with and without the module in the Traffic dataset. Additionally we compare different well-
performing neural forecasting architectures" 0 augmented with the CrossSeriesMLPP including (1) LSTM (Sak
ctt al.) 2014), (2) NBEATS (Oreshkin et al.) 2020 Olivares et al. 2022a), (3) NHITS (Challu et al. 2023 (4)
TFT (Lim et al. 2021), (5) FCGAGA (Oreshkin el al. 2021) a network specialized in spatio-temporal lorecasting.
We use the default implementations available the NeuralForecast library (Olivares et al., 2022b) 7 Table
shows that convolution- based architectures continue to deliver saate-of-the-ate results., more importantly using the CrossSeriesMLP improves sCRPS upon the alternative (without) by 66 percent. The technique bridges the gap to the HierE2E (Rangapuram et al. 2021 which previously outperformed all alternative methods by over 50 percent. We attribute the improvements to the heavy presence of Granger causal relationships between the traffic lanes, as they carry lag historical information that influences each other.
In our second ablation study, we explore learning objective alternatives to Eqn. 17 For this purpose, we replace the last layer of CLDVER with different distribution outputs, including the normal, Student-t, and Poisson mixture distributions (Olivares et al. 2023). In addition we also compare with our own Factor model approach, we can see in Table 6 and Figure 3 the CRPS optimization of the Factor Model improves upon the negative log likelihood by 60 percent in the mcan sCRPS in the validation set. The difference is highly driven by outlier runs, but it is expected as the CRPS objective has much more convenient numerical properties, starting with its bounded gradients.
The Factor Model can augment models from the Neural Forecast library OOlivares et al. 2022b Because of this, it can readily augment AutoFormer, BitCN, DeepAR, DeepNPTS, DilateddNN, DLinear, FedForner, FCGaga, GRU. Informer, Transformer, KAN, LSTM, MLP, NBEATS, NBEATSx, NHITS, NLincar, PStchTST, RMOK, RNN, SOFTS, TCN, TFT, Tide, TimeLLM, TimeMixer, and TemesNet. For this ablation study, we focus on NBEATS, NHITS, LSTM, FCGaga, and TFT.
7 "Optiizzation of the neural forecasting networks follows details from Appendix we increased four times the early stopping patience on NHITS/MBEATS /FGAAA to account. for the gradient variance compared to MQCMN's gradient based forking sequences.
24

---

Published in Transactions on Machine Learning Research (12/2024)
Table 7: Empirical evaluation of mean hierarchical forecasts. Relative squared error (relSE) averaged over 5 runs, cach aggregation level, the best result is highlighted (lower values are preferred). We report 95% confidence intervals, the methods without standard deviation have deterministic solutions.
The ARTM--ER reeuts for Teuriee-t, differ from Rangaparam al.
2021
Ju we improved the numerteal stability of thair implemantscion.
CLOVER (crps)
CLOVER
ARIMA (not hier.
SWaive (not hier.))
DATA
LEVRL.
(energy)
DPHN-GroupBU
ARINA-ERM
ARTHA-MinT-ols
ARIMA-BU
Overall
0.3882400.501
0.381610.00111
0.7896=0.2822
0.5671
0.3136
0.4527
0.3289
5.0683
1abouz
(gco.)
1.14831100014
0.142111..002
0.3307=0.3012
0.0395
0.0320
0.1455
0.0369
5.9572
[goo.]
0.4965-+0. 0888
0.4889+00.18 1322
0.3206-0.4084
0.4274
0.3590
0.5927
0.3600
5.8649
(geo.)
0.66155+0.005 0528
0.6580=00.07G
1.2675=0.2674 2674
0.9561
0.6807
0.8742
0.6723
4.0696
(geo.)
0.7301=#0.0238
0.7234=0.0288
1.1693=0,3414
1.8568
0.76U5
0.8364
0.8364
2.6208
Overall
10008=00.0044
0.0021 0031
11750=0.0090
0.0199
0.0125
0.0217
0.0433
0.0709
Traffic
(geo.)
0.00011000002
0.001320.000 0.027
0.1619=0.0099
0.0133
0,0314
0.0168
0.0302
0.0547
(geo.)
0.0001
0. 0.0014+0.003
0, 0.1835=0.0011
0.0135
0,0380
0.0180
0.0392
0.0676
(geo.)
0.0005
0,0017 00335
0.181900.0100
0.0373
0.0647
0.0295
0.0850
0.0989
(geo.)
0.135422.0025
0.145410.0680
0.9964=0.0430
0. 0.6355
0.5876
0.5669
0.5669
1.3118
Overall
0.1078:0.000 01229
0.1104+0.014. 01131
0.114100.0211
0 0.3781
0.1772
0.09941
0,1886
0.2198
(geo.)
0.1044110.05 01889
0,0874100 0202
0.1362=0.0359
0.4174
0.1722
0.0696
0.1649
0.2596
Tours
(gro.)
0.0969
0.1004++0011
0.0190
0.367
0.1605
0.0819
0.2066
0.1741
(geo.)
0.1381=000094
0.1642++00099
0.1241-0.0200
0.3033
0,2015
0.1684
0819
0.2163
(geo.)
0.1570=:00077
0.1872+0.00 0206
0.1397=0.0371
0.3369
0,2111
0,2218
0,2218
0.2557
Overall
0.0951=0.0145
0.095220.011 01188
0.1113=0.0158
0.1178
0.1251
0.2979
0,1414
0.1306
(goo.)
0.044711.0071
0.040511.0332
0.0597
0.0596
0.0472
0,4002
0.0343
0.0582
Tourisn-1
[geo.)
0.1014++0.0180
0.1016+0.0300
11121-0.0122
1.93
0.76
0.3340
0.2530
0.16
(geo.)
0.2309 0124
0, 2.229990.00 0253
0.2250=0.0196
0 0.2529
0,3556
0,4238
0,1429
0.3695
(geo.)
0.3075=00.033
0226320.0296
0.2980=0.0197
0.3236
0,4288
0.4012
0,4835
0.4766
(prp.)
0.0596.100 0195
0.0560+10.007
0.0196=0.0195
0.0895
0,0856
0.1703
0.0973
0.0615
(prp.)
0.1199110.0115
0.1232+-000085
0.140000
0.1466
0.1537
0.1986
0616
0.1577
7 (prp.)
0.2484
0. 0204440.0190
0. 0.265-10002122
0.2705
0.3017
0.3151
0.2914
0.3699
(prp.)
0.3432++0.0157
0.3830
0.3302=10.0235
0.3543
0.3970
0.3769
0.3769
0.4969
Overall
0.5301
0.763910.1517
0.9270=0.1237
0.9768
1.0253
0.9621
1.0419
0,9288
(gro.)
0.1980110.0220
0.4673 9-233-0-0 0,243
0.6294=0.1769
0.8347
1.0528
0.8271
1.1357
0.6555
N1XI
2 (geo.)
00063++0.0048
0.7627 1917
0. 0.9268=0. 1056 1056
L.0342
0.9913
0.9615
0.8300
1.0672
(geo.)
0.8272
0.919220.366 3670
0.009=000997
1.0844
1.0097
1.1105
1.0398
1,1411
(geo.)
0.8152:1007322
0.933920.3577
1.1789=0.0961
1.0744
1.003U
1.0858
1.0259
1.1095
(prp.)
0.8820:10.0698
0.5834
1.1958=0.0732
1.0881
1.0294
1.0504
1.0504
1.1080
Overall
0.588510.0291
0.751740.022
0.7363-0.0713
U. 0.8163
0,9465
0.8276
0.9665
1.1420
(geo.)
0.6109=:000100
08786440.00 0291
0.79-110.00500 0.068
0. 0,8362
0,8999
0.8415
0.9217
1,1269
2 (geo.)
0.6188:0.0265
0.718320.02 0216
0.7355=0. 1057
0.7830
1.0057
0.8050
1.0451
1.1078
(geo.)
0.5819:10.0256
0.725610.022
0.7303=0.1053
0.7986
1.0418
0.8192
1.0881
1.1315
(gno.)
0.5053
0.7083+000274
0.6770-0.0357
0.8199
0.8808
0.8228
0.8228
1,2815
G Mean Forecast Accuracy Evaluation
To complement the probabilistic and L1-based results in Section 4.1
We also evaluate mean forecasts denoted by РμHM,, (F N1.),..) through the relative squared error relSE (Hyndman & Koehler 2006 that considers the ratio between squared error across forecasts in all levels over squared error of the Naive forecast (i.e., a point forecast using the last observation y11,e) as described by
No [5 el ((0)
relSE (yel)a MOJ (((0))
(25)
No NJJ Vi.e 1m 1(0)
1=
As discussed in Section 3.1 our factor model naturally defines a probabilistic coherent system, where hicrarchi- cal coherence emerges as a consequence. In this section, we compare our method with the following coherent
mean approaches: (1) OPMN-Grouphe
(2) ARIMA-ERY (Ben Taieb & Koo| 2019).
(3) ARIMA-MinT (Wickrama
suriya et al. 2019), (4) ARIMA-BU, (5) ARIMA, and (6) Scasonal Naive. The statistical methods are implemented
using the StatsForecast and HierarchicalForecaa libraries Olivares et al. 2022c Garza et al. 2022). In
particular we use available code on the hierarchical baselines repository and hierarchical datasets repository
The CLOVER model improves relSE accuracy in five out of six datascts by an average of 32%, while degrading accuracy on the smallest dataset, Labour, by 32%. These changes can be attributed to squared error metrics' sensitivity to outliers, though the results are generally consistent with sCRPS outcomes.
25

---

Published in Transactions on Machine Learning Research (12/2024)
H
Efficient Estimation of Scoring Rules
CRPS samples
CRPS 10 samples
100 samples
samples
CRPS Train Loss
Analytical CRPS
10a
400
Gradient Steps
(a) Training Trajectory
CRPS 5 samples
CRPS 10 samples
SCRPS Validation Loss
CRPS 100 samples
CRPS 1000 samples
Analytical CRPS
10-1
200
400
600
800
1000
Gradient Steps
(b) Validation Trajectory
Figure 10: Training and validation Continuous Ranked Probability Score curves on the SF Bay Area Traffic
dataset. We show curves for CLOVER varying as a function of the gradient steps. Fach curve uses a different. number of Monte Carlo samples estimate CRPS.
For the network's optimization we implement the following CRPS Monte Carlo estimator for Eqn. 5
N
N - 35
CRPS (y. Y()) =
y; yl
Σ
(26)
N
2N(N
Where ye Uj
= = 1,..,, are independent samples of the random variable Y()), and a possible observation.
To assess the effect of sample size on the CRPS Monte Carlo estimator, we ran diagnostics throughout the model training process. Figure 10 shows the training and validation loss trajectories for varying Monte Carlo sample sizes, using model configurations from Table 4 For reference, we include the analytical CRPS of
truncated normals (Thorarinsdottir & Gneiting 2009 despite minor distributional differences compared to Eqn.
10 For comprehensive review of alternative CRPS estimators, we refer to Zamo & Naveau 2018
We found that reducing the number of Monte Carlo samples had minimal impact on both training and validation scores. Given GPU memory constraints, decreasing the sample size effectively reduces computational overhead without significantly affecting model performance. We observed similar phenomena with the energy score learning objective and its relationship to Monte Carlo samples.
26

---

Published in Transactions on Machine Learning Research (12/2024)
Note on the Reparameterization Trick
In this section, we discuss the importance of the reparametcrization trick in the context of hicrarchical scries and our multivariate factor model. We explain how it enables us to maintain the network's diflerentiability while optimizing arbitrary learning objectives, regardless of the complexity of our probabilistic model.
Let disaggregated scries be
Y(0) - u(0) + Diag(o F(O)E RR.
(27)
Let the hierarchical series under aggregation constraints be
(0))= = b (0)++
(28)
with e ~ (0, IN.)) and z ~ (O, IN,) the reparameterization random variables. For simplicity of exposition,
we consider a single entry in the reconciled multivariate series (0) that we denote Y(0) = S (0)
For a differentiable learning objective L IR + R, the loss gradient with respect to the network parameters 0 can be obtained as
VeEy [C((()))] Vo I,, C(Y
J VoL(((())P(Y))))) / [(Y())))
(29)
= Ey (Vo(( +
Note that the expectation's gradient is not the expectation of the gradient due to the second term that captures the dependence of the probability of (0) on the parameters of the network. The second term may face challenges when the probability of Y()) does not have an analytical form or its gradient VYP((Y poses challenges for automatic differentiation tools. The solution proposed by Kingma & Welling 2013 reparamcterizes the target variable to recover the differentiability of the gradient loss in Equ. 29
Using the chain rule and the reparameterization trick, the loss' gradient can be recovered:
Vo (F((()))
VolE,,
4 (S (
((0))
Ez,e
VoL
8L
(Ver (2)(())) (0))
(30)
Ez,e
8Y
8L
Ez,
S..6
(Vop Voo VoFe]]
8Y
As highlighted in Section 5 the key insight from Eqn. 30 is that the reparameterizatioon trick provides a highly floxible framework. It is compatible with any differentiable learning objective and can accommodato differentiable constraints that extend beyond traditional aggregation constraints, such as those Eqn. 28 Furthermore, the techniques explored in this paper can be generalized continuous distributions (Figurnov et al. 2018 Ruiz al. 2016 Jankowiak & Obermeyer 2018 beyond Gaussian random variables, opening up exciting avenues for future research.
27


Reviews:
[{"review_id": "X2wZJKmf0p", "reviewer": "Reviewer_RBzB", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* The proposed method is novel and interesting.\n\n* The proposed approach shows promising performance.\n\n* The paper is generally well written.\n\nWeakness:\n* The paper only considers aggregate time-series which are linear combinations of the base time-series (Eq. 1). Can the proposed approach deal with time-series which are composed on non-linear aggregates of the base time-series?\n\n* The paper says that the HierE2E method of Rangapuram et al. (2021) is “too general\". Could the HierE2E model be better regularized with more data? How does the proposed approach compare to HierE2E  with increasing amounts of available training data?\n\n* Real-world datasets: The paper only considers small scale datasets such as Traffic. It would be interesting to see if the proposed approach scales to larger datasets such as nuScenes.\n\n* The paper could consider simpler transformer based approaches as baselines.\n\n* Few typos: Eq. 10.\n\nOverall:\n* A great paper with a novel approach and good results.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The focus of the paper is to perform hieratical probabilistic forecasting of time series. To take advantage of the hierarchy between the forecasted quantities, the proposed approach takes a MQForecaster neural network architecture and augments it with a  deep Gaussian factor forecasting model. The proposed approach achieves promising performance when compared to the state of the art.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See weaknesses above.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "z26K8ylWBW", "reviewer": "Reviewer_3wgV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "My general opinion of the paper is positive, I think that the topic is interesting and the claims are likely true given the empirical results. However, there are some concerns mostly with writing and experimental design. If enough of my concerns are addressed I would lean towards accepting the paper.\n\nStrengths:\n- Interesting problem, relevant to TMLR community.\n- The proposed architecture does what it promises, it models hierarchical time series and aggregates them coherently by design.\n- Strong experimental results. The metrics often improve above 5%.\n- Notation and writing are clear, though there are occasional typos (see Rephrase / Typos section in \"Requested Changes\" part of this review). Note this statement is qualified by the fact that I have a single strong objection to the writing in the introduction section. Please see Additional questions under \"Requested Changes\" part of the review. All other parts of the paper are clear and well written.\n\nWeaknesses:\n- Limited datasets: It seems that the datasets that the authors compare on are only a subset of the datasets explored from the baselines. Why is this? \n- Selection of forecast metrics: CRPS is an excellent metric to report. However, I have my concerns about the usefulness of reporting rSE. In general, squared error based metrics can overestimate the importance of outliers. Just because a method does better wrt to outliers does not mean that the method is better in general. Please see requested change for additional metrics.\n- Baseline selection: It seems that a new model \"HiPerformer\" [https://arxiv.org/pdf/2305.08073] compares on similar datasets and baselines.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes a new method for hierarchical time series forecasting. The paper makes the following contributions:\n- Proposal of Deep Coherent Factor Model Neural Network (DeepCoFactor) model, a method for probabilistic coherent forecasting. DeepCoFactor can directly optimize metrics of interest such as CRPS, MSE, or others rather than naively maximizing the log likelihood which might not always provide the best generalization under model misspecification. This makes the proposed approach versatile.\n- The paper has basic grounding in theory (Lemma A.2) and offers strong experimental results on 3 datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "## Dataset Choice\nCan the authors justify the dataset choice or provide experiments with more datasets? Even a controlled toy experiment might prove useful to explain where the performance improvement of their method comes from and what generalizable insights we can draw about model performance beyond the selected datasets. Some dataset ideas come from the baselines, e.g., in HIER-E2E (see paper for full bibliographical citation): Labour, Tourism, Wiki.\n\n## Baseline selection\nCan the authors explain why they do not compare to \"HiPerformer\" [https://arxiv.org/pdf/2305.08073]? Is a comparison not applicable? If it is can they revise the paper to compare to this approach?\n\n## Additional Metrics:\nAs mentioned in weaknesses section, I have my concerns about reporting rSE. In general, squared error based metrics can overestimate the importance of outliers. Just because a method does better wrt to outliers does not mean that the method is better overall. A very useful metric that is less prone to noise from outliers is symmetric mean absolute percent error sMAPE. This metric would 100% convince me that the proposed method is indeed better, at least for point predictions. For distributional predictions, CRPS is a good metric and I am convinced that the proposed method is strong here. An nice-to-have metric to check would be 95% C.I. coverage (how often is the predicted interval between the 2.5th and 97.5th percentile covering the ground truth)? How does this compare to baselines?\n\n## Notation\nNotation is generally a strong point of the paper. However I have one question.\n- Definition 1.1: \"[...] for any set \\mathcal{B} \\in \\mathbb{F}_{[b]}\". What is \\mathbb{F}_{[b]}? My understanding is that \\mathbb{P} maps elements of the sigma algebra to real valued scalar / vector. Should the notation be \\mathcal{B} \\in \\mathcal{F}_{[b]} instead?\n\n## Rephrase / Typos:\n- Some phrases that need editing: \"DeepCoFactor leverages coherently aggregates the samples of [...]\".\n- Some typos, e.g., in first bullet of contributions, last sentence: \"[...] specializes in mult-step probabilistic forecast\". \nPlease edit the text to fix these.\n\n## Additional questions:\n- In Table 1 the authors claim that baseline method HIER-E2E (see paper for full bibliographical citation) is does not model a joint multivariate probability distribution. However, HIER-E2E claims to be a multivariate method. Is something unclear to me or do these two statements contradict?\n- The paper borrows a lot of definitions and examples at the end of the introduction from DPMN baseline (see paper for full bibliographical citation). Figures 1 and 2 are very similar to DPMN as well. While I do not consider this plagiarism, if definitions and the example of Fig. 1 are directly taken from DPMN it should be made clear that the definitions are being quoted from that paper with proper attribution.\nFor example, \"We can formalize this through Definition 1.1, which essentially tells us that the distribution of a given aggregate random variable is exactly the distribution defined as [...]\". I would rephrase this sentence to point to the fact that definitions are taken from the literature as it might mislead the audience into thinking that this definition is original. Additionally, to make the paper more novel I would generate I completely new example in figure 1 as it is basically identical to DPMN example. Furthermore, I would explain that Fig. 2 mirrors style of DPMN architecture figure to highlight differences between the model inputs, internals and outputs. Overall the similarity in Fig. 2 makes the novelty of this paper seem smaller than it actually is.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ME5GmvaIGY", "reviewer": "Reviewer_jJ2e", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The (long) introduction  is well-written\n- Methods are clearly decribed, down to hyperparameters and including some well-formated code.\n\nWeaknesses\n- The rest of the paper, apart from the introduction but including the abstract, could use some polish in writing style. While is did not find too many errors, it was hard for me to follow sometimes. Additional points see below.\n- The experiments do not show everything that is being claimed:\n  - 3rd contribution: while you show better accuracy than the state of the art, you do not demonstrate \"the adaptability to fit complex temporal relationships between known data and future predictions, while also modeling correlations between elements in the hierarchy\" later in the paper\n  - I think Table 3 and Figure 3.a (and more tables and figures later) are misleading. You validate on CRPS. Of course, this validation works best, if you train on CRPS and not on different criteria. (I know that e.g. Table 4 shows additional results that do not have this problem.)\n  - I do not understand how much of the gain by this method comes from the loss function and how much from certain parts of the model architecture. While you make some ablations, that part could be clearer (and be described more clearly). Why using Convolutions? Why \n- Just a combination of known techniques, less fancy (but ok for TMLR)", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a probabilistic model for hierarchical time series prediction. The model uses a parametrization matrix to predict all (linear) aggregations from the base predictions. Additionally, a new network structure based on temporal convolutions is introduced and the classical, but not well-known loss function CRPS is being used. Experiments are conducted that show the advantage of this approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "You should clearly improve your experiments!\n\n- In general, the introduction seems excessively long with 3.5 pages. You might split everything after your contributions into a new section. Or you cold write the introduction shorter.\n- 2nd contribution: aren't these multiple contributions combined in one long paragraph?\n- Is there no joint probability between time steps in (3)? I.e. are your predictions conditionally independent predictions after conditioning on the past?\n- While notations (p.3) is mostly clear, somewhat more context would help, like writing some more invariants like $[i] = [a]\\uplus [b]$ or $|[a]|=N_a$ etc.\n- What is a probabilistic forecast space? Is this a probability space (as vagely defined in Definition 1.1)? (If it is the latter case, one usually writes \"(Omega,F,P)\" instead of \"Omega,F,P\".) Since probabilistic forecast space appears again later, it should be formally be defined.\n- Definition 1.1: What is a \"set’s image\"? Under which map? The projection on the relevant coordinates?\n- Please introduce CRPS and clearly motivate its main idea and advantage.\n- Various tables: how many standard deviations are you giving in your error estimations?\n- Table 3: regarding the missing results: can you extend the implementation, that prevented these results? Or can you produce the missing results for at least on hierarchy?\n- I do not understand Figure 3. Which data sets are being used? How many repetitions?\n- In Section 3.3, please directly reference the tables and figures with the corresponding results and include these results in the main paper. (You only reference to APpendix F for details, not results.) I find this very confusing.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a novel approach that is based on interesting ideas. The methodology experiments are sufficient to support the claims.

The only minor change that is requested is a moderation of the claim in the introduction that the approach can be “based on any generic deep learning univariate forecasting model”. This should at least be moderated by the condition that the forecasting model can output factor model parameters. But in addition to this, the experiments are conducted with only a single forecasting model, so there is not really empirical evidence that the method can successfully employ any univariate forecasting model.  While “successfully” is not claimed, it is implied by the current text. 

In making their recommendations, all reviewers acknowledged that the method is novel and stated that it would be of interest to the research community. One reviewer commented that further experiments on additional datasets would further solidify the work.
Decision: Accept (poster)