Solution (source code)

= Solution

For $\sigma^2>0$, the <Gaussian linear model> has <likelihood function>
$$
L(\beta,\sigma^2;y)=(2\pi\sigma^2)^{-n/2}\exp\!\left[-\frac{(y-X\beta)^T(y-X\beta)}{2\sigma^2}\right].
$$
Maximizing over $\beta$ minimizes the residual sum of squares. Full column rank makes $X^TX$ positive definite, and differentiating gives the unique <least-squares estimator>
$$
\boxed{\widehat\beta=(X^TX)^{-1}X^TY.}
$$
Put $P=X(X^TX)^{-1}X^T$ and $R=\|(I-P)Y\|^2$. If $R>0$, differentiating the profile log-likelihood in $\sigma^2$ gives
$$
\boxed{\widehat\sigma^2=R/n.}
$$
The denominator is $n$, not the residual degrees of freedom $n-p$; the latter gives the unbiased estimator.

A linear transformation of the <multivariate normal distribution> remains normal, so
$$
\boxed{\widehat\beta\sim N_p(\beta,\sigma^2(X^TX)^{-1}).}
$$
Its <covariance> with the residual vector is $\sigma^2(X^TX)^{-1}X^T(I-P)=0$. Jointly Gaussian vectors with zero cross-covariance are independent, hence $\widehat\beta$ is independent of $(I-P)Y$, and therefore of $\widehat\sigma^2$, a function of that residual. Orthogonal coordinates also give $R/\sigma^2\sim\chi^2_{n-p}$.

\b[The allowed boundary case $p=n$ needs qualification.] Then $P=I$, $R=0$ for every observation, and the likelihood increases without bound as $\sigma^2\downarrow0$ after fitting $Y$ exactly. There is no maximum in the positive-variance model; $R/n=0$ is only a formal boundary estimate. When $n>p$, $R>0$ almost surely under a genuinely positive <variance>, and the displayed MLEs apply. A zero residual in any exceptional exact-fit sample has the same nonexistence issue. The estimator $\widehat\beta$ and its <normal distribution> remain well-defined in the square full-rank case.