def residuals(
        self,
        series: Union[TimeSeries, Sequence[TimeSeries]],
        past_covariates: Optional[Union[TimeSeries, Sequence[TimeSeries]]] = None,
        future_covariates: Optional[Union[TimeSeries, Sequence[TimeSeries]]] = None,
        historical_forecasts: Optional[
            Union[TimeSeries, Sequence[TimeSeries], Sequence[Sequence[TimeSeries]]]
        ] = None,
        forecast_horizon: int = 1,
        num_samples: int = 1,
        train_length: Optional[int] = None,
        start: Optional[Union[pd.Timestamp, float, int]] = None,
        start_format: Literal["position", "value"] = "value",
        stride: int = 1,
        retrain: Union[bool, int, Callable[..., bool]] = True,
        overlap_end: bool = False,
        last_points_only: bool = True,
        metric: METRIC_TYPE = metrics.err,
        verbose: bool = False,
        show_warnings: bool = True,
        predict_likelihood_parameters: bool = False,
        enable_optimization: bool = True,
        data_transformers: Optional[
            dict[str, Union[BaseDataTransformer, Pipeline]]
        ] = None,
        metric_kwargs: Optional[dict[str, Any]] = None,
        fit_kwargs: Optional[dict[str, Any]] = None,
        predict_kwargs: Optional[dict[str, Any]] = None,
        sample_weight: Optional[Union[TimeSeries, Sequence[TimeSeries], str]] = None,
        values_only: bool = False,
        random_state: Optional[int] = None,
    ) -> Union[TimeSeries, list[TimeSeries], list[list[TimeSeries]]]:
        """Compute the residuals that the model produced for historical forecasts on (potentially multiple) `series`.

        This function computes the difference (or one of Darts' "per time step" metrics) between the actual
        observations from `series` and the fitted values obtained by training the model on `series` (or using a
        pre-trained model with `retrain=False`). Not all models support fitted values, so we use historical forecasts
        as an approximation for them.

        In sequence this method performs:

        - use pre-computed `historical_forecasts` or compute historical forecasts for each series (see
          :meth:`~darts.models.forecasting.forecasting_model.ForecastingModel.historical_forecasts` for more details).
          How the historical forecasts are generated can be configured with parameters `num_samples`, `train_length`,
          `start`, `start_format`, `forecast_horizon`, `stride`, `retrain`, `last_points_only`, `fit_kwargs`, and
          `predict_kwargs`.
        - compute a backtest using a "per time step" `metric` between the historical forecasts and `series` per
          component/column and time step (see
          :meth:`~darts.models.forecasting.forecasting_model.ForecastingModel.backtest` for more details). By default,
          uses the residuals :func:`~darts.metrics.metrics.err` (error) as a `metric`.
        - create and return `TimeSeries` (or simply a np.ndarray with `values_only=True`) with the time index from
          historical forecasts, and values from the metrics per component and time step.

        This method works for single or multiple univariate or multivariate series.
        It uses the median prediction (when dealing with stochastic forecasts).

        Parameters
        ----------
        series
            A (sequence of) target time series used to successively train (if `retrain` is not ``False``) and compute
            the historical forecasts.
        past_covariates
            Optionally, a (sequence of) past-observed covariate time series for every input time series in `series`.
            This applies only if the model supports past covariates.
        future_covariates
            Optionally, a (sequence of) future-known covariate time series for every input time series in `series`.
            This applies only if the model supports future covariates.
        historical_forecasts
            Optionally, the (or a sequence of / a sequence of sequences of) historical forecasts time series to be
            evaluated. Corresponds to the output of :meth:`historical_forecasts()
            <darts.models.forecasting.forecasting_model.ForecastingModel.historical_forecasts>`. The same `series` and
            `last_points_only` values must be passed that were used to generate the historical forecasts. If provided,
            will skip historical forecasting and ignore all parameters except `series`, `last_points_only`, `metric`,
            and `reduction`.
        forecast_horizon
            The forecast horizon for the predictions.
        num_samples
            Number of times a prediction is sampled from a probabilistic model. Use values ``>1`` only for probabilistic
            models.
        train_length
            Optionally, use a fixed length / number of time steps for every constructed training set (rolling window
            mode). Only effective when `retrain` is not ``False``. The default is ``None``, where it uses all time
            steps up until the prediction time (expanding window mode). If larger than the number of available time
            steps, uses the expanding mode. Needs to be at least `min_train_series_length`.
        start
            Optionally, the first point in time at which a prediction is computed. This parameter supports:
            ``float``, ``int``, ``pandas.Timestamp``, and ``None``.
            If a ``float``, it is the proportion of the time series that should lie before the first prediction point.
            If an ``int``, it is either the index position of the first prediction point for `series` with a
            `pd.DatetimeIndex`, or the index value for `series` with a `pd.RangeIndex`. The latter can be changed to
            the index position with `start_format="position"`.
            If a ``pandas.Timestamp``, it is the time stamp of the first prediction point.
            If ``None``, the first prediction point will automatically be set to:

            - the first predictable point if `retrain` is ``False``, or `retrain` is a Callable and the first
              predictable point is earlier than the first trainable point.
            - the first trainable point if `retrain` is ``True`` or ``int`` (given `train_length`),
              or `retrain` is a ``Callable`` and the first trainable point is earlier than the first predictable point.
            - the first trainable point (given `train_length`) otherwise

            Note: If `start` is not within the trainable / forecastable points, uses the closest valid start point that
              is a round multiple of `stride` ahead of `start`. Raises a `ValueError`, if no valid start point exists.
            Note: If the model uses a shifted output (`output_chunk_shift > 0`), then the first predicted point is also
              shifted by `output_chunk_shift` points into the future.
            Note: If `start` is outside the possible historical forecasting times, will ignore the parameter
              (default behavior with ``None``) and start at the first trainable/predictable point.
        start_format
            Defines the `start` format.
            If set to ``'position'``, `start` corresponds to the index position of the first predicted point and can
            range from `(-len(series), len(series) - 1)`.
            If set to ``'value'``, `start` corresponds to the index value/label of the first predicted point. Will raise
            an error if the value is not in `series`' index. Default: ``'value'``.
        stride
            The number of time steps between two consecutive predictions.
        retrain
            Whether and/or on which condition to retrain the model before predicting.
            This parameter supports 3 different types: ``bool``, (positive) ``int``, and ``Callable`` (returning a
            ``bool``).
            In the case of ``bool``: retrain the model at each step (`True`), or never retrain the model (`False`).
            In the case of ``int``: the model is retrained every `retrain` iterations.
            In the case of ``Callable``: the model is retrained whenever callable returns `True`.
            The callable must have the following positional arguments:

            - `counter` (int): current `retrain` iteration
            - `pred_time` (pd.Timestamp or int): timestamp of forecast time (end of the training series)
            - `train_series` (TimeSeries): train series up to `pred_time`
            - `past_covariates` (TimeSeries): past_covariates series up to `pred_time`
            - `future_covariates` (TimeSeries): future_covariates series up to `min(pred_time + series.freq *
              forecast_horizon, series.end_time())`

            Note: if any optional `*_covariates` are not passed to `historical_forecast`, ``None`` will be passed
            to the corresponding retrain function argument.
            Note: some models require being retrained every time and do not support anything other than
            `retrain=True`.
            Note: also controls the retraining of the `data_transformers`.
        overlap_end
            Whether the returned forecasts can go beyond the series' end or not.
        last_points_only
            Whether to return only the last point of each historical forecast. If set to ``True``, the method returns a
            single ``TimeSeries`` (for each time series in `series`) containing the successive point forecasts.
            Otherwise, returns a list of historical ``TimeSeries`` forecasts.
        metric
            Either one of Darts' "per time step" metrics (see `here
            <https://unit8co.github.io/darts/generated_api/darts.metrics.html>`_), or a custom metric that has an
            identical signature as Darts' "per time step" metrics, uses decorators
            :func:`~darts.metrics.metrics.multi_ts_support` and :func:`~darts.metrics.metrics.multi_ts_support`,
            and returns one value per time step.
        verbose
            Whether to print the progress.
        show_warnings
            Whether to show warnings related to historical forecasts optimization, or parameters `start` and
            `train_length`.
        predict_likelihood_parameters
            If set to `True`, the model predicts the parameters of its `likelihood` instead of the target. Only
            supported for probabilistic models with a likelihood, `num_samples = 1` and `n<=output_chunk_length`.
            Default: ``False``.
        enable_optimization
            Whether to use the optimized version of `historical_forecasts` when supported and available.
            Default: ``True``.
        data_transformers
            Optionally, a dictionary of `BaseDataTransformer` or `Pipeline` to apply to the corresponding series
            (possibles keys; "series", "past_covariates", "future_covariates"). If provided, all input series must be
            in the un-transformed space. For fittable transformer / pipeline:

            - if `retrain=True`, the data transformer re-fit on the training data at each historical forecast step.
            - if `retrain=False`, the data transformer transforms the series once before all the forecasts.

            The fitted transformer is used to transform the input during both training and prediction.
            If the transformation is invertible, the forecasts will be inverse-transformed.
            Only effective when `historical_forecasts=None`.
        metric_kwargs
            Additional arguments passed to `metric()`, such as `'n_jobs'` for parallelization, `'m'` for scaled
            metrics, etc. Will pass arguments only if they are present in the corresponding metric signature. Ignores
            reduction arguments `"series_reduction", "component_reduction", "time_reduction"`, and parameter
            `'insample'` for scaled metrics (e.g. mase`, `rmsse`, ...), as they are handled internally.
        fit_kwargs
            Optionally, some additional arguments passed to the model `fit()` method.
        predict_kwargs
            Optionally, some additional arguments passed to the model `predict()` method.
        sample_weight
            Optionally, some sample weights to apply to the target `series` labels for training. Only effective when
            `retrain` is not ``False``. They are applied per observation, per label (each step in
            `output_chunk_length`), and per component.
            If a series or sequence of series, then those weights are used. If the weight series only have a single
            component / column, then the weights are applied globally to all components in `series`. Otherwise, for
            component-specific weights, the number of components must match those of `series`.
            If a string, then the weights are generated using built-in weighting functions. The available options are
            `"linear"` or `"exponential"` decay - the further in the past, the lower the weight. The weights are
            computed per time `series`.
        values_only
            Whether to return the residuals as `np.ndarray`. If `False`, returns residuals as `TimeSeries`.
        random_state
            Controls the randomness of probabilistic predictions.

        Returns
        -------
        TimeSeries
            Residual `TimeSeries` for a single `series` and `historical_forecasts` generated with
            `last_points_only=True`.
        List[TimeSeries]
            A list of residual `TimeSeries` for a sequence (list) of `series` with `last_points_only=True`.
            The residual list has length `len(series)`.
        List[List[TimeSeries]]
            A list of lists of residual `TimeSeries` for a sequence of `series` with `last_points_only=False`.
            The outer residual list has length `len(series)`. The inner lists consist of the residuals from
            all possible series-specific historical forecasts.
        """
        # `residuals()` should return metrics per series, component and time step (no reduction)
        metric_kwargs = copy.deepcopy(metric_kwargs) or {}
        metric_kwargs["series_reduction"] = None
        metric_kwargs["component_reduction"] = None
        metric_kwargs["time_reduction"] = None

        historical_forecasts = historical_forecasts or self.historical_forecasts(
            series=series,
            past_covariates=past_covariates,
            future_covariates=future_covariates,
            num_samples=num_samples,
            train_length=train_length,
            start=start,
            start_format=start_format,
            forecast_horizon=forecast_horizon,
            stride=stride,
            retrain=retrain,
            last_points_only=last_points_only,
            verbose=verbose,
            show_warnings=show_warnings,
            predict_likelihood_parameters=predict_likelihood_parameters,
            enable_optimization=enable_optimization,
            data_transformers=data_transformers,
            fit_kwargs=fit_kwargs,
            predict_kwargs=predict_kwargs,
            overlap_end=overlap_end,
            sample_weight=sample_weight,
            random_state=random_state,
        )

        # remember input series type
        series_seq_type = get_series_seq_type(series)
        # validate historical forecasts and convert to multiple series with multiple forecasts case
        series, historical_forecasts = _process_historical_forecast_for_backtest(
            series=series,
            historical_forecasts=historical_forecasts,
            last_points_only=last_points_only,
        )

        # optionally, add nans to end of series to get residuals of same shape for each forecast
        if overlap_end:
            series = _extend_series_for_overlap_end(
                series=series, historical_forecasts=historical_forecasts
            )

        residuals = self.backtest(
            series=series,
            historical_forecasts=historical_forecasts,
            last_points_only=False,
            metric=metric,
            reduction=None,
            data_transformers=data_transformers,
            metric_kwargs=metric_kwargs,
            random_state=random_state,
        )

        # sanity check residual output
        q, q_interval = metric_kwargs.get("q"), metric_kwargs.get("q_interval")
        try:
            series_, res, fc = series[0], residuals[0][0], historical_forecasts[0][0]
            _ = np.reshape(res, (len(fc), -1, 1))
        except Exception as err:
            raise_log(
                ValueError(
                    f"`metric` function did not yield expected output. Make sure "
                    f"to use one of Darts 'per time step' metrics, or a similar "
                    f"custom metric. The following exception was raised: "
                    f"{type(err).__name__}('{err}')"
                ),
                logger=logger,
            )

        # process residuals
        residuals_out = []
        for series_, fc_list, res_list in zip(series, historical_forecasts, residuals):
            res_list_out = []
            if q is not None:
                q = [q] if isinstance(q, float) else q
                # multi-quantile metrics yield more components
                comp_names = likelihood_component_names(
                    components=series_.components,
                    parameter_names=quantile_names(q),
                )
            # `q` and `q_interval` are mutually exclusive
            elif q_interval is not None:
                # multi-quantile metrics yield more components
                q_interval = (
                    [q_interval] if isinstance(q_interval, tuple) else q_interval
                )
                comp_names = likelihood_component_names(
                    components=series_.components,
                    parameter_names=quantile_interval_names(q_interval),
                )
            else:
                comp_names = None
            for fc, res in zip(fc_list, res_list):
                # make sure all residuals have shape (n time steps, n components * n quantiles, n samples=1)
                if len(res.shape) != 3:
                    res = np.reshape(res, (len(fc), -1, 1))
                if values_only:
                    res = res
                elif (q is None and q_interval is None) and res.shape[
                    1
                ] == fc.n_components:
                    res = fc.with_values(res)
                else:
                    # quantile (interval) metrics created different number of components;
                    # create new series with unknown components
                    res = TimeSeries.from_times_and_values(
                        times=fc._time_index,
                        values=res,
                        columns=comp_names,
                    )
                res_list_out.append(res)

            residuals_out.append(res_list_out)

        # if required, reduce to `series` input type
        if series_seq_type == SeriesType.SINGLE:
            return residuals_out[0][0] if last_points_only else residuals_out[0]

        return (
            [res for res_list in residuals_out for res in res_list]
            if last_points_only
            else residuals_out
        )
