Solution (source code)

= Solution

Assume the known <covariance matrix> $\Sigma$ is <positive-definite matrix>[positive definite]. The <Gaussian likelihood>, up to constants, is
$$
\ell(\beta)=-\tfrac12(Y-X\beta)^T\Sigma^{-1}(Y-X\beta).
$$
Its gradient is $X^T\Sigma^{-1}(Y-X\beta)$. The negative <Hessian matrix> $X^T\Sigma^{-1}X$ is <positive-definite matrix>[positive definite] because $X$ has <full column rank>. Thus the unique <maximum-likelihood estimator> is the <generalized least squares> estimator
$$
\boxed{\widehat\beta_{\mathrm{GLS}}=(X^T\Sigma^{-1}X)^{-1}X^T\Sigma^{-1}Y.}
$$
It is unbiased, and multiplication of its linear coefficient matrix by $\Sigma$ gives
$$
\boxed{\operatorname{Cov}(\widehat\beta_{\mathrm{GLS}})=(X^T\Sigma^{-1}X)^{-1}.}
$$
In fact $\widehat\beta_{\mathrm{GLS}}\sim N_p(\beta,(X^T\Sigma^{-1}X)^{-1})$. Equivalently, a <whitening transformation> premultiplies $Y$ and $X$ by $\Sigma^{-1/2}$, after which <ordinary least squares> applies. There is no extra $\sigma^2$ factor when $\Sigma$ is the complete known covariance; such a factor appears if only a shape matrix $\Omega$ is known and $\Sigma=\sigma^2\Omega$.