# Cross-Temporal Forecast Reconciliation at Digital Platforms with Machine
Learning111Correspondence to Ines Wilms, Maastricht University, School of
Business and Economics, P.O. Box 616, 6200 MD Maastricht, The Netherlands,
Email<EMAIL_ADDRESS>We are very grateful to Benjamin
Wolter, Pablo Perez Piskunow, Roger Caminal and Afonso Rodrigues for expert
advice, to Rob Hyndman, Artem Prokhorov and Roberto Renò for comments provided
on earlier versions of the paper, and to the participants of the IMS
International Conference on Statistics and Data Science (ICSDS 2023) for
helpful discussions. The last author was financially supported by the Dutch
Research Council (NWO) under grant number VI.Vidi.211.032.

Jeroen Romboutsa, Marie Ternesb, Ines Wilmsb
aEssec Business School, France
bMaastricht University, School of Business and Economics, The Netherlands

Abstract. Platform businesses operate on a digital core and their decision
making requires high-dimensional accurate forecast streams at different levels
of cross-sectional (e.g., geographical regions) and temporal aggregation
(e.g., minutes to days). It also necessitates coherent forecasts across all
levels of the hierarchy to ensure aligned decision making across different
planning units such as pricing, product, controlling and strategy. Given that
platform data streams feature complex characteristics and interdependencies,
we introduce a non-linear hierarchical forecast reconciliation method that
produces cross-temporal reconciled forecasts in a direct and automated way
through the use of popular machine learning methods. The method is
sufficiently fast to allow forecast-based high-frequency decision making that
platforms require. We empirically test our framework on a unique, large-scale
streaming dataset from a leading on-demand delivery platform in Europe.

Keywords: Hierarchical time series, Forecast reconciliation, Machine learning,
Cross-temporal aggregation, Demand forecasting, Platform econometrics

##  1 Introduction

Time series to be forecasted are oftentimes naturally part of a hierarchical
structure, where higher frequency and granular series are added together to
form lower frequency aggregated series. Separate forecasts of each series
rarely conserve this hierarchy, and forecast reconciliation methods are
therefore required. In this paper, we consider novel forecast reconciliation
methods for on-demand delivery platforms that require forecasts at different
levels of cross-sectional and temporal aggregation. We introduce non-linear
forecast reconciliation based on popular machine learning methods to capture
the complex interdependencies of platform data.

Platforms such as Uber, Lyft, GrubHub, UberEats or DoorDash are nowadays
omnipresent in the global economy. They operate in high frequency and on a
large number of verticals (regions, product categories, etc.). Indeed, their
market place is typically split up in different geographical regions, so a
natural cross-sectional aggregation scheme arises from many individual
delivery areas over zones towards few market places. Moreover, a temporal
aggregation scheme naturally arises since, on the granular end, fast
operational decisions (think in terms of minutes) are needed to ensure the
platform’s service couriers are at the right time and location to serve
consumer demand promptly and to determine compensation schemes for couriers
through dynamic pricing. On the coarser end, strategic business decisions also
require long-term planning since the budget available for each delivery area
is set typically using daily demand forecasts. Accurate and coherent, i.e.
reconciled, demand forecasts across all levels of the cross-sectional and
temporal hierarchy are therefore key to the business’ success and to support
aligned decision making across different planning units.

The literature to reconcile forecasts has been pioneered by Hyndman et al.,
(2011) and Athanasopoulos et al., (2009), see Hyndman and Athanasopoulos,
(2021) for a textbook introduction and Hyndman and Athanasopoulos, (2014) for
an introduction to practitioners given the strong interest from industry.
Typical initial applications concern tourism where demand series are subject
to a cross-sectional hierarchy across different geographical regions, though
many other cases appear more recently, see for instance Caporin et al., (2023)
for a financial application to realized volatility forecasting. While forecast
reconciliation methods initially focused on producing coherent forecasts
across cross-sectional hierarchies, the literature on forecast reconciliation
in the temporal direction (Athanasopoulos et al.,, 2017) as well as
combinations of both in cross-temporal frameworks (Kourentzes and
Athanasopoulos,, 2019; Di Fonzo and Girolimetto, 2023a, ; Di Fonzo and
Girolimetto, 2023b, ) is much sparser (Petropoulos et al.,, 2022) though in
recent years it has proliferated. We refer to Athanasopoulos et al., (2023)
for an extensive review on forecast reconciliation in cross-sectional,
temporal and cross-temporal hierarchies, to Hollyman et al., (2021) for the
link with the forecast combination literature and Athanasopoulos et al.,
(2024) for recent innovations in hierarchical forecasting.

Methodological considerations to forecast reconciliations are discussed in
Hyndman et al., (2016), whereas theoretical results are first established by
Wickramasuriya et al., (2019) by putting forward an optimal linear forecast
reconciliation method, minimizing the sum of forecast error variances across
the collection of time series; see also Van Erven and Cugliari, (2015) for a
game-theoretic approach to forecast reconciliation and Panagiotelis et al.,
(2021) for a unification of the former two works. Traditionally, forecast
reconciliation is thus typically achieved by linearly combining base forecasts
of all the time series in the hierarchy. Recently, however, applications of
machine learning algorithms surface to reconcile forecasts in a non-linear
way.

Machine learning (ML) algorithms nowadays have become an integral component of
every forecaster’s toolkit, boasting a remarkable track record as demonstrated
by various forecasting competitions, such as the M4 and M5 competitions
(Makridakis et al.,, 2020; Makridakis et al., 2022b, ), and successful
economic applications (e.g., Medeiros et al.,, 2021). ML methods have become
standard choices for obtaining forecasts as their algorithms demand minimal
user intervention– typically limited to tuning parameters –and their
sophisticated off-the-shelf implementations are readily available,
facilitating seamless operationalization in production environments. Moreover,
ML methods offer flexibility by imposing no constraints on the number of
features for a given target, and they come with built-in mechanisms for
feature selection and for capturing complex dependencies such as non-
linearities.

In the forecast reconciliation literature, a cross-sectional ML-based
reconciliation scheme is proposed by Spiliotis et al., (2021), on which our
work is mainly built. They first collect base forecasts for all series in the
cross-sectional hierarchy and then train a tree-based ML algorithm (random
forest or XGBoost) with the bottom level series as response and one-step-ahead
forecasts of all the time series in the hierarchy as input features. The ML
fitted algorithm is hereby used to reconcile the bottom level forecasts in a
non-linear way, given all base forecasts, and coherent forecasts across the
complete hierarchy are obtained via bottom-up aggregation. In contrast to
standard linear reconciliation approaches that reconcile forecasts by
minimizing the reconciliation error between the base and reconciled forecasts,
the ML-based reconciliation approach is directly geared towards obtaining the
minimum forecast combination error to improve out-of-sample forecast accuracy
while still adhering to coherence. This results in substantial forecast
improvements as demonstrated on two applications to tourism and retail
industry. Additionally, ML-based methods are more flexible than the standard
ones as they allow for non-linear combinations of the base forecasts and
selectively combine them in a direct and automated way without requiring that
the information in the entire hierarchy must be used for producing reconciled
forecasts.

Anderer and Li, (2022) use deep learning and LightGBM tree-based methods to
propose a variant of the bottom-up method proposed by Spiliotis et al.,
(2021), accounting for biases that are difficult to observe at the bottom
level. The latter forecasts are adjusted to obtain a higher forecasting
accuracy on the upper levels of the hierarchy. Their approach works well on
the M5 competition dataset consisting of over ten thousand hierarchical retail
sales time series (see Makridakis et al., 2022c,  and Makridakis et al.,
2022a,  for details). Abolghasemi et al., (2022) use tree-based methods and
lasso regression in their dynamic top-down bottom-up method applied to fast-
moving consumer goods and the M5 competition dataset. Finally, a new stream of
research uses neural networks for both forecasting and reconciliation;
examples of such end-to-end modeling approaches are Theodosiou and Kourentzes,
(2021) and Wang et al., (2022), which bring promising evidence on various
datasets.

Our paper also proposes an ML-based reconciliation approach, extending the
work of Spiliotis et al., (2021) in the following ways. From a methodological
point of view, we generalize their ML-based forecast reconciliation method for
cross-sectional hierarchies to cross-temporal hierarchies. Additionally
considering reconciliation across the temporal direction requires important
design choices with respect to the ML model used for reconciliation since the
inputs consist of mixed-frequency time series, namely all temporal frequencies
in the temporal hierarchy. Besides, since the ML algorithm is a key backbone
in our forecast reconciliation approach, we carefully investigate the effect
of (i) the ML algorithm used, thereby considering not only random forest and
XGBoost as in Spiliotis et al., (2021) but also LightGBM due to its successful
application in the M5 competition; and (ii) tuning hyperparameters to make
careful considerations on the trade-off between improvements in forecast
accuracy and computational time. Finally, we consider an extensive set of base
forecasts, ranging from academic standards to industry standards (for our
platform application) as well as a combination of several base forecasts
methods.

From an application viewpoint, we showcase the performance of our forecast
reconciliation method on platform demand data streams, thereby bringing novel
insights into the performance of forecast reconciliation methods on a modern
and unique dataset. We use UK demand data (September 2022-September 2023) from
a last mile logistics platform where the bottom level series are 30-minute
demand series for 141 delivery areas in London and in Manchester, their two
largest markets. We aggregate to hourly and daily frequency in the temporal
dimension, and towards zone and market level in the cross-sectional dimension.
Platform data have several features that are distinct from the traditional
datasets considered in the forecast reconciliation literature. First, they are
high-frequency with, our top, hence lowest frequency, level (i.e. daily)
consisting of the bottom level considered in earlier work such as recent M5
competition (other works even consider lower frequency settings such as
monthly to yearly data). Indeed, computational speed and fast decision making
at digital platforms are key to their success, thereby requiring forecasts at
the intra-day level. A second distinctive feature of platform data is their
streaming character which means that they face non-stationarities or shifts
for example due to expansion, competition among platforms, sales of certain
clients, changing customer’s preferences or non-regular events such as sports
events. The data on Manchester highlight the performance of the ML forecast
reconciliation method in presence of data shifts. We compare our approach with
state-of-the art reconciliation benchmarks from the recent literature.

We obtain several insights based on our application. First, our main empirical
finding is that ML-based forecast reconciliation for platform data can result
in substantial forecast accuracy improvements compared to existing linear
reconciliation methods. However, ML-based reconciliation is not uniformly
superior to linear methods. In fact, for the most important series, i.e. high
frequency delivery area level, random forest based reconciliation typically
performs best, for city level series LightGBM based reconciliation yields the
smallest forecast errors. On the other hand, linear reconciliation methods
outperform XGBoost based reconciliation on the bottom level series. The latter
result hold for XGBooost without tuning its hyperparameters to keep computing
time low. However, ignoring the huge computational cost, we also study the
impact of hyperparameter tuning of the ML algorithms and find for example that
XGBoost performs in this case at par with random forest in the lowest levels
of the hierarchy. A second empirical finding is that ML-based reconciliation
seems to work equally well irrespective of the quality of the base forecasts.
Hence, exploring a large model space for obtaining base forecasts is
unnecessary. Finally, data streams are prone to data shifts so we compare the
forecast reconciliation methods on the Manchester market that is impacted by a
major shift in our evaluation period. The results show that, while forecast
performance of all methods suffer during a data shift, ML-based reconciliation
is able to revert quickly back to pre-shift forecast accuracies.

The rest of the paper is organized as follows. Section 2 introduces the
platform data we analyze in this paper. Section 3 presents our ML-based
forecast reconciliation approach for cross-temporal hierarchies. Section 4
describes the forecast set-up we use to evaluate the out-of-sample forecast
performance of our ML-based reconciliation approach versus state-of-the-art
benchmarks whereas Section 5 presents the results. Finally, Section 6
concludes.

