Is it possible to use GridSearchCV for hierarchical data with Reconciler?

[anthonygiorgio97]

Hi,

I have a hierarchical time series dataframe and for each time series I want to find the best forecaster from a list running a GridSearchCV using a params grid space. Then I want to reconcile the result with the ReconcilerForecaster class.

If I insert the reconciler forecaster into the forecaster parameter of ForecastingGridSearchCV I have a message that this is not supported.

Another way, I think, is to fit the ForecastingGridSearchCV for each series separately. But then how can I assign the best forecaster I found for each time series to the ReconcilerForecaster ?

[fkiraly]

Strange, this should work - in theory (and my expectation), you should be able to combine the above-mentioned wrappers and any sktime forecaster arbitrarily.

Metrics, grid search, reconciler, individual forecasters, all should work out-of-the-box with hierarchical data.

There was a bug with ReconcilerForecaster which @ciaran-g fixed in the latest release 0.15.1.

In case something breaks for you, we would appreciate some short and self-contained code with dummy data that has the error you are experiencing, so we can debug.

FYI some people who come to my mind as having worked recently on hierarchical functionality: @ciaran-g, @danbartl, @KishManani.

[fkiraly]

Regarding your second question:

Another way, I think, is to fit the ForecastingGridSearchCV for each series separately. But then how can I assign the best forecaster I found for each time series to the ReconcilerForecaster ?

If you want to do a grid search for each series separately, you need to wrap the grid search in ForecastByLevel, otherwise the grid search will compute the aggregate score and will try to find a single parameter setting that is best for all series, together (in mean aggregation if you use metrics with default values).

If you have the ReconcilerForecaster inside the grid search, the errors will be computed after reconciliation, if it is outside, the error to optimize will be before reconciling.

I'm not sure whether there is currently a way to both (a) fit grid search parameter by series and (b) tune by using metrics that are computed after reconciliation.

It will be possible though once @VyomkeshVyas finishes the HierarchyEnsembleForecaster here: #3905
Then, you would use the HierarchyEnsembleForecaster to construct an ensemble where parameters can be accessed by node, wrap it in the ReconcilerForecaster (or pipeline with the Reconciler transformer), then wrap in grid search.

[anthonygiorgio97]

Here there is an example of my use case:

from sktime.utils._testing.hierarchical import _make_hierarchical
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.trend import PolynomialTrendForecaster
from sktime.forecasting.model_selection import ForecastingGridSearchCV, ExpandingWindowSplitter
from sktime.transformations.hierarchical.aggregate import Aggregator
from sktime.forecasting.reconcile import ReconcilerForecaster
from sktime.forecasting.compose import TransformedTargetForecaster

y = _make_hierarchical()
agg = Aggregator()
y_agg = agg.fit_transform(y)

param_grid = [{"forecaster": [ExponentialSmoothing()],
               "forecaster__trend": ['add','mul']
              },
              {"forecaster": [PolynomialTrendForecaster()],
               "forecaster__degree": [1,2]}
              ]

pipe = TransformedTargetForecaster(steps=[
    ("forecaster", ExponentialSmoothing())])

N_cv_fold = 2
step_cv = 1
fh = [1,2]

initial_window_cv_len = len(y_agg.index.get_level_values(2).unique()) - (N_cv_fold - 1) * step_cv - fh[-1]

cv = ExpandingWindowSplitter(
            initial_window = initial_window_cv_len,
            step_length = step_cv,
            fh = fh)

reconciler = ReconcilerForecaster(pipe, method="ols")

gscv = ForecastingGridSearchCV(
        forecaster=reconciler,
        param_grid=param_grid,
        cv=cv,
        n_jobs=-1,
        verbose = 1
        )

gscv.fit(y_agg)

However in this case, as you said, the grid search will compute the aggregate score and will try to find a single parameter setting that is best for all series, together and not the best forecaster for each series in the hierarchy (which is my goal)

Another approach would be to use a for cicle for the most granular hierarchy index and store best forecaster for each series in a dictionary (however it will not use the benefits of vectorization) as the following code. It would then be ideal if we could configure the forecasters for each series, for example in the dataframe of the image below extract from a tutorial on sktime site:

