Hierarchical forecasting with unequal (short and long) length of time series.

[ggjx22]

Hello sktime Community,

I am currently in the midst of exploring the usage and end-to-end set up of Hierarchical Forecasting and I am met with the following obstacles which I will like to seek guidance here.

This is what I want to achieve:

• To have an end-to-end pipeline that implement Hierarchical Forecasting using a single and tuned global ensemble tree model.
• To forecast each and every time series at lowest node in the hierarchy by a certain forecast horizon regardless the length of the time series.
• To reconcile forecast at each level of the hierarchy.

The nature of my dataset:

• Hierarchy structure: Group companies -> country -> branch -> department -> product.
• forecast target: transactional values of product (all positives and negatives and mix of both).
• Project requirements are to forecast on product level and later on aggregate by sum up to different levels of the hierarchy and measure percentage error between forecast and actual.

Obstacles:

1. Some time series do not have at least 24 observations. They are as short as 12 observations.

Questions about my problem statement:

1. If my dataset suffer from lack of data from some time series, does that mean this approach is not suitable till I have enough data?
2. Aggregating target to a node higher can give me a collective time series with longer length, but are there still other approaches to my problem?

Questions about the code below:

1. If a Hierarchical dataset contains time series of different seasonality, does TransformedTargetForecaster applies the same sp parameter in Deseasonalizer to all time series?
2. Does the usage of PluginParamsTransformer solve this? And then again, there is still a need to fit the training data and manually inspect the value of sp via .get_params()["sp"].

Below is the reproducible example:

# imports
import pandas as pd
import numpy as np
from sktime.utils._testing.hierarchical import _make_hierarchical
from sktime.forecasting.base import ForecastingHorizon
from sklearn.ensemble import RandomForestRegressor
from sktime.pipeline import make_pipeline
from sktime.forecasting.compose import make_reduction, ForecastingPipeline, TransformedTargetForecaster
from sktime.transformations.series.detrend import Detrender, Deseasonalizer
from sktime.transformations.series.difference import Differencer
from sktime.transformations.series.summarize import WindowSummarizer
from sktime.transformations.series.date import DateTimeFeatures
from sktime.transformations.compose import TransformIf
from sktime.param_est.stationarity import StationarityADF
from sklearn.preprocessing import OrdinalEncoder
from sktime.transformations.series.adapt import PandasTransformAdaptor
from sktime.param_est.compose import ParamFitterPipeline
from sktime.param_est.seasonality import SeasonalityACF

RANDOM_STATE = 2

# create hierarchical dataset
data = _make_hierarchical(
    hierarchy_levels=(3, 4),
    max_timepoints=100,
    min_timepoints=12,
    all_positive=False,
    index_type="period",
    random_state=RANDOM_STATE,
)

# simulate exogenous variables out from hierarchy names, assuming they are useful
data['feature_1'] = data.index.get_level_values('h0')
data['feature_2'] = data.index.get_level_values('h1')

# create train and holdout dataset with last 12 observations
test = pd.DataFrame(columns=['c0'])

for name, group in data.groupby(level=['h0', 'h1']):
    test = pd.concat([test, group.tail(12)])
    
test.index = pd.MultiIndex.from_tuples(test.index, names=('h0', 'h1', 'time'))
train = data.drop(test.index)

# split X and y
X_train = train.drop(columns=['c0'])
y_train = train[['c0']]
X_test = test.drop(columns=['c0'])
y_test = test[['c0']]

# doing this because check raise TypeError: input should not have column of 'object' dtype
group_encoder = OrdinalEncoder()
group_encoder.fit(X_train)
train_encoding = group_encoder.transform(X_train)
X_train['feature_1_enc'] = train_encoding[:, 0]
X_train['feature_2_enc'] = train_encoding[:, 1]
X_train.drop(columns=['feature_1', 'feature_2'], inplace=True)

test_encoding = group_encoder.transform(X_test)
X_test['feature_1_enc'] = test_encoding[:, 0]
X_test['feature_2_enc'] = test_encoding[:, 1]
X_test.drop(columns=['feature_1', 'feature_2'], inplace=True)

# define forecast horizon
fh = ForecastingHorizon(np.arange(1, 12+1, 1), is_relative=True)

