Recursive Bivariate Probit With An Endogenous Binary Regressor

Why Model the Two Binary Decisions Jointly?

Suppose an outcome is binary and one of its regressors is also binary, but that regressor is endogenous. For example, an analyst might study whether a patient experiences an outcome and whether the patient receives a treatment. If unobserved factors that affect treatment also affect the outcome, a probit regression that treats treatment as exogenous does not represent the assumed data-generating process.

A recursive bivariate probit addresses this dependence by specifying a latent-variable equation for each binary variable and allowing the two latent errors to be correlated. “Recursive” means the outcome equation contains the realized binary regressor, while the regressor equation does not contain the outcome. The model’s joint likelihood uses the observed outcome and regressor together rather than substituting a fitted value from a separate first stage.

Definition: For observation $i$, let $Y_i$ be a binary outcome, $D_i$ an endogenous binary regressor, $x_i$ a vector of outcome covariates, and $z_i$ a vector of regressor-equation covariates. The recursive bivariate probit is $$ Y_i=\mathbf{1}\{x_i'\beta+\delta D_i+\varepsilon_i>0\}, \qquad D_i=\mathbf{1}\{z_i'\gamma+u_i>0\}, $$ where $\mathbf{1}\{\cdot\}$ equals one when its condition is true and zero otherwise. Conditional on the covariates, assume $$ \begin{pmatrix}\varepsilon_i\\u_i\end{pmatrix} \sim N\left( \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}1&\rho\\\rho&1\end{pmatrix} \right), \qquad -1<\rho<1. $$

The coefficient $\delta$ shifts the outcome’s latent index when $D_i$ changes from zero to one. It is not, by itself, a probability-point effect. The correlation $\rho$ captures dependence between the two latent errors after conditioning on the specified covariates. This bivariate normal specification fixes both latent-error variances at one, as is conventional in probit models.

Assumptions and Exclusion

The model’s likelihood has a structural interpretation only under assumptions about the equations, errors, and covariates. Observations are assumed independent across $i$, conditional on the model specification. For each observation, the latent errors have the stated standard bivariate normal distribution conditional on $x_i,z_i$, with a common correlation $\rho$. The treatment equation is correctly specified as a probit index, and the outcome equation is correctly specified as a probit index that includes $D_i$.

The covariates $z_i$ should include at least one excluded instrument: a variable that shifts the probability of $D_i=1$, conditional on the included covariates, but does not enter the outcome equation directly. If some covariates are common to both equations, their coefficients are still equation-specific. The exclusion restriction is that the excluded variable has no direct effect on $Y_i$ other than through $D_i$. Its relevance means it provides variation in the treatment index after conditioning on the other regressors. The model also requires the instrument to be appropriately exogenous to the structural errors; this is a substantive assumption, not something established simply by fitting the likelihood.

Exclusion and relevance matter for distinguishing the treatment effect from latent-error dependence. Without useful excluded variation, estimates can rely heavily on the assumed bivariate normal form, and the separation of $\delta$ from $\rho$ may be weak. A fitted joint model does not, on its own, validate an instrument or make an observational association causal.

Deriving the Four Joint Probabilities

Write the outcome index at the observed treatment value as $m_i=x_i'\beta+\delta D_i$, and the treatment index as $q_i=z_i'\gamma$. Let $\Phi$ denote the standard normal cumulative distribution function, and let $\Phi_2(a,b;\rho)$ denote the cumulative distribution function of a standard bivariate normal pair with correlation $\rho$, evaluated at $(a,b)$. The observed pair $(Y_i,D_i)$ has four possible values.

Both indicators equal one.
Here $\varepsilon_i>-m_i$ and $u_i>-q_i$. By symmetry of the centered bivariate normal distribution, this probability is $\Phi_2(m_i,q_i;\rho)$.

Outcome zero, regressor one.
Start with $P(D_i=1)=\Phi(q_i)$, then subtract the probability that both latent variables exceed their outcome and treatment thresholds. The result is $\Phi(q_i)-\Phi_2(m_i,q_i;\rho)$.

Outcome one, regressor zero.
Start with $P(D_i=0)=\Phi(-q_i)$, then subtract the event $Y_i=0,D_i=0$. The result is $\Phi(-q_i)-\Phi_2(-m_i,-q_i;\rho)$.

Both indicators equal zero.
The event is $\varepsilon_i\le -m_i$ and $u_i\le -q_i$, so its probability is $\Phi_2(-m_i,-q_i;\rho)$.

