Repository: sktime/sktime
Pull request: [ENH][BUG] clean up `Detrender` and implement correct multivariate behaviour
This PR improves the `Detrender`:
* it fixes https://github.com/sktime/sktime/issues/4051 and also simplifies the logic in the `Detrender`. As the `BaseForecasting` interface supports loop over hierarchy instances and variables natively, all the case distinctions internally can go, and with it goes the incorrect multivariate implementation.
* adds support for detrending with hierarchical forecasters
* adds support for detrending with forecasters that require `fh` in `fit` already (this was not possible previously)

As a consequence, the test `test_polynomial_detrending` had to be updated since the internal mtype of the `Detrender` changed, from `pd.Series` based to `pd.DataFrame` based (requiring a minor change in how objects are being compared).

Problem statement:
[BUG] incorrect multivariate detrending in `Detrender`
Currently, the `Detrender` transformer implements multivariate detrending incorrectly, and in contrast to what its docstring says.

The difference between "should" and "is" happens if the forecaster passed is genuinely multivariate - in that case, the `Detrender` coerces it to univariate, and fits an univariate version per variable, instead of just applying the multivariate forecaster and detrending via it.

This has the status of a bug since the docstring does not agree with content, and the content is contrary to user expectation.

Hints from issue discussion:

Commit message: update test

Patch to review:
diff --git a/sktime/transformations/series/detrend/_detrend.py b/sktime/transformations/series/detrend/_detrend.py
index 1fd3f2fbcb0..3179f687d23 100644
--- a/sktime/transformations/series/detrend/_detrend.py
+++ b/sktime/transformations/series/detrend/_detrend.py
@@ -4,10 +4,11 @@
 """Implements transformations to detrend a time series."""
 
 __all__ = ["Detrender"]
-__author__ = ["mloning", "SveaMeyer13", "KishManani"]
+__author__ = ["mloning", "SveaMeyer13", "KishManani", "fkiraly"]
 
 import pandas as pd
 
+from sktime.datatypes import update_data
 from sktime.forecasting.base._fh import ForecastingHorizon
 from sktime.forecasting.trend import PolynomialTrendForecaster
 from sktime.transformations.base import BaseTransformer
@@ -73,10 +74,11 @@ class Detrender(BaseTransformer):
         "scitype:transform-output": "Series",
         # what scitype is returned: Primitives, Series, Panel
         "scitype:instancewise": True,  # is this an instance-wise transform?
-        "X_inner_mtype": ["pd.DataFrame", "pd.Series"],
+        "X_inner_mtype": ["pd.DataFrame", "pd-multiindex", "pd_multiindex_hier"],
         # which mtypes do _fit/_predict support for X?
-        "y_inner_mtype": "pd.DataFrame",  # which mtypes do _fit/_predict support for y?
-        "univariate-only": True,
+        "y_inner_mtype": ["pd.DataFrame", "pd-multiindex", "pd_multiindex_hier"],
+        # which mtypes do _fit/_predict support for y?
+        "univariate-only": False,
         "fit_is_empty": False,
         "capability:inverse_transform": True,
         "transform-returns-same-time-index": True,
@@ -85,14 +87,18 @@ class Detrender(BaseTransformer):
     def __init__(self, forecaster=None, model="additive"):
         self.forecaster = forecaster
         self.model = model
-        self.forecaster_ = None
+
         super(Detrender, self).__init__()
 
-        # whether this transformer is univariate depends on the forecaster
-        #  this transformer is univariate iff the forecaster is univariate
-        if forecaster is not None:
-            fc_univ = forecaster.get_tag("scitype:y", "univariate") == "univariate"
-            self.set_tags(**{"univariate-only": fc_univ})
+        # default for forecaster - written to forecaster_ to not overwrite param
+        if self.forecaster is None:
+            self.forecaster_ = PolynomialTrendForecaster(degree=1)
+        else:
+            self.forecaster_ = forecaster.clone()
+
+        allowed_models = ("additive", "multiplicative")
+        if model not in allowed_models:
+            raise ValueError("`model` must be 'additive' or 'multiplicative'")
 
     def _fit(self, X, y=None):
         """Fit transformer to X and y.
@@ -110,29 +116,31 @@ def _fit(self, X, y=None):
         -------
         self: a fitted instance of the estimator
         """
