Repository: sktime/sktime
Pull request: [BUG] fix erroneous `int` coercion of `TrendForecaster` and `PolynomialTrendForecaster` on `DatetimeIndex`
Fixes https://github.com/sktime/sktime/issues/4131 and a parallel issue for `TrendForecaster` which is the cause of #3230 (also fixes #3230).

Erroneous coercion caused huge arrays to be allocated when the input data had a `DatetimeIndex`.

Also, previously, any time stamp based fits and coefficients were related to a reset index (naive replacement with `RangeIndex`) rather than the time stamp.

This has been changed to periods since start of 1970 (with correct time stamp conversion). I've now added a clarifying docstring on what the coefficients mean, and changed it consistently to periods since start of 1970.

Since this was subject to a bug previously (and, likely, confusion), we don't need deprecation for this change.

Contains new and updated tests:
*  `test_trendforecaster_with_datetimeindex`, a test for the failure case in #4131, and parallel case.
* The existing test `test_trend` was failing with the bugfix due to the erroneous conversion to `int`. The test has been modified to cover the time stamp case corretly, by using time stamp converted indices in the Vandermonde matrix rather than the `RangeIndex` based indices.
* a test for the time stamp to integer/float conversion logic now used under the hood

Problem statement:
[BUG] `TrendForecaster` gives incorrect output
**Describe the bug**
The TrendForecaster produces unusual outputs on levels different from the target series.

**To Reproduce**

```python
from sktime.datasets import load_airline
from sktime.forecasting.trend import TrendForecaster
from sktime.utils.plotting import plot_series

y = load_airline()
forecaster = TrendForecaster()

y_preds = []
for last_train in [20, 40, 60, 80]:
    forecaster.fit(y[:last_train])
    y_pred = forecaster.predict(fh=[1,2,3,4,5])
    y_preds.append(y_pred)

plot_series(y, *y_preds)
```

**Expected behavior**
The forecasts should be straight lines extrapolating the linear historical trend.

**Additional context**
Would expect that the default behaviour is similar to the default behaviour of the `PolynomialTrendForecaster`. In `TrendForecaster`, the `X` is calculated as
```python
X = y.index.astype("int").to_numpy().reshape(-1, 1)
```
whereas in the polynomial version, it is calculated as
```python
n_timepoints = _get_duration(self._y.index, coerce_to_int=True) + 1
X = np.arange(n_timepoints).reshape(-1, 1)
```



**Versions**
<details>

System:
    python: 3.8.13 (default, Apr 19 2022, 00:53:22)  [GCC 7.5.0]
executable: /usr/bin/python3.8
   machine: Linux-5.4.72-microsoft-standard-WSL2-x86_64-with-glibc2.27

Python dependencies:
          pip: 22.0.4
   setuptools: 56.0.0
      sklearn: 1.1.2
       sktime: 0.13.0
  statsmodels: 0.13.2
        numpy: 1.22.4
        scipy: 1.8.1
       pandas: 1.4.3
   matplotlib: 3.5.2
       joblib: 1.1.0
        numba: 0.56.0
     pmdarima: None
      tsfresh: None

</details>

[BUG] `PolynomialTrendForecaster` stalls when fit on timeseries with datetime index
**Describe the bug**
My computer stalls when running `fit` on a dataframe with a `DatetimeIndex` with the `PolynomialTrendForecaster`.

**To Reproduce**
```python
from sktime.datasets import load_airline
from sktime.forecasting.trend import PolynomialTrendForecaster  

df = load_airline()
df.index = df.index.to_timestamp()

trafo = PolynomialTrendForecaster()
trafo.fit(df)
```

**Expected behavior**
The fitting time shouldn't be orders of magnitude different from using a `PeriodIndex`. 


**Versions**
<details>

System:
    python: 3.8.7 (default, May  6 2021, 21:53:45)  [Clang 12.0.0 (clang-1200.0.32.29)]
executable: /Users/kishan_manani/.pyenv/versions/3.8.7/envs/udemy-ts/bin/python3.8
   machine: macOS-11.6.6-x86_64-i386-64bit

Python dependencies:
          pip: 22.3.1
   setuptools: 49.2.1
      sklearn: 1.2.0
       sktime: 0.15.0
  statsmodels: 0.13.2
        numpy: 1.21.6
        scipy: 1.8.1
       pandas: 1.5.2
   matplotlib: 3.6.2
       joblib: 1.2.0
        numba: 0.56.4
     pmdarima: 1.8.5
      tsfresh: 0.19.0

</details>

