Solution (source code)

= Solution

Write the one-column design as $x\in\mathbb R^n\setminus\{0\}$ and $S=x^Tx$. Maximizing the <normal linear model> likelihood first gives $\widehat\beta=x^TY/S$. Its minimized <residual sum of squares> is $R=\|Y-x\widehat\beta\|^2$. For variance $v>0$, the remaining log likelihood is $-\tfrac n2\log v-R/(2v)$ plus a constant, whose maximum is at
$$
\boxed{\widehat\sigma^2=\frac Rn.}
$$
Here $R>0$ almost surely because $n>1$ and the true variance is positive. In the exceptional data case $R=0$, the likelihood is unbounded as $v\downarrow0$, so there is no positive-variance maximum.

Let $P=xx^T/S$, the <orthogonal projection matrix> onto the design column. Then $\widehat\beta-\beta=x^T\varepsilon/S$ and $Y-x\widehat\beta=(I-P)\varepsilon$. These are orthogonal linear projections of an isotropic <multivariate normal distribution>, hence independent. An orthonormal coordinate system adapted to $x$ makes the residual squared length a sum of $n-1$ independent squared standard normals. Thus the <joint distribution of least-squares and variance estimators> is the product law
$$
\boxed{\widehat\beta\sim N\!\left(\beta,\frac{\sigma^2}{S}\right),\qquad \frac{n\widehat\sigma^2}{\sigma^2}\sim\chi^2_{n-1},\qquad\widehat\beta\ \text{and}\ \widehat\sigma^2\ \text{are independent}.}
$$
In particular the maximum likelihood variance estimator is biased: $\mathbb E\widehat\sigma^2=(n-1)\sigma^2/n$; the unbiased residual estimator uses denominator $n-1$.