HIERARCHICALLY REGULARIZED DEEP FORECASTING
Year: 2021

Abstract:
Hierarchical forecasting is a key problem in many practical multivariate forecasting applications the goal is to simultaneously predict a large number of correlated time series that are arranged in a pre-specified aggregation hierarchy. The main challenge is to exploit the hierarchical correlations to simultaneously obtain good prediction accuracy for time series at different levels of the hierarchy. In this paper, we propose a new approach for hierarchical forecasting which consists of two components. First, decomposing the time series along a global set of basis time series and modeling hierarchical constraints using the coefficients of the basis decomposition. And second, using a linear autoregressive model with coefficients that vary with time. Unlike past methods, our approach is scalable (inference for a specific time series only needs access to its own history) while also modeling the hierarchical structure via (approximate) coherence constraints among the time series forecasts. We experiment on several public datasets and demonstrate significantly improved overall performance on forecasts at different levels of the hierarchy, compared to existing state-of-the-art hierarchical models.

heading: 1 INTRODUCTION
text: Multivariate time series forecasting is a key problem in many domains such as retail demand forecasting (Böse et al., 2017), financial predictions (Zhou et al., 2020), power grid optimization (Hyndman & Fan, 2009), road traffic modeling (Li et al., 2017), and online ads optimization (Cui et al., 2011). In many of these setting, the problem involves simultaneously forecasting a large number of possibly correlated time series for various downstream applications. In the retail domain, the time series may capture sales of items in a product inventory, and in power grids, the time series may correspond to energy consumption in a household. Often, these time series are arranged in a natural multi-level hierarchy - for example in retail forecasting, items are grouped into subcategories and categories, and arranged in a product taxonomy. In the case of power consumption forecasting, individual households are grouped into neighborhoods, counties, and cities. The hierarchical structure among the time series can usually be represented as a tree, with the leaf nodes corresponding to time series at the finest granularity, and the edges representing parent-child relationships. Figure 1 illustrates a typical hierarchy in the retail forecasting domain for time series of product sales.
In such settings, it is often required to obtain good forecasts, not just for the leaf level time-series (fine grained forecasts), but also for the aggregated time-series corresponding to higher level nodes (coarse gained forecasts). Furthermore, for interpretability and business decision making purposes, it is often desirable to obtain predictions that are roughly coherent or consistent (Hyndman et al., 2011) with respect to the hierarchy tree. This means that the predictions for each parent time-series is equal to the sum of the predictions for its children time-series. More importantly, incorporating coherence constraints in a hierarchical forecasting model captures the natural inductive bias in most hierarchical datasets, where the ground truth parent and children time series indeed adhere to additive constraints. For example, total sales of a product category is equal to the sum of sales of all items in that category.
Some standard approaches for hierarchical forecasting include bottom-up aggregation, or reconciliation-based approaches. Bottom-Up aggregation involves training a model to obtain predictions for the leaf nodes, and then aggregate up along the hierarchy tree to obtain predictions for higher-level nodes. Reconciliation methods (Ben Taieb & Koo, 2019; Taieb et al., 2017; Van Erven & Cugliari, 2015; Hyndman et al., 2016; Wickramasuriya et al., 2015; 2020; Panagiotelis et al., 2020) make use of a trained model to obtain predictions for all nodes of the tree, and then, in a separate post-processing phase, reconcile (or modify) them using various optimization formulations to obtain coherent predictions. Both of these approaches suffer from shortcomings in term of either aggregating noise as one moves to higher level predictions (bottom-up aggregation), or not jointly optimizing the forecasting predictions along with the coherence constraints (for instance, reconciliation).
At the same time, there have been several recent advances on using Deep Neural Network models for multivariate forecasting, including Recurrent Neural Network (RNN), Convolutional Neural Network (CNN) architectures (Salinas et al., 2020; Oreshkin et al., 2019; Rangapuram et al., 2018; Benidis et al., 2020), and models designed for multivariate time series based on dimensionality reduction techniques (Sen et al., 2019; Wang et al., 2019; Nguyen & Quanz, 2021; Salinas et al., 2019; Rasul et al., 2020; de Bézenac et al., 2020), that have been shown to outperform classical time-series models such as autoregressive and exponential smoothing models (McKenzie, 1984; Hyndman et al., 2008; Hyndman & Athanasopoulos, 2018), especially for large datasets. However, most of these approaches do not explicitly address the question of how to model the hierarchical relationships in the dataset. Deep forecasting models based on Graph Neural Networks (GNN) (Bai et al., 2020; Cao et al., 2020; Yu et al., 2017; Li et al., 2017; Wu et al., 2020) do offer a general framework for learning on graph-structured data. However it is well known (Bojchevski et al., 2020) that GNNs are hard to scale for learning on graphs with a very large number of nodes, which in real-world settings such as retail forecasting, could involve hundreds of thousands of time series. More importantly, a desirable practical feature for multi-variate forecasting models is to let the prediction of future values for a particular time series only require as input historical data from that time series (along with covariates), without requiring access to historical data from all other time series in the hierarchy. This allows for scalable training and inference of such models using mini-batch gradient descent, without requiring each batch to contain all the time series in the hierarchy. This is often not possible for GNN-based forecasting models, which require batch sizes of the order of the number of time series.
Problem Statement: Based on the above motivations, our goal is to design a hierarchical forecasting model with the following requirements: 1) The model can be trained using a single-stage pipeline on all the time series data, without any separate post-processing, 2) The model captures the additive coherence constraints along the edges of the hierarchy, 3) The model is efficiently trainable on large datasets, without requiring, for instance, batch sizes that scale with the number of time series.
We propose a principled methodology to address all these above requirements for hierarchical forecasting. Our model comprises of two components, both of which can support coherence constraints. The first component is purely a function of the historical values of a time series, without distinguishing between the individual time series themselves in any other way. Coherence constraints on such a model correspond to imposing an additivity property on the prediction function - which constrains the model to be a linear autoregressive model. However, crucially, our model uses timevarying autoregressive coefficients that can themselves be nonlinear functions of the timestamp and other global features (linear versions of time-varying AR have been historically used to deal with non-stationary signals (Sharman & Friedlander, 1984)). We will refer to this component as the time-varying autoregressive model.
The second component focuses on modeling the global temporal patterns in the dataset through identifying a small set of temporal global basis functions. The basis time-series, when combined in different ways, can express the individual dynamics of each time series. In our model, the basis time-series are encoded in a trained seq-2-seq model (Sutskever et al., 2014) model in a functional form. Each time series is then associated with a learned embedding vector that specifies the weights for decomposition along these basis functions. Predicting a time series into the future using this model then just involves extrapolating the global basis functions and combining them using its weight vector, without explicitly using the past values of that time series. The coherence constraints therefore only impose constraints on the embedding vectors of each time series, which can be easily modeled by a hierarchical regularization function. We call this component a basis decomposition model. As we will see, this part of the model is only approximately coherent unless the embedding constraints hold exactly. In particular, in this paper, we focus on improving model accuracy rather than preserving
exact coherency. In Section A.2, we also provide theoretical justification for how such hierarchical regularization using basis decomposition results in improved prediction accuracy.
We experimentally evaluate our model on multiple publicly available hierarchical forecasting datasets. We compare our approach to state-of-the-art (non-hierarchical) deep forecasting models, GNN-based models and reconciliation models, and show that our approach can obtain consistently more accurate predictions at all levels of the hierarchy tree.