<!-- Thanks for contributing! -->

Hints from issue discussion:
This should be fixed now here: https://github.com/sktime/sktime/pull/4133

Cause is the same as bug https://github.com/sktime/sktime/issues/4131 - conversion to integer was erroneous for time stamp indexed time series.
```
Unable to allocate 2.74 TiB for an array with shape (376012800001,) and data type int64
```
???

This is certainly a bug.
I suspect that the large number comes from a hidden coercion from a date format to integer (seconds sine 1970 or similar).

This indicates that we should ensure that `DateTimeIndex` is covered in scenario test cases.
Hope this fixes it: https://github.com/sktime/sktime/pull/4133

Commit message: higher rtol for numerical tests

Patch to review:
diff --git a/sktime/forecasting/tests/test_trend.py b/sktime/forecasting/tests/test_trend.py
index 1034649fc3d..e3cdac78765 100644
--- a/sktime/forecasting/tests/test_trend.py
+++ b/sktime/forecasting/tests/test_trend.py
@@ -2,20 +2,60 @@
 """Test trend forecasters."""
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 
-__author__ = ["mloning"]
+__author__ = ["mloning", "fkiraly"]
 __all__ = ["get_expected_polynomial_coefs"]
 
 import numpy as np
 import pandas as pd
 import pytest
 
-from sktime.forecasting.trend import PolynomialTrendForecaster
+from sktime.datasets import load_airline
+from sktime.forecasting.trend import (
+    PolynomialTrendForecaster,
+    TrendForecaster,
+    _get_X_numpy_int_from_pandas,
+)
 from sktime.utils._testing.forecasting import make_forecasting_problem
 
 
+def test_get_X_numpy():
+    """Test _get_X_numpy_int_from_pandas converts to int/float as expected."""
+    y = load_airline()
+    X_idx = _get_X_numpy_int_from_pandas(y.index)
+
+    # testing pd.PeriodIndex
+    # this should be a 2D nd.nparray, with diffs being 1 (month)
+    # because month is the periodicity
+    assert isinstance(X_idx, np.ndarray)
+    assert X_idx.shape == (len(y), 1)
+    assert (np.diff(X_idx.reshape(-1)) == 1).all()
+
+    y_idx_datetime = y.index.to_timestamp()
+    X_idx_datetime = _get_X_numpy_int_from_pandas(y_idx_datetime)
+
+    # testing pd.DatetimeIndex
+    # this should be a 2D np.ndarray, with diffs being 28, 29, 30, or 31 (days)
+    # because DatetimeIndex time stamps convert to days
+    assert isinstance(X_idx_datetime, np.ndarray)
+    assert X_idx_datetime.shape == (len(y), 1)
+    intdiffs = (np.diff(X_idx_datetime.reshape(-1))).astype(int)
+    assert np.isin(intdiffs, [30, 31, 28, 29]).all()
+
+    # testing integer index
+    int_ix = pd.Index([1, 3, 5, 7])
+    X_idx_int = _get_X_numpy_int_from_pandas(int_ix)
+
+    # testing pd.IntegerIndex
+    # this should be an 2D integer array with entries identical to int_ix
+    assert isinstance(X_idx_int, np.ndarray)
+    assert X_idx_int.shape == (len(int_ix), 1)
+    assert (X_idx_int.reshape(-1) == int_ix.to_numpy()).all()
+
+
 def get_expected_polynomial_coefs(y, degree, with_intercept=True):
     """Compute expected coefficients from polynomial regression."""
-    poly_matrix = np.vander(np.arange(len(y)), degree + 1)
+    t_ix = _get_X_numpy_int_from_pandas(y.index).reshape(-1)
+    poly_matrix = np.vander(t_ix, degree + 1)
     if not with_intercept:
         poly_matrix = poly_matrix[:, :-1]
     return np.linalg.lstsq(poly_matrix, y.to_numpy(), rcond=None)[0]
@@ -31,21 +71,24 @@ def _test_trend(degree, with_intercept):
     # intercept is added in reverse order
     actual = forecaster.regressor_.steps[-1][1].coef_[::-1]
     expected = get_expected_polynomial_coefs(y, degree, with_intercept)
-    np.testing.assert_allclose(actual, expected)
+    np.testing.assert_allclose(actual, expected, rtol=1e-5)
 
 
 @pytest.mark.parametrize("degree", [1, 3])
 @pytest.mark.parametrize("with_intercept", [True, False])
 def test_trend(degree, with_intercept):
