Simultaneous Equations Model with Non-Linear and Linear Dependent Variables on Panel Data

Amélie Adeline¹, Richard K. Moussa²

¹ThEMA, University of Cergy-Pontoise, Cergy-Pontoise, France.
²Ecole Nationale Supérieure de Statistique et d’Economie Appliquée, Abidjan, Cote d’Ivoire.

DOI: 10.4236/tel.[phone]                

Abstract

This paper provides an estimation approach for the multi-equations’ systems in panel data. Multi-equations systems are at the heart of economic modeling. Researchers who want to establish causal links between two outcomes, often need to consider simultaneity between the latter, to overcome endogeneity issues (for instance when considering supply and demand equations). Difficulties arise when considering linear and non-linear outcomes at the same time and this is why Roodman [1] implemented the Stata module cmp for multidimensional models. In this paper, we further develop this technique to allow researchers to implement a simultaneous equations model in a panel dimension setting. Implemented under Stata, our method, xtcmp, is a Full Information Maximum Likelihood (FIML) estimator. This paper explains the associated theory (derivation of the log-likelihood function, the associated gradient and the Hessian matrices of the log-integrand function) and offers an application of t xtcmp, while making comparisons with cmp.

Keywords

Causality, Full Information Maximum Likelihood, Panel Data, Gauss-Hermite Quadrature, Gradient, Hessian

1. Introduction

In empirical economics, a common approach is to consider a linear data-generating process. However, non-linear outcomes are often present and important in research questions. This is due to the structure of the database where interviewers transcribe yes-no answers into binary outcomes. However, when researchers point out a project, they often have to take into account different variables, continuous and categorical, at the same time, while considering simultaneous equations framework in a dynamic setting, in which each dependent variable is endogenous in one (or more) equation(s) of the model. The advantage of simultaneous equations model is to consider the correlation between the error terms of each equation. More specifically, in a dynamic setup, such models allow researchers to consider different individual effects (which are part of error terms, the latter being decomposed into an individual effect fixed across time and an effect which depends on time) across equations. This is of importance since these terms are unobserved, specific to each outcome, and might imply endogeneity issues. For instance in health economics, when investigating causal relationship between health and income, this can run from income to health and from health to income such that both are endogenous to each other [2]. In this way, considering a dynamic simultaneous equation model allows to consider unobserved individual effects such as physical maturity (thanks to genetics for health) or intellectual abilities (for income).

The framework of multi-equations models has been widely used in the literature to address several issues including the case of an endogenous binary outcome. Greene [3] has reformulated the estimation of the impact of an endogenous treatment on a continuous outcome using a multi-equations model. This reformulation has been extended to the analysis of endogenous binary outcome by [4]. Thus, several papers have analyzed the effects of an endogenous treatment on diverse types of outcomes including continuous and count outcomes [5] [6] [7]. Also, some generalizations for the case of noncompliance and nonresponse have been introduced by [8]. However, all these methods focus on cross sectional data and do not account for panel data.

There is almost no automated estimation method in Stata software to estimate parameters in this multi-equations model. An exception is the cmp command which is the first general Stata tool for this class of models, and is written as a Seemingly Unrelated Regressions (SUR) estimator [1]. However, this command does not explicitly consider panel dimension of the data which might be an issue due to the effervescence of databases with a temporal dimension. Moreover, simple relationships among variables at a point in time do not capture adequately the dynamic interaction of changing humans in changing environments. Thus, there is a need to develop a command for simultaneous equations model for panel database.

As a result, we offer an extension of the cmp framework, in a case where there are either two equations (one linear and one binary outcomes) or three equations (either two linear dependent variables with one binary, or one linear and two binary dependent variables), while explicitly considering the panel dimension of the data. In this way, our command xtcmp is a Full Information Maximum Likelihood (FIML) estimator, taking into account time dimension of the data, as well as, linear and non-linear outcomes (which is not feasible with three-stage least squares, because the latter only takes into account linear dependent variables).

As a result, the likelihood function is a multidimensional integral, such that we use the adaptive Gauss-Hermite quadrature method as an approximation (as proposed by Liu and Pierce [9]). For the accuracy of the method and to reduce computing time, we derive the gradient of the log-likelihood, and the Hessian of the respective integrand. The estimation method of xtcmp has been implemented using the d1 method of Stata software (see [10] for further details).

Section 2 derives the likelihood for a FIML estimator in a general setting, as well as, for three specific cases. Section 3 discusses the estimation requirements needed, and details the Hessian matrices, as well as the gradient vector according to the parameters. In Section 4, we give some examples on the use of xtcmp, while making comparisons with Roodman’s [1] cmp results. Section 5 concludes the paper.

2. Likelihood for a FIML Estimator

Let Y denotes a d-dimensional vector of endogenous variables in a simultaneous equations model. Roodman [1] specifically discusses conditions for consistency and identification of such model. Let $Y^k$ denotes the k^th component of Y such that the value of $Y^k$ for individual i at period t is given by $y_{it}^k$. We assume that the first $d_1$ components of Y, where $d_1<d$, are binary outcomes, and the others, $d_2=d-d_1$, are continuous. Let $ϵ$ denotes the vector of associated error terms. By assuming a panel random effects model, $ϵ$ can be decomposed into two terms such that $ϵ=\mu +\nu$, where $\mu$ is time-invariant. In this way, the error term for individual i at period t in the k^th equation is given by ${ϵ}_{it}^k={\mu }_i^k+{\nu }_{it}^k$. The full model can be written as follows:

$\tilde{Y}=X\beta +ϵ$ (1)

where $\tilde{Y}$ contains the related latent variables for the first $d_1$ equation and the original continuous variables for the others. The explanatory matrix X is given by $X=diag\left(Z^k\right)$ where $Z^k$, with $k=1,\cdots ,d$, corresponds to the explanatory variables for the k^th equation. Similarly, the parameter vector is $\beta ={\left({\beta }_1\mathrm{,}\cdots \mathrm{,}{\beta }_d\right)}^{\prime }$ where ${\beta }^k$ is the parameter vector of the k^th equation. We suppose that the classical hypotheses on independence between 1) error components, and 2) the error components and explanatory variables, are satisfied. Furthermore, let us assume that the error components are independent and identically distributed with zero means and covariance matrices ${\Sigma }_{\mu }$ and ${\Sigma }_{\nu }$, the latter being defined as follows:

${\Sigma }_{\nu }=\left(\begin{array}{cc}{\Sigma }_1& {{\Sigma }^{\prime }}_3\\ {\Sigma }_3& {\Sigma }_2\end{array}\right)$ (2)

in which ${\Sigma }_1$ is a $d_1$ dimension matrix with 1 over the diagonal (which corresponds to the covariance matrix structure for simultaneous equations with only binary outcomes), ${\Sigma }_2$ is a $d_2$ dimension matrix, and ${\Sigma }_3$ is a $d_2\ast d_1$ dimension matrix. Thus, the overall individual likelihood is given by:

$L_i={\displaystyle {\int }_{{ℝ}^d}}\left\{{\displaystyle \underset{t=1}{\overset{T_i}{\prod }}}{f}_{\nu }\left({\nu }_{it}^1\mathrm{,}\cdots \mathrm{,}{\nu }_{it}^d\mathrm{<mo>\; |\; </mo>}{\mu }_i^1\mathrm{,}\cdots \mathrm{,}{\mu }_i^d\right)\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}\cdots \mathrm{,}{\mu }_i^d\right)\text{d}{\mu }_i^1\cdots \text{d}{\mu }_i^d$ (3)

where $f_{\mu }\left({\mu }_i^1\mathrm{,}\cdots \mathrm{,}{\mu }_i^d\right)=\frac{1}{{\left(2\pi \right)}^{d/2}\sqrt{\left(\det \left({\Sigma }_{\mu }\right)\right)}}\exp \left(\frac{-1}{2}{\mu }^{\prime }{\Sigma }_{\mu }^{-1}\mu \right)$. The $d_1$ first equations being related to binary outcomes, if we define $q_{it}^k=2\ast y_{it}^k-1$, then, the density function $f_{\nu }\left({\nu }_{it}^1\mathrm{,}\cdots \mathrm{,}{\nu }_{it}^d\mathrm{<mo>\; |\; </mo>}{\mu }_i^1\mathrm{,}\cdots \mathrm{,}{\mu }_i^d\right)$ is given by The density function is denoted by ${\mathcal{l}}_{it}$ for simplification:

