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.


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

  • $\begingroup$ Well, its a singular multinormal distribution, so not that tricky ... $\endgroup$ Oct 13, 2021 at 12:08
  • $\begingroup$ Yes, you're right, it's not really tricky.. I edited the answer. $\endgroup$
    – Pohoua
    Oct 14, 2021 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.