Probabilistic Forecasting: A Level-Set Approach
Year: 2021

Abstract:
Large-scale time series panels have become ubiquitous over the last years in areas such as retail, operational metrics, IoT, and medical domain (to name only a few). This has resulted in a need for forecasting techniques that effectively leverage all available data by learning across all time series in each panel. Among the desirable properties of forecasting techniques, being able to generate probabilistic predictions ranks among the top. In this paper, we therefore present Level Set Forecaster (LSF), a simple yet effective general approach to transform a point estimator into a probabilistic one. By recognizing the connection of our algorithm to random forests (RFs) and quantile regression forests (QRFs), we are able to prove consistency guarantees of our approach under mild assumptions on the underlying point estimator. As a byproduct, we prove the first consistency results for QRFs under the CART-splitting criterion. Empirical experiments show that our approach, equipped with tree-based models as the point estimator, rivals state-of-the-art deep learning models in terms of forecasting accuracy.

heading: 1 Introduction
text: Tree-based methods are known to be robust, general-purpose and high-accuracy methods for general machine learning and data science tasks, particularly those that are not image, video or text based. This is exemplified by their popularity in competitions. In a recent interview, Anthony Goldbloom, the CEO of Kaggle, called the prominence of gradient-boosted trees on Kaggle the “most glaring difference” between what is used on Kaggle and what is “fashionable in academia,” a fact also reflected in the Kaggle Data Science and Machine Learning surveys [3, 4]. In time series forecasting in particular, tree-based methods have performed solidly in public competitions over the years [10], but have had a recent boost in attention via the M5 accuracy competition [23], arguably the most influential forecasting competition, which was dominated by tree-based methods. However, this recent surge in interest is not accompanied by methodological advances. In particular, while being able to generate probabilistic predictions ranks among the top desirable properties of forecasting techniques, the main tree-based method that outputs predictions for multiple quantiles remains Quantile Regression Forests (QRFs), which does not take advantage of gradient boosted trees algorithms.
In this manuscript we provide such an advancement and show its theoretical and empirical soundness: we introduce a novel algorithm, Level Set Forecaster (LSF), which can turn any point estimator algorithm into a probabilistic one. Applied to XGBoost [11] (in which case we refer to LSF as XLSF), and with a simple processing component that turns time series data into tabular data, this yields one of the first notable methodological advancements towards creating tree-based probabilistic forecasts
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
since the introduction of Quantile Regression Forests [27]. That said, the LSF algorithm itself is much more general: it applies to any tabular data, and it can take any point estimator algorithm.
Our contributions can be summarized as follows:
1. We propose a novel algorithm, Level Set Forecaster (LSF), which can turn any point estimator algorithm into a probabilistic one. At a high level, it groups training data whose predictions are sufficiently close, and then uses the resulting bins of true values in the training set as the predicted distributions.
2. We prove the consistency of LSF (Theorem 1); and as a byproduct, prove the consistency of QRFs (Theorem 3). By building on the methods in [32], we introduce a general framework for consistency results that goes beyond random forests, and apply it to LSF and QRF as special cases.
3. We compare LSF with the state-of-the-art models in both tabular (see Section 6.1) and forecasting tasks (See Section 6.2), and the empirical results verify the effectiveness of the proposed approach.
The rest of the article is structured as follows. In Section 2 we introduce the Level Set Forecaster algorithm, and present the consistency of the resulting estimator in Section 3. We continue by describing the method of proof of the result in Section 4, as well as presenting the novel result (Theorem 3) for QRFs. We discuss related work in Section 5, followed by the empirical studies for tabular and time series data in Section 6.

