Deriving likelihood function for IV-probit

Question

So I have a binary model where $y_1^*$ is the latent unobserved variable and $y_1 \in \{0,1\}$ the observed. $y_2$ determines $y_1$ and $z_2$ is thus my instrument. So in short the model is. \begin{eqnarray} y_1^*&=& \delta_1 z_1 + \alpha_1 y_2 + u_1 \\ y_2 &=& \delta_{21} z_1 + \delta_{22}z_2 + v_2 = \textbf{z}\delta + v_2 \\ y_1 &=& \text{1}[y^*>0] \end{eqnarray} Since the error terms aren't independent but, \begin{eqnarray} \begin{pmatrix} u_1 \\ v_2 \end{pmatrix} \sim \mathcal{N} \left(\textbf{0} \; , \begin{bmatrix} 1 &\eta \\ \eta &\tau^2 \end{bmatrix} \right). \nonumber \end{eqnarray} I make use of an IV-probit model.

I'm having trouble deriving the likelihood function. I get that I can write one of the error terms as a linear function of the other so, \begin{eqnarray} u_{1} = \frac{\eta}{\tau^2}v_{2} + \xi, \qquad \text{where} \quad \xi \sim \mathcal{N}(0, 1-\eta^2). \end{eqnarray}
and that $\xi$ should be used in order to impose a normal CDF.

I've looked in the Stata manual (http://www.stata.com/manuals13/rivprobit.pdf) for IV-probit and they suggest using the definition of the conditional density \begin{eqnarray} f(y_1, y_2 \mid \textbf{z}) = f(y_1 \mid y_2, \textbf{z}) f(y_2 \mid \textbf{z}) \end{eqnarray}
in order to derive the likelihood function but I really don't use it (and yea I end up with the wrong result...). My attempt so far is,

\begin{eqnarray} \mathcal{L}(y_1) &=& \prod_{i=1}^n \Pr(y_1=0 \mid y_2, \textbf{z} )^{1-y_1} \Pr(y_1=1 \mid y_2, \textbf{z} )^{y_1} \nonumber \\ &=& \prod_{i=1}^n \Pr(y_1^* \leq 0)^{1-y_1} \Big(\Pr(y_1^* > 0) f(y_2 \mid \textbf{z}) \Big)^{y_1} \nonumber \\ \text{[standardizing]} &=& \prod_{i=1}^n \Pr \Big( \frac{\xi}{\sqrt{1-\eta^2}} \leq - \frac{\delta_1 z_1 + \alpha_1 y_2 + \frac{\eta}{\tau^2}(y_2 - \textbf{z})}{\sqrt{1-\eta^2}}\Big)^{1-y_1} \\ &\cdot& \Big(\Pr \Big( \frac{\xi}{\sqrt{1-\eta^2}} < \frac{\delta_1 z_1 + \alpha_1 y_2 + \frac{\eta}{\tau^2}(y_2 - \textbf{z})}{\sqrt{1-\eta^2}}\Big) f(y_2 \mid \textbf{z}) \Big)^{y_1} \nonumber \\ &=& [1-\Phi(w)]^{1-y_i} \left[ \Phi(w)f(y_2 \mid \textbf{x}) \right]^{y_1} \end{eqnarray} As I said, I haven't used the definition for the joint density function as stated above. Moreover, I end up with also $f(y_2 \mid \textbf{z})$ being raised to $y_1$ which appears to be wrong. Can somebody give me a hint on how to derive the correct (log-) likelihood function or where I went wrong?

Answer 1

Remember that for a bivariate normal variable $$\begin{pmatrix}X \\ Y\end{pmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mu_X\\\mu_Y\end{bmatrix}, \begin{bmatrix}\sigma_X^2 & \rho\sigma_X\sigma_Y\\\rho\sigma_X\sigma_Y & \sigma_Y^2\end{bmatrix}\right),$$ the conditional distribution of $Y$ given $X$ is $$Y\mid X \sim \mathcal{N}\left(\mu_Y+\rho\sigma_Y\frac{X-\mu_X}{\sigma_X},\sigma_Y\left[1-\rho^2\right]\right).$$

In the present case, we have \begin{align} u_1 \mid v_2 &\sim \mathcal{N}\left(0+\frac{\eta}{1\cdot\tau}\cdot1\frac{v_2-0}{\tau}, 1\cdot\left[1-\left(\frac{\eta}{1\cdot\tau}\right)^2\right] \right) \\ &= \mathcal{N}\left(\frac{\eta}{\tau^2}v_2, 1-\frac{\eta^2}{\tau^2} \right), \end{align} which means that $$u_1=\frac{\eta}{\tau^2}v_2+\xi$$ where (and this was your first mistake) $$\xi\sim\mathcal{N}\left(0,1-\frac{\eta^2}{\tau^2}\right).$$