##  2 Platform Data

We use a unique dataset of demand for deliveries from a leading on-demand
last-mile logistics platform in Europe, Stuart, which connects businesses to a
fleet of geographically independent couriers.222The data are provided to us in
the context of ongoing research collaborations. Due to a confidentiality
agreement, we are not allowed to distribute or report actual demand data.
Demand data across all figures are therefore normalized between 0-100. All
analyses were carried out with the original data. We analyze UK London demand
30-minute frequency data from September 5, 2022 to September 24, 2023.333In
addition to the London market, we also analyze the case of Manchester in
Section 5.4. To efficiently organize parcel deliveries, London is split into
18 zones which are constituted by the 117 bottom-level delivery areas.
Planning at the delivery area level is of particular importance because it
determines the remuneration paid to the couriers and and it therefore directly
helps to optimize the platform’s efficiency.

Deliveries are carried out on a daily basis between 7am and 11:30pm in all
areas.444Occasional overnight demand is integrated in the last 30-min of the
business day. This creates a balanced dataset of $13,090$ observations for
each of the 136 (=1+18+117) time series at the 30-minute frequency. For
business planning purposes, the series also require temporarily aggregation to
the hourly and the daily frequency, thereby resulting in temporally aggregated
time series that consist of $13,090/2=6,545$ and $13,090/34=385$ observations
respectively. The forecast reconciliation therefore consists of three cross-
sectional and three temporal layers implying 408 base time series.

Figure 1: Time series plots for full sample of the area City of London
belonging to the zone London Central and the market London.

Figure 1 plots demand for the City of London delivery area, London Central
zone, and London market at 30-minute, hourly and daily frequencies over the
sample period. We observe some typical characteristics of platform data.
First, intra-day area data is spiky because of peak demand periods triggered
by events such as flash sales, food promotions, bad weather. These peaks are
largely tempered at daily frequency. Holidays can also cause large drops in
delivery demand, this is particularly noticeable in the end-of-year period.
Second, there are strong seasonality patterns at all levels of the hierarchy.
Figure A.1 in Appendix A displays these seasonality patterns. We clearly see
that the intra-day seasonality (visible from both the 30-minute and hourly
time series) is mainly driven by food deliveries. The bottom row of the plots
traces daily seasonality, and illustrates that there can be substantial
differences between weekends and the rest of the week. In general, London
faces higher weekend demand while the opposite is true for the delivery area
City of London. This implies that there exists heterogeneity in demand
patterns at the delivery area level. In fact, unreported plots for other
delivery areas reveal that very different dynamics appear.

Figure 2: Left: Aggregated demand in each area (circle) in the London market.
Right: Average demand (top) and proportion of zeros (bottom) across areas; the
median is displayed in solid black, 10% and 90% quantiles in dotted black, 25%
and 75% in dashed black.

Figure 2 illustrates this heterogeneity across areas. In the left panel, we
visualize the aggregated demand over the sample period for each delivery area.
Apart from larger activity above River Thames, there is no specific spatial
pattern in demand volume, low- and high-volume delivery areas co-exist. The
largest volume area is Romford, outside the city center. The City of London
and Westminster area volumes respectively represent about 15% and 5% of
Romford. In the right top panel of Figure 2, we plot the 30-minute average
demand throughout the business day for each delivery area. The solid black
median line visualizes well the lunch and dinner periods, the colored area
curves stress the large variability, except for early morning demand which is
overall low. The right bottom panel of Figure 2 zooms into the proportion of
zeros across areas. Between 10am and 9pm, the proportion of 30-minute zero
demand spells is low for most delivery areas, Stirling Road (upper green
curve) however has mostly zero demand throughout the entire business day.
Large diversity appears in the morning and evening hours with delivery areas
that are idle and others that are active.

##  3 Machine Learning Based Forecast Reconciliation

In Section 3.1, we review the cross-temporal framework using the general
notation from Athanasopoulos et al., (2023). In Section 3.2 we introduce our
ML-based forecast reconciliation approach for cross-temporal hierarchies.

###  3.1 Cross-Temporal Hierarchies

We start with some general notation so that our methodology of the next
section can be widely applied. We denote $y_{b,t}$ for $b=1,\ldots,n_{b}$ and
$t>0$ the high-frequency time series observation for bottom level $b$ at time
$t$. The cross-sectional hierarchy is simply obtained by summing over elements
in the set of the $n_{b}$ bottom level time series. This yields $n_{a}$
additional high-frequency time series resulting in a total of $n=n_{a}+n_{b}$
series.

Regarding the temporal hierarchy, note that in a platform setting data is
streaming but for simplicity let us consider $t=1,\ldots,m$ with $m$ the
number of observations to sum over when going from the highest to the lowest
sampling frequency in the hierarchy. If we define $k=k_{1},...,k_{p}$ the $p$
temporal aggregation levels, with $k_{1}=1$ and $k_{p}=m$, then the non-
overlapping temporally aggregated series, for each series $i=1,\ldots,n$ and
$k=k_{1},...,k_{p}$, are given by

$x_{i,j}^{[k]}=\sum_{t=(j-1)k+1}^{jk}y_{i,t}\qquad j=1,...,m/k.$

Next, denoting $\tau$ as the observation index of the lowest frequency (top
temporal level), we can stack these observations in ($m_{k}\times 1$) vectors
$\boldsymbol{x}_{i,\tau}^{[k]}=\left[x_{i,m_{k}(\tau-1)+1}^{[k]},x_{i,m_{k}(\tau-1)+2}^{[k]},\ldots,x_{i,m_{k}\tau}^{[k]}\right]^{\top},$
where $\top$ denotes the transpose of a vector and $m_{k}=m/k$.

Finally, the cross-temporal hierarchy is obtained by considering for each
variable $i=1,\ldots,n$ of the cross-sectional hierarchy, the entire temporal
hierarchy. Notation wise, we collect in a ($\sum_{l=1}^{p}m_{k_{l}}\times 1$)
vector all temporally aggregated time series
$\boldsymbol{x}_{i,\tau}=\left[x_{i,\tau}^{[m]},\ldots,{\boldsymbol{x}_{i,\tau}^{[1]}}^{\top}\right]^{\top},$
and then stacking all $n$ variables into the vector
$\boldsymbol{x}_{\tau}=\left[\boldsymbol{x}_{1,\tau}^{\top},\ldots,\boldsymbol{x}_{n,\tau}^{\top}\right]^{\top}$.

$x_{i,1}^{[1]}$$x_{i,2}^{[1]}$$x_{i,1}^{[2]}$$x_{i,1}^{[34]}$……30-minhourlydaily$x_{i,33}^{[1]}$$x_{i,34}^{[1]}$$x_{i,17}^{[2]}$$A_{1}$$Z_{118}$$M_{136}$……$A_{2}$$A_{3}$$Z_{135}$$A_{116}$$A_{117}$AreaZoneMarkettemporal
hierarchycross-sectional hierarchy……AreaZoneMarketcross-temporal hierarchy

Figure 3: Cross-sectional (left), temporal (middle) and cross-temporal (right)
hierarchy in the platform application.

In our application, regarding the cross-sectional hierarchy the bottom-level
consists of $n_{b}=117$ delivery areas that are aggregated into the upper
levels consisting of 18 zones (middle) and finally one market (top level), as
visualized in the left panel of Figure 3. The total number of variables is
thus $n=117+18+1=136$. Regarding the temporal hierarchy, we have three
frequencies $p=3$, (30-min - hourly - daily), as visualized in the middle
panel of Figure 3, $m=34$ and $k_{1}=1,\ k_{2}=2$ and $k_{3}=34$ so that
$\boldsymbol{x}_{i,\tau}^{[1]}$ ($\boldsymbol{x}_{i,\tau}^{[2]}$) collects the
$m_{1}=34$ ($m_{2}=17$) 30-min (hourly) demand observations on day $\tau$, and
$x_{i,\tau}^{[34]}$ corresponds to the single daily observation for variable
$i$. Finally, the combined cross-temporal hierarchy is displayed in the right
panel of Figure 3.

###  3.2 Forecast Reconciliation via Machine Learning

We are now ready to introduce our machine learning-based forecast
reconciliation approach for cross-temporal hierarchies, thereby extending the
work of Spiliotis et al., (2021) for cross-sectional hierarchies.

The general idea behind ML-based reconciliation is to model the relationship
between each bottom-level time series $\boldsymbol{x}_{b,\tau}^{[1]}$
($b=1,\ldots,n_{b}$) of the cross-temporal hierarchy and the base forecasts
across the entire cross-temporal tree
$\widehat{\boldsymbol{x}}_{i,\tau}^{[k]}$ ($i=1,\ldots,n$ and
$k=k_{1},\ldots,k_{p}$). We use a time series rolling-window procedure to
construct these base forecasts as visualized in Figure 4. The rolling-window
size of the estimation sample is denoted by $Q$ (length of the blue bars) and
the validation sample with forecast horizon by $H$ (length of the green bars),
both quantities are expressed in number of observations for the top temporal
aggregation level (i.e. the daily time series in our application). Then, for
each cross-sectional series at temporal aggregation level $k$ and rolling
window iteration $r=1,\ldots,R$ (i.e. $R=4$ in Figure 4 as used in our
application), we estimate base forecast models on the estimation sample,
consisting of $m_{k}\times Q$ observations, and obtain $m_{k}\times H$-step
ahead forecasts.

$N$${\color[rgb]{0.29,0.56,0.89}\definecolor[named]{pgfstrokecolor}{rgb}{0.29,0.56,0.89}Q}$$Q+R\times
H$$Q+H$$N+H$ Figure 4: Rolling-window procedure for machine-learning based
forecast reconciliation.

Next, from the observations in the validation sample $\mathcal{V}$ consisting
of observations $\mathcal{V}=\\{Q+1,\ldots,Q+R\times H\\}$, we construct the
response vector and features matrix for the ML-model which can be generally
written as

$\displaystyle\boldsymbol{y}_{b,\mathcal{V}}=f_{b}(\widehat{\boldsymbol{X}}_{b,\mathcal{V}})+\boldsymbol{\varepsilon}_{b,\mathcal{V}},$
(1)

where $f_{b}(\cdot)$ denotes the forecast function to be trained for bottom-
level series $b$, see Section 4.2 for specifications, and
$\boldsymbol{\varepsilon}_{b,\mathcal{V}}$ represents the error term. Note
that every bottom-level series in the cross-temporal hierarchy (i.e. 30-min
delivery area demand in our application) is predicted by a separate ML model,
thereby allowing the forecast reconciliation approach to adapt to different
patterns in each series. The response variable
$\boldsymbol{y}_{b,\mathcal{V}}$ in (1) is simply given by the $m_{1}\times
R\times H$ observations of the bottom-level series for variable $b$; hence

${\boldsymbol{y}}_{b,\mathcal{V}}=\left[{\boldsymbol{x}_{b,Q+1}^{[1]}}^{\top},\dots,{\boldsymbol{x}_{b,Q+RH}^{[1]}}^{\top}\right]^{\top}.$

As features matrix to the ML-model, we use
$\widehat{\boldsymbol{X}}_{b,\mathcal{V}}=$