heading: 2 Level Set Forecaster
text: Given a dataset D := {(xi, yi)}i=1,...,n ⇢ Rd ⇥ R, we train a given point prediction algorithm A on D to arrive at a model f : Rd ! R. Our goal is to output a probabilistic predictor based on f , such that one could query an arbitrary quantile at any testing point x 2 Rd. In more precise terms, assuming (xi, yi)’s are sampled IID from the joint distribution (X,Y ), we wish to find an estimator for P (Y  y|X). In Algorithm 1, we propose a conceptually simple yet effective approach to turn f into a probabilistic forecaster. At a high level, Algorithm 1 creates groupings of training data points such that their predictions are “sufficiently close.” These groupings are used to partition the feature space. With a new testing data point x, we find the partition cell whose predictions are “close to” f(x), and the empirical samples that belong to the same partition cell yield the predictive distribution.
Forming the partition of the feature space. Algorithm 1 has a natural geometric interpretation. First, one trains the algorithm A on the data D to get a model f : Rd ! R. The level sets {f 1(f(xi))} are disjoint, but they do not in general satisfy that their union covers Rd, because it may be that {f(x1), ..., f(xn)} is not equal to the entire image of f . We therefore create Voronoi sets for {f(x1), ..., f(xn)}, and look at inverse images of these. Namely, writing {f(x1), ..., f(xn)} as {v1, ..., vk} without repetition and with v1 < · · · < vk, the jth Voronoi set is:
Vj = {a 2 R|vj = argmin{v1,...,vk}(|vi a|)}. (1)
We now define B̃j = f 1(Vj), with each B̃j being a disjoint union of level sets of f . The inverse Voronoi sets B̃j’s are disjoint, and their union covers all of Rd. This is not yet, however, our final partition: with the basic set-up described above, the cardinality of each bin |{yi 2 {y1, ..., yn}|xi 2 B̃j}| may be too small (often equal to 1) to ensure the desired property of consistency. In the final part of the algorithm we merge the different B̃j together, by going sequentially over B̃1, ..., B̃k (associated with the ascending values v1, ..., vk), grouping them until a sufficient bin size is ensured (repeated yi values are not ignored in this computation), and then proceeding to a new grouping, arriving at the final partition Rd = S · Bj , where each Bj is a union of B̃i’s. The intuition behind the partitioning algorithm is that for a test sample (X,Y ) and a value y 2 R we expect P (Y  y|X) to not vary wildly as X runs over a particular partition cell Bj .
Algorithm 1: Level Set Partitioning Algorithm Input: D := {(xi, yi)}i=1,...,n ⇢ Rd ⇥ R, a natural number min_bin_size, a point prediction
algorithm A (e.g., XGBoost). Train A on D, and call the resulting model f : Rd ! R. Create a dictionary pred_to_bin whose keys are {f(xi)}i, such that the value of f(xi) is the
list of yj satisfying f(xj) = f(xi). Sort the distinct elements of {f(x1), ..., f(xn)} in ascending order, and denote them [v1, ..., vk]. Initialize an empty dictionary res_pred_to_bin, and an empty list current_bin. # If the size of a bin is smaller than the minimum required size, successively merge it with the
bins on its right (higher value) until the minimum size constraint is satisfied.
for v in [v1, ..., vk] do Concatenate pred_to_bin[v] at the end of current_bin. res_pred_to_bin[v] current_bin (pass current_bin by reference) if len(current_bin) min_bin_size then
Point current_bin to a new empty list object.
# If the last bin is too small, merge to the left. if len(res_pred_to_bin[vk]) < min_bin_size then
concatenate res_pred_to_bin[vi] to res_pred_to_bin[vk], where i is the maximal index satisfying that res_pred_to_bin[vi] 6= res_pred_to_bin[vk].
pred_to_bin res_pred_to_bin
Output: A partition such that w1, w2 2 Rd are said to be in the same partition if pred_to_bin maps argmin{v1,...,vk}(|vi f(w1)|) and argmin{v1,...,vk}(|vi f(w2)|) to the same bin.
Using the partition to make inferences. In the notation of Algorithm 1, the algorithm finds the closest element vj in [v1, ..., vk] to f(x) (by binary search), and then computes and caches the quantile of pred_to_bin[vj ]. In more mathematical terms, Algorithm 1 outputs a partition Rd = S · Bj of the feature space into finitely many disjoint sets. Queried at a feature vector x 2 Bl ⇢ Rd, and for any value y 2 R, our probabilistic forecasting algorithm returns an estimator for P (Y  y|X) at X = x:
⌘̂LSF (x) := |{ 2 Ll|  y}|
|Ll| , (2)
where Ll := {yi 2 {y1, ..., yn}|xi 2 Bl} is the associated true target values for the partition cell Bl. We shall study the consistency property of the estimator ⌘̂LSF in the next section. Remark 1. Algorithm 1 resembles the classic Quantile Regression Forecasts (QRFs) ([27]) in the following sense. In QRFs, the random forests are trained in the regular sense, i.e., one splits the nodes according to the CART-splitting criterion. The probabilistic predictions are generated by “opening” the leaf nodes of each tree, and empirically sampling the predictions instead of predicting with the average values in the leaves. Here, the second step is exactly the same while the “leaf nodes” are replaced by partitions induced by the level sets of an arbitrary point estimator. In fact, if one used a decision tree for A and let min_bin_size be 1 then Algorithm 1 reduces to a QRF with one tree.
We refer to the estimator described via Equation 2 with the partitioning algorithm as in Algorithm 1 as the Level Set Forecaster (LSF)1 associated to A. If A is XGBoost, we refer to it as XLSF.

