# Does an explicit expression exist for the moments of the residuals in least squares regression?

Consider the linear regression model is $$Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i,$$ where $$X$$ is a random variable and the error has finite variance $$\sigma^2$$. When we solve this with least squares we find $$\hat \beta_0$$ and $$\hat \beta_1$$. These variables can be used to define $$\hat Y$$ which gives us an estimator of $$Y$$: $$\hat Y_i = \hat\beta_0 + \hat \beta_1 X_i.$$ The residuals are then defined as $$\hat e_i = \hat Y_i - Y_i.$$

I want to compute the expectation of a Taylor expansion of a function $$f$$ of the residuals. Since the residuals have an expected value of zero I get \begin{align} E[f(\hat e_i)] &= E[f(0) + f'(0)\hat e_i + \frac{f''(0)}{2}\hat e_i^2 + \frac{f'''(0)}{3!}\hat e_i^2 + \dots] \\ &= f(0) + \frac{f''(0)}{2}E[\hat e_i^2] + \frac{f'''(0)}{3!}E[\hat e_i^3] + \dots \end{align} To ensure this expansion is valid I need to know the moments of residuals. Are there explicit expressions for these moments? I have looked at many books but could not find anything.

$$y_i = \boldsymbol{x}_i^T\beta + \varepsilon$$ where $$\varepsilon_1, ..., \varepsilon_n \overset{IID}{\sim}\mathcal{N}(0, \sigma^2)$$ and $$\boldsymbol{x}_i^T = (1, x_{i1}, ...x_{ip})$$.
This model can be written in matrix form as: $$Y = X\beta + \boldsymbol{\varepsilon}$$ where $$Y\in\mathbb{R}^n$$ is the vector of the responses, $$X\in\mathbb{R}^{n \times p}$$ is the design matrix and $$\boldsymbol{\varepsilon} \sim\mathcal{N}(0, \sigma^2 I_n)$$ is a multivariate normal vector. The least square estimator is given by $$\hat\beta = (X^T X)^{-1}X^TY$$ and the residual $$\hat{e}_i$$, as you defined it, is given by $$\begin{array}{ccl} \hat{e_i} & = & \boldsymbol{x}_i^T\hat\beta - y_i\\ & = & \boldsymbol{x}_i^T(X^T X)^{-1}X^TY - y_i\\ & = & \boldsymbol{x}_i^T(X^T X)^{-1}X^T(X\beta + \boldsymbol{\varepsilon}) - y_i\\ & = & \boldsymbol{x}_i^T(X^T X)^{-1}X^TX\beta + \boldsymbol{x}_i^T(X^T X)^{-1}X^T\boldsymbol{\varepsilon} - y_i\\ & = & \boldsymbol{x}_i^T\beta - y_i +\boldsymbol{x}_i^T(X^T X)^{-1}X^T\boldsymbol{\varepsilon}\\ & = & -\varepsilon_i + \boldsymbol{x}_i^T(X^T X)^{-1}X^T\boldsymbol{\varepsilon}\\ & = & (-b_i^T + \boldsymbol{x}_i^T(X^TX)^{-1}X^T)\boldsymbol\varepsilon \end{array}$$ where $$b_i$$ is the vector of $$\mathbb{R}^n$$ made of zeros and a 1 at the $$i-th$$ position.
Now, as you know that $$\varepsilon \sim\mathcal{N}(0, \sigma^2 I_n)$$, using the property that for any full rank matrix $$M$$, if $$Z \sim\mathcal{N}(\boldsymbol{\mu}, \Sigma)$$, then $$MZ\sim\mathcal{N}(M\boldsymbol{\mu}, M\Sigma M^T)$$, you get that $$\hat{e}_i \sim{N}(0, s^2)$$ where $$\begin{array}{ccl} s^2 & = & \sigma^2(-b_i^T + \boldsymbol{x}_i^T(X^TX)^{-1}X^T)(-b_i^T + \boldsymbol{x}_i^T(X^TX)^{-1}X^T)^T\\ & = & \sigma^2 (1 - h_{ii}) \end{array}$$ where $$h_{ii} = \boldsymbol{x}_i^T(X^TX)\boldsymbol{x}_i$$ is the leverage of $$\boldsymbol{x}_i$$, between 0 and 1. From that, you can get the moments of the residuals using the moments of the normal distribution.
Getting the joint distribution of the vector of residuals $$\hat{\boldsymbol{e}}$$ is also possible since $$\hat{\boldsymbol{e}} = (I - H)\boldsymbol{\varepsilon}$$ where $$H = X(X^TX)^{-1}X^T$$ is the hat matrix, but it's a bit tricky since the matrix $$(I - H)$$ is singular.