# create transformation & forecasting pipeline for target
kwargs = {
    "lag_feature": {
        "lag": [1],
        "mean": [[1, 3]],
        "std": [[1, 3]],
    }
}

tree_regressor = make_pipeline(RandomForestRegressor(random_state=RANDOM_STATE))

tree_forecaster = make_reduction(
    estimator=tree_regressor,
    window_length=None,
    transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
    pooling='global'
)

detrender = Detrender()
nan_remover = PandasTransformAdaptor(method='dropna')
deseasonalizer = Deseasonalizer(sp=12)

auto_difference = TransformIf(
    if_estimator=StationarityADF(),
    param='stationary',
    condition='==',
    condition_value=False,
    then_trafo=Differencer(),
    else_trafo=None
)

forecast_pipe = TransformedTargetForecaster(
    steps=[
        ('difference_if', auto_difference),
        ('remove_trend', detrender), # creates nan values
        ('remove_nan', nan_remover),
        ('remove_season', deseasonalizer),
        ('forecaster', tree_forecaster)
    ]
)

# define transformers to create exogenous features
date_features = DateTimeFeatures(ts_freq='M', keep_original_columns=True)

# create end to end pipeline
pipe = ForecastingPipeline(
    steps=[
        ('date_features', date_features),
        ('forecaster', forecast_pipe)
    ]
)

pipe.fit(y=y_train, X=X_train, fh=fh) # errors here

error:
ValueError: x must have 2 complete cycles requires 24 observations. x only has 15 observation(s)

package versions:

System:
    python: 3.10.9 | packaged by Anaconda, Inc. | (main, Mar  1 2023, 18:18:15) [MSC v.1916 64 bit (AMD64)]