param_grid = [{"forecaster": [ExponentialSmoothing()],
               "forecaster__trend": ['add','mul']
              },
              {"forecaster": [PolynomialTrendForecaster()],
               "forecaster__degree": [1,2]}
              ]

N_cv_fold = 2
step_cv = 1
fh = [1,2]

initial_window_cv_len = len(y_agg.index.get_level_values(2).unique()) - (N_cv_fold - 1) * step_cv - fh[-1]

hierarchy = list(set(list(zip(y_agg.index.get_level_values(0),
                              y_agg.index.get_level_values(1)))))

forecaster_dict = {}

for ts in hierarchy:

    y_ts = y_agg[(y_agg.index.get_level_values(0) == ts[0]) &
                 (y_agg.index.get_level_values(1) == ts[1])]

    pipe = TransformedTargetForecaster(steps=[
                    ("forecaster", ExponentialSmoothing())])

    cv = ExpandingWindowSplitter(
            initial_window = initial_window_cv_len,
            step_length = step_cv,
            fh = fh)
    
    gscv = ForecastingGridSearchCV(
                forecaster=pipe,
                param_grid=param_grid,
                cv=cv,
                n_jobs=-1,
                verbose = 1
                )

    gscv.fit(y_ts)
    forecaster_dict[ts] = gscv.best_forecaster_

[anthonygiorgio97]

It is odd that it appears only 4 times, and I wouldn't be able to say why the number 4 is special here. Are you sure you just get 4 printouts? If yes, it might be worth investigating and/or reporting as a bug.

I checked again. By restarting the kernel I do not have this issue again.

---

So this would be my final solution for my problem:

from sktime.utils._testing.hierarchical import _make_hierarchical
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.trend import PolynomialTrendForecaster
from sktime.forecasting.model_selection import ForecastingGridSearchCV, ExpandingWindowSplitter
from sktime.transformations.hierarchical.aggregate import Aggregator
from sktime.forecasting.reconcile import ReconcilerForecaster
from sktime.forecasting.compose import TransformedTargetForecaster
from sktime.forecasting.compose import ForecastByLevel

y = _make_hierarchical()
agg = Aggregator()
y_agg = agg.fit_transform(y) 

param_grid = [{"forecaster": [ExponentialSmoothing()],
               "forecaster__trend": ['add','mul']
              },
              {"forecaster": [PolynomialTrendForecaster()],
               "forecaster__degree": [1,2]}
              ]

pipe = TransformedTargetForecaster(steps=[
    ("forecaster", ExponentialSmoothing())])

N_cv_fold = 2
step_cv = 1
fh = [1,2]

initial_window_cv_len = len(y_agg.index.get_level_values(2).unique()) - (N_cv_fold - 1) * step_cv - fh[-1]

cv = ExpandingWindowSplitter(
            initial_window = initial_window_cv_len,
            step_length = step_cv,
            fh = fh)

gscv = ForecastingGridSearchCV(
        forecaster=pipe,
        param_grid=param_grid,
        cv=cv,
        n_jobs=-1,
        verbose = 1
        )

gscv_bylevel = ForecastByLevel(gscv, 'local')

reconciler = ReconcilerForecaster(gscv_bylevel, method="ols")

reconciler.fit(y_agg)

Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits
Fitting 2 folds for each of 4 candidates, totalling 8 fits

---

To access the best fitted forecaster for each time series hierarchy I can run the following code:

best_forecasters = {}
for ts in reconciler.get_fitted_params()['forecaster__ForecastByLevel__forecasters'].index:
    best_forecasters[ts] = reconciler.get_fitted_params()['forecaster__ForecastByLevel__forecasters'] \
                                       .loc[(ts[0], ts[1]),'forecasters'] \
                                       .get_fitted_params()["forecaster__best_forecaster__forecaster"]
best_forecasters