heading: 3 Main Theorem: LSF is Consistent
text: In this section, we aim to provide a quantitative understanding of the Level-set Forecaster (LSF), in particular, the consistency of LSF. The key insight that motivated the LSF algorithm was the observation that the literature regarding the consistency of Random Forests (RFs) can be extended to apply in much more general contexts.
1LSF is implemented in GluonTS: https://github.com/awslabs/gluon-ts/blob/master/src/ gluonts/model/rotbaum/README_LSF.ipynb
We remark that it is precisely because the partitions of trees grown by the CART-splitting algorithm are data-dependent that consistency of RFs has only been proven recently [32], while most other work ([7, 36, 6, 27, 31]) has focused on easier algorithms that are not used in practice. We can leverage these methods to prove consistency of LSF as well (Theorem 1); and as a byproduct we also prove consistency of QRFs under the CART-splitting algorithm (Theorem 3).
In order to have a consistency assurance, we require mild assumptions on the data generating process, and some intuitive assumptions on the point forecasting algorithm A. Assumption 1.
1. Assumption on the data generating process:
(a) There is a uniformly equicontinuous family of functions py(x) (as y varies) so that py(x) integrates to P (Y  y|X) for all y. Further assume that E(Y |X) is bounded.
(b) For every " > 0 there is a > 0 so that for every x in the image of X the probability that X lies in the ball of radius " centered at x is at least .
(c) There exists a function m(x) that integrates to E(Y |X) so that 8" > 0 9 > 0 such that for all x, x0 2 Rd and 0  y  1: |m(x) m(x0)| < implies |py(x) py(x0)| < ".
2. Assumptions on the point estimator A: Let Wn := {f(X1), ..., f(Xn)}, with k := |Wn|. Assume that A has a choice (for every n) of hyperparameters that would satisfy:
(a) Collisions among the training data f(Xi) = f(Xj) are rare in the precise sense that there exists a positive number C such that:
P ✓ C < k
n
◆ ! 1.
(b) For X an independent variable following the distribution of the Xi’s the value f(X) does not tend to be extremal among the f(Xi)’s. To be precise, for any sequence dn ! 0 of positive numbers we have that
P
✓ |{v 2Wn|f(X) < v}|
k > dn
◆ ! 1, P ✓ |{v 2Wn|v < f(X)}|
k > dn
◆ ! 1.
(c) The image of f is “dense in probability”. To be precise,
(ln(n))2 max i=2,...,k 1 Vol(Vi)! 0
in probability, where the Vi’s are the Voronoi sets from Equation 1.
(d) There’s a constant c so that 8j, 8x 2 f 1(Vj)8" > 0 the set f 1(Vj) \ B"(x) a.s. contains a ball of radius c · ".
(e) The base algorithm A is a mean square consistent estimator of the conditional mean:
E(|f(X) E(Y |X)|2)! 0. Theorem 1. Under Assumption 1, letting min_bin_size= (ln(n))2, LSF is mean square consistent. That is, for any value y 2 R, we have
E(|⌘̂LSF (X) P (Y  y|X)|2)! 0, (3)
where the convergence is uniform in y.
As a byproduct of proving this result, we also give the first consistency result (Theorem 3) for QRFs grown under the CART-splitting algorithm under an additive regression model assumption. Remark 2. Assumptions 1a and 1b are much weaker than the additive regression assumptions in [32]. Assumption 1c is what allows the base algorithm to be informative about the conditional pdfs. In light of the empirical success of LSF (see Section 6), it appears that Assumption 1c holds quite generally in real-world data. We expect that Assumptions 2a, 2b, 2c, and 2d, which simply ensure non-degenerate behavior, hold for any reasonable choice of A, but the precise statement is left as a future work. (If f is locally constant on hyperrectangles then Assumption 2d holds with c = 1p
d+1 .)

