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:


*

*Are there generalized projections (for higher-order tensors) such that the nice theory we have for second-moments/variances flow over?

*What is the distribution of $\text{Trace}(\mathcal{S}(\hat{\epsilon}))$? How does it depend on $X$?

 A: 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.