heading: 2 RELATED WORK ON DEEP HIERARCHICAL MODELS
text: In addition to the works referenced in the previous section, we now discuss a few papers that are more relevant to our approach. Specifically, we discuss some recent deep hierarchical forecasting methods that do not require a post-processing reconciliation step. Hierarchical forecasting methods can be roughly divided into two categories: point forecasters and probabilistic forecasters. Mishchenko et al. (2019) propose a point-forecasting approach which imposes coherency on a base model via `2 regularization on the predictions. Gleason (2020) extend the idea further to impose the hierarchy on an embedding space rather than the predictions directly. SHARQ (Han et al., 2021) follows a similar `2 regularization based approach as Mishchenko et al. (2019), and also extends the idea to probabilistic forecasting. Their model is trained separately for each of the hierarchical levels starting from the leaf level, thus requiring a separate prediction model for each level.
Probabilistic forecasting methods include, Hier-E2E (Rangapuram et al., 2021) which produces perfectly coherent forecasts by using a projection operation on base predictions from a DeepVAR model (Salinas et al., 2019). It requires the whole hierarchy of time series to be fed as input to the model leading to a large number of parameters, and hence does not scale well to large hierarchies. Yanchenko et al. (2021) take a fully Bayesian approach by modelling the hierarchy using conditional distributions.

heading: 3 PROBLEM SETTING
text: We are given a set of N coherent time series of length T , arranged in a pre-defined hierarchy consisting of N nodes. At time step t, the time series data can be represented as a vector yt ∈ RN denoting the time series values of all N nodes. We compactly denote the set of time series for all T steps as a matrix Y = [y1, · · · ,yT ]> ∈ RT×N . Also define y(i) as the ith column of the matrix Y denoting all time steps of the i th time series, and y(i)t as the t th value of the i th time series. We compactly denote the H-step history of Y by YH = [yt−H , · · · ,yt−1]> ∈ RH×N and the H-step history of y(i) by y(i)H = [y (i) t−H , · · · ,y (i) t−1] ∈ RH . Similarly define the F -step future of Y as YF = [yt, · · · ,yt+F−1]> ∈ RF×N . We use the ·̂ notation to denote predicted values, for example ŶF , ŷF and ŷt.
Time series forecasts can often be improved by using features as input to the model along with historical time series. The features often evolve with time, for example, categorical features such as type of holiday, or continuous features such as time of the day. We denote the matrix of such features byX ∈ RT×D, where the t th row denotes the D-dimensional feature vector at the t time step. For simplicity, we assume that the features are global, meaning that they are shared across all time series. We similarly defineXH andXF as above.
Hierarchically Coherent Time Series. We assume that the time series data are coherent, that is, they satisfy the sum constraints over the hierarchy. The time series at each node of the hierarchy is the equal to the sum of the time series of its children, or equivalently, equal to the sum of the leaf time series of the sub-tree rooted at that node. Figure 1 shows an example of a sub-tree rooted at a node.
As a result of aggregation, the data can have widely varying scales with the values at higher level nodes being magnitudes higher than the leaf level nodes. It is well known that neural network training is more efficient if the data are similarly scaled. Hence, in this paper, we work with rescaled time series data. The time series at each node is downscaled by the number of leaves in the sub-tree rooted at the node, so that now they satisfy mean constraints rather than sum constraints described above. Denote by L(p), the set of leaf nodes of the sub-tree rooted at p. Hierarchically coherent data satisfy
the following data mean property,
y(p) = 1 |L(p)| ∑ i∈L(p) y(i) (Data Mean Property). (1)

