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?