heading: 4 Method of Proof: Extend Results on RF Consistency to More General Settings
text: Our strategy for proving Theorem 1 is to generalize and strengthen existing results on the consistency of RFs, in particular through the work of [32]. We propose a unified framework that subsumes both the standard RFs (and QRFs) and the proposed LSF estimator, allowing us to arrive at the consistency results. We begin by introducing additional notation and nomenclature.

heading: 4.1 Problem Setup and Existing Results
text: Define a data-based partitioning algorithm B as being any algorithm taking a training dataset, i.e., a finite subset D = {(X1, Y1), ..., (Xn, Yn)} of Rd ⇥ R and outputting a (potentially random) partition2of Rd into disjoint sets S · Bj satisfying that there’s a constant c > 0 so that 8j8x 2 Bj8" > 0 the set Bj \ B"(x) a.s. contains a ball of radius c". (This holds if the Bj’s are hyperrectangles using c = 1p
d+1 .) Examples include Algorithm 1, as well as any tree-growing algorithm such as the
CART-splitting algorithm (see Algorithm 2, recalled below).
We define a mean-regression (resp. quantile-regression) B-estimator, trained on D as an estimator of E(⇠(Y )|X), where ⇠ is the identity map ⇠(Y ) := Y (resp. ⇠(Y ) := IYy) given by first applying B to get a partition Rd = S · Bj ; and then at inference taking a feature vector x 2 Bl to: P 2Ll ⇠( )
|Ll| ,
where Ll := {Yi 2 {Y1, ..., Yn}|Xi 2 Bl}. Now we are ready to introduce the central object of our study, generalized RF/QRF,
Definition 1. A Generalized RF (resp. Generalized QRF) grown using the partitioning algorithm B with M estimators is defined as:
⌘̂(x) :=
PM j=1 ⌘̂j(x)
M , (4)
where ⌘̂j is a mean-regression (resp. quantile-regression) B-estimator trained on only part of the data: for the jth estimator subsample uniformly without replacement an many datapoints from D to train on, where an is a (non-random) sequence of natural numbers satisfying 1  an  n.
We remark that LSF is simply a generalized QRF grown using Algorithm 1 with one estimator (M = 1) and full subsampling (an = n); and that RFs (resp. QRFs) are simply generalized RFs whose partitioning algorithm is induced by a tree-growing algorithm such as the CART splitting algorithm (see Algorithm 2).
To fix some notation, we let tn (a random variable that can depend on the data) be the number of partitions in the first B-estimator, and we let Rx,n be the partition in the first B-estimator that contains a feature vector x 2 Rd. The main question we consider is:
Problem 1. Under what assumptions is the generalized RF/QRF estimator ⌘̂(·) constructed in (4) L2 (mean square) consistent? i.e.,
E(|⌘̂(X) E(⇠(Y )|X)|2)! 0,
where ⇠ is defined either as ⇠(T ) = T or ⇠(T ) = ITy .
Our investigation starts with the recent work on the consistency of random forests. In Theorem 1 of [32] the authors proved mean square consistency in the mean regression case (⇠(T ) = T ) under some mild assumptions, where the partitioning algorithm is the one induced by a tree grown via the CART splitting criterion.
We briefly recall the tree-growing algorithm for the CART splitting criterion:
Algorithm 2: Mean and Quantile Regression Forests with the CART-Splitting Criterion; CART loss is recalled in Appendix A. Note that unlike Algorithm 1 the number of partition cells tn here is a hyper-parameter and therefore not random. Input: We fix a natural number mtry  d; and sequences of natural numbers {tn}n, {an}n such
that tn  an  n.; while not reaching tn many leaves do
For the subsampling of each tree, choose an points uniformly without replacement from the training set. At each step, and for each leaf, choose the split that minimizes CART loss, considering only mtry split dimensions sampled uniformly without replacement from {1, ..., d}.
Mean prediction: output the average value of the leaf nodes. (⇠(T ) = T ); Quantile prediction: output the average value of the leaf nodes. (⇠(T ) = ITy)
In [32], the authors make the following assumptions: Assumption 2 (Assumptions in [32]).
1. Subsamples and tree growth are same as in the CART-splitting algorithm.
2. Additive Regression Model Assumption: Y is normal with fixed variance 2 and with mean Pd i=1 mi(pi(X)), where pi is the projection onto the i th coordinate, where each
mj : [0, 1] ! R is a continuous function, and where the marginal distribution of X is uniform in the hypercube [0, 1]d.
3. Assume that limn!1 tn ln(an)9
an = 0, limn!1 tn =1.
The result of primary interest in [32] is the following. Theorem 2. ([32], Theorem 1) Under Assumptions 2 we have mean square consistency
E h |⌘̂1(X) E(Y |X)|2 i ! 0,
where ⇠(T ) := T is the identity transformation.
Remark 3. The authors [32] phrase their result as the consistency of E(⌘̂1(X)|X, {Xi, Yi}i), i.e. the case of an infinite number of trees, but they in fact reduce first to the stronger statement that a single tree is consistent.
In [27], the quantile regression case (⇠(T ) = ITy) was addressed, but only for the case of labelindependent weights, which precludes the case of the partitioning algorithm being induced by trees grown via the CART-splitting algorithm.