heading: 4 HIERARCHICALLY REGULARIZED DEEP FORECASTING - HIRED
text: We now introduce the two components in our model, namely the time-varying AR model and the basis decomposition model. As mentioned in the introduction a combination of these two components satisfy the three requirements in our problem statement. In particular, we shall see that the coherence property plays a central role in both the components. For simplicity, in this section, all our equations will be for forecasting one step into the future (F = 1), even though all the ideas trivially extend to multi-step forecasting. The defining equation of our model can be written as,
ŷ (i) F = f(y (i) H ,XH,XF ,ZH, θi) = 〈 y
(i) H , a(XH,XF ,ZH) 〉︸ ︷︷ ︸ Time varying AR (TVAR) + 〈 θi, b(XH,XF ,ZH) 〉︸ ︷︷ ︸ Basis decomposition (BD) . (2)
In the above equation, ZH is a latent state vector that contains some summary temporal information about the whole dataset, and θi is the embedding/weight vector for time-series i in the basis decomposition model. ZH can be a relatively low-dimensional temporally evolving variable that represents some information about the global state of the dataset at a particular time. We use the Non-Negative Matrix Factorization (NMF) algorithm by Gillis & Vavasis (2013) to select a small set of representative time series that encode the global state. If the indices of the selected representative time-series is denoted by {i1, · · · , iR} (R denotes the rank of the factorization), then we define Z = [Y (i1), · · · ,Y (iR)] ∈ RT×R. Note that we only feed the past values ZH as input to the model, since future values are not available during forecasting. Also, note that the final basis time-series can be a non-linear function of ZH. In our experiments, R is tuned but is always much much less than N . a and b are functions not dependent on y(i)H and θi. We will provide more details as we delve into individual components.
Time-Varying AR (TVAR): The first part of the expression in equation 2 denoted by Time Varying AR (TVAR) resembles a linear auto-regressive model with coefficients a(XH,XF ,ZH) ∈ RH , that are a function of the input features, and thus can change with time. The AR parameters of this model are shared across all time series and hence do not encode any time-series specific information, a drawback that is overcome by the Basis Decomposition part of our model. This component is coherent by design because it is a shared linear AR model. However, even though the AR weights are shared across all the time-series at a given time-point, they can crucially change with time, thus lending more flexibility to the model.
Implementation: In order to model the sequential nature of the data, we use an LSTM encoder to encode the past XH and ZH. Then, we use a fully connected (FC) decoder for predicting the auto-regressive weights. Similar to Wen et al. (2017)’s multi-horizon approach, we use a different head of the decoder for each future time step resulting in a F -headed decoder producing F -step predictions for TVAR weights. The decoder also takes as input the future covariatesXF if available. The produced weights are then multiplied (inner product) to the history to produce the final TVAR predictions. We illustrate this architecture in Figure 2 (right).
Basis Decomposition (BD) with Hierarchical Regularization: Now we come to the second part of our model in equation 2. As discussed before, this part of the model has per time-series adaptivity, as different time-series can have different embeddings. It resembles an expansion of the time series on a set of basis functions b(XH,XF ,ZH) ∈ RK , with the basis weights/embedding for time series i denoted by θi ∈ RK . Both the basis functions and the time series specific weights are learned from the data, rather than fixing a specific form such as Fourier or Wavelet basis.
The idea of using a basis has also been recently invoked in the time-series literature (Sen et al., 2019; Wang et al., 2019). The basis recovered in the implementation of Wang et al. (2019) is allowed to vary for each individual time-series and therefore is not a true basis. Sen et al. (2019) do explicitly recover an approximate basis in the training set through low-rank matrix factorization regularized
by a deep global predictive model alternatingly trained on the basis vectors, thus not amenable to end-to-end optimization. We shall see that our model can be trained in an end-to-end manner.
Embedding Regularization for Approximate Coherency: The TVAR part of our model is coherent by design due to its linearity. The BD model however requires the embeddings of the time-series to satisfy the mean property along the hierarchy. This directly translates to coherency of the predictions due to linearity with respect to θ.
θp = 1 |L(p)| ∑ i∈L(p) θi (Embedding Mean Property), (3)
We impose this constraint approximately via an `2 regularization on the embedding.
Ereg(θ) = N∑ p=1 ∑ i∈L(p) ‖θp − θi‖22. (4)
The purpose of this regularizer is two fold. Firstly, we observe that, when the leaf embeddings are kept fixed, the regularizer is minimized when the embeddings satisfy the mean property (3), thus encouraging coherency in the predictions. Secondly, it also encodes the inductive bias present in the data corresponding to the hierarchical additive constraints. We provide some theoretical justification for this hierarchical regularization in Section 5.
Implementation: As before, we use an LSTM encoder to encode the pastXH and ZH. Then, we use the encoding from the encoder along with the future featuresXF (sequential in nature) and pass them through an LSTM decoder to yield the F -step basis predictions which are then multiplied with the embedding (inner product) to produce the final BD predictions. We illustrate this architecture in the top right of Figure 2. Thus, a functional representation of the basis time-series is implicitly maintained within the trained weights of the basis generating seq-2-seq model. Note that the embeddings are also trained in our end-to-end model. We illustrate this architecture in Figure 2 (left).
We emphasize that the main ideas in our model are agnostic to the specific type of neural network architecture used. For our experiments, we specifically use an LSTM architecture (Hochreiter & Schmidhuber, 1997) for the encoder and decoder. Other types of architectures including transformers (Vaswani et al., 2017) and temporal convolution networks (Borovykh et al., 2017) can also be used.
Loss Function: During training, we minimize the mean absolute error (MAE) of the predictions along with the embedding regularization term introduced above (our method generalizes to other losses too, such as mean square error, or mean absolute percentage error). For regularization weight λE , and ŷ (i) F defined as Eq (2), and Θ denoting the trainable parameters of a, b, our training loss function is, `(Θ,θ) =
∑ i ∑ F |y(i)F − ŷ
(i) F |︸ ︷︷ ︸
Prediction loss
+ λEEreg(θ)︸ ︷︷ ︸ Embedding regularization . (5)
Note that the time-series dependent part of the loss function can be easily mini-batched and the embeddings are not memory intensive.

