Model selection in reconciling hierarchical time series

Forecasting is essential for supporting decision making, especially in applications that involve a lot of uncertainty. For instance, accurately forecasting the future demand of stock keeping units (SKUs) can significantly improve supply chain management [20], reduce inventory costs [54], and increase service levels [45], particularly under the presence of promotions [7, 3]. In order to obtain more accurate forecasts, forecasters typically try to identify the most appropriate forecasting model for each series from a variety of alternatives. This task, although it can provide significant improvements under perfect foresight [14], is difficult to effectively perform in practice due to model, parameter, and data uncertainty [44]. Thus, many strategies have been proposed in the literature to effectively perform forecasting model selection [15], most of which are based on the in-sample and out-of-sample accuracy of the forecasting models [55], their complexity [27], and the features that time series display [43, 40].

However, in business forecasting applications, data is typically grouped based on its context and characteristics, thus structuring cross-sectional hierarchies. For example, although the demand of an SKU can be reported at a store level, it can be also aggregated (summed) at a regional or national level. Similarly, demand can be aggregated for various SKUs of the same type (e.g., dairy products) or category (e.g., foods). As a result, hierarchical time series introduce additional complexity to the whole forecasting process since, apart from selecting the most appropriate forecasting model for each series, forecasters have also to account for coherence, i.e., make sure that the forecasts produced at the lower hierarchical levels will sum up to those produced at the higher ones [9]. In fact, coherence is a prerequisite in hierarchical forecasting (HF) applications as it ensures that different decisions made across different hierarchical levels will be aligned.

Naturally, the demand recorded at lower hierarchical levels will always add up to the observed demand at higher levels. However, this is rarely the case for forecasts which are usually produced for each series separately and are therefore incoherent. To achieve coherence, various HF methods can be used for reconciling the individual, base forecasts [50]. The most basic HF method is probably the bottom-up (BU), according to which base forecasts are produced just for the series at the lowest level of the hierarchy, being then aggregated to provide forecasts for the series at the higher levels [12]. Top-down (TD) is another option which suggests forecasting just the series at the highest level of the hierarchy and then using proportions to disaggregate these forecasts and predict the series at the lower levels [21, 8]. Middle-out (MO) mixes the above-mentioned methods, producing base forecasts for a middle level of the hierarchy and then aggregating or disaggregating them to forecast the higher and lower levels, respectively. Finally, a variety of HF methods that combine (COM) the forecasts produced at all hierarchical levels have been proposed in the literature succeeding more often than not, along with coherence, better forecasting performance [28, 58, 30].

From the HF methods found in the literature, a COM method, called minimum trace (MinT; [58]), has been distinguished due to the strong theory supporting it and the results of many empirical studies highlighting its merits over other alternatives. However, there are still circumstances under which MinT might fail to provide the most accurate forecasts. For instance, since MinT is based on the estimation of the one-step-ahead error covariance matrix, the method might be proven inappropriate when the in-sample errors of the baseline forecasting models do not represent post-sample accuracy, the assumption that the multi-step forecast error covariance is proportional to the one-step forecast error covariance is unrealistic, or the required estimations are computationally too hard to make. Moreover, since MinT treats all levels equally, it cannot be optimized with respect to certain hierarchical levels of interest. Finally, given that medians are not additive, there is no reason to expect that MinT will always improve forecasting accuracy, at least when measures that are based on averages of forecast errors are used for evaluating forecasting performance.

In such cases, simpler HF methods like the BU and the TD can be proven useful. However, there might still be inadequate evidence about which of the two methods to use [28]. For example, the BU method is typically regarded more suitable for short-term forecasts and for hierarchies which bottom series are not highly correlated but less robust to randomness. On the other hand, the TD method is usually regarded more appropriate for long-term forecasts but less accurate for predicting the series at the lower aggregation levels due to information loss. It seems that no reconciliation method can fit all kinds of HF problems and that, similarly to forecasting model selection, the appropriateness of the different HF methods depends on various factors, including the particularities of the time series [41] and the structure of the hierarchy [16, 17, 21, 1], among others.

