Repository: sktime/sktime
Pull request: [ENH] Add a `CurveFitForecaster` based on `scipy` `optimize_curve`
<!--
Thanks for contributing a pull request! Please ensure you have taken a look
at our contribution guide: https://github.com/sktime/sktime/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
Fixes #5238

<!--
Example: Fixes #1234. See also #3456.

Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->


#### What does this implement/fix? Explain your changes.

Adds a `CurveFitForecaster` to `trend.py`. This forecaster takes as
parameters a function and initial parameters. This function is fitted on
`y` and the dateindex using `curve_fit` of scipy.

#### Does your contribution introduce a new dependency? If yes, which one?
No
<!--
If your contribution does add a new hard dependency, we may suggest to initially develop your contribution in a separate companion package in https://github.com/sktime/ to keep external dependencies of the core sktime package to a minimum.
-->

#### What should a reviewer concentrate their feedback on?
Should this integrated in the TrendForecaster?
<!-- This section is particularly useful if you have a pull request that is still in development. You can guide the reviews to focus on the parts that are ready for their comments. We suggest using bullets (indicated by * or -) and filled checkboxes [x] here -->

#### Did you add any tests for the change?
Only test_params
<!-- This section is useful if you have added a test in addition to the existing ones. This will ensure that further changes to these files won't introduce the same kind of bug. It is considered good practice to add tests with newly added code to enforce the fact that the code actually works. This will reduce the chance of introducing logical bugs.
-->

#### Any other comments?
<!--
Please be aware that we are a loose team of volunteers so patience is necessary; assistance handling other issues is very welcome. We value all user contributions, no matter how minor they are. If we are slow to review, either the pull request needs some benchmarking, tinkering, convincing, etc. or more likely the reviewers are simply busy. In either case, we ask for your understanding during the review process.
-->