+    """Test PolynomialTrendForecaster coefficients."""
     _test_trend(degree, with_intercept)
 
 
 # zero trend does not work without intercept
 def test_zero_trend():
+    """Test PolynomialTrendForecaster with degree zero."""
     _test_trend(degree=0, with_intercept=True)
 
 
 def test_constant_trend():
+    """Test expected output from constant trend."""
     y = pd.Series(np.arange(30))
     fh = -np.arange(30)  # in-sample fh
 
@@ -53,3 +96,15 @@ def test_constant_trend():
     y_pred = forecaster.fit(y).predict(fh)
 
     np.testing.assert_array_almost_equal(y, y_pred)
+
+
+def test_trendforecaster_with_datetimeindex():
+    """Test PolyonmialTrendForecaster with DatetimeIndex, see #4131."""
+    df = load_airline()
+    df.index = df.index.to_timestamp()
+
+    f = PolynomialTrendForecaster()
+    f.fit(df)
+
+    f = TrendForecaster()
+    f.fit(df)
diff --git a/sktime/forecasting/trend.py b/sktime/forecasting/trend.py
index 10a9924db6e..cfa7efd2915 100644
--- a/sktime/forecasting/trend.py
+++ b/sktime/forecasting/trend.py
@@ -3,10 +3,9 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Implements trend based forecasters."""
 
-__author__ = ["tensorflow-as-tf", "mloning", "aiwalter"]
+__author__ = ["tensorflow-as-tf", "mloning", "aiwalter", "fkiraly"]
 __all__ = ["TrendForecaster", "PolynomialTrendForecaster", "STLForecaster"]
 
-import numpy as np
 import pandas as pd
 from sklearn.base import clone
 from sklearn.linear_model import LinearRegression
@@ -14,13 +13,36 @@
 from sklearn.preprocessing import PolynomialFeatures
 
 from sktime.forecasting.base import BaseForecaster
-from sktime.utils.datetime import _get_duration
+
+
+def _get_X_numpy_int_from_pandas(x):
+    """Convert pandas index to an sklearn compatible X, 2D np.ndarray, int type."""
+    if isinstance(x, (pd.DatetimeIndex)):
+        x = x.astype("int64") / 864e11
+    else:
+        x = x.astype("int64")
+    return x.to_numpy().reshape(-1, 1)
 
 
 class TrendForecaster(BaseForecaster):
-    """Trend based forecasts of time series data.
+    r"""Trend based forecasts of time series data, regressing values on index.
+
+    Uses a `sklearn` regressor `regressor` to regress values of time series on index:
+
+    In `fit`, for input time series :math:`(v_i, t_i), i = 1, \dots, T`,
+    where :math:`v_i` are values and :math:`t_i` are time stamps,
+    fits an `sklearn` model :math:`v_i = f(t_i) + \epsilon_i`, where `f` is
+    the model fitted when `regressor.fit` is passed `X` = vector of :math:`t_i`,
+    and `y` = vector of :math:`v_i`.
+
+    In `predict`, for a new time point :math:`t_*`, predicts :math:`f(t_*)`,
+    where :math:`f` is the function as fitted above in `fit`.
 
-    Default settings train a linear regression model.
+    Default for `regressor` is linear regression = `sklearn` `LinearRegression` default.
+
+    If time stamps are `pd.DatetimeIndex`, fitted coefficients are in units
+    of days since start of 1970. If time stamps are `pd.PeriodIndex`,
+    coefficients are in units of (full) periods since start of 1970.
 
     Parameters
     ----------
@@ -71,11 +93,13 @@ def _fit(self, y, X=None, fh=None):
         else:
             self.regressor_ = clone(self.regressor)
 
-        # transform data
-        X = y.index.astype("int64").to_numpy().reshape(-1, 1)
+        # we regress index on series values
+        # the sklearn X is obtained from the index of y
+        # the sklearn y can be taken as the y seen here
+        X_sklearn = _get_X_numpy_int_from_pandas(y.index)
 
         # fit regressor
-        self.regressor_.fit(X, y)
+        self.regressor_.fit(X_sklearn, y)
         return self
 
     def _predict(self, fh=None, X=None):
@@ -94,9 +118,9 @@ def _predict(self, fh=None, X=None):
             Point predictions for the forecast
         """
         # use relative fh as time index to predict
-        fh = self.fh.to_absolute_int(self._y.index[0], self.cutoff)
-        X_pred = fh.to_numpy().reshape(-1, 1)
-        y_pred = self.regressor_.predict(X_pred)
+        fh = self.fh.to_absolute(self.cutoff)
+        X_sklearn = _get_X_numpy_int_from_pandas(fh.to_pandas())
+        y_pred = self.regressor_.predict(X_sklearn)
         return pd.Series(y_pred, index=self.fh.to_absolute(self.cutoff))
 
     @classmethod
