Blending gradient boosted trees and neural networks for point and
  probabilistic forecasting of hierarchical time series
In this paper we tackle the problem of point and probabilistic forecasting by
describing a blending methodology of machine learning models that belong to
gradient boosted trees and neural networks families. These principles were
successfully applied in the recent M5 Competition on both Accuracy and
Uncertainty tracks. The keypoints of our methodology are: a) transform the task
to regression on sales for a single day b) information rich feature engineering
c) create a diverse set of state-of-the-art machine learning models and d)
carefully construct validation sets for model tuning. We argue that the
diversity of the machine learning models along with the careful selection of
validation examples, where the most important ingredients for the effectiveness
of our approach. Although forecasting data had an inherent hierarchy structure
(12 levels), none of our proposed solutions exploited that hierarchical scheme.
Using the proposed methodology, our team was ranked within the gold medal range
in both Accuracy and the Uncertainty track. Inference code along with already
trained models are available at
https://github.com/IoannisNasios/M5_Uncertainty_3rd_place
Ioannis Nasios, Konstantinos Vogklis
Blending gradient boosted trees and neural networks for point and probabilistic forecasting of hierarchical time series

###### Abstract

In this paper we tackle the problem of point and probabilistic forecasting by describing a blending methodology of machine learning models that belong to gradient boosted trees and neural networks families. These principles were successfully applied in the recent M5 Competition on both Accuracy and Uncertainty tracks. The keypoints of our methodology are: a) transform the task to regression on sales for a single day b) information rich feature engineering c) create a diverse set of state-of-the-art machine learning models and d) carefully construct validation sets for model tuning. We argue that the diversity of the machine learning models along with the careful selection of validation examples, where the most important ingredients for the effectiveness of our approach. Although forecasting data had an inherent hierarchy structure (12 levels), none of our proposed solutions exploited that hierarchical scheme. Using the proposed methodology, our team was ranked within the gold medal range in both Accuracy and the Uncertainty track. Inference code along with already trained models are available at [https://github.com/IoannisNasios/M5_Uncertainty_3rd_place](https://github.com/IoannisNasios/M5_Uncertainty_3rd_place)

keywords: M5 Competition, Point forecast, Probabilistic forecast, Regression models, Gradient Boosted Trees, Neural Networks, Machine learning +
Footnote †: journal: International Journal of Forecasting

## 1 Introduction

Machine Learning (ML) methods have been well established in the academic literature as alternatives to statistical ones for time series forecasting [1; 2; 3; 4; 5]. The results of the recent M-competitions also reveal such a trend from the classical Exponential Smoothing and ARIMA methods [6] towards more data-driven generic ML models such as Neural Networks [7] and Gradient Boosted Trees [8; 9].

One of the most important results of the recent M5 Competition [1] was the superiority of Machine Learning methods, particularly the Gradient Boosted Trees and Neural Networks, against statistical time series methods (e.g., exponential smoothing, ARIMA etc.). This realization is in the opposite direction
of the findings of previous M competitions that classical time series methods were more accurate, and establishes a milestone in the domination of Machine Learning in yet another scientific domain.

The M5 forecasting competition was designed to empirically evaluate the accuracy of new and existing forecasting algorithms in a real-world scenario of hierarchical unit sales series. The dataset provided contains \(42,840\) hierarchical sales data from Walmart. It covers stores in three US states (California Texas and Wisconsin) and includes sales data per item level, department and product categories and store details ranging from February 2011 to April 2016. The products have a (maximum) selling history of 1941 days. The competition was divided into two tracks, one requiring point forecasts (Accuracy track), and one requiring the estimation of the uncertainty distribution (Uncertainty track). The time horizon for both tracks was 28 days ahead. For the Accuracy track, the task was to predict the sales for each one of the \(42,840\) hierarchical time series following day 1941. For the Uncertainty track, the task was to provide probabilistic forecasts for the corresponding median and four prediction intervals (50%, 67%, 95%, and 99%).

Table 1 presents all hierarchical groupings of the data. Level 12, containing \(30,490\) unique combinations of product per store, is the most disaggregated level. Following the competition rules, only these \(30,490\) series needed to be submitted; the forecasts of all higher levels would be automatically calculated by aggregating (summing) the ones of this lowest level. So, regardless the way we used the hierarchy information to produce 28-day ahead predictions, all we needed to actually submit was level 12 series only.

\begin{table}
\begin{tabular}{l l l r} \hline Level & Level Description & Aggr. Level & \#of series \\ \hline
1 & Unit sales of all products, aggregated for all stores/states & Total & 1 \\
2 & Unit sales of all products, aggregated for each State & State & 3 \\
3 & Unit sales of all products, aggregated for each store & Store & 10 \\
4 & Unit sales of all products, aggregated for each category & Category & 3 \\
5 & Unit sales of all products, aggregated for each department & Department & 7 \\
6 & Unit sales of all products, aggregated for each State and category & State/Category & 9 \\
7 & Unit sales of all products, aggregated for each State and department & State/Department & 21 \\
8 & Unit sales of all products, aggregated for each store and category & Store/Category & 30 \\
9 & Unit sales of all products, aggregated for each store and department & Store/Department & 70 \\
10 & Unit sales of product x, aggregated for all stores/s-tates & Product & 3,049 \\
11 & Unit sales of product x, aggregated for each State & Product/State & 9,147 \\
12 & Unit sales of product x, aggregated for each store & Product/Store & 30,490 \\ Total & & & 42,840 \\ \end{tabular}
\end{table}
Table 1: M5 series aggregation levels
## 2 Machine Learning based forecasting

Machine Learning methods are designed to learn patterns from data and they make no assumptions about their nature. Time series forecasting can be easily formulated as a supervised learning task. The goal is to approximate a function \(f(\cdot;\theta):\mathbf{R}^{d}\rightarrow\mathbf{R}\) controlled by a set of parameters \(\theta\), which corresponds to the relation between a vector of input variables (features) and a target quantity. The machine learning setup is completed once we define: a) a (training) dataset \(\mathcal{D}=(x^{(k)},y^{(k)}),\quad k=1...n\) consisting of a set of tuples containing features \(x^{(k)}\in\mathbf{R}^{d}\) that describe a target \(y^{(k)}\in\mathbf{R}\) (number of daily sales in the context of M5 competition) and b) a suitable loss function to be minimized \(\mathcal{L}(Y,f(x;\theta))\). The parameters of \(f(\cdot;\theta)\) are then tuned to minimize the loss function \(\mathcal{L}\) using the tuples of the training dataset \(\mathcal{D}\).

