def multivariate_support(func) -> Callable[..., METRIC_OUTPUT_TYPE]:
    """
    This decorator transforms a metric function that takes as input two univariate TimeSeries instances
    into a function that takes two equally-sized multivariate TimeSeries instances, computes the pairwise univariate
    metrics for components with the same indices, and returns a float value that is computed as a function of all the
    univariate metrics using a `component_reduction` subroutine passed as argument to the metric function.
    """

    @wraps(func)
    def wrapper_multivariate_support(*args, **kwargs) -> METRIC_OUTPUT_TYPE:
        params = signature(func).parameters
        # we can avoid checks about args and kwargs since the input is adjusted by the previous decorator
        actual_series = args[0]
        pred_series = args[1]
        num_series_in_args = 2

        q, q_comp_names = kwargs.get("q"), None
        if q is None:
            # without quantiles, the number of components must match
            if actual_series.n_components != pred_series.n_components:
                raise_log(
                    ValueError(
                        f"Mismatch between number of components in `actual_series` "
                        f"(n={actual_series.width}) and `pred_series` (n={pred_series.width})."
                    ),
                    logger=logger,
                )
            # compute median for stochastic predictions
            if pred_series.is_stochastic:
                q = np.array([0.5])
        else:
            # `q` is required to be a tuple (handled by `multi_ts_support` wrapper)
            if not isinstance(q, tuple) or not len(q) == 2:
                raise_log(
                    ValueError(
                        "`q` must be of tuple of `(np.ndarray, Optional[pd.Index])` "
                        "where the (quantile values, optional quantile component names). "
                        f"Received `q={q}`."
                    ),
                    logger=logger,
                )
            q, q_comp_names = q
            if not pred_series.is_stochastic:
                # quantile component names are required if the predictions are not stochastic (as for stochastic
                # predictions, the quantiles can be retrieved from the sample dimension for each component)
                if q_comp_names is None:
                    q_comp_names = pd.Index(
                        likelihood_component_names(
                            components=actual_series.components,
                            parameter_names=quantile_names(q=q),
                        )
                    )
                if not q_comp_names.isin(pred_series.components).all():
                    raise_log(
                        ValueError(
                            f"Computing a metric with quantile(s) `q={q}` is only supported for probabilistic "
                            f"`pred_series` (num samples > 1) or `pred_series` containing the predicted "
                            f"quantiles as columns / components. Either pass a probabilistic `pred_series` or "
                            f"a series containing the expected quantile components: {q_comp_names.tolist()} "
                        ),
                        logger=logger,
                    )

        if "q" in params:
            kwargs["q"] = (q, q_comp_names)

        # handle `insample` parameters for scaled metrics
        input_series = (actual_series, pred_series)
        if "insample" in params:
            insample = args[2]
            if actual_series.n_components != insample.n_components:
                raise_log(
                    ValueError(
                        f"Mismatch between number of components in `actual_series` "
                        f"(n={actual_series.width}) and `insample` (n={insample.width}."
                    ),
                    logger=logger,
                )
            input_series += (insample,)
            num_series_in_args += 1

        vals = func(*input_series, *args[num_series_in_args:], **kwargs)
        # bring vals to shape (n_time, n_comp, n_quantile)
        if not 2 <= len(vals.shape) <= 3:
            raise_log(
                ValueError(
                    "Metric output must have 2 dimensions (n components, n quantiles) "
                    "for aggregated metrics (e.g. `mae()`, ...), "
                    "or 3 dimension (n times, n components, n quantiles)  "
                    "for time dependent metrics (e.g. `ae()`, ...)"
                ),
                logger=logger,
            )
        if len(vals.shape) == 2:
            vals = np.expand_dims(vals, TIME_AX)

        time_reduction = _get_reduction(
            kwargs=kwargs,
            params=params,
            red_name="time_reduction",
            axis=TIME_AX,
            sanity_check=False,
        )
        if time_reduction is not None:
            # -> (1, n_comp, n_quantile)
            vals = np.expand_dims(time_reduction(vals, axis=TIME_AX), axis=TIME_AX)

        component_reduction = _get_reduction(
            kwargs=kwargs,
            params=params,
            red_name="component_reduction",
            axis=COMP_AX,
            sanity_check=False,
        )
        if component_reduction is not None:
            # -> (*, n_quantile)
            vals = component_reduction(vals, axis=COMP_AX)
        else:
            # -> (*, n_comp * n_quantile), with order [c0_q0, c0_q1, ... c1_q0, c1_q1, ...]
            vals = vals.reshape(vals.shape[0], -1)
        return vals

    return wrapper_multivariate_support

    This decorator transforms a metric function that takes as input two univariate TimeSeries instances
    into a function that takes two equally-sized multivariate TimeSeries instances, computes the pairwise univariate
    metrics for components with the same indices, and returns a float value that is computed as a function of all the
    univariate metrics using a `component_reduction` subroutine passed as argument to the metric function.
    
darts/metrics/metrics.py