Solution (source code)

= Solution

The <regression residual> vector is $e=Y-X\widehat\beta=(I-H)Y$, where $H=X(X^TX)^{-1}X^T$ is the <hat matrix>. Both $H$ and $I-H$ are symmetric <idempotent> matrices, and $(I-H)X=0$. Therefore $e=(I-H)\varepsilon$ and the <linear image of a multivariate normal vector> gives
$$
\boxed{e\sim N_n\bigl(0,\sigma^2(I-H)\bigr).}
$$
This <multivariate normal distribution> is singular, with rank $n-p$: <fitted-residual orthogonality> forces $X^Te=0$. Specifically, $\operatorname{Var}(e_i)=\sigma^2(1-h_{ii})$ and $\operatorname{Cov}(e_i,e_j)=-\sigma^2h_{ij}$ for distinct observations. Thus <regression residuals> are generally neither independent nor identically distributed, even though the model errors are. Also $\operatorname{RSS}/\sigma^2\sim\chi^2_{n-p}$ and is independent of $\widehat\beta$, which justifies the <standard errors> and exact <Student t-tests> below.