### Feature engineering

One of the key ingredients in any Machine Learning method is the creation of representative and information-rich input features \(x^{(\cdot)}\in\mathbf{R}^{d}\). The task of feature engineering is a largely ad hoc procedure, considered equal parts science and art, and it crucially depends on the experience of the practitioner on similar tasks. In the context of the M5 Competition, we worked with the following feature groups:

1. Categorical id-based features: This is a special type of features that take only discrete values. Their numerical encoding can take many forms and depends on the methodology chosen. 1. Categorical variables encoding via mean target value: It is a process of encoding the target in a predictor variable perfectly suited for categorical variables. Each category is replaced with the corresponding probability of the target value in the presence of this category [10]. 2. Trainable embedding encoding: A common paradigm that gained a lot of popularity since its application to Natural Language Processing [11; 12] is to project the distinct states of a categorical feature to a real-valued, low-dimensional latent space.This is usually implemented as a Neural Network layer (called Embedding Layer) and it is used to compactly encode all the discrete states of a categorical feature. Contrary to one hot encoding which are binary, sparse, and very high-dimensional, trainable embeddings are low-dimensional floating-point vectors (see Fig. 1):
2. Price related: Sell prices, provided on a week level for each combination of store and product. Prices are constant at weekly basis, although they may change through time. Using this information we calculate statistical features such as maximum, minimum, mean, standard deviation and also number of unique historical prices for each combination of store and product (level 12).
3. Calendar related
* Special events and holidays (e.g. Super Bowl, Valentine's Day, and Orthodox Easter), organized into four classes, namely Sporting, Cultural, National, and Religious.
* Supplement Nutrition Assistance Program (SNAP) activities that serve as promotions. This is a binary variable (0 or 1) indicating whether the stores of CA, TX or WI allow SNAP purchases on the examined date.
4. Lag related features: * Lag only: These features are based the historical sales for each store/product combination (level 12) for \(28+k\) days before a given date \(t\) with \(k\) ranging from 1 to 14, thus spanning two weeks. * Rolling only: Rolling mean and standard deviation of historical sales for each store/product (level 12) ending 28 days before a given date \(t\). * Lag and rolling: Rolling mean and standard deviation until a lag date in the past.

In total we devised around 80 input features, each one used to predict the unit sales of a specific product/store (Level 12) for one specific date.

### Cross Validation

The inherent ordering of time series forecasting (i.e. the time component) forces ML practitioners to define special cross validation schemes and avoid k-fold validation [14] random splits back and forth in time. We chose the following different training/validation splits.

* Validation split 1:

Figure 1: Trainable embeddings vs one hot encoding. Image taken from [13]
* Training days \([d_{1},d_{2},\ldots,d_{1940}]\)
* Validation days \([d_{1914},d_{1915},\ldots d_{1941}]\) (last 28 days)
* Validation split 2:
* Training days \([d_{1},d_{2},\ldots,d_{1885}]\)
* Validation days \([d_{1886},d_{1887},\ldots,d_{1913}]\) (28 days before last 28 days)
* Validation split 3:
* Training days \([d_{1},d_{2},\ldots,d_{1577}]\)
* Validation days \([d_{1578},d_{1579},\ldots,d_{1605}]\) (exactly one year before)

Modeling and blending was based on improving the mean of the competition metric on these splits. For each split we followed a two step modeling procedure:

**Tuning**: Use all three validation sets and the competition metric to select the best architecture and parameters of the models.

**Full train**: Use the fine-tuned parameters of the previous step to perform a full training run until day \(d_{1941}\).

The biggest added value of using multiple validation sets was the elimination of any need for external adjustments on the final prediction(see Finding 5 in [1]).

## 3 Point Forecast Methodology

Point forecast is based on ML models that predict the number of daily sales for a specific date and a specific product/store combination. In order to forecast the complete 28 days horizon we apply a recursive multi-step scheme [15]. This involves using the prediction of the model in a specific time step (day) as an input in order to predict the subsequent time step. This process is repeated until the desired number of steps have been forecasted. We found this methodology to be superior to a direct multi-step scheme were the models would predict at once 28 days in the future.

The performance measure selected for this forecasting task was a variant of Mean Absolute Scaled Error (MASE) [16] called Root Mean Squared Scaled Error (RMSSE). The measure for the 28 day horizon is defined in Eq. 1

\[RMSSE=\sqrt{\frac{\frac{1}{h}\sum\limits_{t=n+1}^{n+h}\left(y_{t}-\hat{y}_{t} \right)^{2}}{\frac{1}{n-1}\sum\limits_{t=2}^{n}\left(y_{t}-y_{t-1}\right)^{2}}} \tag{1}\]

where \(y\) is the actual future value of the examined time series (for a specific aggregation level \(l\)) at point \(t\) and \(\hat{y}\) the predicted value, \(n\) the number of historical observations, and \(h\) the 28 day forecasting horizon. The choice of
this metric is justified from the intermittency of forecasting data that involve sporadic unit sales with lots of zeros.

The overall accuracy of each forecasting method at each aggregation level is computed by averaging the RMSSE scores across all the series of the dataset using appropriate price related weights. The measure, called weighted RMSSE (WRMSSE) by the organizers, is defined in Eq. 2.

\[WRMSSE=\sum_{i=0}^{42,840}w_{i}\cdot RMSSE_{i} \tag{2}\]

where \(w_{i}\) is a weight assigned on the \(i_{th}\) series. This weight is computed based on the last 28 observations of the training sample of the dataset, i.e., the cumulative actual dollar sales that each series displayed in that particular period (sum of units sold multiplied by their respective price) [1]. The weights \(w_{i}\) are computed once and kept constant throughout the analysis.

### Models

Here we describe briefly the specific instances of GBM and neural network models used.

#### 3.1.1 LightGBM models

LightGBM [17] is an open source Gradient Boosting Decision Tree (GBDT) [9; 18] implementation by Microsoft. It uses a histogram-based algorithm to speed up the training process and reduce memory. LightGBM models have proven to be very efficient in terms of speed and quality in many practical regression problems. For the Accuracy track, we trained two variations of LightGBM models:

1. lgb_cos: Single lightGBM regression model for all available data (see Table 2)
2. lgb_nas: A different lightGBM regression model for each store (10 models in total) (see Table 3)

The target output of each LightGBM models was the sales count of a specific product. Since the target quantity (sales count) is intermittent and has a lot of zeros (\(\approx 68\%\) of daily sales are zeros), using the MSE loss function may lead to suboptimal solutions. For this reason, we implemented a special loss function that follows the Tweedie-Gaussian distribution [19]. Tweedie regression [20] is designed to deal with right-skewed data where most of the target values are concentrated around zero. In Figure 2 we show the histogram of sales for all available training data; it is obviously right skewed with a lot of concentration around zero..

The formula of the Tweedie loss function given a predefined parameter \(p\) is shown in Eq. 3. In our implementation we used the default value of \(p=1.5\) which is a good balance between the two terms:

\[\mathcal{L}(y,\hat{y})=-\sum_{k}y^{(k)}\cdot\frac{\left(\hat{y}^{(k)}\right)^ {1-p}}{1-p}+\frac{\left(\hat{y}^{(k)}\right)^{2-p}}{2-p},\ \ \hat{y}^{(k)}=f(x^{(k)};\theta) \tag{3}\]
\begin{table}
\begin{tabular}{|l|l|} \hline
**parameter** & **value** \\ \hline boosting\_type & gbdt \\ \hline  &
where \(f(\cdot;\theta):\mathbf{R}^{d}\rightarrow\mathbf{R}\) is a regression model controlled by a set of parameters \(\theta\) that maps a set of multidimensional features \(x^{(k)}\in\mathbf{R}^{d}\) to a target value \(y^{(k)}\) (daily sales count). \(\hat{y}^{(k)}\) is the output of the regression model for input \(x^{(k)}\).

#### 3.1.2 Neural Network models

We implemented the following two classes of Neural Network models:

* Keras MLP models. Keras [13] is a model-level library, providing high-level building blocks for developing neural network models. It has been adopted by Google to become the standard interface for Tensorflow [21], its flagship Machine Learning library. Using the highly intuitive description of Keras building modules one can define a complex Neural Network architectures and experiment on training and inference. In total we trained 15 slightly different Keras models and we averaged their predictions. These models were grouped in 3 slightly different architecture groups, and within each group we kept the final 5 weights during training. This strategy helped reduce variance among predictions and stabilized the result, both in every validation split and in final training. The loss function chosen was mean squared error between the actual and predicted sale count. All model groups share the shame architecture depicted in Figure 3 and are presented in detail at [https://github.com/IoannisNasios/M5_Uncertainty_3rd_place](https://github.com/IoannisNasios/M5_Uncertainty_3rd_place).
* FastAI MLP models. FastAI [22] is a Pytorch [23] based deep learning library which provides high-level components for many machine learning tasks. To tackle the regression task at hand, we incorporated the _tabular_ module that implements state of the art deep learning models for tabular data. One important thing about FastAI tabular module is again the use of embedding layers for categorical data. Similarly to the Keras case, using the embedding layer introduced good interaction for the categorical variables and leveraged deep learning's inherent mechanism of automatic feature extraction. Our modelling approach includes implementing a special Tweedie loss function to be used during training. Roughly speaking, the FastAI tabular model is the Pytorch equivalent of Keras / Tensorflow keras_nas model with the difference of using a specially implentned objective function. Empirical results on different regression contexts, mostly from Kaggle competitions, support this decision of using similar modelling methodology on complete different Deep Learning frameworks.

### Ensembling

The final part of our modelling approach was to carefully blend the predictions of the diverge set of models. We divide our models in 4 groups:
* lgb_cos: 1 LightGBM model, using all available data
* lgb_nas: 10 LightGBM models (1 model per store), using all available data for every store.
* keras_nas: 3 Keras models using only last \(17\times 28\) days of data and simple averaging.
* fastai_cos: 1 FastAI model, using all available data

All these 4 model groups were finetuned to perform best on the average of the three validation splits described in Section 2.2 above. After fine-tuning, and using the best parameters, all 4 model groups were retrained using the information available until \(d_{1941}\) and then used to produce forecasts for days \(d_{1942}\) until \(d_{1969}\).

The individual predictions were blended using geometric averaging shown in

Figure 3: Basic Keras model architecture
Eq. 4.

\[\text{blend}=\begin{cases}\left(\text{lgb\_nas}^{3.5}\cdot\text{lgb\_cos}^{1} \cdot\text{keras\_nas}^{1}\cdot\text{fastai\_cos}^{0.5}\right)^{\frac{1}{6}}& \text{for days 1-27}\\ \\ \left(\text{lgb\_nas}^{3}\cdot\text{lgb\_cos}^{0.5}\cdot\text{fastai\_cos}^{1.5 }\right)^{\frac{1}{6}}&\text{for day 28}\end{cases} \tag{4}\]

By reviewing several predictions by the keras_nas group of models, we noticed an unexpected large peak on the last day of the private set (day \(d_{1969}\)). After confirming this behavior to be dominant in almost all cases, we decided to exclude that group's prediction for the last day. This is the reasoning for not using keras_nas predictions on the lower branch of Eq 4. In Table 3.2 we present the weights for each model component as well as the validation WRMSSE score for each validation split.

### Post processing

We tried several post processing smoothing techniques to improve forecasting accuracy. In the final solution we employed a simple exponential smoothing (\(\alpha=0.96\)) per product and store id (Level 12 - \(30,490\)) that lead to substantial improvement in both validation sets and the final evaluation (private leaderboard).

We notice that, although this post-processing should be performed for both competition tracks, we were able to only use it for the Uncertainty track as it was a last minute finding (2 hours before competitions closing); as it turned out, had we applied it also for the Accuracy track submissions we would have ended 3 positions higher in the Accuracy track leaderboard.

## 4 Probabilistic Forecast Methodology

The performance measure selected for this competition track was the Scaled Pinball Loss (SPL) function. The measure is calculated for the 28 days horizon

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline  & \multicolumn{2}{|c|}{\(\text{lgb\_nas}\)} & \multicolumn{2}{|c|}{\(\text{lgb\_cos}\)} & \multicolumn{2}{|c|}{\(\text{keras\_nas}\)} & \multicolumn{2}{|c|}{\(\text{fastai\_cos}\)} & \multicolumn{2}{|c|}{\(\text{Ensemble}\)} \\ \hline Val. split 1 & 0.474 & 0.470 & 0.715 & 0.687 & 0.531 \\ \hline Val. split 2 & 0.641 & 0.671 & 0.577 & 0.631 & 0.519 \\ \hline Val. split 3 & 0.652 & 0.661 & 0.746 & 0.681 & 0.598 \\ \hline Mean / std & 0.589 & 0.08 & 0.602 & 0.09 & 0.679 & 0.05 & 0.667 & 0.02 & 0.549 & 0.03 \\ \hline \end{tabular}
\end{table}
Table 6: Competition scores for each component

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline  & \multicolumn{2}{|c|}{\(\text{lgb\_nas}\)} & \multicolumn{2}{|c|}{\(\text{lgb\_cos}\)} & \multicolumn{2}{|c|}{\(\text{keras\_nas}\)} & \multicolumn{2}{|c|}{\(\text{fastai\_cos}\)} & \multicolumn{2}{|c|}{\(\text{Ensemble}\)} \\ \hline weights \(d_{1}-d_{27}\) & 3.5 & 1.0 & 1.0 & 0.5 & \\ \hline weight \(d_{28}\) & 3.0 & 0.5 & 0.0 & 1.5 & \\ \hline \end{tabular}
\end{table}
Table 5: Final weights selected for blending
for each series and quantile, as shown in Eq. 5:

\[SPL=\frac{1}{h}\frac{\sum_{t=n+1}^{n+h}\left(Y_{t}-Q_{t}(u)\right)u\mathbf{1}\{Q_{ t}(u)\leq Y_{t}\}+\left(Q_{t}(u)-Y_{t}\right)(1-u)\mathbf{1}\{Q_{t}(u)>Y_{t}\}}{ \frac{1}{n-1}\sum_{t=2}^{2}|Y_{t}-Y_{t-1}|} \tag{5}\]

where \(Y_{t}\) is the actual future value of the examined time series at point \(t\), \(Q_{t}(u)\) the generated forecast for quantile \(u\), \(h\) the forecasting horizon, \(n\) is the number of historical observations, and \(\mathbf{1}\) the indicator function (being 1 if Y is within the postulated interval and 0 otherwise). The values \(u\) were set to \(u_{1}=0.005\), \(u_{2}=0.025\), \(u_{3}=0.165\), \(u_{4}=0.25\), \(u_{5}=0.5\), \(u_{6}=0.75\), \(u_{7}=0.835\), \(u_{8}=0.975\), and \(u_{9}=0.995\), so that they correspond to the requested median, 50%, 67%, 95%, and 99% prediction intervals.

After estimating the SPL for all time series and all the requested quantiles of the competition, the Uncertainty track competition entries were be ranked using the Weighted SPL (WSPL) shown in Eq. 6:

\[WSPL=\sum_{i=1}^{42,840}w_{i}*\frac{1}{9}\sum_{j=1}^{9}SPL(u_{j}) \tag{6}\]

where weights \(w_{i}\) were same described in Section 3.

### Quantile estimation via optimization

Using our best point forecast as median, we proceed on optimizing WSPL objective function on validation split 1 (last 28 days) and calculate the factors in Table 7. Due to time restrictions of the competition we could not extend our analysis to cover all three validation splits. These factors were used to multiply median solution (quantile \(u=0.50\)) and produce the remaining upper and lower quantiles. We assumed symmetric distributions on levels 1-9 and skewed distributions on levels 10-12. Furthermore, due to right-skewness of our sales data (zero-bounded on the left) for every level, on last quantile (99.5%) distributions proposed factor was multiplied by either 1.02 or 1.03. These factors were determined so as to minimize WSPL on validation split 1.

\begin{table}
\begin{tabular}{l l l l l l l l l l l} \hline
**Level Aggr** & \# & **0.005** & **0.025** & **0.165** & **0.25** & **0.5** & **0.75** & **0.835** & **0.975** & **0.995** \\ \hline
1 & Total & 1 & 0.890 & 0.922 & 0.963 & 0.973 & 1.000 & 1.027 & 1.037 & 1.078 & 1.143 \\
2 & State & 3 & 0.869 & 0.907 & 0.956 & 0.969 & 1.000 & 1.031 & 1.043 & 1.093 & 1.166 \\
3 & Store & 10 & 0.848 & 0.893 & 0.950 & 0.964 & 1.000 & 1.036 & 1.049 & 1.107 & 1.186 \\
4 & Category & 3 & 0.869 & 0.907 & 0.951 & 0.969 & 1.000 & 1.031 & 1.043 & 1.093 & 1.166 \\
5 & Dept & 7 & 0.827 & 0.878 & 0.943 & 0.960 & 1.000 & 1.040 & 1.057 & 1.123 & 1.209 \\
6 & State/Cat. & 9 & 0.827 & 0.878 & 0.943 & 0.960 & 1.000 & 1.040 & 1.057 & 1.123 & 1.209 \\
7 & State/Dept. & 21 & 0.787 & 0.850 & 0.930 & 0.951 & 1.000 & 1.048 & 1.070 & 1.150 & 1.251 \\
8 & Store/Cat. & 30 & 0.767 & 0.835 & 0.924 & 0.947 & 1.000 & 1.053 & 1.076 & 1.166 & 1.272 \\
9 & Store/Dept. & 70 & 0.707 & 0.793 & 0.905 & 0.934 & 1.000 & 1.066 & 1.095 & 1.208 & 1.335 \\
10 & Product & 3.049 & 0.249 & 0.416 & 0.707 & 0.705 & 1.000 & 1.218 & 1.323 & 1.720 & 2.041 \\
11 & Product/State & 9.147 & 0.111 & 0.254 & 0.590 & 0.708 & 1.000 & 1.336 & 1.504 & 2.158 & 2.662 \\
12 & Product/Store & 30.490 & 0.005 & 0.055 & 0.295 & 0.446 & 1.000 & 1.884 & 2.328 & 3.548 & 4.066 \\ \hline \end{tabular}
\end{table}
Table 7: Quantile factors for all levels
### Quantile correction via statistical information

Level 12 (the only non-aggregated one) was the most difficult for accurate estimations, and the previously calculated factors could be further improved. Statistical information on past days played a major role for this level. For every product/store id, eight statistical sales quantiles (four intervals excluding median) were calculated over the last \(13\times 28\) days (1 year) and over the last 28 days. These quantiles were first averaged and then used to correct the quantile estimation of Section 4.1 for the corresponding level. For the same reason we calculated _w_eekly sales quantiles over the last \(13\times 28\) days and \(3\times 28\) days. The final formula for estimating level 12 quantiles, which led to minimum WSPL score for validation split 1, is given below:

level 12 \[= 0.2\cdot\text{quantile estimation}\] \[+ 0.7\cdot\left(\frac{\text{statistical correction}_{13\times 28\text{ days}}+1.75\cdot\text{statistical correction}_{28\text{ days}}}{2.75}\right)\] \[+ 0.1\cdot\left(\frac{\text{weekly statistical correction}_{13\times 28 \text{ days}}+\text{weekly statistical correction}_{3\times 28\text{ days}}}{2}\right)\]

Level 11 was also corrected in a similar manner and the final formula is given below:

level 11 \[= 0.91\cdot\text{quantile estimation}\] \[+ 0.09\cdot\left(\frac{\text{statistical correction}_{13\times 28\text{ days}}+1.75\cdot\text{statistical correction}_{28\text{ days}}}{2.75}\right)\]

No corrections were applied for levels other than 12 and 11, so the respective quantile factors for levels 1-10 remain as shown in Table 7.

The evolution of our submission attempts (SPL score / ranking) for the Uncertainty track of the competition is shown in Fig. 4. This scatter plot highlights:

Figure 4: Leaderboard ranking vs. score
* the importance of ensembling diverse models, as our rank increased by almost 90 places only by ensembling,
* the contribution of statistical correction of levels 11 and 12 that led us to a winning placement.

The overall procedure of ensembling, optimizing, and correcting for the probabilistic forecasting is depicted in Fig. 5.

## 5 Discussion

We have presented a Machine Learning solution for point and probabilistic forecasting of hierarchical time series that represent daily unit sales of retail products. The methodology involves two state-of-the-art Machine Learning approaches, namely Gradient Boosting Trees and Neural Networks, tuned and combined using carefully selected training and validation sets. The proposed methodology was applied successfully on the recent M5 Competition, yielding a prize position placement in the Uncertainty track.

### Point Forecasting score breakdown

In order to get a deeper insight on the point forecasting task, we broke down the WRMSSE calculation of Eq. 2 on each hierarchical level. The resulting per-level WRMSSE for validation split 1 is shown in Fig. 6. It is obvious that

Figure 5: Flowchart of probabilistic forecasting pipeline
the performance varies significantly among different aggregation levels, with levels 10, 11 and 12 being the harder to predict. An interesting observation is that, although Level 1 is simply the aggregation of all Level 12 predictions, the Level 1 loss is less than half the magnitude of the Level 12 loss. It would seem that aggregation _cancels out_ the poor predictive capability on Level 12. This is stressed out even more by the fact that, even though mixture of traditional forecasting techniques (ARIMA, Exponential smoothing) achieved lower error on levels 10, 11, 12, than the one shown in Figure 6, the overall score was really poor.

The mean validation WRMSSE score using the final blend of Eq. 4 is 0.549, while our final submission score was 0.552 - only 0.003 units apart. This close tracking of the unseen test data error by the validation error is generally sought after by ML practitioners in real world problems, and it serves as an extra indication for the soundness of our methodology.

### Probabilistic Forecasting Factors visualized

In Fig. 7 we present a graphical plot of the factors sorted by increasing level. We notice an unsymmetrical widening of the calculated factors as the number of series of the corresponding level increases. This delta-like shape is correlated to the increasing WRMSSE error on levels 10-12 shown in Fig. 6.

### Takeaways

It is worth mentioning that our point and probabilistic forecasting results are highly correlated since we use the point forecasting as the starting point for our probabilistic analysis. It is expected that, starting from better point forecasts, will result to better probabilistic forecasts as well.

Model diversity and selected training/validation split was crucial for the overall performance and choices eliminated any need for external adjustments on the final prediction in both tracks. This is also highlighted in M5 Competition

Figure 6: Loss per Level for predictions of validation split 1
he magic multiplier 0.97 we would end up first in Accuracy Track.

Finally, although the M5 Competition dataset was hierarchical in nature, we have not used this information explicitly via a reconciliation procedure as described in [24; 6]. Hierarchy was only implicitly considered via the nature of WRMSSE objective function.

### Extensions

The competition data were carefully curated and clean providing a rich variety of potential features. There were also a multitude of ML regression methods that could be tested. However we feel that the most crucial decision was not the selection of features or ML models but the selection of a representative validation sets. Although the three validation splits described in Section 2.2 were enough to stabilize our point forecasting results we believe that using more validation splits could be beneficial and would even eliminate the need of _magic_ external multipliers (0.97) to reach the top place.

One obvious extension to the probabilistic forecast methodology was to use a weighted averaged on three validation splits instead of validation split 1 (last 28 days).