@@ -125,10 +149,30 @@ def get_test_params(cls, parameter_set="default"):
 
 
 class PolynomialTrendForecaster(BaseForecaster):
-    """Forecast time series data with a polynomial trend.
+    r"""Forecast time series data with a polynomial trend.
+
+    Uses a `sklearn` regressor `regressor` to regress values of time series on index,
+    after extraction of polynomial features.
+    Same `TrendForecaster` where `regressor` is pipelined with transformation step
+    `PolynomialFeatures(degree, with_intercept)` applied to time, at the start.
+
+    In `fit`, for input time series :math:`(v_i, p(t_i)), i = 1, \dots, T`,
+    where :math:`v_i` are values, :math:`t_i` are time stamps,
+    and :math:`p` is the polynomial feature transform with degree `degree`,
+    and with/without intercept depending on `with_intercept`,
+    fits an `sklearn` model :math:`v_i = f(p(t_i)) + \epsilon_i`, where `f` is
+    the model fitted when `regressor.fit` is passed `X` = vector of :math:`p(t_i)`,
+    and `y` = vector of :math:`v_i`.
+
+    In `predict`, for a new time point :math:`t_*`, predicts :math:`f(p(t_*))`,
+    where :math:`f` is the function as fitted above in `fit`,
+    and :math:`p` is the same polynomial feature transform as above.
+
+    Default for `regressor` is linear regression = `sklearn` `LinearRegression` default.
 
-    Default settings train a linear regression model with a 1st degree
-    polynomial transformation of the feature.
+    If time stamps are `pd.DatetimeIndex`, fitted coefficients are in units
+    of days since start of 1970. If time stamps are `pd.PeriodIndex`,
+    coefficients are in units of (full) periods since start of 1970.
 
     Parameters
     ----------
@@ -195,12 +239,13 @@ def _fit(self, y, X=None, fh=None):
             regressor,
         )
 
-        # transform data
-        n_timepoints = _get_duration(self._y.index, coerce_to_int=True) + 1
-        X = np.arange(n_timepoints).reshape(-1, 1)
+        # we regress index on series values
+        # the sklearn X is obtained from the index of y
+        # the sklearn y can be taken as the y seen here
+        X_sklearn = _get_X_numpy_int_from_pandas(y.index)
 
         # fit regressor
-        self.regressor_.fit(X, y)
+        self.regressor_.fit(X_sklearn, y)
         return self
 
     def _predict(self, fh=None, X=None):
@@ -219,9 +264,9 @@ def _predict(self, fh=None, X=None):
             Point predictions for the forecast
         """
         # use relative fh as time index to predict
-        fh = self.fh.to_absolute_int(self._y.index[0], self.cutoff)
-        X_pred = fh.to_numpy().reshape(-1, 1)
-        y_pred = self.regressor_.predict(X_pred)
+        fh = self.fh.to_absolute(self.cutoff)
+        X_sklearn = _get_X_numpy_int_from_pandas(fh.to_pandas())
+        y_pred = self.regressor_.predict(X_sklearn)
         return pd.Series(y_pred, index=self.fh.to_absolute(self.cutoff))
 
     @classmethod


Code review comments:
diff_hunk: @@ -2,20 +2,60 @@
 """Test trend forecasters."""
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 
-__author__ = ["mloning"]
+__author__ = ["mloning", "fkiraly"]
 __all__ = ["get_expected_polynomial_coefs"]
 
 import numpy as np
 import pandas as pd
 import pytest
 
-from sktime.forecasting.trend import PolynomialTrendForecaster
+from sktime.datasets import load_airline
+from sktime.forecasting.trend import (
+    PolynomialTrendForecaster,
+    TrendForecaster,
+    _get_X_numpy_int_from_pandas,
+)
 from sktime.utils._testing.forecasting import make_forecasting_problem
 
 
+def test_get_X_numpy():
+    """Test _get_X_numpy_int_from_pandas converts to int/float as expected."""
+    y = load_airline()
+    X_idx = _get_X_numpy_int_from_pandas(y.index)
+
+    # testing pd.PeriodIndex
+    # this should be a 2D nd.nparray, with diffs being 1 (month)
path: sktime/forecasting/tests/test_trend.py
text: @user1:
Nit: there is a typo here `nd.nparray`. 

@author:
fixed