${\mathcal{l}}_{it}=f_{\nu }\left({\nu }_{it}^1,\cdots ,{\nu }_{it}^d|{\mu }_i^1,\cdots ,{\mu }_i^d\right)={\displaystyle \underset{-\infty }{\overset{q_{it}^1{z}_{it}^1}{\int }}}\cdots {\displaystyle \underset{-\infty }{\overset{q_{it}^{d_1}{z}_{it}^{d_1}}{\int }}}{\varphi }_d\left({\nu }_{it}^1,\cdots ,{\nu }_{it}^d\right)\text{d}{\nu }_i^1\cdots \text{d}{\nu }_i^{d_1}$ (4)

where $z_{it}^k=Z_{it}^k{\beta }^k+{\mu }_i^k$, with $k=1,\cdots ,d$ and ${\nu }_{it}^k={\tilde{y}}_{it}^k-z_{it}^k$, with $k=d_1+1,\cdots ,d$. The idiosyncratic error is $\nu =\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\mathrm{,}{\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)\sim \mathcal{N}\left(\mathrm{0,}{\Sigma }_{\nu }\right)$, such that $\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)\sim \mathcal{N}\left(\mathrm{0,}{\Sigma }_2\right)$ and $\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\right)\mathrm{<mo>\; |\; </mo>}\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)\sim \mathcal{N}\left(m_{\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\right)\mathrm{<mo>\; |\; </mo>}\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)}\mathrm{,}{\Sigma }_{\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\right)\mathrm{<mo>\; |\; </mo>}\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)}\right)$ with:

$m_{d_1|d_2}=m_{\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\right)\mathrm{<mo>\; |\; </mo>}\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)}={{\Sigma }^{\prime }}_3{\left({\Sigma }_2\right)}^{-1}{\left({\nu }_{d_1+1},\cdots ,{\nu }_d\right)}^{\prime }$ (5)

and

${\Sigma }_{d_1|d_2}={\Sigma }_{\left({\nu }_1\mathrm{,}\cdots \mathrm{,}{\nu }_{d_1}\right)\mathrm{<mo>\; |\; </mo>}\left({\nu }_{d_1+1}\mathrm{,}\cdots \mathrm{,}{\nu }_d\right)}={\Sigma }_1-{{\Sigma }^{\prime }}_3{\left({\Sigma }_2\right)}^{-1}{\Sigma }_3$ (6)

Thus, we have:

${\mathcal{l}}_{it}={\varphi }_{d_2}\left({\nu }_{it}^{d_1+1},\cdots ,{\nu }_{it}^d\right){\Phi }_{d_1}\left(\left(q_{it}^1{z}_{it}^1,\cdots ,q_{it}^{d_1}{z}_{it}^{d_1}\right),m_{d_1|d_2},{\Sigma }_{d_1|d_2}\right)$ (7)

in which ${\Phi }_{d_1}\left(\left(q_{it}^1{z}_{it}^1\mathrm{,}\cdots \mathrm{,}{q}_{it}^{d_1}{z}_{it}^{d_1}\right)\mathrm{,}{m}_{d_1\mathrm{<mo>\; |\; </mo>}{d}_2}\mathrm{,}{\Sigma }_{d_1\mathrm{<mo>\; |\; </mo>}{d}_2}\right)$ denotes the cumulative distribution function of a multivariate normal function with mean $m_{d_1\mathrm{<mo>\; |\; </mo>}{d}_2}$ and a covariance matrix ${\Sigma }_{d_1\mathrm{<mo>\; |\; </mo>}{d}_2}$.

We now focus on two cases. The first one, related to a simultaneous equations model with two outcomes, one binary and the other one continuous, is treated in Section 2.1. Then, we focus on a case with three outcomes, composed of either one binary and two continuous variables (that is developed in Section 2.2), or two binary and one continuous variables (analyzed in Section 2.3).

2.1. Case with Two Outcomes: One Binary and One Continuous

Let us consider the two following equations:

${\tilde{y}}_{it}^1=Z_{it}^1{\beta }^1+{ϵ}_{it}^1$ (8)

${\tilde{y}}_{it}^2=Z_{it}^2{\beta }^2+{ϵ}_{it}^2$ (9)

in which $y_{it}^1$ is a binary variable equal to 1 if ${\tilde{y}}_{it}^1>0$, and ${\tilde{y}}_{it}^2=y_{it}^2$ is a linear outcome. The associated variance/covariance matrices of error components are:

${\Sigma }_{\nu }=\left(\begin{array}{cc}1& {\rho }_1\sigma \\ {\rho }_1\sigma & {\sigma }^2\end{array}\right)$

${\Sigma }_{\mu }=\left(\begin{array}{cc}{\sigma }_1^2& {\rho }_{1,2}{\sigma }_1{\sigma }_2\\ {\rho }_{1,2}{\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right)$

By identification, we have ${\Sigma }_1=1$, ${\Sigma }_2={\sigma }^2$, and ${\Sigma }_3={\rho }_1\sigma$. Thus, $m_{{\nu }_1|{\nu }_2}=\frac{{\rho }_1}{\sigma }{\nu }_2$ and ${\Sigma }_{{\nu }_1|{\nu }_2}=1-{\rho }_1^2$. In this way, the likelihood has the following form:

$\begin{array}{c}{L}_i={\displaystyle {\int }_{{ℝ}^2}}\left\{{\displaystyle \underset{t=1}{\overset{T_i}{\prod }}}{f}_{\nu }\left({\nu }_{it}^1,{\nu }_{it}^2|\mu \right)\right\}{f}_{\mu }\left({\mu }_i^1,{\mu }_i^2\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\\ ={\displaystyle {\int }_{{ℝ}^2}}\left\{{\displaystyle \underset{t}{\prod }}{\mathcal{l}}_{it}\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}{\mu }_i^2\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\end{array}$

in which ${\mathcal{l}}_{it}$ is the individual likelihood:

${\mathcal{l}}_{it}={\varphi }_1\left({\nu }_{it}^2\mathrm{,0,}{\sigma }^2\right){\Phi }_1\left(\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_1}{\sigma }{\nu }_{it}^2\right)}{\sqrt{1-{\rho }_1^2}}\right)$ (10)

2.2. Case with Three Outcomes: One Binary and Two Continuous

Let us consider the two previous equations (Equation (8) and Equation (9)), and a new one:

${\tilde{y}}_{it}^3=Z_{it}^3{\beta }^3+{ϵ}_{it}^3$ (11)

also corresponding to a linear equation with ${\tilde{y}}_{it}^3=y_{it}^3$. The associated variance/covariance matrices of error components are:

${\Sigma }_{\nu }=\left(\begin{array}{ccc}1& {\rho }_1{\sigma }_a& {\rho }_2{\sigma }_b\\ {\rho }_1{\sigma }_a& {\sigma }_a^2& {\rho }_3{\sigma }_a{\sigma }_b\\ {\rho }_2{\sigma }_b& {\rho }_3{\sigma }_a{\sigma }_b& {\sigma }_b^2\end{array}\right)$

${\Sigma }_{\mu }=\left(\begin{array}{ccc}{\sigma }_1^2& {\rho }_{1,2}{\sigma }_1{\sigma }_2& {\rho }_{1,3}{\sigma }_1{\sigma }_3\\ {\rho }_{1,2}{\sigma }_1{\sigma }_2& {\sigma }_2^2& {\rho }_{2,3}{\sigma }_2{\sigma }_3\\ {\rho }_{1,3}{\sigma }_1{\sigma }_3& {\rho }_{2,3}{\sigma }_2{\sigma }_3& {\sigma }_3^2\end{array}\right)$

