Vec operator and covariance matrix You have a matrix containing $T$ observations of each of $K$ random variables
\begin{align}
   U =
   \begin{bmatrix}
       u_{11} & \dots  & u_{1T} \\
       \vdots & \ddots & \vdots \\
       u_{K1} & \dots  & u_{KT}
   \end{bmatrix}.
\end{align} The covariance matrix of the $u_{kt}$'s is given by $\Sigma_u := E(u_t u_t')$ for $t = 1, \dots, T$ and where $ u_t = (u_{1t}, \dots, u_{Kt})'$ is a column vector by convention.
Now, we introduce the $vec$ operator which stacks columns of a matrix one over the next, i.e.
\begin{align}
    vec(U) = 
    \begin{bmatrix}
       u_{11} \\
       \vdots \\
       u_{1T} \\
       \vdots \\
       u_{K1} \\
       \vdots \\
       u_{KT}
    \end{bmatrix} =
    \begin{bmatrix}
       u_{1}' \\
       \vdots \\
       u_{K}'
    \end{bmatrix}
\end{align}
again using as convention that $u_k = ( u_{k1}, \dots, u_{kT} )'$ are column vectors for $k = 1, \dots, K$. It turns out that the covariance matrix of $vec(U)$ is
\begin{align}
    E\left( vec(U) vec(U)' \right) &= I_T \otimes \Sigma_u
\end{align}
which is to say, a block diagonal $TK \times TK$ matrix with entries $E( u_t u_t' )$. While I can see the intuition, I don't know how to formally establish this result. Any help would be appreciated.
 A: I think the easiest way is just index chasing, not anything elegant.
I'll assume the $u$ are all mean zero; if not, the identity does not hold, and the left-hand side isn't the covariance matrix.
A. $E[u_{it}u_{js}]=\sigma_{ij}$ if $t=s$ and zero otherwise.
B. Now write $S$ for the matrix on the right-hand side, and $s_{itjs}$ for the entry corresponding to $u_{it}$ and $u_{js}$. Clearly $E[u_{it}u_{js}]=s_{itjs}$; that's just what A said.
C. The entries of the right-hand side matrix (call it $S$) are also either $\sigma_{ij}$ or 0 for some $(i,j)$. In fact, $T$ of them are $\sigma_{ij}$ for each $(i,j)$, and $T(T-1)K^2$ of them are zero
D. So we just need to work out if the indices match up.   The $i$ index varies fastest, so the top left block of the left-hand side has $u_{i1,j1}$, then the next block to the right has $u_{i1,j2}$ and the first block on the second row has $u_{i2,j1}$ and so on. That is, the $(it,jt)$ entries -- the diagonal blocks -- have entries $\sigma_{ij}$ and the $(it,js)$ entries for $s\neq t$ are zero.  That's exactly $I_T\otimes \Sigma_u$.
More formally, if we write $(p,q)$ for the indices of the big matrices, then $p=i+(t-1)K$ and $p=j+(s-1)K$, and that works for both sides.
