# Effect of omitting interaction term on OLS estimators

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!

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