How do you calculate a confidence interval for a parameter in multiple linear regression? Suppose $Y = Xβ + ε$, where $Y$ is a random $n$-vector, $X$ is a random $k × n$ matrix, $β$ is a fixed $k$-vector, and $ε$ is a random $n$-vector where the coordinates are independent (of each other and of $X$) and distributed as $N(0, σ^2)$. Given least-squares estimates $\hat β_1,\; \hat β_2,\; …,\; \hat β_k$ of the parameters, how do I compute a $(1 - α)$ confidence interval for one of the parameters $β_i$? Please provide a reference or a proof.
Although this is a pretty fundamental question, I couldn't find a previous CV question that exactly answers it, and poking around textbooks and websites seems to show some disagreement (e.g., should you use the quantile of a normal distribution or a $t$ distribution when multiplying by a standard error?). If there are several reasonable approaches, I'd appreciate a discussion of the differences.
 A: The covariance matrix of the parameters is the residual variance multiplied by the inverse of the quadratic form of the design matrix. It is a result of the Gauss-Markov theorem (or, specifically, such a variance is a Best Linear Unbiased Estimator: BLUE).
Using matrix notation $\hat{\beta} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^TY$.
PROOF: The variance is:
$\text{var}(\hat{\beta}) = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T var(Y|X) \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}$
$= (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T diag(\sigma^2) \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}$
$= \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}$
$= \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}$.
One result which does not require normal residuals: The limiting distribution of the $\hat{\beta}$ is normal.
Thus a 95% CI can be constructed: $\hat{\beta} \pm \mathcal{Z}_{1-\alpha/2} \sqrt{diag(\sigma^2 (\mathbf{X}^T\mathbf{X})^{-1})}$
Another result when we know residuals are normal but their variance unknown: the estimator is a linear combination of normals and thus is normal. To obtain finite sample inference, the values follow a T-distribution with degrees of freedom equal to the sample size minus the number of free parameters (including the intercept) so that the CI is expressed as:
$\hat{\beta} \pm \mathcal{T}_{1-\alpha/2, n-p} \sqrt{diag(\sigma^2 (\mathbf{X}^T\mathbf{X})^{-1})}$
Notationally, beta-hat a vector, the SEs also a vector and the $\pm$ a $\mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R} \times \mathbb{R} $ operator, basically a sum with an outer product with the [-1, 1] vector.
Seber Lee, Linear Regression Analysis 2nd ed
