Solution (source code)

= Solution

The <fitted values> are $\widehat Y=X\widehat\beta$, so the <hat matrix> is
$$
\boxed{H=X(X^TX)^{-1}X^T.}
$$
It is symmetric, and
$$
H^2=X(X^TX)^{-1}(X^TX)(X^TX)^{-1}X^T=H.
$$
It fixes every vector in the column space of $X$, and annihilates its orthogonal complement. Thus $H$ is the <orthogonal projection> onto that column space, has <matrix rank> and <trace> $p$, and $I-H$ is the complementary <orthogonal projection>. Consequently
$$
\boxed{\operatorname{Cov}(\widehat Y)=\sigma^2HH^T=\sigma^2H,\qquad
\operatorname{Cov}(\widehat e)=\sigma^2(I-H).}
$$
Since $(I-H)X=0$, the <covariance matrix> between the <regression residual> vector and the coefficient estimate is
$$
\operatorname{Cov}(\widehat e,\widehat\beta)
=\sigma^2(I-H)X(X^TX)^{-1}=0.
$$
These quantities have a joint <normal distribution>, so they are also independent, a stronger conclusion than zero <covariance>.

Now $\widehat e=(I-H)\varepsilon$ has mean zero. Using the expectation of a <quadratic form>,
$$
\mathbb E\,\mathrm{RSS}
=\mathbb E(\widehat e^T\widehat e)
=\operatorname{tr}\operatorname{Cov}(\widehat e)
=\sigma^2(n-p).
$$
Thus the <maximum-likelihood estimator> has
$$
\mathbb E\widehat\sigma^2_{\mathrm{ML}}=\frac{n-p}{n}\sigma^2,
\qquad \operatorname{Bias}(\widehat\sigma^2_{\mathrm{ML}})=-\frac pn\sigma^2.
$$
The unbiased alternative is
$$
\boxed{s^2=\frac{\mathrm{RSS}}{n-p}.}
$$
The <orthogonal projection of a Gaussian vector> also yields $\mathrm{RSS}/\sigma^2\sim\chi^2_{n-p}$, explaining the <residual degrees of freedom>.