Large collections of time series organized into structures at different aggregation levels often require their forecasts to follow their aggregation constraints, 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 classic hierarchical reconciliation.

In this notebook we will show how to use the StatsForecast library to produce base forecasts, and use HierarchicalForecast package to perform hierarchical reconciliation.

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

​
1. Libraries

# !pip install hierarchicalforecast
# !pip install -U numba statsforecast datasetsforecast

​
2. Load Data

In this example we will use the TourismSmall dataset. The following cell gets the time series for the different levels in the hierarchy, the summing matrix S 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, HierarchicalInfo

group_name = 'TourismSmall'
group = HierarchicalInfo.get_group(group_name)
Y_df, S_df, tags = HierarchicalData.load('./data', group_name)
Y_df['ds'] = pd.to_datetime(Y_df['ds'])

Y_df.head()
unique_iddsy
0total1998-03-3184503
1total1998-06-3065312
2total1998-09-3072753
3total1998-12-3170880
4total1999-03-3186893
S_df.iloc[:5, :5]
nsw-hol-citynsw-hol-noncityvic-hol-cityvic-hol-noncityqld-hol-city
total1.01.01.01.01.0
hol1.01.01.01.01.0
vfr0.00.00.00.00.0
bus0.00.00.00.00.0
oth0.00.00.00.00.0
tags

{'Country': array(['total'], dtype=object),
 'Country/Purpose': array(['hol', 'vfr', 'bus', 'oth'], dtype=object),
 'Country/Purpose/State': array(['nsw-hol', 'vic-hol', 'qld-hol', 'sa-hol', 'wa-hol', 'tas-hol',
        'nt-hol', 'nsw-vfr', 'vic-vfr', 'qld-vfr', 'sa-vfr', 'wa-vfr',
        'tas-vfr', 'nt-vfr', 'nsw-bus', 'vic-bus', 'qld-bus', 'sa-bus',
        'wa-bus', 'tas-bus', 'nt-bus', 'nsw-oth', 'vic-oth', 'qld-oth',
        'sa-oth', 'wa-oth', 'tas-oth', 'nt-oth'], dtype=object),
 'Country/Purpose/State/CityNonCity': array(['nsw-hol-city', 'nsw-hol-noncity', 'vic-hol-city',
        'vic-hol-noncity', 'qld-hol-city', 'qld-hol-noncity',
        'sa-hol-city', 'sa-hol-noncity', 'wa-hol-city', 'wa-hol-noncity',
        'tas-hol-city', 'tas-hol-noncity', 'nt-hol-city', 'nt-hol-noncity',
        'nsw-vfr-city', 'nsw-vfr-noncity', 'vic-vfr-city',
        'vic-vfr-noncity', 'qld-vfr-city', 'qld-vfr-noncity',
        'sa-vfr-city', 'sa-vfr-noncity', 'wa-vfr-city', 'wa-vfr-noncity',
        'tas-vfr-city', 'tas-vfr-noncity', 'nt-vfr-city', 'nt-vfr-noncity',
        'nsw-bus-city', 'nsw-bus-noncity', 'vic-bus-city',
        'vic-bus-noncity', 'qld-bus-city', 'qld-bus-noncity',
        'sa-bus-city', 'sa-bus-noncity', 'wa-bus-city', 'wa-bus-noncity',
        'tas-bus-city', 'tas-bus-noncity', 'nt-bus-city', 'nt-bus-noncity',
        'nsw-oth-city', 'nsw-oth-noncity', 'vic-oth-city',
        'vic-oth-noncity', 'qld-oth-city', 'qld-oth-noncity',
        'sa-oth-city', 'sa-oth-noncity', 'wa-oth-city', 'wa-oth-noncity',
        'tas-oth-city', 'tas-oth-noncity', 'nt-oth-city', 'nt-oth-noncity'],
       dtype=object)}

We split the dataframe in train/test splits.

Y_test_df = Y_df.groupby('unique_id').tail(group.horizon)
Y_train_df = Y_df.drop(Y_test_df.index)

​
3. Base forecasts

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

from statsforecast.core import StatsForecast
from statsforecast.models import AutoARIMA, Naive

fcst = StatsForecast(
    df=Y_train_df, 
    models=[AutoARIMA(season_length=group.seasonality), Naive()], 
    freq=group.freq, 
    n_jobs=-1
)
Y_hat_df = fcst.forecast(h=group.horizon)

​
4. Hierarchical reconciliation

The following cell makes the previous forecasts coherent using the HierarchicalReconciliation class. The used methods to make the forecasts coherent are:

  • BottomUp: The reconciliation of the method is a simple addition to the upper levels.
  • TopDown: The second method constrains the base-level predictions to the top-most aggregate-level serie and then distributes it to the disaggregate series through the use of proportions.
  • MiddleOut: Anchors the base predictions in a middle level.
from hierarchicalforecast.core import HierarchicalReconciliation
from hierarchicalforecast.methods import BottomUp, TopDown, MiddleOut

reconcilers = [
    BottomUp(),
    TopDown(method='forecast_proportions'),
    MiddleOut(middle_level='Country/Purpose/State', 
              top_down_method='forecast_proportions')
]
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)

​
5. 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.set_index('unique_id'),
        tags=tags, benchmark='Naive'
)
evaluation.filter(like='ARIMA', axis=1).T
levelOverallCountryCountry/PurposeCountry/Purpose/StateCountry/Purpose/State/CityNonCity
metricmse-scaledmse-scaledmse-scaledmse-scaledmse-scaled
AutoARIMA0.1689570.1328360.1938020.1823510.199594
AutoARIMA/BottomUp0.1001880.085180.0770970.1557710.199594
AutoARIMA/TopDown_method-forecast_proportions0.1442390.1328360.119760.2044260.227395
AutoARIMA/MiddleOut_middle_level-Country/Purpose/State_top_down_method-forecast_proportions0.1342190.1457760.0912760.1823510.215141

​
References

InstallIntroduction