Assume we have the following system of equations
Eq. 1: $z = x\gamma + \epsilon$
Eq. 2: $y = z\theta + x\beta + \eta$
where the error terms fulfill standard assumptions. Simulating data from Eq. 1 and Eq. 2 and estimating Eq. 2 using OLS gives an unbiased estimate of $\theta$ because we condition on $x$. Similarly, we can use a different parametrization of the problem and write
Eq. 3: $y = (z - x\gamma)\theta + x\rho + \eta$
where we made explicit that $\theta$ is the effect of the exogenous variation in $z$ on $y$. Using Eq. 3 to simulate the data and Eq. 2 to estimate the parameters results in the same estimate of $\theta$ as before. So far, so good.
I am interested in including an interaction term between $x$ and $z$ in the model like follows
Eq. 4: $y = z\theta + zx\delta + x\tau + \kappa.$
Simulating data from Eq. 4 and estimating the parameters using OLS gives unbiased estimates of $\delta$ and $\theta$. Based on these results, I was expecting that simulating from
Eq. 5: $y = (z - x\gamma)\theta + (z - x\gamma)x\delta + x\nu + \zeta$
and estimating a specification like in Eq. 4 using OLS will also give unbiased results for $\theta$ and $\delta$. However, while the estimate of $\theta$ is unbiased, the estimate of $\delta$ is now heavily biased. I do not really have an explanation for this, but I assume that it has to do with the fact that $z$ is exogenous conditional on $x$, while $z*x$ is not exogenous conditional on $x$. Any thoughts are appreciated. Here is my simulation code in R.