def temporal_autocorrelation_multivariate(
    x, d=0, mask=None, window="gaussian", window_radius=np.inf
):
    """
    For a :math:`q`-variate time series
    :math:`\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_n`, compute the lag-l
    correlation matrices :math:`\mathbf{\Gamma}_l`, where
    :math:`\Gamma_{l,i,j}=\gamma_{l,i,j}` and
    :math:`\gamma_{l,i,j}=\mbox{corr}(x_i(t),x_j(t-l))` for
    :math:`i,j=1,2,\dots,q` and :math:`l=0,1,\dots,n-1`.

    Parameters
    ----------
    x: array_like
        Array of shape (n, q, ...) containing the time series :math:`\mathbf{x}_i`.
        The inputs are assumed to be in increasing order with respect to time,
        and the time step is assumed to be regular. All inputs are required to
        have finite values. The remaining dimensions after the second one are
        flattened before computing the correlation coefficients.
    d: {0,1}
        The order of differencing. If d=1, the input time series is differenced
        before computing the correlation coefficients. In this case, a time
        series of length n+1 is needed for computing the n-1 coefficients.
    mask: array_like
        Optional mask to use for computing the correlation coefficients. Input
        elements with mask==False are excluded from the computations. The shape
        of the mask is expected to be x.shape[2:].
    window: {"gaussian", "uniform"}
        The weight function to use for the moving window. Applicable if
        window_radius < np.inf. Defaults to 'gaussian'.
    window_radius: float
        If window_radius < np.inf, the correlation coefficients are computed in
        a moving window. Defaults to np.inf (i.e. the correlations are computed
        over the whole domain). If window is 'gaussian', window_radius is the
        standard deviation of the Gaussian filter. If window is 'uniform', the
        size of the window is 2*window_radius+1.

    Returns
    -------
    out: list
        List of correlation matrices :math:`\Gamma_0,\Gamma_1,\dots,\Gamma_{n-1}`
        of shape (q,q). If window_radius<np.inf, the elements of the list are
        arrays of shape (x.shape[2:],q,q). In this case, nan values are assigned,
        when the sample size for computing the correlation coefficients is too
        small.

    References
    ----------
    :cite:`CP2002`

    """
    if len(x.shape) < 3:
        raise ValueError("the dimension of x must be >= 3")
    if mask is not None and mask.shape != x.shape[2:]:
        raise ValueError(
            "dimension mismatch between x and mask: x.shape[2:]=%s, mask.shape=%s"
            % (str(x.shape[2:]), str(mask.shape))
        )
    if np.any(~np.isfinite(x)):
        raise ValueError("x contains non-finite values")

    if d == 1:
        x = np.diff(x, axis=0)

    p = x.shape[0] - 1
    q = x.shape[1]

    gamma = []
    for k in range(p + 1):
        if window_radius == np.inf:
            gamma_k = np.empty((q, q))
        else:
            gamma_k = np.empty(np.hstack([x.shape[2:], [q, q]]))
        for i in range(q):
            x_i = x[-1, i, :]
            for j in range(q):
                x_j = x[-(k + 1), j, :]
                if window_radius == np.inf:
                    gamma_k[i, j] = np.corrcoef(x_i.flatten(), x_j.flatten())[0, 1]
                else:
                    gamma_k[..., i, j] = _moving_window_corrcoef(
                        x_i, x_j, window_radius, window=window, mask=mask
                    )

        gamma.append(gamma_k)

    return gamma
For a :math:`q`-variate time series
:math:`\mathbf{x}_1,\mathbf{x}_2,\dots,\mathbf{x}_n`, compute the lag-l
correlation matrices :math:`\mathbf{\Gamma}_l`, where
:math:`\Gamma_{l,i,j}=\gamma_{l,i,j}` and
:math:`\gamma_{l,i,j}=\mbox{corr}(x_i(t),x_j(t-l))` for
:math:`i,j=1,2,\dots,q` and :math:`l=0,1,\dots,n-1`.

Parameters
----------
x: array_like
Array of shape (n, q, ...) containing the time series :math:`\mathbf{x}_i`.
The inputs are assumed to be in increasing order with respect to time,
and the time step is assumed to be regular. All inputs are required to
have finite values. The remaining dimensions after the second one are
flattened before computing the correlation coefficients.
d: {0,1}
The order of differencing. If d=1, the input time series is differenced
before computing the correlation coefficients. In this case, a time
series of length n+1 is needed for computing the n-1 coefficients.
mask: array_like
Optional mask to use for computing the correlation coefficients. Input
elements with mask==False are excluded from the computations. The shape
of the mask is expected to be x.shape[2:].
window: {"gaussian", "uniform"}
The weight function to use for the moving window. Applicable if
window_radius < np.inf. Defaults to 'gaussian'.
window_radius: float
If window_radius < np.inf, the correlation coefficients are computed in
a moving window. Defaults to np.inf (i.e. the correlations are computed
over the whole domain). If window is 'gaussian', window_radius is the
standard deviation of the Gaussian filter. If window is 'uniform', the
size of the window is 2*window_radius+1.

Returns
-------
out: list
List of correlation matrices :math:`\Gamma_0,\Gamma_1,\dots,\Gamma_{n-1}`
of shape (q,q). If window_radius<np.inf, the elements of the list are
arrays of shape (x.shape[2:],q,q). In this case, nan values are assigned,
when the sample size for computing the correlation coefficients is too
small.

References
----------
:cite:`CP2002`
pysteps/timeseries/correlation.py