Solution (source code)

= Solution

Write $A=(X^TX)^{-1}X^T$. Since $AX=I$, the estimator is $\widehat\theta=\theta+A\epsilon$. Therefore
$$
\boxed{\mathbb E\widehat\theta=\theta,\qquad
\operatorname{Var}(\widehat\theta)=A(\sigma^2I)A^T=\sigma^2(X^TX)^{-1}.}
$$
Its <expectation> proves unbiasedness. It is also jointly normal, being an affine transformation of the normal error <vector>.