Multivariate probit model

In the ordinary probit model, there is only one binary dependent variable ${\displaystyle Y}$ and so only one latent variable ${\displaystyle Y^{*}}$ is used. In contrast, in the bivariate probit model there are two binary dependent variables ${\displaystyle Y_{1}}$ and ${\displaystyle Y_{2}}$, so there are two latent variables: ${\displaystyle Y_{1}^{*}}$ and ${\displaystyle Y_{2}^{*}}$. It is assumed that each observed variable takes on the value 1 if and only if its underlying continuous latent variable takes on a positive value:

${\displaystyle Y_{1}={\begin{cases}1&{\text{if }}Y_{1}^{*}>0,\\0&{\text{otherwise}},\end{cases}}}$

${\displaystyle Y_{2}={\begin{cases}1&{\text{if }}Y_{2}^{*}>0,\\0&{\text{otherwise}},\end{cases}}}$

with

${\displaystyle {\begin{cases}Y_{1}^{*}=X_{1}\beta _{1}+\varepsilon _{1}\\Y_{2}^{*}=X_{2}\beta _{2}+\varepsilon _{2}\end{cases}}}$

and

${\displaystyle {\begin{bmatrix}\varepsilon _{1}\\\varepsilon _{2}\end{bmatrix}}\mid X\sim {\mathcal {N}}\left({\begin{bmatrix}0\\0\end{bmatrix}},{\begin{bmatrix}1&\rho \\\rho &1\end{bmatrix}}\right)}$

Fitting the bivariate probit model involves estimating the values of ${\displaystyle \beta _{1},\ \beta _{2},}$ and ${\displaystyle \rho }$. To do so, the likelihood of the model has to be maximized. This likelihood is

${\displaystyle {\begin{aligned}L(\beta _{1},\beta _{2})={\Big (}\prod &P(Y_{1}=1,Y_{2}=1\mid \beta _{1},\beta _{2})^{Y_{1}Y_{2}}P(Y_{1}=0,Y_{2}=1\mid \beta _{1},\beta _{2})^{(1-Y_{1})Y_{2}}\\[8pt]&{}\qquad P(Y_{1}=1,Y_{2}=0\mid \beta _{1},\beta _{2})^{Y_{1}(1-Y_{2})}P(Y_{1}=0,Y_{2}=0\mid \beta _{1},\beta _{2})^{(1-Y_{1})(1-Y_{2})}{\Big )}\end{aligned}}}$

Substituting the latent variables ${\displaystyle Y_{1}^{*}}$ and ${\displaystyle Y_{2}^{*}}$ in the probability functions and taking logs gives

${\displaystyle {\begin{aligned}\sum &{\Big (}Y_{1}Y_{2}\ln P(\varepsilon _{1}>-X_{1}\beta _{1},\varepsilon _{2}>-X_{2}\beta _{2})\\[4pt]&{}\quad {}+(1-Y_{1})Y_{2}\ln P(\varepsilon _{1}<-X_{1}\beta _{1},\varepsilon _{2}>-X_{2}\beta _{2})\\[4pt]&{}\quad {}+Y_{1}(1-Y_{2})\ln P(\varepsilon _{1}>-X_{1}\beta _{1},\varepsilon _{2}<-X_{2}\beta _{2})\\[4pt]&{}\quad {}+(1-Y_{1})(1-Y_{2})\ln P(\varepsilon _{1}<-X_{1}\beta _{1},\varepsilon _{2}<-X_{2}\beta _{2}){\Big )}.\end{aligned}}}$

After some rewriting, the log-likelihood function becomes:

${\displaystyle {\begin{aligned}\sum &{\Big (}Y_{1}Y_{2}\ln \Phi (X_{1}\beta _{1},X_{2}\beta _{2},\rho )\\[4pt]&{}\quad {}+(1-Y_{1})Y_{2}\ln \Phi (-X_{1}\beta _{1},X_{2}\beta _{2},-\rho )\\[4pt]&{}\quad {}+Y_{1}(1-Y_{2})\ln \Phi (X_{1}\beta _{1},-X_{2}\beta _{2},-\rho )\\[4pt]&{}\quad {}+(1-Y_{1})(1-Y_{2})\ln \Phi (-X_{1}\beta _{1},-X_{2}\beta _{2},\rho ){\Big )}.\end{aligned}}}$

Note that ${\displaystyle \Phi }$ is the cumulative distribution function of the bivariate normal distribution. ${\displaystyle Y_{1}}$ and ${\displaystyle Y_{2}}$ in the log-likelihood function are observed variables being equal to one or zero.

For the general case, ${\displaystyle \mathbf {y_{i}} =(y_{1},...,y_{j}),\ (i=1,...,N)}$ where we can take ${\displaystyle j}$ as choices and ${\displaystyle i}$ as individuals or observations, the probability of observing choice ${\displaystyle \mathbf {y_{i}} }$ is

${\displaystyle {\begin{aligned}\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int _{A_{J}}\cdots \int _{A_{1}}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )dy_{1}^{*}\dots dy_{J}^{*}\\\Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )=&\int \mathbb {1} _{y^{*}\in A}f_{N}(\mathbf {y} _{i}^{*}|\mathbf {X_{i}\beta } ,\Sigma )d\mathbf {y} _{i}^{*}\end{aligned}}}$

Where ${\displaystyle A=A_{1}\times \cdots \times A_{J}}$ and,

${\displaystyle A_{j}={\begin{cases}(-\infty ,0]&y_{j}=0\\(0,\infty )&y_{j}=1\end{cases}}}$

The log-likelihood function in this case would be ${\displaystyle \sum _{i=1}^{N}\log \Pr(\mathbf {y_{i}} |\mathbf {X_{i}\beta } ,\Sigma )}$

Except for ${\displaystyle J\leq 2}$ typically there is no closed form solution to the integrals in the log-likelihood equation. Instead simulation methods can be used to simulated the choice probabilities. Methods using importance sampling include the GHK algorithm (Geweke, Hajivassilou, McFadden and Keane),[3] AR (accept-reject), Stern's method. There are also MCMC approaches to this problem including CRB (Chib's method with Rao-Blackwellization), CRT (Chib, Ritter, Tanner), ARK (accept-reject kernel), and ASK (adaptive sampling kernel).[4] A variational approach scaling to large datasets is proposed in Probit-LMM (Mandt, Wenzel, Nakajima et al.).[5]

• Greene, William H., Econometric Analysis, seventh edition, Prentice-Hall, 2012.