-6.8em-7em  $\begin{bmatrix}\hat{\boldsymbol{x}}_{1,Q+1}^{[1]}&\dotsc&\hat{\boldsymbol{x}}_{b-1,Q+1}^{[1]}&\hat{\boldsymbol{x}}_{b+1,Q+1}^{[1]}&\dotsc&\hat{\boldsymbol{x}}_{n,Q+1}^{[1]}&\hat{\boldsymbol{x}}_{b,Q+1}^{[1]}&\hat{\boldsymbol{x}}_{b,Q+1}^{[k_{2}]}\otimes\boldsymbol{1}_{k_{2}}&\dotsc&\hat{x}_{b,Q+1}^{[k_{p}]}\otimes\boldsymbol{1}_{k_{p}}\\\ \vdots&\ddots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\ \hat{\boldsymbol{x}}_{1,Q+RH}^{[1]}&\dotsc&\hat{\boldsymbol{x}}_{b-1,Q+RH}^{[1]}&\hat{\boldsymbol{x}}_{b+1,Q+RH}^{[1]}&\dotsc&\hat{\boldsymbol{x}}_{n,Q+RH}^{[1]}&\hat{\boldsymbol{x}}_{b,Q+RH}^{[1]}&\hat{\boldsymbol{x}}_{b,Q+RH}^{[k_{2}]}\otimes\boldsymbol{1}_{k_{2}}&\dotsc&\hat{x}_{b,Q+RH}^{[k_{p}]}\otimes\boldsymbol{1}_{k_{p}}\end{bmatrix},$ which contains various base forecasts on the validation sample, where $\otimes$ denotes the Kronecker product and $\boldsymbol{1}$ is a column vector of ones. More specifically, the first set of predictors, highlighted in yellow, exploits information in the cross-sectional hierarchy and contains the $n_{b}$ highest-frequency (i.e. 30-min in our application) base forecasts for all bottom-level variables, as well as the $n_{a}$ highest-frequency base forecasts for the upper level series in the cross-temporal hierarchy (i.e. zone and market in our application). The second set of predictors highlighted in blue exploits information in the temporal hierarchy and contains all temporal (i.e. 30-min, hourly and daily in our application) base forecasts for the bottom-level variable $b$ for which we are constructing our ML-model. Notice that due to the frequency mismatch, each of the lower-frequency (i.e. hourly and daily in our application) at level $k_{a}$ for $a=2,\ldots,p$ are repeated respectively $k_{a}$ times. We hereby make the assumption that temporal gain in reconciling forecasts mainly comes from the focal variable’s own temporally aggregated series to keep this second set of predictors compact. Finally, note the central role for the bottom level variable $b$ itself (highlighted in yellow and blue) whose base forecast we aim to improve upon in this reconciliation step.

Finally, to reconcile the actual out-of-sample demand forecasts of horizon $h$
for forecast evaluation purposes, we use the complete sample up to
$N=Q+R\times H$ to obtain the $m_{k}\times H$ step ahead forecasts for each
cross-sectional series at temporal aggregation level $k$. This is visualized
in the bottom part of Figure 4, where the gray bar represents the training
sample (of length $N$) on which the base forecast models are estimated and the
red bar (of length $H$) represents the test set on which we evaluate out-of-
sample forecast performance. The out-of-sample base forecasts
$\widehat{\boldsymbol{x}}_{i,\tau}^{[k]}$ ($i=1,\ldots,n$ and
$k=k_{1},\ldots,k_{p}$) are used as inputs to the previously estimated ML
models (one for each bottom-level series) to provide the reconciled bottom-
level forecasts $\widetilde{\boldsymbol{x}}_{b,\tau}^{[1]}$
($b=1,\ldots,n_{b}$) of the cross-temporal hierarchy. As a last step, these
reconciled bottom-level forecasts are then aggregated as explained in Section
3.1 to obtain coherent forecasts across the complete hierarchy.

##  4 Design of the Forecast Study

We start by discussing in Section 4.1 the forecast set-up we use to evaluate
out-of-sample forecast performance. In Section 4.2, we summarize the base
forecast models all forecast reconciliation methods rely on, together with the
machine learning methods used by our procedure to obtain the reconciled
forecasts. Finally, Section 4.3 outlines the benchmark forecast methods
against which we compare the performance of our procedure.

###  4.1 Forecast Set-up

To compare the forecast performance of our forecast reconciliation method
against its benchmarks, we use a standard rolling-window approach using the
period February 20, 2023 to September 24, 2023 (31 weeks) as test sample. In
each rolling window iteration, we use the most recent six months of data to
estimate the base models on the estimation set (first five months), and the
most recent month of data to estimate the ML model on the validation set.
Hence, in line with the notation used in Figure 4, we take $Q=140$ days, $H=7$
days and $R=4$. Note that since the demand data in our platform application
consists of non-negative integers we ensure this by rounding the base
forecasts prior to using them as inputs in the ML model, i.e.
$\hat{x}_{i,t}^{{[k]}^{\star}}=\text{max}(0,\lfloor\hat{x}_{i,t}^{[k]}\rceil)$,
we similarly round the bottom-level reconciled forecasts before aggregating
them.

In a next step, we combine the estimation and validation sample ($N=168$ days)
to estimate the base models and obtain final out-of-sample base forecasts as
well as reconciled forecasts on the test set. As forecast horizon, we consider
one week of data (i.e. horizon $m_{k}\times H$ for temporal aggregation level
$k$) as this corresponds to the planning horizon for the compensation schemes
of couriers through dynamic pricing which crucially relies on demand forecasts
as inputs. Hence, for the 30-min series, we obtain $34\times 7=238$ forecasts;
for the hourly series we obtain $17\times 7=119$ forecasts and for the daily
series we obtain $H=7$ forecasts. This concludes one iteration in the outer
rolling window used for forecast evaluation purposes. In the next iteration,
we move one-week ahead and repeat the same steps. We roll forward until we
reach the end of our sample, thereby resulting in 31 outer rolling windows.

Figure B.1 in the Appendix B provides a visualization of iteration $r$ and
$r+1$ in the outer rolling window. We hereby highlight the estimation (blue
bars) and validation (green bars) sets for the inner rolling window used to
construct the inputs to the ML model (see Section 3.2). Due to our inner
rolling-window set-up, the training set (blue bar) and forecast horizon (green
bar) of the last three (inner) rolling windows in iteration $r$ overlap with
the first three (inner) rolling windows in iteration $r+1$, as indicated with
the shading. This offers a clear computational advantage as beyond iteration
$r=1$, we can re-use the base forecasts for the first three weeks from the
previous iteration and need to fit the forecasting model only once to obtain
the base forecasts for the fourth week to construct the validation set.

To evaluate the forecast performance of the various forecast methods, we use
the Weighted Absolute Percentage Error (WAPE), also known as normalized
version of the mean absolute error (MAE). The WAPE loss function for each
temporal factor $k$ and cross-sectional time series $i=1,\ldots,n$ is given by

$\text{WAPE}^{[k]}_{i}=\frac{\sum_{j=1}^{T_{\text{test}}/k}|A_{i,j}^{[k]}-F_{i,j}^{[k]}|}{\sum_{j=1}^{T_{\text{test}}/k}|A_{i,j}^{[k]}|},$

which measures the overall deviation of the forecasted values $F_{i,j}^{[k]}$
from the actual values $A_{i,j}^{[k]}$, for all time points $j$ in the test
set of size $T_{\text{test}}/k$ where $T_{\text{test}}$ denotes the number of
high-frequency (i.e. 30-min) observations in the test set. We then average
these WAPEs across all cross-sectional units belonging to the same cross-
temporal level; for instance the overall WAPE for the bottom-level series
(i.e. 30-min for the $n_{b}$ areas) is then given by
$\text{WAPE}^{[1]}=(1/n_{b})\sum_{b=1}^{n_{b}}\text{WAPE}^{[1]}_{b}$.

The WAPE is a popular forecast metric to evaluate performance at delivery
platforms (i) because of its scale-independence thereby facilitating
comparisons across heterogeneous areas as well as different time/geographical
aggregation levels, (ii) it is agnostic to zeros as it is very unlikely that
the denominator is zero, thereby making it a suitable metric even when the
demand in an area is low or intermittent, and (iii) it emphasizes accuracy for
items with larger demand, thereby prioritizing accurate forecasting of demand
during busy periods. Note that as robustness check we also report results for
the root mean squared scaled error as well as the symmetric mean absolute
percentage error.

###  4.2 Base Forecasts and Machine Learning Methods

Our procedure allows practitioners to choose their preferred base forecast
models as well as their preferred ML method to perform the ML-based
reconciliation. In this paper, we consider four base forecast models and three
popular ML methods to this end.

Base Forecasts. We consider one popular industry forecast model for platform
data, two all-round time series models for forecasting univariate time series
and a forecast combination of the former three as ensemble, base forecast
model. The first base forecast model is a simple yet often-used model in the
platform industry which we will label as “Naive” since the forecast for a
specific slot (30-min, hourly or daily) of next week is simply the value
observed in the same slot and day of the previous week. The following two base
forecast models are oftentimes used in the hierarchical forecasting literature
to obtain base forecasts, namely the class of SARIMA and exponential smoothing
models. Full details on the base forecast methods are available in Appendix B.

Machine Learning Methods. We consider three popular ML models, namely random
forest, XGBoost and LightGBM. The former two are also considered in Spiliotis
et al., (2021) whereas the latter has recently shown great success on the M5
competition (Makridakis et al., 2022a, ) consisting of an application with
hierarchical retail sales time series, thereby making all suitable candidates
for our platform data. More details on the ML models are available in Appendix
B.

For all three ML methods, we report main results based on their standard
implementations with default tuning parameters (see Appendix B). We do this,
first to avoid excessive computing times on our streaming platform data, but
second, and importantly, also because off-the-shelf implementations of tree-
based methods often attain excellent performance across a variety of settings
(e.g, Januschowski et al.,, 2022). In Section 5.3, we investigate the
improvements in forecast performance one can obtain by tuning these ML models.

Finally, when reporting our main results, we use the sum of squared errors as
loss function to train all ML models. However, our framework allows
practitioners to use another loss function as they see suited for the data at
hand. We also investigated the performance of the ML-based reconciliation
methods when using the Tweedie loss function instead of the standard squared
loss, since the former showed superior performance in the M5 competition
(Januschowski et al.,, 2022; Makridakis et al., 2022a, ) on forecasting
hierarchical time series as it is specifically geared towards sparse (zero-
inflated) target data, which some of the considered delivery areas also
display (see Figure 2). However, using a Tweedie loss instead of the regular
squared loss function did not improve forecast results of our proposed ML
reconciliation methods on the considered data, and are hence omitted but
available from the authors upon request.

###  4.3 Benchmark Reconciliation Methods

We compare our new ML forecast reconciliation procedure against the base
forecasts and five state-of-the-art linear reconciliation benchmarks for
cross-temporal hierarchies, as summarized in Table 1 and implemented in the
FoReco package (Girolimetto and Di Fonzo,, 2023) in `R` (R Core Team,, 2023).
Apart from the Bottom Up method, the reconciliation step requires specifying
how to estimate the covariance matrix of the in-sample fit errors since the
benchmark methods are geared towards minimizing reconciliation errors as
opposed to our procedure which is geared towards minimizing forecast
combination errors. For the tcs, cst and ite benchmark methods, the temporal
reconciliation covariance matrix estimator is the series variance scaling
matrix (thf_comb = “wlsv”) which is the diagonal matrix that contains the
estimated variances of the in-sample residuals across each level. The cross-
sectional reconciliation covariance matrix estimator is a shrinkage covariance
matrix of the in-sample residuals (hts_comb = “shr”), where the off-diagonal
elements are shrunk towards zero. For the oct method cross-temporal based
covariance matrix estimator is simply containing the in-sample residual
variance on the main diagonal (csts_comb = “wlsv”). Full details on the cross-
temporal benchmark methods are available in Di Fonzo and Girolimetto, 2023a .

Table 1: Overview of cross-temporal forecast reconciliation benchmarks.

