Large collections of time series organized into structures at different aggregation levels often require their forecasts to follow their aggregation constraints and to be nonnegative, which poses the challenge of creating novel algorithms capable of coherent forecasts.

The HierarchicalForecast package provides a wide collection of Python implementations of hierarchical forecasting algorithms that follow nonnegative hierarchical reconciliation.

In this notebook, we will show how to use the HierarchicalForecast package to perform nonnegative reconciliation of forecasts on Wiki2 dataset.

You can run these experiments using CPU or GPU with Google Colab.

Open In Colab
!pip install hierarchicalforecast statsforecast datasetsforecast

1. Load Data

In this example we will use the Wiki2 dataset. The following cell gets the time series for the different levels in the hierarchy, the summing dataframe S_df which recovers the full dataset from the bottom level hierarchy and the indices of each hierarchy denoted by tags.

import numpy as np
import pandas as pd

from datasetsforecast.hierarchical import HierarchicalData
Y_df, S_df, tags = HierarchicalData.load('./data', 'Wiki2')
Y_df['ds'] = pd.to_datetime(Y_df['ds'])
Y_df.head()
unique_iddsy
0Total2016-01-01156508
1Total2016-01-02129902
2Total2016-01-03138203
3Total2016-01-04115017
4Total2016-01-05126042
S_df.iloc[:5, :5]
de_AAC_AAG_001de_AAC_AAG_010de_AAC_AAG_014de_AAC_AAG_045de_AAC_AAG_063
Total11111
de11111
en00000
fr00000
ja00000
tags
{'Views': array(['Total'], dtype=object),
 'Views/Country': array(['de', 'en', 'fr', 'ja', 'ru', 'zh'], dtype=object),
 'Views/Country/Access': array(['de_AAC', 'de_DES', 'de_MOB', 'en_AAC', 'en_DES', 'en_MOB',
        'fr_AAC', 'fr_DES', 'fr_MOB', 'ja_AAC', 'ja_DES', 'ja_MOB',
        'ru_AAC', 'ru_DES', 'ru_MOB', 'zh_AAC', 'zh_DES', 'zh_MOB'],
       dtype=object),
 'Views/Country/Access/Agent': array(['de_AAC_AAG', 'de_AAC_SPD', 'de_DES_AAG', 'de_MOB_AAG',
        'en_AAC_AAG', 'en_AAC_SPD', 'en_DES_AAG', 'en_MOB_AAG',
        'fr_AAC_AAG', 'fr_AAC_SPD', 'fr_DES_AAG', 'fr_MOB_AAG',
        'ja_AAC_AAG', 'ja_AAC_SPD', 'ja_DES_AAG', 'ja_MOB_AAG',
        'ru_AAC_AAG', 'ru_AAC_SPD', 'ru_DES_AAG', 'ru_MOB_AAG',
        'zh_AAC_AAG', 'zh_AAC_SPD', 'zh_DES_AAG', 'zh_MOB_AAG'],
       dtype=object),
 'Views/Country/Access/Agent/Topic': array(['de_AAC_AAG_001', 'de_AAC_AAG_010', 'de_AAC_AAG_014',
        'de_AAC_AAG_045', 'de_AAC_AAG_063', 'de_AAC_AAG_100',
        'de_AAC_AAG_110', 'de_AAC_AAG_123', 'de_AAC_AAG_143',
        'de_AAC_SPD_012', 'de_AAC_SPD_074', 'de_AAC_SPD_080',
        'de_AAC_SPD_105', 'de_AAC_SPD_115', 'de_AAC_SPD_133',
        'de_DES_AAG_064', 'de_DES_AAG_116', 'de_DES_AAG_131',
        'de_MOB_AAG_015', 'de_MOB_AAG_020', 'de_MOB_AAG_032',
        'de_MOB_AAG_059', 'de_MOB_AAG_062', 'de_MOB_AAG_088',
        'de_MOB_AAG_095', 'de_MOB_AAG_109', 'de_MOB_AAG_122',
        'de_MOB_AAG_149', 'en_AAC_AAG_044', 'en_AAC_AAG_049',
        'en_AAC_AAG_075', 'en_AAC_AAG_114', 'en_AAC_AAG_119',
        'en_AAC_AAG_141', 'en_AAC_SPD_004', 'en_AAC_SPD_011',
        'en_AAC_SPD_026', 'en_AAC_SPD_048', 'en_AAC_SPD_067',
        'en_AAC_SPD_126', 'en_AAC_SPD_140', 'en_DES_AAG_016',
        'en_DES_AAG_024', 'en_DES_AAG_042', 'en_DES_AAG_069',
        'en_DES_AAG_082', 'en_DES_AAG_102', 'en_MOB_AAG_018',
        'en_MOB_AAG_022', 'en_MOB_AAG_101', 'en_MOB_AAG_124',
        'fr_AAC_AAG_029', 'fr_AAC_AAG_046', 'fr_AAC_AAG_070',
        'fr_AAC_AAG_087', 'fr_AAC_AAG_098', 'fr_AAC_AAG_104',
        'fr_AAC_AAG_111', 'fr_AAC_AAG_112', 'fr_AAC_AAG_142',
        'fr_AAC_SPD_025', 'fr_AAC_SPD_027', 'fr_AAC_SPD_035',
        'fr_AAC_SPD_077', 'fr_AAC_SPD_084', 'fr_AAC_SPD_097',
        'fr_AAC_SPD_130', 'fr_DES_AAG_023', 'fr_DES_AAG_043',
        'fr_DES_AAG_051', 'fr_DES_AAG_058', 'fr_DES_AAG_061',
        'fr_DES_AAG_091', 'fr_DES_AAG_093', 'fr_DES_AAG_094',
        'fr_DES_AAG_136', 'fr_MOB_AAG_006', 'fr_MOB_AAG_030',
        'fr_MOB_AAG_066', 'fr_MOB_AAG_117', 'fr_MOB_AAG_120',
        'fr_MOB_AAG_121', 'fr_MOB_AAG_135', 'fr_MOB_AAG_147',
        'ja_AAC_AAG_038', 'ja_AAC_AAG_047', 'ja_AAC_AAG_055',
        'ja_AAC_AAG_076', 'ja_AAC_AAG_099', 'ja_AAC_AAG_128',
        'ja_AAC_AAG_132', 'ja_AAC_AAG_134', 'ja_AAC_AAG_137',
        'ja_AAC_SPD_013', 'ja_AAC_SPD_034', 'ja_AAC_SPD_050',
        'ja_AAC_SPD_060', 'ja_AAC_SPD_078', 'ja_AAC_SPD_106',
        'ja_DES_AAG_079', 'ja_DES_AAG_081', 'ja_DES_AAG_113',
        'ja_MOB_AAG_065', 'ja_MOB_AAG_073', 'ja_MOB_AAG_092',
        'ja_MOB_AAG_127', 'ja_MOB_AAG_129', 'ja_MOB_AAG_144',
        'ru_AAC_AAG_008', 'ru_AAC_AAG_145', 'ru_AAC_AAG_146',
        'ru_AAC_SPD_000', 'ru_AAC_SPD_090', 'ru_AAC_SPD_148',
        'ru_DES_AAG_003', 'ru_DES_AAG_007', 'ru_DES_AAG_017',
        'ru_DES_AAG_041', 'ru_DES_AAG_071', 'ru_DES_AAG_072',
        'ru_MOB_AAG_002', 'ru_MOB_AAG_040', 'ru_MOB_AAG_083',
        'ru_MOB_AAG_086', 'ru_MOB_AAG_103', 'ru_MOB_AAG_107',
        'ru_MOB_AAG_118', 'ru_MOB_AAG_125', 'zh_AAC_AAG_021',
        'zh_AAC_AAG_033', 'zh_AAC_AAG_037', 'zh_AAC_AAG_052',
        'zh_AAC_AAG_057', 'zh_AAC_AAG_085', 'zh_AAC_AAG_108',
        'zh_AAC_SPD_039', 'zh_AAC_SPD_096', 'zh_DES_AAG_009',
        'zh_DES_AAG_019', 'zh_DES_AAG_053', 'zh_DES_AAG_054',
        'zh_DES_AAG_056', 'zh_DES_AAG_068', 'zh_DES_AAG_089',
        'zh_DES_AAG_139', 'zh_MOB_AAG_005', 'zh_MOB_AAG_028',
        'zh_MOB_AAG_031', 'zh_MOB_AAG_036', 'zh_MOB_AAG_138'], dtype=object)}