By identification, we have ${\Sigma }_1=1$, ${\Sigma }_2=\left(\begin{array}{cc}{\sigma }_a^2& {\rho }_3{\sigma }_a{\sigma }_b\\ {\rho }_3{\sigma }_a{\sigma }_b& {\sigma }_b^2\end{array}\right)$, and ${\Sigma }_3={\left({\rho }_1{\sigma }_a,{\rho }_2{\sigma }_b\right)}^{\prime }$. Thus, $m_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}=\frac{\left({\rho }_1-{\rho }_2{\rho }_3\right)\frac{{\nu }_2}{{\sigma }_a}+\left({\rho }_2-{\rho }_1{\rho }_3\right)\frac{{\nu }_3}{{\sigma }_b}}{1-{\rho }_3^2}$ and ${\Sigma }_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}=\frac{{\rho }_1^2+{\rho }_2^2-2{\rho }_1{\rho }_2{\rho }_3}{1-{\rho }_3^2}$. As a result, we can write the following likelihood:

$\begin{array}{c}{L}_i={\displaystyle {\int }_{{ℝ}^3}}\left\{{\displaystyle \underset{t=1}{\overset{T_i}{\prod }}}{f}_{\nu }\left({\nu }_{it}^1\mathrm{,}{\nu }_{it}^2\mathrm{,}{\nu }_{it}^3\mathrm{<mo>\; |\; </mo>}\mu \right)\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}{\mu }_i^2\mathrm{,,}{\mu }_i^3\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\text{d}{\mu }_i^3\\ ={\displaystyle {\int }_{{ℝ}^3}}\left\{{\displaystyle \underset{t}{\prod }}{\mathcal{l}}_{it}\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}{\mu }_i^2\mathrm{,,}{\mu }_i^3\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\text{d}{\mu }_i^3\end{array}$

in which ${\mathcal{l}}_{it}$ is the individual likelihood, defined as:

${\mathcal{l}}_{it}={\varphi }_2\left(\left({\nu }_{it}^2\mathrm{,}{\nu }_{it}^3\right)\mathrm{,0,}{\Sigma }_2\right){\Phi }_1\left(\frac{q_{it}^1\left(z_{it}^1+m_{{\nu }_1\mathrm{<mo>\; |\; </mo>}\left({\nu }_2\mathrm{,}{\nu }_3\right)}\right)}{\sqrt{{\Sigma }_{{\nu }_1\mathrm{<mo>\; |\; </mo>}\left({\nu }_2\mathrm{,}{\nu }_3\right)}}}\right)$ (12)

2.3. Case with Three Outcomes: Two Binary and One Continuous

In order to derive the likelihood function for a case composed of two binary and one continuous outcomes, let us consider Equation (8) for the first binary outcome and Equation (11) for the linear variable. To consider another binary outcome, let us redefine the following:

${\tilde{y}}_{it}^2=Z_{it}^2{\beta }^2+{ϵ}_{it}^2$ (13)

where $y_{it}^2$ is a binary variable, equal to 1 if ${\tilde{y}}_{it}^2>0$. The associated variance/covariance matrices of error components are:

${\Sigma }_{\nu }=\left(\begin{array}{ccc}1& {\rho }_1& {\rho }_2\sigma \\ {\rho }_1& 1& {\rho }_3\sigma \\ {\rho }_2\sigma & {\rho }_3\sigma & {\sigma }^2\end{array}\right)$

${\Sigma }_{\mu }=\left(\begin{array}{ccc}{\sigma }_1^2& {\rho }_{1,2}{\sigma }_1{\sigma }_2& {\rho }_{1,3}{\sigma }_1{\sigma }_3\\ {\rho }_{1,2}{\sigma }_1{\sigma }_2& {\sigma }_2^2& {\rho }_{2,3}{\sigma }_2{\sigma }_3\\ {\rho }_{1,3}{\sigma }_1{\sigma }_3& {\rho }_{2,3}{\sigma }_2{\sigma }_3& {\sigma }_3^2\end{array}\right)$

By identification, we have ${\Sigma }_1=\left(\begin{array}{cc}1& {\rho }_1\\ {\rho }_1& 1\end{array}\right)$, ${\Sigma }_2={\sigma }^2$, and ${\Sigma }_3={\left({\rho }_2\sigma \mathrm{,}{\rho }_3\sigma \right)}^{\prime }$. Thus, $m_{\left({\nu }_1,{\nu }_2\right)|{\nu }_3}={\left(\frac{{\rho }_2}{\sigma }{\nu }_3\mathrm{,}\frac{{\rho }_3}{\sigma }{\nu }_3\right)}^{\prime }$ and ${\Sigma }_{\left({\nu }_1,{\nu }_2\right)|{\nu }_3}=\left(\begin{array}{cc}1-{\rho }_2^2& {\rho }_1-{\rho }_2{\rho }_3\\ {\rho }_1-{\rho }_2{\rho }_3& 1-{\rho }_3^2\end{array}\right)$. The likelihood has the following form:

$\begin{array}{c}{L}_i={\displaystyle {\int }_{{ℝ}^3}}\left\{{\displaystyle \underset{t=1}{\overset{T_i}{\prod }}}{f}_{\nu }\left({\nu }_{it}^1\mathrm{,}{\nu }_{it}^2\mathrm{,}{\nu }_{it}^3\mathrm{<mo>\; |\; </mo>}\mu \right)\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}{\mu }_i^2\mathrm{,,}{\mu }_i^3\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\text{d}{\mu }_i^3\\ ={\displaystyle {\int }_{{ℝ}^3}}\left\{{\displaystyle \underset{t}{\prod }}{\mathcal{l}}_{it}\right\}{f}_{\mu }\left({\mu }_i^1\mathrm{,}{\mu }_i^2\mathrm{,,}{\mu }_i^3\right)\text{d}{\mu }_i^1\text{d}{\mu }_i^2\text{d}{\mu }_i^3\end{array}$

in which ${\mathcal{l}}_{it}$ is the individual likelihood:

${\mathcal{l}}_{it}={\varphi }_1\left({\nu }_{it}^3\mathrm{,0,}{\sigma }^2\right){\Phi }_2\left(\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_2}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_2^2}}\mathrm{,}\frac{q_{it}^2\left(z_{it}^2-\frac{{\rho }_3}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_3^2}}\mathrm{<mo>\; ;\; </mo>}{q}_{it}^1{q}_{it}^2\rho \right)$ (14)

with $\rho =\frac{{\rho }_1-{\rho }_2{\rho }_3}{\sqrt{\left(1-{\rho }_2^2\right)\left(1-{\rho }_3^2\right)}}$.

3. Estimation Requirement

The likelihood function being a d-dimensional integral function, we use the Gauss-Hermite quadrature method (see Moussa and Delattre [11]). Implementing this method requires to 1) compute the mode $\hat{\mu }$ of the log-integrand $\log \left(f\right)=\log \left(\left\{{\displaystyle {\prod }_{t=1}^{T_i}}{\mathcal{l}}_{it}\right\}{f}_{\mu }\left({\mu }_i^1,\cdots ,{\mu }_i^d\right)\right)$ in $\mu =\left({\mu }_i^1,\cdots ,{\mu }_i^d\right)$ and to derive the Hessian matrix H at $\hat{\mu }$ with respect to $\mu$ ; and, 2) to derive the gradient of the overall likelihood function with respect to the parameters.

Let Q denotes the selected number of quadrature points, x denotes the Q dimension vector of quadrature nodes, and w denotes the Q dimension vector of quadrature weight. By applying the adaptive Gauss-Hermite quadrature [9], the likelihood function in Equation (3) can be rewritten as:

$L_i={\displaystyle \underset{k_1=1}{\overset{Q}{\sum }}}\cdots {\displaystyle \underset{k_d=1}{\overset{Q}{\sum }}}{w}_{k_1}^{*}\cdots w_{k_d}^{*}\left\{{\displaystyle \underset{t=1}{\overset{T_i}{\prod }}}{\mathcal{l}}_{it}\right\}{f_{\mu }\left({\mu }_i^1,\cdots ,{\mu }_i^d\right)|}_{\mu =x^{*}}$ (15)

