I have an instrumental variable (IV) problem where both the endogenous variable $y_2$ and the independent variable $y_1$ are binary. Thus, I need a generalised linear model in both stages. Using the probit link function $Φ^{-1}$, the equations for the first and second stage, respectively, might look as follows:
$$ Φ^{-1}(y_2) = z_1 + z_2 + z_3 + \varepsilon_2 $$$$ Φ^{-1}(y_2) = z_1 + z_2 + z_3 $$ $$ Φ^{-1}(y_1) = z_1 + y_2 + z_3 + \varepsilon_1 $$$$ Φ^{-1}(y_1) = z_1 + y_2 + z_3 $$
with $z_1$ and $z_3$ being exogenous control variables and $z_2$ a valid instrument for $y_2$.
I cannot do the standard heckit model as not only my first but also my second stage is probit. "Forbidden regression" (i.e. using standard 2SLS) is also not a great way out.
Does anybody know how to fit such a model? If there is a solution, I am particular interested in a step-by-step explanation which can be implemented relatively easy in R (or Stata).