heading: 4.2 Consistency of Generalized RFs/QRFs
text: To begin with, we state a novel result, the first consistency proof for QRFs under the CART-splitting criterion. Theorem 3. Under Assumption 2, with the assumption 2.3 replaced by the weaker assumption that
limn!1 tn ln(an)
an = 0, limn!1 tn =1, we have that:
E(|⌘̂1(X) P (Y  y|X)|2)! 0, uniformly in y, where ⇠(T ) := ITy .
We remark that convergence in probability (the notion used in [27]) and L2 convergence coincide for QRFs because the estimator and the CDF are both bounded.
Both Theorems 3 and Theorem 1 are corollaries from the following generalization of Theorem 2: 2In order to be rigorous one really ought to use random closed sets. We will ignore this subtlety for convenience, and assume that the clever reader can deduce the proper adjustments.
Theorem 4. A generalized QRF with a single estimator (a “generalized quantile regression tree”) is mean square consistent under the data generating assumptions 1a, 1b in Assumption 1, and the
following conditions:
1. The expected variance of the conditional CDFs in the cell containing X goes to 0 uniformly in the quantile. To be more precise, for (X 0, Y 0) an IID copy of (X,Y ) restricted to the sub probability space where X 0 2 RX,n, assume that E(V (P (Y 0  y|X 0)|X)) goes to 0 uniformly in y.
2. tn ln(an)
an ! 0 in probability.
Namely, under these conditions, we have that:
E h |⌘̂1(X) P (Y  y|X)|2 i ! 0,
where the convergence is uniform in y, and where ⇠(T ) := ITy .
Theorem 1 follows from Theorem 4 as an immediate corollary:
Proof. (Theorem 1) By Assumption 2a, and our choice min_bin_size= (ln(n))2, it follows that in probability tn has the same order of magnitude as n(ln(n))2 . It is therefore clear that Assumption 2 in Theorem 4 is satisfied.
Note that RX,n is the union of at most (ln(n))2 many f 1(Vj). By Assumption 2b, with probability converging to 1 none of these Vj’s are V1 nor Vtn . Now Assumption 2c immediately implies that Vol(f(RX,n)) goes to 0 in probability. Letting (X 0, Y 0) be an IID copy of (X,Y ) restricted to the sub probability space where X 0 2 RX,n, this implies that E(V (f(X 0)|X)) ! 0. By Assumption 2e, we have that V (f(X) E(Y |X)) ! 0, and thus (by the law of total variance) E(V (f(X 0) E(Y 0|X 0)|X)) ! 0. Therefore E(V (E(Y 0|X 0)|X)) ! 0. By an application of Chebyshev’s Inequality and the boundedness of E(Y |X), Assumption 1c implies that E(V (P (Y 0  y|X 0)|X))! 0 uniformly in y.
Note that Assumptions 2a, 2b and 2c could have been replaced simply by E(V (f(X 0)|X)) going to 0 in the notation of the proof above. It is less obvious to see why Theorem 4 implies Theorem 3. The assertion boils down to proving that the assumptions of Theorem 3 imply Assumption 1 of Theorem 4. We refer to Appendix B for the proof of Theorems 3 and 4.