heading: 5 THEORETICAL JUSTIFICATION FOR HIERARCHICAL MODELING
text: In this section, we theoretically analyze the benefits of modeling hierarchical constraints in a much simplified setting, and show how it can result in provably improved accuracy, under some assumptions. Since analyzing our actual deep non-linear model for an arbitrary hierarchical set of time series can be complex, we make some simplifying assumptions to the problem and model. We assume that all the time series in the dataset is a linear combination of a small set of basis time series. That is, Y = Bθ + w, where B ∈ RT×K denotes the set of basis vectors, θ = [θ1, · · · , θN ] ∈ RK×N denotes the set of weight vectors used in the linear combination for each time series, andw ∈ RT×N denotes the noise matrix sampled i.i.d as w ∼ N (0, σ2) for the leaf nodes. A classical example of such a basis set can be a small subset of Fourier or Wavelet basis (Strang, 1993; van den Oord et al., 2016) that is relevant to the dataset. Note that we ignore the TVAR model for the sake of analysis and focus mainly on the BD model which includes the hierarchical regularization.
Algorithm 1: Basis Recovery Input: Observed y, basis dict B̄, regularization parameter λL Output: Estimated basisB α̂0 ← argmin α∈Rn 1 2T ‖y0 − B̄α‖ 2 2 + λL‖α‖1 Estimate support Ŝ = {i | |α̂0| > 0} Estimate true basisB ← B̄Ŝ
Algorithm 2: Parameter Recovery Input: Observed time series y, estimated basisB, regularization parameter λE Output: Estimated parameters θ θ̂0 ← argminθ0 1 T ‖y0 −Bθ0‖ 2 2 for n ∈ L(0) do θ̂n ← argminθn 1 T ‖yn−Bθn‖ 2 2+λE‖θ̂0−θn‖22. end
In this section, we consider a 2-level hierarchy of time-series, consisting of a single root node (indexed by 0) with L children (denoted by L(0)). We will also assume that instead of learning the K basis vectors B from scratch, the K basis vectors are assumed to come from a much larger dictionary B̄ ∈ RT×D of D ( K) vectors that is fixed and known to the model. While the original problem learns the basis and the coefficients θ simultaneously, in this case the goal is to select the basis from among a larger dictionary, and learn the coefficients θ .
We analyze this problem, and show that under the reasonable assumption of the parent embedding θ0 being close to all the children embeddings θn, using the hierarchical constraints can result in a meansquare error at the leaf nodes that is a multiplicative factor L smaller than the optimal mean-square error of any model that does not use the hierarchical constraints. Our proposed HIRED model, when applied in this setting would result in the following (hierarchically) regularized regression problem:
min θ
1
NT ‖y −Bθ‖22 + λ ∑ n∈L(0) ‖θ0 − θn‖22. (6)
For the sake of analysis, we instead consider a two-stage version, described in Algorithm 1 and Algorithm 2: we first recover the support of the basis using Basis Pursuit (Chen et al., 2001). We then estimate the parameters of the root node, which is then plugged-in to solve for the parameters of the children node. We also define the baseline (unregularized) optimization problem for the leaf nodes that does not use any hierarchical information, as
θ̃n = argmin θn
1 T ‖yn −Bθn‖22 ∀n ∈ L(0). (7)
The basis support recovery follows from standard analysis (Wainwright, 2009) detailed in Lemma 1 in the Appendix. We focus on the performance of Algorithm 2 here. The following theorem bounds the error of the unregularized (θ̃n) and the hierarchically-regularized (θ̂n, see Algorithm 2) optimization solutions. A proof of the theorem can be found in Appendix A.2. Theorem 1. Suppose the rows ofB are norm bounded as ‖Bi‖2 ≤ r, and ‖θn − θ0‖2 ≤ β. Define Σ = BTB/T as the empirical covariance matrix. For λE = σ 2K Tβ2 , θ̃n and θ̂n can be bounded as,
E‖θ̃n − θn‖2Σ ≤ σ2K
T , E‖θ̂n − θn‖2Σ ≤ 3
σ2K
T
1
1 + σ 2K
Tr2β2
+ 6 σ2K
TL . (8)
In fixed design linear regression ‖θ̂n − θn‖2Σ = ‖B(θ̂n − θn)‖2 is the population squared error (see Appendix A.5 for a bound on the parameter estimation error). The gains due to the regularization can be understood by considering the case when β is upper bounded by a sufficiently small quantity. Note that an upper bound on β essentially implies that the children time-series have structural similarities as further elaborated in Appendix A.5. We show that the above assumption yields a smaller upper bound on the error. In fact, if β = o( √ K/T ), then the numerator 1 + σ
2K Tr2β2 in Eq. (8) is ω(1)
resulting in E‖θ̂n− θn‖2Σ = o(σ 2K T ) which decays faster than σ2K T . Furthermore, if β is even smaller as β = O( √ K/LT ), then following similar calculations, E‖θ̂n − θn‖2Σ = O(σ 2K LT ) which is again smaller than the unregularized bound.

heading: 6 EXPERIMENTS
text: We implemented our proposed model in Tensorflow (Abadi et al., 2016) and compared against multiple baselines on popular hierarchical time-series datasets.
Datasets. We experimented with three hierarchical forecasting datasets - Two retail forecasting datasets, M5 (M5, 2020) and Favorita (Favorita, 2017); and the Tourism (Tourism, 2019) dataset consisting of tourist count data. The history length and forecast horizon (H,F ) were set to (28, 7), (28, 7) and (24, 4), for Favorita, M5 and Tourism respectively. More information can be found in Appendix B.2. We divide each of the datasets into training, validation and test sets, with details on the splits provided in Appendix B.3.
Baselines. We compare our proposed approach HIRED with the following baseline models: (i) RNN - we use a seq-2-seq model shared across all the time series, (ii) DeepGLO (Sen et al., 2019), (iii) DCRNN (Li et al., 2017), a GNN based approach where we feed the hierarchy tree as the input graph, (iv) Deep Factors (DF) (Wang et al., 2019), (v) L2Emb (Gleason, 2020), which is an improvement over Mishchenko et al. (2019), (vi) SHARQ (Han et al., 2021), (vii) HierE2E (Rangapuram et al., 2021) - the standard implementation produces probabilistic forecasts, we however adapt it to point forecasts by using their projection step on top of a seq-2-seq model. We use code publicly released by the authors for DeepGLO1 and DCRNN2. We implemented our own version of DeepFactors, HierE2E, and L2Emb for a fair comparison, since the official implementations either make rolling probabilistic forecasts, or use a different set of covariates. Additionally, we also compare with the recent ERM (Ben Taieb & Koo, 2019) reconciliation method applied to the base forecasts from the RNN model, denoted as (vii) RNN+ERM. It has been shown in (Ben Taieb & Koo, 2019) to outperform many previous reconciliation techniques such as MinT (Wickramasuriya et al., 2019). For a fair comparison, we use the same Mean Absolute Error (MAE) loss for all the compared methods, and make sure that all the models have approximately the same number of parameters. Further details about the baselines and training parameters can be found in Appendix B.
Metrics. We compare the accuracy of the various approaches with respect to two metrics: (i) weighted absolute percentage error (WAPE), and (ii) symmetric mean absolute percentage error (SMAPE). A description of these metrics can be found in Appendix B.1. We report the metrics on the test data, for each level of the hierarchy (with level 0 denoting the root) in Table 1. As a measure of the aggregate performance across all the levels, we also report the mean of the metrics in all the levels of the hierarchy denoted by Mean.

