Repository: sktime/sktime
Pull request: [ENH] Forecaster convenience method to return forecast residuals
This PR implements a frequently requested convenience method (`predict_residuals`) for forecasters to return forecast residuals. See also #1105.
Tests for `predict_residuals` are also included (checks whether indices are correct).

FYI @ngupta23, @AngelPone, is this what you were looking for?

Problem statement:
Forecaster interface: getting in-sample forecasts, getting residuals
**Is your feature request related to a problem? Please describe.**
I think there should be a method or property to obtain in-sample fitted value or residuals, which are important in many forecasting methods, such as forecast combination.

**Describe the solution you'd like**

Add a method or property to `BaseForecaster` and it should be implemented for every forecaster.

**Additional context**
A traditional forecast combination method based on variance of in-sample residuals.

Hints from issue discussion:
Ok, I think there are two requests in there, and a proposed way of addressing these.

The two requests: residuals and in-sample values.

Regarding residuals, I fully agree that there should be an easy way to return them!
But, I don't think this is a concern of the `BaseForecaster`.

There also already is a way to obtain residuals from a forecaster, it's the `Detrender` - simply detrend with the forecaster you would like to get residuals of - it's a transformer, and the output gives you the residuals.
Perhaps we should say this in the forecasting tutorial, since it's not obvious from the name "detrender"?

Regarding in-sample fitted values: this is definitely a `BaseForecaster` concern in my opinion, but doesn't this already exist? Is this not what you would get when you run `predict` at the values that have already been seen? Simply pass a fixed forecasting horizon with the indices that you've observed.

`forecaster.fit(y).predict(fh=y.index)` or `forecaster.fit(y, fh=y.index).predict()`. Does this work for you?

What I'm wondering, should there perhaps be a `fit_predict(y)`, which, if `fh` is not passed, defaults to in-sample forecasts, i.e., `fh=y.index`? There currently isn't, but it might be the "in-sample" convenience method.

Feedback on whether it is what you need, or whether it's "obvious" would be appreciated.


for forecast combination, you may like to look at the `HedgeExpertEnsemble`, it may do what you want, i.e., re-weight depending on errors or forecast residuals.