executable: [c:\Users\agpgago\data-science\forecast-wizard-test\venv\Scripts\python.exe](file:///C:/Users/agpgago/data-science/forecast-wizard-test/venv/Scripts/python.exe)
   machine: Windows-10-10.0.19042-SP0

Python dependencies:
          pip: 23.3.2
       sktime: 0.26.1
      sklearn: 1.4.1.post1
       skbase: 0.7.2
        numpy: 1.26.4
        scipy: 1.12.0
       pandas: 2.1.4
   matplotlib: 3.8.3
       joblib: 1.3.2
        numba: 0.58.1
  statsmodels: 0.14.1
     pmdarima: 2.0.4
statsforecast: 1.7.3
      tsfresh: 0.20.2
      tslearn: 0.6.3
        torch: None
   tensorflow: 2.15.0
tensorflow_probability: None

[fkiraly]

1. If my dataset suffer from lack of data from some time series, does that mean this approach is not suitable till I have enough data?

That depends on the forecaster. If you use panel forecasts via reduction to tabular regressor, e.g., gradient boosting, you do not necessarily need long individual time series. That depends on the capabilities of the forecaster. Having said that, many off-shelf implementations and classical approaches might assume equal length and minimum length, these would not be suitable. You may like to try make_reduction and DirectReductionForecaster (feedback appreciated).

Aggregating target to a node higher can give me a collective time series with longer length, but are there still other approaches to my problem?

Instead of aggregating, you could try reconciliation. I am not sure how well ReconcilerForecaster behaves if some time series are short or unequal length - FYI @ciaran-g.

If a Hierarchical dataset contains time series of different seasonality, does TransformedTargetForecaster applies the same sp parameter in Deseasonalizer to all time series?

Yes - you've set it to 12 explicitly as far as I can see, so it will always be 12 periods.

I'm not sure where your error comes from - details on the traceback would be helpful. Does it come from Deseasonalizer? It is plausible that that needs to complete periods.

Does the usage of PluginParamsTransformer solve this? And then again, there is still a need to fit the training data and manually inspect the value of sp via .get_params()["sp"].

If you use the PluginParamsTransformer, then indeed the seasonality should be automatically determined per series. You do not necessarily need to manually inspect, it will set it automatically. It should also not detect a seasonality that is longer than half the data length, so that issue would also be solved (need to check though if that's true).

[ggjx22]

Hello @fkiraly thank you for your replies and guidance, I very much appreciate it. But there are more doubts to clarify, which led to more confusion on my end. Sorry for the long reply.

Here are some feedbacks to your replies.

That depends on the forecaster. If you use panel forecasts via reduction to tabular regressor, e.g., gradient boosting, you do not necessarily need long individual time series. That depends on the capabilities of the forecaster. Having said that, many off-shelf implementations and classical approaches might assume equal length and minimum length, these would not be suitable.

Yes, I realise most statistical (classical) models require at least 24 months of observations. Not sure about unequal length because I only intend to use ensemble tree models to learn dependency relationships between time series. I think statistical models does not have those capability, unless I'm wrong (happy to be proved wrong).

You may like to try make_reduction and DirectReductionForecaster (feedback appreciated).

Unlike make_reduction which accepts window_length=None, DirectReductionForecaster seems to require this argument to be given. I would expect window_length to be inferred from WindowSummarizer, as explained in the tutorial, but that does not seem to be the case.

tree_regressor = make_pipeline(RandomForestRegressor(random_state=RANDOM_STATE))

# tree_forecaster = make_reduction(
#     estimator=tree_regressor,
#     window_length=None,
#     transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
#     pooling='global'
# )

tree_forecaster = DirectReductionForecaster(
    estimator=tree_regressor,
    window_length=None,
    transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
    pooling='global'
)

TypeError                                 Traceback (most recent call last)
Cell In[37], line 19
     10 tree_regressor = make_pipeline(RandomForestRegressor(random_state=RANDOM_STATE))
     12 # tree_forecaster = make_reduction(
     13 #     estimator=tree_regressor,
     14 #     window_length=None,
     15 #     transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
     16 #     pooling='global'
     17 # )
---> 19 tree_forecaster = DirectReductionForecaster(
     20     estimator=tree_regressor,
     21     window_length=None,
     22     transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
     23     pooling='global'
     24 )
     26 detrender = Detrender()
     27 nan_remover = PandasTransformAdaptor(method='dropna')

File c:\Users\agpgago\data-science\forecast-wizard-test\venv\lib\site-packages\sktime\forecasting\compose\_reduce.py:1906, in DirectReductionForecaster.__init__(self, estimator, window_length, transformers, X_treatment, impute_method, pooling)
   1904 self.impute_method = impute_method
   1905 self.pooling = pooling
-> 1906 self._lags = list(range(window_length))
   1907 super().__init__()
   1909 warn(
...
   1911     "user feedback is appreciated in issue #3224 here: "
   1912     "https://github.com/sktime/sktime/issues/3224"
   1913 )

TypeError: 'NoneType' object cannot be interpreted as an integer

Instead of aggregating, you could try reconciliation.

I meant to say - by aggregating the time series to a node higher and then fitting them to the model can help me solve my issue of having short time series. But doing this does not fulfil my project requirement.

Anyway, hope to have more feedback about the usability of ReconcilerForecaster .

I'm not sure where your error comes from - details on the traceback would be helpful. Does it come from Deseasonalizer? It is plausible that that needs to complete periods.

Yes, that's right. It comes from Deseaonalizer. And that is because of Deseasonalizer(sp=12). This is resolve through the implementation of PluginParamsTransformer which plugs in the estimated sp parameter into Deseasonalizer() for time series which does not follows to sp=12.

sp_estimator = Differencer() * SeasonalityACF()  
deseason = Deseasonalizer()
deseason_w_sp_est = PluginParamsTransformer(sp_estimator, deseason)

If you use the PluginParamsTransformer, then indeed the seasonality should be automatically determined per series. You do not necessarily need to manually inspect, it will set it automatically. It should also not detect a seasonality that is longer than half the data length, so that issue would also be solved (need to check though if that's true).

The issue is resolved. Fitting the end to end pipeline has no problem. There is another issue though. During pipe.predict(X=X_test), PandasTransformAdaptor does not have inverse_transform.

The reason for it to exist in the first place is because I thought I have seen NaN values after applying Detrender(). Anyhow, I rechecked y_train again after applying Detrender() and the NaNs were gone! I removed it from TransformedTargetForecaster, fit the end to end pipeline again, another issue arise. I'm not sure if this is a genuine bug or wrong usage or the case of data quality.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[101], line 9
      1 # create end to end pipeline
      2 pipe = ForecastingPipeline(
      3     steps=[
      4         ('date_features', date_features),
      5         ('forecaster', forecast_pipe)
      6     ]
      7 )
----> 9 pipe.fit(y=y_train, X=X_train, fh=fh)

File c:\Users\agpgago\data-science\forecast-wizard-test\venv\lib\site-packages\sktime\forecasting\base\_base.py:387, in BaseForecaster.fit(self, y, X, fh)
    385 # we call the ordinary _fit if no looping/vectorization needed
    386 if not vectorization_needed:
--> 387     self._fit(y=y_inner, X=X_inner, fh=fh)
    388 else:
    389     # otherwise we call the vectorized version of fit
    390     self._vectorize("fit", y=y_inner, X=X_inner, fh=fh)

File c:\Users\agpgago\data-science\forecast-wizard-test\venv\lib\site-packages\sktime\forecasting\compose\_pipeline.py:509, in ForecastingPipeline._fit(self, y, X, fh)
    507 name, forecaster = self.steps_[-1]
    508 f = forecaster.clone()
--> 509 f.fit(y=y, X=X, fh=fh)
    510 self.steps_[-1] = (name, f)
...
--> 153     raise ValueError("This function does not handle missing values")
    154 if model.startswith("m"):
    155     if np.any(x <= 0):

ValueError: This function does not handle missing values

So, I went to check on my target transformer pipeline if they were fitted per series in y_train. And there were no missing values:

for name, group in y_train.groupby(level=['h0', 'h1']):
    print(name)
    print(group.shape)
    transfo_pipe = auto_difference * detrender * deseason_w_sp_est
    _t = transfo_pipe.fit_transform(group)
    print('Any NaN values: ', _t.isna().sum())
    print('='*50)

('h0_0', 'h1_0')
(40, 1)
Any NaN values:  c0    0
dtype: int64
==================================================
('h0_0', 'h1_1')
(15, 1)
Any NaN values:  c0    0
dtype: int64
==================================================
('h0_1', 'h1_0')
(72, 1)
Any NaN values:  c0    0
dtype: int64
==================================================
('h0_1', 'h1_1')
(22, 1)
Any NaN values:  c0    0
dtype: int64
==================================================

And X_train looks like this:

for name, group in X_train.groupby(level=['h0', 'h1']):
    print(name)
    print(group.head())
    print('='*50)

('h0_0', 'h1_0')
                   feature_1_enc  feature_2_enc
h0   h1   time                                 
h0_0 h1_0 2004-01            0.0            0.0
          2004-02            0.0            0.0
          2004-03            0.0            0.0
          2004-04            0.0            0.0
          2004-05            0.0            0.0
==================================================
('h0_0', 'h1_1')
                   feature_1_enc  feature_2_enc
h0   h1   time                                 
h0_0 h1_1 2006-02            0.0            1.0
          2006-03            0.0            1.0
          2006-04            0.0            1.0
          2006-05            0.0            1.0
          2006-06            0.0            1.0
==================================================
('h0_1', 'h1_0')
                   feature_1_enc  feature_2_enc
h0   h1   time                                 
h0_1 h1_0 2001-05            1.0            0.0
          2001-06            1.0            0.0
          2001-07            1.0            0.0
          2001-08            1.0            0.0
          2001-09            1.0            0.0
==================================================
('h0_1', 'h1_1')
                   feature_1_enc  feature_2_enc
h0   h1   time                                 
h0_1 h1_1 2005-07            1.0            1.0
          2005-08            1.0            1.0
          2005-09            1.0            1.0
          2005-10            1.0            1.0
          2005-11            1.0            1.0
==================================================

I am more confused than I was before. The thing is if PandasTransformAdaptor(method='dropna') exist in the target transformer pipeline, then pipe.fit(y=y_train, X=X_train, fh=fh) fits with no errors. But the end to end pipeline fails at pipe.predict(X=X_test). The forecaster is still using make_reduction and not DirectReductionForecaster

tree_regressor = make_pipeline(RandomForestRegressor(random_state=RANDOM_STATE))

tree_forecaster = make_reduction(
    estimator=tree_regressor,
    window_length=None,
    transformers=[WindowSummarizer(**kwargs, n_jobs=-1, truncate='bfill')],
    pooling='global'

detrender = Detrender()
nan_remover = PandasTransformAdaptor(method='dropna')
sp_estimator = Differencer() * SeasonalityACF()  
deseason = Deseasonalizer()
deseason_w_sp_est = PluginParamsTransformer(sp_estimator, deseason)

auto_difference = TransformIf(
    if_estimator=StationarityADF(),
    param='stationary',
    condition='==',
    condition_value=False,
    then_trafo=Differencer(),
    else_trafo=None
)

forecast_pipe = TransformedTargetForecaster(
    steps=[
        ('difference_if', auto_difference),
        ('remove_trend', detrender),
        # ('remove_nan', nan_remover),
        ('remove_season', deseason_w_sp_est),
        ('forecaster', tree_forecaster)
    ]
)

[fkiraly]

Replying to explicit and implicit questions above.

I'm opening a couple threads as there are multiple diferent comments in the above.

[fkiraly]

I think statistical models does not have those capability, unless I'm wrong (happy to be proved wrong).

If you mean, the ability to deal with unequal length series in a panel or hierarchical data set: this is subtle, but I do think that models like ARIMA can do that. The reason is that these classical models are not panel models but deal with single time series. Any such model is broadcast by sktime to individual time series ("for loop" or parallel), and there it does not matter whether they are unequal length, because ARIMA sees only a single time series.

It is also a limitation that there is no tag for forecasters that tells the user whether it can deal with unequal length time series. This could be added, but I wonder what it should be set to, in the above situation. Probably True?

[ggjx22]

It may be beneficial to include tags that tells user if a model is suitable for panel and local/global forecasting task? E.g. {'panel_hierarchical_task': 'local'} (For ARIMA case). For tag that deals with unequal length time series, it can be as what you suggested.

[fkiraly]

yes, good idea: opened issue here #6070

[fkiraly]

Unlike make_reduction which accepts window_length=None, DirectReductionForecaster seems to require this argument to be given. I would expect window_length to be inferred from WindowSummarizer, as explained in the tutorial, but that does not seem to be the case.

Design-wise, it seems not right to have an estimator reach into its components and get parameters or attributes from there.

Do you mean, if no window_length is given to DirectReductionForecaster, but one of its components is WindowSummarizer, it should default to identical window_length of that component? This is at least not what I would expect.

Could you explain why this is your expectation?

[ggjx22]

No you're right, window_length and lags windown of WindowSummarizer should not be related.

[fkiraly]

I meant to say - by aggregating the time series to a node higher and then fitting them to the model can help me solve my issue of having short time series. But doing this does not fulfil my project requirement.

Anyway, hope to have more feedback about the usability of ReconcilerForecaster.

Any feedback?

[ggjx22]

I tried with dummy data but no avail. I will create an issue.

[ggjx22]

#6068 created.

[fkiraly]

The issue is resolved. Fitting the end to end pipeline has no problem. There is another issue though. During pipe.predict(X=X_test), PandasTransformAdaptor does not have inverse_transform.

Can you explain how the lack of inverse_transform becomes a problem? We could add that, the question is how to expose this to the user.

What I see is that something seems to upset a pandas function.
This may be a bug, in any case to say more we'd need a complete code snippet, perhaps you can reproduce the issue with dummy/random data?

Another option could be using DropNA or Impute, these transformers are dedicated to handling nans, and may be "safer" in that respect, of course I can't promise these will solve the nan problem.

[ggjx22]

For some reason, after applying Detrender(), lots of nans appear. I have removed it from the pipeline and replaced it with Differencer(). I can create a reproducible example and open an issue.

[ggjx22]

What I have experienced is described in the Detrender() documentation.

Depending on time indices, this can generate in-sample or out-of-sample residuals.

I'm not sure why is this happening. Maybe I won't create an issue for this. Below is the reproducible example.

from sktime.utils._testing.hierarchical import _make_hierarchical
from sktime.transformations.series.detrend import Detrender
from sktime.split import temporal_train_test_split

RANDOM_STATE = 2

# create hierarchical dataset
data = _make_hierarchical(
    hierarchy_levels=(2, 2),
    max_timepoints=100,
    all_positive=False,
    index_type="period",
    random_state=RANDOM_STATE,
)

# create hold out set
train, test = temporal_train_test_split(y=data, test_size=12)

# split X and y
X_train = train
y_train = train[['c0']]
X_test = test
y_test = test[['c0']]

y_train_t = Detrender().fit_transform(y_train)

y_train_t.loc[('h0_0', 'h1_0')]

[fkiraly]

yes, separate issue would be appreciated. My guess is that it is maybe sth in combination with the differencer, that creates nans at the start.