Joint distribution of residuals in a simple linear regression model with iid normal errors

Question

In a simple linear regression $Y=X\beta + \varepsilon$, residuals are given by $\hat\varepsilon = M\varepsilon$, where $M = I_n - P$ is the annihilator matrix, and $P = X(X^TX)^{−1}X^T$ is the projection matrix, and $X$ is the design matrix. Assuing that the errors $\varepsilon$ are iid normal with mean $0$ and standard deviation $\sigma$, what is the joint (conditional on $X$) distribution of the residuals $\hat\varepsilon$? It seems to me that it should be multivariate normal, but what is the covariance matrix $\Sigma$?

Answer 1

Note that the residuals $\hat{\varepsilon}$ are linear transformation of $Y$, that are multivariate normal, hence, as $\hat{\varepsilon} = (I-H)Y$, thus $\hat{\varepsilon}$ is multivariate normal with $\mu = (I-H)EY = X\beta - X\beta = 0$ and $\Sigma_{\hat{\varepsilon}} = (I-H)\Sigma_Y(I-H) = \sigma_{Y}^2(I-H) $.