Time series can often be naturally disaggregated by various attributes of interest. For example, the total number of bicycles sold by a cycling manufacturer can be disaggregated by product type such as road bikes, mountain bikes and hybrids. Each of these can be disaggregated into finer categories. For example hybrid bikes can be divided into city, commuting, comfort, and trekking bikes; and so on. These categories are nested within the larger group categories, and so the collection of time series follows a hierarchical aggregation structure. Therefore we refer to these as “hierarchical time series”.
Hierarchical time series often arise due to geographic divisions. For example, the total bicycle sales can be disaggregated by country, then within each country by state, within each state by region, and so on down to the outlet level.
Alternative aggregation structures arise when attributes of interest are crossed rather than nested. For example, the bicycle manufacturer may be interested in attributes such as frame size, gender, price range, etc. Such attributes do not naturally disaggregate in a unique hierarchical manner as the attributes are not nested. We refer to the resulting time series of crossed attributes as “grouped time series”.
More complex structures arise when attributes of interest are both nested and crossed. For example, it would be natural for the bicycle manufacturer to be interested in sales by product type and also by geographic division. Then both the product groupings and the geographic hierarchy are mixed together. We introduce alternative aggregation structures in Section 11.1.
Forecasts are often required for all disaggregate and aggregate series, and it is natural to want the forecasts to add up in the same way as the data. For example, forecasts of regional sales should add up to forecasts of state sales, which should in turn add up to give a forecast for national sales.
In this chapter we discuss forecasting large collections of time series that aggregate in some way. The challenge is that we require forecasts that are coherent across the entire aggregation structure. That is, we require forecasts to add up in a manner that is consistent with the aggregation structure of the hierarchy or group that defines the collection of time series. In this chapter we will rely on StatsForecast (Garza et al. (2022)) and HierarchicalForecast(Olivares et al. (2024)).
11.1 Hierarchical and grouped time series
Hierarchical time series
Figure 11.1 shows a simple hierarchical structure. At the top of the hierarchy is the “Total”, the most aggregate level of the data. The tth observation of the Total series is denoted by y_t for t=1,\dots,T. The Total is disaggregated into two series, which in turn are divided into three and two series respectively at the bottom level of the hierarchy. Below the top level, we use y_{j,t} to denote the tth observation of the series corresponding to node j. For example, y_{\text{A},t} denotes the tth observation of the series corresponding to node A, y_{\text{AB},t} denotes the tth observation of the series corresponding to node AB, and so on.
Figure 11.1: A two level hierarchical tree diagram.
In this small example, the total number of series in the hierarchy is n=1+2+5=8, while the number of series at the bottom level is m=5. Note that n>m in all hierarchies.
For any time t, the observations at the bottom level of the hierarchy will sum to the observations of the series above. For example,
y_{t}=y_{\text{AA},t}+y_{\text{AB},t}+y_{\text{AC},t}+y_{\text{BA},t}+y_{\text{BB},t},
\tag{11.1}
y_{\text{A},t}=y_{\text{AA},t}+y_{\text{AB},t}+y_{\text{AC},t}\qquad \text{and} \qquad y_{\text{B},t}=y_{\text{BA},t}+y_{\text{BB},t}.
\tag{11.2} Substituting Equation 11.2 into Equation 11.1, we also get y_{t}=y_{\text{A},t}+y_{\text{B},t}.
Example: Australian tourism hierarchy
Australia is divided into six states and two territories, with each one having its own government and some economic and administrative autonomy. For simplicity, we refer to both states and territories as “states”. Each of these states can be further subdivided into regions as shown in Figure 11.2 and Table 11.1. In total there are 76 such regions. Business planners and tourism authorities are interested in forecasts for the whole of Australia, for each of the states and territories, and also for the regions.
Figure 11.2: Australian states and tourism regions.
Table 11.1: Australian tourism regions.
State
Region
Australian Capital Territory
Canberra
New South Wales
Blue Mountains, Capital Country, Central Coast, Central NSW, Hunter, New England North West, North Coast NSW, Outback NSW, Riverina, Snowy Mountains, South Coast, Sydney, The Murray.
Northern Territory
Alice Springs, Barkly, Darwin, Kakadu Arnhem, Katherine Daly, Lasseter, MacDonnell.
Queensland
Brisbane, Bundaberg, Central Queensland, Darling Downs, Fraser Coast, Gold Coast, Mackay, Northern Outback, Sunshine Coast, Tropical North Queensland, Whitsundays.
South Australia
Adelaide, Adelaide Hills, Barossa, Clare Valley, Eyre Peninsula, Fleurieu Peninsula, Flinders Ranges and Outback, Kangaroo Island, Limestone Coast, Murraylands, Riverland, Yorke Peninsula.
Tasmania
East Coast, Hobart and the South, Launceston Tamar and the North, North West, Wilderness West.
Victoria
Ballarat, Bendigo Loddon, Central Highlands, Central Murray, Geelong and the Bellarine, Gippsland, Goulburn, Great Ocean Road, High Country, Lakes, Macedon, Mallee, Melbourne, Melbourne East, Murray East, Peninsula, Phillip Island, Spa Country, Upper Yarra, Western Grampians, Wimmera.
Western Australia
Australia’s Coral Coast, Australia’s Golden Outback, Australia’s North West, Australia’s South West, Experience Perth.
The tourism dataset contains data on quarterly domestic tourism demand, measured as the number of overnight trips Australians spend away from home. The key variables State and Region denote the geographical areas, while a further key Purpose describes the purpose of travel. For now, we will ignore the purpose of travel and just consider the geographic hierarchy. To make the graphs and tables simpler, we will recode State to use abbreviations.
Using the aggregate() function, we can create the hierarchical time series with overnight trips in regions at the bottom level of the hierarchy, aggregated to states, which are aggregated to the national total. A hierarchical time series corresponding to the nested structure is created using a parent/child specification.
The new Y_df now has some additional rows corresponding to state and national aggregations for each quarter, while the purpose of travel has been aggregated out. The aggregate() function also provides the summing dataframe S_df, which indicates which bottom level time series are included in each higher level time series, as well as tags, which lists the unique_id values at each level in the hierarchy.
S_df.iloc[:, :5]
unique_id
Australia
ACT
Canberra
Australia
NSW
Blue Mountains
Australia
NSW
Capital Country
Australia
NSW
Central Coast
0
Australia
1
1
1
1
1
Australia
ACT
1
0
0
0
2
Australia
NSW
0
1
1
1
3
Australia
NT
0
0
0
0
...
...
...
...
...
81
Australia
WA
Australia's Golden Outback
0
0
0
0
82
Australia
WA
Australia's North West
0
0
0
0
83
Australia
WA
Australia's South West
0
0
0
0
84
Australia
WA
Experience Perth
0
0
0
0
for (tag, unique_ids) in tags.items():print(tag)print(unique_ids, end='\n\n')
Country
['Australia']
Country/State
['Australia/ACT' 'Australia/NSW' 'Australia/NT' 'Australia/QLD'
'Australia/SA' 'Australia/TAS' 'Australia/VIC' 'Australia/WA']
Country/State/Region
['Australia/ACT/Canberra' 'Australia/NSW/Blue Mountains'
'Australia/NSW/Capital Country' ... "Australia/WA/Australia's North West"
"Australia/WA/Australia's South West" 'Australia/WA/Experience Perth']
Figure 11.3 shows the aggregate total overnight trips for the whole of Australia as well as the states, revealing diverse and rich dynamics. For example, there is noticeable national growth since 2010 and for some states such as the ACT, New South Wales, Queensland, South Australia, and Victoria. There seems to be a significant jump for Western Australia in 2014.
Figure 11.3: Domestic overnight trips from 1998 Q1 to 2017 Q4 aggregated by state.
The seasonal pattern of the northern states, such as Queensland and the Northern Territory, leads to peak visits in winter (corresponding to Q3) due to the tropical climate and rainy summer months. In contrast, the southern states tend to peak in summer (corresponding to Q1). This is highlighted in the seasonal plots shown in Figure 11.4 for Queensland and the Northern Territory (shown in the left column) versus the most southern states of Victoria and Tasmania (shown in the right column).
Figure 11.4: Seasonal plots for overnight trips for Queensland and the Northern Territory, and Victoria and Tasmania highlighting the contrast in seasonal patterns between northern and southern states in Australia.
Code
plot_regions = ["North Coast NSW", "Snowy Mountains", "Hunter","New England North West", "Alice Springs", "Darwin","Kakadu Arnhem", "MacDonnell", "Brisbane","Gold Coast", "Northern Outback", "Sunshine Coast","Tropical North Queensland", "Adelaide Hills", "Murraylands","Yorke Peninsula", "Kangaroo Island", "Ballarat","Great Ocean Road", "High Country", "Goulburn", "Australia's Coral Coast","Australia's Golden Outback", "Australia's North West","Australia's North West",]aus_tourism_rq = aus_tourism.query("Region in @plot_regions")plot_states = ["NSW", "NT", "QLD", "SA", "VIC", "WA"]fig, axs = plt.subplots(2, 3, sharex=True)for state, ax inzip(plot_states, axs.flatten()): plot_state = ( aus_tourism_rq.query("State == @state") .groupby(["Region", "State", "ds"]) .sum(numeric_only=True) .reset_index() ) sns.lineplot(data=plot_state, x="ds", y="y", hue="Region", palette="husl", ax=ax, legend=False) ax.set(title=state, xlabel="", ylabel="")fig.supylabel("Trips ('000)")fig.supxlabel("Quarter [1Q]")fig.suptitle("Australian tourism: by regions nested within states", x=.53)plt.show()
Figure 11.5: Domestic overnight trips from 1998 Q1 to 2017 Q4 for some selected regions.
The plots in Figure 11.5 show data for some selected regions. These help us to visualise the diverse regional dynamics within each state, with some series showing strong trends or seasonality, some showing contrasting seasonality, while some series appear to be just noise.
Grouped time series
With grouped time series, the data structure does not naturally disaggregate in a unique hierarchical manner. Figure 11.6 shows a simple grouped structure. At the top of the grouped structure is the Total, the most aggregate level of the data, again represented by y_t. The Total can be disaggregated by attributes (A, B) forming series y_{\text{A},t} and y_{\text{B},t}, or by attributes (X, Y) forming series y_{\text{X},t} and y_{\text{Y},t}. At the bottom level, the data are disaggregated by both attributes.
Figure 11.6: Alternative representations of a two level grouped structure.
This example shows that there are alternative aggregation paths for grouped structures. For any time t, as with the hierarchical structure,
However, for the first level of the grouped structure,
y_{\text{A},t}=y_{\text{AX},t}+y_{\text{AY},t}\quad \quad y_{\text{B},t}=y_{\text{BX},t}+y_{\text{BY},t}
\tag{11.3} but also
y_{\text{X},t}=y_{\text{AX},t}+y_{\text{BX},t}\quad \quad y_{\text{Y},t}=y_{\text{AY},t}+y_{\text{BY},t}.
\tag{11.4} Grouped time series can sometimes be thought of as hierarchical time series that do not impose a unique hierarchical structure, in the sense that the order by which the series can be grouped is not unique.
Example: Australian prison population
In this example we consider the Australia prison population data introduced in Chapter 2. The top panel in Figure 11.7 shows the total number of prisoners in Australia over the period 2005Q1–2016Q4. This represents the top-level series in the grouping structure. The panels below show the prison population disaggregated or grouped by (a) state (b) legal status (whether prisoners have already been sentenced or are in remand waiting for a sentence), and (c) gender. The three factors are crossed, but none are nested within the others.
Figure 11.7: Total Australian quarterly adult prison population, disaggregated by state, by legal status, and by gender.
Figure 11.8 shows the Australian prison population grouped by all possible combinations of two attributes at a time: state and gender, state and legal status, and legal status and gender.
Figure 11.9: Bottom-level time series for the Australian adult prison population, grouped by state, legal status and gender.
Mixed hierarchical and grouped structure
Often disaggregating factors are both nested and crossed. For example, the Australian tourism data can also be disaggregated by the four purposes of travel: holiday, business, visiting friends and relatives, and other. This grouping variable does not nest within any of the geographical variables. In fact, we could consider overnight trips split by purpose of travel for the whole of Australia, and for each state, and for each region. We describe such a structure as a “nested” geographic hierarchy “crossed” with the purpose of travel. Using aggregate() this can be specified by simply combining the factors.
The aus_tourism dataset contains 425 series, including the 85 series from the hierarchical structure, as well as another 340 series obtained when each series of the hierarchical structure is crossed with the purpose of travel.
Figure 11.11: Australian domestic overnight trips over the period 1998 Q1 to 2017 Q4 disaggregated by purpose of travel and by state.
Figures 11.10 and 11.11 show the aggregate series grouped by purpose of travel, and the series grouped by purpose of travel and state, revealing further rich and diverse dynamics across these series.
11.2 Single level approaches
Traditionally, forecasts of hierarchical or grouped time series involved selecting one level of aggregation and generating forecasts for that level. These are then either aggregated for higher levels, or disaggregated for lower levels, to obtain a set of coherent forecasts for the rest of the structure.
The bottom-up approach
A simple method for generating coherent forecasts is the “bottom-up” approach. This approach involves first generating forecasts for each series at the bottom level, and then summing these to produce forecasts for all the series in the structure.
For example, for the hierarchy of Figure 11.1, we first generate h-step-ahead forecasts for each of the bottom-level series:
\hat{y}_{\text{AA},h},~~\hat{y}_{\text{AB},h},~~\hat{y}_{\text{AC},h},~~ \hat{y}_{\text{BA},h}~~\text{and}~~\hat{y}_{\text{BB},h}.
(We have simplified the previously used notation of \hat{y}_{T+h|T} for brevity.)
Summing these, we get h-step-ahead coherent forecasts for the rest of the series: \begin{align*}
\tilde{y}_{h} & =\hat{y}_{\text{AA},h}+\hat{y}_{\text{AB},h}+\hat{y}_{\text{AC},h}+\hat{y}_{\text{BA},h}+\hat{y}_{\text{BB},h}, \\
\tilde{y}_{\text{A},h} & = \hat{y}_{\text{AA},h}+\hat{y}_{\text{AB},h}+\hat{y}_{\text{AC},h}, \\
\text{and}\quad
\tilde{y}_{\text{B},h} &= \hat{y}_{\text{BA},h}+\hat{y}_{\text{BB},h}.
\end{align*} (In this chapter, we will use the “tilde” notation to indicate coherent forecasts.)
An advantage of this approach is that we are forecasting at the bottom level of a structure, and therefore no information is lost due to aggregation. On the other hand, bottom-level data can be quite noisy and more challenging to model and forecast.
Example: Generating bottom-up forecasts
Suppose we want national and state forecasts for the Australian tourism data, but we aren’t interested in disaggregations using regions or the purpose of travel. So we first create a simple dataset containing only state and national trip totals for each quarter.
However, we want a more general approach that will work with all the forecasting methods discussed in this chapter. So we will use the reconcilers object to specify how we want to compute coherent forecasts.
The reconcilers object can be used in the HierarchicalReconciliation class to produce a set of hierarchically coherent forecasts. The Y_rec_df dataframe contains the AutoETS() forecasts as well as the coherent AutoETS()/BottomUp() forecasts, for the 8 states and the national aggregate. At the state level, these forecasts are identical, but the national AutoETS() forecasts will be different from the national AutoETS()/BottomUp() forecasts.
For bottom-up forecasting, this is rather inefficient as we are not interested in the ETS model for the national total, and the resulting Y_rec_df contains a lot of duplicates. But later we will introduce more advanced methods where we will need models for all levels of aggregation, and where the coherent forecasts are different from any of the original forecasts.
Workflow for forecasting aggregation structures
The code below illustrates the general workflow for hierarchical and grouped forecasts.
# Load datadata = ...# Define aggregation structure and build Y_df objectspec = []Y_df, S_df, tags = aggregate(df=data, spec=spec)# Choose a set of models and create the base forecastssf = StatsForecast(models=[.....], freq=..., n_jobs=-1)Y_hat_df = sf.forecast(h=..., df=Y_df)# Choose a set of reconcilers and reconcile the base forecastsreconcilers = [BottomUp()]hrec = HierarchicalReconciliation(reconcilers=reconcilers)Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df, Y_df=Y_df, S_df=S_df, tags=tags)
Begin with a DataFrame (here labelled data) containing the individual bottom-level series.
Define in spec the aggregation structure and build a Y_df object that also contains the aggregate series.
Specify a set of models for each series, at all levels of aggregation.
Specify in reconcilers how the coherent forecasts are to be generated from the selected models.
Use the HierarchicalReconciliation class with the reconcile() function to generate forecasts for the whole aggregation structure.
Top-down approaches
Top-down approaches involve first generating forecasts for the Total series y_t, and then disaggregating these down the hierarchy.
Let p_1,\dots,p_{m} denote a set of disaggregation proportions which determine how the forecasts of the Total series are to be distributed to obtain forecasts for each series at the bottom level of the structure. For example, for the hierarchy of Figure 11.1, using proportions p_1,\dots,p_{5} we get
\tilde{y}_{\text{AA},t}=p_1\hat{y}_t,~~~\tilde{y}_{\text{AB},t}=p_2\hat{y}_t,~~~\tilde{y}_{\text{AC},t}=p_3\hat{y}_t,~~~\tilde{y}_{\text{BA},t}=p_4\hat{y}_t~~~\text{and}~~~~~~\tilde{y}_{\text{BB},t}=p_5\hat{y}_t.
Once the bottom-level h-step-ahead forecasts have been generated, these are aggregated to generate coherent forecasts for the rest of the series.
Top-down forecasts can be generated using TopDown() within the reconcilers object.
There are several possible top-down methods that can be specified. The two most common top-down approaches specify disaggregation proportions based on the historical proportions of the data. These performed well in the study of Gross & Sohl (1990).
Average historical proportions
p_j=\frac{1}{T}\sum_{t=1}^{T}\frac{y_{j,t}}{{y_t}}
for j=1,\dots,m. Each proportion p_j reflects the average of the historical proportions of the bottom-level series y_{j,t} over the period t=1,\dots,T relative to the total aggregate y_t.
This approach is implemented in the TopDown() class by setting method ="average_proportions".
Proportions of the historical averages
p_j={\sum_{t=1}^{T}\frac{y_{j,t}}{T}}\Big/{\sum_{t=1}^{T}\frac{y_t}{T}}
for j=1,\dots,m. Each proportion p_j captures the average historical value of the bottom-level series y_{j,t} relative to the average value of the total aggregate y_t.
This approach is implemented in the TopDown() class by setting method ="proportion_averages".
A convenient attribute of such top-down approaches is their simplicity. One only needs to model and generate forecasts for the most aggregated top-level series. In general, these approaches seem to produce quite reliable forecasts for the aggregate levels and they are useful with low count data. On the other hand, one disadvantage is the loss of information due to aggregation. Using such top-down approaches, we are unable to capture and take advantage of individual series characteristics such as time dynamics, special events, different seasonal patterns, etc.
Forecast proportions
Because historical proportions used for disaggregation do not take account of how those proportions may change over time, top-down approaches based on historical proportions tend to produce less accurate forecasts at lower levels of the hierarchy than bottom-up approaches. To address this issue, proportions based on forecasts rather than historical data can be used (Athanasopoulos et al., 2009).
Consider a one level hierarchy. We first generate h-step-ahead forecasts for all the series. We don’t use these forecasts directly, and they are not coherent (they don’t add up correctly). Let’s call these “initial” forecasts. We calculate the proportion of each h-step-ahead initial forecast at the bottom level, to the aggregate of all the h-step-ahead initial forecasts at this level. We refer to these as the forecast proportions, and we use them to disaggregate the top-level h-step-ahead initial forecast in order to generate coherent forecasts for the whole of the hierarchy.
For a K-level hierarchy, this process is repeated for each node, going from the top to the bottom level. Applying this process leads to the following general rule for obtaining the forecast proportions:
p_j=\prod^{K-1}_{\ell=0}\frac{\hat{y}_{j,h}^{(\ell)}}{\hat{S}_{j,h}^{(\ell+1)}}
where j=1,2,\dots,m, \hat{y}_{j,h}^{(\ell)} is the h-step-ahead initial forecast of the series that corresponds to the node which is \ell levels above j, and \hat{S}_{j,h}^{(\ell)} is the sum of the h-step-ahead initial forecasts below the node that is \ell levels above node j and are directly connected to that node. These forecast proportions disaggregate the h-step-ahead initial forecast of the Total series to get h-step-ahead coherent forecasts of the bottom-level series.
We will use the hierarchy of Figure 11.1 to explain this notation and to demonstrate how this general rule is reached. Assume we have generated initial forecasts for each series in the hierarchy. Recall that for the top-level “Total” series, \tilde{y}_{h}=\hat{y}_{h}, for any top-down approach. Here are some examples using the above notation:
Moving down the farthest left branch of the hierarchy, coherent forecasts are given by
\tilde{y}_{\text{A},h} = \Bigg(\frac{\hat{y}_{\text{A},h}}{\hat{S}_{\text{A},h}^{(1)}}\Bigg) \tilde{y}_{h} =
\Bigg(\frac{\hat{y}_{\text{AA},h}^{(1)}}{\hat{S}_{\text{AA},h}^{(2)}}\Bigg) \tilde{y}_{h}
and
\tilde{y}_{\text{AA},h} = \Bigg(\frac{\hat{y}_{\text{AA},h}}{\hat{S}_{\text{AA},h}^{(1)}}\Bigg) \tilde{y}_{\text{A},h}
=\Bigg(\frac{\hat{y}_{\text{AA},h}}{\hat{S}_{\text{AA},h}^{(1)}}\Bigg) \Bigg(\frac{\hat{y}_{\text{AA},h}^{(1)}}{\hat{S}_{\text{AA},h}^{(2)}}\Bigg)\tilde{y}_{h}.
Consequently,
p_1=\Bigg(\frac{\hat{y}_{\text{AA},h}}{\hat{S}_{\text{AA},h}^{(1)}}\Bigg) \Bigg(\frac{\hat{y}_{\text{AA},h}^{(1)}}{\hat{S}_{\text{AA},h}^{(2)}}\Bigg).
The other proportions can be obtained similarly.
This approach is implemented in the TopDown() class by setting method ="forecast_proportions". This approach tends to work better than other top-down methods.
One disadvantage of all top-down approaches, is that they may not produce unbiased coherent forecasts (Hyndman et al., 2011) even if the base forecasts are unbiased.
Middle-out approach
The middle-out approach combines bottom-up and top-down approaches. Again, it can only be used for strictly hierarchical aggregation structures.
First, a “middle” level is chosen and forecasts are generated for all the series at this level. For the series above the middle level, coherent forecasts are generated using the bottom-up approach by aggregating the “middle-level” forecasts upwards. For the series below the “middle level”, coherent forecasts are generated using a top-down approach by disaggregating the “middle level” forecasts downwards.
This approach is implemented in the MiddleOut() class by specifying the appropriate middle level via the middle_level argument and selecting the top-down approach with the top_down_method argument.
11.3 Forecast reconciliation
Warning: the rest of this chapter is more advanced and assumes a knowledge of some basic matrix algebra.
Matrix notation
Recall that Equations 11.1 and 11.2 represent how data, that adhere to the hierarchical structure of Figure 11.1, aggregate. Similarly Equation 11.3 and Equation 11.4 represent how data, that adhere to the grouped structure of Figure 11.6, aggregate. These equations can be thought of as aggregation constraints or summing equalities, and can be more efficiently represented using matrix notation.
For any aggregation structure we construct an n\times m matrix \bm{S} (referred to as the “summing matrix”) which dictates the way in which the bottom-level series aggregate.
For the hierarchical structure in Figure 11.1, we can write
\begin{bmatrix}
y_{t} \\
y_{\text{A},t} \\
y_{\text{B},t} \\
y_{\text{AA},t} \\
y_{\text{AB},t} \\
y_{\text{AC},t} \\
y_{\text{BA},t} \\
y_{\text{BB},t}
\end{bmatrix}
=
\begin{bmatrix}
1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 \\
1 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
y_{\text{AA},t} \\
y_{\text{AB},t} \\
y_{\text{AC},t} \\
y_{\text{BA},t} \\
y_{\text{BB},t}
\end{bmatrix}
or in more compact notation
\bm{y}_t=\bm{S}\bm{b}_{t},
\tag{11.5} where \bm{y}_t is an n-dimensional vector of all the observations in the hierarchy at time t, \bm{S} is the summing matrix, and \bm{b}_{t} is an m-dimensional vector of all the observations in the bottom level of the hierarchy at time t. Note that the first row in the summing matrix \bm{S} represents Equation 11.1, the second and third rows represent Equation 11.2. The rows below these comprise an m-dimensional identity matrix \bm{I}_m so that each bottom-level observation on the right hand side of the equation is equal to itself on the left hand side.
Similarly for the grouped structure of Figure 11.6 we write
This matrix notation allows us to represent all forecasting methods for hierarchical or grouped time series using a common notation.
Suppose we forecast all series ignoring any aggregation constraints. We call these the base forecasts and denote them by \hat{\bm{y}}_h where h is the forecast horizon. They are stacked in the same order as the data \bm{y}_t.
Then all coherent forecasting approaches for either hierarchical or grouped structures can be represented as1
\tilde{\bm{y}}_h=\bm{S}\bm{G}\hat{\bm{y}}_h,
\tag{11.7} where \bm{G} is a matrix that maps the base forecasts into the bottom level, and the summing matrix \bm{S} sums these up using the aggregation structure to produce a set of coherent forecasts\tilde{\bm{y}}_h.
The \bm{G} matrix is defined according to the approach implemented. For example if the bottom-up approach is used to forecast the hierarchy of Figure 11.1, then
\bm{G}=
\begin{bmatrix}
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\
\end{bmatrix}.
Notice that \bm{G} contains two partitions. The first three columns zero out the base forecasts of the series above the bottom level, while the m-dimensional identity matrix picks only the base forecasts of the bottom level. These are then summed by the \bm{S} matrix.
If any of the top-down approaches were used then
\bm{G}=
\begin{bmatrix}
p_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
p_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
p_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
p_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
p_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
\end{bmatrix}.
The first column includes the set of proportions that distribute the base forecasts of the top level to the bottom level. These are then summed up by the \bm{S} matrix. The rest of the columns zero out the base forecasts below the highest level of aggregation.
For a middle out approach, the \bm{G} matrix will be a combination of the above two. Using a set of proportions, the base forecasts of some pre-chosen level will be disaggregated to the bottom level, all other base forecasts will be zeroed out, and the bottom-level forecasts will then be summed up the hierarchy via the summing matrix.
Forecast reconciliation
Equation 11.7 shows that pre-multiplying any set of base forecasts with \bm{S}\bm{G} will return a set of coherent forecasts.
The traditional methods considered so far are limited in that they only use base forecasts from a single level of aggregation which have either been aggregated or disaggregated to obtain forecasts at all other levels. Hence, they use limited information. However, in general, we could use other \bm{G} matrices, and then \bm{S}\bm{G} combines and reconciles all the base forecasts in order to produce coherent forecasts.
In fact, we can find the optimal \bm{G} matrix to give the most accurate reconciled forecasts.
The MinT optimal reconciliation approach
Wickramasuriya et al. (2019) found a \bm{G} matrix that minimises the total forecast variance of the set of coherent forecasts, leading to the MinT (Minimum Trace) optimal reconciliation approach.
Suppose we generate coherent forecasts using Equation 11.7. First we want to make sure we have unbiased forecasts. If the base forecasts \hat{\bm{y}}_h are unbiased, then the coherent forecasts \tilde{\bm{y}}_h will be unbiased provided2\bm{S}\bm{G}\bm{S}=\bm{S}. This provides a constraint on the matrix \bm{G}. Interestingly, no top-down method satisfies this constraint, so all top-down approaches result in biased coherent forecasts.
Next we need to find the errors in our forecasts. Wickramasuriya et al. (2019) show that the variance-covariance matrix of the h-step-ahead coherent forecast errors is given by \begin{equation*}
\bm{V}_h = \text{Var}[\bm{y}_{T+h}-\tilde{\bm{y}}_h]=\bm{S}\bm{G}\bm{W}_h\bm{G}'\bm{S}'
\end{equation*} where \bm{W}_h=\text{Var}[(\bm{y}_{T+h}-\hat{\bm{y}}_h)] is the variance-covariance matrix of the corresponding base forecast errors.
The objective is to find a matrix \bm{G} that minimises the error variances of the coherent forecasts. These error variances are on the diagonal of the matrix \bm{V}_h, and so the sum of all the error variances is given by the trace of the matrix \bm{V}_h. Wickramasuriya et al. (2019) show that the matrix \bm{G} which minimises the trace of \bm{V}_h such that \bm{S}\bm{G}\bm{S}=\bm{S}, is given by
\bm{G}=(\bm{S}'\bm{W}_h^{-1}\bm{S})^{-1}\bm{S}'\bm{W}_h^{-1}.
Therefore, the optimally reconciled forecasts are given by
\tilde{\bm{y}}_h=\bm{S}(\bm{S}'\bm{W}_h^{-1}\bm{S})^{-1}\bm{S}'\bm{W}_h^{-1}\hat{\bm{y}}_h.
\tag{11.8} We refer to this as the MinT (or Minimum Trace) optimal reconciliation approach. MinT is implemented by the MinTrace() class.
To use this in practice, we need to estimate \bm{W}_h, the forecast error variance of the h-step-ahead base forecasts. This can be difficult, and so we provide four simplifying approximations that have been shown to work well in both simulations and in practice.
Set \bm{W}_h=k_h\bm{I} for all h, where k_{h} > 0.3 This is the most simplifying assumption to make, and means that \bm{G} is independent of the data, providing substantial computational savings. The disadvantage, however, is that this specification does not account for the differences in scale between the levels of the structure, or for relationships between series.
Setting \bm{W}_h=k_h\bm{I} in Equation 11.8 gives the ordinary least squares (OLS) estimator we introduced in Section 7.9 with \bm{X}=\bm{S} and \bm{y}=\hat{\bm{y}}. Hence this approach is usually referred to as OLS reconciliation. It is implemented in MinTrace() by setting method ="ols".
Set \bm{W}_{h} = k_{h}\text{diag}(\hat{\bm{W}}_{1}) for all h, where k_{h} > 0, \hat{\bm{W}}_{1} = \frac{1}{T}\sum_{t=1}^{T}\bm{e}_{t}\bm{e}_{t}', and \bm{e}_{t} is an n-dimensional vector of residuals of the models that generated the base forecasts stacked in the same order as the data.
This specification scales the base forecasts using the variance of the residuals and it is therefore referred to as the WLS (weighted least squares) estimator using variance scaling. The approach is implemented in MinTrace() by setting method ="wls_var".
Set \bm{W}_{h}=k_{h}\bm{\Lambda} for all h, where k_{h} > 0, \bm{\Lambda}=\text{diag}(\bm{S}\bm{1}), and \bm{1} is a unit vector of dimension m (the number of bottom-level series). This specification assumes that the bottom-level base forecast errors each have variance k_{h} and are uncorrelated between nodes. Hence each element of the diagonal \bm{\Lambda} matrix contains the number of forecast error variances contributing to each node. This estimator only depends on the structure of the aggregations, and not on the actual data. It is therefore referred to as structural scaling. Applying the structural scaling specification is particularly useful in cases where residuals are not available, and so variance scaling cannot be applied; for example, in cases where the base forecasts are generated by judgmental forecasting (Chapter 6). The approach is implemented in MinTrace() by setting method ="wls_struct".
Set \bm{W}_h = k_h \bm{W}_1 for all h, where k_h>0. Here we only assume that the error covariance matrices are proportional to each other, and we directly estimate the full one-step covariance matrix \bm{W}_1. The most obvious and simple way would be to use the sample covariance. This is implemented in MinTrace() by setting method ="mint_cov".
However, for cases where the number of bottom-level series m is large compared to the length of the series T, this is not a good estimator. Instead we use a shrinkage estimator which shrinks the sample covariance to a diagonal matrix. This is implemented in MinTrace() by setting method ="mint_shrink".
In summary, unlike any other existing approach, the optimal reconciliation forecasts are generated using all the information available within a hierarchical or a grouped structure. This is important, as particular aggregation levels or groupings may reveal features of the data that are of interest to the user and are important to be modelled. These features may be completely hidden or not easily identifiable at other levels.
For example, consider the Australian tourism data introduced in Section 11.1, where the hierarchical structure followed the geographic division of a country into states and regions. Some areas will be largely summer destinations, while others may be winter destinations. We saw in Figure 11.4 the contrasting seasonal patterns between the northern and the southern states. These differences will be smoothed at the country level due to aggregation.
11.4 Forecasting Australian domestic tourism
We will compute forecasts for the Australian tourism data that was described in Section 11.1. We use the data up to the end of 2015 as a training set, withholding the final two years (eight quarters, 2016Q1–2017Q4) as a test set for evaluation. The code below demonstrates the full workflow for generating coherent forecasts using the bottom-up, OLS and MinT methods. The workflow for this example is:
Begin with a DataFrame (here labelled aus_tourism) containing the individual bottom-level series.
Define in spec the aggregation structure and build a Y_df object that also contains the aggregate series.
Specify a set of models for each series, at all levels of aggregation. In this example we use AutoETS().
Specify in reconcilers how the coherent forecasts are to be generated from the selected models. In this example, we use BottomUp() and two versions of MinTrace().
Use the HierarchicalReconciliation class with the reconcile() function to generate forecasts for the whole aggregation structure.
# 2. Define the spec and aggregate the dataspec = [ ["Country"], ["Country", "State"], ["Country", "Purpose"], ["Country", "State", "Region"], ["Country", "State", "Purpose"], ["Country", "State", "Region", "Purpose"],]Y_df, S_df, tags = aggregate( aus_tourism.drop(columns=["unique_id"]), spec)Y_train_df = Y_df.loc[lambda x: x["ds"] <"2016"]Y_test_df = Y_df.loc[lambda x: x["ds"] >="2016"]# 3. Specify a set of models and forecastsf = StatsForecast(models=[AutoETS(season_length=4)], freq="Q", n_jobs=-1)Y_hat_df = sf.forecast(h=8, df=Y_train_df, fitted=True)Y_fitted_df = sf.forecast_fitted_values()# 4. Specify the reconcilersreconcilers = [ BottomUp(), MinTrace(method="ols"), MinTrace(method="mint_shrink"),]# 5. Reconcile the forecastshrec = HierarchicalReconciliation(reconcilers=reconcilers)Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df, Y_df=Y_fitted_df, S_df=S_df, tags=tags)
Figures 11.12 and 11.13 plot the four point forecasts for the overnight trips for the Australian total, the states, and the purposes of travel, along with the actual observations of the test set.
Code
# Helper for simplifying column namesdef rename_models(df, labels):"""Apply at most one string replacement to each column"""return df.rename(columns=lambda c: next( (c.replace(k, v) for k, v in labels.items() if k in c), c))renames = {"AutoETS/MinTrace_method-mint_shrink": "MinT","AutoETS/MinTrace_method-ols": "OLS","AutoETS/BottomUp": "BottomUp","AutoETS": "Base",}Y_rec_df = rename_models(Y_rec_df, renames)plot_tags = [*tags["Country"], *tags["Country/State"]]plot_Y_df = Y_df.loc[lambda x: x["unique_id"].isin(plot_tags)& (x["ds"] >="2011")]plot_Y_rec_df = Y_rec_df.loc[lambda x: x["unique_id"].isin(plot_tags)]fig, axes = plt.subplots(3, 3, figsize=(8, 7), sharex=True)plot_series(plot_Y_df, plot_Y_rec_df, xlabel="", ylabel="", palette="black_and_4color", ax=axes)for ax in axes[-1]: plt.setp(ax.get_xticklabels(), rotation=90, ha="right")handles, labels = axes.flat[-1].get_legend_handles_labels()fig.legend(handles, labels, bbox_to_anchor=(1.02, 0.5), loc='center left', borderaxespad=0, frameon=False)fig.supylabel("Trips ('000)")fig.supxlabel("Quarter [1Q]", x=.55)fig
Figure 11.12: Forecasts of overnight trips for Australia and its states over the test period 2016Q1–2017Q4.
Figure 11.13: Forecasts of overnight trips by purpose of travel over the test period 2016Q1–2017Q4.
To make it easier to see the differences, we have included only the last five years of the training data, and have omitted the prediction intervals. In most panels, the increase in overnight trips, especially in the second half of the test set, is higher than what is predicted by the point forecasts. This is particularly noticeable for the mainland eastern states of ACT, New South Wales, Queensland and Victoria, and across all purposes of travel.
The accuracy of the forecasts over the test set can be evaluated using the evaluate() function. We summarise some results in Table 11.2 using RMSE and MASE.
The scales of the series at different levels of aggregation are quite different, due to aggregation. Hence, we need to be cautious when comparing or calculating scale dependent error measures, such as the RMSE, across levels as the aggregate series will dominate. Therefore, we compare error measures across each level of aggregation, before providing the error measures across all the series in the bottom-row. Notice, that the RMSE increases as we go from the bottom level to the aggregate levels above.
The following code generates the accuracy measures shown in Table 11.2.
# Define a set of evaluation tagseval_tags = {"Total": tags["Country"],"Purpose": tags["Country/Purpose"],"State": tags["Country/State"],"Regions": tags["Country/State/Region"],"Bottom": tags["Country/State/Region/Purpose"],}# Run the evaluationeval_df = Y_rec_df.merge(Y_test_df, on=["unique_id", "ds"])evaluation = evaluate( df=eval_df, tags=eval_tags, train_df=Y_train_df, metrics=[rmse, partial(mase, seasonality=4)],)# Format the evaluation tablemodel_names ="Base", "BottomUp", "MinT", "OLS"evaluation = ( evaluation.pivot_table(index="level", columns="metric", values=model_names, sort=False) .rename(index={"Overall": "All series"}) .reorder_levels([1, 0], axis="columns") .sort_index(axis="columns") .rename(columns={"mase": "MASE", "rmse": "RMSE"}))evaluation
Table 11.2: Accuracy of forecasts for Australian overnight trips over the test set 2016Q1–2017Q4.
metric
MASE
RMSE
Base
BottomUp
MinT
OLS
Base
BottomUp
MinT
OLS
level
Total
1.53
3.09
1.93
1.60
1713
2988
2009
1780
Purpose
1.30
2.17
1.41
1.22
524
784
553
502
State
1.31
1.85
1.32
1.19
298
407
306
284
Regions
1.11
1.18
0.99
0.99
51
54
46
46
Bottom
0.99
0.99
0.94
1.02
19
19
18
18
All series
1.02
1.06
0.96
1.02
41
49
39
38
Reconciling the base forecasts using OLS and MinT results in more accurate forecasts compared to the bottom-up approach. This result is commonly observed in applications as reconciliation approaches use information from all levels of the structure, resulting in more accurate coherent forecasts compared to the older traditional methods which use limited information. Furthermore, reconciliation usually improves the incoherent base forecasts for almost all levels.
11.5 Reconciled distributional forecasts
So far we have only discussed the reconciliation of point forecasts. However, we are usually also interested in the forecast distributions so that we can compute prediction intervals.
Panagiotelis et al. (2023) present several important results for generating reconciled probabilistic forecasts. We focus here on two fundamental results.
If the base forecasts are normally distributed, i.e.,
\hat{\bm{y}}_h\sim N(\hat{\bm\mu}_h,\hat{\bm\Sigma}_h),
then the reconciled forecasts are also normally distributed,
\tilde{\bm{y}}_h \sim N(\bm{S}\bm{G}\hat{\bm{\mu}}_h,\bm{S}\bm{G}\hat{\bm{\Sigma}}_{h}\bm{G}'\bm{S}').
To generate prediction intervals in this way, we simply set intervals_method ='normality' in the reconcile() function. This is the default choice.
If it is unreasonable to assume normality for the base forecasts, we can use bootstrapping. Bootstrapped prediction intervals were introduced in Section 5.5. The same idea can be used here. We can simulate future sample paths from the model(s) that produce the base forecasts, and then reconcile these sample paths. Coherent prediction intervals can be computed from the reconciled sample paths.
Suppose that (\hat{\bm{y}}_h^{[1]},\dots,\hat{\bm{y}}_h^{[B]}) are a set of B simulated sample paths, generated independently from the models used to produce the base forecasts. Then (\bm{S}\bm{G}\hat{\bm{y}}_h^{[1]},\dots,\bm{S}\bm{G}\hat{\bm{y}}_h^{[B]}) provides a set of reconciled sample paths, from which percentiles can be calculated in order to construct coherent prediction intervals.
To generate bootstrapped prediction intervals in this way, we simply set intervals_method ='bootstrap' in the reconcile() function.
11.6 Forecasting Australian prison population
Returning to the Australian prison population data (Section 11.1), we will compare the forecasts from bottom-up and MinT methods applied to base ETS models, using a test set comprising the final two years or eight quarters 2015Q1–2016Q4 of the available data. We will also forecast using the Naive model, to compare the results against later on.
Figure 11.14: Forecasts for the total Australian quarterly adult prison population for the period 2015Q1–2016Q4.
Figure 11.14 shows the three sets of forecasts for the aggregate Australian prison population. The base and bottom-up forecasts from the ETS models seem to underestimate the trend over the test period. The MinT approach combines information from all the base forecasts in the aggregation structure; in this case, the base forecasts at the top level are adjusted upwards.
The MinT reconciled prediction intervals are much tighter than the base forecasts, due to MinT being based on an estimator that minimises variances. The base forecast distributions are also incoherent, and therefore carry with them the extra uncertainty of the incoherency error.
We exclude the bottom-up forecasts from the remaining plots in order to simplify the visual exploration. However, we do revisit their accuracy in the evaluation results presented later.
Figures 11.15–11.18 show the MinT and base forecasts at various levels of aggregation. To make it easier to see the effect, we only show the last five years of training data. In general, MinT adjusts the base forecasts in the direction of the test set, hence improving the forecast accuracy. There is no guarantee that MinT reconciled forecasts will be more accurate than the base forecasts for every series, but they will be more accurate on average (see Panagiotelis et al., 2021).
Figure 11.15: Forecasts for the Australian quarterly adult prison population, disaggregated by state.
Figure 11.15 shows forecasts for each of the eight states. There is a general upward trend during the test set period across all the states. However, there appears to be a relatively large and sudden surge in New South Wales and Tasmania, which means the test set observations are well outside the upper bound of the forecast intervals for both these states. Because New South Wales is the state with the largest prison population, this surge will have a substantial impact on the total. In contrast, Victoria shows a substantial dip in 2015Q2–2015Q3, before returning to an upward trend. This dip is not captured in any of the Victorian forecasts.
Figure 11.18: Forecasts for bottom-level series of the Australian quarterly adult prison population, disaggregated by state, by legal status and by gender.
Figure 11.18 shows the forecasts for some selected bottom-level series of the Australian prison population. The four largest states are represented across the columns, with legal status and gender down the rows. These allow for some interesting analysis and observations that have policy implications. The large increase observed across the states during the 2015Q1–2016Q4 test period appears to be driven by large increases in the remand prison population. These increases seem to be generally missed by both forecasts. In contrast to the other states, for New South Wales there is also a substantial increase in the sentenced prison population. In particular, the increase in numbers of sentenced males in NSW contributes substantially to the rise in state and national prison numbers.
Using the evaluate() function, we evaluate the forecast accuracy across the grouped structure. The results are shown in Table 11.3, summarising the accuracy of the base, bottom-up and the MinT reconciled forecasts over the 2015Q1–2016Q4 test period across each of the levels of the grouped aggregation structure as well as all the levels.
Code
# Define a set of evaluation tagseval_tags = {"Total": tags["Country"],"State": tags["Country/State"],"Legal status": tags["Country/Legal"],"Gender": tags["Country/Gender"],"Bottom": tags["Country/State/Gender/Legal"],}# Run the evaluation and reshapeeval_df = Y_rec_df.merge(Y_test_df, on=["unique_id", "ds"])evaluation = evaluate( df=eval_df, tags=eval_tags, train_df=Y_train_df, metrics=[partial(mase, seasonality=4), mqloss], level=levels,)evaluation = ( evaluation.pivot_table(index="level", columns="metric", values=model_names, sort=False) .rename(index={"Overall": "All series"}) .reorder_levels([1, 0], axis="columns") .sort_index(axis="columns"))# Calculate CRPS skill score using Naive's CRPS score as the basecrps_ets = evaluation['mqloss'].filter(like="AutoETS").to_numpy()crps_naive = evaluation['mqloss'].filter(like="Naive").to_numpy()crps_skill_score =100* (crps_naive - crps_ets) / crps_naive# Keep only AutoETS entries and prettify tablerenames_final = {f"AutoETS/{name}": namefor name in ("Base", "BottomUp", "MinT")}evaluation_final = ( evaluation.filter(like="AutoETS") .rename(columns={"mqloss": "Skill Score (CRPS)", "mase": "MASE"}) .reorder_levels([1, 0], axis="columns").sort_index(axis="columns") .rename(columns=renames_final) .reorder_levels([1, 0], axis="columns").sort_index(axis="columns"))evaluation_final.loc[:, "Skill Score (CRPS)"] = crps_skill_score
Table 11.3: Accuracy of Australian prison population forecasts for different groups of series.
metric
MASE
Skill Score (CRPS)
Base
BottomUp
MinT
Base
BottomUp
MinT
level
Total
1.64
1.51
1.54
43.9
47.9
47.3
State
1.67
1.64
1.67
25.7
26.5
25.8
Legal status
2.55
2.57
2.49
26.7
28.8
28.9
Gender
1.37
1.30
1.34
48.8
48.8
47.5
Bottom
2.11
2.11
2.10
14.6
14.6
14.9
All series
2.01
2.00
1.99
31.0
32.7
32.1
We use scaled measures because the numbers of prisoners vary substantially across the groups. The MASE gives a scaled measure of point-forecast accuracy (see Section 5.8), while the CRPS skill score gives a scaled measure of distributional forecast accuracy (see Section 5.9). A low value of MASE indicates a good forecast, while a high value of the skill score indicates a good forecast.
The results show that the MinT reconciled forecasts improve on the accuracy of the base forecasts and are also more accurate than the bottom-up forecasts. As the MinT optimal reconciliation approach uses information from all levels in the structure, it generates more accurate forecasts than the traditional approaches (such as bottom-up) which use limited information.
11.7 Exercises
Consider the PBS data which has aggregation structure ATC1/ATC2 * Concession * Type.
Produce plots of the aggregated Scripts data by Concession, Type and ATC1.
Forecast the PBS Scripts data using ETS, ARIMA and Seasonal Naïve models, applied to all but the last three years of data.
Reconcile each of the forecasts using MinT.
Which type of model works best on the test set?
Does the reconciliation improve the forecast accuracy?
Why doesn’t the reconciliation make any difference to the Seasonal Naïve forecasts?
Repeat the tourism example from Section 11.4, but also evaluate the forecast distribution accuracy using CRPS scores. Which method does best on this measure?
Repeat the prison example from Section 11.6, but using a bootstrap to generate the forecast distributions rather than assuming normality. Does it make much difference to the CRPS scores?
11.8 Further reading
There are no other textbooks which cover hierarchical forecasting in any depth, so interested readers will need to tackle the original research papers for further information.
Gross & Sohl (1990) provide a good introduction to the top-down approaches.
A recent survey of forecast reconciliation is provided by Athanasopoulos et al. (2020).
The reconciliation methods were developed in a series of papers. The later papers summarise previous results and present the most general theory: Wickramasuriya et al. (2019), Panagiotelis et al. (2021), Panagiotelis et al. (2023).
Athanasopoulos et al. (2017) extends the reconciliation approach to deal with temporal hierarchies.
The tourism example is discussed in more detail in Athanasopoulos et al. (2009), Wickramasuriya et al. (2019), and Kourentzes & Athanasopoulos (2019).
AutoETS()model - For automatic exponential smoothing
HierarchicalForecast
HierarchicalReconciliationclass - For reconciling hierarchical forecasts
BottomUp(), TopDown(), MiddleOut(), MinTrace()methods - Different reconciliation approaches
aggregate()utility - For creating hierarchical structures
HierarchicalPlot()utility - For visualising hierarchical forecasts
evaluate()utility - For evaluating hierarchical forecasts
UtilsForecast
plot_series()utility - For creating time series visualisations
rmse(), mae(), mase(), scaled_crps()metrics - For forecast evaluation
Athanasopoulos, G, Ahmed, RA, & Hyndman, RJ. (2009). Hierarchical forecasts for Australian domestic tourism. International Journal of Forecasting, 25, 146–166.
Athanasopoulos, G, Gamakumara, P, Panagiotelis, A, Hyndman, RJ, & Affan, M. (2020). Hierarchical forecasting. In P Fuleky (Ed.), Macroeconomic forecasting in the era of big data (pp. 689–719). Springer.
Athanasopoulos, G, Hyndman, RJ, Kourentzes, N, & Petropoulos, F. (2017). Forecasting with temporal hierarchies. European Journal of Operational Research, 262(1), 60–74.
Garza, A, Mergenthaler Canseco, M, Challú, C, & Olivares, KG. (2022). StatsForecast: Lightning fast forecasting with statistical and econometric models. PyCon Salt Lake City, Utah, US 2022.
Gross, CW, & Sohl, JE. (1990). Disaggregation methods to expedite product line forecasting. Journal of Forecasting, 9, 233–254.
Hyndman, RJ, Ahmed, RA, Athanasopoulos, G, & Shang, HL. (2011). Optimal combination forecasts for hierarchical time series. Computational Statistics and Data Analysis, 55(9), 2579–2589.
Kourentzes, N, & Athanasopoulos, G. (2019). Cross-temporal coherent forecasts for Australian tourism. Annals of Tourism Research, 75, 393–409.
Olivares, KG, Garza, A, Luo, D, Challú, C, Mergenthaler, M, Ben Taieb, S, Wickramasuriya, SL, & Dubrawski, A. (2024). HierarchicalForecast: A reference framework for hierarchical forecasting in python.
Panagiotelis, A, Gamakumara, P, Athanasopoulos, G, & Hyndman, RJ. (2021). Forecast reconciliation: A geometric view with new insights on bias correction. International Journal of Forecasting, 37(1), 343–359.
Panagiotelis, A, Gamakumara, P, Athanasopoulos, G, & Hyndman, RJ. (2023). Probabilistic forecast reconciliation: Properties, evaluation and score optimisation. European J Operational Research, 306(2), 693–706.
Wickramasuriya, SL, Athanasopoulos, G, & Hyndman, RJ. (2019). Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. Journal of the American Statistical Association, 114(526), 804–819.
Actually, some recent nonlinear reconciliation methods require a slightly more complicated equation. This equation is for general linear reconciliation methods.↩︎
This “unbiasedness preserving” constraint was first introduced in Hyndman et al. (2011). Panagiotelis et al. (2021) show that this is equivalent to \bm{S}\bm{G} being a projection matrix onto the m-dimensional coherent subspace for which the aggregation constraints hold.↩︎
Note that k_{h} is a proportionality constant. It does not need to be estimated or specified here as it gets cancelled out in Equation 11.8.↩︎