in which $x^{*}=\hat{\mu }+\sqrt{2}{H}^{-1/2}x$ and $w^{*}={\left(w_1^{*}\cdots w_d^{*}\right)}^{\prime }=2^{d/2}\det \left(H^{-1/2}\right)\cdot \text{diag}\left(w^{\prime }\cdot \exp \left(x^2\right)\right)$.

The derivation of the Hessian matrix is explained in Section 3.1 while the gradient of the overall likelihood function is derived in Section 3.2.

3.1. Hessian Matrix at $\hat{\mu }$

Based on the expressions of ${\mathcal{l}}_{it}$ for each case described in Section 2 (Equation (10), Equation (12), and Equation (14)), we first need to write the associated log-integrand $\log \left(f\right)$ corresponding to each three cases. Then, we can focus on

the calculation of the Hessian matrices, where we need to derive $\frac{{\partial }^2}{\partial {\eta }_i^k\partial {\eta }_i^j}\log \left(f\right)$, in which $k,j=1,\cdots ,d$.

Focusing on the first case with two equations, we have the following log-integrand:

$\log \left(f\right)=\log \left(f_{\mu }\left({\mu }_i^1,{\mu }_i^2\right)\right)+{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\varphi }_1\left({\nu }_{it}^2,0,{\sigma }^2\right)\right)+{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\Phi }_1\left(\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_1}{\sigma }{\nu }_{it}^2\right)}{\sqrt{1-{\rho }_1^2}}\right)\right)$

with the notation $b_{it}=\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_1}{\sigma }{\nu }_{it}^2\right)}{\sqrt{1-{\rho }_1^2}}$, we find:

$\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)=-\frac{1}{{\sigma }_1^2\left(1-{\rho }_{1,2}^2\right)}+\frac{1}{1-{\rho }_1^2}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)-{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}$

$\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)=-\frac{1}{{\sigma }_2^2\left(1-{\rho }_{1,2}^2\right)}-\frac{T_i}{{\sigma }^2}+\frac{{\rho }_1^2}{{\sigma }^2\left(1-{\rho }_1^2\right)}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)-{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}$

$\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)=\frac{{\rho }_{1,2}}{{\sigma }_1{\sigma }_2\left(1-{\rho }_{1,2}^2\right)}+\frac{{\rho }_1}{\sigma \left(1-{\rho }_1^2\right)}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}{q}_{it}^1\frac{-b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}$

Thus, the Hessian matrix is given by:

$H=\left(\begin{array}{cc}-\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)\\ -\frac{{\partial }^2}{\partial {\eta }_i^1\delta {\eta }_i^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)\end{array}\right)$ (16)

Then, focusing on the case of three equations with one binary outcome, the log-integrand is given by:

$\begin{array}{c}\log \left(f\right)=\log \left(f_{\mu }\left({\mu }_i^1,{\mu }_i^2,{\mu }_i^3\right)\right)+{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\varphi }_2\left(\left({\nu }_{it}^2,{\nu }_{it}^3\right),0,{\Sigma }_2\right)\right)\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\Phi }_1\left(\frac{q_{it}^1\left(z_{it}^1+m_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}\right)}{\sqrt{{\Sigma }_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}}}\right)\right)\end{array}$

The associated derivatives, assuming $b_{it}=\frac{q_{it}^1\left(z_{it}^1+m_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}\right)}{\sqrt{{\Sigma }_{{\nu }_1|\left({\nu }_2,{\nu }_3\right)}}}$, $a_1=\sqrt{1-{\rho }_{1,3}^2}$, $a_2=\sqrt{1-{\rho }_{2,3}^2}$, and $ra=\frac{{\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}}{a_1{a}_2}$, are given by:

$\begin{array}{l}\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)=-\frac{1}{{\left({\sigma }_1{a}_1\sqrt{1-r{a}^2}\right)}^2}\\ -\frac{1}{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)=-\frac{1}{{\left({\sigma }_2{a}_2\sqrt{1-r{a}^2}\right)}^2}-\frac{T_i}{\left(1-{\rho }_3^2\right){\sigma }_a^2}\\ -\frac{{\left({\rho }_1-{\rho }_2{\rho }_3\right)}^2}{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right){\left({\sigma }_a\left(1-{\rho }_3^2\right)\right)}^2}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\left({\eta }_i^3\right)}^2}\log \left(f\right)=-\frac{1-{\rho }_{1,2}^2}{{\left({\sigma }_3{a}_1{a}_2\sqrt{1-r{a}^2}\right)}^2}-\frac{T_i}{\left(1-{\rho }_3^2\right){\sigma }_b^2}\\ -\frac{{\left({\rho }_2-{\rho }_1{\rho }_3\right)}^2}{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right){\left({\sigma }_b\left(1-{\rho }_3^2\right)\right)}^2}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)=\frac{{\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}}{{\sigma }_1{\sigma }_2{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ -\frac{{\rho }_1-{\rho }_2{\rho }_3}{{\sigma }_a\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)\left(1-{\rho }_3^2\right)}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}{q}_{it}^1\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^3}\log \left(f\right)=\frac{{\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}}{{\sigma }_1{\sigma }_3{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ -\frac{{\rho }_2-{\rho }_1{\rho }_3}{{\sigma }_b\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)\left(1-{\rho }_3^2\right)}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}{q}_{it}^1\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)=\frac{{\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}}{{\sigma }_2{\sigma }_3{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}+\frac{{\rho }_3{T}_i}{\left(1-{\rho }_3^2\right){\sigma }_a{\sigma }_b}\\ -\frac{\left({\rho }_2-{\rho }_3{\rho }_3\right)\left({\rho }_1-{\rho }_2{\rho }_3\right)}{{\sigma }_a{\sigma }_b\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right){\left(1-{\rho }_3^2\right)}^2}{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{b_{it}{\varphi }_1\left(b_{it}\right){\Phi }_1\left(b_{it}\right)+{\left({\varphi }_1\left(b_{it}\right)\right)}^2}{{\left({\Phi }_1\left(b_{it}\right)\right)}^2}\end{array}$

Thus, the Hessian matrix is given by:

$H=\left(\begin{array}{ccc}-\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^3}\log \left(f\right)\\ -\frac{{\partial }^2}{\partial {\eta }_i^1\delta {\eta }_i^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)\\ -\frac{{\partial }^2}{\partial {\eta }_i^1\delta {\eta }_i^3}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\left({\eta }_i^3\right)}^2}\log \left(f\right)\end{array}\right)$ (17)

Finally, for the three equations with two binary outcomes case, the log-integrand can be written as:

$\begin{array}{c}\log \left(f\right)=\log \left(f_{\mu }\left({\mu }_i^1,{\mu }_i^2,{\mu }_i^3\right)\right)+{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\varphi }_1\left({\nu }_{it}^3,0,{\sigma }^2\right)\right)\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\log \left({\Phi }_2\left(\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_2}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_2^2}},\frac{q_{it}^2\left(z_{it}^2+\frac{{\rho }_3}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_3^2}};q_{it}^1{q}_{it}^2\rho \right)\right)\end{array}$

Then, considering $a_1$, $a_2$ and $ra$ previously defined and the following notations:

$b_{it}^1=\frac{q_{it}^1\left(z_{it}^1+\frac{{\rho }_2}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_2^2}}$

$b_{it}^2=\frac{q_{it}^2\left(z_{it}^2+\frac{{\rho }_3}{\sigma }{\nu }_{it}^3\right)}{\sqrt{1-{\rho }_3^2}}$

$r_n=\sqrt{1-{\rho }_2^2}$

$r_m=\sqrt{1-{\rho }_3^2}$

$\rho =\frac{{\rho }_1-{\rho }_2{\rho }_3}{r_n{r}_m}$

$p_{it}^1=q_{it}^1\varphi \left(b_{it}^1\right)\Phi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)$

$p_{it}^2=q_{it}^2\varphi \left(b_{it}^2\right)\Phi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)$