#### PR checklist
<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [X] I've added myself to the [list of contributors](https://github.com/sktime/sktime/blob/main/CONTRIBUTORS.md) with any new badges I've earned :-)
  How to: add yourself to the [all-contributors file](https://github.com/sktime/sktime/blob/main/.all-contributorsrc) in the `sktime` root directory (not the `CONTRIBUTORS.md`). Common badges: `code` - fixing a bug, or adding code logic. `doc` - writing or improving documentation or docstrings. `bug` - reporting or diagnosing a bug (get this plus `code` if you also fixed the bug in the PR).`maintenance` - CI, test framework, release.
  See here for [full badge reference](https://allcontributors.org/docs/en/emoji-key)
- [ ] Optionally, I've added myself and possibly others to the [CODEOWNERS](https://github.com/sktime/sktime/blob/main/CODEOWNERS) file - do this if you want to become the owner or maintainer of an estimator you added.
  See here for further details on the [algorithm maintainer role](https://www.sktime.net/en/latest/get_involved/governance.html#algorithm-maintainers).
- [X] The PR title starts with either [ENH], [MNT], [DOC], or [BUG]. [BUG] - bugfix, [MNT] - CI, test framework, [ENH] - adding or improving code, [DOC] - writing or improving documentation or docstrings.

##### For new estimators
- [X] I've added the estimator to the API reference - in `docs/source/api_reference/taskname.rst`, follow the pattern.
- [x] I've added one or more illustrative usage examples to the docstring, in a pydocstyle compliant `Examples` section.
- [ ] If the estimator relies on a soft dependency, I've set the `python_dependencies` tag and ensured
  dependency isolation, see the [estimator dependencies guide](https://www.sktime.net/en/latest/developer_guide/dependencies.html#adding-a-soft-dependency).


<!--
Thanks for contributing!
-->


Problem statement:
[ENH] `scipy.curve_fit` based forecaster
**Is your feature request related to a problem? Please describe.**
In scipy there exist the `curve_fit` method. This takes a function, x, and y data and some initial parameters. If x are time series indexes this can be used to perform forecasting: 


**Describe the solution you'd like**
Example usage would looks like that:
```python
def sine_3(x,amplitude, phase, offset, amplitude2, amplitude3, phase2):
    sbase = np.sin(x * (1 / 365 / 24 * np.pi * 4) + phase) * amplitude + offset 
    s1 = amplitude2 * x * (np.sin(x* (1 / 365 / 24 * np.pi * 2) + (365* 24* np.pi)+ phase2) + 1)
    s2 = amplitude3 * x * (np.sin(x* (1 / 365 / 24 * np.pi * 2) + phase2) - 1)
    return sbase + s1 + s2


function = FunctionForecaster(function=sine_3, initial_params=[3 * np.std(rolling_mean) / (2 ** 0.5), 
                                         0, 
                                         rolling_mean.mean(),
                                         10,
                                         1,
                                         1])

function.fit(rolling_mean)
function.predict(fh=ForecastingHorizon(y_test.index, is_relative=False))
```

During the `fit`, `curve_fit` is used to detect the best parameters. During `predict`, the found parameters and the fh are used to forecast the time series. 



**Describe alternatives you've considered**
* Implement a new forecaster.
* Extend TrendForecaster (idea from @fkiraly)


**Additional context**

The fit code would look as follows: 

```python
def fit(y, X, fh):
        t = ForecastingHorizon(y.index, is_relative=False).to_relative(pd.DatetimeIndex([self._y.index[0]], name=self.cutoff.name, freq=self.cutoff.freq))
    
        self.params_ = curve_fit(self.function, list(t), y.values, self.initial_params)
        return self
```
and the predict code as:
```python

def predict(fh, X):
        t = fh.to_relative(pd.DatetimeIndex([self._y.index[0]], name=self.cutoff.name, freq=self.cutoff.freq))
        return pd.Series(self.function(np.array(t), 
                             *self.params_[0]), index=pd.DatetimeIndex(fh.to_absolute(self.cutoff)))

```

Note using `[self._y.index[0]]` as cutoff  in `to_relative` has the effect that the x values used in the function are natural numbers.

Hints from issue discussion:
> Describe alternatives you've considered
> Implement a new forecaster.
> Extend `TrendForecaster` (idea from @fkiraly)

FYI, looking at your PR and the code, I think writing a new estimator is probably the better solution.

Commit message: Fix doc string issues

Patch to review:
diff --git a/docs/source/api_reference/forecasting.rst b/docs/source/api_reference/forecasting.rst
index 3506e711f10..7715dbd5adc 100644
--- a/docs/source/api_reference/forecasting.rst
+++ b/docs/source/api_reference/forecasting.rst
@@ -149,6 +149,7 @@ Trend forecasters
     TrendForecaster
     PolynomialTrendForecaster
     STLForecaster
+    FunctionForecaster
 
 .. currentmodule:: sktime.forecasting.statsforecast
 
diff --git a/sktime/forecasting/trend/__init__.py b/sktime/forecasting/trend/__init__.py
index b7724ab8c04..f1df869127b 100644
--- a/sktime/forecasting/trend/__init__.py
+++ b/sktime/forecasting/trend/__init__.py
@@ -6,8 +6,10 @@
     "TrendForecaster",
     "PolynomialTrendForecaster",
     "STLForecaster",
+    "CurveFitForecaster",
 ]
 
+from sktime.forecasting.trend._curve_fit_forecaster import CurveFitForecaster
 from sktime.forecasting.trend._polynomial_trend_forecaster import (
     PolynomialTrendForecaster,
 )
diff --git a/sktime/forecasting/trend/_curve_fit_forecaster.py b/sktime/forecasting/trend/_curve_fit_forecaster.py
new file mode 100644
index 00000000000..1aa588a5542
--- /dev/null
+++ b/sktime/forecasting/trend/_curve_fit_forecaster.py
@@ -0,0 +1,177 @@
+# !/usr/bin/env python3 -u
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Implements FunctionForecaster."""
+
+__author__ = ["benheid"]
+__all__ = ["CurveFitForecaster"]
+
+import numpy as np
+import pandas as pd
+from scipy.optimize import curve_fit
+
+from sktime.forecasting.base import BaseForecaster, ForecastingHorizon
+
+
+class CurveFitForecaster(BaseForecaster):
+    """The CurveFitForecaster takes a function and tits it by using scipy curve_fit.
+
+    The CurveFitForecaster applies the scipy curve_fit method to find the best
+    parameter of a function. This function maps a list of integers to real values.
+    The integers are derrived from the index of the input time series.
+
+    In `fit`
+    1. The index of the input time series is transformed to a list of Integers.
+    2. The scipy curve_fit is called using the list of integers as x values,
+       and the time series values as y values.
+    In `predict`
+    1. The ForecastingHorizon is transformed to a list of integers.
+    2. The list of integers is passed together with the fitted parameters to the
+       function to provide the forecast.
+
+
+    Parameters
+    ----------
+    function: Callable
+        The function that should be fitted and used to make forecasts.
+        The signature of the functions is `function(x, ...)`.
+        It takes the independet variables as first argument and the parametrs
+        to fit as separate remaining arguments.
+        See scipy.optimize.curve_fit for more information.
+    curve_fit_params: dict, default=None
+        Additional parameters that should be passed to the curve_fit method.
+        See scipy.optimize.curve_fit for more information.
+
+
+    Examples
+    --------
+    >>> from sktime.forecasting.trend import CurveFitForecaster
+    >>> from sktime.datasets import load_airline
+    >>> y = load_airline()
+    >>> def linear_function(x, a, b):
+    ...     return a * x + b
+    >>> forecaster = CurveFitForecaster(function=linear_function,
+    ...                                 curve_fit_params={"p0":[-1, 1]})
+    >>> forecaster.fit(y)
+    CurveFitForecaster(...)
+    >>> y_pred = forecaster.predict(fh=[1, 2, 3])
+    """
+
+    _tags = {
+        "scitype:y": "univariate",
+        "y_inner_mtype": "pd.Series",
+        "X_inner_mtype": "pd.DataFrame",
+        "ignores-exogeneous-X": True,
+        "requires-fh-in-fit": False,
+    }
+
+    def __init__(
+        self,
+        function,
+        curve_fit_params=None,
+    ):
+        self.function = function
+        self.curve_fit_params = curve_fit_params
+        self._curve_fit_params = (
+            curve_fit_params if curve_fit_params is not None else {}
+        )
+        super().__init__()
+
+    def _fit(self, y, X=None, fh=None):
+        """Fit forecaster to training data.
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        y : pd.DataFrame
+            time series to which to fit the forecaster.
+        fh : guaranteed to be ForecastingHorizon or None, optional (default=None)
+            The forecasting horizon with the steps ahead to predict.
+            Ignored in fit
+        X : pd.DataFrame, optional (default=None)
+            Ignored in fit
+
+        Returns
+        -------
+        self : reference to self
+        """
+        t = ForecastingHorizon(y.index, is_relative=False).to_relative(self.cutoff)
+        self.params_ = curve_fit(
+            self.function, np.array(t), y.values, **self._curve_fit_params
+        )
+        return self
+
+    def _predict(self, fh, X=None):
+        """Forecast time series at future horizon.
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
+            self.cutoff
+
+        Parameters
+        ----------
+        fh : guaranteed to be ForecastingHorizon,
+            The forecasting horizon with the steps ahead to to predict.
+        X : pd.DataFrame, optional (default=None)
+            Ignored in predict
+
+        Returns
+        -------
+        y_pred : pd.DataFrame
+            Point predictions
+        """
+        t = fh.to_relative(self.cutoff)
+        return pd.Series(
+            self.function(np.array(t), *self.params_[0]),
+            index=list(fh.to_absolute(self.cutoff)),
+            name=self._y.name,
+        )
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for forecasters.
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class.
+        """
+        params1 = {
+            "function": _test_function,
+            "curve_fit_params": {
+                "p0": [1, 1],
+                "sigma": None,
+                "absolute_sigma": True,
+                "check_finite": True,
+                "bounds": [(-100, -100), (100, 100)],
+                "method": "dogbox",
+            },
+        }
+
+        params2 = {
+            "function": _test_function_2,
+            "curve_fit_params": {
+                "method": "trf",
+            },
+        }
+
+        return [params1, params2]
+
+
+def _test_function(x, a, b):
+    return a * x + b
+
+
+def _test_function_2(x, a, b, c):
+    return a * x**2 + b * x + c


Code review comments:
diff_hunk: @@ -0,0 +1,177 @@
+# !/usr/bin/env python3 -u
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Implements FunctionForecaster."""
+
+__author__ = ["benheid"]
+__all__ = ["CurveFitForecaster"]
+
+import numpy as np
+import pandas as pd
+from scipy.optimize import curve_fit
+
+from sktime.forecasting.base import BaseForecaster, ForecastingHorizon
+
+
+class CurveFitForecaster(BaseForecaster):
+    """The CurveFitForecaster takes a function and tits it by using scipy curve_fit.
path: sktime/forecasting/trend/_curve_fit_forecaster.py
text: @user1:
typo, fits

@author:
Done

diff_hunk: @@ -0,0 +1,177 @@
+# !/usr/bin/env python3 -u
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+"""Implements FunctionForecaster."""
+
+__author__ = ["benheid"]
+__all__ = ["CurveFitForecaster"]
+
+import numpy as np
+import pandas as pd
+from scipy.optimize import curve_fit
+
+from sktime.forecasting.base import BaseForecaster, ForecastingHorizon
+
+
+class CurveFitForecaster(BaseForecaster):
+    """The CurveFitForecaster takes a function and tits it by using scipy curve_fit.
+
+    The CurveFitForecaster applies the scipy curve_fit method to find the best
+    parameter of a function. This function maps a list of integers to real values.
+    The integers are derrived from the index of the input time series.
+
+    In `fit`
+    1. The index of the input time series is transformed to a list of Integers.
+    2. The scipy curve_fit is called using the list of integers as x values,
+       and the time series values as y values.
+    In `predict`
+    1. The ForecastingHorizon is transformed to a list of integers.
+    2. The list of integers is passed together with the fitted parameters to the
+       function to provide the forecast.
+
+
+    Parameters
+    ----------
+    function: Callable
+        The function that should be fitted and used to make forecasts.
+        The signature of the functions is `function(x, ...)`.
+        It takes the independet variables as first argument and the parametrs
+        to fit as separate remaining arguments.
+        See scipy.optimize.curve_fit for more information.
+    curve_fit_params: dict, default=None
+        Additional parameters that should be passed to the curve_fit method.
+        See scipy.optimize.curve_fit for more information.
+
+
+    Examples
+    --------
+    >>> from sktime.forecasting.trend import CurveFitForecaster
+    >>> from sktime.datasets import load_airline
+    >>> y = load_airline()
+    >>> def linear_function(x, a, b):
+    ...     return a * x + b
+    >>> forecaster = CurveFitForecaster(function=linear_function,
+    ...                                 curve_fit_params={"p0":[-1, 1]})
+    >>> forecaster.fit(y)
+    CurveFitForecaster(...)
+    >>> y_pred = forecaster.predict(fh=[1, 2, 3])
+    """
+
+    _tags = {
+        "scitype:y": "univariate",
+        "y_inner_mtype": "pd.Series",
+        "X_inner_mtype": "pd.DataFrame",
+        "ignores-exogeneous-X": True,
+        "requires-fh-in-fit": False,
+    }
+
+    def __init__(
+        self,
+        function,
+        curve_fit_params=None,
+    ):
+        self.function = function
+        self.curve_fit_params = curve_fit_params
+        self._curve_fit_params = (
+            curve_fit_params if curve_fit_params is not None else {}
+        )
+        super().__init__()
+
+    def _fit(self, y, X=None, fh=None):
+        """Fit forecaster to training data.
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        y : pd.DataFrame
+            time series to which to fit the forecaster.
+        fh : guaranteed to be ForecastingHorizon or None, optional (default=None)
+            The forecasting horizon with the steps ahead to predict.
+            Ignored in fit
+        X : pd.DataFrame, optional (default=None)
+            Ignored in fit
+
+        Returns
+        -------
+        self : reference to self
+        """
+        t = ForecastingHorizon(y.index, is_relative=False).to_relative(self.cutoff)
path: sktime/forecasting/trend/_curve_fit_forecaster.py
text: @user1:
will this always result in integers?

@author:
I revised this. Now, to extract the integer values, I use partially the `trend.utils._get_X_numpy_int_from_pandas`

However, to enable an easier definition of the function without the need to consider that zero is 1970, I set the zero to the index of the index of the first value of the training data. Moreover, I the distance between two consecutive values to one.