Name |  Description
---|---
Bottom Up |  Simple aggregation of bottom level forecasts for any cross-temporal hierarchy
tcs |  First-temporal-then-cross-sectional cross-temporal forecast reconciliation - Kourentzes and Athanasopoulos, (2019)
cst |  First-cross-sectional-then-temporal cross-temporal forecast reconciliation - Di Fonzo and Girolimetto, 2023a
ite |  Iterative cross-temporal forecast reconciliation - Di Fonzo and Girolimetto, 2023a
oct |  Optimal combination cross-temporal forecast reconciliation in a single reconciliation step - Di Fonzo and Girolimetto, 2023a

Notes: This table summarizes the benchmark cross-temporal forecast
reconciliation methods.

##  5 Results

This section is divided in four parts. In Section 5.1, we discuss the overall
forecast performance of the forecast reconciliation methods, followed by
Section 5.2 where we present detailed insights on our ML-based reconciliation
approach. In Section 5.3, we discuss the trade-off between computing time and
forecast accuracy when incorporating hyperparameter tuning into our ML-based
reconciliation proposal. Finally, in Section 5.4, we discuss how shifts in the
demand data streams affect the performance of the reconciliation methods.

###  5.1 Overall Forecast Performance

Table 2 summarizes the forecast accuracy across the complete cross-temporal
hierarchy for the different combinations of base forecasts and reconciliation
methods.555Note that for the Naive base method, the base and bottom up
forecasts are identical, and no results are reported for the linear forecast
reconciliation methods since the original base forecasts are already
reconciled.  Overall, all methods find that forecasting is easier higher up in
the hierarchy, i.e. lower temporal frequency and larger geographic areas,
which confirms the general findings in the literature. For the set of
considered base forecasts and reconciliation methods, WAPE falls from around
30% to 5% for respectively 30-minute delivery areas and the daily London
hierarchy levels.

Table 2:  London forecast reconciliation results.

Base | Forecast | Delivery Area | Zone | Market
---|---|---|---|---
Forecasts | Method |  |  |  |  |  |  |  |  |
|  | 30min | hour | day | 30min | hour | day | 30min | hour | day
Naive | Base | 0.3576 | 0.2703 | 0.1017 | 0.1475 | 0.1196 | 0.0638 | 0.0793 | 0.0756 | 0.0559
| Random forest | 0.2762 | 0.2162 | 0.0976 | 0.1271 | 0.1081 | 0.0636 | 0.0819 | 0.0782 | 0.0576
| XGBoost | 0.3017 | 0.2338 | 0.0993 | 0.1352 | 0.1132 | 0.0654 | 0.0841 | 0.0796 | 0.0588
| LightGBM | 0.2898 | 0.2252 | 0.0991 | 0.1321 | 0.1110 | 0.0650 | 0.0842 | 0.0796 | 0.0587
ETS | Base | 0.3180 | 0.2654 | 0.0837 | 0.1935 | 0.1769 | 0.0539 | 0.1621 | 0.1521 | 0.0471
| Bottom Up | 0.3180 | 0.2682 | 0.1755 | 0.1842 | 0.1716 | 0.1378 | 0.1523 | 0.1501 | 0.1313
| tcs | 0.3105 | 0.2591 | 0.1587 | 0.1769 | 0.1641 | 0.1290 | 0.1440 | 0.1420 | 0.1229
| cst | 0.3070 | 0.2557 | 0.1536 | 0.1699 | 0.1570 | 0.1191 | 0.1349 | 0.1331 | 0.1117
| ite | 0.3115 | 0.2610 | 0.1621 | 0.1786 | 0.1664 | 0.1321 | 0.1459 | 0.1444 | 0.1256
| oct | 0.3099 | 0.2577 | 0.1545 | 0.1765 | 0.1633 | 0.1267 | 0.1435 | 0.1411 | 0.1206
| Random forest | 0.2809 | 0.2229 | 0.1064 | 0.1360 | 0.1192 | 0.0729 | 0.0959 | 0.0933 | 0.0673
| XGBoost | 0.3161 | 0.2545 | 0.1189 | 0.1439 | 0.1249 | 0.0735 | 0.0920 | 0.0889 | 0.0634
| LightGBM | 0.2960 | 0.2374 | 0.1164 | 0.1387 | 0.1209 | 0.0737 | 0.0926 | 0.0896 | 0.0645
SARIMA | Base | 0.3118 | 0.2584 | 0.0894 | 0.1793 | 0.1656 | 0.0568 | 0.1490 | 0.1445 | 0.0526
| Bottom Up | 0.3118 | 0.2605 | 0.1637 | 0.1808 | 0.1682 | 0.1339 | 0.1511 | 0.1489 | 0.1286
| tcs | 0.3049 | 0.2524 | 0.1486 | 0.1742 | 0.1612 | 0.1253 | 0.1436 | 0.1416 | 0.1201
| cst | 0.3013 | 0.2488 | 0.1424 | 0.1679 | 0.1549 | 0.1164 | 0.1351 | 0.1334 | 0.1100
| ite | 0.3054 | 0.2536 | 0.1511 | 0.1753 | 0.1630 | 0.1279 | 0.1450 | 0.1435 | 0.1225
| oct | 0.3043 | 0.2510 | 0.1443 | 0.1718 | 0.1583 | 0.1210 | 0.1408 | 0.1385 | 0.1164
| Random forest | 0.2801 | 0.2206 | 0.1002 | 0.1305 | 0.1132 | 0.0657 | 0.0883 | 0.0861 | 0.0586
| XGBoost | 0.3195 | 0.2528 | 0.1109 | 0.1385 | 0.1182 | 0.0659 | 0.0845 | 0.0815 | 0.0551
| LightGBM | 0.2933 | 0.2327 | 0.1067 | 0.1320 | 0.1134 | 0.0641 | 0.0832 | 0.0806 | 0.0538
Forecast | Base | 0.2950 | 0.2365 | 0.0850 | 0.1502 | 0.1333 | 0.0538 | 0.1134 | 0.1089 | 0.0485
Combination | Bottom Up | 0.2950 | 0.2358 | 0.1289 | 0.1493 | 0.1333 | 0.0991 | 0.1120 | 0.1097 | 0.0938
| tcs | 0.2925 | 0.2349 | 0.1215 | 0.1470 | 0.1310 | 0.0962 | 0.1093 | 0.1073 | 0.0913
| cst | 0.2909 | 0.2331 | 0.1182 | 0.1438 | 0.1277 | 0.0908 | 0.1046 | 0.1027 | 0.0853
| ite | 0.2932 | 0.2360 | 0.1237 | 0.1484 | 0.1328 | 0.0985 | 0.1112 | 0.1095 | 0.0934
| oct | 0.2917 | 0.2336 | 0.1182 | 0.1454 | 0.1290 | 0.0928 | 0.1070 | 0.1046 | 0.0879
| Random forest | 0.2758 | 0.2161 | 0.1018 | 0.1285 | 0.1108 | 0.0717 | 0.0871 | 0.0844 | 0.0666
| XGBoost | 0.3026 | 0.2355 | 0.1011 | 0.1336 | 0.1129 | 0.0671 | 0.0836 | 0.0803 | 0.0601
| LightGBM | 0.2900 | 0.2260 | 0.0999 | 0.1301 | 0.1103 | 0.0664 | 0.0831 | 0.0799 | 0.0593

Notes: This table shows forecast accuracy measured in weighted absolute
percentage error (WAPE). The forecast horizon is one week-ahead. The best
forecast reconciliation results for each base forecast (Naive, ETS, SARIMA,
Forecast Combination) are highlighted in gray. The linear benchmark methods
are defined in Table 1.

The ML-based forecast reconciliation results using random forest are
consistently outperforming all benchmarks at the delivery area and zone levels
at all frequencies. For example, using 30-minute frequency SARIMA base
forecasts, we obtain a WAPE of 0.3013 for the cst method (i.e. the best linear
forecast reconciliation method) compared to 0.2801 for random forest based
reconciliation (i.e. the best non-linear forecast reconciliation method).
However, ML-based forecast reconciliation does not systematically dominate
linear methods as shown for instance by XGBoost with a WAPE of 0.3195.
Furthermore, ML-based forecast reconciliation seems to be rather insensitive
to the choice of the base forecasts. In fact, the Naive, ETS and SARIMA
results are qualitatively the same. For example in the case of random forest,
the respective WAPEs are very close with 0.2762, 0.2809 and 0.2801
respectively.

To assess the significance of our findings, we report in Figure 5 the results
of the accuracy ranking based multiple comparisons with the best (MCB) test, a
methodology popularized in the forecasting literature since Koning et al.,
(2005). We compute the MCB test across all levels in the cross-temporal
hierarchy and report results for each of the four base forecast methods. To
this end, we use the R-package tsutils (Kourentzes,, 2023) with function
nemenyi() to compute the MCB test. The test ranks the performance of the
examined methods across the series being forecast comparing their average
ranks by considering a critical difference, determined through a confidence
interval. For each base forecast panel in Figure 5, the methods that do not
overlap with the best ranked procedure (as highlighted by the gray zone) are
significantly worse than the best. Visual inspection of the results confirms
our overall finding that random forest reconciliation dominates the other
procedures.

Figure 5:  MCB test with confidence level 0.95, using all levels in the
hierarchy. Methods that do not overlap with the gray zone indicating the
confidence interval of the best method are significantly worse than the best.

The important business question is whether ML-based forecast reconciliation
helps improving the base forecast methods for high-frequency (i.e. 30-minute
and hourly) delivery area levels, because these forecasts are direct inputs
for automated decision making. The answer is yes and the forecast gains are
between 10 and 15%. As an illustration for random forest reconciliation with
ETS base forecasts, the WAPE is 0.2809 compared to 0.3180 for non-reconciled
base forecasts. It also turns out that these high-frequency delivery area
level gains become even stronger at the zone and market levels. In fact, the
forecast gains mount to even 40% at the market level. For example, the London
30-minute WAPEs decrease from 0.1621 to 0.0959 with ETS base forecasts and
random forest reconciliation. Here, the linear reconciliation methods only
deliver marginal accuracy gains, with the best (cst) WAPE being 0.1349.
However, we also note that at the London level, the overall smallest WAPEs are
obtained using Naive base forecasts for 30-minute and hourly frequencies.

Finally, in the bottom panel of Table 2, we provide forecast combination
results, that is by combining with equal weights the base forecasts. We find
this strategy to slightly pay off compared to the best performing base
methods. For example, at the 30-min delivery area level, the random forest ML
method for the base forecasts are respectively 0.2762, 0.2809 and 0.2801 while
the forecast combination yields 0.2758. The gains are overall stronger for the
linear reconciliation methods compared to the non-linear ones, yet the former
are, overall, outperformed by the ML-based reconciliation methods. Still in
practice, we generally do not know a priori which base forecasts method will
perform best out-of-sample, and moreover, this likely varies across the
different levels of the cross-temporal hierarchy. Then, a forecast combination
base forecast approach forms an effective practice to obtain accurate and
robust base forecasts, which can still be further improved upon via ML-based
reconciliation to obtain coherent forecasts.

To sum up, random forest forecast reconciliation performs particularly well
but not all ML methods are doing better than linear methods. The strongest
gains are obtained for the 30-minute and hourly frequencies. Forecast
combination, i.e. averaging the base forecasts, is (marginally) beneficial.
Note that these main findings are qualitatively the same when we measure
forecast accuracy in terms of the symmetric mean absolute percentage error
(SMAPE), a popular asymmetric loss function, or in terms of the root mean
squared scaled error (RMSSE), a popular metric in the forecast reconciliation
literature (see Hyndman and Koehler,, 2006 for a discussion on these metrics)
instead of WAPE. We refer to Tables C.1 and C.2 in Appendix C for respectively
the SMAPE and RMSSE results.

###  5.2 ML-Based Reconciliation Insights

