Solution (source code)

= Solution

Put $C=(X^TX)^{-1}$ and $\nu=n-p$. The normal law of the estimator gives $(\widehat\theta_i-\theta_i)/(\sigma\sqrt{C_{ii}})\sim N(0,1)$. The fitted and residual Gaussian projections are <orthogonal> and hence independent, so this quantity is independent of $\nu s^2/\sigma^2\sim\chi^2_\nu$. Their ratio therefore has <Student's t-distribution>:
$$
\frac{\widehat\theta_i-\theta_i}{s\sqrt{C_{ii}}}\sim t_\nu.
$$
If $t_{\nu,1-\alpha/2}$ is its upper $1-\alpha/2$ <quantile>, rearranging its central <probability> event gives
$$
\boxed{\left[\widehat\theta_i-t_{\nu,1-\alpha/2}s\sqrt{C_{ii}},\quad
\widehat\theta_i+t_{\nu,1-\alpha/2}s\sqrt{C_{ii}}\right].}
$$
This interval has exact coverage $1-\alpha$ under the given normal linear model, rather than merely large-sample coverage.