The above findings reconfirm George Box’s quote that “all models are wrong, but some are useful” [59], and highlight the potential benefits of conditional HF, i.e., the improvements in terms of forecasting accuracy that could be possibly achieved if forecasters were able to select the most appropriate HF method according to the characteristics of the series that form a hierarchy. In this paper we propose an approach for performing such a conditional selection using time series features as leading indicators [31, 51] and Machine Learning (ML) methods for conducting the classification. Essentially, we suggest that the forecasting accuracy of the different HF methods found in the literature is closely related with the characteristics of the individual series and that, based on these relationships, “horses for courses” can be effectively identified [43].

We benchmark the accuracy of the proposed approach against various HF methods, both standard and state-of-the-art, using a large dataset coming from the retail industry. Our results suggest that conditional HF leads to superior forecasts which are significantly better than those of the individual HF methods examined. Thus, we conclude that selection should not be limited to forecasting models, as done till now, but be expanded to HF methods as well.

The remainder of the paper is organised as follows. Section 2 describes the most popular HF methods found in the literature and Section 3 introduces conditional HF. Section 4 presents the dataset used for the empirical evaluation of the proposed approach and describes the experimental set-up. Section 5 presents the results of the experiment and discusses our findings. Finally, Section 6 concludes the paper.

In this section, we discuss the TD, BU, and COM as three well-established HF methods that are widely used in the literature and in practice for reconciling hierarchical base forecasts. These methods are also the ones considered in this study, both as alternatives of the conditional HF approach to be described in the next section and as benchmarks. For a more detailed discussion on the existing HF methods, their advantages, and drawbacks, please refer to the study of [9].

Before proceeding, we introduce the following notations and parameters that will facilitate the discussion of the three methods:

We can express a hierarchical time series as $\bm{Y}_{t}=\bm{S}\bm{Y}_{k,t}$, where $\bm{S}$ is a summing matrix of order $m\times m_{k}$. For example, we can express the three-level hierarchical time series shown in Figure 1 as:

Accordingly, we can express various hierarchical structures with a unified format as $\tilde{\bm{Y}}_{n}(h)=\bm{S}\bm{G}\hat{\bm{Y}}_{n}(h)$, where $\bm{G}$ is a matrix of order $m\times m_{k}$ which elements depend on the type of the reconciliation method used, in our case the BU, TD, and COM methods [26].

Time series often depict different patterns, such as seasonality, randomness, noise, and auto-correlation [31]. As a result, there is no model that can consistently forecast all types of series more accurately than other models, even relatively simple ones [43, 15]. Similarly, although some models may perform better on a time series dataset of particular characteristics, there is no guarantee that this will always be the case [51]. For example, despite that exponential smoothing [19] typically produces accurate forecasts for seasonal series, it might be outperformed by ML methods when a large number of observations is available [48]. Thus, selecting the most appropriate forecasting model for each series becomes a promising, yet challenging task for improving overall forecasting accuracy [40].

Model selection has been extensively studied in the forecasting literature. Although there is no unique way to determine the most appropriate forecasting model for each series, empirical studies have provided effective strategies for performing this task [15]. From these strategies, the approaches that build on time series features are among the most promising ones given that the latter can effectively represent the behaviour of the series in an abstract form and match it with the relative performance of various forecasting models [46, 39, 57, 43, 31, 2].

Expert systems and rule-based forecasting were two of the early approaches to be suggested for forecasting model selection [11, 36]. [11] considered domain knowledge along with 18 time series features and proposed a framework that consisted of 99 rules to select the most appropriate forecasting model from 4 alternatives. In another study, [5] considered 6 features and 64 rules to select the most accurate forecasting model from 3 alternatives. Similarly, [6] proposed an approach to automatically extract time series features and choose the best forecasting model. [43] measured the impact of 7 time series features plus the length of the forecasting horizon on the accuracy of 14 popular forecasting models, while [31] and [51] linked the performance of standard time series forecasting models with that of various indicative features using data from well-known forecasting competitions. More recently, [40] used 42 time series features to determine the weights for optimally combining 9 different forecasting models, winning the second place in the M4 forecasting competition [38].

