Solution (source code)

= Solution

Define the <hat matrix> $H=X(X^TX)^{-1}X^T$. It is symmetric and idempotent, the <orthogonal projection> onto the column space of $X$. Since $HX=X$,
$$
\widehat Y=HY=X\beta+H\epsilon,\qquad
\boxed{\widehat Y\sim N_n(X\beta,\sigma^2H).}
$$
The <covariance> is $\sigma^2HH^T=\sigma^2H$. This <normal distribution> can be singular when $p<n$, since all fitted vectors lie in the $p$-dimensional model subspace.