heading: 6.1 RESULTS
text: Table 1 shows the averaged test metrics for M5, Favorita, and Tourism datasets. We present only the Mean metrics for the Tourism dataset due to lack of space. Complete results with confidence intervals for all the three datasets can be be found in Appendix B.6.
We find that for all three datasets, our proposed model either yields the smallest error or close to the smallest error across most metrics and most levels. In particular, we find that our proposed method achieves the smallest errors in the mean column for all datasets in terms of WAPE and SMAPE, thus indicating good performance generally across all levels. We find that RNN+ERM in general, yields
1https://github.com/rajatsen91/deepglo 2https://github.com/liyaguang/DCRNN/
an improvement over the base RNN predictions for the higher levels closer to the root node (Levels 0 and 1), while, worsening at the lower levels. DCRNN, despite using the hierarchy as a graph also does not perform as well as our approach, especially in Tourism and M5 Datasets - possibly due to the fact that a GNN is not the most effective way to model the tree hierarchies. We notice that HierE2E performs reasonably well for the smaller Tourism while performing badly for the larger M5 and Favorita datasets - a possible explanation being that this is a VAR model that requires much more parameters to scale to large datasets. Therefore, for HierE2E we perform experiments with 50× more parameters for M5 and Favorita and report the results in Table 4 in the appendix, showing that while the results improve, it still performs much worse than our model. Overall, we find that our proposed method consistently works better or at par with the other baselines at all hierarchical levels.
Ablation study. Next, we perform an ablation study of our proposed model to understand the effects of its various components, the results of which are presented in Table 2. We compare our proposed model, to the same model without any regularization (set λE = 0 in Eq (5)), and a model consisting of only TVAR. We find that both these components in our model are important, and result in improved accuracy in most metrics.
Coherence. We also compare the coherence of our predictions to that of the RNN model and an ablated model with λE = 0. Specifically, for each node p we measure the deviation of our forecast from c(p) = 1/L(p) ∑ i∈L(p) ŷ
(i), the mean of the leaf node predictions of the corresponding subtree. Perfectly coherent predictions will have a zero deviation from this quantity. In Table 4, we report the WAPE metric between the predictions from our model ŷ and c, for each of the hierarchical levels. The leaf level predictions are trivially coherent. We find that our proposed model consistently produces more coherent predictions compared to both the models, indicating that our hierarchical regularization indeed encourages coherency in predictions, in addition to improving accuracy.
Basis Visualization. We visualize the basis generated by the BD model for the M5 validation set in Figure 3 (left). We notice that the bases capture various global temporal patterns in the dataset. In particular, most of the bases have a period of 7, indicating that they represent weekly patterns. We also show the predictions made from the various components of our model at all hierarchical levels, in Figure 3 (right). We notice that the predictions from the BD part closely resemble the general patterns of the true time series values, where as the AR model adds further adjustments to the predictions, including a constant bias, for most time series. For the leaf level (Level 3) predictions however, the final prediction is dominated by the AR model indicating that global temporal patterns may be less useful in this case.

heading: 7 CONCLUSION
text: In this paper, we proposed a method for hierarchical time series forecasting, consisting of two components, the TVAR model, and the BD model. The TVAR model is coherent by design, whereas we regularize the BD model to impose approximate coherency. Our model is fully differentiable and is trainable via SGD, while also being scalable with respect to the number of nodes in the hierarchy. Furthermore, it also does not require any additional pre-processing steps.
We empirically evaluated our method on three benchmark datasets and showed that our model consistently improved over state of the art baselines for most levels of the hierarchy. We perform an ablation study to justify the important components of our model and also show empirically that our forecasts are more coherent than the RNN baseline. Lastly, we also visualize the learned basis and observe that they capture various global temporal patterns.
In this work, we aimed at maximizing the overall performance without emphasizing performance at individual hierarchical levels. For future work, we plan to treat this is as a multi-objective problem with the aim of understanding the performance tradeoffs at various levels of the hierarchy. Finally, we also plan to extend our current model to probabilistic forecasting.
Reproducibility statement: We provide code for our model along with instructions in the supplementary file. The instructions are provided in the readme.txt file inside the zipped folder. Our code includes links to public datasets and code for preprocessing them. We also fix the random seed in our code in order for it to be reproducible. However, exact results may depend on the version of the python libraries, type of GPU and compute used.

heading: A THEORY