-        if self.forecaster is None:
-            self.forecaster = PolynomialTrendForecaster(degree=1)
-
-        # univariate: X is pd.Series
-        if isinstance(X, pd.Series):
-            forecaster = self.forecaster.clone()
-            # note: the y in the transformer is exogeneous in the forecaster, i.e., X
-            self.forecaster_ = forecaster.fit(y=X, X=y)
-        # multivariate
-        elif isinstance(X, pd.DataFrame):
-            self.forecaster_ = {}
-            for colname in X.columns:
-                forecaster = self.forecaster.clone()
-                self.forecaster_[colname] = forecaster.fit(y=X[colname], X=y)
+        if not self.forecaster_.get_tag("requires-fh-in-fit", True):
+            self.forecaster_.fit(y=X, X=y)
         else:
-            raise TypeError("X must be pd.Series or pd.DataFrame")
-
-        allowed_models = ("additive", "multiplicative")
-        if self.model not in allowed_models:
-            raise ValueError("`model` must be 'additive' or 'multiplicative'")
-
+            self._X = X
+            self._y = y
         return self
 
+    def _get_fh_from_X(self, X):
+        """Obtain fh from X, which can be simple or hierarchical."""
+        if not isinstance(X.index, pd.MultiIndex):
+            time_index = X.index
+        else:
+            time_index = X.index.get_level_values(-1).unique()
+        return ForecastingHorizon(time_index, is_relative=False)
+
+    def _get_fitted_forecaster(self, X, y, fh):
+        """Obtain fitted forecaster from self."""
+        if self.forecaster_.get_tag("requires-fh-in-fit", True):
+            X = update_data(self._X, X)
+            y = update_data(self._y, y)
+            forecaster = self.forecaster_.clone().fit(y=X, X=y, fh=fh)
+        else:
+            forecaster = self.forecaster_
+        return forecaster
+
     def _transform(self, X, y=None):
         """Transform X and return a transformed version.
 
@@ -150,37 +158,15 @@ def _transform(self, X, y=None):
         Xt : pd.Series or pd.DataFrame, same type as X
             transformed version of X, detrended series
         """
-        fh = ForecastingHorizon(X.index, is_relative=False)
-
-        # univariate: X is pd.Series
-        if isinstance(X, pd.Series):
-            # note: the y in the transformer is exogeneous in the forecaster, i.e., X
-            X_pred = self.forecaster_.predict(fh=fh, X=y)
-            if self.model == "additive":
-                return X - X_pred
-            elif self.model == "multiplicative":
-                return X / X_pred
-        # multivariate: X is pd.DataFrame
-        elif isinstance(X, pd.DataFrame):
-            Xt = X.copy()
-            # check if all columns are known
-            X_fit_keys = set(self.forecaster_.keys())
-            X_new_keys = set(X.columns)
-            difference = X_new_keys.difference(X_fit_keys)
-            if len(difference) != 0:
-                raise ValueError(
-                    "X contains columns that have not been "
-                    "seen in fit: " + str(difference)
-                )
-            for colname in Xt.columns:
-                X_pred = self.forecaster_[colname].predict(fh=fh, X=y)
-                if self.model == "additive":
-                    Xt[colname] = Xt[colname] - X_pred
-                elif self.model == "multiplicative":
-                    Xt[colname] = Xt[colname] / X_pred
-            return Xt
-        else:
-            raise TypeError("X must be pd.Series or pd.DataFrame")
+        fh = self._get_fh_from_X(X=X)
+        forecaster = self._get_fitted_forecaster(X=X, y=y, fh=fh)
+
+        X_pred = forecaster.predict(fh=fh, X=y)
+
+        if self.model == "additive":
+            return X - X_pred
+        elif self.model == "multiplicative":
+            return X / X_pred
 
     def _inverse_transform(self, X, y=None):
         """Logic used by `inverse_transform` to reverse transformation on `X`.
@@ -197,36 +183,15 @@ def _inverse_transform(self, X, y=None):
         Xt : pd.Series or pd.DataFrame, same type as X
             inverse transformed version of X
         """