$p_{it}^3=\frac{{\rho }_3}{\sigma \sqrt{1-{\rho }^2}}\varphi \left(b_{it}^2\right)\Phi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)+\frac{{\rho }_2}{\sigma }\varphi \left(b_{it}^1\right)\Phi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)$

We find the following:

$\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)=-\frac{a_2^2}{{\left({\sigma }_1{a}_1{a}_2\sqrt{1-r_a^2}\right)}^2}+{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{q_{it}^1}{r_n}\left(-b_{it}^1{p}_{it}^1-q_{it}^2\rho \varphi \left(b_{it}^1\right)\varphi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)-{\left(p_{it}^1\right)}^2}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)=-\frac{a_1^2}{{\left({\sigma }_2{a}_1{a}_2\sqrt{1-r_a^2}\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{q_{it}^2}{r_m}\left(-b_{it}^2{p}_{it}^2-q_{it}^1\rho \varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)-{\left(p_{it}^2\right)}^2}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\left({\eta }_i^3\right)}^2}\log \left(f\right)=-\frac{1-{\rho }_1^2}{{\left({\sigma }_3{a}_1{a}_2\sqrt{1-r_a^2}\right)}^2}-\frac{T_i}{{\sigma }^2}-{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{{\left(p_{it}^3\right)}^2}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{{\rho }_3}{\sigma a_2\sqrt{1-{\rho }^2}}\left(\left(\frac{{\rho }_2}{\sigma a_1}-q_{it}^1{q}_{it}^2\rho \frac{{\rho }_3}{\sigma a_2}\right)\varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)-\frac{{\rho }_3}{\sigma }{b}_{it}^2{p}_{it}^2\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{{\rho }_2}{\sigma a_1\sqrt{1-{\rho }^2}}\left(\left(\frac{{\rho }_3}{\sigma a_2}-q_{it}^1{q}_{it}^2\rho \frac{{\rho }_2}{\sigma a_1}\right)\varphi \left(b_{it}^1\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)-\frac{{\rho }_2}{\sigma }{b}_{it}^1{p}_{it}^1\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)=\frac{{\rho }_1-{\rho }_2{\rho }_3}{{\sigma }_1{\sigma }_2{\left(a_1{a}_2\sqrt{1-r_a^2}\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{q_{it}^1{q}_{it}^2}{r_n{r}_m\sqrt{1-{\rho }^2}}\left(\varphi \left(b_{it}^1\right)\varphi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)-p_{it}^1{p}_{it}^2}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^3}\log \left(f\right)=\frac{{\rho }_2-{\rho }_1{\rho }_3}{{\sigma }_1{\sigma }_3{\left(a_1{a}_2\sqrt{1-r_a^2}\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{q_{it}^1}{r_n}\left(\frac{\frac{{\rho }_3}{\sigma r_m}-q_{it}^1{q}_{it}^2\rho \frac{{\rho }_2}{\sigma r_n}}{\sqrt{1-{\rho }^2}}\varphi \left(b_{it}^1\right)\varphi \left(\frac{b_{it}^2-q_{it}^1{q}_{it}^2\rho b_{it}^1}{\sqrt{1-{\rho }^2}}\right)-\frac{{\rho }_2}{\sigma }{b}_{it}^1{p}_{it}^1\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)-p_{it}^1{p}_{it}^3}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\end{array}$

$\begin{array}{l}\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)=\frac{{\rho }_3-{\rho }_1{\rho }_2}{{\sigma }_2{\sigma }_3{\left(a_1{a}_2\sqrt{1-r_a^2}\right)}^2}\\ +{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\frac{q_{it}^2}{r_m}\left(\frac{\frac{{\rho }_2}{\sigma r_n}-q_{it}^1{q}_{it}^2\rho \frac{{\rho }_3}{\sigma r_m}}{\sqrt{1-{\rho }^2}}\varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)-\frac{{\rho }_3}{\sigma }{b}_{it}^2{p}_{it}^2\right){\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)-p_{it}^2{p}_{it}^3}{{\left({\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)\right)}^2}\end{array}$

Thus, the Hessian matrix has the same form as the previous one (matrix 17). In other words, it is given by:

$H=\left(\begin{array}{ccc}-\frac{{\partial }^2}{\partial {\left({\eta }_i^1\right)}^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^1\partial {\eta }_i^3}\log \left(f\right)\\ -\frac{{\partial }^2}{\partial {\eta }_i^1\delta {\eta }_i^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\left({\eta }_i^2\right)}^2}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)\\ -\frac{{\partial }^2}{\partial {\eta }_i^1\delta {\eta }_i^3}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\eta }_i^2\partial {\eta }_i^3}\log \left(f\right)& -\frac{{\partial }^2}{\partial {\left({\eta }_i^3\right)}^2}\log \left(f\right)\end{array}\right)$ (18)

3.2. Gradient Vector with Respect to the Parameters

Based on the likelihood function given by Equation (15), parameters to estimate are ${\beta }^k$, with $k=1,\cdots ,d$, and the associated covariance matrices ${\Sigma }_{\mu }$ and ${\Sigma }_{\nu }$. Thus, the gradient has to be calculated with respect to these parameters. The first order derivative of the log-likelihood function with respect to a parameter $\alpha$, in the set of parameters, is given by:

$\frac{\partial \log \left(L_i\right)}{\partial \alpha }={\displaystyle \underset{k_1=1}{\overset{Q}{\sum }}}\cdots {\displaystyle \underset{k_d=1}{\overset{Q}{\sum }}}\frac{\partial f/\partial \alpha }{L_i}$ (19)

Focusing on the three cases, we apply this formula to compute derivatives with respect to each parameter.

First, considering the two outcomes case, we need to consider the six following parameters: ${\beta }^1$, ${\beta }^2$, $\sigma$, ${\rho }_1$, ${\sigma }_1$, ${\sigma }_2$, ${\rho }_{12}$. As in subsection 3.1, we consider the previously defined $b_{it}$, which is specific to the case with two outcomes, such that we have:

$\frac{\partial f}{\partial {\beta }^1}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{q_{it}^1{\varphi }_1\left(b_{it}\right)}{\sqrt{1-{\rho }_1^2}{\Phi }_1\left(b_{it}\right)}$

