Cross-correlation of multidimensional time-series and Python function for computing it So far I know, the cross-correlation of two time-series $a(t)$ and $b(t)$ for which $N$ observations are available is given by $r^N(\tau)=\frac{1}{N} \sum_{t=\tau+1}^Na(t-\tau)b(t)$, where $\tau$ indicates the time-lag.
However, I am not sure how to cross-correlate two multi-dimensional time-series.
For instance, let $x(t)$ and $y(t)$ two multidimensional time-series of dimension $q$ and $p$, respectively, and assume that I have $N$ observations that I arranged in a matrix form as it follows
\begin{equation}
        x^{N}:=
    \begin{bmatrix}
        x_1(1) & \dots &  x_q(1) \\
        \vdots & \ddots & \vdots \\
        x_1(N) & \dots &  x_q(N)
    \end{bmatrix}\,,
\end{equation}
and
\begin{equation}
        y^{N}:=
    \begin{bmatrix}
        y_1(1) & \dots &  y_p(1) \\
        \vdots & \ddots & \vdots \\
        y_1(N) & \dots &  y_p(N)
    \end{bmatrix}\,.
\end{equation}
My instinct would suggest that the cross-correlation $R_{xy}(\tau)$ between $x^N$ and $y^N$  is given by
\begin{equation}
        R_{xy}(\tau)=
        \begin{bmatrix}
            \frac{1}{N}\sum_{t=\tau+1}^N x_1(t-\tau)y_1(t) & \frac{1}{N}\sum_{t=\tau+1}^N x_1(t-\tau)y_2(t) & \cdots & \frac{1}{N}\sum_{t=\tau+1}^N x_1(t-\tau)y_p(t)] \\
            \frac{1}{N}\sum_{t=\tau+1}^N x_2(t-\tau)y_1(t) & \frac{1}{N}\sum_{t=\tau+1}^N x_2(t-\tau)y_2(t) & \cdots & \frac{1}{N}\sum_{t=\tau+1}^N x_2(t-\tau)y_p(t)] \\
            \vdots & \vdots & \ddots & \vdots \\
            \frac{1}{N}\sum_{t=\tau+1}^N x_q(t-\tau)y_1(t) & \frac{1}{N}\sum_{t=\tau+1}^N x_q(t-\tau)y_2(t) & \cdots & \frac{1}{N}\sum_{t=\tau+1}^N x_q(t-\tau)y_p(t)]
        \end{bmatrix}\,.
    \end{equation}
Is it correct?
Also, is there any Python function that would allow me to compute the cross-correlation between $x^N$ and $y^N$? If so, what's that function name?

The summation indices should be from $\tau+1$ to $N$, assuming 1-indexed series, and the cross-correlation matrix looks fine. At least in statsmodels library, which is a popular tool for time series workloads, there isn't any function to calculate a sample cross-correlation matrix with $\tau=[0,N-1]$ entries, which would be a 3D tensor.
