Residualizing dependent variable and two step linear regression

Question

Assume we have a DGP of the form

$y = \beta_0 + \beta_1 * x_1 + \beta_2 * x_2 + \beta_3 * x_3 + \epsilon$

where $\epsilon$ is a standard i.i.d. error term. Does residualizing $y$ using a linear regression including only an intercept, $x_1$ and $x_2$ and regressing the resulting residuals on $x_3$ give us the same estimate for the coefficient of $x_3$ as in the full specification?

More formally, if you estimate the model in a two step procedure of the form

$y = \bar{\beta_0} + \bar{\beta_1} * x_1 + \bar{\beta_2} * x_2 + \eta$

$\hat{\eta} = \bar{\beta_3} * x_3 + \xi$

where $\eta$ and $\xi$ are error terms and $\hat{\eta}$ are the estimated residuals from the first regression, is it proven that $\bar{\beta_3}$ and $\beta_3$ will be the same?

If I simulate a model in R, it seems to work flawlessly. However, I have a problem when trying to replicate the two step procedure with a real dataset. Now I'm not sure whether I made a coding error or whether this two step procedure is wrong.

Simulation Code

x1 <- rnorm(1000)
x2 <- rnorm(1000)
x3 <- rnorm(1000)

y <- 4 + 3*x1 - 2*x2 + 22*x3 + rnorm(1000)

res <- resid(lm(y~x1+x2))
summary(lm(res~x3-1))

Answer 1

It general these regressions are not the same, but since you have simulated independence it works out.

Consider the full regression \begin{align*} Y = \alpha + X_1\beta_1 + X_2\beta_2 + X_3\beta_3 +\epsilon \end{align*}

Let me define a $n\times 3$ matrix $Z = [1,X_1,X_2]$, a column of ones then the values of $X_1$ and $X_3$. Let $M_z$ be the residual making projection matrix. In other words this matrix when applied to a variable is the same as regression that variable on a column of ones, $X_1$ and $X_2$.

\begin{equation} \hat{\beta_3} = \frac{X_3^TM_zY}{X_3^TM_zX_3} \end{equation}

The above formula follows from the application of the FWL theorem. You can derive the result equally from minimizing the sum of squared residuals, but the matrix notation and FWL theorem make things much cleaner and in this case give us insight into the question you asked.

Recall that in simple regression, the general formula for the OLS estimate is $\hat{\beta} = \frac{X^TY}{X^TX}$. You can show for yourself using this formula that the $\beta_3$ that we derived above is the same as the following simple regression: \begin{align*} M_zY = M_zX_3\beta_3 + \epsilon \end{align*}

(Note: $(M_z)^TM_z = M_z$ by properties of orthogonal projection matrices. Second note, I wrote $\epsilon$ again because these residuals are numerically the same by FWL. In other words this is the exact same regression).

So to answer the question, it would be same as residualizes Y by regression in a constant, $X_1$ and $X_2$ and then regressing these residuals on the residuals of $X_3$ regressed on a constant, $X_1$ and $X_2$.

In you simulation you have all your $Xs$ as independent standard normal variables. In expectation when you regress $X_3$ on the others, the result is the same so the residualizing of $X_3$ doesn't matter. In the real data, my guess is that $X_3$ is not mean 0 and has some real relationship with the other two variables so it will not work.