Reconciliation

This example provides a minimal example of how to reconcile a set of existing forecasts.

This example is based on chapter 11 of:

Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3. Accessed on 05/07/2022.

The data comes from the online book of Hyndman, see link below.

The code for this example can be found at:

Jupyter Notebook: https://github.com/elephaint/hierts/tree/main/docs/notebooks
Python script: https://github.com/elephaint/hierts/tree/main/examples

Import packages & read data

[1]:

#%% Read packages
import pandas as pd
import numpy as np
from hierts.reconciliation import calc_summing_matrix, apply_reconciliation_methods,\
                                  calc_level_method_rmse

[2]:

#%% Read data
df = pd.read_csv("https://OTexts.com/fpp3/extrafiles/prison_population.csv")

[3]:

# Let's look at the data
df.head(10)

[3]:

	Date	State	Gender	Legal	Indigenous	Count
0	2005-03-01	ACT	Female	Remanded	ATSI	0
1	2005-03-01	ACT	Female	Remanded	Non-ATSI	2
2	2005-03-01	ACT	Female	Sentenced	ATSI	0
3	2005-03-01	ACT	Female	Sentenced	Non-ATSI	5
4	2005-03-01	ACT	Male	Remanded	ATSI	7
5	2005-03-01	ACT	Male	Remanded	Non-ATSI	58
6	2005-03-01	ACT	Male	Sentenced	ATSI	5
7	2005-03-01	ACT	Male	Sentenced	Non-ATSI	101
8	2005-03-01	NSW	Female	Remanded	ATSI	51
9	2005-03-01	NSW	Female	Remanded	Non-ATSI	131

As you can see, we have a column indicating the time index (‘Date’), four columns that provide information on the hierarchy/aggregation of the data (‘State’, ‘Gender’, ‘Legal’, ‘Indigenous’), and a target column (‘Count’).

Set aggregations and calculate summing matrix

For this dataset, the aggregations are in the columns: [‘State’, ‘Gender’, ‘Legal’, ‘Indigenous’]. These columns will be in our aggregations.

[4]:

aggregation_cols = ['State', 'Gender', 'Legal', 'Indigenous']

Next, we define the aggregations that we are interested in. In this case, we are interested in the following aggregations:

[5]:

aggregations = [['State'],
                ['State', 'Gender'],
                ['State', 'Legal'],
                ['State', 'Indigenous'],
                ['Gender', 'Legal']]

Don’t include the top (total) level and bottom-level: these will be added automatically later on. Now, we can create a summing matrix, which shows how each of the bottom-level series maps to our chosen aggregations:

[6]:

df_S = calc_summing_matrix(df, aggregation_cols, aggregations)

Let’s have a look at df_S:

[7]:

df_S

[7]:

	bottom_timeseries	ACT-Female-Remanded-ATSI	ACT-Female-Remanded-Non-ATSI	ACT-Female-Sentenced-ATSI	ACT-Female-Sentenced-Non-ATSI	ACT-Male-Remanded-ATSI	ACT-Male-Remanded-Non-ATSI	ACT-Male-Sentenced-ATSI	ACT-Male-Sentenced-Non-ATSI	NSW-Female-Remanded-ATSI	NSW-Female-Remanded-Non-ATSI	...	VIC-Male-Sentenced-ATSI	VIC-Male-Sentenced-Non-ATSI	WA-Female-Remanded-ATSI	WA-Female-Remanded-Non-ATSI	WA-Female-Sentenced-ATSI	WA-Female-Sentenced-Non-ATSI	WA-Male-Remanded-ATSI	WA-Male-Remanded-Non-ATSI	WA-Male-Sentenced-ATSI	WA-Male-Sentenced-Non-ATSI
Aggregation	Value
Total	Total	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	...	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0
State	ACT	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	NSW	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	1.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	NT	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	QLD	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
bottom_timeseries	WA-Female-Sentenced-Non-ATSI	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0	0.0
	WA-Male-Remanded-ATSI	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0	0.0
	WA-Male-Remanded-Non-ATSI	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0	0.0
	WA-Male-Sentenced-ATSI	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0	0.0
	WA-Male-Sentenced-Non-ATSI	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	...	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0

125 rows × 64 columns

As you can see, df_S is a mapping between the bottom-level series (columns of df_S) and the aggregations specified. Some observations:

‘Total’ contains the total across all bottom-level series. Hence, this is a row consisting of ones - the total is the sum of all the bottom-level series.
‘Bottom level’ contains the bottom level series. This is an identity matrix, as each bottom-level series in the rows only maps to a single bottom level series in the columns
Everything in between ‘Bottom level’ and ‘Total’ denotes how we can construct the aggregations.

Create random forecasts for each time series

Now we can create a forecasting model for each time series in the aggregation matrix df_S. In this minimal example, we’ll just use random forecasts.

[8]:

# Set target, time_index and split of train and test.
target = 'Count'
time_index = 'Date'
end_train = '2015-12-31'
start_test = '2016-01-01'

[9]:

# Create random forecasts - in this case we just do random poisson sampling of the actual values
rng = np.random.default_rng(seed=0)
df[f'{target}_predicted'] = rng.poisson(lam=df[f'{target}'])

[10]:

# Add bottom_timeseries identifier and create actuals & forecasts dataframe for all aggregations
df['bottom_timeseries'] = df[aggregation_cols].agg('-'.join, axis=1)
actuals_bottom_timeseries = df.set_index(['bottom_timeseries', time_index])[target]\
                              .unstack(1)\
                              .loc[df_S.columns]
forecasts_bottom_timeseries = df.set_index(['bottom_timeseries', time_index])[f'{target}_predicted']\
                                .unstack(1)\
                                .loc[df_S.columns]
actuals = df_S @ actuals_bottom_timeseries
forecasts = df_S @ forecasts_bottom_timeseries

Reconciliation

We can now reconcile the forecasts.

[11]:

# Create forecast test set and apply a set of reconciliation methods.
forecasts_test = forecasts.loc[:, start_test:]
forecasts_reconciled = apply_reconciliation_methods(forecasts_test, df_S, \
                                                    actuals.loc[:, :end_train],\
                                                    forecasts.loc[:, :end_train],\
                                                    methods=['ols', 'wls_var', 'mint_shrink'])

Method ols, reconciliation time: 0.0024s
Method wls_var, reconciliation time: 0.0006s
Method mint_shrink, reconciliation time: 3.7461s

Finally, we can compute the root mean-squared error on all of our methods.

[12]:

# Calculate error for all levels and methods. We set bottom-up as the base method to compare against
# in the relative rmse
rmse, rel_rmse = calc_level_method_rmse(forecasts_reconciled, actuals, base='base')

Let’s look at the errors:

[13]:

rmse

[13]:

Method	base	ols	wls_var	mint_shrink
Aggregation
Total	200.393987	200.392359	200.392720	200.393987
Gender-Legal	71.493444	71.493115	71.493200	71.493444
State	70.355570	70.355378	70.355157	70.355570
State-Indigenous	50.534487	50.534313	50.534269	50.534487
State-Legal	48.516911	48.516884	48.516749	48.516911
State-Gender	47.778624	47.778658	47.778282	47.778624
bottom_timeseries	24.704899	24.704902	24.704829	24.704899
All series	45.105233	45.105116	45.105014	45.105233

In this case, the reconciliation methods don’t do very much. This example is just to illustrate what constitutes a minimal working example to produce reconciled forecasts on a given set of existing forecasts.