Solution (source code)

= Solution

Assume $n>p\ge1$, the ordinary nonsaturated <Gaussian linear model>. The log <likelihood> is, up to a constant,
$$
\ell(\beta,\sigma^2)=-\frac n2\log\sigma^2-\frac{\|Y-X\beta\|^2}{2\sigma^2}.
$$
The normal equations and full column rank give the <least-squares estimator> $\widehat\beta=(X^TX)^{-1}X^TY$. Write $H=X(X^TX)^{-1}X^T$ and $Q=I-H$, an <orthogonal projection> of rank $n-p$. Maximizing in $\sigma^2$ yields
$$
\boxed{\widehat\sigma^2=\frac{Y^TQY}{n}.}
$$
Since $QX=0$, the residual sum of squares is $\varepsilon^TQ\varepsilon$. Expanding its <expectation> uses only $\mathbb E\varepsilon_i\varepsilon_j=\sigma^2\delta_{ij}$:
$$
\mathbb E(Y^TQY)=\sum_{i,j}Q_{ij}\mathbb E(\varepsilon_i\varepsilon_j)
=\sigma^2\operatorname{tr}Q=(n-p)\sigma^2.
$$
This does not use <Cochran's theorem>. It proves bias $-p\sigma^2/n$ and gives
$$
\boxed{\mathbb E\widehat\sigma^2=(1-p/n)\sigma^2,\qquad
\widetilde\sigma^2=\frac{Y^TQY}{n-p},\qquad\mathbb E\widetilde\sigma^2=\sigma^2.}
$$
In the saturated case $n=p$, residuals vanish identically, the likelihood is unbounded as the positive <variance> tends to zero, and the unbiased residual <estimator> is unavailable. Thus the usual variance-estimation request requires $n>p$.