Finding covariance matrix for MLE of correlated outputs

I generate data using the following model:

$\begin{pmatrix}Y_1\\Y_2\end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix}\mathbf{X}\beta_1\\\mathbf{X}\beta_2\end{pmatrix}, \mathbf{\Sigma} \right)$

I create $N$ test pairs, $Y_1^*, Y_2^*$, I'd like to fit a least-squares model to this data and evaluate its likelihood. So, I set:

$\hat\beta_1 = (\mathbf{X^TX})^{-1}\mathbf{X^T}Y_1^*$ and $\hat\beta_2 = (\mathbf{X^TX})^{-1}\mathbf{X^T}Y_2^*$

I'm having trouble computing $\mathbf{\hat\Sigma}$.

So far, I've tried

$\hat e_1 = Y_1^* - \mathbf{X}\hat\beta_1; \;\hat e_2 = Y_2^* - \mathbf{X}\hat\beta_2 \; \; (\star)$

$\mathbf{\hat\Sigma} = \frac{1}{N}\left( \begin{array}{ c c } e_1^Te_1 & e_1^Te_2\\ e_2^Te_1 & e_2^Te_2 \end{array} \right)$

But I've found that $\mathbf{\hat\Sigma}$ is pretty far from the true $\mathbf{\Sigma}$ that generated the data. Do I need to account somehow for the correlation between the two outputs when doing the estimation at step $(\star)$?

\begin{align*} Cov(\hat\beta_1, \hat\beta_2) &= Cov\left((\mathbf{X^TX})^{-1}\mathbf{X^T}Y_1^*, (\mathbf{X^TX})^{-1}\mathbf{X^T}Y_2^*\right) \\ & = (\mathbf{X^TX})^{-1}\mathbf{X^T} Cov(Y_1^*, Y_2^*) \mathbf{X} (\mathbf{X^TX})^{-1} \\ &= (\mathbf{X^TX})^{-1}\mathbf{X^T} \Sigma_{12} \mathbf{X} (\mathbf{X^TX})^{-1} \end{align*}
And a similar calculation can be used to find $Cov(\hat\beta_1, \hat\beta_1)= Var(\hat\beta_1)$, etc. (One may use the Kronecker product to write this all more compactly.)
Since in practice you don't know $\mathbf \Sigma$, one may plug in the consistent estimator $\boldsymbol{\hat{\Sigma}}$ (cross product of residuals) noted at the end of your question.
• Thanks for the solution and the reading material. I'm simulating data and comparing the performance of least-squares and expectation-maximization, so in this case I do know $\mathbf{\Sigma}$, but that's obviously not true in general. – Michael K Mar 4 '15 at 2:00