$\frac{\partial f}{\partial {\beta }^2}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{{\nu }_{it}^2}{{\sigma }^2}-\frac{{\rho }_1{q}_{it}^1{\varphi }_1\left(b_{it}\right)}{\sigma \sqrt{1-{\rho }_1^2}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log \left(\sigma \right)}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(-1+{\left(\frac{{\nu }_{it}^2}{\sigma }\right)}^2-\frac{{\rho }_1{q}_{it}{\nu }_{it}^2{\varphi }_1\left(b_{it}\right)}{\sigma \sqrt{1-{\rho }_1^2}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_1}{1-{\rho }_1}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left({\rho }_1{q}_{it}{z}_{it}^1+\frac{q_{it}{\nu }_{it}^2}{\sigma }\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_1\right)}=f\ast \left(-1+\frac{{\left(\frac{{\mu }_i^1}{{\sigma }_1}\right)}^2-{\rho }_{1,2}\frac{{\mu }_i^1{\mu }_i^2}{{\sigma }_1{\sigma }_2}}{1-{\rho }_{1,2}^2}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_2\right)}=f\ast \left(-1+\frac{{\left(\frac{{\mu }_i^2}{{\sigma }_2}\right)}^2-{\rho }_{1,2}\frac{{\mu }_i^1{\mu }_i^2}{{\sigma }_1{\sigma }_2}}{1-{\rho }_{1,2}^2}\right)$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{12}}{1-{\rho }_{12}}\right)}^{1/2}}=f{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{\left(q_{it}^1{\rho }_1{z}_{it}^1-\frac{{\nu }_{it}^2}{\sigma }\right){\varphi }_1\left(b_{it}\right)}{\sqrt{1-{\rho }_1^2}{\Phi }_1\left(b_{it}\right)}\\ =f\ast \left({\rho }_{1,2}-\frac{{\rho }_{1,2}{\left(\frac{{\mu }_i^1}{{\sigma }_1}\right)}^2+{\rho }_{1,2}{\left(\frac{{\mu }_i^2}{{\sigma }_2}\right)}^2-\left(1+{\rho }_{1,2}^2\right)\frac{{\mu }_i^1{\mu }_i^2}{{\sigma }_1{\sigma }_2}}{1-{\rho }_{1,2}^2}\right)\end{array}$

Now, focusing on the case of three equations with one binary outcome, we consider notations associated to this case in subsection 3.1 ( $b_{it}$, $a_2$, and $ra$ ). We compute derivatives with respect to ${\beta }^1$, ${\beta }^2$, ${\beta }^3$, ${\sigma }_a$, ${\sigma }_b$, ${\rho }_1$, ${\rho }_2$, ${\rho }_3$, ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$, ${\rho }_{\mathrm{1,2}}$, ${\rho }_{\mathrm{1,3}}$ and ${\rho }_{\mathrm{2,3}}$, such that:

$\frac{\partial f}{\partial {\beta }^1}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{q_{it}^1{\varphi }_1\left(b_{it}\right)}{\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}$

$\frac{\partial f}{\partial {\beta }^2}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{\frac{{\nu }_{it}^2}{{\sigma }_a^2}-{\rho }_3\frac{{\nu }_{it}^3}{{\sigma }_a{\sigma }_b}}{1-{\rho }_3^2}-\frac{q_{it}^1\left({\rho }_1-{\rho }_2{\rho }_3\right){\varphi }_1\left(b_{it}\right)}{{\sigma }_a\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial {\beta }^3}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{\frac{{\nu }_{it}^3}{{\sigma }_b^2}-{\rho }_3\frac{{\nu }_{it}^2}{{\sigma }_a{\sigma }_b}}{1-{\rho }_3^2}-\frac{q_{it}^1\left({\rho }_2-{\rho }_1{\rho }_3\right){\varphi }_1\left(b_{it}\right)}{{\sigma }_b\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_a\right)}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(-1+\frac{{\left(\frac{{\nu }_{it}^2}{{\sigma }_a}\right)}^2-{\rho }_3\frac{{\nu }_{it}^2{\nu }_{it}^3}{{\sigma }_a{\sigma }_b}}{1-{\rho }_3^2}-\frac{q_{it}^1\left({\rho }_1-{\rho }_2{\rho }_3\right){\nu }_{it}^2{\varphi }_1\left(b_{it}\right)}{{\sigma }_a\left(1-{\rho }_3^2\right)\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_b\right)}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(-1+\frac{{\left(\frac{{\nu }_{it}^3}{{\sigma }_b}\right)}^2-{\rho }_3\frac{{\nu }_{it}^2{\nu }_{it}^3}{{\sigma }_a{\sigma }_b}}{1-{\rho }_3^2}-\frac{q_{it}^1\left({\rho }_2-{\rho }_1{\rho }_3\right){\nu }_{it}^3{\varphi }_1\left(b_{it}\right)}{{\sigma }_b\left(1-{\rho }_3^2\right)\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_1}{1-{\rho }_1}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\left(1-{\rho }_1^2\right)\frac{\left(-\frac{q_{it}^1\left({\rho }_3\frac{{\nu }_{it}^3}{{\sigma }_b}-\frac{{\nu }_{it}^2}{{\sigma }_a}\right)}{1-{\rho }_3^2}+b_{it}\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}\frac{{\rho }^2}{\left(1-{\rho }^2\right)\left({\rho }_1-{\rho }_2{\rho }_3\right)}\right){\varphi }_1\left(b_{it}\right)}{\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_2}{1-{\rho }_2}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{\left(\frac{1-{\rho }_2^2}{1-{\rho }_3^2}{q}_{it}^1\left(-{\rho }_3\frac{{\nu }_{it}^2}{{\sigma }_a}+\frac{{\nu }_{it}^3}{{\sigma }_b}\right)+b_{it}\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}\left({\rho }_2+\frac{\left({\rho }_1{\rho }_2-{\rho }_3\right){\rho }^2}{\left(1-{\rho }^2\right)\left({\rho }_1-{\rho }_2{\rho }_3\right)}\right)\right){\varphi }_1\left(b_{it}\right)}{\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\right)$

$\begin{array}{l}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_3}{1-{\rho }_3}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left({\rho }_3\left(1-\frac{{\left(\frac{{\nu }_{it}^2}{{\sigma }_a}\right)}^2+{\left(\frac{{\nu }_{it}^3}{{\sigma }_b}\right)}^2}{1-{\rho }_3^2}\right)+\left(1+{\rho }_3^2\right)\frac{\frac{{\nu }_{it}^2{\nu }_{it}^3}{{\sigma }_a{\sigma }_b}}{1-{\rho }_3^2}\right)\\ +f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{q_{it}^1\left(-\frac{{\rho }_2{\nu }_{it}^2}{{\sigma }_a}-\frac{{\rho }_1{\nu }_{it}^3}{{\sigma }_b}+2{\rho }_3\frac{\left({\rho }_1-{\rho }_2{\rho }_3\right)\frac{{\nu }_{it}^2}{{\sigma }_a}+\left({\rho }_2-{\rho }_1{\rho }_3\right)\frac{{\nu }_{it}^3}{{\sigma }_b}}{1-{\rho }_3^2}+\frac{b_{it}\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}\left({\rho }_1{\rho }_3-{\rho }_2\right){\rho }^2}{\left(1-{\rho }^2\right)\left({\rho }_1-{\rho }_2{\rho }_3\right)}\right){\varphi }_1\left(b_{it}\right)}{\sqrt{\left(1-{\rho }^2\right)\left(1-{\rho }_2^2\right)}{\Phi }_1\left(b_{it}\right)}\end{array}$

