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, but it's a bit tricky since the matrix $(I - H)$ is singular.