Inspired by the work done in the area of forecasting model selection, we posit that HF methods can be similarly selected using time series features and, based on such a selection, improve forecasting accuracy for the case of hierarchical series. In this respect, we proceed by computing various time series features across all hierarchical levels and propose using these features for selecting the HF method that best suites the hierarchy, i.e., produces on average the most accurate forecasts for all the series it consists of [2, 43]. The selection is done by employing a popular ML classification method.

The proposed approach, to be called conditional HF (CHF), is summarised below and presented in the flowchart of Figure 2.

Given that CHF builds on time series features and its accuracy is directly connected with the representativeness of the features used, as well as the capacity of the algorithm employed for selecting the most appropriate HF method, it becomes evident that choosing a set of diverse, yet finite features is a prerequisite for enhancing the performance of the proposed classification approach. There are many features that can be used to describe time series patterns. For example, [18] extracted more than 7,700 features for describing the behaviour of the time series and then summarised them into 22 canonical features, losing just 7% of accuracy [35]. Similarly, [56] and [31] suggested that a relatively small number of features can be used for effectively visualising time series and performing forecasting model selection. Based on the above, we decided to consider 32 features for the CHF algorithm so that the patters of the hierarchical series are effectively captured without exaggerating. These features included entropy, lumpiness, hurst, acffeatures, pacffeatures, nperiods, seasonal-period, trend, spike, linearity, curvature,e-acf1, e-acf10, seasonal-strength, x-acf1, x-acf10, diff1-acf1, diff1-acf10, diff2-acf1, diff2-acf10, seas-acf1, x-pacf5, diff1x-pacf5, diff2x-pacf5, seas-pacf, linearity and non-linearity, arch-test, luctanal-prop-r1, and unitroot-kpss, and were computed using the tsfeatures package for R [23]. For more details about these features, please refer to the studies of [56] and [31].

Note that CHF is flexible in terms of the method that will be employed for performing the classification. That is, users can choose their classification method of choice for identifying the most accurate HF method and reconciling the base forecasts produced for the examined hierarchy. In the preliminary experiment, we considered three methods, namely eXtreme Gradient Boosting (XGB), Support Vector Machines (SVM), and Random Forests (RF). However, since XGB performed best, we decided to exclude the last two from our analysis. This is aligned with the mainstream knowledge in data science [42, 10, 13]. XGB is an implementation of gradient boosted decision trees that is based on ensembling and uses a large number of hyper-parameters for performing the classification. Since XGB results are fairly robust to the chosen values of the hyper-parameters, we set their values manually by considering the literature and users’ experience [53]. More specifically, we set the learning rate (eta) to 0.01, colsample-bytree to 1, min-child weight to 5 , max-depth to 5, sub sample size to 0.7, and the number of iterations to 1000. We trained XGB by setting the objective to multi:softprob and evaluated its performance with mlogloss.

Note that CHF is also flexible in terms of the forecasting models that will be used for producing the base forecasts and the HF methods that will be considered for their reconciliation. For the latter case, in the present study we considered three HF methods (BU, TD, and COM), as described in Section 2. The reasoning is two-fold. First, classification methods tend to perform better when the number of classes is limited and, as a result, the key differences between the classes are easier to identify [22]. Second, we believe that the selected HF methods are diverse enough, each one focusing on different levels of the hierarchy and adopting a significantly different approach for reconciling the base forecasts. Although we could have considered more HF methods of those proposed in the literature, they are mostly variants of the examined three methods (especially the COM method) and are therefore sufficiently covered.