$\frac{\partial f}{\partial \log \left({\sigma }_1\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2-\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}-\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_2\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2-\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}-\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_3\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2-\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}-\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{1,2}}{1-{\rho }_{1,2}}\right)}^{1/2}}=f\ast \left(\frac{1-{\rho }_{1,2}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left({\rho }_{1,2}{\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2+\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}-{\rho }_{2,3}\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}-{\rho }_{1,3}\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}\right)\right)\\ +f\frac{\left(1-{\rho }_{1,2}^2\right)\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(1-\frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ +f\frac{\left(1-{\rho }_{1,2}^2\right)\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(\frac{-\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2+2\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\end{array}$

$+f\frac{\left(1-{\rho }_{1,2}^2\right)\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(2\frac{\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{1,3}}{1-{\rho }_{1,3}}\right)}^{1/2}}=f\ast \left(\frac{\left(1-{\rho }_{1,3}^2\right)\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)+\left({\rho }_{1,3}{\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2-{\rho }_{2,3}\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}-{\rho }_{1,2}\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}\right)a_1^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ -f\left({\rho }_{1,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{1,3}-{\rho }_{2,3}\right)\right)\\ \ast \frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2+\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ +2f\left({\rho }_{1,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{1,3}-{\rho }_{2,3}\right)\right)\\ \ast \frac{\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\end{array}$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{2,3}}{1-{\rho }_{2,3}}\right)}^{1/2}}=f\ast \left(\frac{\left(1-{\rho }_{2,3}^2\right)\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)+\left({\rho }_{2,3}{\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2-{\rho }_{1,3}\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}-{\rho }_{1,2}\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}\right)a_2^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ -f\left({\rho }_{2,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{2,3}-{\rho }_{1,3}\right)\right)\\ \ast \frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2+\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ +2f\left({\rho }_{2,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{2,3}-{\rho }_{1,3}\right)\right)\\ \ast \frac{\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\end{array}$

Finally, we compute derivatives with respect to ${\beta }^1$, ${\beta }^2$, ${\beta }^3$, $\sigma$, ${\rho }_1$, ${\rho }_2$, ${\rho }_3$, ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$, ${\rho }_{\mathrm{1,2}}$, ${\rho }_{\mathrm{1,3}}$ and ${\rho }_{\mathrm{2,3}}$ for the three equations with two binary outcomes case. To do so, we consider notations defined for this case in subsection 3.1 concerning $b_{it}^1$, $b_{it}^2$, $a_1$, $a_2$, $ra$, $r_n$, $r_m$, $\rho$, $p_{it}^1$, and $p_{it}^2$, such that:

$\frac{\partial f}{\partial {\beta }^1}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{p_{it}^1}{{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}$

$\frac{\partial f}{\partial {\beta }^2}=f*{\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{p_{it}^2}{{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}$

$\frac{\partial f}{\partial {\beta }^3}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{{\nu }_{it}^3}{{\sigma }^2}-\frac{{\rho }_3{p}_{it}^2}{\sigma q_{it}^2{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}-\frac{{\rho }_2{p}_{it}^1}{\sigma q_{it}^1{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}\right)$

$\frac{\partial f}{\partial \log \left(\sigma \right)}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(-1+{\left(\frac{{\nu }_{it}^3}{\sigma }\right)}^2-\frac{{\rho }_3{\nu }_{it}^3{p}_{it}^2}{\sigma q_{it}^2{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}-\frac{{\rho }_2{\nu }_{it}^3{p}_{it}^1}{\sigma q_{it}^1{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_1}{1-{\rho }_1}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\frac{q_{it}^1{q}_{it}^2\left(1-{\rho }_1^2\right)\varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)}{r_n{r}_m\sqrt{1-{\rho }^2}{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_2}{1-{\rho }_2}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{q_{it}^1{q}_{it}^2\left({\rho }_1{\rho }_2-{\rho }_3\right)\varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)}{r_n{r}_m\sqrt{1-{\rho }^2}{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}+\frac{\left({\rho }_2{b}_{it}^1-\frac{r_n{\nu }_{it}^3}{\sigma }\right)p_{it}^1}{q_{it}^1{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_3}{1-{\rho }_3}\right)}^{1/2}}=f\ast {\displaystyle \underset{t=1}{\overset{T_i}{\sum }}}\left(\frac{q_{it}^1{q}_{it}^2\left({\rho }_1{\rho }_3-{\rho }_2\right)\varphi \left(b_{it}^2\right)\varphi \left(\frac{b_{it}^1-q_{it}^1{q}_{it}^2\rho b_{it}^2}{\sqrt{1-{\rho }^2}}\right)}{r_n{r}_m\sqrt{1-{\rho }^2}{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}+\frac{\left({\rho }_3{b}_{it}^2-\frac{r_m{\nu }_{it}^3}{\sigma }\right)p_{it}^2}{q_{it}^2{\Phi }_2\left(b_{it}^1,b_{it}^2;q_{it}^1{q}_{it}^2\rho \right)}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_1\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_i^1}{{\sigma }_1}\right)}^2-\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_i^1{\mu }_i^2}{{\sigma }_1{\sigma }_2}-\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_i^1{\mu }_i^3}{{\sigma }_1{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_2\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_i^2}{{\sigma }_2}\right)}^2-\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_i^1{\mu }_i^2}{{\sigma }_1{\sigma }_2}-\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_i^2{\mu }_i^3}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\frac{\partial f}{\partial \log \left({\sigma }_3\right)}=f\ast \left(-1+\frac{\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_i^3}{{\sigma }_3}\right)}^2-\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_i^1{\mu }_i^3}{{\sigma }_1{\sigma }_3}-\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_i^2{\mu }_i^3}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)$

$\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{1,2}}{1-{\rho }_{1,2}}\right)}^{1/2}}=f\ast \left(\frac{1-{\rho }_{1,2}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left({\rho }_{1,2}{\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2+\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}-{\rho }_{2,3}\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}-{\rho }_{1,3}\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}\right)\right)$

