def r2_score(
    actual_series: Union[TimeSeries, Sequence[TimeSeries]],
    pred_series: Union[TimeSeries, Sequence[TimeSeries]],
    intersect: bool = True,
    *,
    q: Optional[Union[float, list[float], tuple[np.ndarray, pd.Index]]] = None,
    component_reduction: Optional[Callable[[np.ndarray], float]] = np.nanmean,
    series_reduction: Optional[Callable[[np.ndarray], Union[float, np.ndarray]]] = None,
    n_jobs: int = 1,
    verbose: bool = False,
) -> METRIC_OUTPUT_TYPE:
    """Coefficient of Determination :math:`R^2` (see [1]_ for more details).

    For the true series :math:`y` and predicted series :math:`\\hat{y}` of length :math:`T`, it is computed per
    component/column and (optional) quantile as:

    .. math:: 1 - \\frac{\\sum_{t=1}^T{(y_t - \\hat{y}_t)^2}}{\\sum_{t=1}^T{(y_t - \\bar{y})^2}},

    where :math:`\\bar{y}` is the mean of :math:`y` over all time steps.

    This metric is not symmetric.

    If :math:`\\hat{y}_t` are stochastic (contains several samples) or quantile predictions, use parameter `q` to
    specify on which quantile(s) to compute the metric on. By default, it uses the median 0.5 quantile
    (over all samples, or, if given, the quantile prediction itself).

    Parameters
    ----------
    actual_series
        The (sequence of) actual series.
    pred_series
        The (sequence of) predicted series.
    intersect
        For time series that are overlapping in time without having the same time index, setting `True`
        will consider the values only over their common time interval (intersection in time).
    q
        Optionally, the quantile (float [0, 1]) or list of quantiles of interest to compute the metric on.
    component_reduction
        Optionally, a function to aggregate the metrics over the component/column axis. It must reduce a `np.ndarray`
        of shape `(t, c)` to a `np.ndarray` of shape `(t,)`. The function takes as input a ``np.ndarray`` and a
        parameter named `axis`, and returns the reduced array. The `axis` receives value `1` corresponding to the
        component axis. If `None`, will return a metric per component.
    series_reduction
        Optionally, a function to aggregate the metrics over multiple series. It must reduce a `np.ndarray`
        of shape `(s, t, c)` to a `np.ndarray` of shape `(t, c)` The function takes as input a ``np.ndarray`` and a
        parameter named `axis`, and returns the reduced array. The `axis` receives value `0` corresponding to the
        series axis. For example with `np.nanmean`, will return the average over all series metrics. If `None`, will
        return a metric per component.
    n_jobs
        The number of jobs to run in parallel. Parallel jobs are created only when a ``Sequence[TimeSeries]`` is
        passed as input, parallelising operations regarding different ``TimeSeries``. Defaults to `1`
        (sequential). Setting the parameter to `-1` means using all the available processors.
    verbose
        Optionally, whether to print operations progress

    Returns
    -------
    float
        A single metric score (when `len(q) <= 1`) for:

        - a single univariate series.
        - a single multivariate series with `component_reduction`.
        - a sequence (list) of uni/multivariate series with `series_reduction` and `component_reduction`.
    np.ndarray
        A numpy array of metric scores. The array has shape (n components * n quantiles,) without component reduction,
        and shape (n quantiles,) with component reduction and `len(q) > 1`.
        For:

        - the same input arguments that result in the `float` return case from above but with `len(q) > 1`.
        - a single multivariate series and at least `component_reduction=None`.
        - a sequence of uni/multivariate series including `series_reduction` and `component_reduction=None`.
    list[float]
        Same as for type `float` but for a sequence of series.
    list[np.ndarray]
        Same as for type `np.ndarray` but for a sequence of series.

    References
    ----------
    .. [1] https://en.wikipedia.org/wiki/Coefficient_of_determination
    """
    y_true, y_pred = _get_values_or_raise(
        actual_series,
        pred_series,
        intersect,
        remove_nan_union=True,
        q=q,
    )
    ss_errors = np.nansum((y_true - y_pred) ** 2, axis=TIME_AX)
    y_hat = np.nanmean(y_true, axis=TIME_AX)
    ss_tot = np.nansum((y_true - y_hat) ** 2, axis=TIME_AX)
    return 1 - ss_errors / ss_tot
Coefficient of Determination :math:`R^2` (see [1]_ for more details).

    For the true series :math:`y` and predicted series :math:`\hat{y}` of length :math:`T`, it is computed per
    component/column and (optional) quantile as:

    .. math:: 1 - \frac{\sum_{t=1}^T{(y_t - \hat{y}_t)^2}}{\sum_{t=1}^T{(y_t - \bar{y})^2}},

    where :math:`\bar{y}` is the mean of :math:`y` over all time steps.

    This metric is not symmetric.

    If :math:`\hat{y}_t` are stochastic (contains several samples) or quantile predictions, use parameter `q` to
    specify on which quantile(s) to compute the metric on. By default, it uses the median 0.5 quantile
    (over all samples, or, if given, the quantile prediction itself).

    Parameters
    ----------
    actual_series
        The (sequence of) actual series.
    pred_series
        The (sequence of) predicted series.
    intersect
        For time series that are overlapping in time without having the same time index, setting `True`
        will consider the values only over their common time interval (intersection in time).
    q
        Optionally, the quantile (float [0, 1]) or list of quantiles of interest to compute the metric on.
    component_reduction
        Optionally, a function to aggregate the metrics over the component/column axis. It must reduce a `np.ndarray`
        of shape `(t, c)` to a `np.ndarray` of shape `(t,)`. The function takes as input a ``np.ndarray`` and a
        parameter named `axis`, and returns the reduced array. The `axis` receives value `1` corresponding to the
        component axis. If `None`, will return a metric per component.
    series_reduction
        Optionally, a function to aggregate the metrics over multiple series. It must reduce a `np.ndarray`
        of shape `(s, t, c)` to a `np.ndarray` of shape `(t, c)` The function takes as input a ``np.ndarray`` and a
        parameter named `axis`, and returns the reduced array. The `axis` receives value `0` corresponding to the
        series axis. For example with `np.nanmean`, will return the average over all series metrics. If `None`, will
        return a metric per component.
    n_jobs
        The number of jobs to run in parallel. Parallel jobs are created only when a ``Sequence[TimeSeries]`` is
        passed as input, parallelising operations regarding different ``TimeSeries``. Defaults to `1`
        (sequential). Setting the parameter to `-1` means using all the available processors.
    verbose
        Optionally, whether to print operations progress

    Returns
    -------
    float
        A single metric score (when `len(q) <= 1`) for:

        - a single univariate series.
        - a single multivariate series with `component_reduction`.
        - a sequence (list) of uni/multivariate series with `series_reduction` and `component_reduction`.
    np.ndarray
        A numpy array of metric scores. The array has shape (n components * n quantiles,) without component reduction,
        and shape (n quantiles,) with component reduction and `len(q) > 1`.
        For:

        - the same input arguments that result in the `float` return case from above but with `len(q) > 1`.
        - a single multivariate series and at least `component_reduction=None`.
        - a sequence of uni/multivariate series including `series_reduction` and `component_reduction=None`.
    list[float]
        Same as for type `float` but for a sequence of series.
    list[np.ndarray]
        Same as for type `np.ndarray` but for a sequence of series.

    References
    ----------
    .. [1] https://en.wikipedia.org/wiki/Coefficient_of_determination
    
darts/metrics/metrics.py