Hierarchical Time Series Forecasting with Robust Reconciliation
Authors: Shuhei Aikawa
Aru Suzuki
Kei Yoshitake
Kanata Teshigawara
Iwabuchi Akira
Ken Kobayashi
Kazuhide Nakata
Venue: Accepted by TMLR ()

Abstract: This paper focuses on forecasting hierarchical time-series data, where each higher-level observation equals the sum of its corresponding lower-level time series. In such contexts, the forecast values should be coherent, meaning that the forecast value of each parent series exactly matches the sum of the forecast values of its child series. Existing hierarchical forecasting methods typically generate base forecasts independently for each series and then apply a reconciliation procedure to adjust them so that the resulting forecast values are coherent across the hierarchy. These methods generally derive an optimal reconciliation, using a covariance matrix of the forecast error. In practice, however, the true covariance matrix is unknown and has to be estimated from finite samples in advance. This gap between the true and estimated covariance matrix may degrade forecast performance. To address this issue, we propose a robust optimization framework for hierarchical reconciliation that accounts for uncertainty in the estimated covariance matrix. We first introduce an uncertainty set for the estimated covariance matrix and formulate a reconciliation problem that minimizes the worst-case average of weighted squared residuals over this uncertainty set. We show that our problem can be cast as a semidefinite optimization problem. Numerical experiments demonstrate that the proposed robust reconciliation method achieved better forecast performance than existing hierarchical forecasting methods, which indicates the effectiveness of integrating uncertainty into the reconciliation process.

Published in Transactions on Machine Learning Research (02/2026)
Hierarchical Time Series Forecasting with Robust Reconcil- iation
Shuhei Aikawa
rikovo...66470m.isttt.
Deportment of Industrial Engineering and Economics Institute of Science Tokyo
Aru Suzuki
Benuk....3822mmmsscc ac.j.
Department Industrial Engineering and Economics Institute of Science Tokyo
Kei Yoshitake
woshitake.k.661100 www.issctt.... ascp.
Department of Industrial Engineering and Economics Institute of Science Tokyo
Kanata Teshigawara
Reshigrwarokk./6c, Darec4...iisee @m..iscti.......... ac.jp
Department Industrial Engineering and Economics Institute of Science Tokyo
Akira Iwabuchi
twabuchi. 3.55 &mm. isct. ac.jp
Department of Industrial Engineering and Economics Institute of Science Tokyo
Ken Kobayashi
bobayashi.küie..ce inw.bee..engis.s........ ac.jp
Department of Industrial Engineering and Economics Institute of Science Tokyo
Kazuhide Nakata
nakato. k. ac am www.ttittccc.... ac.jp
Department Industrial Engineering and Economics Institute of Science Tokyo
Reviewed on OpenReview: openreviee. net/ pert/formm? id= XHPL)F52gy
Abstract
This paper focuses on forecasting hierarchical time-series data, where each higher-level ob- servation equals the sum of its corresponding lower-level time series. In such contexts, the forecast values should be coherent, meaning that the forecast value of each parent series exactly matches the sum of the forecast values of its child series. Existing hierarchical fore- casting methods typically generate base forecasts independently for each series and then apply a reconciliation procedure to adjust them so that the resulting forecast values are co- herent across the hierarchy. These methods generally yield an optimal reconciliation, using covariance matrix of the forecast errors. In practice, however, the true covariance matrix unknown and has to be estimated from finite samples in advance. This gap between the true and estimated covariance matrix may degrade forecast performance. To address this issue, we propose robust optimization framework for hicrarchical reconciliation that accounts for uncertainty in the estimated covariance matrix. We first. induuee an uncertainly set for the estimated covariance matrix and formulate a reconciliation problem that minimizes the worst-case average of weighted squared residuals over this uncortainty set. We show that our problem can be cast as a semidefinite optimization problem. Numerical experiments demonstrate that the proposed robust reconciliation method achieved better forecast. per-

---

Published in Transactions on Machine Learning Research (02/2026)
formance than existing hicrarchical forecasting methods, which indicates the effectiveness of integrating uncertainty into the reconciliation process.
1 Introduction
Time-series forecasting is indispensable across diverse fields, including sales planning and inventory manage- ment (Aviv] 2003 Ramos al. 2015), energy supply planning (Suganthi & Samuel, 2012 Hernandez et al,, 2014), and economic analysis and stock investment decision-making (Krollner al. 2010). For instance, in the retail sector, accurate sales predictions based on historical data are crucial for optimizing inventory levels and preventing both overstocking and shortages. Similarly, for electric power companies, forecasting electricity consumption enables efficient facility operations and effective supply-demand balance manage- ment. Moreover, at both the individual and national levels, leveraging forecasts of economic indicators and stock prices can significantly contribute to wealth creation. Conversely, low forecast accuracy can lead to substantial losses and missed opportunities for individuals, businesses, and society as whole. Consequently, extensive research has focused on developing various time-series forecasting methods, and there remains a strong demand for more precise techniques (Mahalakshmi et al. 2016 Liu et all. 2021 Wen et al. 2022).
Many real-world datasets inherently possess hicrarchical structures. Examples include sales data organized by region or demographic statistics categorized by gender or age group, which are commonly recorded across multiple levels of aggregation. In practice, the appropriate hierarchical level for forecasting depends on the specific application, and this choice can significantly influence prediction outcomes. Generally, as one descends lower hicrarchical levels, the data bocomes more gramular but also more suscoptible to individual variations and noise, leading to increased uncertainly. This often means that aggregated data higher levels tends to be more stable and achieve greater forecast accuracy (Grunfeld & Griliches) 1960). However, higher- level aggregated data can obscure fine-grained patterns and individual variation factors. Therefore, it has also been suggested that utilizing detailed data from lower levels, if appropriately modeled, can potentially yield superior forecast accuracy (Orcutt et al. 1968 Edwards & Orcutt 1969).
Given this context, time-series forecasting methods that explicitly account for hierarchical structures have garnered increasing attention (Athanasopoulos et al.] 2009: Hyndman et al.) 2011 Wickramasuriya et al., 2019 Shiratori et al., 2020: Panagiotelis et al.] 2021 Hyndman & Athanasopoulos 2021p. These approaches
aim to adjust forecasts across both lower and higher levels to ensure coherence when aggregating forocast val- ues within the hierarchy. Such methods are expected to enhance forecast accuracy compared to conventional time-series forecasting based on a single level.
Despite these advancements, existing hicrarchical time-scries forecasting methods face certain challenges. Traditional approaches lypically aim to minimize the expected forecast error al a given target. time point, which necessitates covariance matrix of the forecast errors. This matrix is commonly estimated from the residuals between obscrved and forecast values. However, if underlying data trends shift the forecasting model is inaccurate, discrepancies the estimated covariance matrix itself carries inherent uncertainty, which must be addressed. Prior work by Møller et al. (2024) focused on this issue, decomposing covariance matrix estimation into parameter estimation errors and stochastic irreducible errors to quantify uncertainty and improve forecast accuracy. Nevertheless, even with their method, the true covariance matrix cannot be perfectly determined, leaving room for further improvements in forccast accuracy.
arise between the estimated and true covariance matrices. Thus,
To address the uncertainty inherent in estimators, robust optimization has emerged as a powerful technique. This methodology is designed to yield solutions that remain effective even when the underlying data fluctuates within a defined uncertain range. Specifically, it involves establishing a range for uncertain parameters or data and then seeking an optimal solution that performs best under the worst-case scenario within that range. Since its inception by Ben- Tal & Nemirovski (1998), robust optimization has been extensively researched in both thcorctical and applicd domains (Bertsimas et al. 2011). Notably, models that incorporate covariance matrix uncertainty have been developed and applied to various problems, such as portfolio optimization
Lobo & Boyd 2000: Halldórsson & Tütüncū 2003 Tütūneñ & Koenig 2004).
2