-        fh = ForecastingHorizon(X.index, is_relative=False)
-
-        # univariate: X is pd.Series
-        if isinstance(X, pd.Series):
-            # note: the y in the transformer is exogeneous in the forecaster, i.e., X
-            X_pred = self.forecaster_.predict(fh=fh, X=y)
-            if self.model == "additive":
-                return X + X_pred
-            elif self.model == "multiplicative":
-                return X * X_pred
-        # multivariate: X is pd.DataFrame
-        if isinstance(X, pd.DataFrame):
-            X = X.copy()
-            # check if all columns are known
-            X_fit_keys = set(self.forecaster_.keys())
-            X_new_keys = set(X.columns)
-            difference = X_new_keys.difference(X_fit_keys)
-            if len(difference) != 0:
-                raise ValueError(
-                    "X contains columns that have not been "
-                    "seen in fit: " + str(difference)
-                )
-            for colname in X.columns:
-                X_pred = self.forecaster_[colname].predict(fh=fh, X=y)
-                if self.model == "additive":
-                    X[colname] = X[colname] + X_pred
-                elif self.model == "multiplicative":
-                    X[colname] = X[colname] * X_pred
-
-            return X
+        fh = self._get_fh_from_X(X=X)
+        forecaster = self._get_fitted_forecaster(X=X, y=y, fh=fh)
+
+        X_pred = forecaster.predict(fh=fh, X=y)
+
+        if self.model == "additive":
+            return X + X_pred
+        elif self.model == "multiplicative":
+            return X * X_pred
 
     def _update(self, X, y=None, update_params=True):
         """Update the parameters of the detrending estimator with new data.
@@ -247,24 +212,11 @@ def _update(self, X, y=None, update_params=True):
         -------
         self : an instance of self
         """
-        # multivariate
-        if isinstance(X, pd.DataFrame):
-            # check if all columns are known
-            X_fit_keys = set(self.forecaster_.keys())
-            X_new_keys = set(X.columns)
-            difference = X_new_keys.difference(X_fit_keys)
-            if len(difference) != 0:
-                raise ValueError(
-                    "Z contains columns that have not been "
-                    "seen in fit: " + str(difference)
-                )
-            for colname in X.columns:
-                self.forecaster_[colname].update(
-                    y=X[colname], X=y, update_params=update_params
-                )
-        # univariate
-        else:
+        if not self.forecaster_.get_tag("requires-fh-in-fit", True):
             self.forecaster_.update(y=X, X=y, update_params=update_params)
+        else:
+            self._X = update_data(self._X, X)
+            self._y = update_data(self._y, y)
         return self
 
     @classmethod
diff --git a/sktime/transformations/series/detrend/tests/test_detrend.py b/sktime/transformations/series/detrend/tests/test_detrend.py
index 50ff9e676ae..970f8242d56 100644
--- a/sktime/transformations/series/detrend/tests/test_detrend.py
+++ b/sktime/transformations/series/detrend/tests/test_detrend.py
@@ -24,7 +24,7 @@ def y_dataframe():
 
 
 def test_polynomial_detrending():
-
+    """Test that transformer results agree with manual detrending."""
     y = pd.Series(np.arange(20) * 0.5) + np.random.normal(0, 1, size=20)
     forecaster = PolynomialTrendForecaster(degree=1, with_intercept=True)
     transformer = Detrender(forecaster)
@@ -38,9 +38,11 @@ def test_polynomial_detrending():
     np.testing.assert_array_almost_equal(actual_coefs, expected_coefs)
 
     # check trend