{('__total', '__total'): ExponentialSmoothing(trend='add'),
 ('h0_0', '__total'): PolynomialTrendForecaster(),
 ('h0_0', 'h1_0'): ExponentialSmoothing(trend='mul'),
 ('h0_0', 'h1_1'): ExponentialSmoothing(trend='mul'),
 ('h0_0', 'h1_2'): ExponentialSmoothing(trend='add'),
 ('h0_0', 'h1_3'): PolynomialTrendForecaster(),
 ('h0_1', '__total'): PolynomialTrendForecaster(),
 ('h0_1', 'h1_0'): ExponentialSmoothing(trend='mul'),
 ('h0_1', 'h1_1'): PolynomialTrendForecaster(),
 ('h0_1', 'h1_2'): ExponentialSmoothing(trend='add'),
 ('h0_1', 'h1_3'): ExponentialSmoothing(trend='add')}

[jpss95]

This solution seemed to be suitable for my problem, where I wanted to take a local forecasting approach by finding the best model for each time series. When I tried it on a small scale, it worked perfectly, but when attempting to use it on a problem with many time series, it didn't scale.

I'm using a machine with multiple cores, and I'm providing the argument n_jobs=-1 to ForecastingGridSearchCV to make use of them. However, in reality, it doesn't utilize them, and all the work is done sequentially.

cc: @fkiraly

[fkiraly]

However, in reality, it doesn't utilize them, and all the work is done sequentially.

ForecastingGridSearchCV parallelizes parameters, not time series instances in the hierarchy - we could write parallelization for the multiple series though.

[fkiraly]

or is this a more general issue with the tuners? Could you perhaps try the same without hierarchical data, and check your cores?

[jpss95]

Thanks for the clarification, Franz. It makes sense that ForecastingGridSearchCV focuses on hyperparameter optimization. But I believe that introducing a parallelization feature at the time series level could greatly enhance the efficiency of the library for users like me who deal with multiple time series.

[fkiraly]

Hm, what you did here is very interesting:

best_forecasters = {}
for ts in gscv_bylevel.get_fitted_params()["forecasters"].index:
    best_forecasters[ts] = gscv_bylevel.get_fitted_params()["forecasters"] \
                                       .loc[(ts[0], ts[1]),'forecasters'] \
                                       .get_fitted_params()["forecaster__best_forecaster__forecaster"]
best_forecasters

that is, the data frame is filled with nested parameters if called with the name of the nested parameter.

Perhaps this should just happen by default if you call gscv_bylevel.get_fitted_params()["best_forecaster"] or similar, don't you think?

[anthonygiorgio97]

If I try this I get an error. However it will be very helpful if implemented

[fkiraly]

yes, I know. That's what I'm saying, it should just produce the data frame you made! Currently it doesn't.

[fkiraly]

i.e., I got the idea that it would be very helpful just to happen by default, in general, from looking at what you did. Not claiming that it currently behaves that way.

[fkiraly]

here is a prototype, @anthonygiorgio97:
#4107

would appreciate some feedback what you think!

[fkiraly]

It would then be ideal if we could configure the forecasters for each series, for example in the dataframe of the image below extract from a tutorial on sktime site:

Can you explain what you mean here by "it would be ideal"?

• does the ForecastByLevel approach work for you? Or is that not what you want?
• or, do you think it's too "clunky", and you would prefer a parameter, say, in ForecastingGridSearchCV?

[fkiraly]

oh, I mentioned it here in my first reply, it might have gotten lost in all the threads: #4101 (comment)

[fkiraly]

confused - you also used it here in your own code?? #4101 (reply in thread)

[anthonygiorgio97]

Yes. Sorry, my mistake: I have posted my first reply and after I documented and tried the ForecastByLevel approach

[fkiraly]

ah, I see, so "was not aware" refers to your state of knowledge 4 hours ago. Makes sense. Would appreciate feedback in case you decide to try it out.

[fkiraly]

might also be good if we add some of the newer stuff to the hierarchical forecasting tutorial.

[daniel-torresc]

Hello,

I have tried to reproduce the code on this thread (see below) and had an issue when trying to access the reconciler.forecaster.get_fitted_params() method of the reconciler when fitted. The output is the following:

NotFittedError                            Traceback (most recent call last)
Cell In[17], line 1
----> 1 reconciler.forecaster.get_fitted_params()

File ~\AppData\Roaming\Python\Python39\site-packages\sktime\forecasting\base\_base.py:1256, in BaseForecaster.get_fitted_params(self, deep)
1254 # if self is not vectorized, run the default get_fitted_params
1255 if not getattr(self, "_is_vectorized", False):
-> 1256     return super(BaseForecaster, self).get_fitted_params(deep=deep)
1258 # otherwise, we delegate to the instances' get_fitted_params
1259 # instances' parameters are returned at dataframe-slice-like keys
1260 fitted_params = {}