Would be nice to hear from you whether you think that's what you need or not, and whether it should be made clearer in the tutorial.
@fkiraly and @AngelPone outside of the forecast combination case, I could see it being useful to be able to expect the residuals from a model. For example, to see if they have a clear pattern (suggesting the model doesn't fit well) or run a statistical test on them (e.g. are they normally distributed in the case of models like Arima that make this assumption).

Adding a convenience method seems straight-forward (but I could be missing something) -- something like:

```python
class BaseForecaster(BaseEstimator):
    def insample_residuals(self):
        self.check_is_fitted()
        if self.get_tag("requires-fh-in-fit", False):
             # Not sure if we should use `_fit` or `fit` here
             # Maybe `_fit` since it doesn't override the stored `fh`?
            return self._fit(self._y, fh=self._y.index).predict()
        else:
            return self.predict(fh=self._y.index)
```

Regarding residuals:
from a [discussion on Slack](https://sktime-group.slack.com/archives/CQVH0NJ1M/p1628935128001400) with @ngupta23 : As some models e.g. from `statsmodels` have residuals internally calculated, the best would be to have a default function `def _get_residuals(self)` or sth in the parent which is called then in `fit()` after` _fit()` was called that calculates residuals based on in-sample predictions. So we could overwrite this function in the child class for some models that provide this internally. So the parent class would get an additional attribute `residuals_`. Not sure if this should be optionally or always calculated. The reason for optionally would be that for large data sets, we maybe dont want to store the residuals in `self`.
What does everybody think?
alternative suggestion: a `ResidualsTrafo` transformer that wrape a forecaster? That way, it can be used as a transformer in pipelines.

We could do both, of course - though let's think carefully before we proliferate base class methods.
I think you can just run:
```python
transformer = Detrender(MyForecaster())
residuals = transformer.fit_transform(y)
```
oh yes, we already have the `ResidualsTrafo` and it's called `Detrender`.

Perhaps we should introduce an alias `Residuals` and/or add a "computing residuals" section to the tutorial?
So just to clarify, `Detrender` will give the true in-sample residuals (in-sample prediction MINUS in-sample actual)?
I think so, yes.
I see [it now](https://github.com/alan-turing-institute/sktime/blob/3c54d8d1e66ec7863d12738e64c1f733d0e6a22e/sktime/transformations/series/detrend/_detrend.py#L141). 

> Perhaps we should introduce an alias Residuals and/or add a "computing residuals" section to the tutorial?

Yes, that will be a useful addition (from a clarity perspective) even though it will be an alias.

Commit message: typos in error msg

Patch to review:
diff --git a/sktime/forecasting/base/_base.py b/sktime/forecasting/base/_base.py
index cc4b5397592..17b505a3178 100644
--- a/sktime/forecasting/base/_base.py
+++ b/sktime/forecasting/base/_base.py
@@ -8,11 +8,17 @@ class name: BaseForecaster
 Scitype defining methods:
     fitting            - fit(y, X=None, fh=None)
     forecasting        - predict(fh=None, X=None)
+    updating           - update(y, X=None, update_params=True)
+
+Convenience methods:
     fit&forecast       - fit_predict(y, X=None, fh=None)
+    update&forecast    - update_predict(cv=None, X=None, update_params=True)
+    forecast residuals - predict_residuals(y, X=None, fh=None)
+    forecast scores    - score(y, X=None, fh=None)
+
+Optional, special capability methods (check capability tags if available):
     forecast intervals - predict_interval(fh=None, X=None, coverage=0.90)
     forecast quantiles - predict_quantiles(fh=None, X=None, alpha=[0.05, 0.95])
-    updating           - update(y, X=None, update_params=True)
-    update&forecast    - update_predict(cv=None, X=None, update_params=True)
 
 Inspection methods:
     hyper-parameter inspection  - get_params()
@@ -631,6 +637,48 @@ def update_predict_single(
             alpha=alpha,
         )
 
+    def predict_residuals(self, y, X=None, fh=None):
+        """Return residuals of time series forecasts.
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_".
+            self.cutoff, self._is_fitted
+
+        Writes to self:
+            Stores fh to self.fh if fh is passed and has not been passed previously.
+
+        Parameters
+        ----------
+        y : pd.Series, pd.DataFrame, or np.ndarray (1D or 2D)
+            Time series with ground truth observations, to compute residuals to.
+            Must have same type, dimension, and indices as expected return of predict.
+        X : pd.DataFrame, or 2D np.ndarray, optional (default=None)
+            Exogeneous time series to predict from
+            if self.get_tag("X-y-must-have-same-index"), X.index must contain fh.index
+        fh : int, list, np.ndarray or ForecastingHorizon
+            Forecasting horizon
+
+        Returns
+        -------
+        y_res : pd.Series, pd.DataFrame, or np.ndarray (1D or 2D)
+            Forecast residuals at fh, with same index as fh
+            y_pred has same type as y passed in fit (most recently)
+        """
+        y_pred = self.predict(fh=fh, X=X)
+
+        if not type(y_pred) == type(y):
+            raise TypeError(
+                "y must have same type, dims, index as expected predict return. "
+                f"expected type {type(y_pred)}, but found {type(y)}"
+            )
+
+        y_res = y - y_pred
+
+        return y_res
+
     def score(self, y, X=None, fh=None):
         """Scores forecast against ground truth, using MAPE.
 


Code review comments:
diff_hunk: @@ -631,6 +637,48 @@ def update_predict_single(
             alpha=alpha,
         )
 
+    def predict_residuals(self, y, X=None, fh=None):
path: sktime/forecasting/base/_base.py
text: @user1:
Shouldn't `y` be optional for in-sample predictions?

@author:
this is not necessarily in-sample, but we could make it that if no `y` is passed, in-sample residuals are computed, and the already seen `_y` are used. If that's convenient?

@user1:
So just confirming --- 

## For in sample residuals

```python
# Option 1: All In-sample residuals
forecaster.predict_residuals()  

# Last "n" (n=4 in this case) in-sample residuals
forecaster.predict_residuals(fh=[0, -1, -2, -3]) 
```

## For out of sample, it would be 

```python
# Option 1: # fh would be same length as Sth
forecaster.predict_residuals(y=Sth)  

# Option 2
# Would non-matching indices (between y and predictions) be ignored or raise error or warning?
forecaster.predict_residuals(y=Sth1, fh=Sth2)  
```




@author:
ok, here's with what I'm currently thinking.

if you want in-sample residuals:
```python
# if fh must be passed in fit
forecaster.fit(y, fh=y.index).predict_residuals()
# if fh can be passed in predict
forecaster.fit(y).predict_residuals()
```

for out-of-sample residuals:
```python
# if fh must be passed in fit
forecaster.fit(y, fh=my_fh).predict_residuals(y_true)
# if fh can be passed in predict
forecaster.fit(y).predict_residuals(y_true)
```

The `fh` is inferred in either case, it's the inde of `y` in `fit` (in-sample) or `y_true` (out-of-sample); if no `y` is passed, the `y` in `fit` is used as the argument in `predict_residuals`.



@user1:
Hi @author Sorry for the delayed response. I was out over the holidays. I think your alternate solution is cleaner and we can work with this. 

@author:
alright!

@author:
I'll write tests then so it's ready to go in.

@author:
(done)

diff_hunk: @@ -631,6 +637,48 @@ def update_predict_single(
             alpha=alpha,
         )
 
+    def predict_residuals(self, y, X=None, fh=None):
+        """Return residuals of time series forecasts.
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_".
+            self.cutoff, self._is_fitted
+
+        Writes to self:
+            Stores fh to self.fh if fh is passed and has not been passed previously.
+
+        Parameters
+        ----------
+        y : pd.Series, pd.DataFrame, or np.ndarray (1D or 2D)
+            Time series with ground truth observations, to compute residuals to.
+            Must have same type, dimension, and indices as expected return of predict.
+        X : pd.DataFrame, or 2D np.ndarray, optional (default=None)
+            Exogeneous time series to predict from
+            if self.get_tag("X-y-must-have-same-index"), X.index must contain fh.index
+        fh : int, list, np.ndarray or ForecastingHorizon
+            Forecasting horizon
+
+        Returns
+        -------
+        y_res : pd.Series, pd.DataFrame, or np.ndarray (1D or 2D)
+            Forecast residuals at fh, with same index as fh
+            y_pred has same type as y passed in fit (most recently)
+        """
+        y_pred = self.predict(fh=fh, X=X)
path: sktime/forecasting/base/_base.py
text: @user1:
@author I do not see the handling for `fh=None` (i.e. return all in-sample residuals). Is that something that will be added later?

@author:
we could do this in this PR if you'd prefer that? I wouldn't mind.
Also not sure what the other core devs think what the default behaviour for `fh=None` should be.

@author:
ok, done - see below (if no argument is passed)