Frisch-Waugh-Lovell Theorem: Partialing out a set of regressors I am trying to understand the result of the Frisch-Waugh-Lovell Theorem that we can partial out a set out regressors.
The model I am looking at is $y=X_1\beta_1 + X_2\beta_2 +u$
So the first step would be to regress $X_2$ on $X_1$:
$$
\begin{align}
X_2&=X_1\hat{\gamma}_1+\hat{w}\\
&=X_1\hat{\gamma}_1+M_{X_1}X_2
\end{align}
$$ with $M_X$ being the orthogonal projection matrix ($M_X=I-P_X$). The second step is then to regress $y$ on $X_1$:
$$
\begin{align}
y&=X_1\hat{\gamma}_2+\hat{v}\\
&=X_1\hat{\gamma}_2+M_{X_1}y
\end{align}
$$
And the third step would be to regress the residuals of the previous regressions on each other:
$$
\begin{align}
\hat{w}&=\hat{v}\beta_2+u\\
M_{X_1}y&=M_{X_1}X_2\beta_2+u
\end{align}
$$
Can somebody tell me whether this would be a correct way of describing this result and in particular whether the residuals in the last equation $u$ are equal to the residuals in the model?
 A: To see the residuals are the same consider that by definition of the residuals
$$y = X_1\hat \beta_1 + X_2\hat\beta_2 + \hat u$$
which is the first regression and its residuals $\hat u$. Now multiply with $M_{X_1}$ to get
$$M_{X_1}y =  M_{X_1}X_2\hat\beta_2 + \hat u,$$
because $M_{X_1} X_1 = \mathbf 0$ and $M_{X_1}\hat u = \hat u - P_{X_1}\hat u = \hat u$ because $P_{X_1}\hat u = \mathbf 0$. But the second equation is the regression of the residuals $M_{X_1}y$ from regression $y$ on  $X_1$, regressed on the residuals $M_{X_1}X_2$ being residuals from regression of $X_2$ on $X_1$. The residuals in the second equation are however also $\hat u$ so the residuals from the two regressions must be the same.
Heres R code to illustrate:
N <- 1000
z <- rnorm(N) # just a way to make x1 and x2 correlated
x1 <- rnorm(N) + z
x2 <- rnorm(N) + z
cov(x1,x2)
y <- 2 + 2*x1 - 2*x2 + rnorm(N) 

e_y_x1 <- lm(y~x1)$residuals
e_x2_x1 <- lm(x2~x1)$residuals

u1 <- lm(e_y_x1~e_x2_x1)$residuals
u2 <- lm(y ~ x1 + x2)$residuals

plot(u1,u2)