From the previous results, we find that random forest based forecast
reconciliation outperforms the considered benchmarks, especially at the more
granular levels in the hierarchy. To gain more insights, we investigate the
link between forecast accuracy and size of the delivery areas. In addition, we
investigate what drives forecast performance by computing SHAP values which
are considered to be the current state-of-the-art method for interpreting ML
models, see Molnar, (2022) for an introduction.

Figure 6 plots forecast performance (WAPE, vertical axis) of each delivery
area (dot) based on its average demand size (horizontal axis). The left panel
shows this relationship using un-reconciled base forecasts, and we find a
smirk shape pattern with many WAPEs above 0.5. The Naive base forecasts are,
overall, worse than ETS and SARIMA when average demand is below 20. The middle
figure shows the random forest forecast reconciliation results and we see that
the smirk flattens out, with hardly any noticeable difference between the base
forecast methods, and few areas having WAPEs above 0.5. The right figure
highlights that random forest based forecast reconciliation has a capacity to
“correct” bad quality base forecasts as the relative WAPEs of the ML-
reconciled forecasts over the base forecasts are mostly below 1.

Figure 6: Forecast performance measured in WAPE for base forecasts (left), for
ML-reconciled forecasts via random forest (middle) and relative improvement in
WAPE of the latter over the former versus average 30-min demand for the
different delivery areas (dots).

Next, for each random forest model with Naive base forecasts in the (outer)
rolling-window setup, we compute SHAP values (Shapley et al.,, 1953; Lundberg
and Lee,, 2017) for all observations in the test set. To this end, we use the
package treeshap (Komisarczyk et al.,, 2023) in R, a fast implementation for
tree ensemble models (Lundberg et al.,, 2018). Figure 7 shows the resulting
relative variable importance plots by variable group for City of London as an
illustration. Given the large number of features used in the ML-model, we
focus on the 36 highest ranked features in terms of variable importance (which
together account for more than 80% of the overall variable importance).
Overall, the proportions are stable over the test sample. While City of
London’s own 30-min and hourly streams count for a sizable portion in the
20-30% range, the other areas 30-minute streams count together for almost 40%.
Zone and market level, or own daily information contributes only marginally.

Figure 7:  Variable importance in random forest ML model with Naive base
forecasts for the delivery area City of London.

###  5.3 Computational Aspects and Hyperparameter Tuning

Our main results above are obtained by training the ML algorithms with their
default hyperparameters as provided by the software packages. By doing so, we
save a substantial amount of computing time which is an important
consideration when deploying forecast reconciliation in a production
environment. However, in this section, we investigate whether the performance
of the ML-based forecast reconciliation methods can be further improved by
hyperparameter tuning. To this end, we use the common practice of grid search
to tune the hyperparameters of the random forest algorithm, since there are
only few hyperparameters, whereas we use Bayesian optimization to tune the
hyperparameters of the XGBoost and LightGBM algorithms as the latter two have
a large amount of hyperparameters. Details on the tuning process are provided
in Appendix B.

Table 3 reports results for the ML-based reconciliation methods– default and
tuned –with SARIMA base forecasts. The random forest algorithm with default
tuning parameters, these were best performing in most of our previous results,
performs equally well as its tuned version. This is a great advantage when
implementing the forecasting procedure in a high-dimensional streaming data
environment with frequent updates. In contrast to random forest, the XGBoost
algorithm greatly benefits from tuning at the high-frequency delivery area
level. Here, the loss decreases from 0.3196 to 0.2790, a reduction of about
13% and XGBoost is now at the same accuracy level as random forest. For the
zone and market hierarchies, the differences between default and tuning are
marginal. Finally, the LightGBM algorithm is only mildly impacted, with small
tuning improvements at the delivery area level.

Table 3:  London forecast reconciliation results with hyperparameter tuning.

Forecast Method | Time | Delivery Area | Zone | Market
---|---|---|---|---
|  | (seconds) |  |  |  |  |  |  |  |  |
|  |  | 30min | hour | day | 30min | hour | day | 30min | hour | day
Random forest | default | 7 | 0.2801 | 0.2206 | 0.1002 | 0.1305 | 0.1132 | 0.0657 | 0.0883 | 0.0861 | 0.0586
Random forest | tuned | 400 | 0.2784 | 0.2198 | 0.1033 | 0.1324 | 0.1155 | 0.0692 | 0.0911 | 0.0889 | 0.0621
XGBoost | default | 1 | 0.3195 | 0.2528 | 0.1109 | 0.1385 | 0.1182 | 0.0659 | 0.0845 | 0.0815 | 0.0551
XGBoost | tuned | 600 | 0.2790 | 0.2209 | 0.1032 | 0.1314 | 0.1144 | 0.0652 | 0.0886 | 0.0864 | 0.0568
LightGBM | default | 1 | 0.2933 | 0.2327 | 0.1067 | 0.1320 | 0.1134 | 0.0641 | 0.0832 | 0.0806 | 0.0538
LightGBM | tuned | 900 | 0.2801 | 0.2221 | 0.1034 | 0.1334 | 0.1165 | 0.0668 | 0.0917 | 0.0895 | 0.0589

Notes: This table shows forecast accuracy measured in weighted absolute
percentage error (WAPE) for SARIMA base forecasts. The forecast horizon is one
week-ahead. Computing times are approximate, measured in seconds and recorded
on a Windows 10 Enterprise LTSC machine with Intel Core i5 2.4GHz processor.

While hyperparameter tuning does result, overall, in a marginal gain in
forecast accuracy, it is important to investigate the computational burden
that comes along with it. Table 3 reports in the third column the approximate
computing times (in seconds) for one ML fit with fixed tuning parameters as
well as one round of hyperparameter tuning, and this for random forest,
XGBoost and LightGBM. As well known, both the XGBoost and LightGBM algorithms
are fast with default tuning parameters, requiring only about one second to
train. Random forest is in comparison relatively slow with seven seconds to
train. However, the prohibitive cost of tuning is appalling. All algorithms
need large multiples of computing time to tune the hyperparameters. In
particular the XGBoost and LightGBM algorithms with many tuning
hyperparameters are long to optimize.666We use Bayesian optimization to tune
the hyperparameters of XGBoost and LightGBM. Other methods or only subsets of
parameters could be tuned. This would make the tuning faster.

To conclude, among the considered ML algorithms, the random forest algorithm
is preferred because the default tuning parameters yield good accuracy levels
close to the tuned ones. This is a great advantage when faced with
computational constraints such as given by the platform environments we
consider in this paper.

###  5.4 Data Shifts: The Case of Manchester

The streaming platform data for the London market are reasonably stable over
time. In this section, we investigate how our ML-based forecast reconciliation
procedure performs when facing shifts in the data streams. In fact, breaks or
data shifts are regularly observed in platform data because of new competitors
in the marketplace and companies moving some of their locations to other on-
demand platforms.

Figure 8: Time series plots for full sample of area Bolton Central (top)
belonging to the market Manchester (bottom), with a data shift for Bolton
Central at June 14 2023.

To study the impact of breaks on forecast reconciliation, we focus on
Manchester for which we have a cross-temporal hierarchy that consists cross-
sectionally of 24 delivery area and Manchester as the entire market (there is
no intermediate level consisting of zones in Manchester), and temporally of
30-minute, hourly and daily frequencies for the same sample period as London.
Figure A.2 in Appendix A displays aggregated demand over the sample period in
each delivery area in Manchester. Figure 8, top panel, illustrates how severe
the data shift is on June 14, 2023 in all considered temporal frequencies for
the Bolton Central delivery area situated North-West of Manchester.

Table 4:  Manchester forecast reconciliation results.

Period | Forecast | Delivery Area | Market
---|---|---|---
| Method |  |  |  |  |  |
|  | 30min | hour | day | 30min | hour | day
Before shift | Base | 0.3549 | 0.2922 | 0.1100 | 0.1884 | 0.1832 | 0.0815
| Bottom Up | 0.3549 | 0.2945 | 0.1879 | 0.1949 | 0.1905 | 0.1630
| tcs | 0.3459 | 0.2837 | 0.1734 | 0.1847 | 0.1794 | 0.1533
| cst | 0.3435 | 0.2807 | 0.1680 | 0.1794 | 0.1735 | 0.1460
| ite | 0.3463 | 0.2841 | 0.1743 | 0.1855 | 0.1800 | 0.1549
| oct | 0.3470 | 0.2846 | 0.1731 | 0.1848 | 0.1797 | 0.1520
| Random forest | 0.3226 | 0.2523 | 0.1146 | 0.1216 | 0.1141 | 0.0764
| XGBoost | 0.3650 | 0.2885 | 0.1315 | 0.1205 | 0.1121 | 0.0743
| LightGBM | 0.3418 | 0.2706 | 0.1300 | 0.1196 | 0.1118 | 0.0759
During shift | Base | 0.8124 | 0.7687 | 0.4634 | 0.2091 | 0.2034 | 0.0615
| Bottom Up | 0.8124 | 0.7587 | 0.5923 | 0.2092 | 0.2059 | 0.1805
| tcs | 0.7991 | 0.7683 | 0.5926 | 0.1965 | 0.1926 | 0.1642
| cst | 0.7986 | 0.7651 | 0.5747 | 0.1885 | 0.1841 | 0.1520
| ite | 0.8114 | 0.7730 | 0.6021 | 0.1999 | 0.1958 | 0.1679
| oct | 0.7885 | 0.7606 | 0.5789 | 0.1910 | 0.1875 | 0.1581
| Random forest | 0.8004 | 0.7346 | 0.5885 | 0.1110 | 0.1022 | 0.0646
| XGBoost | 0.8330 | 0.7534 | 0.5683 | 0.1128 | 0.1017 | 0.0633
| LightGBM | 0.7516 | 0.6818 | 0.5230 | 0.1101 | 0.1002 | 0.0658
After shift | Base | 0.7166 | 0.6402 | 0.1610 | 0.1776 | 0.1740 | 0.0524
| Bottom Up | 0.7166 | 0.6464 | 0.4506 | 0.1825 | 0.1786 | 0.1526
| tcs | 0.6969 | 0.6446 | 0.3694 | 0.1682 | 0.1631 | 0.1364
| cst | 0.7023 | 0.6484 | 0.3381 | 0.1615 | 0.1558 | 0.1256
| ite | 0.7049 | 0.6531 | 0.3617 | 0.1699 | 0.1647 | 0.1381
| oct | 0.6925 | 0.6402 | 0.3755 | 0.1650 | 0.1601 | 0.1322
| Random forest | 0.4206 | 0.3300 | 0.1457 | 0.0992 | 0.0904 | 0.0465
| XGBoost | 0.4863 | 0.3874 | 0.1464 | 0.1007 | 0.0896 | 0.0485
| LightGBM | 0.4386 | 0.3488 | 0.1551 | 0.0980 | 0.0880 | 0.0474

Notes: This table shows forecast accuracy measured in weighted absolute
percentage error (WAPE) for SARIMA base forecasts. The forecast horizon is one
week-ahead. The best forecast reconciliation results for the periods Before
shift (prior to June 14, 2023), After shift (after July 14, 2023), and During
shift (period in between) are highlighted in gray. The linear benchmark
methods are defined in Table 1.