We should also highlight that the rolling origin evaluation of the off-line phase can be adjusted to any desirable set-up that might be more suitable to the user. For example, if computational cost is not an issue, instead of updating the forecast origin by $h$ periods at a time, a step of one period could be considered to further increase the size of the set used for training the classification method and facilitate learning. The main motivation for considering an $h$-step-ahead update is that this practice suits the way the retail firms operate when forecasting their sales and making their plans, creating also a rich set on which the ML classification method can be effectively trained, without exaggerating in terms of computational cost. Finally, although we chose to train the classification method so that the average accuracy of CHF is minimized across the entire hierarchy, this objective can be easily adjusted in order for CHF to provide more accurate forecasts for a specific level of interest. This choice depends on the decision-makers and can vary based on their focus and objectives. However, we do believe that our choice to optimize forecasting accuracy across the entire hierarchy, weighting equally all hierarchical level, is realistic when dealing with demand forecasting and supply chain management given that in such settings each level supports very different, yet equally important decisions. A similar weighting scheme was adopted in the latest M competition, M5¹¹1https://mofc.unic.ac.cy/m5-competition/, whose objective was to produce the most accurate point forecasts for 42,840 time series that represent the hierarchical unit sales of ten Walmart stores.

Table 2 displays the forecasting accuracy (MASE and RMSSE) of the three HF methods considered in this study as benchmarks as we all as classes for training the CHF algorithm in the retail dataset presented in Section 4.1. The accuracy is reported both per hierarchical level and on average (arithmetic mean of the three levels), while CHF is implemented using the XGB classifier, as described in section 4.2.

The results indicate that, on average, CHF is the best HF method according to both accuracy measures used. More specifically, CHF provides about 4%, 9%, and 2% more accurate forecasts than BU, TD, and COM, respectively, indicating that the XGB method has effectively managed to classify the hierarchies based on the feautures that their series display.

The improvements are similar if not better for the middle and bottom levels of the hierarchy. However, CHF fails to outperform TD and COM for the top hierarchical level, being about 6% and 2% less accurate, despite being still 8% better than the BU method. This finding confirms our initial claim that, depending on the hierarchical level of interest, different HF methods may be more suitable. In this study, we focused on the average accuracy of the hierarchical levels and optimised CHF with such an objective. However, as described in section 3, different objectives could be considered in order to explicitly optimise the top, middle or bottom level of the hierarchy.

In order to better evaluate the performance of the proposed approach, we proceed by investigating the distribution of the error ratios reported between CHF and the three benchmark methods examined in our study across all the 55 hierarchies of the dataset. The results, presented per hierarchical level and accuracy measure, are visualised in Figure 5 where box-plots are used to display the minimum, 1^st quantile, median, 3^rd quantile, and maximum values of the ratios, as well as any possible outliers. Values lower than unity indicate that CHF provides more accurate forecasts and vice versa. As observed, in most of the cases, CHF outperforms the rest of the HF methods, having a median ratio value lower than unity. The only exception is when forecasting level 0 and using RMSSE for measuring forecasting accuracy, where the TD method tends to provide superior forecasts. Thus, we conclude that CHF does not only provide the most accurate forecasts on average across all the 55 hierarchies, but also the most accurate forecasts across the individual ones.

To validate this finding, we also examine the significance of the differences reported between the various HF methods using the multiple comparisons with the best (MCB) test, as proposed by [33]. According to MCB, the methods are first ranked based on the accuracy they display for each series of the hierarchy and then their average ranks are compared considering a critical difference, $r_{\alpha,K,N}$, as follows:

The results of the MCB test are presented in Figure 6. If the intervals of two methods do not overlap, this indicates a statistically different performance. Thus, methods that do not overlap with the grey zone of Figure 6 are significantly worse than the best, and vice versa. As seen, our results indicate that CHF provides significantly better forecasts than the rest of the HF methods, both in terms of MASE and RMSSE. Moreover, we find that CHF is followed by COM and then by TD and BU. Interestingly, the performance of TD is not significantly different than the one of BU according to RMSSE or MASE. In this regard, we conclude that CHF performs better than the state-of-the-art HF methods found in the literature, being also significantly more accurate than the standards used for HF, such as BU and TD. As such, CHF can be effectively used for selecting the most appropriate HF method from a set of alternatives and improving the overall forecasting accuracy of various hierarchies.

