def qr(
    actual_series: Union[TimeSeries, Sequence[TimeSeries]],
    pred_series: Union[TimeSeries, Sequence[TimeSeries]],
    intersect: bool = True,
    *,
    q: Union[float, list[float], tuple[np.ndarray, pd.Index]] = 0.5,
    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:
    """Quantile Risk (QR)

    QR is a metric that quantifies the accuracy of a specific quantile :math:`q` from the predicted value
    distribution of a stochastic/probabilistic `pred_series` containing N samples.

    The main difference to the Quantile Loss (QL) is that QR computes the quantile and loss on the aggregate of all
    sample values summed up along the time axis (QL computes the quantile and loss per time step).

    For the true series :math:`y` and predicted stochastic/probabilistic series (containing N samples) :math:`\\hat{y}`
    of of shape :math:`T \\times N`, it is computed per column/component and quantile as:

    .. math:: 2 \\frac{QL(Z, \\hat{Z}_q)}{Z},

    where :math:`QL` is the Quantile Loss (:func:`~darts.metrics.metrics.ql`), :math:`Z = \\sum_{t=1}^{T} y_t` is
    the sum of all target/actual values, :math:`\\hat{Z} = \\sum_{t=1}^{T} \\hat{y}_t` is the sum of all predicted
    samples along the time axis, and :math:`\\hat{Z}_q` is the quantile :math:`q` of that sum.

    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
        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.
    """
    if not pred_series.is_stochastic:
        raise_log(
            ValueError(
                "quantile risk (qr) should only be computed for stochastic predicted TimeSeries."
            ),
            logger=logger,
        )

    z_true, z_hat = _get_values_or_raise(
        actual_series,
        pred_series,
        intersect,
        q=None,
        remove_nan_union=True,
    )
    z_true = np.nansum(z_true, axis=TIME_AX)
    z_hat = np.nansum(
        z_hat, axis=TIME_AX
    )  # aggregate all individual sample realizations
    # quantile loss
    q, _ = q
    z_hat_rho = np.quantile(
        z_hat, q=q, axis=1
    ).T  # get the quantile from aggregated samples

    errors = z_true - z_hat_rho
    losses = 2 * np.maximum((q - 1) * errors, q * errors)
    return losses / z_true
Quantile Risk (QR)

    QR is a metric that quantifies the accuracy of a specific quantile :math:`q` from the predicted value
    distribution of a stochastic/probabilistic `pred_series` containing N samples.

    The main difference to the Quantile Loss (QL) is that QR computes the quantile and loss on the aggregate of all
    sample values summed up along the time axis (QL computes the quantile and loss per time step).

    For the true series :math:`y` and predicted stochastic/probabilistic series (containing N samples) :math:`\hat{y}`
    of of shape :math:`T \times N`, it is computed per column/component and quantile as:

    .. math:: 2 \frac{QL(Z, \hat{Z}_q)}{Z},

    where :math:`QL` is the Quantile Loss (:func:`~darts.metrics.metrics.ql`), :math:`Z = \sum_{t=1}^{T} y_t` is
    the sum of all target/actual values, :math:`\hat{Z} = \sum_{t=1}^{T} \hat{y}_t` is the sum of all predicted
    samples along the time axis, and :math:`\hat{Z}_q` is the quantile :math:`q` of that sum.

    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
        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.
    
darts/metrics/metrics.py