Compute the residuals that the model produced for historical forecasts on (potentially multiple) `series`.

        This function computes the difference (or one of Darts' "per time step" metrics) between the actual
        observations from `series` and the fitted values obtained by training the model on `series` (or using a
        pre-trained model with `retrain=False`). Not all models support fitted values, so we use historical forecasts
        as an approximation for them.

        In sequence this method performs:

        - use pre-computed `historical_forecasts` or compute historical forecasts for each series (see
          :meth:`~darts.models.forecasting.forecasting_model.ForecastingModel.historical_forecasts` for more details).
          How the historical forecasts are generated can be configured with parameters `num_samples`, `train_length`,
          `start`, `start_format`, `forecast_horizon`, `stride`, `retrain`, `last_points_only`, `fit_kwargs`, and
          `predict_kwargs`.
        - compute a backtest using a "per time step" `metric` between the historical forecasts and `series` per
          component/column and time step (see
          :meth:`~darts.models.forecasting.forecasting_model.ForecastingModel.backtest` for more details). By default,
          uses the residuals :func:`~darts.metrics.metrics.err` (error) as a `metric`.
        - create and return `TimeSeries` (or simply a np.ndarray with `values_only=True`) with the time index from
          historical forecasts, and values from the metrics per component and time step.

        This method works for single or multiple univariate or multivariate series.
        It uses the median prediction (when dealing with stochastic forecasts).

        Parameters
        ----------
        series
            A (sequence of) target time series used to successively train (if `retrain` is not ``False``) and compute
            the historical forecasts.
        past_covariates
            Optionally, a (sequence of) past-observed covariate time series for every input time series in `series`.
            This applies only if the model supports past covariates.
        future_covariates
            Optionally, a (sequence of) future-known covariate time series for every input time series in `series`.
            This applies only if the model supports future covariates.
        historical_forecasts
            Optionally, the (or a sequence of / a sequence of sequences of) historical forecasts time series to be
            evaluated. Corresponds to the output of :meth:`historical_forecasts()
            <darts.models.forecasting.forecasting_model.ForecastingModel.historical_forecasts>`. The same `series` and
            `last_points_only` values must be passed that were used to generate the historical forecasts. If provided,
            will skip historical forecasting and ignore all parameters except `series`, `last_points_only`, `metric`,
            and `reduction`.
        forecast_horizon
            The forecast horizon for the predictions.
        num_samples
            Number of times a prediction is sampled from a probabilistic model. Use values ``>1`` only for probabilistic
            models.
        train_length
            Optionally, use a fixed length / number of time steps for every constructed training set (rolling window
            mode). Only effective when `retrain` is not ``False``. The default is ``None``, where it uses all time
            steps up until the prediction time (expanding window mode). If larger than the number of available time
            steps, uses the expanding mode. Needs to be at least `min_train_series_length`.
        start
            Optionally, the first point in time at which a prediction is computed. This parameter supports:
            ``float``, ``int``, ``pandas.Timestamp``, and ``None``.
            If a ``float``, it is the proportion of the time series that should lie before the first prediction point.
            If an ``int``, it is either the index position of the first prediction point for `series` with a
            `pd.DatetimeIndex`, or the index value for `series` with a `pd.RangeIndex`. The latter can be changed to
            the index position with `start_format="position"`.
            If a ``pandas.Timestamp``, it is the time stamp of the first prediction point.
            If ``None``, the first prediction point will automatically be set to:

            - the first predictable point if `retrain` is ``False``, or `retrain` is a Callable and the first
              predictable point is earlier than the first trainable point.
            - the first trainable point if `retrain` is ``True`` or ``int`` (given `train_length`),
              or `retrain` is a ``Callable`` and the first trainable point is earlier than the first predictable point.
            - the first trainable point (given `train_length`) otherwise

            Note: If `start` is not within the trainable / forecastable points, uses the closest valid start point that
              is a round multiple of `stride` ahead of `start`. Raises a `ValueError`, if no valid start point exists.
            Note: If the model uses a shifted output (`output_chunk_shift > 0`), then the first predicted point is also
              shifted by `output_chunk_shift` points into the future.
            Note: If `start` is outside the possible historical forecasting times, will ignore the parameter
              (default behavior with ``None``) and start at the first trainable/predictable point.
        start_format
            Defines the `start` format.
            If set to ``'position'``, `start` corresponds to the index position of the first predicted point and can
            range from `(-len(series), len(series) - 1)`.
            If set to ``'value'``, `start` corresponds to the index value/label of the first predicted point. Will raise
            an error if the value is not in `series`' index. Default: ``'value'``.
        stride
            The number of time steps between two consecutive predictions.
        retrain
            Whether and/or on which condition to retrain the model before predicting.
            This parameter supports 3 different types: ``bool``, (positive) ``int``, and ``Callable`` (returning a
            ``bool``).
            In the case of ``bool``: retrain the model at each step (`True`), or never retrain the model (`False`).
            In the case of ``int``: the model is retrained every `retrain` iterations.
            In the case of ``Callable``: the model is retrained whenever callable returns `True`.
            The callable must have the following positional arguments:

            - `counter` (int): current `retrain` iteration
            - `pred_time` (pd.Timestamp or int): timestamp of forecast time (end of the training series)
            - `train_series` (TimeSeries): train series up to `pred_time`
            - `past_covariates` (TimeSeries): past_covariates series up to `pred_time`
            - `future_covariates` (TimeSeries): future_covariates series up to `min(pred_time + series.freq *
              forecast_horizon, series.end_time())`

            Note: if any optional `*_covariates` are not passed to `historical_forecast`, ``None`` will be passed
            to the corresponding retrain function argument.
            Note: some models require being retrained every time and do not support anything other than
            `retrain=True`.
            Note: also controls the retraining of the `data_transformers`.
        overlap_end
            Whether the returned forecasts can go beyond the series' end or not.
        last_points_only
            Whether to return only the last point of each historical forecast. If set to ``True``, the method returns a
            single ``TimeSeries`` (for each time series in `series`) containing the successive point forecasts.
            Otherwise, returns a list of historical ``TimeSeries`` forecasts.
        metric
            Either one of Darts' "per time step" metrics (see `here
            <https://unit8co.github.io/darts/generated_api/darts.metrics.html>`_), or a custom metric that has an
            identical signature as Darts' "per time step" metrics, uses decorators
            :func:`~darts.metrics.metrics.multi_ts_support` and :func:`~darts.metrics.metrics.multi_ts_support`,
            and returns one value per time step.
        verbose
            Whether to print the progress.
        show_warnings
            Whether to show warnings related to historical forecasts optimization, or parameters `start` and
            `train_length`.
        predict_likelihood_parameters
            If set to `True`, the model predicts the parameters of its `likelihood` instead of the target. Only
            supported for probabilistic models with a likelihood, `num_samples = 1` and `n<=output_chunk_length`.
            Default: ``False``.
        enable_optimization
            Whether to use the optimized version of `historical_forecasts` when supported and available.
            Default: ``True``.
        data_transformers
            Optionally, a dictionary of `BaseDataTransformer` or `Pipeline` to apply to the corresponding series
            (possibles keys; "series", "past_covariates", "future_covariates"). If provided, all input series must be
            in the un-transformed space. For fittable transformer / pipeline:

            - if `retrain=True`, the data transformer re-fit on the training data at each historical forecast step.
            - if `retrain=False`, the data transformer transforms the series once before all the forecasts.

            The fitted transformer is used to transform the input during both training and prediction.
            If the transformation is invertible, the forecasts will be inverse-transformed.
            Only effective when `historical_forecasts=None`.
        metric_kwargs
            Additional arguments passed to `metric()`, such as `'n_jobs'` for parallelization, `'m'` for scaled
            metrics, etc. Will pass arguments only if they are present in the corresponding metric signature. Ignores
            reduction arguments `"series_reduction", "component_reduction", "time_reduction"`, and parameter
            `'insample'` for scaled metrics (e.g. mase`, `rmsse`, ...), as they are handled internally.
        fit_kwargs
            Optionally, some additional arguments passed to the model `fit()` method.
        predict_kwargs
            Optionally, some additional arguments passed to the model `predict()` method.
        sample_weight
            Optionally, some sample weights to apply to the target `series` labels for training. Only effective when
            `retrain` is not ``False``. They are applied per observation, per label (each step in
            `output_chunk_length`), and per component.
            If a series or sequence of series, then those weights are used. If the weight series only have a single
            component / column, then the weights are applied globally to all components in `series`. Otherwise, for
            component-specific weights, the number of components must match those of `series`.
            If a string, then the weights are generated using built-in weighting functions. The available options are
            `"linear"` or `"exponential"` decay - the further in the past, the lower the weight. The weights are
            computed per time `series`.
        values_only
            Whether to return the residuals as `np.ndarray`. If `False`, returns residuals as `TimeSeries`.
        random_state
            Controls the randomness of probabilistic predictions.

        Returns
        -------
        TimeSeries
            Residual `TimeSeries` for a single `series` and `historical_forecasts` generated with
            `last_points_only=True`.
        List[TimeSeries]
            A list of residual `TimeSeries` for a sequence (list) of `series` with `last_points_only=True`.
            The residual list has length `len(series)`.
        List[List[TimeSeries]]
            A list of lists of residual `TimeSeries` for a sequence of `series` with `last_points_only=False`.
            The outer residual list has length `len(series)`. The inner lists consist of the residuals from
            all possible series-specific historical forecasts.
        
darts/models/forecasting/forecasting_model.py