PERMBU/Normality probabilistic reconciliation fail with `level=0` prediction interval

[melopeo]

Reconciliation fails when intervals_method="permbu"

import pandas as pd

# compute base forecast no coherent
from statsforecast.core import StatsForecast
from statsforecast.models import AutoARIMA, Naive

# obtain hierarchical reconciliation methods and evaluation
from hierarchicalforecast.core import HierarchicalReconciliation
from hierarchicalforecast.methods import BottomUp, TopDown, MiddleOut

# obtain hierarchical datasets
from datasetsforecast.hierarchical import HierarchicalData

# Load TourismSmall dataset
Y_df, S, tags = HierarchicalData.load("./data", "TourismSmall")
Y_df["ds"] = pd.to_datetime(Y_df["ds"])

level = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90]

# Compute base level predictions
sf = StatsForecast(df=Y_df, models=[AutoARIMA(season_length=12)], freq="M", n_jobs=1)

forecasts_df = sf.forecast(h=12, fitted=True, level=level)

# Reconcile the base predictions
reconcilers = [BottomUp()]

hrec = HierarchicalReconciliation(reconcilers=reconcilers)

reconciled_forecasts = hrec.reconcile(
    Y_hat_df=forecasts_df,
    Y_df=sf.forecast_fitted_values(),
    S=S,
    tags=tags,
    level=level,
    intervals_method="permbu",
)

Thrown error:

/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/core.py:158: RuntimeWarning: invalid value encountered in divide
  sigmah = sign * (sigmah - y_hat_model) / z
/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/sklearn/preprocessing/_encoders.py:808: FutureWarning: `sparse` was renamed to `sparse_output` in version 1.2 and will be removed in 1.4. `sparse_output` is ignored unless you leave `sparse` to its default value.
  warnings.warn(
Traceback (most recent call last):
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3442, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-c88d133fa83d>", line 1, in <module>
    reconciled_forecasts = hrec.reconcile(
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/core.py", line 197, in reconcile
    fcsts_model = reconcile_fn(y_hat=y_hat_model, **kwargs)
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/methods.py", line 92, in reconcile
    return _reconcile(S, P, W, y_hat, level=level,
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/methods.py", line 48, in _reconcile
    res = sampler.get_prediction_levels(res=res, level=level)
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/probabilistic_methods.py", line 340, in get_prediction_levels
    samples = self.get_samples(y_hat=res['mean'])
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/hierarchicalforecast/probabilistic_methods.py", line 324, in get_samples
    parent_samples = np.einsum('ab,bhs->ahs', Agg, children_samples)
  File "<__array_function__ internals>", line 180, in einsum
  File "/Users/XXX/anaconda3/envs/YYY/lib/python3.8/site-packages/numpy/core/einsumfunc.py", line 1371, in einsum
    return c_einsum(*operands, **kwargs)
ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (36,44)->(36,newaxis,newaxis,44) (18,12,35)->(12,35,18) 
/Applications/PyCharm with Anaconda plugin .app/Contents/plugins/python/helpers/pydev/_pydevd_bundle/pydevd_utils.py:462: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.
  for item in s.iteritems():

[kdgutier]

Hi @melopeo,

From the logs, is it the case that the 0 level breaks the following division?

sigmah = sign * (sigmah - y_hat_model) / z

Would you be able to try level 0.001?
In particular, I believe that combining the Normality assumption on PERMBU with level=0 translates into an unbounded prediction.

If this is the case, we will add a Raise.Exception protection to the levels domain.

[kdgutier]

Yes just tested the case level = [10, 20, 30, 40, 50, 60, 70, 80, 90]
And it works, I added the domain protection the level list to live in (0,100].

Thanks for reporting the bug.

[melopeo]

Thank you!