Will MSE-based estimator generate symmetric residuals if the error has got symmetric support (not distribution)? This question is more specific than :my old question
Take follow regression model:
$y=f(x)+e$
Where $e\sim D$ with a such symmetric support $A=(-a,a)$, not symmetric distribution. Now given a data set $\{x_i,y_i\}_{i=1}^{N}$ and a estimator $g(.,\mathbf{\beta})$ to $f$.
Will I always get symmetric residuals if:
$\beta=argmin_{\beta}(\sum_{i=1}^{N}(g(x_i,\beta)-y_i)^2)$
I've simulated a regression problem where $y=x^2+e$ with $e$ having the follow distribution:
$f(e)=\frac{1}{exp(1)-exp(-1)}exp(e) I(e)_{\{-1,1\}}$
That is not symmetric at all:

But using random forest with MSE impurity to train $\{x_i,y_i\}_{i=1}^{N}$ I got these residuals:

Where are clearly symmetric. I think it'll always happen but I'm not able to proof it.

STILL NO
I have a similar argument as I gave in your previous question. Let's do the same model but with an error term that is an asymmetric beta distribution that is shifted to have support on $[-0.5, 0.5]$. I then take this a step further and multiply by a scaling factor to expand the support to $[-50, 50]$. I encourage readers to play with the a, b, and m values to see what happens.
set.seed(2021)
N <- 10000
a <- 1
b <- 3
m <- 100
x <- seq(0, 10, 10/(N - 1))
e <- rbeta(N, a, b) - 0.5 # subtract 0.5 so support is [-0.5, 0.5], not [0, 1]
e <- m * e # Multiply by m to stretch or compress the support to [-m/2, m/2]
y <- x + e
plot(x, y)
L <- lm(y ~ x)
par(mfrow = c(2, 1))
hist(resid(L))
hist(e)
par(mfrow = c(1, 1))

The empirical range of my e values is $-49.99677$ to $45.6901$, so the support looks pretty symmetric about $0$ (which we expect from the construction). However, neither the residuals nor the error term are symmetric, as the histograms below show.

