# Higher moments of linear regression residuals?

I previously asked this on Math StackExchange, with no success, but this post will add to that with some simulations.

Background

In the following linear regression with i.i.d $$\epsilon_i$$ $$(i = 1, \cdots, n)$$ with mean 0 finite moments of all order (say, $$\mu_2 = \sigma^2$$ for the centralized second moment, $$\mu_3$$ for centralized third moment, etc), \begin{align*} Y_i = X_i^\intercal\beta + \epsilon_i \end{align*} we know the least-squares estimator for $$\beta$$ is \begin{align*} \hat{\beta}=(X^\intercal X)^{-1}X^\intercal Y \end{align*} where $$X = (X_1, \cdots, X_n)^\intercal \in \mathbb{R}^{n\times p}$$, $$Y = (Y_1,\cdots, Y_n)^\intercal \in \mathbb{R}^n$$. Now, we can easily derive the variance of $$\hat{\epsilon} = Y - X\hat{\beta}$$ with \begin{align*} \text{Var}(\hat{\epsilon}) &= \text{Var}((I-H)Y) \\ &\overset{\heartsuit}{=} (I-H)\text{Var}(Y)(I-H)^\intercal \\ &= \sigma^2(I-H)(I-H)^\intercal \\ &\overset{\spadesuit}{=} \sigma^2(I-H) \end{align*} where $$H = X(X^\intercal X)^{-1}X^\intercal$$ is a projection matrix of rank $$p$$, thus justifying equality $$(\spadesuit)$$. Hence, this allows us to find an unbiased estimator for $$\sigma^2$$, since \begin{align*} \text{Trace}(\text{Var}(\hat{\epsilon})) = \sigma^2\text{Trace}(I-H) = \sigma^2(n-p) \implies \widehat{\sigma^2}=(n-p)^{-1}\|Y - \hat{Y}\|_2^2 \end{align*} An attempt for higher moments

How would we generalize this computation in coming up with unbiased estimators for the third centralized moments \begin{align*} \mu_3 \overset{\text{def}}{=} \mathbb{E}\epsilon^3_i \end{align*} or beyond? Here's my attempt. Define \begin{align*} \mathcal{S}(\hat{\epsilon}) = \mathbb{E}(\hat{\epsilon}_i - \mathbb{E}\hat{\epsilon}_i)^{\otimes 3} \end{align*} Expanding the tensor product, we get \begin{align*} \mathcal{S}(\hat{\epsilon}) &= [(\mathcal{S}(Y)\times_1 M)\times_2 M]\times_3 M \\ &=\mu_3[(I_{n\times n\times n}\times_1 M)\times_2 M]\times_3 M \end{align*} where $$M = I - H$$, and $$\times_i$$ is defined by the $$n$$-mode product. This just means that \begin{align*} [\mathcal{S}(\hat{\epsilon})]_{ijk} = \mu_3 \sum_{v=1}^{n}M_{iv}M_{jv}M_{kv} \end{align*} So the trace is \begin{align*} \text{Trace}(\mathcal{S}(\hat{\epsilon})) = \sum_{i=1}^{n}[\mathcal{S}(\hat{\epsilon})]_{iii} = \sum_{i=1}^{n}\sum_{v=1}^{n}M_{iv}^3 \end{align*} This all looks reasonable, considering that \begin{align*} \text{Trace}(\text{Var}(\hat{\epsilon})) = \sum_{i=1}^{n}\sum_{v=1}^{n}M_{iv}^2 \end{align*} I was curious what the third moment trace might equal, so I conducted some simulations:

p = 5
n = 100
X = matrix(rnorm(p*n), nrow = n, ncol = p)
H = X%*%solve(t(X)%*%X)%*%t(X)
M = diag(n) - H

#Trace of variance (second-moment outer product)
sum(M^2)             #Always = n - p
#Trace of third-moment outer product
sum(M^3)             #Not only doesn't equal n - p, but changes for each X!


This is somewhat disappointing, since we no longer have a degree of freedom interpretation for higher moments.

Future

Some questions remain:

1. Are there generalized projections (for higher-order tensors) such that the nice theory we have for second-moments/variances flow over?
2. What is the distribution of $$\text{Trace}(\mathcal{S}(\hat{\epsilon}))$$? How does it depend on $$X$$?
• It seems to me that under the (implicitly) general assumptions you have made--namely, only that $\epsilon_i$ are iid of zero mean and variance $\sigma^2$--it is impossible to obtain an unbiased estimator of any higher moment. After all, those moments might not exist or they might be infinite. For you to have any chance at a "nice theory," I think you will need to specify a parametric family of distributions for the $\epsilon_i.$
– whuber
Feb 1 '19 at 20:02
• Sure, I've edited. Assume finite moments of all orders. I shouldn't try and specify a fully parametric distribution for $\epsilon_i$, because linear model theory all follows without a parametric assumption on $\epsilon_i$ for the second moments, and I'd like to see if something similar can follow for third. If this is too restrictive, I might consider assuming $\epsilon_i$ are normal. Feb 1 '19 at 20:10
• Permit me to repeat then: for a nonparametric model it's unlikely you can find any unbiased estimator. If you assume normality the theory becomes almost trivial because the first two moments determine all the others.
– whuber
Feb 1 '19 at 20:18
• Even assuming normality doesn't make the theory trivial. The difficulty does not lie with $\mu_3$ or whatever higher moments, but rather the $3$-mode products from $M = I - H$, which is a function of $X$. Even if $X$, along with $\epsilon_i$, were assumed normal, the result isn't obvious. Feb 1 '19 at 20:29
• The best unbiased estimator of the odd central moments in the Normal case is $0.$ That looks pretty trivial to me. I'll agree, before you post any rejoinder, that this does not imply it is trivial to find unbiased estimators for the non-central moments: but it clearly points the way.
– whuber
Feb 1 '19 at 20:32

Your method here involves an attempt to derive the higher-order moments of the residual vector $$\hat{\boldsymbol{\varepsilon}}$$ from its relationship to the error vector $$\hat{\boldsymbol{\varepsilon}}$$. If you want to do that then you will need to prescribe a distributional family for the error vector. So, for example, the distribution of the trace is determined by the underlying distribution of the error vector --- you cannot derive its distribution without this. If you assume that the errors are normally distributed then it is possible to derive the higher-order moments of the residual vector, but the results are going to be quite complicated (and are probably best computed by simulation).
On the other hand, if you wish to proceed without an assumption about the error distribution, then you would not attempt to determine the moments of interest from the underlying error vector at all. Instead, you would simply use standard methods to estimate the higher-order moments empirically from the residuals. For large $$n$$ (and under mild assumptions about the regressors) the residuals become asymptotically uncorrelated, which should allow you to use standard sample estimators and have good convergence properties.