Table 4 summarizes the forecast accuracy for the cross-temporal hierarchy for
SARIMA base forecasts and all reconciliation methods; the other base forecast
results are similar and available on request. We report results for the period
before the data shift (i.e. prior to June 14), the period after the data shift
(i.e. after July 14), and the one month period in between characterizing the
data shift. There are several interesting findings for the delivery area level
of the hierarchy. First, for “Before shift” period we see that random forest
reconciliation for delivery areas outperforms other methods with the same
relative magnitudes as found in the London market where no data shifts occur.
For example, random forest and cst WAPEs are 0.3226 and 0.3435, respectively.
Second, we see that the volatile “During shift” period WAPEs at the delivery
area level all more than double compared to the stable period before. The
LightGBM ML reconciliation method performs best with a WAPE of 0.7516. Third,
the “After shift” period is generally characterized by better accuracies than
during the shift, but the size of the gains are much stronger for the ML-based
reconciliation methods. To illustrate the magnitude of this, consider cst with
WAPEs of 0.7986 and 0.7023 and random forest WAPEs of 0.8004 and 0.4206,
respectively for the “During” and “After shift” periods. Fourth, the un-
reconciled “After shift” SARIMA base forecasts are of low quality with WAPEs
still double the size compared to the “Before shift” period. This can somehow
be expected given that parameter estimates are based on historical data
subject to the shift. However, our ML-based reconciliation method is able to
transform these low quality projections into forecasts that yield WAPEs much
closer to the period prior to the shift. At the market level part of the
hierarchy, we observe only small differences in terms of forecasting accuracy
over the three periods. This is expected since the aggregation of delivery
areas not all being subject to a data shift, results in a smoother market
level time series, see Figure 8, bottom panel. The LightGBM method has
consistently the smallest WAPEs at the 30-minute and hourly frequency levels.

##  6 Conclusion

Forecast reconciliation based on linear methods has been shown to be very
useful and successful in a wide variety of applications. Recently, ML-based
forecast reconciliation methods are proposed for cross-sectional time series
hierarchies. We extend the latter approach to cross-temporal hierarchies. We
provide a general description of the methodology and an application to a
unique platform dataset from an on-demand logistics company. In particular,
the data consists of demand streams that constitute a rich hierarchy in
spatial and time dimensions, with bottom level time series defined as delivery
areas at 30-minute frequency.

Our key empirical finding is that ML-based forecast reconciliation for our
platform data can result in substantial forecast accuracy improvements
compared to existing linear reconciliation methods. However, unless one can
face the huge computational cost of hyperparameter tuning, ML-based
reconciliation is not uniformly superior to linear methods. Another key
finding is that, as platform data streams are potentially impacted by data
shifts, our approach is able to react swiftly to such instability when
compared with linear approaches.

Several extensions should be considered. First, applications to other datasets
consisting of different cross-temporal hierarchies. Second, our ML algorithms
in the platform application consist of tree-based methods only, though our
approach can be used with any ML algorithm. Taking into account computing
time, it would be interesting to test alternative neural network based methods
using our approach. Third, we currently retrain the ML algorithms at every
iteration of the rolling window, it would be interesting to investigate if
online ML methods that incrementally update the forecast function as new data
becomes available can lead to further computational efficiencies.

