Why are residuals Normally distributed?

When they assume the $$\varepsilon_i$$ have a normal distribution, with mean $$0$$ and variance $$\sigma^2,$$ then the estimated errors $$\hat\varepsilon$$ have a normal distribution with mean $$0$$ and variance $$\operatorname{Var}(\hat\varepsilon).$$

But what is the connection? Why, when errors are normally distributed, the estimated residuals also are normally distributed?

• Because this is a statistics site and we have been around a fairly long time, you can find a great deal written about this with a site search. The key phrase "hat matrix" is especially helpful.
– whuber
Jul 13 at 13:36

The connection is made through a specific procedure: namely, ordinary least squares. When you suppose a response $$Y$$ is a random variable related to other variables $$X$$ in the form

$$Y = X\beta + \varepsilon$$

where $$\varepsilon$$ is assumed to have a joint Normal distribution (plus some other assumptions that don't matter here), then the least squares solution is

$$\hat\beta = (X^\prime X)^{-}X^\prime Y = (X^\prime X)^{-}X^\prime (X\beta + \epsilon) = \beta + J\varepsilon$$

for a matrix $$J = (X^\prime X)^{-}X^\prime$$ that is computed only from $$X.$$ The residuals, or estimated errors, are the differences

$$\hat\varepsilon = Y - \hat Y = (X\beta + \varepsilon) - X\hat\beta = \varepsilon - XJ\varepsilon = (\mathbb I - H)\varepsilon$$

where $$H = XJ = X(X^\prime X)^{-}X^\prime$$ is the "hat matrix." This exhibits the residuals as linear combinations of $$\varepsilon$$ (with coefficients given by $$\mathbb I - H$$). Linear combinations of jointly Normal variables are Normal, QED.