A Trainable Reconciliation Method
for Hierarchical Time-Series

For simplicity, we only consider the task of forecasting the one-step-ahead time-series values. To train the forecasting models we consider two alternative strategies, i.e. estimating an individual model for each time-series, or a global model for all of them. For the individual strategy, we consider linear autoregressive models ($AR(p)$), taking the lagged time-series values as input. For the datasets which present a high number of time-series ($>500$), we estimate a single global predictive model to forecast all the time-series. This approach allows us to exploit the time-series similarities and estimating a more complex model [11]. We consider the Light Gradient Boosting (LightGBM) model, taking as input the scaled lagged values, time-series specific features (such as the level in the hierarchy), and temporal features.

For both modeling strategies, we choose the best hyperparameters combination by performing ten-folds blocked cross-validation, i.e. with validation sets belonging to the same time-window. This validation technique was shown to be effective for time-series tasks [12, 13]. For the individual strategy with $AR(p)$, we perform a grid search to optimize the number of predictive lagged values $p$. For the global strategy with LighGBM, we perform a grid search to optimize the number of leaves and minimum number of observations in each leaf hyperparameters, and we keep the default values for all the others. Once the best configuration is found, we keep the ten models trained during cross-validation in order to validate the reconciler parameters, as shown in Section 4.4.

For the reconciler’s encoder, we consider feed-forward networks with ReLU activation functions and we set the size of the output layers to the number of bottom level time-series. We consider two architecture alternatives. The first is a standard fully connected network. The second is a “shrunk” version of a fully connected network such that the output for a given bottom-level time-series depends only on forecasts for itself and his parents in the hierarchy at all levels. As an example, consider the hierarchy in Figure 1. Let us consider the reconciled predictions for bottom level time-series $\tilde{D}_{t}$. In the fully connected case, we have that $\tilde{D}_{t}$ is influenced by forecasts for all the time series at all levels: $\tilde{D}_{t}=f(\hat{A}_{t},...,\hat{G}_{t})$. In the shrunk case instead, $\tilde{D}_{t}$ depends only on the forecast for itself, its parent and its grand-parent: $\tilde{D}_{t}=f(\hat{A}_{t},\hat{B}_{t},\hat{D}_{t})$. The “shrunk” version does not have the same representative power of the fully connected one. However note that it still includes bottom-up, top-down, and middle-out representation spaces, and any mix of them. Furthermore, for a given number of hierarchy levels, the number of trainable parameters grows with the number of time-series quadratically in the fully connected case and linearly in the shrunk case. To reduce the overfitting risk, we consider the shrunk architecture when the number of time-series is greater than ten times their length. However, the shrunk option can be generally considered as a tunable hyperparameter. An illustration of the two encoder architectures is given in Figure 3.

Applying the model directly to the raw time-series presents a challenge: hierarchical time-series naturally exhibit large amplitude differences, which could compromise the model training since neural networks are sensitive to the scale of the data [14]. To deal with the scale problem, we perform a normalization step. We divide each time-series by a scaling factor before feeding the data to the encoder, and we perform the inverse scaling on the encoded output. The scaling factor is set to one plus the time-series average computed in the training period.

Given the short size of the training samples, learning a multi-dimensional transformation could be challenging for the network. To simplify the training process, we want to initialize the network to the BU method, meaning that the untrained encoder should map the raw predictions to the raw bottom-level predictions. We explore two different ways to achieve this, i.e. using residual connections [15] or using a custom initialization of the network weights. We choose the second option due to its slightly higher performance, and we refer to this method as “bottom-up initialization”. As an example, consider the two architectures shown in Figure 3. With both these architectures we initialize the weights connecting each bottom level input to the corresponding output to one (i.e. the connections $\hat{D}_{t}\rightarrow{\tilde{D}_{t}},\hat{E}_{t}\rightarrow{\tilde{E}_{t}},\hat{F}_{t}\rightarrow{\tilde{F}_{t}},\hat{G}_{t}\rightarrow{\tilde{G}_{t}}$) and all the other weights and biases to zero. This technique can be trivially extended to multiple layers, assuming that the size of the hidden layers is at least equal to the output size.

To improve the stability of the model, we ensemble ten networks and we take the average prediction.

To demonstrate the flexibility of our approach, we consider two evaluation metrics. The first one is the Mean Absolute Scaled Error (MASE) [16], which scales the predictions errors by the average in-sample naïve error and has the nice properties of being symmetric, stable for values close to zero, and scale-free. These, in general, are desirable statistical properties. However, the scale-free property sometimes may not be aligned with the business point of view. If, on the one hand, scale-dependent metrics such as Mean Absolute Error (MAE) may be misleading since the prediction error for small time-series would become irrelevant in presence of much larger ones, on the other hand, higher forecast errors typically correspond to higher costs, regardless of the time-series scale. Furthermore, the MASE is undefined or unstable for constant or almost constant time-series, respectively. For these reasons, we also consider the Mean Logarithm of Absolute Error (MLAE). The MLAE depends on the time-series scale, but at the same time its dependency is “reduced” via the natural logarithmic transformation $x\mapsto log(x+1)$, and is stable for constant and almost constant time-series.

We trained the models with different loss functions based on the two target metrics: the MASE loss for the MASE metric, and the MLAE loss for the MLAE metric. Note that using the performance metric as a loss function is a convenient way to optimize for our target and is proven effective in our experiments. However, there are cases in which this may not be the optimal solution. For instance, using the MASE as a loss may be unstable in case of almost constant time series; if this is the case we suggest to use a different loss or a regularized version of it, e.g. by scaling the errors by one plus the average in-sample naïve error.

To choose the hyperparameters for the reconciliation model we exploit, again, ten folds blocked cross-validation. Specifically, we divide the training time interval into ten time-windows, and at each iteration, we perform the following steps:

1. 1.

   Set one time window as validation period, and the remaining ones as training period.

2. 2.

   Select the forecasting models trained on the training period and use them to generate forecasts for both periods. These are the reconciler inputs in the next step.

3. 3.

   Train the reconciler with forecasts in the training period and evaluate its performance with forecasts in the validation period.

Note that, for each fold, we produce the reconciler input with the forecasting models trained on the reconciler training period. One may wonder why we do not simply use the predictions of the forecasting models trained on the whole training period. This is because in doing so, the validation labels for the reconciler would be used twice, first to train the forecasting models and then to evaluate the reconciler, and this would lead to an over-optimistic evaluation.

We select the best combination of hyperparameters via random search; details about the training procedure and the hyperparameters search are shown in the Appendix.

We train all the models via backpropagation, using the AdamW [21] optimizer and a batch size equal to 128 . For the fully connected networks, we set the size of the hidden layers to the number of bottom level time-series. For shrunk networks, we set the size of the hidden layers for each sub-network relative to each bottom level time-series to 8, and therefore the number of nodes in each hidden layer is equal to 8 times the number of bottom level time-series.

We optimize the reconciler hyperparameters by evaluating 100 random combinations via the procedure described in Section 4.4. Combinations were uniformly sampled from the following values:

• •

  dropout : $[0,0.1,0.2]$

• •

  learning rate : $[1^{-3},1^{-4},1^{-5}]$

• •

  weight decay : $[1^{-1},3^{-2},1^{-2}]$

• •

  number of training epochs : $[50,100,200,500]$

• •

  number of hidden layers : $[0,1,2,3]$

In Table 3 we show the performances on the test set for each dataset, at each level of the hierarchy. Level 0 denotes the top level.