Assume the 2x1 random vector, $E_t$, has a multivariate normal distribution with a mean vector of 0 and covariance matrix of $Q_t$. Note here that $E_t$ is a time series variable and the covariance matrix is differs across time
$$E_t = [e_{1,t}, e_{2,t}]^T\sim MVN(0, Q_t), \quad t=0,..,T$$
Now let $E$ be the matrix of all the time series observations stacked on top of each other:
$$E = \begin{pmatrix} e_{0,0} & e_{1, 0}\\ e_{0,1} & e_{1, 1}\\... & ...\\e_{0,T} & e_{1,T}\end{pmatrix} $$
Lastly, let $\hat{E}$ be the vectorization of $E$:
$$\hat{E} = vec(E) = \begin{pmatrix} e_{0,0} & e_{0,1} & ... & e_{0,T}& e_{1,0} & e_{1,1} & ... & e_{1,T}\end{pmatrix}^T$$
My question is how to specify the covariance function of the vectorized $\hat{E}$
That is $$vec(E) \sim MVN(0, ???)$$
What I've tried:
It's clear that $E$ has a Matrix Normal distribution with some row and column covariance matrices $U$ and $V$, which implies that $\hat{E} \sim MVN(0, V\bigotimes U)$ but I don't see how to make this work since the column-wise covariance matrix is time dependent. If $Q_t = Q$ for all $t$ and the observations are IID than we would have $$\hat{E}\sim MVN(0, Q\bigotimes I)$$ but in this case $Q$ is time varying and there is potentially dependence across time