We split the dataframe in train/test splits.

Y_test_df = Y_df.groupby('unique_id').tail(7)
Y_train_df = Y_df.drop(Y_test_df.index)
Y_test_df = Y_test_df.set_index('unique_id')
Y_train_df = Y_train_df.set_index('unique_id')

2. Base Forecasts

The following cell computes the base forecast for each time series using the ETS and naive models. Observe that Y_hat_df contains the forecasts but they are not coherent.

from statsforecast.models import ETS, Naive
from statsforecast.core import StatsForecast
fcst = StatsForecast(
    df=Y_train_df, 
    models=[ETS(season_length=7, model='ZAA'), Naive()], 
    freq='D', 
    n_jobs=-1
)
Y_hat_df = fcst.forecast(h=7)

Observe that the ETS model computes negative forecasts for some series.

Y_hat_df.query('ETS < 0')
dsETSNaive
unique_id
de_AAC_AAG_0012016-12-25-487.601532340.0
de_AAC_AAG_0012016-12-26-215.634201340.0
de_AAC_AAG_0012016-12-27-173.175613340.0
de_AAC_AAG_0012016-12-30-290.836060340.0
de_AAC_AAG_0012016-12-31-784.441040340.0
…………
zh_AAC_AAG_0332016-12-31-86.52642137.0
zh_MOB2016-12-26-199.5348821036.0
zh_MOB2016-12-27-69.5272601036.0
zh_MOB_AAG2016-12-26-199.5348821036.0
zh_MOB_AAG2016-12-27-69.5272601036.0

