# Not recovering true coefficient with recursive bivariate probit model on simulated data

I have built a simulated dataset to try to build my intuition about the recursive bivariate probit model. The challenge I'm running into is that I'm unable to recover the true coefficient in my simulated data in the presence of an unobserved confounder, despite having an instrumental variable.

The true model for observation $$i$$ is as follows:

\begin{align*} X^*_i &= \beta_0+\beta_1Z_i+\beta_2O_i+\beta_3U_i+\zeta_i \\ X_i &= \mathbb{I}[X^*_i > 0] \\ Y^*_i &= \beta_4+\beta_5X_i+\beta_6O_i+\beta_7U_i+\epsilon_i \\ Y_i &= \mathbb{I}[Y^*_i > 0] \end{align*}

I want to learn the relationship between $$X_i$$ and $$Y_i$$ (so the goal is to accurately estimate $$\beta_5$$); I have instrumental variable $$Z_i$$ (which affects $$X_i$$ but not $$Y_i$$), as well as observed factor $$O_i$$ and unobserved factor $$U_i$$, both of which impact both $$X_i$$ and $$Y_i$$.

I simulate $$Z_i$$, $$O_i$$, $$U_i$$, $$\zeta_i$$, and $$\epsilon_i$$ as IID standard normal random variables, and I use coefficients $$\beta_0=\beta_1=\beta_4=\beta_5=-1$$ and $$\beta_2=\beta_3=\beta_6=\beta_7=1$$. With the full model specification (including access to the unobserved confounder $$U_i$$), I can perfectly estimate the coefficients using the recursive bivariate probit model:

bivariateProbit(cbind(b0=1, b1=Z, b2=O, b3=U), cbind(b4=1, b5=X, b6=O, b7=U), X, Y)
#            b0            b1            b2            b3            b4
# -1.0004121410 -0.9883708493  0.9916502262  1.0030378739 -0.9986777626
#            b5            b6            b7           rho
# -1.0003586241  0.9967899947  1.0005921044  0.0007548157


However, once I drop $$U_i$$, my estimate for $$\hat\beta_5$$ becomes -0.717, quite far from the true value of -1:

bivariateProbit(cbind(b0=1, b1=Z, b2=O), cbind(b4=1, b5=X, b6=O), X, Y)
#         b0         b1         b2         b4         b5         b6        rho
# -0.7077195 -0.6983809  0.7046800 -0.7028858 -0.7172548  0.7091427  0.5080552


What is causing the estimate to vary so far from the true value once I drop my unobserved confounder? Are there better approaches to obtain an estimate closer to the true value of -1?

R code to construct the example (note that bivariateProbit is a function I implemented to make sure I understood how the bivariate probit model worked):

set.seed(144)
N <- 100000  # Observations
Z <- rnorm(N) ; O <- rnorm(N) ; U <- rnorm(N) ; zeta <- rnorm(N) ; eps <- rnorm(N)
b0 <- -1 ; b1 <- -1 ; b2 <- 1 ; b3 <- 1 ; b4 <- -1 ; b5 <- -1 ; b6 <- 1 ; b7 <- 1
X <- as.numeric(b0 + b1*Z + b2*O + b3*U + zeta > 0)
Y <- as.numeric(b4 + b5*X + b6*O + b7*U + eps > 0)
library(pbivnorm)
bivariateProbit <- function(X1, X2, y1, y2) {
optim(setNames(rep(0, ncol(X1)+ncol(X2)+1), c(colnames(X1), colnames(X2), "rho")),
function(beta) -sum(log(pbivnorm((2*y1-1)*as.numeric(X1 %*% head(beta, ncol(X1))), (2*y2-1)*as.numeric(X2 %*% beta[seq(ncol(X1)+1, ncol(X1)+ncol(X2))]), (2*y1-1)*(2*y2-1)*tail(beta, 1)))),
lower=rep(c(-Inf, -1), c(ncol(X1)+ncol(X2), 1)), upper=rep(c(Inf, 1), c(ncol(X1)+ncol(X2), 1)), method="L-BFGS-B")\$par
}