Effect of omitting interaction term on OLS estimators

Question

Say we know that $Y$ follows the model $$Y = \beta_0 + \beta_1X_1 + \beta_2X_2+\beta_3X_1X_2+\epsilon$$ Suppose that $X_1$ is a categorical binary predictor, while $X_2$ is a continuous predictor that has been standardized to have sample mean $0$ and variance $1$.

We fit a model for $Y$ on $X_1$ and $X_2$, without the interaction term $X_1X_2$.

How does one derive an expression for the bias of the least squares estimators $\hat{\beta_1}$ and $\hat{\beta_2}$ in this case? I'm having trouble determining what the $X'X$ matrix is supposed to look like. I would appreciate any help or hints in the right direction!

Answer 1

I would start by demeaning the model

$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1X_2 + \epsilon$$

to get

$$\mathbb E[Y] = \beta_0 + \beta_1 \mathbb E[X_1] + \beta_2 \mathbb E[X_2] + \beta_3 \mathbb E[X_1X_2] + \mathbb E[\epsilon]$$

substract from the model to get demeaned model

$$\tilde Y = \beta_1\tilde X_1 + \beta_2 \tilde X_2 + \beta_3 \widetilde{X_1X_2} + \tilde \epsilon$$

define $\tilde X = (\tilde X_1,\tilde X_2)^\top$ and get

$$\tilde X\tilde Y = \tilde X \tilde X^\top \beta + \beta_3 \tilde X\widetilde{X_1X_2} + \tilde X\tilde \epsilon$$

then take expectations

$$\mathbb E[\tilde X\tilde Y] = \mathbb E[\tilde X \tilde X^\top] \beta + \beta_3\mathbb E[ \tilde X\widetilde{X_1X_2}]$$

and multiply with $\mathbb E[\tilde X \tilde X^\top]^{-1}$ to get

$$plim (\hat \beta_{OLS}) = \mathbb E[\tilde X \tilde X^\top]^{-1}\mathbb E[\tilde X\tilde Y] = \beta + \beta_3\mathbb E[\tilde X \tilde X^\top]^{-1}\mathbb E[ \tilde X\widetilde{X_1X_2}]$$,

such that $bias = \beta_3\mathbb E[\tilde X \tilde X^\top]^{-1}\mathbb E[ \tilde X\widetilde{X_1X_2}]$. Also, note that $Var(X)^{-1}=\mathbb E[\tilde X \tilde X^\top]^{-1}$ and $Cov(X,Y) = \mathbb E[\tilde X \tilde Y]$.

Here is some R code to verify.

N <- 10000

x1 <- as.numeric(runif(N)<0.5)
x2 <- rnorm(N)
e <- rnorm(N)

b0 <- 1
b1 <- -1
b2 <- 1
b3 <- 2

y <- b0 + b1*x1 + b2*x2 + b3*x1*x2 + e
lm(y~x1+x2)

# OLS estimator
ols <- solve(var(cbind(x1,x2)))%*%cov(cbind(x1,x2),y)
ols
# Bias
bias <- b3 * solve(var(cbind(x1,x2)))%*%cov(cbind(x1,x2),x1*x2)

ols - bias