File ~\AppData\Roaming\Python\Python39\site-packages\sktime\base\_base.py:417, in BaseEstimator.get_fitted_params(self, deep)
    386 """Get fitted parameters.
    387 
    388 State required:
(...)
    414       e.g., `[componentname]__[componentcomponentname]__[paramname]`, etc
    415 """
    416 if not self.is_fitted:
--> 417     raise NotFittedError(
    418         f"estimator of type {type(self).__name__} has not been "
    419         "fitted yet, please call fit on data before get_fitted_params"
    420     )
    422 # collect non-nested fitted params of self
    423 fitted_params = self._get_fitted_params()

NotFittedError: estimator of type ForecastByLevel has not been fitted yet, please call fit on data before get_fitted_params

I have come up with a "partial solution" which is fitting again the forecaster, in addition to the previous reconciler fit. This is:

Step 1: Calling reconciler.fit(y_train_agg) (here when trying to access reconciler.forecaster.get_fitted_params() it throws the exception)
Step 2: Additionally, execute reconciler.forecaster.fit(y_train_agg) (ok when accessing get_fitted_params() after executing this fit)

I have tried this with sktime 0.17.1, 0.17.2 and the latest 0.18 versions and all had the same issue.

Is there any way to access get_fitted_params() without needing to explicitly fit the forecaster (step 2) again? I would like to access it once the reconciler is fitted.

Has anyone been able to execute this code without having errors? I understand the solution to my question is what is proposed in the thread, but I am having the issue listed above.

The code I have tried is the one below:

from sktime.utils._testing.hierarchical import _make_hierarchical
from sktime.forecasting.exp_smoothing import ExponentialSmoothing
from sktime.forecasting.trend import PolynomialTrendForecaster
from sktime.forecasting.model_selection import ForecastingGridSearchCV, ExpandingWindowSplitter
from sktime.transformations.hierarchical.aggregate import Aggregator
from sktime.forecasting.reconcile import ReconcilerForecaster
from sktime.forecasting.compose import TransformedTargetForecaster
from sktime.forecasting.compose import ForecastByLevel

y = _make_hierarchical()
agg = Aggregator()
y_agg = agg.fit_transform(y) 

param_grid = [{"forecaster": [ExponentialSmoothing()],
               "forecaster__trend": ['add','mul']
              },
              {"forecaster": [PolynomialTrendForecaster()],
               "forecaster__degree": [1,2]}
              ]

pipe = TransformedTargetForecaster(steps=[
    ("forecaster", ExponentialSmoothing())])

N_cv_fold = 2
step_cv = 1
fh = [1,2]

initial_window_cv_len = len(y_agg.index.get_level_values(2).unique()) - (N_cv_fold - 1) * step_cv - fh[-1]

cv = ExpandingWindowSplitter(
            initial_window = initial_window_cv_len,
            step_length = step_cv,
            fh = fh)

gscv = ForecastingGridSearchCV(
        forecaster=pipe,
        param_grid=param_grid,
        cv=cv,
        n_jobs=-1,
        verbose = 1
        )

gscv_bylevel = ForecastByLevel(gscv, 'local')

reconciler = ReconcilerForecaster(gscv_bylevel, method="ols")

reconciler.fit(y_agg)

# reconciler.forecaster.fit(y_agg)  # If this line is uncommented, the error disappears

reconciler.forecaster.get_fitted_params()  # here we get the error when accessing ".get_fitted_params()"

Thank you.

[fkiraly]

apologies, just saw you pasted it. Let me try.

[fkiraly]

ok, here's the thing: that reconciler.forecaster.get_fitted_params fails is not a bug, because the forecaster is just a template, a clone is fitted.

However, if you call reconciler.get_fitted_params, it fails - that is a bug!

[fkiraly]

reported here: #4574

[daniel-torresc]

Thank you so much for the quick response! I will keep an eye on it and will check again when it's fixed.

[fkiraly]

fixed in 0.18.1, @daniel-torresc