heading: 5 Related Work
text: The classical conformal prediction algorithm [33], whose underlying principle is bootstrapping the residuals on a validation set, is another approach for turning point estimator algorithms into probabilistic ones. The natural setting in this approach is to create prediction intervals, rather than to make predictions for particular quantiles. Unlike LSF, the prediction intervals from conformal predictions are always symmetric about the point prediction with a fixed length. More sophisticated conformal prediction algorithms exist that allow for variable length prediction intervals (e.g. [29], which takes as input two quantile regression estimators), but they no longer satisfy that they take a single point estimator algorithm to a probabilistic one. We refer to Appendix C for more details.
The main available probabilistic tree-based algorithm other than LSF (applied to a tree-based point estimator algorithm) that does not require a list of quantiles during training is QRFs [27]. We remark that LSF applied to a decision tree is equivalent to a QRF with a single tree, and that therefore QRFs with multiple trees are an ensemble of these simple LSF algorithms.
Another approach for creating tree-based probabilistic predictions is to learn the parameters of a family of conditional distributions (e.g., [9]). For tree-based methods, this has been proposed [25, 26]. We do not further compare against this method in our experiments because the implementations are not available and instead focus on comparisons with state-of-the-art methods for forecasting that use a similar approach (e.g., DeepAR [30]).
For probabilistic time series forecasting, a number of global [18] models have been proposed (e.g., [28, 12, 22, 14, 37]). For a comprehensive review for neural network models for forecasting, please refer to [16, 5]. In contrast to classical forecasting models (e.g., [17]) which are local (parameters are estimated per time series independently), tree-based methods are best employed as global models (learning parameters over the entire panel of time series, see e.g., [19] for a discussion and [35] for a concurrent contribution). In our empirical comparisons, we chose DeepAR [30] and CNN-QR [38] as global, deep-learning based models for their robust state-of-the-art performance.

heading: 6 Experiments
text: In Section 6.1 we compare XLSF (LSF with XGBoost) against QRFs (as a leading tree-based probabilistic algorithm) and Conformalized Predictions (as a baseline for turning point estimators into probabilistic ones). In Section 6.2 we apply LSF to time series data. All experiments are done using Amazon Sagemaker [21] with instance type ml.m4.16xlarge. All hyperparameters used are specified in Appendix F.

heading: 6.1 Experiments on Tabular Data
text: A main competitor for LSF for using tree-based methods to create probabilistic predictions is QRFs ([27]). In the experiments we use the skgarden implementation ([2]). QRFs are made up of multiple trees, but XLSF is one generalized tree made up of an XGBoost ensemble. It is therefore interesting to see how the two compare. By design, conformal predictions are phrased in terms of prediction intervals (a pair of quantiles) rather than a range of quantiles. For example for an ↵ = 0.1, conformal predictions aim to estimate P05 and P95 so as to have an accuracy of the prediction intervals (percent of times that the true value is in the prediction interval) that is at least 90%. We therefore chose only the pair of quantiles P05 and P95 in our experiments, as opposed to 3 quantiles as in Section 6.2.
The results are shown in Table 1, and we see that XLSF outperforms conformal predictions by a wide margin on quantile prediction, though conformal predictions’ prediction interval accuracy was more faithfully close to 90%. We remark that to reap the benefits of both algorithms, one can feed XLSF into the “Conformalized Quantile Prediction” algorithm; see [29], and the discussion in Appendix C. We also remark that QRF is much slower compared to XLSF, as shown in Table 1. The datasets in this section were small, but in Section 6.2 we will see that this makes classical QRF impractical for time series predictions, while XLSF is quite competitive.
cqr/tree/master/datasets. Accuracy is percent of times the true value was in the prediction interval. (It should revolve around 90%.) P05 and P95 are weighted quantile losses.