We also investigate the performance of the XGB classifier in terms of precision, recall, and $\text{F}_{1}$ score. The precision metric measures the number of correct predictions in the total number of predictions made for each class. Recall (also known as sensitivity) informs us about the number of times that the classifier has successfully chosen the best HF method for each class. Finally, the $\text{F}_{1}$ score is the harmonic mean of precision and recall, computed by $F_{1}=2\frac{p*r}{p+r}$, and is used to combine the information provided by the other two metrics [49].

Before presenting the results, we should first note that the TD, BU, and COM methods have been identified as the most accurate HF method in 1,860, 2,745, and 2,820 cases, respectively. This indicates that the dataset used for training the classification method was equally populated, displaying a uniform probability distribution. Having a balanced training sample is important for our experiment since it facilitates the training of the ML algorithm (enough observations from each class are available and biases can be effectively mitigated) and provides more opportunities for accuracy improvements (if no HF method is dominant, selecting between different HF methods becomes promising).

The performance of the XGB classification method is presented in Table 3. Observe that, according to the precision metric, XGB manages to select quite accurately the BU and COM methods but finds it difficult to make appropriate selections when the TD method is the most preferable one. This indicates that, although XGB identifies the conditions under which the BU or COM methods are more preferable than the TD method, the opposite is not true.

As a final step in our analysis we investigate the significance of the time series features used by the XGB classifier, i.e., the number of times that each feature was considered by the method for making a prediction. Figure 7 shows the top-5 selected features of the classifier. As seen, the non-linearity at level 2 is the most frequently used feature and, therefore, the most critical variable for selecting a HF method. This is because non-linearity is one of the strongest features in our dataset, with promotions frequently affecting the sales of the products strongly and changing the volatility of their demand not only over promotional periods but over a number of periods [2]. Thus, the non-linearity at level 1 and level 0 are also among the top-selected features. Maximum variance shift at level 0 is another feature that is frequently selected. This may attribute to the sudden changes and spikes in sales data set caused by promotion. Finally, e-acf-a.L1, which contains the first autocorrelation coefficient of the error terms of series at level 1, is also among top selected features. One explanation is because sales are affected by promotions and have a high level of entropy.

Note that when conducting this experiment we considered another training set-up where for the classifiers, apart from time series features, we also provided information about the correlation of the series, both across levels and within each level separately. The results were similar to those reported in Table 2 and therefore we decided to exclude those features from our models for reasons of simplicity. However, we do believe that features related with the structure and the relationships of the series could help further improving the accuracy of CHF and leave this for future research.

Another extension of the proposed approach that could be also examined in a future study, is that it focuses on selecting the most appropriate hierarchical forecasting method per hierarchy. However, numerous empirical studies have shown that combining forecasts from multiple forecasting methods can improve forecasting accuracy [37]. Thus, replacing classifiers with other methods that would combine various HF methods using appropriate weights becomes a promising alternative to CHF. Simple, equal-weighted combinations of standard HF methods have already been proven useful under some settings [4], while feature-based forecast model averaging has demonstrated its potential to generate robust and accurate forecasts [40].

This paper introduced conditional hierarchical forecasting, a dynamic approach for effectively selecting the most accurate method for reconciling incoherent hierarchical forecasts. Inspired by the work done in the area of forecasting model selection and the advances reported in the field of Machine Learning, the proposed approach computes various features for the time series of the examined hierarchy and relates their values to the forecasting accuracy achieved by different hierarchical forecasting methods, like bottom-up, top-down, and combination methods, using an appropriate classification method. Based on the lessons learned, and depending on the characteristics of time series in the hierarchy, the most suitable hierarchical forecasting method can be chosen and used to enhance overall forecasting performance.

We exploited various time series features at different levels of the hierarchy that represent their behavior and trained an extreme gradient boosting classification model to choose the most appropriate type of hierarchical forecasting method for a hierarchical time series with the selected features. The accuracy of the proposed approach was evaluated using a large dataset coming from the retail industry and compared to that of three popular hierarchical forecasting methods. Our results indicate that conditional hierarchical forecasting can produce significantly more accurate forecasts than the benchmarks considered, especially at lower hierarchical levels. Thus, we suggest that, when dealing with hierarchical forecasting applications, selection should be expanded from forecasting model to reconciliation methods as well.