## References

  * Abolghasemi et al., (2022) Abolghasemi, M., Tarr, G., and Bergmeir, C. (2022).  Machine learning applications in hierarchical time series forecasting: Investigating the impact of promotions.  International Journal of Forecasting, forthcoming.
  * Anderer and Li, (2022) Anderer, M. and Li, F. (2022).  Hierarchical forecasting with a top-down alignment of independent-level forecasts.  International Journal of Forecasting, 38(4):1405–1414.
  * Athanasopoulos et al., (2009) Athanasopoulos, G., Ahmed, R. A., and Hyndman, R. J. (2009).  Hierarchical forecasts for Australian domestic tourism.  International Journal of Forecasting, 25(1):146–166.
  * Athanasopoulos et al., (2023) Athanasopoulos, G., Hyndman, R. J., Kourentzes, N., and Panagiotelis, A. (2023).  Forecast reconciliation: A review.  International Journal of Forecasting, forthcoming.
  * Athanasopoulos et al., (2024) Athanasopoulos, G., Hyndman, R. J., Kourentzes, N., and Panagiotelis, A. (2024).  Editorial: Innovations in hierarchical forecasting.  International Journal of Forecasting, forthcoming.
  * Athanasopoulos et al., (2017) Athanasopoulos, G., Hyndman, R. J., Kourentzes, N., and Petropoulos, F. (2017).  Forecasting with temporal hierarchies.  European Journal of Operational Research, 262(1):60–74.
  * Bischl et al., (2016) Bischl, B., Lang, M., Kotthoff, L., Schiffner, J., Richter, J., Studerus, E., Casalicchio, G., and Jones, Z. M. (2016).  mlr: Machine learning in R.  Journal of Machine Learning Research, 17(170):1–5.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45:5–32.
  * Caporin et al., (2023) Caporin, M., Di Fonzo, T., and Girolimetto, D. (2023).  Exploiting intraday decompositions in realized volatility forecasting: A forecast reconciliation approach.  arXiv preprint arXiv:2306.02952.
  * Chen and Guestrin, (2016) Chen, T. and Guestrin, P. C. (2016).  XGBoost: A scalable tree boosting system.  in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. New York: ACM, 0(0):785–794.
  * Chen et al., (2023) Chen, T., He, T., Benesty, M., Khotilovich, V., Tang, Y., Cho, H., Chen, K., Mitchell, R., Cano, I., Zhou, T., Li, M., Xie, J., Lin, M., Geng, Y., Li, Y., and Yuan, J. (2023).  xgboost: Extreme Gradient Boosting.  R package version 1.7.6.1.
  * (12) Di Fonzo, T. and Girolimetto, D. (2023a).  Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives.  International Journal of Forecasting, 39(1):39–57.
  * (13) Di Fonzo, T. and Girolimetto, D. (2023b).  Spatio-temporal reconciliation of solar forecasts.  Solar Energy, 251:13–29.
  * Efron and Hastie, (2016) Efron, B. and Hastie, T. (2016).  Computer Age Statistical Inference.  Cambridge University Press, Cambridge.
  * Elliott and Timmermann, (2016) Elliott, G. and Timmermann, A. (2016).  Economic Forecasting.  Princeton University Press, Princeton, New Jersey.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  The Annals of Statistics, 29(5):1189–1232.
  * Girolimetto and Di Fonzo, (2023) Girolimetto, D. and Di Fonzo, T. (2023).  FoReco: Point Forecast Reconciliation.  R package version 2.6.
  * Hollyman et al., (2021) Hollyman, R., Petropoulos, F., and Tipping, M. E. (2021).  Understanding forecast reconciliation.  European Journal of Operational Research, 294(1):149–160.
  * Hyndman et al., (2023) Hyndman, R., Athanasopoulos, G., Bergmeir, C., Caceres, G., Chhay, L., O’Hara-Wild, M., Petropoulos, F., Razbash, S., Wang, E., and Yasmeen, F. (2023).  forecast: Forecasting functions for time series and linear models.  R package version 8.21.1.
  * Hyndman et al., (2011) Hyndman, R. J., Ahmed, R. A., Athanasopoulos, G., and Shang, H. L. (2011).  Optimal combination forecasts for hierarchical time series.  Computational Statistics & Data Analysis, 55(9):2579–2589.
  * Hyndman and Athanasopoulos, (2014) Hyndman, R. J. and Athanasopoulos, G. (2014).  Optimally reconciling forecasts in a hierarchy.  Foresight: The International Journal of Applied Forecasting, 35:42–48.
  * Hyndman and Athanasopoulos, (2021) Hyndman, R. J. and Athanasopoulos, G. (2021).  Forecasting: principles and practice.  3rd. Melbourne, Australia: OTexts.
  * Hyndman and Khandakar, (2008) Hyndman, R. J. and Khandakar, Y. (2008).  Automatic time series forecasting: The forecast package for R.  Journal of Statistical Software, 26(3):1–22.
  * Hyndman and Koehler, (2006) Hyndman, R. J. and Koehler, A. B. (2006).  Another look at measures of forecast accuracy.  International Journal of Forecasting, 22(4):679–688.
  * Hyndman et al., (2002) Hyndman, R. J., Koehler, A. B., Snyder, R. D., and Grose, S. (2002).  A state space framework for automatic forecasting using exponential smoothing methods.  International Journal of forecasting, 18(3):439–454.
  * Hyndman et al., (2016) Hyndman, R. J., Lee, A. J., and Wang, E. (2016).  Fast computation of reconciled forecasts for hierarchical and grouped time series.  Computational Statistics & Data Analysis, 97:16–32.
  * Januschowski et al., (2022) Januschowski, T., Wang, Y., Torkkola, K., Erkkilä, T., Hasson, H., and Gasthaus, J. (2022).  Forecasting with trees.  International Journal of Forecasting, 38(4):1473–1481.
  * Komisarczyk et al., (2023) Komisarczyk, K., Kozminski, P., Maksymiuk, S., and Biecek, P. (2023).  treeshap: Compute SHAP values for your tree-based models using the ‘TreeSHAP’ algorithm.  R package version 0.2.5.
  * Koning et al., (2005) Koning, A. J., Franses, P. H., Hibon, M., and Stekler, H. (2005).  The M3 competition: Statistical tests of the results.  International Journal of Forecasting, 21(3):397–409.
  * Kourentzes, (2023) Kourentzes, N. (2023).  tsutils: Time Series Exploration, Modelling and Forecasting.  R package version 0.9.4.
  * Kourentzes and Athanasopoulos, (2019) Kourentzes, N. and Athanasopoulos, G. (2019).  Cross-temporal coherent forecasts for Australian tourism.  Annals of Tourism Research, 75:393–409.
  * Liaw and Wiener, (2002) Liaw, A. and Wiener, M. (2002).  Classification and regression by randomforest.  R News, 2(3):18–22.
  * Lundberg et al., (2018) Lundberg, S. M., Erion, G. G., and Lee, S.-I. (2018).  Consistent individualized feature attribution for tree ensembles.  arXiv preprint arXiv:1802.03888.
  * Lundberg and Lee, (2017) Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  Advances in neural information processing systems, 30.
  * Makridakis et al., (2020) Makridakis, S., Spiliotis, E., and Assimakopoulos, V. (2020).  The M4 competition: 100,000 time series and 61 forecasting methods.  International Journal of Forecasting, 36(1):54–74.
  * (36) Makridakis, S., Spiliotis, E., and Assimakopoulos, V. (2022a).  M5 accuracy competition: Results, findings, and conclusions.  International Journal of Forecasting, 38(4):1346–1364.
  * (37) Makridakis, S., Spiliotis, E., and Assimakopoulos, V. (2022b).  The M5 competition: Background, organization, and implementation.  International Journal of Forecasting, 38(4):1325–1336.
  * (38) Makridakis, S., Spiliotis, E., and Assimakopoulos, V. (2022c).  Predicting/hypothesizing the findings of the M5 competition.  International Journal of Forecasting, 38(4):1337–1345.
  * Medeiros et al., (2021) Medeiros, M. C., Vasconcelos, G. F. R., Álvaro Veiga, and Zilberman, E. (2021).  Forecasting inflation in a data-rich environment: The benefits of machine learning methods.  Journal of Business & Economic Statistics, 39(1):98–119.
  * Molnar, (2022) Molnar, C. (2022).  Interpretable Machine Learning.  2 edition.
  * Panagiotelis et al., (2021) Panagiotelis, A., Athanasopoulos, G., Gamakumara, P., and Hyndman, R. J. (2021).  Forecast reconciliation: A geometric view with new insights on bias correction.  International Journal of Forecasting, 37(1):343–359.
  * Petropoulos et al., (2022) Petropoulos, F., Apiletti, D., Assimakopoulos, V., Babai, M. Z., Barrow, D. K., Taieb, S. B., Bergmeir, C., Bessa, R. J., Bijak, J., Boylan, J. E., et al. (2022).  Forecasting: theory and practice.  International Journal of Forecasting, 38(3):705–871.
  * R Core Team, (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  * Rostami-Tabar and Hyndman, (2024) Rostami-Tabar, B. and Hyndman, R. J. (2024).  Hierarchical time series forecasting in emergency medical services.  Journal of Service Research, forthcoming.
  * Shapley et al., (1953) Shapley, L. S. et al. (1953).  A value for n-person games.  Contributions to the theory of games, 2:307–317.
  * Shi et al., (2023) Shi, Y., Ke, G., Soukhavong, D., Lamb, J., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.-Y., and Titov, N. (2023).  lightgbm: Light Gradient Boosting Machine.  R package version 3.3.5.
  * Spiliotis et al., (2021) Spiliotis, E., Abolghasemi, M., Hyndman, R. J., Petropoulos, F., and Assimakopoulos, V. (2021).  Hierarchical forecast reconciliation with machine learning.  Applied Soft Computing, 112:107756.
  * Theodosiou and Kourentzes, (2021) Theodosiou, F. and Kourentzes, N. (2021).  Forecasting with deep temporal hierarchies.  Available at SSRN 3918315.
  * Van Erven and Cugliari, (2015) Van Erven, T. and Cugliari, J. (2015).  Game-theoretically optimal reconciliation of contemporaneous hierarchical time series forecasts.  In Modeling and stochastic learning for forecasting in high dimensions, pages 297–317. Springer.
  * Wang et al., (2022) Wang, S., Zhou, F., Sun, Y., Ma, L., Zhang, J., and Zheng, Y. (2022).  End-to-end modeling of hierarchical time series using autoregressive transformer and conditional normalizing flow-based reconciliation.  In 2022 IEEE International Conference on Data Mining Workshops (ICDMW), pages 1087–1094. IEEE.
  * Wickramasuriya et al., (2019) Wickramasuriya, S. L., Athanasopoulos, G., and Hyndman, R. J. (2019).  Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization.  Journal of the American Statistical Association, 114(526):804–819.
  * Yan, (2021) Yan, Y. (2021).  rBayesianOptimization: Bayesian Optimization of Hyperparameters.  R package version 1.2.0.

##  Appendix A Platform Data

Figure A.1: Time series plots for 30-min, hourly and daily demand (in the
rows) for the area City of London belonging to the zone London Central and the
market London (in the columns). The black solid line is the pointwise median.
The shaded blue area displays the 5% and 95% pointwise quantiles. Figure A.2:
Aggregated demand in each delivery area in Manchester.

##  Appendix B Design of the Forecast Study

###  B.1 Rolling-Window Set-Up

iteration $\displaystyle r$ iteration $\displaystyle r+1$ Figure B.1: Outer
rolling window in forecast set-up and inner rolling windows procedure in the
machine-learning based forecast reconciliation.

###  B.2 Base Forecasts

We consider four base forecast models: (i) Naive, (ii) SARIMA, (iii) ETS and a
(iv) forecast combination of the former three.

Naive. The first base forecast model for temporal aggregation level $k$ of a
specific cross-sectional time series is simply
$x_{j}^{[k]}=x_{j-7m/k}^{[k]}+\varepsilon^{[k]}_{j}$, which we call “Naive”
since it simply implies that the forecast of a specific time slot (30-min,
hourly or daily) is the value observed in the same time slot and day of the
previous week. This way of forecasting is popular in the industry since it is
ultra fast, i.e. it does not require parameter estimation, and it works well
in the case of strong seasonality patterns as observed in platform data.

SARIMA. Secondly, we consider ARIMA models, as implemented in the forecast
package (Hyndman et al.,, 2023; Hyndman and Khandakar,, 2008) in R. For the
daily time series, we search for each series the optimal model in the space of
seasonal ARIMA models, thereby allowing for weekly seasonality (seasonal
period equal to 7 for daily data) and a maximal MA order $q=3$, and maximal AR
order of $p=3$. For the seasonal components we set the maximal MA order to
$Q=1$ and maximal AR order to $P=1$. We restrict the maximum number of non-
seasonal differences to $d=1$. All other arguments are kept at their default
values.

For the hourly and 30-min time series, we need to additionally account for the
strong intra-day seasonal patterns in the platform data. Since the seasonal
period is long in both cases ($m_{1}=34$ for 30-min data and $m_{2}=17$ for
hourly data), we follow the common practice based on Fourier series to capture
the seasonality. Specifically, for the 30-min time series $y_{t}$, we consider
the regression model

$y_{t}=\mu+\sum_{s=1}^{S}\left[\alpha_{s}\text{sin}(2\pi
st/m_{1})+\beta_{s}\text{cos}(2\pi st/m_{1})\right]+\varepsilon_{t},$ (B.1)

where $\varepsilon_{t}$ is subsequently modeled as an ARIMA process.
Analogously for the hourly data thereby taking the seasonal period $m_{2}=17$.
We select the optimal value of $S$ via the Bayesian Information Criterion.

ETS. Thirdly, we consider exponential smoothing (ETS; Hyndman et al.,, 2002)
as another popular model to produce the base forecasts. To this end, we use
the ets function (with default arguments) as implemented in the forecast
package. As above, for the 30-min and hourly series we first estimate the
intra-day seasonality through the Fourier series approach and subsequently
model the error term in equation (B.1) using ETS.

Forecast Combination. Our final base forecast model consists of a simple
forecast combination of the previous three base forecast models. Forecast
combination approaches have shown to perform well as base forecast methods in
the hierarchical forecasting literature, thereby oftentimes providing an
effective practice for improving forecast accuracy; see for instance Rostami-
Tabar and Hyndman, (2024) for an application on emergency medical services
data. We opt for equal weights in the forecast combination, mainly because of
its simplicity, ease in implementation and its established track-record of
good performance in the forecast combination literature (e.g., Elliott and
Timmermann,, 2016, Chapter 14).

###  B.3 Machine Learning Methods

We consider three machine learning methods: (i) random forest, (ii) XGBoost
and (iii) LightGBM.

Random Forest. Random forests, proposed by Breiman, (2001), produce forecasts
by combining regression trees. A regression tree is a nonparametric method
that partitions the feature space to compute local averages as forecasts, see
Efron and Hastie, (2016) for a textbook treatment. The tuning parameters are
the number of trees that are used in the forecast combination, the number of
features to randomly select when constructing each regression tree split, and
the minimum number of observations in each terminal node to compute the local
forecasts. We use the standard implementation of the randomForest package
(Liaw and Wiener,, 2002) in `R` with default settings for the hyperparamters
(i.e. number of trees: ntree = 500, number of variables sampled at each split:
$\texttt{mtry}=\text{\\# of features}/3$, minimum size of terminal nodes:
nodesize = 5) when reporting our main results.

In Section 5.3, we also investigate the performance of the random forest based
forecast reconciliation method when the hyperparameters are tuned. To this
end, we follow Spiliotis et al., (2021) and use the package mlr in R (Bischl
et al.,, 2016) which implements a random grid search with cross-validation to
find the optimal hyperparameters. In the cross-validation procedure, we leave
each of the four subsequent weeks in the validation set once out as test
sample and use the root mean squared error as cross-validation score to
compute the optimal hyperparameters. Note that we tune the hyperparameters
only once every four iterations (i.e. once every month) in the outer rolling
window (so in 8 out of the 31 outer rolling windows) to keep the computational
burden low; in between we use the optimal hyperparameters from the previous
tuning round as fixed values. The lower and upper bounds for the
hyperparameters are set to $(2,50)$ for mtry and $(5,50)$ for nodesize. We set
the bounds for ntree between 50 and 500 on interval steps of 10.

XGBoost. Gradient boosting, proposed by Friedman, (2001), constructs forecasts
by sequentially fitting small regression trees, i.e. weak learners, to the
residuals by the ensemble of the previous trees. This procedure results in a
one final tree, constructed as a sum of trees, used for forecasting. Extreme
gradient boosting (XGBoost), introduced by Chen and Guestrin, (2016),
optimizes the implementation of the gradient boosting framework in terms of
speed and flexibility. We use the xgboost package (Chen et al.,, 2023) in `R`
with fixed choices of the tuning parameters when reporting our main results.
We fix the hyperparameters to the default values as follows: 100 boosting
iterations (nrounds), 6 as max tree depth (max_depth), 0.3 as learning rate
(eta), 1 as subsample ratio (subsample), 1 as subsample ratio of columns
(colsample_bytree), 1 as minimum sum of instance weight (hessian)
(min_child_weight) and 0 as minimum loss reduction (gamma).

Subsequently, to tune the hyperparameters, we use the same procedure as
discussed for random forest but this time use Bayesian optimization to tune
the hyperparameters since grid search is computationally very expensive in
this case. To this end, we use the rBayesianOptimization (Yan,, 2021) in R and
consider the following intervals with lower and upperbounds for each
hyperparameter, in line with Spiliotis et al., (2021): we set the values of
max_depth between (2, 10), the learning rate (eta) between (0.01, 0.05),
subsample values between (0.3, 1), colsample_bytree values between (0.3, 1),
min_child_weight between (0, 10) and gamma between (0, 5). The values for the
maximum number of boosting iterations (nrounds) are over the range of 50 and
200.

LightGBM. LightGBM, put forward by Microsoft in 2016, is a gradient boosting
framework that uses tree-based learning algorithms like XGBoost but as the
name suggests has computational advantages with respect to training speed,
memory usage and parallelization. The implementation of gradient-based one-
side sampling and exclusive feature bundling techniques allows handling large
training datasets. We use the lightgbm package (Shi et al.,, 2023) in `R` with
fixed, default hyperparameters when reporting our main results: 100 boosting
iterations (nrounds), 31 as maximum number of leaves (num_leaves), 0.1 as
learning rate (eta), 1 as subsample ratio (subsample), 1 as subsample ratio of
columns (colsample_bytree), 0.001 as minimum sum of instance weight (hessian)
(min_child_weight) and 0 as $\ell_{1}$-regularization (lambda_l1) and no limit
to the max depth (max_depth) for the tree model.

Finally, to tune the hyperparameters, we also use Bayesian optimization with
the following similar prior values for the hyperparameters as with XGBoost.
The maximum number of leaves (num_leaves) fixed to the range of 5 and 31. The
learning rate (eta) is set between (0.01, 0.05), subsample values between
(0.3, 1), colsample_bytree values are set between (0.3, 1), min_child_weight
between (0, 10), max_depth between (2, 10) and lambda_l1 between (0, 5). The
values for the maximum number of boosting iterations (nrounds) are over the
range of 50 and 200.

##  Appendix C Platform Application: Additional Results

This Appendix contains the forecast accuracy results on the London data for
two additional forecast error metrics, namely the symmetric mean absolute
percentage error (SMAPE) and the root mean squared scaled error (RMSSE), see
Hyndman and Koehler, (2006) for an overview.

The two forecast error metrics, for temporal factor $k$ and cross-sectional
time series $i=1,\ldots,n$, are defined as

$\text{SMAPE}^{[k]}_{i}=\frac{1}{T_{\text{test}}/k}\sum_{j=1}^{T_{\text{test}}/k}\frac{|A_{i,j}^{[k]}-F_{i,j}^{[k]}|}{|A_{i,j}^{[k]}|+|F_{i,j}^{[k]}|},$

and

$\text{RMSSE}_{i,r}^{[k]}=\sqrt{\left[\frac{1}{m_{k}H}\sum\limits_{j=m_{k}N+1}^{m_{k}(N+H)}(A_{i,j}^{[k]}-F^{[k]}_{i,j})^{2}\right]/\left[\frac{1}{m_{k}(N-7)}\sum\limits_{j=7m_{k}+1}^{m_{k}N}(A_{i,j}^{[k]}-A_{i,j-7m_{k}}^{[k]})^{2}\right]},$

for outer rolling window $r=1,\ldots,31$, and
$\text{RMSSE}_{i}^{[k]}=\frac{1}{31}\sum_{r=1}^{31}\text{RMSSE}_{i,r}^{[k]}$
given the 31 outer rolling windows in our forecast set-up. Table C.1 presents
the results for the SMAPE, Table C.2 for the RMSSE.

Table C.1:  London forecast reconciliation results for SMAPE loss.

Base | Forecast | Delivery Area | Zone | London
---|---|---|---|---
Forecasts | Method |  |  |  |  |  |  |  |  |
|  | 30min | hour | day | 30min | hour | day | 30min | hour | day
naive | Base | 22.5751 | 16.9011 | 5.2106 | 8.5895 | 6.7946 | 3.2439 | 4.2526 | 4.0147 | 2.8348
| Random forest | 18.0685 | 13.9355 | 5.3208 | 7.3302 | 6.0777 | 3.2188 | 4.4706 | 4.2126 | 2.9165
| XGBoost | 19.2921 | 14.8049 | 5.2141 | 7.7537 | 6.3345 | 3.3164 | 4.5151 | 4.2266 | 2.9741
| LightGBM | 18.6969 | 14.4084 | 5.3153 | 7.5802 | 6.2178 | 3.2874 | 4.5408 | 4.2495 | 2.9651
ETS | Base | 20.7684 | 17.2927 | 4.2926 | 11.4502 | 10.2841 | 2.7382 | 9.4232 | 8.4256 | 2.3900
| Bottom Up | 20.7684 | 17.3651 | 9.1734 | 10.5947 | 9.6581 | 6.8178 | 8.2864 | 8.0427 | 6.4905
| tcs | 20.2981 | 16.7859 | 7.9347 | 10.0560 | 9.1228 | 6.3790 | 7.7067 | 7.5057 | 6.0820
| cst | 20.1352 | 16.6636 | 7.7101 | 9.6296 | 8.7143 | 5.8954 | 7.1204 | 6.9684 | 5.5321
| ite | 20.4516 | 17.0330 | 8.1244 | 10.2236 | 9.3600 | 6.5371 | 7.8553 | 7.7311 | 6.2203
| oct | 20.2486 | 16.6548 | 7.7160 | 10.0888 | 9.1167 | 6.2570 | 7.7516 | 7.4954 | 5.9571
| Random forest | 18.4454 | 14.4383 | 5.7600 | 7.9979 | 6.8666 | 3.6242 | 5.4889 | 5.2527 | 3.3360
| XGBoost | 20.2685 | 16.1830 | 6.2633 | 8.4360 | 7.1863 | 3.6507 | 5.2181 | 4.9604 | 3.1321
| LightGBM | 19.1597 | 15.2385 | 6.2077 | 8.1262 | 6.9507 | 3.6457 | 5.2457 | 4.9954 | 3.1812
SARIMA | Base | 20.0513 | 16.4029 | 4.5864 | 10.1443 | 9.1848 | 2.8773 | 8.0568 | 7.6200 | 2.6608
| Bottom Up | 20.0513 | 16.4626 | 8.5675 | 10.2068 | 9.2690 | 6.6271 | 8.2048 | 7.9624 | 6.3712
| tcs | 19.7353 | 16.0791 | 7.4439 | 9.7834 | 8.8115 | 6.2130 | 7.6081 | 7.3829 | 5.9539
| cst | 19.5960 | 15.9584 | 7.1493 | 9.4172 | 8.5024 | 5.7717 | 7.0327 | 6.8877 | 5.4519
| ite | 19.8005 | 16.2166 | 7.5748 | 9.8460 | 8.9688 | 6.3424 | 7.6368 | 7.5063 | 6.0710
| oct | 19.6591 | 15.9352 | 7.2266 | 9.6154 | 8.6084 | 5.9998 | 7.4799 | 7.2344 | 5.7753
| Random forest | 18.1941 | 14.1153 | 5.3504 | 7.5216 | 6.3654 | 3.3066 | 4.8442 | 4.6770 | 2.9467
| XGBoost | 20.0262 | 15.6200 | 5.6302 | 7.9671 | 6.6464 | 3.3281 | 4.5853 | 4.3811 | 2.7765
| LightGBM | 18.8722 | 14.7949 | 5.6537 | 7.5975 | 6.3711 | 3.2282 | 4.5371 | 4.3372 | 2.7075
Forecast Combination | Base | 19.1283 | 15.0267 | 4.3534 | 8.6499 | 7.4707 | 2.7334 | 6.2835 | 5.8143 | 2.4571
| Bottom Up | 19.1283 | 15.1227 | 6.7835 | 8.5620 | 7.4380 | 4.9205 | 6.0844 | 5.8620 | 4.6645
| tcs | 19.0291 | 14.9696 | 6.1030 | 8.3911 | 7.2776 | 4.7630 | 5.8247 | 5.6441 | 4.5212
| cst | 18.9747 | 14.9280 | 5.9454 | 8.2179 | 7.1169 | 4.4956 | 5.5536 | 5.4126 | 4.2237
| ite | 19.1102 | 15.0884 | 6.2119 | 8.4948 | 7.4385 | 4.8716 | 5.9654 | 5.8388 | 4.6208
| oct | 18.9531 | 14.8599 | 5.9487 | 8.2957 | 7.1462 | 4.6062 | 5.7255 | 5.5037 | 4.3655
| Random forest | 17.9998 | 13.9051 | 5.4693 | 7.4518 | 6.2836 | 3.5961 | 4.8270 | 4.6312 | 3.3328
| XGBoost | 19.2742 | 14.8543 | 5.2093 | 7.7242 | 6.3802 | 3.3615 | 4.5575 | 4.3324 | 2.9994
| LightGBM | 18.6192 | 14.3817 | 5.2450 | 7.4903 | 6.2069 | 3.3260 | 4.5124 | 4.2996 | 2.9680

Notes: This table shows forecast accuracy measured in symmetric mean absolute
percentage error (SMAPE). The forecast horizon is one week-ahead. The best
forecast reconciliation results for each base forecast (Naive, ETS, SARIMA,
Forecast Combination) are highlighted in gray.

Table C.2:  London forecast reconciliation results for RMSSE loss.

Base | Forecast | Delivery Area | Zone | London
---|---|---|---|---
Forecasts | Method |  |  |  |  |  |  |  |  |
|  | 30min | hour | day | 30min | hour | day | 30min | hour | day
naive | Base | 0.9651 | 0.9353 | 0.7624 | 0.8704 | 0.8256 | 0.6609 | 0.7267 | 0.7111 | 0.6211
| Random forest | 0.7537 | 0.7590 | 0.7056 | 0.7533 | 0.7478 | 0.6485 | 0.7400 | 0.7242 | 0.6296
| XGBoost | 0.8242 | 0.8202 | 0.7333 | 0.8044 | 0.7856 | 0.6652 | 0.7647 | 0.7434 | 0.6414
| LightGBM | 0.7929 | 0.7923 | 0.7232 | 0.7830 | 0.7683 | 0.6609 | 0.7584 | 0.7367 | 0.6382
ETS | Base | 0.8920 | 0.9665 | 0.6217 | 1.1561 | 1.2359 | 0.5479 | 1.4382 | 1.4153 | 0.5102
| Bottom Up | 0.8920 | 0.9741 | 1.3343 | 1.1120 | 1.2014 | 1.4209 | 1.3668 | 1.3848 | 1.4444
| tcs | 0.8721 | 0.9438 | 1.2217 | 1.0748 | 1.1574 | 1.3197 | 1.3054 | 1.3230 | 1.3374
| cst | 0.8602 | 0.9277 | 1.1709 | 1.0385 | 1.1147 | 1.2210 | 1.2393 | 1.2559 | 1.2152
| ite | 0.8766 | 0.9520 | 1.2566 | 1.0886 | 1.1766 | 1.3688 | 1.3305 | 1.3515 | 1.3915
| oct | 0.8705 | 0.9396 | 1.1900 | 1.0700 | 1.1496 | 1.2943 | 1.2968 | 1.3124 | 1.3142
| Random forest | 0.7683 | 0.7848 | 0.7606 | 0.7961 | 0.8095 | 0.7215 | 0.8307 | 0.8284 | 0.7087
| XGBoost | 0.8663 | 0.8977 | 0.8670 | 0.8421 | 0.8475 | 0.7367 | 0.7987 | 0.7916 | 0.6807
| LightGBM | 0.8120 | 0.8379 | 0.8385 | 0.8106 | 0.8203 | 0.7351 | 0.7995 | 0.7949 | 0.6904
SARIMA | Base | 0.8694 | 0.9383 | 0.6629 | 1.0799 | 1.1639 | 0.5829 | 1.3213 | 1.3339 | 0.5648
| Bottom Up | 0.8694 | 0.9411 | 1.2200 | 1.0796 | 1.1623 | 1.3347 | 1.3242 | 1.3403 | 1.3666
| tcs | 0.8574 | 0.9219 | 1.1428 | 1.0560 | 1.1340 | 1.2666 | 1.2881 | 1.3041 | 1.3026
| cst | 0.8474 | 0.9083 | 1.0933 | 1.0261 | 1.0992 | 1.1819 | 1.2329 | 1.2491 | 1.1998
| ite | 0.8607 | 0.9282 | 1.1690 | 1.0661 | 1.1489 | 1.3036 | 1.3065 | 1.3263 | 1.3440
| oct | 0.8521 | 0.9130 | 1.1008 | 1.0386 | 1.1116 | 1.2087 | 1.2582 | 1.2721 | 1.2418
| Random forest | 0.7657 | 0.7773 | 0.7225 | 0.7741 | 0.7809 | 0.6561 | 0.7828 | 0.7788 | 0.6266
| XGBoost | 0.8630 | 0.8795 | 0.8090 | 0.8199 | 0.8130 | 0.6598 | 0.7516 | 0.7394 | 0.5874
| LightGBM | 0.8037 | 0.8205 | 0.7740 | 0.7800 | 0.7795 | 0.6381 | 0.7355 | 0.7279 | 0.5718
Forecast | Base | 0.8122 | 0.8371 | 0.6317 | 0.8963 | 0.9300 | 0.5502 | 1.0049 | 1.0034 | 0.5216
Combination | Bottom Up | 0.8122 | 0.8395 | 0.9569 | 0.8903 | 0.9246 | 0.9950 | 0.9930 | 0.9998 | 1.0022
| tcs | 0.8079 | 0.8321 | 0.9219 | 0.8821 | 0.9151 | 0.9704 | 0.9803 | 0.9879 | 0.9805
| cst | 0.8025 | 0.8241 | 0.8896 | 0.8641 | 0.8929 | 0.9136 | 0.9441 | 0.9509 | 0.9111
| ite | 0.8104 | 0.8365 | 0.9412 | 0.8903 | 0.9266 | 0.9999 | 0.9960 | 1.0057 | 1.0140
| oct | 0.8045 | 0.8263 | 0.8929 | 0.8709 | 0.8998 | 0.9302 | 0.9582 | 0.9637 | 0.9356
| Random forest | 0.7518 | 0.7576 | 0.7360 | 0.7574 | 0.7585 | 0.7098 | 0.7627 | 0.7558 | 0.6977
| XGBoost | 0.8215 | 0.8213 | 0.7444 | 0.7907 | 0.7766 | 0.6720 | 0.7437 | 0.7311 | 0.6403
| LightGBM | 0.7902 | 0.7914 | 0.7286 | 0.7679 | 0.7569 | 0.6594 | 0.7340 | 0.7203 | 0.6291

Notes: This table shows forecast accuracy measured in root mean squared scaled
error (RMSSE). The forecast horizon is one week-ahead. The best forecast
reconciliation results for each base forecast (Naive, ETS, SARIMA, Forecast
Combination) are highlighted in gray.