---

Published in Transactions on Machine Learning Research (02/2026)
In this paper, propose novel method that frames hierarchical time-series forecasting robust op- timization problem. Our approach aims to minimize the average of forecast residuals over an observation period under the worst-case scenario within an uncertainty set for the covariance matrix of forecast errors. We demonstrate through duality that this robust optimization problem can be formulated as semidefinite optimization problem, which theoretically solvable efficiently. Furthermore, we present numerical experi- ments on five real-world datasets, showing that our proposed method achieves more accurate forecasts than existing hierarchical time-series forecasting techniques.
2
Hierarchical time-series forecasting
2.1
Notation
This section defines the notation used throughout the paper, which is consistent with prior studies on hierarchical time-series forecasting (Athanasopoulos et al. 2009: Hyndman et al. 2011 Wickramasuriya
et al. 2019 Panagiotelis et al., 2021 Hyndman & Athanasopoulos 2021).
A hicrarchical structure is defined by series of nested levels. Level is the fully aggregated series. Level consists of the series obtained by disaggregating the Level 0 series, and Level 2 contains the series that further disaggregate each Level series. This process continues until the bottom level, denoted as Level K, where its series can no longer be disaggregated.
Let 2 IR denote the observation of series X at time t. The label X is series of labels representing the indices of each level. For example, a series X that belongs to series i at Level 1, series j at Level 2, and series k at Level 3 can be denoted by ijk. The series at Level 0 is simply written as (()), without series name X. A key property of hierarchical data is that, at any given time point, the value of a series at specific level equals the sum of the values of the series nested directly below it:
(
Σ (E) U1
(t) Vig Vig
(1)
Σ
(4)
Viyk
To simplify the notation of the hierarchical structure, matrix and vector expression is often used. Let n be the total number of series and m be the number of bottom-level series, which satisfy the vector of all series observations at time t as y E R" and the vector of bottom-leve observations as B(O) E Rm. With the summing matrix E TR that dictates the way in which the bottom-level series aggregate, the hierarchical structure can written as:
We denote
y (t) = (t)
(1)
When Equation 1 holds for the values all series and the bottom-level series at each time t, said that the hierarchy satisfied.
Total
B
AA
AB
AC
BA
BB
Figure 1: An example of hierarchical structure
Figure D provides simple example of hierarchical structure. In this case, K = 2, 5, n 8, and the following aggregation relationship holds:
(4) (0) (1) (1) (6) (1) (0)
V UAA VAB + VAC> UB VBA UBB
3

---

Published in Transactions on Machine Learning Research (02/2026)
/0) yハ
(2)
(9
(9) (6) (2) (8)
+
Vis
UAA WAB + MAC UBA VBB
The corresponding matrix and vector representation is defined by:
y
y
yA
(8)
yB
(0)
MAA
VAB VAC UB VBB
(6)
(0)
609)
(8)
(0)
(0)
(0)
(0)
V
VAA
VAB
VAC
VBA
VBB
1
1
1
0
0
O
0
0 0
S
1
0
0
0
1
0
0
0
0
1
0
0
0
0
0
1
With these definitions, the hierarchical structure is fully captured by Equation
2.2 Reconciliation methods
Hierarchical time-series forecasting is process that adjusts, "reconciles," " forecast values to ensure they satisfy the hierarchical aggregation constraints. We assume that observations for each series available for an observation period t = I,...,T, and the objective is to forecast values for the forecast period 7 = 7 7 ++.....,T+77?
Forecasting, ignoring any aggregation constraints, is called a base forecast. Let g()) R" denote the vector of base forecasts for all series time 7. A forecast that satisfies the hierarchical aggregation constraints is called coherent forecast. Let g(r) E R² denote the vector of coherent forecasts time т. To transform base forecasts into coherent forecasts, hierarchical forecasting methods estimate reconciliation matrix P E Rmxn that maps the base forecasts to the boltom level. All hierarchical forecasting approaches can be expressed in the general form:
(P
where S is the summing matrix defined in Equation [
The bottom-up and top-down approaches are two traditional reconciliation methods. We again use the example from Figure [ to illustrate these concepts. In the bottom-up approach, coherent forecasts are derived by summing the base forecasts of the bottom-level series. This corresponds the reconciliation matrix:
0
0
0
0
0
P
0
0
0
0
0
0
0
0
0
0
1
0
0
0
1
Conversely, in the top-down approach, coherent forecasts are obtained by disaggregating the base forecast of the top-level series. If px the proportion that allocates the total values to each bottom-level series X, the reconciliation matrix becomes:
PAA
0
0
0
0
0
0
PAB
0
0
0
P
PAC
0
PBA
0
0
0
0
PBB
0
0
0
A common way to determine these proportions is based on the average historical proportions of the data.
4

---

Published in Transactions on Machine Learning Research (02/2026)
Because the bottom-up and top-down methods utilize base forecasts from only single level of aggregation, they rely on limited information. To overcome this, subsequent research has proposed methods that use base forecasts from all series to estimate a reconciliation matrix, thus producing more comprehensive, coherent forecasts.
Hyndman et al. (2011) proposed the generalized least. squares (GL.) reconciliation, regression-based ap- proach. In this method, the reconciliation matrix is chosen minimize errors between the base and coherent forecasts. Specifically, consider the regression model for base forecasts
g(T) = SB(r) + +
Let ((r) [ 1 ........ y) Rm be the expectation of base forecasts the bottom level, and let the error term e have zero mean with covariance matrix Σ = Var (e(*) were known, the minimumn variance unbiased estimator of ((r) would be the GLS estimator
y()) E Raxx, If (()
( ( S (ST)()))))) S g(ggg)
where E(r) is the generalized inverse of (r), Comparing this expression with Equation
2
gives
p EST S s(1!
In practice, however, the covariance matrix E(r) is unknown and cannot be estimated. It represents the covariance of the reconciliation errors at time T, but the errors are not observable until coherent forecasts have already been produced. [Hyndman et al. (2011) showed that, under the assumption that the errors themselves satisfy the aggregation constraints, (T can be replaced with the n-dimensional identity matrix I, This replacement is equivalent to computing the ordinary least squares (OLS) estimator instead of the GLS estimator.
Wickramasuriya al. (2019) proposed the minimum-trace (MinT) reconciliation, which determines the roconciliation matrix P by minimizing the total variance of crrors between the observed values and the coherent forecasts. Let
b+) Var - :(()) y (y
,....,,,))
E mRxx
denole the covariance matrix of those errors, and define the covariance matrix of the errors between observed values and base forecasts
(r) = Var (y(+) gi (1)
yr) E "x
Under the linear transformation in Equation the error covariance matrix of the coherent forecasts is expressed as (Tr) SPW(") Consequently, the optimal reconciliation matrix that minimizes the total error variance tr (V(r)) is given by
P ST (w(T)) $ -1 ST (w(rt)
(ST(W()))
-1 -s)
-1
A notable property of this reconciliation that the resulting coherent forecasts g()) are guaranteed least as accurate as the base forecasts (+) in terms of the weighted squared error using (W(+))
(3 )()) (w(+)) 3 (T) - g(d) < (y(1)) -g(t)) T (w(t)) (1(1) g(t)).
As noted by Panagiotelis et al. (2021 when the base forecast. (r) is unbiased, the trace minimization MinT is equivalent to minimizing the expected squared Euclidean norm of the errors between the observed values and the coherent forecasts:
min P
l (7) - Sp)(
(4)
5

---

Published in Transactions on Machine Learning Research (02/2026)
where ||-I] is the Euclidean norm. Furthermore, they provided a key insight that extends this Euclidean norm to more general weighted norm. Specifically, they showed that if we know the true covariance matrix W the solution in Equation 3 is invariant to the choice of the weight matrix M in the following generalized optimization problem:
(y SSP(n) M - (
min
(5)
where M is any positive definite matrix. This property implies that the specific weighting of errors does not alter the optimal reconciliation matrix P.
3 Robust reconciliation
The invariance property of Min' is a key advantage. At the same time, applying this to real-world data often involves practical considerations. For example, base forecasts may contain some empirical biases. Further- morc, because the true crror distribution is unknown, we typically replace the expectation Equation (5) with the average of residuals over specific observation period 1,..., T. In such cases, the choice of the weight matrix M can influence the performance of reconciliation.
In Min framework, it is natural to use (W(r))- as the weight matrix M to account for the different scales of errors across series. Although the true covariance matrix W(() is unknown, it is common to estimate it using the observations and base forecasts from the observation period = 1,.,... T. This estimate is usually assumed to constant for all forecast periods, so we consider using an estimated covariance matrix W as an alternative to ((r) for all T= + 1...., T
However, this single estimate W may not always be accurate due to limited data changes in the underlying error structure. To address this uncertainty, propose robust optimization framework. We treat W-1 an uncertain parameter and introduce uncertainly set. By formulating problem that mininizes the worst-case error over this uncertainty set, we aim to ensure stable forecast performance even when the estimate of W not perfect.
3.1 Formulation
As we focus on the empirical average over T time steps from the observation period, the minimization problem of weighted squared residuals is formulated as follows:
T
min P
Σ T
yot - - SPg
w- (y'( SPG
(6)
where ())) is the base forecast time in the observation period. To incorporate robustness against the uncertainty in estimating W, introduce uncertainty set M C R=*, Then, consider the following robust optimization problem:
Σ E - T M -
min max
MEM
where we omit. the coefficient: 1/T from the objective function, since it does not affect the optimal solution.
In robust optimization problems with uncertainty in the covariance matrix, box uncertainty set is often used (oboo & Boyd 2000 Halldórsson & Tütüncü, 2003). The box uncertainty set places upper and lower bounds on the inverse of covariance matrix, so that.
M M M 2
where M. M Rxx are the upper and lower bounds of M, respectively, and the inequalities M < M M represent element-wise inequalities, meaning My Mj s Mij for all i, 1, n. M 2 O denotes that
6

---

Published in Transactions on Machine Learning Research (02/2026)
M is positive semidefinite. Therefore, the proposed method determines the reconciliation matrix by solving the following robust optimization problem:
I
Σ
y1)) - -SPOl) M (v -v0
min max M
==1
(7)
s.t.
M M S M,
M 2 0.
3.2 Equivalent reformulation
We show that the min-ma max problem 7 is equivalently reformulated as a semidefinite optimization problem following the approach of Lobo & Boyd 2000) and Halldórsson & Tütüneū (2003). Let I and O be the identity matrix and zero matrix of appropriate dimensions, respectively. For two symmetric matrices and B of the same size, A B denotes the standard inner product of the two matrices, defined A B = tr (AB), which is the sum of the element-wist product of A and B.
Proposition 1. Assume that there cvists positive definite matric M' satisfying M < M M. Then, the problem con be equivalently reformulated as the following semidefinite optimization problem:
min X.XEE,P
M X MoX
X
X E
s.t.
=
E
I
(8)
E ( -SP .....,y")) "...,""SSPp77)
X X20,
where X and X are n-dimensional symmetric matric variables, and E is on X T-dimensional matric variable whose columns are the error vectors between the observed values and coherent forecasts over the observation period.
Proof. Fix P and consider the inner maximization of the problem
T
Σ ( - M -30'SSp'"
max M
(9)
s.t.
M S M M,
M 1=0 0.
Its dual problem is formulated as
min
M X M
X
8.6. X X Σ E(y0) SP (y0) spg 00,
(10)
X. X.
where X, X are dual variables. Let us define the matrix E = ( SPt ....,,, -7)-577 -P Then, by the Schur complement, the semidefinite constraint in the problem 10 follows:
is equivalently transformed as
X X- X
Σ y - SP 1-00SPp'') 20
X
X
EE
00
X
X
E I
O.
(11)
E

---

Published in Transactions on Machine Learning Research (02/2026)
Therefore, replace the semidefinite constraint of the problem the following equivalent dual problem:
10 with Equation
11
which yields
min X
M X
M
X
X
X X E E
(12)
s.t.
O,
1
X20.
Finally, we show the dual problem 12 has a strictly feasible solution, which implies that the strong duality between the problems (9] and 12 holds. Let J E RRxn be matrix whose clements are all onc. For the
fixed P, we consider the following solution:
=E= +EE++JJ+ +7J+1, X'=
where y > 0 is a sufficiently large constant such that EE +7J > O. X X' O. Focusing on the semidefinite constraint of the problem 12
Since we take >> sufficiently large,
we have
X
X
E
EE T
+1 I
E I
> O,
E
I
thus X X''
is strictly feasible. Since we assume that the primal problem (9) is strictly feasible, the strong
duality theorem holds and the optimal values of the two problems 19 and 12 are equal (Vandenberghe & Boyd 1996 As we took P arbitrarily, the above argument holds for all P. Therefore, the min-max problem is equivalently reformulated as the semidefinite optimization problem
Hence, the approach of the proposed method is to solve the semidefinite optimization problem and use the optimal solution for P the reconciliation matrix for hierarchical time-series forecasts.
3.3 Uncertainty set
In order to deal with the robust optimization problem described in the previous sections, it is necessary to determine the uncertainty set, ie., M and M in advance. We set the upper and lower bounds of the uncertainty set from the observation period data using bootstrap with reference to Bertsimas et al. (0018).
We summarize the method for determining the uncertainty set in Algorithm [ In the first step, we calculate parameter À for the shrinkage approach, similar to the existing hierarchical time-series forecasting described in Section 2.2 This shrinkage approach is applied to the inverse of covariance matrix estimated by unbiased variance from the observation period data. In the next step, estimate the inverse of covariance matrix of each sample obtained by sampling the data of the observation period. For sampling, we select the same number of time points as the observation period T with replacement. Then, using the shrinkage intensity parameter obtained in the first step, the shrinkage approach is applied to the sampled inverse covariance matrix. This sampling is repeated NB 1 times obtain NB inverse covariance matrices. In the last step, determine upper and lower bounds from each element. of the sampled inverse covariance matrices. Let 0 < 1, then set the width of the uncertainty set to be times the width of the maximum and minimum sample values of the inverse of covariance matrix.
4 Numerical experiments
To verify the effectiveness the proposed method, we compared its forecast accuracy with that of cxisting hicrarchical timo-serics forccasting methods across multiple real-world datascts. Code used in these experiments is available http::///gttuub..om//tttmaaaeeee
8

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm Designing uncertainly set
Input: n. T, (y (E JE.
Output: M, M
M estimate inverse of covariance matrix from (y Ji= and (o
À calculate shrinkage intensity parameter of M
for ===,, l,...., do
+ sample T data points from (0 and (0 with replacement (s) + estimate inverse of covariance matrix from (x ==== and (0))7
(0(")
Fp=
(a) t update (s by shrinkage estimator with X end for
for i = 1, , n, j j= 1....,77 n do
Mj f calculate 100 (1 a2/2 peercenteeee of {M(s)
Mij calculate 100 (1 peercentte of ()) )
end for
return M, M
4.1 Datasets
The numerical experiments used four datasets from prior studies. An overview each dataset provided below, including the number of hierarchical levels and series, as well the lengths of the observation and forecast periods. These four datasets allowed real-world conditions, inchuding hierarchies with different scales of levels and series.
assess the proposed method's performance under various
The first dataset is the Australian births data (Hyndman & Athanasopoulos 2021). It records the number of births in Australia every month from January 1975 to December 2022. As summarized in Table D the hierarchy consisted of single disaggregation level (K = 1) with nine bottom-level series (m = 9) and ten series in total (n 10). Tevel 1, the bottom level, disaggregated the national series into nine states and territorics: ACT, AUS, NSW, NT, QLD, SA, TAS, VIC, and WA. For this experiment, the first 516 months (January 1975 to December 2017) served the observation period, while the following 60 months (January 2018 to November 2022) comprised the forecast period.
Table 1: Hierarchy for Australian births data
Level
Number of series
Total series per level
lotal
state
The second dataset is the Australian tourism data AAthanasopoulos et al. 2009 which records the number of domestic travelers every quarter from Q1 1998 to Q4 2017. As Table shows, the hierarchy was structured with K = 2, m = 27, and n = 35. Level 1 disaggregated the national total into seven states: NSW, NT, QLD, SA, TAS, VIC, and WA. Note that this state grouping differed from the one used for the births datasct, as we followed prior work. Level 2, the bottom level, further broke down each state into finer geographic zones. Here, the observation period consisted the first 68 quarters (Q1 1998 Q4 2014) and the forecast period covered the subsequent 12 quarters (Q1 2015 Q4 2017).
Table 2: Hierarchy for Australian tourism data
Level
Number of series
Total series per level
total
state
zone
6-2-4-4-3-5-3
27

---

Published in Transactions on Machine Learning Research (02/2026)
The third dataset is the U.S. Walmart sales data (Mancuso al. 2021 J. It tracks weekly sales from Jamuary 3, 2011, to May 29, 2016. Table 3 shows the hierarchy as K 2, m 10, and n 14. Level 1 split the national total into three states: CA, TX, and WI. Level 2, the bottom level, further subdivided each state into its constituent stores: four in CA, three in TX, and three in WI. In this experiment, we used the first 261 weeks (January 3, 2011, to January 3, 2016) as the observation period and the next weeks (January 4, 2016, to May 29, 2016) for the forecast period.
Table 3: Hicrarchy for Walmart sales data
Level Number of series
Total series per level
total
state
store
4-3-3
10
The fourth dalaset is the Swiss electricity demand data (Nespoli al. 2020), recording electricity supply every ten minutes from Jamuary 13, 2018, to January 19, 2019. As outlined in Table 4 the hierarchy
K = 3, m = 24, and = 31. Level 1 disaggregated the grid into two synthetic meter aggregations: S1 and S2. Level 2 further subdivided each aggregation into two synthetic sub-aggregations: S11 and S12 for S1, and S21 and S22 for S2. Level 3, the bottom level, then separated each aub-aggregation into six individual meters. After converting the data to a daily unit, we used the first 353 days (January 13, 2018, to December 31, 2018) for the observation period and the subsequent 19 days (January 1, 2019, January 19, 2019) for the forecast period.
Table 4: Hierarchy for Swiss electricity demand data
Level
Number of series
Total series per level
grid
1
agg.
2
2
sub-agg. 2-2
4
meter
6-6-6-6
24
4.2 Experimental setup
This subsection details the benchmark methods, our proposed robust methods, and evaluation metrics.
Hierarchical forecasting first requires a set of base forecasts (Base), which do not account for the hierarchical structure. We generated these base forecasts using Prophet [Taylor & Letham 2018 an open-source time- series library from Meta, with its default settings. As a preliminary experiment, we also performed base forccasts using other timo-scrics forecasting methods such as ARIMA and LightGBM, but Prophot showed the best results. Our comparative methods included the bottom-up (BU), top-down (TD), GLS reconeiliation (GLS), and MinT reconciliation (MinT) approaches introduced in Section 2
For our proposed method (Robust), the optimization problem (MOSEK ApS. 2025). Moreover, we needed to set two parameters for the boolstrap to design the uncertainly set, the number of sampling Ng and the width of the uncertainty set a. In these experiments, Ng was fixed at 5000 and a was determined by validation from among the candidates specified in advance. For validation, we divided the observation period data into train and validation data in a ratio of 9:1, and used the a with the smallest. RMSE in the validation data. The candidates for were 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0.
8 was solved using the MOSEK solver
10

---

Published in Transactions on Machine Learning Research (02/2026)
As evaluation metrics, used the mean absolute error (MAE) and the root mean squared error (RMSE) for each series. For given series X, these are defined as:
T+T Σ
I+T Σ
(
MAE T
UX
RMSE
UX
T
T=741
Note that RMSE gives harsh evaluation when the prediction deviates significantly from the MAE.
4.3 Results and discussion
The experimental results for all datasets are summarized in Table To make the results casier understand, calculated the ratio of MAF and RMSE for each hierarchical time-series forecasting method when the MAE and RMSE of the base forecasts are set to 1, and then calculated the mcan and standard deviation for series included in the same hierarchical level. We defined these metrics as relative MAE and relative RMSE, respectively. The original MAE and RMSE of all series are listed in the appendix. In the tables, the format is "mean = standard deviation." Since the top level of hicrarchical structure (Level 0) includes only one series, the standard deviation for the series not shown. The underlined values indicate the best forecast. accuracy in terms of the mean for the corresponding hierarchical level.
For the Australian births dataset, the validation decide a parameter for bootstrap resulted in a 1.0. The proposed method achieved the highest accuracy in most cases. In the top level, our proposed method was the only approach that achieved higher accuracy than the base forecast. In the bottom level, considering the standard deviation, some series showed significant improvement in forecast accuracy with the proposed method.
Table 5: Forecast accuracy for Australian births data
(a) Relative MAE
Base
BU
TD
GLS
MinT
Robust
total
1.000
1.034
1.000
1.003
0.999
0.923
state
1.000
1.000 + 0.000
1.513 + 1.302
0.955 + 0.136
0.978 + 0.043
0.944 = 0.167
(b) Relative RMSE
Base
BU
TD
GLS
MinT
Robust.
total
1.000
1.014
1.000
1.001
1.000
0.968
state
1.000
1.000
0.000
1.297 + 0.831
0.970 + 0.090 0.989 + 0.026 0.973 = 0.103
In the Australian tourism dataset, the bootstrap parameter was determined via validation to be 0.5. Our proposed method yielded lower forecast accuracy than the existing GLS and MinT reconciliation methods, except for the top level. This underperformance seems to be due to a slight discrepancy between the estimated covariance matrix and the true one in the forecast period, which reduced the benefit of explicitly accounting for covariance uncertainty.
For the Walmart sales, the validation procodure for selecting the bootstrap parameter yielded a 0.9. Our method achieved the highest forecast accuracy for all levels, with substantial performance gains over existing methods.
On the Swiss electricity demand datascts, determined by validation, the bootstrap parameter was set to a = 0.7. Our proposed method achieved the highest forecast accuracy at the upper levels, but the accuracy deteriorated at the bottom level compared to the base forecast. In other words, the improvement at the upper levels came at the expense of accuracy the bottom level. As with the Australian tourism data, the discrepancy between the true and estimated covariance matrices may have been small, limiting the overall impact of robustification.
11

---

Published in Transactions on Machine Learning Research (02/2026)
Table 6: Forecast accuracy for Australian tourism data
(a) Relative MAE
Base
BU
TD
GLS
MinT
Robust
total
1.000
2.315
1.000
1,125
0.663
0.518
state
1.000
1.277 + 0.247
0.968
0.400
0.635 + 0.181
0.573 + 0.135
0.793 0.460
zone
1.000
1.000 = 0.000
0.925 0.380
0.623 + 0.166
0.568 0.158
0.843 0.387
(b) Relative RMSE
Base
BU
TD
GLS
MinT
Robust
total
1.000
2.132
1.000
1.105
0.698
0.519
state
1.000
1.227
0.193
0.989 + 0.388
0.680 + 0.185
0.617 + 0.135
0.810 + 0.428
zone
1.000
1.000 0.000
0.959 + 0.367
0.657 0.154
0.606 + 0.141
0.844 + 0.339
Table 7: Forecast accuracy for Walmart sales data.
(a) Relative MAE
Base
BU
TD
GLS
MinT
Robust
total
1.000
1,135
1.000
1.011
0.966
0.810
state
1.000
1.123 + 0.117
1.263
0.986
0.999 + 0.029
0.956
0.049
0.801 - 0.036
store
1.000
1.000 + 0.000
1.618
1.001
0.882 + 0.069
0.868 0.076
0.741 0.117
(b) Relative RMSE
Base
BC
TD
GLS
MinT
Robust
total
1.000
1.129
1.000
1.011
0.967
0.817
state
1.000
1.118
+ 0.116
1.241 + 0.948
0.998 + 0.030
0.955 + 0.050
0.807 + 0.038
store
1.000
1.000 0.000
1.580 + 0.945
0.885 0.067
0.872 + 0.070
0.751 + 0.111
Table 8: Forecast accuracy for Swiss electricity demand data
(a) Relative MAE
Base
BU
TD
GLS
MinT
Robust
grid
1.000
1.041
1.000
1.022
1.110
0.888
agg.
1.000
1.045
0.034
0.983 = 0.422
0.983 = 0.020
= 1995=00.011
0.808 + 0,277
sub-agg.
1.000
1.072
0.071
1.390 + 0.549
0.999
0.010
1.014 0.099
0.867 + 0.232
meter
1.000
1.000 + 0.000
1.810 0.921
0.987 + 0.102
1.023 + 0.177
1.036 + 0.275
(b) Relative RMSE
Base
BU
TD
GLS
Min'T
Robust
grid
1.000
1.022
1.000
1.012
1.050
0.946
agg.
1.000
1.013 + 0.045
0.970 = 0.213
0.986 = 0.008
1.005 = 0.023
0.895 + 0.171
sub-agg.
1.000
1.032 + 0.058
1.170 = 0.359
0.995 = 0.003
1.025
0.071
0.929 + 0.158
meter
1.000
1.000 + 0.000
1.717 0.781
0.981 0.086
1.021 + 0.141
1.013 + 0.227
12

---

Published in Transactions on Machine Learning Research (02/2026)
The numerical experiments on these four datasets demonstrated the effectiveness of our proposed hierarchical time-series forecasting method. Across most target series and hierarchical levels, the proposed approach achieved higher forecast accuracy than existing techniques. However, in a few cases, the proposed method's prediction performance was slightly inferior to that of the GLS and MinT reconciliation methods.
As discussed in the Australian tourism data results, when the discrepancy between the estimated and true covariance matrices is small, incorporating the uncertainty may lead to overly conservative forecasts, resulting in lower predictive performance. This is further supported by the parameter selection during validation: a larger was chosen for the Australian births and Walmart sales datasets, where method was effective. Conversely, smaller was selected for the Australian tourism and Swiss electricity demand datasets, where the method was less effective. This suggests that for these latter datasets, the estimated covariance matrix closely approximated the true covariance matrix, reducing the need to account for its uncertainty in the reconciliation process.
5 Conclusion
In this paper, we proposed a robust hierarchical time-series forecasting method. This approach introduces an uncertainty for the inverse of covariance matrix of forecast errors and minimizes the forecast error between the observation values and the coherent forecasts. Through numerical experiments, we demonstrated that the proposed method often provides more accurate forecasts than existing hierarchical timo-series forecasting vary across different datasets.
methods, although its performance
The limitations of our method are twofold: its accuracy may not always surpass existing methods, and its scalability is limited due to the computational demands of the optimization problem. The first issuc is as stated in the discussion of the experimental results. It arises when the discrepancy between the true and estimated covariance matrices small, in which case the robust approach offers little advantage. The second limitation concerns the computational time required to obtain reconeiliation matrix, as our method relies solving a semidefinite optimization problem. The size of the optimization problem depends on the total number of series and the length of observation periods. While a solution was achievable in tens of seconds for the dataset sizes used in our experiments, this approach would not be practical for very large-scale predictions.
Acknowledgments
This study was conducted part of the Data Analysis Competition hosted by the Joint Association Study Group of Management Science. The authors would like to thank the organizers and Rescarch and Innovation Co., Ltd.
References
George Athanasopoulos, Roman A Ahmed, and Rob J Hyndman. Hierarchical forecasts for australian domestic tourism. International Journal of Forecusting, 25(1):146-IM) 1666,666666 2009.
Yossi Aviv. A time-series framework for supply-chain inventory management. Operations Research, 51 (2): 210 227, 2003.
Aharon Ben-Tal and Arkadi Nemirovski. Robust convex optimization. Mathematics of Operations Research, 7344)77698805 1998.
Dimitris Bertsimas, David B Brown, and Constantine Caramanis. Theory and applications of robust opti- mization. SIAM Review, 5364:4664- 2011.
Dimitris Bertsimas, Vishal Gupta, and Nathan Kallus. Data-driven robust optimization. Mathematical Programming, 167:235-22 292, 2018.
John B Edwards and Guy Orcutt. Should aggregation prior to estimation be the rule? The Review of Economics and Statistics, pp. 409-420, 420, 1969.
13

---

Published in Transactions on Machine Learning Research (02/2026)
Yehuda Grunfeld and Zvi Griliches. Is aggregation necessarily bad? The Review of Economics and Statistics, pp. 1-13, 1960.
Bjarni V Halldórsson and Reha Tülüneü. An interior-point method for a class of saddle-point problems. Journal of Optimization Theory and Applications, 116:559 590, 2003.
Luis Hernandez, Carlos Baladron, Javier M Aguiar, Belen Carro, Antomio J Sanchez-Esguevilne Esguevillas, Jaime Lloret, and Joaquim Massana. A survey on electric power demand forecasting: Future trends in smart grids, microgrids and smart buildings. IEEE Communicatione Surveys & Tulorials, 16(3):1460 1495, 2014.
Rob J Hyndman and George Athanasopoulos. Porecosting: Principles and practice. OTexts, Melbourne, Australia, 3rd edition, 2021. URL
Accessed on 2025-05-06.
Rob J Hyndman, Roman Ahmed, George Athanasopoulos, and Han Lin Shang. Optimal combination forecasts for hierarchical time series. Computational Statistics & Data Analysis, 55(9): :2799 2589, 2011.
Bjoern Krollner, Bruce Vanstone, and Gavin Finnie. Financial time series forecasting with machine learning techniques: A survey. In European Symposium on Artificial Newral Networks: Computational Intelligence and Machine Learning, pp. 25 30, 2010.
Zhenyu Liu, Zhengtong Zhu, Jing Gao, and Cheng Xu. Forecast methods for time series data: A survey, IEEE Access, 9-91896-91912, 2021.
Miguel Sousa Lobo and Stephen Boyd. The worst-case risk of a portfolio. Technical report, Stanford University, 2000. URL https://web. stanford. http:///too///ppppp///////////// 2025-05-06.
bndiiisss___d___ee bnd. pdf Accessed on
Ganapathy Mahalakshmi, S Sridevi, and Shyamsundar Rajaram. A survey on forecasting of time series data. 2016 International Conference on Computing Technologies and Intelligent Data Engineering (ICCTIDE'66) pp. 8-8. TEEE, 2016.
Paolo Mancuso, Veronica Piccialli, and Antonio M Sudoso. A machine learning approach for forecasting hierarchical time series. Expert Systems with Applications, 182, 2021.
Jan Kloppenborg Moller, Peter Nystrup, Poul G Hjorth, and Henrik Madsen. Optimal forecast reconciliation with uncertainty quantification. arXiv preprint ar Xrx: 06480, 2024.
MOSEK ApS. The MOSEK Python Fusion APT manual. Version 11.0., 2025. URL https:://ocss.ooot .oose.kk......k........... Accessed on 2025-05-06.
Lorenzo Nespoli, Vasco Medici, Kristijan Lopatichki, and Fabrizio Sossan. Hierarchical demand forecasting benchmark for the distribution grid. Electric Power Systems Research, 189:106755, 2020.
Guy H Orcutt, Harold W Watts, and John Edwards. Data aggregation and information loss. The American Economic Review, 58(4):773 787, 1968.
Anastasios Panagiotelis, George Athanasopoulos, Puwasala Gamakumara, and Rob J Hyndman. Forecast reconciliation: A geometric view with new insights on bias correction. International Journal of Forecasting, 37(1))343 359, 2021.
Patrícia Ramos, Nicolau Santos, and Rui Rebelo. Performance of state space and arima models for consumer retail sales forecasting. Robotics and Computer-integnated Manafacturingg 34:151-163, 2015.
Tomokaze Shiratori, Ken Kobayashi, and Yuichi Takano. Prediction of hicrarchical time scries using struc- tured regularization and its application to artificial neural networks. PLOS ONE, 15(1)))e0242099, 2020.
L Suganthi and Anand A Samuel. Energy models for demand forecasting a review. Renewable and Sus- tainable Energy Reviews, 16(2):1223-122) -2440 2012.
Sean J Taylor and Benjamin Letham. Forecasting at scale. The American Statistician, 72(1):37 -4,7-45,55 2018.
14

---

Published in Transactions on Machine Learning Research (02/2026)
Reha H Tütüncü and Mark Koenig. Robust asset allocation. Annals of Operations Research, 132:157-187, 2004.
Lieven Vandenberghe and Stephen Boyd. Semidefinite programming. SIAM review, 38(1):49-99) 95, 1996.
Qingsong Wen, Tian Zhou, Chaoli Zhang, Weiqi Chen, Ziqing Ma, Junchi Yan, and Liang Sun. Transformers in time series: A survey. arXiv preprint ar 221:222000722 2022.
Shanika L Wickramasuriya, George Athanasopoulos, and Rob J Hyndman. Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. Journal of the American Statistical Association, :84 819, 2019.
A Full results of numerical experiments
This section reports the original MAE and RMSE for all series from the numerical experiments that could not be included in Section 4.3] Table and Table I0 present the MAE and RMSE for the Australian births data, respectively. Similarly, Table [I and Table 12 report the forecast accuracy for the Australian tourism data, Table 13 and Table 14 for the Walmart sales data, and Table 15 and Table I6 for the Swiss electricity demand data. Within each lable, series separated by horizontal lines correspond to the same hierarchical level.
Table 9: MAE for Australian births data
Base
BU
TD
GLS
MinT
Robust.
total
3697.93
3822.66
3697.93
3710.04
3695.48
3413.01
ACT
79.34
79.34
23.95
65.71
77.18
71.65
AUS
1863.94
1863.94
1848.74
1851.83
1819.48
1706.96
NSW
433.36
433.36
934.11
423.74
410.68
372.41
NT
41.86
41.86
43.48
30.60
39.88
36.31
QLD
222.39
222.39
231.03
217.53
217.92
202.75
SA
107.33
107.33
328.02
97.16
100.02
83.49
TAS
40.31
40.31
163.39
49.53
43.70
55.08
VIC
711.71
711.71
430.98
699.45
700.14
675.65
WA
479.92
479.92
180.03
466.29
468.86
451.56
Table 10:
RMSE for Australian births data.
Base
BU
TD
GLS
MinT
Robust
total
6419.88
6507.75
6419.88
6428.59
6418.09
6213.25
ACT
93.99
93.99
44.91
81.93
92.04
87.08
AUS
3220.69
3220.69
3209.85
3211.97
3188.33
3106.97
NSW
782.22
782.22
1165.31
774.65
764.60
736.15
NT
55.47
55.47
57.62
45.46
53.67
50.50
QLD
452.41
452.41
462.21
447.36
447.76
433.83
SA
166.40
166.40
354.25
157.74
160.18
146.77
TAS
54.89
54.89
169.39
62.80
57.77
67.82
VIC
1216.24
1216.24
1004.19
1207.09
1207.60
1190.46
WA
626.98
626.98
383.23
615.56
617.71
603.22
15

---

Published in Transactions on Machine Learning Research (02/2026)
Table 11: MAE for Australian tourism data
Base
BU
TD
GL.S
MinT
Robust
total
1507.73
3490.86
1507.73
1696.30
999.98
780.57
NSW
706.38
1056.55
336.42
594.78
445.01
678.46
NT
122.88
122.99
151.28
57.50
82.99
71.12
QLD
565.88
626.18
503.91
427.08
206.32
219.97
SA
213.29
278.99
90.94
82.64
83.73
78.47
TAS
132.18
132.50
158.77
67.01
84.52
62.70
VIC
509.36
828.63
532.14
416.27
347.01
640.93
WA
331.08
462.09
498.96
222.41
207.11
503.88
NSW ACT
131.54
131.54
68.63
63.37
68.05
102.13
NSW Nth
139.01
139.01
92.79
75.93
82.49
57.30
NSW Sth
131.20
131.20
115.03
77.85
94.87
165.99
NSW Metro
313.69
313.69
161.47
241.31
142.86
323.77
NSW Nth Coast
237.56
237.56
119.50
160.60
130.42
190.16
NSW_St Coast
118.41
118.41
136.72
53.19
44.83
110.12
NT Central
72.85
72.85
84.78
37.56
51.82
32.97
NT_Nth _ Coast
50.14
50.14
67.35
30.65
31.72
76.22
QLD Metro
360.78
360.78
209.51
318.68
197.90
173.06
QLD Central Coast
53.20
53.20
69.45
33.78
38.48
60.19
QLD Inland
150.29
150.29
172.18
107.54
64.82
84.05
QLD Nth Coast
80.10
80.10
140.75
58.20
63.47
128.42
SA Metro
127.55
127.55
34.11
76.69
36.13
38.15
SA Inland
78.95
78.95
30.53
34.33
32.78
38.62
SA. West West__Coast
27.14
27.14
24.60
31.01
16.42
22.17
SA Sth Coast
52.22
52.22
70.85
37.05
40.02
57.76
TAS _Nth East
50.76
50.76
61.78
29.15
32.34
39.28
TAS Sth
53.34
53.34
70.71
39.37
41.03
46.34
TAS Nth West
39.50
39.50
33.94
23.07
25.23
21.32
VIC Nth West
100.44
100.44
64.15
53.47
53.19
59.12
VIC Nth East
209.01
209.01
86.61
127.61
98.84
87.13
VIC Metro
335.23
335.23
339.66
250.59
176.29
335.11
VIC East Coast
113.59
113.59
126.28
73.29
79.66
154.80
VIC West. Coast.
85.57
85.57
93.60
39.71
43.26
59.43
WA West Coast
237.14
237.14
326.91
171.87
199.24
392.71
WA Sth
96.10
96.10
56.24
30.70
30.76
70.90
WA Nth
130.13
130.13
115.82
50.24
35.39
57.09
16

---

Published in Transactions on Machine Learning Research (02/2026)
Table 12: RMSE for Australian tourism data
Base
BU
TD
GL.
MinT
Robust
total
1726.07
3680.00
1726.07
1906.48
1204.56
895.21
NSW
814.62
1153.92
390.45
709.89
556.17
767.23
NT
135.97
136.13
177.71
62.33
101.17
90.27
QLD
633.93
691.06
581.55
503.91
264.39
267.44
SA
226.56
291.26
105.01
98.52
99.75
94.40
TAS
148.29
148.68
185.51
86.92
96.55
78.14
VIC
642.69
938.29
697.81
553.19
475.28
725.78
WA
362.15
481.64
513.88
272.77
233.58
567.47
NSW ACT
146.72
146.72
86.05
79.23
83.41
117.14
NSW Nth
158.32
158.32
110.01
93.05
99.64
73.21
NSW Sth
145.49
145.49
128.61
90.52
104.94
186.20
NSW Metro
370.82
370.82
192.91
300.52
178.58
365.35
NSW Nth Coast
265.32
265.32
145.86
193.64
162.47
222.41
NSW BSW_Sth__ Stt_Coastt
131.95
131.95
165.04
64.79
54.05
122.21
NT Central
87.08
87.08
102.50
48.88
69.71
41.73
_T_Nth ____tt__Co_t Coast
58.57
58.57
81.15
36.06
43.98
81.02
QLD Metro
432.75
432.75
278.40
389.27
243.99
207.95
QLD Central_Coass
62.53
62.53
78.83
42.89
48.47
64.68
QLD Inland
174.46
174.46
194.94
131.45
78.24
109.60
QLD Nth Coast
101.77
101.77
183.26
73.57
70.65
145.45
SA Metro
134.08
134.08
41.14
84.94
45.04
43.86
SA Inland
87.19
87.19
36.99
43.59
42.36
46.13
SA West. Coast
31.63
31.63
33.36
35.80
21.76
31.26
SA Sth Coast
66.34
66.34
82.00
42.47
46.66
69.52
TAS Nth East
58.96
58.96
70.44
35.35
39.60
47.58
TAS Sth
61.49
61.49
84.42
47.51
48.14
55.55
TAS Nth West
44.88
44.88
42.50
30.26
33.01
25.17
VIC Nth West
124.89
124.89
86.11
76.51
74.90
81.43
VIC Nth East
232.47
232.47
108.39
153.51
116.29
106.73
VIC Metro
385.23
385.23
397.61
307.86
238.76
370.23
VIC East Coast
142.78
142.78
164.21
87.36
97.22
177.92
VIC West Coast
101.05
101.05
122.62
43.97
50.62
68.44
WA West Coast
275.32
275.32
364.37
214.80
225.39
446.61
WA Sth
101.13
101.13
63.93
41.77
42.26
79.92
WA Nth
135.48
135.48
134.23
63.10
48.77
66.16
17

---

Published in Transactions on Machine Learning Research (02/2026)
Table 13: MAE for Walmart sales data
Base
BU
TD
GLS
MinT
Robust
total
2343.98
2659.41
2343.98
2370.30
2264.97
1898.06
CA
1237.39
1397.73
473.01
1248.40
1195.54
993.95
TX
586.12 12
587.64 64
1364.16
566.09
529.02
447.96
WT
551.73
682.22
595.33
563.94
551.06
460.93
CA
217.28
217.28
509.66
182.14
196.40
169.73
CA
707.61
707.61
626.00
670.27
556.24
426.06
CA
205.31
205.31
371.79
167.98
193.55
174.84
CA
268.06
268.06
153.50
230.73
250.44
224.32
TX
139.57
139.57
496.24
132.39
122.58
95.30
TX 2
210.50
210.50
505.23
203.52
190.62
161.45
TX 3
239.12
239.12
362.70
232.15
216.86
193.64
WI
266.21
266.21
60.84
226.78
197.83
129.05
2
178.76
178.76
171.12
139.94
135.31
138.43
WT
238.70
238.70
456.46
199.22
220.13
195.16
Table 14: RMSE for Walmart sales data
Base
BU
TD
GLS
MinT
Robust
total
2479.68
2799.50
2479.68
2507.51
2398.58
2026.11
CA
1315.76
1466.86
502.11
1323.84
1267.95
1076.92
TX
613.93
615.60
1386.53
592.66
553.26
469.26
WI
586.42
724.49
635.82
599.93
586.67
491.14
CA 1
228.58
228.58
516.72
191.94
206.43
178.76
CA
761.58
761.58
630.90
727.71
623.03
518.46
CA
212.33
212.33
377.04
174.67
200.11
180.73
CA
278.59
278.59
163.45
241.39
260.40
233.43
TX 1
147.41
147.41
503.43
139.73
129.54
101.01
TX 2
219.32
219.32
511.96
211.67
198.63
168.66
T
250.99
250.99
371.73
243.43
227.41
202.47
WI
283.85
283.85
90.82
241.77
211.29
136.84
WI 191.39
191.39
183.70
151.72
147.07
151.58
WI
252.62
252.62
464.93
211.08
233.04
207.25
18

---

Published in Transactions on Machine Learning Research (02/2026)
Table 15: MAE for Swiss electricity demand data
Base
BU
TD
GLS
Min
Robust
grid
39.471
41.091
39.471
40.343
43.814
35.068
S1
18.902
19.298
24.210
18.845
19.039
18.977
S2
32.668
34.938
22.365
31 .6.633
32.098
19.997
S11
7.783
8.606
16.944
7.880
7.647
7.629
S12
11.358
11.356
14.279
11.362
11.437
11.348
S21
24.839
25.529
22.372
24.545
22.737
12.900
S22
8.388
9.680
10.263
8.347
9.649
8.115
S11 1
1.087
1.087
1.864
0.783
0.478
0.489
S11 2
2.265
2.265
4.394
2.627
2.666
2.000
S11 3
1.710
1.710
2.321
2.078
2.142
1.877
S11 4
3.008
3.008
3.335
2.968
2.909
2.994
S11 5
6.053
6.053
15.819
5.771
5.806
5.135
$11 6
2.619
2.619
3.602
2.724
2.753
2.616
S12
4.558
4.558
8.822
4.563
4.760
5.088
S12 2
1.758
1.758
1.667
1.756
2.023
2.144
S12 3
4.519
4.519
7.814
4.515
5.325
4.572
S12 4
4.263
4.263
6.355
4.268
4.233
3.121
S12_5
4.499
4.499
4.249
4.504
4.169
5.501
S12 6
2.505
2.505
2.161
2.502
2.884
2.610
$21 1
1.359
1.359
2.935
1.470
1.174
1.490
S21 2
1.124
1.124
1.135
1.180
L.106
1.297
$21 3
0.872
0.872
2.793
0.848
0.841
1.385
821 4
2.751
2.751
10.583
2.812
2.745
3.352
$21 5
6.562
6.562
4.230
6.398
7.161
5.460
$21 6
17.541
17.541
18.843
17.377
13.768
9.224
S22 1
2.173
2.173
5.764
1.866
2.493
2.296
S22 2
1.928
1.928
7.793
1.687
2.118
2.091
S22 3
3.880
3.880
8.993
3.711
3.702
4.317
S22
1.100
1.100
1.983
1.198
1.093
1.101
$22 5
1.225
1.225
1.180
1.088
1.266
1.066
$22 6
1.168
1.168
1.998
0.994
1.581
2.003
19

---

Published in Transactions on Machine Learning Research (02/2026)
Table 16: RMSE for Swiss electricity demand data
Base
BU
TD
GLS
Min
Robust
grid
64.301
65.685
64.301
65.102
67.498
60.821
S1
30.940
30.341
34.684
30.683
31.597
31.432
S2
39.665
41.439
32.503
38.863
39.198
30.729
S11
12.023
11.953
20.456
11.951
12.321
12.442
S12
19.037
18.966
17.871
18.976
19.450
19.160
S21
28.347
28.919
27.276
28.105
26.616
19.665
S22
11.817
13.200
12.733
11.781
13.164
11.596
S11 1
1.166
1.166
2.088
0.877
0.649
0.656
S11 2
2.800
2.800
5.079
3.117
3.152
2.526
S11 3
1.968
1.968
2.775
2.304
2.365
2.111
S11 4
3.466
3.466
4.334
3.409
3.596
3.862
S11 5
6.399
6.399
16.354
6.105
6.142
5.449
$11 6
4.682
4.682
5.583
4.823
4.850
4.666
S12
6.568
6.568
12.222
6.571
6.566
6.621
S12 2
2.350
2.350
2.105
2.348
2.586
2.688
S12 3
5.037
5.037
8.329
5.034
5.782
5.077
S12 4
4.726
4.726
7.384
4.730
4.707
4.183
S12_5
5.247
5.247
5.082
5.252
5.062
6.016
S12 6
2.744
2.744
2.989
2.741
3.099
2.844
$21 1
1.739
1.739
3.783
1.800
1.655
1.811
S21 2
1.566
1.566
1.614
1.592
1.560
1.648
$21 3
1.022
1.022
3.715
0.996
0.993
1.625
821 4
4.406
4.406
11.238
4.360
4.415
4.228
$21 5
7.181
7.181
5.095
7.034
7.725
6,216
$21 6
18.548
18.548
19.619
18.400
15.199
11.585
S22 1
2.559
2.559
6.203
2.256
2.858
2.672
S22 2
2.339
2.339
8.260
2.046
2.590
2.557
S22 3
6.345
6.345
10.247
6.159
5.731
5.548
S22
1.326
1.326
2.465
1.443
1.310
1.320
$22 5
1.717
1.717
1.497
1.522
1.758
1,484
$22 6
1.424
1.424
2.417
1.226
1.838
2.248
20


Reviews:
[{"review_id": "lEpE8iYNUp", "reviewer": "Reviewer_z2xd", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper is proposed a way to improve upon the reconciliation method for hierarchical time series problems. Its main starting point is the minT method, which it improves upon. Whereas in minT, the co-variance used for least squares estimation is a single step process consisting of all the samples, in the current method, they take bootstrap samples from the dataset and come up with min and max (or rather, low and high percentile) values of these. They then solve a min-max problem (min P max W Error ), so that it can be performant even under the worst case scenarios of W - this is different from minT which does it over a single pass. \n\nThe minmax problem is solved using positive semi-definite programming. I am not very familiar with the methods used so I cannot comment on the SDP approach with confidence. \n\nThey carry out evaluations over three different hierarchical time series datasets (Australian tourism, Walmart sales and Swiss electricity demand), and show comparisons with other methods (topdown, bottom up, ordinary least squares, minT). The numbers show that they with the robust method, they are able to obtain lower MAE/RMSE compared to the other methods. \n\nI would like to emphasize that I am not an expert in this area and my confidence estimate of the review is therefore is low.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The motivation is very clear in that the method overcomes limitations of minT which can be sensitive to the choice of W. Results are better than minT.", "audience_yn": "Yes", "audience_explanation": "This is an important addition to the time series reconciliation problem, overcoming some of the limitations of the well known minT method.", "requested_changes": "Could the authors provide some insight into why the method underperforms in the tourism dataset?", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "utQCBse4TZ", "reviewer": "Reviewer_4M2K", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The main technical problem in hierarchical time series forecasting is the estimation of the reconciliation matrix that encodes the hierarchical aggregation constrains. In state-of-the-art techniques this reconciliation matrix is obtained by estimating the covariance matrix of the errors between observation values and base forecasts. In this paper they propose a novel robust method that estimates the reconciliation matrix as the min-max solution over an uncertainty set of possible covariance matrices. The resulting min-max optimization problem is solved by formulating it as a semidefinite program. The method is tested in several dataset and results in improved MAE and RMSE metrics for in most of them when compared to state-of-the-art hierarchical time series forecasting methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The problem statement, modelling steps and the experimental setup are explained with sufficient detail and clarity.", "audience_yn": "Yes", "audience_explanation": "The paper may be interesting for both the time series forecasting and robust optimization communities.", "requested_changes": "Your method is basically is a tool to convert a base forecast into a hierarchical forecast, and in theory, it is agnostic to the base forecasting technique used. In the experiments you mention the use of the Prophet library for base forecasting. Have you compared your method to other reconciliation methods using different base forecasts? In principle, it could happen that different base forecasting techniques interacted in unexpected ways with the different reconciliation methods you analyze. From what you mention in subsection 4.2, it wouldn't be hard to add some extra experiments to verify that your method is consistently superior over different base forecasting techniques. If there are theoretical reasons to expect that, they also deserve some discussion.", "broader_impact_concerns": "No", "additional_comments": "N/A", "extra_scores": {}, "extra_text": {}}, {"review_id": "2QxTBGFjqZ", "reviewer": "Reviewer_DVTz", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors suggest using a robust optimization procedure for the reconciliation of hierarchical time-series forecasts. The authors start with the MinT method proposed by Wickramasuriya et al. (2019). They modify an optimization problem to a robust optimization problem to account for errors in covariance matrix estimation. The efficiency of the proposed method is demonstrated in experiments with time-series data.\n\n**Strengths**\n1. The motivation for the underlying problem is clear. Time-series forecasting is an important problem for practical applications.\n2. The writing is clear.\n3. The proposed approach seems novel.\n\n**Weaknesses**\n1. The motivation for the exact research question of the paper is vague (see below).\n2. The paragraph about the limitations of the method at the end of Section 4.3 does not directly lead to a practical rule to choose between MinT and the new method.", "claims_and_evidence_yn": "No", "claims_explanation": "1. I am confused about the optimization problem. Panagiotelis et al. (2021) indeed showed that the MinT method is equivalent to the optimization of Equation (3). However, the canonical weight matrix in this problem is $W = (W^{(\\tau)})^{-1}$, not $W = W^{(\\tau)}$.\n2. Similarly, Theorem 3.3 in Panagiotelis et al. (2021) shows that the MinT solution is optimal for any matrix. It suggests that the exact choice of matrix $W$ might not be so crucial for the optimization. Thus, the rationale for applying robust methods to the matrix $W$ is unclear.\n3. The current theoretical results in Section 3 do not demonstrate the advantage of the proposed method in terms of expected error.", "audience_yn": "Yes", "audience_explanation": "I think the paper might be of interest to communities that perform time-series forecasting.", "requested_changes": "1. Please explain the motivation behind the approach.\n2. Explain why the expected error is supposed to decrease.\n3. (minor) Please expand on the limitations and provide a practical guide for using your method.", "broader_impact_concerns": "I am not concerned about the broader ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)