-    expected_trend = expected_coefs[0] + np.arange(len(y)) * expected_coefs[1]
-    actual_trend = transformer.forecaster_.predict(-np.arange(len(y)))
-    np.testing.assert_array_almost_equal(actual_trend, expected_trend)
+    n = len(y)
+    expected_trend = expected_coefs[0] + np.arange(n) * expected_coefs[1]
+    expected_trend_2D = np.reshape(expected_trend, (n, 1))
+    actual_trend = transformer.forecaster_.predict(-np.arange(n))
+    np.testing.assert_array_almost_equal(actual_trend, expected_trend_2D)
 
     # check residuals
     actual = transformer.transform(y)


Code review comments:
diff_hunk: @@ -197,36 +183,15 @@ def _inverse_transform(self, X, y=None):
         Xt : pd.Series or pd.DataFrame, same type as X
             inverse transformed version of X
         """
-        fh = ForecastingHorizon(X.index, is_relative=False)
-
-        # univariate: X is pd.Series
-        if isinstance(X, pd.Series):
-            # note: the y in the transformer is exogeneous in the forecaster, i.e., X
-            X_pred = self.forecaster_.predict(fh=fh, X=y)
-            if self.model == "additive":
-                return X + X_pred
-            elif self.model == "multiplicative":
-                return X * X_pred
-        # multivariate: X is pd.DataFrame
-        if isinstance(X, pd.DataFrame):
-            X = X.copy()
-            # check if all columns are known
-            X_fit_keys = set(self.forecaster_.keys())
-            X_new_keys = set(X.columns)
-            difference = X_new_keys.difference(X_fit_keys)
-            if len(difference) != 0:
-                raise ValueError(
-                    "X contains columns that have not been "
-                    "seen in fit: " + str(difference)
-                )
-            for colname in X.columns:
-                X_pred = self.forecaster_[colname].predict(fh=fh, X=y)
-                if self.model == "additive":
-                    X[colname] = X[colname] + X_pred
-                elif self.model == "multiplicative":
-                    X[colname] = X[colname] * X_pred
-
-            return X
+        fh = self._get_fh_from_X(X=X)
+        forecaster = self._get_fitted_forecaster(X=X, y=y, fh=fh)
path: sktime/transformations/series/detrend/_detrend.py
text: @user1:
If `X` is the de-trended data, wouldn't passing this to `_get_fitted_forecaster` return a newly fitted forecaster when `if self.forecaster_.get_tag("requires-fh-in-fit", True)`, resulting in the forecaster being fitted on the de-trended data? Is there an example forecaster you have in mind where this is needed to add the trend back to the data?

@author:
ah, good catch! Didn't think carefully about this. I suppose fitting on the detrended data does not really make sense here.

Should it then just raise an error, if `self.forecaster_.get_tag("requires-fh-in-fit", True)`? I.e., one more line with either exception or warning under that condition?

@author:
here's I think what should happen in the `inverse_transform`: the data should *not* be updated with the new data, simply the data in `fit` should be remembered. Changed to this.

@user1:
The only downside is that now we'll be carrying around a copy of the data inside the transformer. Perhaps this is not much of an issue because it's just reference rather than a deep copy of the data. My concern was If I wanted to just use something like the `PolynomialTrendForecaster` to do de-trending and re-trending, then I don't actually require the data after `fit`, I certainly wouldn't want the burden of carrying around a copy of the data unnecessarily. But we don't have the memory burden because it's not a deep copy, so hopefully should be fine.

@author:
well, memorization will not happen for the `PolynomialTrendForecaster` (not even a pointer copy), because its tag is `"requires-fh-in-fit": False`.

The memorization happens *only* in cases where this tag is `True` - a case where the transformer previously would just crash.

In the `False` case, behaviour is precisely as before for univariate forecasters, and the "bugfixed" behaviour for multivariate or hierarchical ones.