We can thus rewrite the first equation \begin{align} y_1^* &= \delta_1 z_1 + \alpha_1 y_2 + u_1 \\ &= \delta_1 z_1 + \alpha_1 y_2 + \frac{\eta}{\tau^2}v_2+\xi \\ &= \delta_1 z_1 + \alpha_1 y_2 + \frac{\eta}{\tau^2}(y_2-\textbf{z}\delta)+\xi. \end{align}

Now, remember that the conditional probability density function of $X=x$ given $Y=y$ is $$f_{X}(x \mid y)=\frac{f_{XY}(x,y)}{f_{Y}(y)}.$$

In the present case, we have $$f_{1}(y_1 \mid y_2, \mathbf{z})=\frac{f_{12}(y_1,y_2 \mid \mathbf{z})}{f_{2}(y_2 \mid \mathbf{z})},$$ which can be rearranged to your expression $$f_{12}(y_1, y_2 \mid \mathbf{z})= f_{1}(y_1 \mid y_2, \mathbf{z})f_{2}(y_2 \mid \mathbf{z}).$$

Then, we can write the likelihood as a function of the densities of the two independent shocks $v_1,\xi_1$: \begin{align} \mathcal{L}(y_1,y_2\mid \mathbf{z}) &= \prod_i^n f_{1}(y_{1i} \mid y_{2i}, \mathbf{z}_i)f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(y_{1i}=1\right)^{y_{1i}}\Pr\left(y_{1i}=0\right)^{1-y_{1i}}f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(y_{1i}^*>0\right)^{y_{1i}}\Pr\left(y_{1i}^*\leq0\right)^{1-y_{1i}}f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_{i}\delta)+\xi_i>0\right)^{y_{1i}}\\ &\qquad\quad \Pr\left(\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)+\xi_i\leq0\right)^{1-y_{1i}}\\ &\qquad\quad f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(\xi_i>-\left[\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)\right]\right)^{y_{1i}}\\ &\qquad\quad \Pr\left(\xi_i\leq-\left[\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)\right]\right)^{1-y_{1i}}\\ &\qquad\quad f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(\frac{\xi_i-0}{\sqrt{1-\frac{\eta^2}{\tau^2}}}>-\frac{\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)+0}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\right)^{y_{1i}}\\ &\qquad\quad \Pr\left(\frac{\xi_i-0}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\leq-\frac{\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)+0}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\right)^{1-y_{1i}}\\ &\qquad\quad f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \Pr\left(\frac{\xi_i}{\sqrt{1-\frac{\eta^2}{\tau^2}}}>-w_i\right)^{y_{1i}} \Pr\left(\frac{\xi_i}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\leq-w_i\right)^{1-y_{1i}} f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i^n \left[1-\Pr\left(\frac{\xi_i}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\leq-w_i\right)\right]^{y_{1i}} \Pr\left(\frac{\xi_i}{\sqrt{1-\frac{\eta^2}{\tau^2}}}\leq-w_i\right)^{1-y_{1i}} f_{2}(y_{2i} \mid \mathbf{z}_i) \\ &= \prod_i \left[1-\Phi(-w_i)\right]^{y_{1i}} \Phi(-w_i)^{1-y_{1i}} \varphi\left(\frac{y_{2i}-\mathbf{z}_i\delta}{\tau}\right) \\ &= \prod_i^n \Phi(w_i)^{y_{1i}} \left[1-\Phi(w_i)\right]^{1-y_{1i}} \varphi\left(\frac{y_{2i}-\mathbf{z}_i\delta}{\tau}\right) \\ &= \Phi(w)^{y_{1}} \left[1-\Phi(w)\right]^{1-y_{1}} \varphi\left(\frac{y_{2}-\mathbf{z}\delta}{\tau}\right) \\ \end{align} where \begin{align} w_i = \frac{\delta_1 z_{1i} + \alpha_1 y_{2i} + \frac{\eta}{\tau^2}(y_{2i}-\textbf{z}_i\delta)}{\sqrt{1-\frac{\eta^2}{\tau^2}}}. \end{align} $\Phi(z)$ and $\varphi(z)$ are the cumulative density function and probability density function of the standard normal distribution.