Why a regression of OLS residuals on regressors, yields a $R^2$ of 0? In a regression of residuals
on X, the OLS estimator would be $\hat{\gamma} = (X'X)^{−1}X'\hat{\varepsilon}$ which is a vector of zeros ($\hat{\gamma} = 0$)
since $X'\hat{\varepsilon} = 0$. Could anyone tell me how did we arrive to the above conclusion about the OLS estimator?
 A: $$
\hat\beta_{ols} = \left(X^TX\right)^{-1}X^Ty\\\Downarrow\\
\left(X^TX\right)\hat\beta_{ols} = X^Ty\\\Downarrow\\
\left(X^TX\right)\hat\beta_{ols} - X^Ty = 0\\\Downarrow\\
X^T\left(X\hat\beta_{ols} - y\right) = 0\\\Downarrow\\
X^T\hat\epsilon = 0
$$
(When you go through the derivation of the OLS estimator, $\hat\beta_{ols}$, it is the second line that implies the first. Let's set that aside.)
Since $X^T\hat\epsilon = 0$, then $\left(X^TX\right)^{-1}X^T\hat\epsilon = \left(X^TX\right)^{-1}0 = 0$.
To get into why $R^2=0$ for this regression on the residuals, note that this means all of the coefficients (including the intercept) are zero. Also note that the mean of the residuals is zero. As $R^2$ compares the performance of your model (which always predicts zero) to the performance of a model that always predicts the mean of the outcome (which will be the residuals in the case, so a mean of zero), the two models have identical performance.
Take the $y_i$ to be observed residuals from the original model, $\hat y_i$ to be the predictions made by regressing on the residuals from the original model, and $\bar y$ to be the mean of the residuals (which will be zero).
$$
R^2 = 1-\left(
\dfrac{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -\hat y_i
\right)^2
}{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -\bar y
\right)^2
}
\right)
= 1-\left(
\dfrac{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -0
\right)^2
}{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -0
\right)^2
}
\right) = 1-1 = 0
$$
This is reflected in an R simulation.
set.seed(2023)
N <- 10                 # sample size
X <- rnorm(N)           # simulate a feature
Y <- X + rnorm(N)       # simulate a response
L <- lm(Y ~ X)          # fit a regression
y <- resid(L)           # extract the residuals of that regression
new_L <- lm(y ~ X)      # regress those residuals on the original feature
y_hat <- predict(new_L) # make predictions from that second regression "new_L"

1 - (sum((y - y_hat)^2))/(sum((y - mean(y))^2)) # This is R^2. I get 0.
```

A: The OLS estimates are solutions to
$$ \hat{\beta} = \min_{\beta\in\mathbb{R^n}} (y-X\beta)^T(y-X\beta) $$
You an show that the estimates also satisfy
$$ 0 = 2X^T(y-X\beta) $$
by differentiating the first equation with respect to $\beta$ and setting the resulting gradient to 0.
Now, note that $\varepsilon = (y-X\beta)$.  The consequence is that the residual is orthogonal to the columns of $X$.
Note also that $\beta = (X^TX)^{-1}X'y$.  If we regress the data onto the residual (yielding estimates $\beta_{\varepsilon}$) , then the orthogonality we're discovered means $\beta_{\varepsilon}=0$, which means the predicted residual $X\beta_{\varepsilon}=0$.
There is then no variation in $X\beta_{\varepsilon}$ and hence no variation in the residuals explained by $X$.
A: $\renewcommand{\epsilon}{\varepsilon}$
To clarify, let me first restate your question as follows:

Given a response vector $y \in \mathbb{R}^n$ and a design matrix $X \in \mathbb{R}^{n \times p}$.  After fitting $(y, X)$ to a linear model, the residual vector is $\hat{\epsilon} = (I - H)y$, where $H = X(X'X)^{-1}X'$ is the hat matrix.  Show that the $R^2 = 1 - \frac{SSE}{SST}$ of the linear model fitted by $(\hat{\epsilon}, X)$ is $0$.

To prove it, just note that when the response vector is $\hat{\epsilon}$ and design matrix is $X$, the average of new responses $\hat{\epsilon}_1, \ldots, \hat{\epsilon}_n$, denoted by $\bar{\hat{\epsilon}}$, equals to  $n^{-1}1'\hat{\epsilon} = 0$ (as the sum of residuals in a regression is $0$ -- under the condition that the regression model contains an intercept), and the new residual vector is $(I - H)\hat{\epsilon} = (I - H)(I - H)y = (I - H)y$, which coincides the old residual vector, thanks to the idempotency of $I - H$.  It then follows by definitions of $SST$ and $SSE$ that
\begin{align}
& SST = (\hat{\epsilon} - \bar{\hat{\epsilon}})'(\hat{\epsilon} - \bar{\hat{\epsilon}}) = \hat{\epsilon}'\hat{\epsilon} = y'(I - H)y, \\
& SSE = ((I - H)\hat{\epsilon})'(I - H)\hat{\epsilon} = 
((I - H)y)'(I - H)y = y'(I - H)y.
\end{align}
Therefore, $SST = SSE$, whence $R^2 = 0$.
Intuitively, since $\hat{\epsilon}$ is orthogonal to the space spanned by columns of $X$, $X$ does not provide any information to predict $\hat{\epsilon}$.  Therefore, $R^2 = 0$.
