Solution (source code)

= Solution

The <least-squares estimator> minimizes $\|Y-Xb\|^2$. Differentiating in $b$ gives $-2X^T(Y-Xb)=0$, hence $\boxed{X^TX\widehat\beta=X^TY}$. Full column rank makes $X^TX$ positive definite, so
$$
\widehat\beta=(X^TX)^{-1}X^TY=\beta+(X^TX)^{-1}X^T\epsilon.
$$
This is a <linear transformation> of a <multivariate normal distribution>. Its mean is $\beta$ and its <covariance> is $\sigma^2(X^TX)^{-1}$. Therefore $\boxed{\widehat\beta\sim N_p(\beta,\sigma^2(X^TX)^{-1})}$.