heading: 6.2 LSF for Time Series Prediction
text: In order to apply LSF to time series data, a few choices need to be made. First, the method for preprocessing the data into tabular data; second, whether to use a single model for the entire forecast
horizon, or one model for each timestep; and finally, one has to decide on what point estimator algorithm A to use and its hyperparameters. The winning solution to the M5 accuracy competition has already made all of these decisions, and so we can simply wrap it with LSF, and see how it compares with the top 5 solutions for the M5 uncertainty competition. We also propose a simple default choice (using only lag features, one model for each timestep, and using XGBoost) and apply it to multiple datasets.
M5 competition The M-competitions are a series of competitions led by Spyros Makridakis intended to evaluate and compare the accuracy of different forecasting methods. The most recent competition, the M5 competition [23], was in fact a twin competitions: an accuracy competition (point forecasting) and an uncertainty competition (probabilistic forecasting) with the same training data. The winning solution (by YeonJun-IN, available in [1]) to the M5 accuracy competition is a tree-based solution that trains 220 lightgbm ([20]) models (with different feature engineering set ups and/or trained on different parts of the data) and combines them to create a point forecast. The preprocessing involved more than mere lag features.
We remark that the M5 uncertainty competition requires forecasts not only of the individual time series, but also of hierarchical aggregates. In our reported metrics we only consider the bottom level of the hierarchy, since there is no trivial way to produce aggregate forecasts using the LSF-wrapping of the winning M5 accuracy competition’s solution.
We can see that wrapping the M5 accuracy competition’s winning solution is competitive among the top 5 solutions for the M5 uncertainty competition, showing the out-of-the-box we were able to turn quality point-forecasts into quality probabilistic forecasts.
Benchmarking datasets. We now turn to some benchmarking with certain default choices of preprocessing and inference on common time series datasets. To be precise, we preprocess the time series data into tabular data by sampling P context windows with only lag features, and then training a separate probabilistic forecaster for each time-step in the future. (See Algorithm 3 in the appendix for more details; see also Appendix F.)
We refer to this set-up as Tree-based Time Series Wrapper (TTSW). We remark that by design, TTSW with XLSF never predicts above (resp. below) the maximal (resp. minimal) value observed. We also include TTSW with QRFs instead of XLSF, and with lightgbm ([20]) with quantile loss, which we refer to as Quantile Regression. If one chooses XLSF or QRF there is no need to specify quantiles at training time, but for Quantile Regression, one specifies quantiles before training, and for each timestep in the forecast horizon it trains as many models as the number of quantiles requested.
By default P = 1000000 and the context window is equal to the forecast horizon. In Table 3, since QRF is much slower to train, we had to reduce P to 10000, which means it had 100 times less datapoints to train on. Even with this adjustment, it took considerably longer than the other methods. Our results in Table 3 show that TTSW (both XLSF and with Quantile Regression) perform competitively with state-of-the-art deep learning forecasting methods ([30, 38]), while maintaining advantages of tree-based methods such as interpretability. For detailed results, please refer to Appendix E.

Review summaries:
The authors propose a general method to derive probabilistic / conditional distribution predictions for any given trained base prediction model f(x) - i.e., CDF of y|x, by using an approach somewhat related to quantile regression forests (QRFs).  In this approach prediction function f(x_i) outputs are grouped with nearby values, so that groups have a minimum size from the training data.  Each group's true values {y_i} define a sample for that group that can be used to compute any CDF value - and, for instance provide a prediction interval via a set of quantiles.  A prediction for a new sample x is transformed into the probabilistic prediction by finding the nearest group to f(x), and using its samples to derive the distribution statistics.  The authors prove consistency under fairly reasonable assumptions of the prediction algorithm, as well as proving consistency for QRFs with CART splitting.  The authors also perform extensive experiments showing the proposed approach is effective at estimating distribution statistics / prediction intervals for a variety of datasets, compared to other approaches to deriving distributions from point predictions (conformal approach) as well as other probabilistic approaches (QRF, quantile regression approaches, DeepAR), showing it is both efficient and competitive in terms of accuracy, even out-performing other methods.  These include experiments with forecasting datasets showing how their approach applied to the point forecasting winner of M5 made it competitive for probabilistic forecasting with the top 5 probabilistic forecasting methods from the competition.  

The authors propose a method for turning any prediction model into a probabilistic forecaster. Briefly, the approach works by grouping training data with similar outputs, and then considers the binned true values as a prediction distribution. They prove that their algorithm is consistent by extending proofs about the consistency of random Forests to more general contexts. To evaluate their method, the authors analyze publicly available tabular datasets and time series datasets. They also adapt use their method to adapt the M5 forecasting competition’s winner to perform probabilistically, and compare their results to the top entrants in the M5 uncertainty competition. 

This paper proposes a level-set approach (LSF) to model the uncertainty in forecasting problems. It can process the result of any point forecaster into probabilistic forecasts. The main idea of LSF is to group the predictions of training data that are sufficiently close, which is also used to partition the feature space. One then obtains different bins of true values in the training set, and this can be used as the predicted distributions. The authors have proved the mean-square consistency of LSF and further extended to random forests in more general settings. Empirical evaluations on both tabular and time series competition datasets demonstrate the effectiveness of this method.