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.
 A: Let's take a classical linear regression model:
$$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: $\hat{\boldsymbol{e}}$ follows a singular multivariate normal distribution (singular since its variance matrix $(I - H)$ is singular).