heading: A.1 SUPPORT RECOVERY
text: Lemma 1. SupposeB satisfies the lower eigenvalue condition (Assumption 1 in Appendix A.3) with parameterCmin and the mutual incoherence condition (Assumption 2 in Appendix A.3) with parameter γ. Also assume that the columns of the basis pool B̄ are normalized so that maxj∈Sc ‖B̄(j)‖ ≤ √ T , and the true parameter θ0 of the root satisfies
‖θ0‖∞ ≥ λL [∣∣∣∣∣∣Σ−1∣∣∣∣∣∣∞ + 4σ√LCmin ] , (9)
where |||A|||∞ = maxi ∑ j |Aij | denotes matrix operator `∞ norm, and Σ = BTB/T denotes the
empirical covariance matrix. Then for λL ≥ 2γ √ 2σ2 log d LT , with probability least 1− 4 exp(−c1Tλ 2)
(for some constant c1), the support Ŝ = {i | |α̂0| > 0} recovered from the Lasso solution (see Algorithm 1) is equal to the true support S.
Proof. We are given a pool of basis vectors B̄ from which the observed data is generated using a subset of K columns which we have denoted by B in the text. We denote the correct subset of columns by S and recover them from the observed data using basis pursuit - also known as the support recovery problem in the literature. Given the observed data and the pool of basis vectors B̄, we recover the support from the following regression problem corresponding to the root node time series. y0 = B̄α+ w0, w0 ∼ N (0, σ2I/L), (10) where α is K-sparse with the non-zero indices at S, and the non-zero values equal θ0 - the true parameters of the root node. Here we have used the fact that the root node has a 1/L times smaller variance due to aggregation. The true support S can be recovered from the observed data y0, by solving the sparse regression problem (Lasso) given in Algorithm 2. A number of standard Lasso assumptions are needed to ensure that the support is identifiable, and that the non-zero parameters are large enough to be estimated. Assuming that B̄S (= B) and α satisfy all the assumptions of Theorem 2, the theorem ensures that the true support S is recovered with high probability.

heading: A.2 PROOF OF THEOREM 1 - ERROR BOUNDS FOR REGULARIZED ESTIMATORS
text: For this proof, we assume that the support S is recovered and the true basis functionsB are known with high probability (see Section A.1). We divide the proof into multiple steps.
Step I: By Corollary 1, the OLS estimate θ̂0 (see Algorithm 2) of parameters of the root node and the OLS estimate θ̃n (see Eq. (7)) can be bounded as,
E[‖θ̂0 − θ0‖2Σ] ≤ σ2K
TL , E[‖θ̃n − θn‖2Σ] ≤
σ2K
T ∀n ∈ L(0). (11)
Step II: Next, using change of variables, we notice that the ridge regression loss for the child nodes (see Algorithm 2) is equivalent to the following.
ψ̂n = argmin ψn
1 T ‖yn −Bθ̂0 −Bψn‖22 + λ‖ψn‖22 ∀n ∈ L(0), (12)
where ψn = θn − θ̂0. The final estimate for the child parameters can be written as a sum of the ψn estimate and the root node estimate, θ̂n = ψ̂n + θ̂0. We also consider a related problem that will help us in computing the errors bounds.
ψ̃n = argmin ψn
1 T ‖yn −Bθ0 −Bψn‖22 + λ‖ψn‖22 ∀n ∈ L(0). (13)
Here we have replaced θ̂0 with the true value θ0. Note that this regression problem cannot be solved in practice since we do not have access to the true value of θ0. We will only use it to assist in the
analysis. Now we will bound the difference between the estimates ψ̂n and ψ̃n. The closed form solution for ridge regression is well known in the literature.
ψ̂n = T −1(Σ + λI)−1BT (yn −Bθ̂0)
ψ̃n = T −1(Σ + λI)−1BT (yn −Bθ0),
where Σ = BTB/T , as defined earlier. The norm of the difference of the estimates can be bounded as
ψ̂n − ψ̃n = T−1(Σ + λI)−1BTB(θ̃0 − θ̂0) =⇒ ‖ψ̂n − ψ̃n‖2Σ = (θ̃0 − θ̂0)TΣ(Σ + λI)−1Σ(Σ + λI)−1Σ(θ̃0 − θ̂0)
= (θ̃0 − θ̂0)TΣ(Σ + λI)−1Σ(Σ + λI)−1Σ(θ̃0 − θ̂0) = (θ̃0 − θ̂0)TV D [
λ3i (λi + λ)2
] V T (θ̃0 − θ̂0).
Here we have used an eigen-decomposition of the symmetric sample covariance matrix as Σ = V D[λi]V
T . We use the notation D[λi] to denote a diagonal matrix with values λi on the diagonal. The above can be further upper bounded using the fact that λi ≤ λi + λ.
‖ψ̂n − ψ̃n‖2Σ ≤ (θ̃0 − θ̂0)TV D[λi]V T (θ̃0 − θ̂0) = ‖θ̃0 − θ̂0‖2Σ. (14)
Step III: Now we will bound the error on ψ̃n and finally use it in the next step with triangle inequality to prove our result. Note that yn −Bθ0 = B(θn − θ0) +wn. Therefore, we can see from Eq. (13) that ψ̃n is an estimate for θn − θ0. Using the fact that ‖θn − θ0‖2 ≤ β and corollary 1, ψ̃n can be bounded as,
E[‖ψ̃ − (θn − θ0)‖2Σ] ≤ r2β2σ2K
Tr2β2 + σ2K . (15)
Finally using triangle inequality, we bound the error of our estimate θ̂n.
‖θ̂n − θn‖2Σ = ‖ψ̂n + θ̂0 − θn‖2Σ (Using the decomposition from Step II). ≤ 3 ( ‖ψ̂n − ψ̃n‖2Σ + ‖ψ̃n − (θn − θ0)‖2Σ + ‖θ̂0 − θ0‖2Σ ) (Using triangle and Cauchy-Schwartz inequality)
≤ 3 ( ‖ψ̃n − (θn − θ0)‖2Σ + 2‖θ̂0 − θ0‖2Σ ) (Using Eq. (14)).
Taking the expectation of the both sides, and using Eq. (11) and (15), we get the desired result.
E‖θ̂n − θn‖2Σ ≤ 3 r2β2σ2K
Tr2β2 + σ2K + 6
σ2K
TL .

heading: A.3 REVIEW OF SPARSE LINEAR REGRESSION
text: We consider the following sparse recovery problem. We are given data (X, y) ∈ Rn×d × Rn following the observation model y = Xθ∗ + w, where w ∼ N (0, σ2I), and θ∗ is supported in the indices indexed by a set S (S-sparse). We estimate θ∗ using the following Lagrangian Lasso program,
θ̂ ∈ argmin θ∈Rn
{ 1
2n ‖y −Xθ‖22 + λn‖θ‖1
} (16)
We consider the fixed design setting, where the matrix X is fixed and not sampled randomly. Following (Wainwright, 2009), we make the following assumptions required for recovery of the true support S of θ∗. Assumption 1 (Lower eigenvalue). The smallest eigenvalue of the sample covariance sub-matrix indexed by S is bounded below:
Λmin ( XTSXS n ) ≥ Cmin > 0 (17)
Assumption 2 (Mutual incoherence). There exists some γ ∈ (0, 1] such that∣∣∣∣∣∣XTScXS(XTSXS)−1∣∣∣∣∣∣∞ ≤ 1− γ, (18) where |||A|||∞ = maxi ∑ j |Aij | denotes matrix operator `∞ norm. Theorem 2 (Support Recovery, Wainwright (2009)). Suppose the design matrix satisfies assumptions 1 and 2. Also assume that the design matrix has its n-dimensional columns normalized so that maxj∈Sc ‖Xj‖2 ≤ √ n. Then for λn satisfying,
λn ≥ 2
γ
√ 2σ2 log d
n , (19)
the Lasso solution θ̂ satisfies the following properties with a probability of at least 1−4 exp(−c1nλ2n):
1. The Lasso has a unique optimal solution θ̂ with its support contained within the true support S(θ̂) ⊆ S(θ∗) and satisfies the `∞ bound
‖θ̂S − θ∗S‖∞ ≤ λn [∣∣∣∣∣ ∣∣∣∣∣ ∣∣∣∣∣ ( XTSXS n )−1∣∣∣∣∣ ∣∣∣∣∣ ∣∣∣∣∣ ∞ + 4σ√ Cmin ] ︸ ︷︷ ︸
g(λn)
, (20)
2. If in addition, the minimum value of the regression vector θ∗ is lower bounded by g(λn), then it recovers the exact support.

heading: A.4 REVIEW OF RIDGE REGRESSION
text: In this section we review the relevant background from (Hsu et al., 2012) on fixed design ridge regression. As usual, we assume data (X, y) ∈ Rn×d × Rn following the observation model y = Xθ∗ + w, where w ∼ N (0, σ2I). Define the ridge estimator θ̂ as the minimizer of the `2 regularized mean squared error,
θ̂ ∈ argmin θ∈Rn
{ 1
n ‖y −Xθ‖22 + λ‖θ‖22
} (21)
We denote the sample covariance matrix by Σ = XTX/n. Then for any parameter θ, the expected `2 prediction error is given by, ‖θ − θ∗‖2Σ = ‖X(θ − θ∗)‖22/n. We also assume the standard ridge regression setting of bounded ‖θ∗‖2 ≤ B. We have the following proposition from Hsu et al. (2012) on expected error bounds for ridge regression. Proposition 1 (Hsu et al. (2012)). For any regularization parameter λ > 0, the expected prediction loss can be upper bounded as
E[‖θ̂ − θ∗‖2Σ] ≤ ∑ j λj (λj/λ+ 1)2 θ∗j 2 + σ2 n ∑ j ( λj λj + λ )2 , (22)
where λi denote the eigenvalues of the empirical covariance matrix Σ.
Using the fact that λj ≤ tr (Σ), and x/(x+ c) is increasing in x for x ≥ 0, the above bound can be simplified as,
E[‖θ̂ − θ∗‖2Σ] ≤ tr (Σ) (tr (Σ) /λ+ 1)2 ‖θ∗‖2 + σ
2d
n
( tr (Σ)
tr (Σ) + λ )2 ≤ tr (Σ)B 2λ2 + tr (Σ) 2 σ2d/n
(tr (Σ) + λ)2
Assuming that the covariate vectors Xi are norm bounded as ‖Xi‖2 ≤ r, and using the fact that tr (Σ) ≤ r2, gives us the following corollary. Corollary 1. When choosing λ = σ
2d nB2 , the prediction loss can be upper bounded as,
E[‖θ̂ − θ∗‖2Σ] ≤ r2B2σ2d
nr2B2 + σ2d . (23)
The usual ordinary least squares bound of σ 2d n can be derived when considering the limit B →∞, corresponding to λ = 0.

heading: A.5 FURTHER DISCUSSION ON THEOREM 1
text: Upper bound ‖θn−θ0‖2 ≤ β: The upper bound β essentially bounds the distance between sibling leaf embeddings belonging to the same parent. This is directly related to an upper bound on the distance between the parent embedding θ0 and the leaf embeddings θn, as θn is essentially the mean of the leaf nodes (mean property). In many practical scenarios, the children time series of a parent may not have too different seasonal trends (for example power consumption of houses in the same neighborhood, or sales of items under the same category) resulting in the parent time series following similar trends as well.
Bounding ‖θn − θ0‖2: In most theoretical analyses of linear regression (Wainwright, 2019), the main quantity of interest is the prediction error ‖θn − θ0‖Σ rather than the parameter estimation error ‖θn−θ0‖2, as the former is directly related to the performance metric of the model. However, a bound on the parameter estimation error can be easily established using the property that ‖θn − θ0‖2 ≤ ‖θn − θ0‖Σ/ √ Cmin, where Cmin is the lower bound on the smallest eigenvalue of the sample covariance matrix as defined in Assumption 1 in Section A.3.

heading: B FURTHER EXPERIMENTAL DETAILS

heading: B.1 ACCURACY METRICS
text: In this section we define the evaluation metrics used in this paper. Denote the true values by y and the predicted values by ŷ, both n-dimensional vectors.
1. Symmetric mean absolute percent error SMAPE = 2n ∑ i |ŷi−yi| |yi|+|ŷi| .
2. Weighted absolute percentage error WAPE = ∑
i |ŷi−yi|∑ i |yi| .

heading: B.2 DATASET DETAILS
text: We use three publicly available benchmark datasets for our experiments.
1. The M5 dataset3 consists of time series data of product sales from 10 Walmart stores in three US states. The data consists of two different hierarchies: the product hierarchy and store location hierarchy. For simplicity, in our experiments we use only the product hierarchy consisting of 3k nodes and 1.8k time steps. The validation scores are computed using the predictions from time steps 1843 to 1877, and test scores on steps 1878 to 1913.
2. The Favorita dataset4 is a similar dataset, consisting of time series data from Corporación Favorita, a South-American grocery store chain. As above, we use the product hierarchy, consisting of 4.5k nodes and 1.7k time steps. The validation scores are computed using the predictions from time steps 1618 to 1652, and test scores on steps 1653 to 1687.
3. The Australian Tourism dataset5 consists of monthly domestic tourist count data in Australia across 7 states which are sub-divided into regions, sub-regions, and visit-type. The data consists of around 500 nodes and 230 time steps. The validation scores are computed using the predictions from time steps 122 to 156, and test scores on steps 157 to 192.
For the three datasets, all the time-series (corresponding to both leaf and higher-level nodes) of the hierarchy that we used are present in the training data.

heading: B.3 TRAINING DETAILS
text: Table 3 presents all the hyperparameters used in our proposed model. All models were trained via SGD using the Adam optimizer (Kingma & Ba, 2014), and the training data was standardized to mean zero and unit variance. The datasets were split into train, validation, and test sets, the sizes of which
3https://www.kaggle.com/c/m5-forecasting-accuracy/ 4https://www.kaggle.com/c/favorita-grocery-sales-forecasting/ 5https://robjhyndman.com/publications/mint/
are given in the Table 3. We used learning rate decay and early stopping using the Mean WAPE score on the validation set, with a patience of 10 for all models. We tuned the model hyper-parameters using the same metric. The various model hyper-parameters are given in Table 3.
All our experiments were implemented in Tensorflow 2, and run on a Titan Xp GPU with 12GB of memory. The computing server we used, had 256GB of memory, and 32 CPU cores, however, our code did not seem to use more than 10GB of memory and 4 CPU cores.
Mini-Batching: During each training iteration, we sample a minibatch of time series for computing the gradients. For constructing a minibatch, first a time window (of length H + F ) is uniformly randomly selected from the training time steps, after which a subset of nodes is sampled from the hierarchy tree. The time series data corresponding to this subset of nodes and the sampled time window constitutes our minibatch.
Predicting for new time steps: Once trained, our model can be used to predict for new datapoints without any retraining. However, in practice it may be beneficial to retrain the model with newer data to improve performance, even though this is not a constraint for our proposed approach.
Knowledge of the hierarchy: Our proposed approach requires the hierarchy to be known during training to be able to regularize the embeddings. In scenarios where new aggregations are presented during test time, a reasonable prediction can be produced using the following strategy: the embedding of the new aggregation is set to the mean of the embeddings of nodes in that aggregation - the rest of the model remains the same. Predictions for unseen aggregations are made in the same way as with seen aggregations. This idea is beyond the scope of the current paper, thus left for future exploration.

heading: B.4 FURTHER DETAILS ABOUT GLOBAL STATE Z
text: In many practical scenarios, the evolving global state of the set of time series may not be captured by the global covariates X only. For instance, when there is an overall increase/decrease in sales across all time series, it is captured in the past values Y rather than X . As a result, it may be required to feed in past values of time series to the BD model. However, since we cannot feed in the whole set of time series (order of 1000s) without leading to scalability issues, we choose a small set of representative time series using Non-negative Matrix Factorization, thus approximating the global state. NMF assumes that the columns of the time series matrix Y ∈ RT×N lies in the convex set spanned by a small subset of columns. However, this assumption may not hold true for many datasets, and therefore, most NMF algorithms compute an approximate factorization. We use this aforementioned subset of columns as an approximation to the global state. We would also like to
emphasize that our prediction model only sees the past valuesZH as input since the future time series values are unknown.
While NMF is one of the choices for Z, it is definitely not the only choice. Another option, PCA (or equivalently SVD), may lead to latent vectors which do not have any temporal dependencies thus requiring additional temporal regularization (Sen et al., 2019; Wilson et al., 2008). One may also use more sophisticated methods such as Temporal Latent Auto-Encoders (Nguyen & Quanz, 2021). We leave this idea for future exploration.

heading: B.5 BASELINES
text: Further details about the baselines we compare with are provided as below. For a fair comparison, we use the the same Mean Absolute Error loss functions for training, and ensure that all the baseline models have similar number of parameters.
1. RNN: We use an LSTM decoder and encoder to implement a seq-2-seq model shared across all the time series, trained using mean absolute error loss.
2. DeepGLO (Sen et al., 2019): We use the implementation released by the authors on Github6. We modify the loss function, data handling and evaluation to adapt to our setting.
3. DCRNN (Li et al., 2017): DCRNN being a GNN based approach requires a correlation graph as input. We use the official implementation released by the authors7 and provide the hierarchy tree as the input graph. The implementation uses MAE loss by default.
4. Deep Factors (DF) (Wang et al., 2019): The original implementation released by the authors makes rolling probabilistic forecasts. We implement our own version in Tensorflow using an LSTM encoder for the global model (producing point predictions) as in the original implementation, while leaving out the probabilistic local model. We manually tune the hyper-parameters for each of the datasets on the validation set.
5. L2Emb (Gleason, 2020): We implement this model using an LSTM decoder and encoder with MAE as the main training loss. In addition, we also use node embeddings which are fed as input to the encoder and decoder, and regularized according to the hierarchy as described by Gleason (2020).
6. SHARQ (Han et al., 2021): We were not able to find an official release of SHARQ. We implemented it using an LSTM based seq-to-seq model, with layer-wise training as described in the paper. We used MAE as the data fit loss function and the default squared error reconciliation loss regularizer (See Han et al. (2021) for terminology).
7. HierE2E (Rangapuram et al., 2021): The official implementation produces rolling probabilistic forecasts. We adapt it to point forecasts by using the proposed projection step on outputs from a seq-2-seq model.
8. RNN+ERM (Ben Taieb & Koo, 2019): This approach involves learning a sparse projection matrix from data, resulting in coherent predictions. We tune the sparsity parameter using the model accuracy on the validation set.

heading: B.6 MORE RESULTS
text: Table 4 show the test metrics averaged over 10 independent runs on the three datasets along with the standard deviations.
6https://github.com/rajatsen91/deepglo 7https://github.com/liyaguang/DCRNN/

Review summaries:
The paper considers point forecasting of hierarchical time series, i.e. multivariate time series with hierarchical aggregation constraints. The authors propose a new approach based on decomposing the series along a global set of basis time series where (approximate) hierarchical constraints are applied on the coefficients of the basis decomposition. Forecasts are produced using a dynamic linear autoregressive model. Compared to existing state-of-the-art hierarchical models, the proposed approach improved overall performance on forecasts at different levels of the hierarchy on several public datasets. 

The paper introduces a method for hierarchical time series forecasting. The problem setting is: given historical hierarchical univariate time series  data and given historical and future features (like holidays, etc..), try to  predict future values for all time series, while keeping the coherence  constraints of the hierarchy. The forecasting is broken down into an autoregressive part with shared  parameters for all time series, and a basis decomposition with different  weights for each time series (regularized by the hierarchy structure). The method is then being compared to different baseline methods on three  datasets.