In steps 2 and 3, $m_i$ is defined using the observed $D_i$. Thus, each likelihood contribution evaluates the outcome index at the treatment value actually observed for that individual. Write $p_{yd,i}$ for the probability of observing $Y_i=y,D_i=d$. Under independent observations, the likelihood and log-likelihood are

L(\beta,\delta,\gamma,\rho) =\prod_{i=1}^n p_{Y_iD_i,i}, \qquad \ell(\beta,\delta,\gamma,\rho) =\sum_{i=1}^n\log p_{Y_iD_i,i}.

The parameters are estimated jointly by maximizing this likelihood. In practice, evaluating the bivariate normal probabilities and optimizing the resulting objective are numerical tasks. The four probabilities are not four unrelated regressions: they are linked through the same indices and the same correlation parameter.

Worked Example: A Joint Likelihood Contribution

Consider a hypothetical observation with $Y_i=1,D_i=1$. At a trial parameter value, suppose $m_i=0$, $q_i=0$, and $\rho=0.4$. Calculate this observation’s likelihood contribution.

For zero thresholds, the standard bivariate normal probability has the identity $$ \Phi_2(0,0;\rho)=\frac14+\frac{\arcsin(\rho)}{2\pi}. $$ Therefore, $$ p_{11,i}=\Phi_2(0,0;0.4) =\frac14+\frac{\arcsin(0.4)}{2\pi} \approx 0.3155. $$

Since the observed pair is $(1,1)$, its likelihood contribution is $p_{11,i}\approx0.3155$, and its log-likelihood contribution is $\log(0.3155)\approx-1.154$. This is the contribution at the stated trial parameter values; it is not a fitted estimate based on one observation.

Worked Example: Log-Likelihood for Hypothetical Counts

Suppose ten hypothetical observations all have zero outcome and treatment indices at a particular trial parameter value, with $\rho=0.4$. The observed counts are $n_{11}=4,n_{10}=2,n_{01}=3,n_{00}=1$, where the first index is $Y$ and the second is $D$. Calculate the log-likelihood at these trial values.

The preceding calculation gives $p_{11}=0.3155$. By symmetry, $p_{00}=0.3155$. The marginal probability of either treatment value is $0.5$, so each off-diagonal cell has probability $0.5-0.3155=0.1845$. The probabilities sum to one: $2(0.3155)+2(0.1845)=1.000$.

Using the counts, $$ \ell =4\log(0.3155)+2\log(0.1845)+3\log(0.1845)+\log(0.3155) =5\log(0.3155)+5\log(0.1845) \approx -14.22. $$

The likelihood is approximately $\exp(-14.22)$, or $6.70\times10^{-7}$. The numerical likelihood is small because it multiplies ten probabilities; its absolute size is not usually interpreted on its own. Estimation compares the log-likelihood across parameter values.

Why the Correlation Signals Endogeneity

When $\rho\ne0$, individuals selected into $D_i=1$ tend to have a different distribution of $\varepsilon_i$ from individuals with $D_i=0$, even after conditioning on the covariates. This is the model’s representation of conditional dependence between the treatment and outcome errors. The connection can be shown directly from the treatment threshold.

Given $z_i$, $D_i=1$ is equivalent to $u_i>-q_i$. For a standard normal variable, the mean above this threshold is $\phi(q_i)/\Phi(q_i)$, where $\phi$ is the standard normal density. Since the conditional mean of $\varepsilon_i$ given $u_i$ is $\rho u_i$, averaging over the selected values of $u_i$ gives

E[\varepsilon_i\mid D_i=1,x_i,z_i] =\rho\frac{\phi(q_i)}{\Phi(q_i)}, \qquad E[\varepsilon_i\mid D_i=0,x_i,z_i] =-\rho\frac{\phi(q_i)}{\Phi(-q_i)}.

For example, the first equality follows by iterated expectation: $E[\varepsilon_i\mid D_i=1,z_i]=E[E[\varepsilon_i\mid u_i]\mid u_i>-q_i]$. The inner conditional expectation is $\rho u_i$, and the mean of a standard normal variable truncated below at $-q_i$ is $\phi(q_i)/\Phi(q_i)$. The second equality uses the corresponding mean below the threshold.

If $\rho=0$, the bivariate normal errors are independent, so selection into $D_i$ does not change the conditional distribution of $\varepsilon_i$. If $\rho\ne0$, the selected groups differ in their outcome errors, conditional on the specified covariates. Thus, under the model, testing $H_0:\rho=0$ is a test of this latent-error dependence. It is not a direct test of whether an instrument is valid, and a nonzero estimate does not establish which omitted factor caused the dependence.