$\begin{array}{c}+f\frac{\left(1-{\rho }_{\mathrm{1,2}}^2\right)\left({\rho }_{\mathrm{1,2}}-{\rho }_{\mathrm{1,3}}{\rho }_{\mathrm{2,3}}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(1-\frac{\left(1-{\rho }_{\mathrm{2,3}}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{\mathrm{1,3}}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ +f\frac{\left(1-{\rho }_{\mathrm{1,2}}^2\right)\left({\rho }_{\mathrm{1,2}}-{\rho }_{\mathrm{1,3}}{\rho }_{\mathrm{2,3}}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(\frac{-\left(1-{\rho }_{\mathrm{1,2}}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2+2\left({\rho }_{\mathrm{1,2}}-{\rho }_{\mathrm{1,3}}{\rho }_{\mathrm{2,3}}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ +f\frac{\left(1-{\rho }_{\mathrm{1,2}}^2\right)\left({\rho }_{\mathrm{1,2}}-{\rho }_{\mathrm{1,3}}{\rho }_{\mathrm{2,3}}\right)}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\left(2\frac{\left({\rho }_{\mathrm{1,3}}-{\rho }_{\mathrm{1,2}}{\rho }_{\mathrm{2,3}}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{\mathrm{2,3}}-{\rho }_{\mathrm{1,2}}{\rho }_{\mathrm{1,3}}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\end{array}$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{1,3}}{1-{\rho }_{1,3}}\right)}^{1/2}}=f\ast \left(\frac{\left(1-{\rho }_{1,3}^2\right)\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)+\left({\rho }_{1,3}{\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2-{\rho }_{2,3}\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}-{\rho }_{1,2}\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}\right)a_1^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ -f\left({\rho }_{1,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{1,3}-{\rho }_{2,3}\right)\right)\\ \ast \frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2+\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ +2f\left({\rho }_{1,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{1,3}-{\rho }_{2,3}\right)\right)\\ \ast \frac{\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\end{array}$

$\begin{array}{c}\frac{\partial f}{\partial \log {\left(\frac{1+{\rho }_{2,3}}{1-{\rho }_{2,3}}\right)}^{1/2}}=f\ast \left(\frac{\left(1-{\rho }_{2,3}^2\right)\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)+\left({\rho }_{2,3}{\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2-{\rho }_{1,3}\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}-{\rho }_{1,2}\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}\right)a_2^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\right)\\ -f\left({\rho }_{2,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{2,3}-{\rho }_{1,3}\right)\right)\\ \ast \frac{\left(1-{\rho }_{2,3}^2\right){\left(\frac{{\mu }_1^i}{{\sigma }_1}\right)}^2+\left(1-{\rho }_{1,3}^2\right){\left(\frac{{\mu }_2^i}{{\sigma }_2}\right)}^2+\left(1-{\rho }_{1,2}^2\right){\left(\frac{{\mu }_3^i}{{\sigma }_3}\right)}^2}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}\\ +2f\left({\rho }_{2,3}{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2+\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\left({\rho }_{1,2}{\rho }_{2,3}-{\rho }_{1,3}\right)\right)\end{array}$

$\ast \frac{\left({\rho }_{1,2}-{\rho }_{1,3}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_2^i}{{\sigma }_1{\sigma }_2}+\left({\rho }_{1,3}-{\rho }_{1,2}{\rho }_{2,3}\right)\frac{{\mu }_1^i{\mu }_3^i}{{\sigma }_1{\sigma }_3}+\left({\rho }_{2,3}-{\rho }_{1,2}{\rho }_{1,3}\right)\frac{{\mu }_2^i{\mu }_3^i}{{\sigma }_2{\sigma }_3}}{{\left(a_1{a}_2\sqrt{1-r{a}^2}\right)}^2}$

4. Examples and Comparisons with Roodman’s Command

In order to shed light on the advantages and the consistency of our method (xtcmp), we decide to implement examples using a dataset, previously used for xtsur¹, in Stata software. This database is an unbalanced panel database of 1672 observations, corresponding to 142 individuals followed between 1990 and 2003. All explanatory variables used are quantitative and do not contain any missing values.

We implement two cases: 1) a system of two equations with one linear and one binary dependent variables; 2) a system of three equations with two binary and one continuous outcomes. Indeed, let us consider the three following equations:

${\tilde{y}}_{it}^1={\beta }_0+x{1}_{it}{\beta }_1+x{2}_{it}{\beta }_2+x{3}_{it}{\beta }_3+x{4}_{it}{\beta }_4+{\mu }_i^1+{\nu }_{it}^1$ (20)

${\tilde{y}}_{it}^2={\gamma }_0+x{4}_{it}{\gamma }_1+x{6}_{it}{\gamma }_2+x{7}_{it}{\gamma }_3+{\mu }_i^2+{\nu }_{it}^2$ (21)

${\tilde{y}}_{it}^3={\alpha }_0+x{7}_{it}{\alpha }_1+x{9}_{it}{\alpha }_2+{\mu }_i^3+{\nu }_{it}^3$ (22)

where $y_{it}^1$ and $y_{it}^2$ are binary variables equal to 1 if ${\tilde{y}}_{it}^1>0$ and ${\tilde{y}}_{it}^2>0$ ; and where ${\tilde{y}}_{it}^3=y_{it}^3$ is a linear outcome.

When considering the first case (a system of two equations) with one continuous and one binary outcomes, we focus on the simultaneous estimation of Equation (20) and Equation (22). In this case, the associated variance/covariance matrices of the error components are given by:

${\Sigma }_{\nu }=\left(\begin{array}{cc}1& \rho \sigma \\ \rho \sigma & {\sigma }^2\end{array}\right)\text{and}{\Sigma }_{\mu }=\left(\begin{array}{cc}{\sigma }_1^2& {\rho }_1{\sigma }_1{\sigma }_2\\ {\rho }_1{\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right)$

We consider four estimation techniques: 1) each equation is estimated separately as a single panel equation; 2) the two equations are estimated with cmp, while considering a pooled equation; 3) the two equations are estimated with cmp, with a posterior estimate of random effects; and, 4) the two equations are estimated with our method (xtcmp) presented in Section 2. Results are presented in Table 1.

Results suggest that first, estimating equations separately is misleading since the covariance between idiosyncratic errors and individual random effects is not considered; and, the significance of the coefficients appears to be false. On the other hand, considering equations simultaneously aims at obtaining consistent estimates among the three tests performed (coefficients are closer to each others in the last three columns). However, one should notice that, when using cmp with random effects, the estimation of the covariance matrices for both individual random effects and idiosyncratic errors is done after the estimation of the

Table 1. Two equations system with both a linear and a non-linear outcome.

***: significant at the 1% level; **: significant at the 5% level; *: significant at the 10% level.

coefficients (post-estimation). Comparing with our results (xtcmp, last column), we can see that the variance of the individual effects seems to be overestimated in cmp’s case.

Then, we offer an example for the second case, a three equations system with two binary outcomes and one continuous dependent variable, such that we consider Equations (20)-(22). In this case, the associated variance/covariance matrices of the error components are given by:

${\Sigma }_{\nu }=\left(\begin{array}{ccc}1& {\rho }_1& {\rho }_2\sigma \\ {\rho }_1& 1& {\rho }_3\sigma \\ {\rho }_2\sigma & {\rho }_3\sigma & {\sigma }^2\end{array}\right)\text{and}{\Sigma }_{\mu }=\left(\begin{array}{ccc}{\sigma }_1^2& {\rho }_{\mathrm{1,2}}{\sigma }_1{\sigma }_2& {\rho }_{\mathrm{1,3}}{\sigma }_1{\sigma }_3\\ {\rho }_{\mathrm{1,2}}{\sigma }_1{\sigma }_2& {\sigma }_2^2& {\rho }_{\mathrm{2,3}}{\sigma }_2{\sigma }_3\\ {\rho }_{\mathrm{1,3}}{\sigma }_1{\sigma }_3& {\rho }_{\mathrm{2,3}}{\sigma }_2{\sigma }_3& {\sigma }_3^2\end{array}\right)$

For this example, we cannot provide the third estimation technique where we used cmp command with a posterior estimate of the random effects, because this procedure does not converge. However, we provide the other estimation techniques: 1) each equation is estimated separately as a single panel equation; 2) the three equations are estimated with cmp, considering a pooled database; and, 3) the three equations are estimated with our method (xtcmp). Results are provided in Table 2.

Table 2. Three equations system with two binary outcomes.

***: significant at the 1% level; **: significant at the 5% level; *: significant at the 10% level.

Results suggest that, as before, estimating equations separately leads to errors in the significance of coefficients, especially for the first outcome. Moreover, such method does not consider the covariance between idiosyncratic errors and individual random effects. On the other hand, considering equations simultaneously allows to obtain more consistent estimates (coefficient estimates seem closer in the last two columns, and significance of the latter is persistent along the two last columns).

5. Conclusions

xtcmp is a command implemented under Stata software. We focus on three main cases: 1) a simultaneous equations model with two equations (including one linear and one binary outcomes); 2) a case with three equations composed of two linear and one binary outcomes; and, 3) a three equations case with one linear and two binary dependent variables. This command further develops Roodman’s [1] command cmp which does not explicitly consider the panel dimension of the data, nor simultaneous equations model since it is written as a SUR estimator. This technical note gives detailed description of the computations, namely likelihood functions, log-integrand associated, Hessian matrices and gradient vectors with respect to each parameter of interest, specific to the three cases described above.

xtcmp’s estimation framework could be further developed in order to consider a broader range of non-linear outcomes (such as ordered probit, multinomial probit or truncated framework, for instance) or to consider much more equations simultaneously, in a dynamic setup.

Still, as it stands, xtcmp represents a significant development in Stata’s commands. Indeed, it allows researchers to resolve endogeneity issues in a panel dimension context by analyzing correlation in the error terms of the equations and thus specific individual effects depending on the outcomes.

Acknowledgements

We thank E. Delattre for his helpful comments.

NOTES

¹See Nguyen [12] for more detailed information.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Roodman, D. (2011) Fitting Fully Observed Recursive Mixed-Process Models with cmp. The Stata Journal, 11, 159-206. https://doi.org/10.1177/1536867X1101100202
[2]	Adeline, A. and Delattre, E. (2018) Health and Income: Testing for Causality on European Elderly People. THEMA Working Paper No. 2018-07.
[3]	Williams Greene, H. (2012) Econometric Analysis. Seventh Edition, Prentice Hall, Upper Saddle River, NJ.
[4]	Jeffrey Wooldridge, M. (2010) Econometric Analysis of Cross Section and Panel Data. Second Edition, MIT Press, Cambridge, MA.
[5]	Deb, P. and Trivedi, P.K. (2006) Maximum Simulated Likelihood Estimation of a Negative Binomial Regression Model with Multinomial Endogenous Treatment. The Stata Journal, 6, 246-255. https://doi.org/10.1177/1536867X0600600206
[6]	Kenkel, D.S. and Terza, J.V. (2001) The Effect of Physician Advice on Alcohol Consumption: Count Regression with an Endogenous Treatment Effect. Journal of Applied Econometrics, 16, 165-184. https://doi.org/10.1002/jae.596
[7]	Bratti, M. and Miranda, A. (2011) Endogenous Treatment Effects for Count Data Models with Endogenous Participation or Sample Selection. Health Economics, 20, 1090-1109. https://doi.org/10.1002/hec.1764
[8]	Esterling, K.M., Neblo, M.A. and Lazer, D.M.J. (2011) Estimating Treatment Effects in the Presence of Noncompliance and Nonresponse: The Generalized Endogenous Treatment Model. Political Analysis, 19, 205-226. https://doi.org/10.1093/pan/mpr005
[9]	Liu, Q. and Pierce, D. (1994) A Note on Gauss-Hermite Quadrature. Biometrika, 83, 624-629. https://doi.org/10.2307/2337136
[10]	Gould, W., Pitblado, J. and Poi, B. (2010) Maximum Likelihood Estimation with Stata. Fourth Edition, Stata Press, College Station, Texas.
[11]	Moussa, R. and Delattre, E. (2018) On the Estimation of Causality in a Bivariate Dynamic Probit Model on Panel Data with Stata Software. A Technical Review. Theoretical Economics Letter, 8, 1257-1278. https://doi.org/10.4236/tel.[phone]
[12]	Nguyen, M.C. (2008) XTSUR: Stata Module to Estimate Seemingly Unrelated Regression Model on Unbalanced Panel Data. Statistical Software Components S456953, Boston College Department of Economics.