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

Question

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.

Answer 1

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).