Worked Example: Testing for Latent-Error Dependence

Suppose a hypothetical joint maximum-likelihood fit gives $\hat{\rho}=0.36$ with $SE(\hat{\rho})=0.14$. Test $H_0:\rho=0$ against $H_a:\rho\ne0$ using a large-sample Wald test at $\alpha=0.05$.

The Wald statistic is $$ z=\frac{\hat{\rho}-0}{SE(\hat{\rho})} =\frac{0.36}{0.14} \approx 2.571. $$ Under the null and the large-sample approximation, this statistic is compared with a standard normal distribution. Its two-sided p-value is $$ 2P(Z\ge 2.571)\approx 0.01013, \qquad Z\sim N(0,1). $$

Because $0.01013<0.05$, reject $H_0$ at the five-percent level. Under the specified model, the data provide evidence that the latent errors are correlated after conditioning on the included covariates. This conclusion depends on the model and its assumptions; it does not verify the exclusion restriction or identify a particular source of endogeneity.

Interpreting the Outcome Coefficient

The coefficient $\delta$ shifts the latent outcome index, not the observed outcome probability by a constant amount. When $\rho=0$, the outcome error is independent of treatment selection under the model. Conditional on $x_i$ and $D_i=d$, the outcome probability is then $\Phi(x_i'\beta+\delta d)$. When $\rho\ne0$, conditioning on the observed treatment also conditions on a truncated part of $u_i$, and hence changes the distribution of $\varepsilon_i$. The conditional outcome probability must then account for the joint distribution; the simple single-index expression does not apply.

Worked Example: A Latent-Index Shift Is Not a Probability-Point Effect

Suppose $\rho=0$, and for a hypothetical individual the outcome index without treatment is $x_i'\beta=-0.2$. Let $\delta=0.5$. Calculate the two model-implied outcome probabilities and their difference under the zero-correlation case.

Without treatment, the probability is $\Phi(-0.2)\approx0.4207$. With treatment, the index is $-0.2+0.5=0.3$, so the probability is $\Phi(0.3)\approx0.6179$. The difference is $0.6179-0.4207=0.1972$, or about 19.72 percentage points for this index value under the stated assumptions.

The change is not equal to $\delta=0.5$: that coefficient is measured in latent-index units. Nor is 19.72 percentage points a universal effect for every covariate value. If $\rho\ne0$, these calculations would also fail to account for selection-induced dependence.

Common Pitfalls and Practical Interpretation

Treating $\rho$ as a causal effect. It is a correlation between latent errors under the model, not the effect of treatment on the outcome. The treatment-index coefficient $\delta$ and the correlation parameter answer different questions.
Omitting the exclusion restriction from the argument. Joint normality is a distributional assumption, not a substitute for credible excluded variation. Relevance and the exclusion restriction need substantive justification.
Reading $\delta$ as a probability difference. It shifts a latent index. Probability changes depend on the index and, when $\rho\ne0$, on the selection structure.
Using a product of separate marginal likelihoods. When $\rho\ne0$, the joint probabilities are not products of the outcome and treatment marginal probabilities. The joint likelihood must use the four bivariate-normal cell probabilities.
Overstating a test of $\rho=0$. Rejecting that null is evidence of latent-error dependence under the specified model. It neither proves the instrument valid nor rules out misspecification.

Key takeaway: The recursive bivariate probit models a binary outcome and an endogenous binary regressor through two latent indices with correlated standard normal errors. Its joint likelihood respects that dependence; a nonzero $\rho$ signals conditional latent-error dependence under the assumptions, while credible exclusion and relevance remain essential for interpretation.

Check Your Understanding

Try these questions without looking back.

In the recursive specification, which equation contains the observed binary regressor, and which equation contains the excluded instrument?
For $Y_i=0,D_i=0$, write the joint probability in terms of $\Phi_2$, $m_i$, and $q_i$.
Explain why $\rho\ne0$ makes the outcome error distribution differ between the $D_i=1$ and $D_i=0$ groups.
What does the exclusion restriction require of a variable that appears in the treatment equation but not the outcome equation?
If $\hat{\rho}=0.20$ and $SE(\hat{\rho})=0.10$, calculate the Wald statistic for testing $H_0:\rho=0$. What additional assumption is needed to interpret its reference p-value?
Why is $\delta$ not generally equal to the change in outcome probability when treatment changes from zero to one?