3. Non-Negative Reconciliation

The following cell makes the previous forecasts coherent and nonnegative using the HierarchicalReconciliation class.

from hierarchicalforecast.methods import MinTrace
from hierarchicalforecast.core import HierarchicalReconciliation
reconcilers = [
    MinTrace(method='ols'),
    MinTrace(method='ols', nonnegative=True)
]
hrec = HierarchicalReconciliation(reconcilers=reconcilers)
Y_rec_df = hrec.reconcile(Y_hat_df=Y_hat_df, Y_df=Y_train_df,
                          S=S_df, tags=tags)

Observe that the nonnegative reconciliation method obtains nonnegative forecasts.

Y_rec_df.query('`ETS/MinTrace_method-ols_nonnegative-True` < 0')
dsETSNaiveETS/MinTrace_method-olsNaive/MinTrace_method-olsETS/MinTrace_method-ols_nonnegative-TrueNaive/MinTrace_method-ols_nonnegative-True
unique_id

The free reconciliation method gets negative forecasts.

Y_rec_df.query('`ETS/MinTrace_method-ols` < 0')
dsETSNaiveETS/MinTrace_method-olsNaive/MinTrace_method-olsETS/MinTrace_method-ols_nonnegative-TrueNaive/MinTrace_method-ols_nonnegative-True
unique_id
de_DES2016-12-25-2553.932861495.0-3468.745214495.02.262540e-15495.0
de_DES2016-12-26-2155.228271495.0-2985.587125495.01.356705e-30495.0
de_DES2016-12-27-2720.993896495.0-3698.680055495.06.857413e-30495.0
de_DES2016-12-29-3429.432617495.0-2965.207609495.02.456449e+02495.0
de_DES2016-12-30-3963.202637495.0-3217.360371495.03.646790e+02495.0
……………………
zh_MOB_AAG_0362016-12-2675.298317115.0-165.799776115.03.207772e-14115.0
zh_MOB_AAG_0362016-12-2772.895554115.0-134.340626115.02.308198e-14115.0
zh_MOB_AAG_1382016-12-2594.79662365.0-47.00981365.03.116938e-1465.0
zh_MOB_AAG_1382016-12-2671.29398365.0-169.80411065.00.000000e+0065.0
zh_MOB_AAG_1382016-12-2762.04974465.0-145.18643665.00.000000e+0065.0

4. Evaluation

The HierarchicalForecast package includes the HierarchicalEvaluation class to evaluate the different hierarchies and also is capable of compute scaled metrics compared to a benchmark model.

from hierarchicalforecast.evaluation import HierarchicalEvaluation
def mse(y, y_hat):
    return np.mean((y-y_hat)**2)

evaluator = HierarchicalEvaluation(evaluators=[mse])
evaluation = evaluator.evaluate(
        Y_hat_df=Y_rec_df, Y_test_df=Y_test_df, 
        tags=tags, benchmark='Naive'
)
evaluation.filter(like='ETS', axis=1).T
levelOverallViewsViews/CountryViews/Country/AccessViews/Country/Access/AgentViews/Country/Access/Agent/Topic
metricmse-scaledmse-scaledmse-scaledmse-scaledmse-scaledmse-scaled
ETS1.0115850.73581.1903541.1036571.0895151.397139
ETS/MinTrace_method-ols0.9791630.6983551.0625211.1432771.1133491.354041
ETS/MinTrace_method-ols_nonnegative-True0.9450750.6778921.0046391.1847191.1414421.158672

Observe that the nonnegative reconciliation method performs better that its unconstrained counterpart.

References