On the Estimation of Causality in a Bivariate Dynamic Probit Model on Panel Data with Stata Software: A Technical Review

In order to assess causality between binary economic outcomes, we consider the estimation of a bivariate dynamic probit model on panel data that has the particularity to account the initial conditions of the dynamic process. Due to the intractable form of the likelihood function that is a two dimensions integral, we use an approximation method: The adaptative Gauss-Hermite quadrature method. For the accuracy of the method and to reduce computing time, we derive the gradient of the log-likelihood and the Hessian of the integrand. The estimation method has been implemented using the d1 method of Stata software. We made an empirical validation of our estimation method by applying on simulated data set. We also analyze the impact of the number of quadrature points on the estimations and on the estimation process duration. We then conclude that when exceeding 16 quadrature points on our simulated data set, the relative differences in the estimated coefficients are around 0.01% but the computing time grows up exponentially.

KEYWORDS

1. Introduction

Testing Granger causality has generated a large set of paper in the literature. The larger part of this literature concerns the case where we have continuous dependent variables. For binary outcomes, there is also a way to consider the causality problem. As described by [1] for a vector of dependant variables, the one order Granger causality can be analyse as a probability conditional independence given a set of exogenous variables and the first order lagged dependent variables. And for a binary outcome in the dependent vector, one can use a probit probability that implies the use of latent variable.

For panel data case, as the one way fix effects model estimated on a finite sample has necessarily inconsistent estimators [2] , the random effect model is used. Due to the fact that we aim to test for one order Granger causality, lagged dependent variables are included as explanatory variables. For the first wave of the panel, we do not have previous values for the dependent variables, and treating them casually or as exogenous leads to inconsistent estimators [2] . So we specify an other equation for initial conditions as described by [3] . The equation is allowed to have different explanatory variables and different idiosyncratic error terms from the dynamic equation.

This specification leads to a likelihood function with an intractable form that is a two dimensions integral with a large set of parameters to be estimated. The estimation of this likelihood function requires the use of numerical approximation of integral function such as maximum simulated likelihood (see [4] for more details) or Gauss-Hermite quadrature (for more details see [5] [6] [7]).

The main goal of this paper is to propose and to test a method for estimating a two equations system where the explanatory variables are binary in a panel data framework. To the extent of our knowledge, there is no program to do so, especially as we propose the calculation of the Hessian matrix and the gradient vector of our maximisation program.

In this paper, we discuss on the problem of testing Granger causality with a bivariate dynamic probit model taking into account the initial conditions. The organization of this paper is the following one. In Section 2 we explain the causality test method for bivariate probit model with panel data. In Section 3, we describe the estimation method available when the likelihood function has an intractable form (two dimensions integral in our case). Section 4 presents the calculation of the gradient with respect to the model parameters and the calculation of the Hessian matrix with respect to the random effects vector. In Section 5, we present a robustness analysis of our selected estimation method by doing some simulations1.

2. Testing Causality with a Bivariate Dynamic Probit Model

This section aims to describe causality test method in the case of binary variables. We start by presenting the general approach in time series before introducing panel data case. We end this section by a discussion on the initial conditions problem.

2.1. Testing Causality: General Approach

Causality concept was introduced by [8] as a better predictability of a variable Y by the use of it lag values, the lag value of an other variable Z and some controls X. In his paper, [8] distinguishes instantaneous causality that means Zt is causing Yt (if Zt is included in the model it improves the predictability of Yt than if not) from lag causality that means lag values of Z improve the predictability of Yt. In this section, we rule out the instantaneous causality and deal with lag causality of one period.

The one period Granger causality can be rephrase in terms of conditional independence. Without lost of generality, we present the univariate case for time series. Let’s Yt and Zt denote some dependent variables and Xt denote a set of controls variables. One period Granger non-causality from Z to Y is the conditional independence of Yt from ${Z}_{t-1}$ conditionally to Xt and ${Y}_{t-1}$ . More clearly, Granger non-causality from Z to Y is:

$f\left({Y}_{t}|{Y}_{t-1},{X}_{t},{Z}_{t-1}\right)=f\left({Y}_{t}|{Y}_{t-1},{X}_{t}\right)$ (1)

Note that the same kind of relationship can be written for Granger non-causality from Y to Z. As Yt and Zt are binary outcome variables, we can use latent variables ( ${Y}^{*}$ and ${Z}^{*}$ respectively) and make the assumption that Y and Z have positive outcomes (equals to 1) if their latent variables are positive. The latent variables are defined as follows:

For the left side of the Equation (1) ( $f\left({Y}_{t}|{Y}_{t-1},{X}_{t},{Z}_{t-1}\right)$ ):

${Y}_{t}^{*}={X}_{t}{\beta }_{1}+{\delta }_{11}{Y}_{t-1}+{\delta }_{12}{Z}_{t-1}+{ϵ}_{t}^{1}$ (2)

${Z}_{t}^{*}={X}_{t}{\beta }_{2}+{\delta }_{21}{Y}_{t-1}+{\delta }_{22}{Z}_{t-1}+{ϵ}_{t}^{2}$ (3)

For the right side of the Equation (1) ( $f\left({Y}_{t}|{Y}_{t-1},{X}_{t}\right)$ ):

${Y}_{t}^{*}={X}_{t}{\beta }_{1}+{\delta }_{11}{Y}_{t-1}+{ϵ}_{t}^{1}$ (4)

${Z}_{t}^{*}={X}_{t}{\beta }_{2}+{\delta }_{21}{Z}_{t-1}+{ϵ}_{t}^{2}$ (5)

where

$\left(\begin{array}{c}{ϵ}^{1}\\ {ϵ}^{2}\end{array}\right)⇝N\left(0,{\Sigma }_{ϵ}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}{\Sigma }_{ϵ}=\left(\begin{array}{cc}1& {\rho }_{ϵ}\\ {\rho }_{ϵ}& 1\end{array}\right)$

To fit the joint distribution of Y and Z conditionally to X (meaning that we estimate a bivariate model), we need to analyze four available situations that are $\left(Y=Z=1\right)$ , $\left(Y=Z=0\right)$ , $\left(Y=1;Z=0\right)$ and $\left(Y=0;Z=1\right)$ . For each of these situations, we have:

$\begin{array}{l}P\left({Y}_{t}=1,{Z}_{t}=1|{X}_{t}\right)\\ =P\left({ϵ}_{t}^{1}>-{X}_{t}{\beta }_{1}-{\delta }_{11}{Y}_{t-1}-{\delta }_{12}{Z}_{t-1},{ϵ}_{t}^{2}>-{X}_{t}{\beta }_{2}-{\delta }_{21}{Y}_{t-1}-{\delta }_{22}{Z}_{t-1}\right)\end{array}$

$\begin{array}{l}P\left({Y}_{t}=0,{Z}_{t}=0|{X}_{t}\right)\\ =P\left({ϵ}_{t}^{1}<-{X}_{t}{\beta }_{1}-{\delta }_{11}{Y}_{t-1}-{\delta }_{12}{Z}_{t-1},{ϵ}_{t}^{2}<-{X}_{t}{\beta }_{2}-{\delta }_{21}{Y}_{t-1}-{\delta }_{22}{Z}_{t-1}\right)\end{array}$

$\begin{array}{l}P\left({Y}_{t}=1,{Z}_{t}=0|{X}_{t}\right)\\ =P\left({ϵ}_{t}^{1}>-{X}_{t}{\beta }_{1}-{\delta }_{11}{Y}_{t-1}-{\delta }_{12}{Z}_{t-1},{ϵ}_{t}^{2}<-{X}_{t}{\beta }_{2}-{\delta }_{21}{Y}_{t-1}-{\delta }_{22}{Z}_{t-1}\right)\end{array}$

$\begin{array}{l}P\left({Y}_{t}=0,{Z}_{t}=1|{X}_{t}\right)\\ =P\left({ϵ}_{t}^{1}<-{X}_{t}{\beta }_{1}-{\delta }_{11}{Y}_{t-1}-{\delta }_{12}{Z}_{t-1},{ϵ}_{t}^{2}>-{X}_{t}{\beta }_{2}-{\delta }_{21}{Y}_{t-1}-{\delta }_{22}{Z}_{t-1}\right)\end{array}$

As we can see, by assuming ${q}_{t}^{1}=2{Y}_{t}-1$ and ${q}_{t}^{2}=2{Z}_{t}-1$ , we can rewrite the probabilities above as:

$\begin{array}{l}P\left({Y}_{t},{Z}_{t}|{X}_{t}\right)\\ ={\Phi }_{2}\left({q}_{t}^{1}\left({X}_{t}{\beta }_{1}+{\delta }_{11}{Y}_{t-1}+{\delta }_{12}{Z}_{t-1}\right),{q}_{t}^{2}\left({X}_{t}{\beta }_{2}+{\delta }_{21}{Y}_{t-1}+{\delta }_{22}{Z}_{t-1}\right),{q}_{t}^{1}{q}_{t}^{2}{\rho }_{ϵ}\right)\end{array}$ (6)

where ${\Phi }_{2}\left(\text{ }\right)$ stands for the bivariate normal c.d.f.

Then testing Granger non-causality in this specification is testing $H0:{\delta }_{12}=0$ for Z is not causing Y and testing $H0:{\delta }_{21}=0$ for Y is not causing Z.

2.2. Testing Causality: Panel Data Case

For panel data case, two major approaches can be used. The first one is to consider that causal effect is not the same for all individuals in the panel ( [9] ). This approach is useful when individuals are heterogeneous or when the causal effect is not homogeneous. The specification for latent variables are:

${Y}_{it}^{*}={X}_{t}^{1}{\beta }_{1}+{\delta }_{11,i}{Y}_{i,t-1}+{\delta }_{12,i}{Z}_{i,t-1}+{\eta }_{i}^{1}+{\zeta }_{it}^{1}$ (7)

${Z}_{it}^{*}={X}_{t}^{2}{\beta }_{2}+{\delta }_{21,i}{Y}_{i,t-1}+{\delta }_{22,i}{Z}_{i,t-1}+{\eta }_{i}^{2}+{\zeta }_{it}^{2}$ (8)

where ${\left({\eta }_{i}^{1},{\eta }_{i}^{2}\right)}^{\prime }$ denotes the individual random effects which are zero mean and covariance matrix ${\Sigma }_{\eta }$ and ${\left({\zeta }_{it}^{1},{\zeta }_{it}^{2}\right)}^{\prime }$ denote the idiosyncratic shocks which are zero mean and covariance matrix ${\Sigma }_{\zeta }$ with

${\Sigma }_{\eta }=\left(\begin{array}{cc}{\sigma }_{1}^{2}& {\sigma }_{1}{\sigma }_{2}{\rho }_{\eta }\\ {\sigma }_{1}{\sigma }_{2}{\rho }_{\eta }& {\sigma }_{2}^{2}\end{array}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }{\Sigma }_{\zeta }=\left(\begin{array}{cc}1& {\rho }_{\zeta }\\ {\rho }_{\zeta }& 1\end{array}\right)$

In this approach, testing Granger non-causality is equivalent to test ${\delta }_{12,i}=0,i=1,\cdots ,N$ for Z is not causing Y and to test ${\delta }_{21,i}=0,i=1,\cdots ,N$ for Y is not causing Z.

The second approach (that is used in this paper) is to assume the causal effects, if they exist, are the same for all individuals in the panel. With the same notation that the previous case, the latent variables are:

${Y}_{it}^{*}={X}_{t}^{1}{\beta }_{1}+{\delta }_{11}{Y}_{i,t-1}+{\delta }_{12}{Z}_{i,t-1}+{\eta }_{i}^{1}+{\zeta }_{it}^{1}$ (9)

${Z}_{it}^{*}={X}_{t}^{2}{\beta }_{2}+{\delta }_{21}{Y}_{i,t-1}+{\delta }_{22}{Z}_{i,t-1}+{\eta }_{i}^{2}+{\zeta }_{it}^{2}$ (10)

Then testing Granger non-causality is equivalent to test $H0:{\delta }_{12}=0$ for Z is not causing Y and to test $H0:{\delta }_{21}=0$ for Y is not causing Z.

Finally, Equations (9) and (10) are the core of our problem. Since Y and Z are binary panel outcomes and each equation includes lag dependent variables, estimating jointly these two equations can be viewed as the estimation of a bivariate dynamic probit model.

2.3. Dealing with Initial Conditions

For the first wave of the panel (initial conditions), due to the fact that we do not have data for the previous state on Y and Z (no values for ${Y}_{i,0}$ and ${Z}_{i,0}$ ) we are not able to evaluate $P\left({Y}_{i1},{Z}_{i1}|{Y}_{i,0},{Z}_{i,0},{X}_{i}\right)$ . By ignoring it in the individual likelihood, researchers also ignore the data generation process for the first wage of the panel. This means that they assume the data generating process of the first wave of the panel to be exogenous or to be in equilibrium. These assumptions hold only if the individual random effects are degenerated. If this assumption is not fulfilled, the initial conditions (the first wave of the panel) are explained by the individual random effects and ignoring them leads to inconsistent parameter estimates [2] .

The solution proposed by [2] for the univariate case and generalized by [3] is to estimate a static equation for the first wave of the panel (meaning that we do not introduce lagged dependent variables). In this static equation, the random effects are a linear combination of the random effects in the next wave of the panel and idiosyncratic error terms may have different structure from the idiosyncratic error terms in the dynamic equation. Formally, the latent variables for the first wave of the panel are defined as follows:

${Y}_{i,1}^{*}={X}_{i}^{1}{\gamma }_{1}+{\lambda }_{11}{\eta }_{i}^{1}+{\lambda }_{12}{\eta }_{i}^{2}+{ϵ}_{i}^{1}$ (11)

${Z}_{i,1}^{*}={X}_{i}^{2}{\gamma }_{2}+{\lambda }_{21}{\eta }_{i}^{1}+{\lambda }_{22}{\eta }_{i}^{2}+{ϵ}_{i}^{2}$ (12)

where ${\left({ϵ}_{i}^{1},{ϵ}_{i}^{2}\right)}^{\prime }$ denotes the vector of idiosyncratic shocks which are zero mean and covariance matrix ${\Sigma }_{ϵ}$ with ${\Sigma }_{ϵ}=\left(\begin{array}{cc}1& {\rho }_{ϵ}\\ {\rho }_{ϵ}& 1\end{array}\right)$ .

As ${\eta }^{1}$ and ${\eta }^{2}$ are individual random effects respectively on Y and Z, ${\lambda }_{12}$ and ${\lambda }_{21}$ can be interpreted as the influence of the Y random individual effects (respectively Z random individual effects) on Z (respectively on Y) at the first wave of the panel.

3. Estimation Methods

Due to the fact that the likelihood function has an intractable form (an integral function), it is impossible to estimate this likelihood by usual methods. We then deal with numerical integration methods that are numerical approximation method for an integral. In this section we describe two major methods and argue for one of them to estimate our likelihood function.

The Gauss-Hermite quadrature is a numerical approximation method use to close the value of an integral function. The default approach is related to an univariate integral of the form:

${\int }_{ℝ}\text{ }\text{ }f\left(x\right)\mathrm{exp}\left(-{x}^{2}\right)\text{d}x$ (13)

where $\mathrm{exp}\left(-{x}^{2}\right)$ denotes the Gaussian factor2. Then the integral above can be approximated using:

${\int }_{ℝ}\text{ }\text{ }f\left(x\right)\mathrm{exp}\left(-{x}^{2}\right)\text{d}x=\underset{q=1}{\overset{Q}{\sum }}\text{ }\text{ }{w}_{q}\ast f\left({x}_{q}\right)$ (14)

where ${x}_{q},q=1,\cdots ,Q$ are nodes from the Hermite polynomial and ${w}_{q},q=1,\cdots ,Q$ are corresponding weights.

This approximation supposes that the integrand can be well approximated by an $2Q+1$ order polynomial and that the integrand is sampled on a symmetric range centered on zero. So, for suitable results, these two assumptions must be taken into account.

We first assume that finding the optimal number of quadrature points can be achieved numerically. For the accuracy of the approximation, it is required to choose the optimal number of quadrature points. To do this, one can start with a number $\stackrel{¯}{q}$ of quadrature points and increase it to assess if it significantly changes the result, and repeat this process until convergence in terms of overall likelihood value variation and estimated coefficients variation. But, it is also important to take into account the fact that increasing number of quadrature points also increases the computing time. An example of the impact of number of quadrature points on estimated results is given in Section 5.

For the problem of suitable sampling range, the solution of using the adaptative Gauss-Hermite quadrature was proposed by [5] and [6] . In this approach, instead of using $\mathrm{exp}\left(-{x}^{2}\right)$ as a Gaussian factor, they use a Gaussian density $\varphi \left(\mu ,{\sigma }^{2}\right)$ of mean μ and variance σ2. That means (see [5] ):

${\int }_{ℝ}\text{ }f\left(x\right)\text{d}x=\underset{q=1}{\overset{Q}{\sum }}\text{ }\text{ }{w}_{q}^{\ast }f\left({x}_{q}^{*}\right)$ (15)

Then the sampling range is transformed and the new nodes are ${x}_{q}^{*}=\mu +\sqrt{2}\sigma {x}_{q}$ and weights are ${w}_{q}^{*}=\sqrt{2}\sigma {w}_{q}\mathrm{exp}\left({x}_{q}^{2}\right)$ . For [5] , one can choose the normal density with posterior mean and variance equal respectively to μ and σ2. For the implementation, we can start with $\mu =0$ and $\sigma =1$ and at each iteration of the likelihood maximization process, calculate the posterior weighted mean and variance of the quadrature points and use them to calculate the nodes and weights for the next iteration. For [6] , one can choose μ to be the mode of the integrand $f\left(x\right)$ and σ to be the square of the Hessian of the log of integrand taken in the mode.

$\sigma ={\left(-{\frac{{\partial }^{2}}{\partial {x}^{2}}\mathrm{log}\left(f\left(x\right)\right)|}_{x=\stackrel{^}{x}}\right)}^{-1/2}$ (16)

For the multivariate integral case, the same approach is used. Without lost of generality, we discuss the bivariate case that can be apply to others multivariate cases. The function to approximate is written as follows:

${\int }_{{ℝ}^{2}}\text{ }f\left(x,y\right)\text{d}x\text{d}y$ (17)

With the assumption of independence between x and y (that can be overcome by using a Cholesky decomposition ${x}^{\prime }=x$ and ${y}^{\prime }=\rho {x}^{\prime }+y$ , see [5] or [7] for more precision on these Cholesky transformation or other transformations that can lead to similar results) the integral above can be approximated by:

${\int }_{{ℝ}^{2}}\text{ }f\left(x,y\right)\text{d}x\text{d}y=\underset{{q}_{1}=1,{q}_{2}=1}{\overset{Q}{\sum }}{w}_{{q}_{1}}^{*}{w}_{{q}_{2}}^{*}f\left({x}_{{q}_{1}}^{*},{y}_{{q}_{1}}^{*}\right)$ (18)

And in this case, the nodes and weights are derived as follows:

$\left(\begin{array}{c}{x}_{{q}_{1}}^{*}\\ {y}_{{q}_{1}}^{*}\end{array}\right)=\stackrel{^}{x}+\sqrt{2}\ast {\left(-{\frac{{\partial }^{2}}{\partial {x}^{2}}\mathrm{log}\left(f\left(x,y\right)\right)|}_{x,y=\stackrel{^}{x}}\right)}^{-1/2}\ast \left(\begin{array}{c}{x}_{{q}_{1}}\\ {y}_{{q}_{1}}\end{array}\right)$ (19)

and

$\left(\begin{array}{c}{w}_{{q}_{1}}^{*}\\ {w}_{{q}_{2}}^{*}\end{array}\right)=2\ast {|-{\frac{{\partial }^{2}}{\partial {x}^{2}}\mathrm{log}\left(f\left(x,y\right)\right)|}_{\left(x,y\right)=\stackrel{^}{x}}|}^{-1/2}\ast \left(\begin{array}{c}{w}_{{q}_{1}}\mathrm{exp}\left({x}_{{q}_{1}}^{2}\right)\\ {w}_{{q}_{2}}\mathrm{exp}\left({x}_{{q}_{2}}^{2}\right)\end{array}\right)$ (20)

where $|A|$ denotes the determinant of matrix A, and $\stackrel{^}{x}$ denotes the mode of the integrand $f\left(x,y\right)$ .

Jackel (2005) also suggests that for the nodes with low weights (when contributions to the integral value are not significant) we can prune the range from those nodes in order to save calculation time. That means to set a scalar

$\tau =\frac{{w}_{1}{w}_{|\left(Q+1\right)/2|}}{Q}$ and drop all nodes with weights lower than this scalar.

3.2. Maximum Simulated Likelihood Method

Maximum Simulated Likelihood method was introduced by [4] as a solution to maximization problems that have an integral as objective function. In this approach, the likelihood function is supposed to be defined as:

$f\left(x,y\right)={\int }_{{ℝ}^{2}}\text{ }{f}^{*}\left(x,y,{u}_{1},{u}_{2}\right)g\left({u}_{1},{u}_{2}\right)\text{d}{u}_{1}\text{d}{u}_{2}$ (21)

where $g\left({u}_{1},{u}_{2}\right)$ is a probability distribution function, ${f}^{*}\left(x,y,{u}_{1},{u}_{2}\right)$ is called simulator and denotes the function from which the mean value at some draws u1 and u2 gives an approximation of the overall likelihood. Without lost of generality, we only define the two dimensions case that can be generalized to fewer or larger dimensions integral. For this kind of likelihood function, [4] proposed as simulator the function ${f}^{*}\left(x,y,{u}_{1},{u}_{2}\right)$ with u1 and u2 drawn from the same probability distribution function g (the probability distribution function of the individual random effects). Then the overall likelihood function can be approximated by (u1d denotes the dth draw from u1; the same definition holds for u2d):

$f\left(x,y\right)=\frac{1}{D}\underset{d=1}{\overset{D}{\sum }}\text{ }\text{ }{f}^{*}\left(x,y,{u}_{1d},{u}_{2d}\right)$ (22)

where D denotes the number of draws.

To implement this method, we start by simulating a bivariate normal draw $N\left(0,{I}_{2}\right)$ and we give them the $\left({u}_{1},{u}_{2}\right)$ covariance matrix structure. Then we calculate the value of the simulator at these transformed draws and we repeat D times. The overall likelihood is the mean of the simulator value at each transformed draw. At each iteration, once the random effects covariance matrix is calculated, we apply it to the simulated first normal draws to transform them in draws of the random effects and use them to calculate the likelihood. This process is repeated until convergence.

The simulated likelihood estimator is consistent and asymptotically equivalent to the likelihood estimator ( [4] ) if the number of draws tend to infinity faster than $\sqrt{N}$ .

3.3. GHQ or MSL: What Method to Choose?

As described above, they are two main methods to estimate our likelihood function. To choose which method to implement, we deal with the accuracy and the computing time requirement.

For our estimations, we choose the adaptative Gauss-Hermite quadrature proposed by [7] for three main reasons.

• Our dataset is an unbalanced panel data with 10,569 individuals observed in mean over 26 years, that leads 255,206 observations. Due to the fact that the simulated likelihood method requires that the number of draw D be larger than the square of the number of observations, we do not use it to avoid waste of time in computing process.

• The Gauss-Hermite quadrature requires that we find the best number of quadrature Q that is the one for whom the integrand can be well approximated by an $2Q+1$ order polynomial. If Q is small, that reduces computing time. Our estimations are achieved in general for Q between 8 and 14. It means that at each iteration, for the likelihood value calculation, we do a weighted sum of between ${8}^{2}=64$ and ${14}^{2}=196$ terms.

• Using the Gauss-Hermite quadrature method reduces computing time but this computing time remains very long if the integrand is not sampled at the suitable range (meaning that the adaptative method has not been used). And in this case, the maximization process spends between two and three weeks before achieving convergence on an Intel Core i7 computer at 3.4 GHz with 8 GB of RAM memory. By applying the adaptative Gauss-Hermite quadrature, the computing time is significantly reduced and then, we spend between two and three days for achieving convergence on the same computer.

Note that the reduced convergence time mentioned above is in part due to the implementation of the first order derivatives of the likelihood function. Using the overall log-likelihood approximated by the Liu and Pierce adaptative Gauss-Hermite quadrature method, we can get derivatives with respect of all model parameters. The implementation of these derivative in the maximization process allows us to used the Stata’s d1 method. The convergence time saved by this method is clearly huge. On our overall data set, with 8 quadratures points, when we use a non adaptative quadrature method, the convergence is not achieved: after 3 weeks of computation, the model underflows. When we use the [6] adaptative Gauss-Hermite quadrature, but without implementing the first order derivatives, the estimation process takes 11 days and 10 hours to achieve convergence. When we use the adaptative Gauss-Hermite quadrature in [6] with implemented the first order derivatives, the estimation process achieve convergence only after 1 day and 17 hours.

4. Chosen Method Requirements

In this section we describe some requirements of the selected method that is the adaptative Gauss-Hermite Quadrature. The first one is the fact that the adaptative Gauss-Hermite quadrature requires to derive the Hessian of the log of the integrand ( [6] ). The second one is that we derive the gradient of the overall likelihood function in order to use Stata’s d1 method (see [10] ) for more accuracy and more speed in the calculations.

The gradient of the overall log-likelihood function has been calculated to speed up the maximization process. This will allow us to use the Stata’s d1 method that requires the implementation of the gradient vector in addition to the log-likelihood. The likelihood function for an individual i is:

${L}_{i}={\int }_{{ℝ}^{2}}\text{ }{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\underset{t=2}{\overset{{T}_{i}}{\prod }}\text{\hspace{0.17em}}{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\varphi \left({\eta }_{i},{\Sigma }_{\eta }\right)\text{d}{\eta }_{i}^{1}\text{d}{\eta }_{i}^{2}$ (23)

where

${q}_{it}^{1}=2{y}_{it}^{1}-1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall \text{\hspace{0.17em}}i,t$

${q}_{it}^{2}=2{y}_{it}^{2}-1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall \text{\hspace{0.17em}}i,t$

${h}_{i}^{0}={Z}_{i}^{1}{\gamma }_{1}+{\lambda }_{11}{\eta }_{i}^{1}+{\lambda }_{12}{\eta }_{i}^{2}$

${w}_{i}^{0}={Z}_{i}^{2}{\gamma }_{2}+{\lambda }_{21}{\eta }_{i}^{1}+{\lambda }_{22}{\eta }_{i}^{2}$

${\stackrel{¯}{h}}_{it}={X}_{it}^{1}{\beta }_{1}+{\delta }_{11}{h}_{i,t-1}+{\delta }_{12}{w}_{i,t-1}+{\eta }_{i}^{1}$

${\stackrel{¯}{w}}_{it}={X}_{it}^{2}{\beta }_{2}+{\delta }_{21}{h}_{i,t-1}+{\delta }_{22}{w}_{i,t-1}+{\eta }_{i}^{2}$

Using the adaptative Gauss-Hermite quadrature method by [6] , the overall likelihood function is given by (we use the same notation that those used in Section 3):

$\begin{array}{l}{L}_{i}=\underset{k=1,j=1}{\overset{Q}{\sum }}{w}_{k}^{*}{w}_{j}^{*}{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}×\underset{t=2}{\overset{{T}_{i}}{\prod }}\text{ }{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){\varphi \left({\eta }_{i},{\Sigma }_{\eta }\right)|}_{{\eta }_{i}^{1}={x}_{k}^{*},{\eta }_{i}^{2}={x}_{j}^{*}}\end{array}$ (24)

To get the gradient vector, the log-likelihood above must be derived with respect to 13 parameters that are: ${\stackrel{¯}{\beta }}_{1}={\left({\beta }_{1},{\delta }_{11},{\delta }_{12}\right)}^{\prime }$ , ${\stackrel{¯}{\beta }}_{2}={\left({\beta }_{2},{\delta }_{21},{\delta }_{22}\right)}^{\prime }$ , ${\gamma }_{1}$ , ${\gamma }_{2}$ , ${\lambda }_{11}$ , ${\lambda }_{12}$ , ${\lambda }_{21}$ , ${\lambda }_{22}$ , ${\sigma }_{1}$ , ${\sigma }_{2}$ , ${\rho }_{\eta }$ , ${\rho }_{\zeta }$ , and ${\rho }_{ϵ}$ .

Let’s ${l}_{kj}$ denotes:

${l}_{kj}={\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\underset{t=2}{\overset{{T}_{i}}{\prod }}\text{ }{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){\varphi \left({\eta }_{i},{\Sigma }_{\eta }\right)|}_{{\eta }_{i}^{1}={x}_{k}^{*},{\eta }_{i}^{2}={x}_{j}^{*}}$

Then, the first order derivatives with respect to each α of the 13 parameters is given by:

$\frac{\partial \mathrm{log}\left({L}_{i}\right)}{\partial \alpha }=\underset{k=1,j=1}{\overset{Q}{\sum }}\frac{\partial {l}_{kj}/\partial \alpha }{{L}_{i}}$

With respect to ${\stackrel{¯}{\beta }}_{1}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\stackrel{¯}{\beta }}_{1}}={l}_{kj}\underset{t=2}{\overset{{T}_{i}}{\sum }}\frac{{q}_{it}^{1}\varphi \left({q}_{it}^{1}{\stackrel{¯}{h}}_{it}\right){\Phi }_{1}\left(\frac{{q}_{it}^{2}{\stackrel{¯}{w}}_{it}-{q}_{it}^{2}{\rho }_{\zeta }{\stackrel{¯}{h}}_{it}}{\sqrt{1-{\rho }_{\zeta }^{2}}}\right)}{{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}$

With respect to ${\stackrel{¯}{\beta }}_{2}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\stackrel{¯}{\beta }}_{2}}={l}_{kj}\underset{t=2}{\overset{{T}_{i}}{\sum }}\frac{{q}_{it}^{2}\varphi \left({q}_{it}^{2}{\stackrel{¯}{w}}_{it}\right){\Phi }_{1}\left(\frac{{q}_{it}^{1}{\stackrel{¯}{h}}_{it}-{q}_{it}^{1}{\rho }_{\zeta }{\stackrel{¯}{w}}_{it}}{\sqrt{1-{\rho }_{\zeta }^{2}}}\right)}{{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}$

With respect to ${\gamma }_{1}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\gamma }_{1}}={l}_{kj}\frac{{q}_{i0}^{1}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{q}_{i0}^{2}{\rho }_{ϵ}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\gamma }_{2}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\gamma }_{2}}={l}_{kj}\frac{{q}_{i0}^{2}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{q}_{i0}^{1}{\rho }_{ϵ}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\lambda }_{11}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\lambda }_{11}}={l}_{kj}\frac{{q}_{i0}^{1}{x}_{k}^{*}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{q}_{i0}^{2}{\rho }_{ϵ}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\lambda }_{12}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\lambda }_{12}}={l}_{kj}\frac{{q}_{i0}^{1}{x}_{j}^{*}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{q}_{i0}^{2}{\rho }_{ϵ}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\lambda }_{21}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\lambda }_{21}}={l}_{kj}\frac{{q}_{i0}^{2}{x}_{k}^{*}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{q}_{i0}^{1}{\rho }_{ϵ}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\lambda }_{22}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial {\lambda }_{22}}={l}_{kj}\frac{{q}_{i0}^{2}{x}_{j}^{*}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{q}_{i0}^{1}{\rho }_{ϵ}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

With respect to ${\sigma }_{1}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial \mathrm{log}\left({\sigma }_{1}\right)}={l}_{kj}\ast \left(-1+\frac{{\left({x}_{k}^{*}/{\sigma }_{1}\right)}^{2}-{\rho }_{\eta }{x}_{k}^{*}{x}_{j}^{*}/\left({\sigma }_{1}{\sigma }_{2}\right)}{1-{\rho }_{\eta }^{2}}\right)$

With respect to ${\sigma }_{2}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial \mathrm{log}\left({\sigma }_{2}\right)}={l}_{kj}\ast \left(-1+\frac{{\left({x}_{j}^{*}/{\sigma }_{2}\right)}^{2}-{\rho }_{\eta }{x}_{k}^{*}{x}_{j}^{*}/\left({\sigma }_{1}{\sigma }_{2}\right)}{1-{\rho }_{\eta }^{2}}\right)$

With respect to ${\rho }_{\eta }$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial \mathrm{log}{\left(\frac{1+{\rho }_{\eta }}{1-{\rho }_{\eta }}\right)}^{1/2}}={l}_{kj}\ast \left({\rho }_{\eta }-\frac{{\rho }_{\eta }\left({\left({x}_{k}^{*}/{\sigma }_{1}\right)}^{2}+{\left({x}_{j}^{*}/{\sigma }_{2}\right)}^{2}\right)-\left(1+{\rho }_{\eta }^{2}\right){x}_{k}^{*}{x}_{j}^{*}/\left({\sigma }_{1}{\sigma }_{2}\right)}{1-{\rho }_{\eta }^{2}}\right)$

With respect to ${\rho }_{\zeta }$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial \mathrm{log}{\left(\frac{1+{\rho }_{\zeta }}{1-{\rho }_{\zeta }}\right)}^{1/2}}={l}_{kj}\underset{t=2}{\overset{{T}_{i}}{\sum }}\frac{{q}_{it}^{1}{q}_{it}^{2}\varphi \left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}$

With respect to ${\rho }_{ϵ}$ the first order derivative is:

$\frac{\partial {l}_{kj}}{\partial \mathrm{log}{\left(\frac{1+{\rho }_{ϵ}}{1-{\rho }_{ϵ}}\right)}^{1/2}}={l}_{kj}\frac{{q}_{i0}^{1}{q}_{i0}^{2}\varphi \left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}$

Remarks:

• For ${\sigma }_{1}$ , ${\sigma }_{2}$ , ${\rho }_{\eta }$ , ${\rho }_{\zeta }$ , and ${\rho }_{ϵ}$ , we used some transformations on parameters to insure that in the maximization process, each σ remains positive and each ρ remains between −1 and 1 at all iteration. For σ we use exponential transformation then in the derivation, we derive with respect to

$\mathrm{log}\left(\sigma \right)$ . For ρ we use arc-tangency transformation (i.e. $\frac{\mathrm{exp}\left(2\rho \right)-1}{\mathrm{exp}\left(2\rho \right)+1}$ ) then in the derivation, we derive with respect to $\mathrm{log}{\left(\frac{1+\rho }{1-\rho }\right)}^{1/2}$ .

• To easily derive a bivariate normal probability with zero mean, variance one and correlation ρ, we can transform it into an integral where the integrand is a product of an univariate normal density and an univariate normal probability as follows:

${\Phi }_{2}\left(x,y,\rho \right)={\int }_{-\infty }^{y}\text{ }\text{ }\varphi \left(v\right)\Phi \left(\frac{x-\rho v}{\sqrt{1-{\rho }^{2}}}\right)\text{d}v={\int }_{-\infty }^{x}\text{ }\text{ }\varphi \left(u\right)\Phi \left(\frac{y-\rho u}{\sqrt{1-{\rho }^{2}}}\right)\text{d}u.$

• Given the transformation above, the first order derivatives of ${\Phi }_{2}\left(x,y,\rho \right)$ with respect to x and y are respectively given by:

$\frac{\partial {\Phi }_{2}\left(x,y,\rho \right)}{\partial x}=\varphi \left(x\right)\Phi \left(\frac{y-\rho x}{\sqrt{1-{\rho }^{2}}}\right)$

$\frac{\partial {\Phi }_{2}\left(x,y,\rho \right)}{\partial y}=\varphi \left(y\right)\Phi \left(\frac{x-\rho y}{\sqrt{1-{\rho }^{2}}}\right)$

4.2. Hessian Matrix Calculation

For the requirement of the adaptative Gauss-Hermite quadrature method, we need to derive the Hessian matrix of the log of the integrand function with respect to the random effects vector3. From the individual likelihood function defined in Equation 23, the log of the integrand is:

$\begin{array}{l}\mathrm{log}\left(g\left({\eta }_{i}^{1},{\eta }_{i}^{2}\right)\right)\\ =\mathrm{log}\left({\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\underset{t=2}{\overset{{T}_{i}}{\prod }}\text{\hspace{0.17em}}{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\varphi \left({\eta }_{i},{\Sigma }_{\eta }\right)\right)\end{array}$ (25)

We derive from the log of the integrand in Equation (25) the Hessian matrix by calculating:

$-\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{1}\right)}^{2}}\mathrm{log}\left(g\left({\eta }_{i}^{1},{\eta }_{i}^{2}\right)\right)$

$-\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{2}\right)}^{2}}\mathrm{log}\left(g\left({\eta }_{i}^{1},{\eta }_{i}^{2}\right)\right)$

$-\frac{{\partial }^{2}}{\partial {\eta }_{i}^{1}\partial {\eta }_{i}^{1}}\mathrm{log}\left(g\left({\eta }_{i}^{1},{\eta }_{i}^{2}\right)\right)$

The first order derivatives are given by:

$\begin{array}{c}-\frac{\partial }{\partial {\eta }_{i}}\mathrm{log}\left(g\right)=-\frac{{{\Phi }^{\prime }}_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}-\underset{t=2}{\overset{{T}_{i}}{\sum }}\frac{{{\Phi }^{\prime }}_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{\eta }_{i}^{1}/{\sigma }_{1}^{2}-\rho {\eta }_{i}^{2}/\left({\sigma }_{1}{\sigma }_{2}\right)}{1-{\rho }_{\eta }^{2}}\end{array}$

With respect to ${\eta }_{i}^{1}$ we have:

$\begin{array}{c}{{\Phi }^{\prime }}_{2{\eta }_{i}^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)={q}_{i0}^{1}{\lambda }_{11}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{q}_{i0}^{2}{\rho }_{ϵ}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{q}_{i0}^{2}{\lambda }_{21}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{q}_{i0}^{1}{\rho }_{ϵ}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\end{array}$

${{\Phi }^{\prime }}_{2{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)={q}_{it}^{1}\varphi \left({q}_{it}^{1}{\stackrel{¯}{h}}_{it}\right){\Phi }_{1}\left(\frac{{q}_{it}^{2}{\stackrel{¯}{w}}_{it}-{q}_{it}^{2}{\rho }_{\zeta }{\stackrel{¯}{h}}_{it}}{\sqrt{1-{\rho }_{\zeta }^{2}}}\right)$

And with respect to ${\eta }_{i}^{2}$ we have:

$\begin{array}{c}{{\Phi }^{\prime }}_{2{\eta }_{i}^{2}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)={q}_{i0}^{1}{\lambda }_{12}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{q}_{i0}^{2}{\rho }_{ϵ}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{q}_{i0}^{2}{\lambda }_{22}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{q}_{i0}^{1}{\rho }_{ϵ}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\end{array}$

${{\Phi }^{\prime }}_{2{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)=\varphi \left({q}_{it}^{2}{\stackrel{¯}{w}}_{it}\right){\Phi }_{1}\left(\frac{{q}_{it}^{1}{\stackrel{¯}{h}}_{it}-{q}_{it}^{1}{\rho }_{\zeta }{\stackrel{¯}{w}}_{it}}{\sqrt{1-{\rho }_{\zeta }^{2}}}\right)$

The second order derivatives are given by:

$\begin{array}{c}-\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{1}\right)}^{2}}\mathrm{log}\left(g\right)=-\frac{{{\Phi }^{″}}_{2{\eta }_{i}^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right){\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{{\Phi }^{\prime }}_{2{\eta }_{i}^{1}}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\underset{t=2}{\overset{{T}_{i}}{\sum }}\left(\frac{{{\Phi }^{″}}_{2{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{{{\Phi }^{\prime }}_{2{\eta }_{i}^{1}}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\right)+\frac{1}{{\sigma }_{1}^{2}\left(1-{\rho }_{\eta }^{2}\right)}\end{array}$ (26)

$\begin{array}{c}-\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{2}\right)}^{2}}\mathrm{log}\left(g\right)=-\frac{{{\Phi }^{″}}_{2{\eta }_{2}^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right){\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{{\Phi }^{\prime }}_{2{\eta }_{i}^{2}}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\underset{t=2}{\overset{{T}_{i}}{\sum }}\left(\frac{{{\Phi }^{″}}_{2{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{{{\Phi }^{\prime }}_{2{\eta }_{i}^{2}}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\right)+\frac{1}{{\sigma }_{1}^{2}\left(1-{\rho }_{\eta }^{2}\right)}\end{array}$ (27)

$\begin{array}{c}-\frac{{\partial }^{2}}{\partial {\eta }_{i}^{1}\delta {\eta }_{i}^{2}}\mathrm{log}\left(g\right)=-\frac{{{\Phi }^{″}}_{2{\eta }_{i}^{1}{\eta }_{i}^{2}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{\epsilon }\right){\Phi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{{\Phi }^{\prime }}_{2{\eta }^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right){{\Phi }^{\prime }}_{2{\eta }^{2}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}{{\Phi }_{2}^{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\underset{t=2}{\overset{{T}_{i}}{\sum }}\left(\frac{{{\Phi }^{″}}_{2{\eta }_{i}^{1}{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){\Phi }_{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\end{array}$

$\begin{array}{c}-\frac{{{\Phi }^{\prime }}_{2{\eta }^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right){{\Phi }^{\prime }}_{2{\eta }^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}{{\Phi }_{2}^{2}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)}\right)\\ -\frac{{\rho }_{\eta }}{{\sigma }_{1}{\sigma }_{2}\left(1-{\rho }_{\eta }^{2}\right)}\end{array}$ (28)

where

$\begin{array}{l}{{\Phi }^{″}}_{2{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\\ =-{\stackrel{¯}{h}}_{it}{{\Phi }^{\prime }}_{2{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)-{\rho }_{\zeta }{\varphi }_{{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\end{array}$

$\begin{array}{l}{{\Phi }^{″}}_{2{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\\ =-{\stackrel{¯}{w}}_{it}{{\Phi }^{\prime }}_{2{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)-{\rho }_{\zeta }{\varphi }_{{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)\end{array}$

${{\Phi }^{″}}_{2{\eta }_{i}^{1}{\eta }_{i}^{2}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)={q}_{it}^{1}{q}_{it}^{2}\rho \zeta {\varphi }_{{\eta }_{i}^{1}}\left({q}_{it}^{1}{\stackrel{¯}{h}}_{it},{q}_{it}^{2}{\stackrel{¯}{w}}_{it},{q}_{it}^{1}{q}_{it}^{2}{\rho }_{\zeta }\right)$

$\begin{array}{l}{{\Phi }^{″}}_{2{\eta }_{i}^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ =\left(2{\lambda }_{11}{\lambda }_{21}-{\rho }_{ϵ}\left({\lambda }_{11}^{2}+{\lambda }_{21}^{2}\right)\right)\varphi \left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{11}^{2}{h}_{i}^{0}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{\rho }_{\epsilon }{q}_{i0}^{2}{h}_{i}^{0}}{\sqrt{1-{\rho }_{\epsilon }^{2}}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{21}^{2}{w}_{i}^{0}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{\rho }_{ϵ}{q}_{i0}^{1}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\end{array}$

$\begin{array}{l}{{\Phi }^{″}}_{2{\eta }_{i}^{2}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ =\left(2{\lambda }_{12}{\lambda }_{22}-{\rho }_{ϵ}\left({\lambda }_{12}^{2}+{\lambda }_{22}^{2}\right)\right)\varphi \left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{12}^{2}{h}_{i}^{0}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{\rho }_{ϵ}{q}_{i0}^{2}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{22}^{2}{w}_{i}^{0}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{\rho }_{\epsilon }{q}_{i0}^{1}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\end{array}$

$\begin{array}{l}{{\Phi }^{″}}_{2{\eta }_{i}^{1}{\eta }_{i}^{1}}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ ={q}_{i0}^{1}{q}_{i0}^{2}\left({\lambda }_{11}{\lambda }_{22}+{\lambda }_{12}{\lambda }_{21}-{\rho }_{ϵ}\left({\lambda }_{11}{\lambda }_{12}+{\lambda }_{21}{\lambda }_{22}\right)\right){\varphi }_{2}\left({q}_{i0}^{1}{h}_{i}^{0},{q}_{i0}^{2}{w}_{i}^{0},{q}_{i0}^{1}{q}_{i0}^{2}{\rho }_{ϵ}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{11}{\lambda }_{12}{h}_{i}^{0}\varphi \left({q}_{i0}^{1}{h}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{2}{w}_{i}^{0}-{\rho }_{ϵ}{q}_{i0}^{2}{h}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\lambda }_{21}{\lambda }_{22}{w}_{i}^{0}\varphi \left({q}_{i0}^{2}{w}_{i}^{0}\right){\Phi }_{1}\left(\frac{{q}_{i0}^{1}{h}_{i}^{0}-{\rho }_{ϵ}{q}_{i0}^{1}{w}_{i}^{0}}{\sqrt{1-{\rho }_{ϵ}^{2}}}\right)\end{array}$

Then, the Hessian matrix is given by:

$H=\left(\begin{array}{cc}-\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{1}\right)}^{2}}\mathrm{log}\left(g\right)& -\frac{{\partial }^{2}}{\partial {\eta }_{i}^{1}\partial {\eta }_{i}^{2}}\mathrm{log}\left(g\right)\\ -\frac{{\partial }^{2}}{\partial {\eta }_{i}^{1}\delta {\eta }_{i}^{2}}\mathrm{log}\left(g\right)& -\frac{{\partial }^{2}}{\partial {\left({\eta }_{i}^{2}\right)}^{2}}\mathrm{log}\left(g\right)\end{array}\right)$ (29)

As described in Section 3.1, after having derived this Hessian matrix, we calculate its value at the mode of the integrand and use it to re-sample the integrand.

5. Robustness Analysis Based on Simulations

This section aims to insure that the implemented method gives suitable results. We consider that the implemented method give us suitable results if for a given relationship between variables, by applying the estimation method on these variables we find approximatively the same coefficients. To reach this goal, we perform a robustness analysis on the estimation method. This robustness analysis is an empirical one based on simulations. We use two different approaches for that.

The first approach is to simulate bivariate binary variables by specifying a relationship between some explanatory variables (it means that we set coefficients of explanatory variables) and estimate this relationship with the implemented method in order to compare the results with the relationship specified before. In the second approach, we introduce new variables (that were not used in the data generating process) when estimating the relationship with the implemented method and compare the new results with the first ones. The implemented method is robust when it is able to correctly estimate the relationship specified even if we introduce other variables and also to estimate non significant coefficients to those other variables. Finally, the method we make use of to check for the robustness is the same that in [11] .

As the estimation method implemented is a numerical approximation method, the results will depend on the selected number of quadrature points. We deal with the incidence of number of quadrature points on results in the last part of this section. For a better analysis of the results we also add the standard errors of each estimated coefficients.

5.1. Simulated Relationship between Real Variables

In this section, we use variables from the French SIP (Santé et Itinéraire Professionnel) survey data set and we simulate error terms and a relationship between some selected variables. The subset of the database use for this section is an unbalanced panel of 1202 individuals with total waves per individual between 5 and 10 waves.

We set the error terms parameters as ${\sigma }_{1}=2.1$ , ${\sigma }_{2}=3.1$ , ${\rho }_{\eta }=0.7$ , ${\rho }_{\zeta }=0.5$ and ${\rho }_{ϵ}=0.4$ .

We simulate idiosyncratic errors vectors $\zeta ={\left({\zeta }_{1},{\zeta }_{2}\right)}^{\prime }$ and $ϵ={\left({ϵ}_{1},{ϵ}_{2}\right)}^{\prime }$ as bivariate normal variables with zero mean, variance equal to 1 and covariances respectively equal to ${\rho }_{\zeta }$ and ${\rho }_{ϵ}$ . We also simulate individual random effects as bivariate normal variables with zero mean, covariance equals to ${\rho }_{\eta }$ and variance equals to ${\sigma }_{1}^{2}$ for the first component of the random effects vector and equals to ${\sigma }_{2}^{2}$ for the second component of the random effects vector. It has been done as follows:

${ϵ}_{1}=rnormal\left(0,1\right)$

${ϵ}_{2}=rnormal\left(0,1\right)\ast \sqrt{1-{\rho }_{ϵ}^{2}}+{\rho }_{ϵ}{ϵ}_{1}$

${\zeta }_{1}=rnormal\left(0,1\right)$

${\zeta }_{2}=rnormal\left(0,1\right)\ast \sqrt{1-{\rho }_{\zeta }^{2}}+{\rho }_{\zeta }{\zeta }_{1}$

where $rnormal\left(\mu ,\sigma \right)$ denote the random normal density with mean μ and standard deviation σ. As individuals effects are time invariant, we simulate η as follows:

${\eta }_{1}=rnormal\left(0,{\sigma }_{1}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{ }\text{ }t=1$

${\eta }_{2}=rnormal\left(0,{\sigma }_{2}\right)\ast \sqrt{1-{\rho }_{\eta }^{2}}+{\rho }_{\eta }\frac{{\sigma }_{2}}{{\sigma }_{1}}{\eta }_{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{if}\text{ }\text{\hspace{0.17em}}t=1$

${\eta }_{1}={\eta }_{1}\left[t=1\right]\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{ }\text{\hspace{0.17em}}t\ne 1$

${\eta }_{2}={\eta }_{2}\left[t=1\right]\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{if}\text{\hspace{0.17em}}t\ne 1$

For the initial conditions ( $t=1$ ), the simulated relationship is the following:

${y}_{1}^{*}=-0.2+0.3ill-0.2unemp+0.4{\eta }_{1}-0.5{\eta }_{2}+{ϵ}_{1}$

${y}_{2}^{*}=2-0.2ill-0.08age+0.3{\eta }_{1}+0.5{\eta }_{2}+{ϵ}_{2}$

${y}_{1}=\mathbb{I}\left({y}_{1}^{*}>0\right)$

${y}_{2}=\mathbb{I}\left({y}_{2}^{*}>0\right)$

For $t>1$ , we specify the following relationship:

${y}_{1t}^{*}=1.9+0.3{y}_{1,t-1}+0.1{y}_{2,t-1}-0.05Mal{e}_{t}-0.2unem{p}_{t}+{\eta }_{1}+{\zeta }_{1t}$

${y}_{2t}^{*}=-0.4-0.1{y}_{1,t-1}+0.4{y}_{2,t-1}+0.05Mal{e}_{t}-0.5den{s}_{t}+{\eta }_{2}+{\zeta }_{2t}$

${y}_{1t}=\mathbb{I}\left({y}_{1t}^{*}>0\right)$

${y}_{2t}=\mathbb{I}\left({y}_{2t}^{*}>0\right)$

The variable ill denotes having an illness episode in the year, unemp denotes being out of labour marking during the year, age denotes the age of individual, and Male is 1 if individual is male and 0 otherwise. Estimation results for 16 quadrature points are displayed in Table 1. For all equations, we give the coefficients that are used in the DGP and those that are estimated by our program. As we can see, all the coefficients from the DGP are very closed from the estimates ones.

5.2. Simulated Relationship with Additional Variables

In this section, we keep the same DGP than in Section 5.1 and we add other variables in the model that we estimate in order to evaluate the robustness of the estimation method by the fact that all estimated coefficients for variables in the DGP should remain the same and the added variables coefficients should not

Table 1. Simulated data set estimation’s results.

Estimated standard deviations for estimated coefficients are given within parenthesis. ***: significant at the 1% level, **: significant at the 5% level.

significant. We introduce two variables rural and nationality (not French) in the dynamic equations of the regression.

Results are in Table 2. Columns 1 and 2 in Table 2 are the same than corresponding columns in Table 1. We provide in Table 2, column 3, the new results with the additional variables in order to compare with previous estimates4. As we can see in the Table 2, the coefficients estimated (using again 16 quadrature points) for those variables are not significant and all initial coefficients in the model remain approximately the same.

Table 2. Simulated data set with added variables estimation’s results.

Estimated standard deviations for estimated coefficients are given within parenthesis. ***: significant at the 1% level. **: significant at the 5% level.

5.3. Impact of Number of Quadrature Points on Estimated Results

As the accuracy of the method depends on the number of quadrature points used for the likelihood calculation, we propose an assessment of how it affects the results when this number increases. For doing so, we fit the same model with different numbers of quadrature points and we calculate the relative difference in log-likelihood and in estimated parameters.

We fit some models by using the same simulated relationship between variables as in Section 5.1.

The results are displayed in the Table 3 for dynamic equations and in the Table 4 for initial conditions equations and errors terms covariance matrix structure.

As we can see from Table 3 and Table 4, by increasing the number of quadrature points the changes in results decline and the relative differences are around 0.01% for significant coefficients and 0.1% or at most 1% for non significant coefficients. After 16 quadrature points, the relative differences in log-likelihood and in estimated coefficients become fewer as we increase the number of quadrature points. The estimations with 22 quadrature points are closer to those with 24 quadrature points than the others. So when we increase the number of quadrature points the changes in estimated coefficients are not significant but the computing time grows up exponentially. For these models, estimation time on an i5 core computer at 2.5 GHz with 6 GB of RAM memory for the different number of quadrature points are given in Table 5.

6. Conclusions

This paper describes the bivariate dynamic probit model with endogenous initial conditions starting by justifying the econometric specification of the model, giving the estimation method and its requirements and ending by presenting a robustness analysis. We calculate the derivatives of the log-likelihood function (the gradient) with respect to the 13 parameters in the model. This is the main contribution of our research as many programs use numerical computation of the gradient vector instead of encoding the mathematically derived expression of

Table 3. Impact of the number of quadrature points on estimation results. Part A.

Estimated standard deviations for estimated coefficients are given within parenthesis. ***: significant at the 1% level. **: significant at the 5% level.

Table 4. Impact of the number of quadrature points on estimation results. Part B.

Estimated standard deviations for estimated coefficients are given within parenthesis. ***: significant at the 1% level. **: significant at the 5% level. *: significant at the 10% level.

Table 5. Computing time for different number of quadrature points.

the gradient. Furthermore, for the use of the adaptative Gauss-Hermite quadrature, we also calculate the Hessian matrix with respect to individual random effects vector.

The implementation has been done using Stata software. We wrote 2 ado-files for this purpose. We use Stata’s d1 method for the maximization process. For the use of this method, we implement the gradient vector for the 13 parameters and we also implement the Hessian matrix with respect the random effects vector in order to use the adaptative Gauss-Hermite quadrature. We also wrote two others ado-files for the estimation of the bivariate probit for panel data and the bivariate dynamic probit without initial conditions for panel data. These ado-files are written using the same method (Stata’s d1 method) with the adaptative Gauss-Hermite quadrature. These ado-files are available upon request.

Due to the fact that the integration is bi-dimensional, estimation time is very high and stills increasing when the number of quadrature points or the number of observation or the number of explanatory variable increase. For an estimated model, one should insure that when increasing the number of quadrature point, the computed results don’t change significantly before using them. It means that the relative difference in the results must be around 0.1% or fewer, and if so, we can conclude that the results remain stable when increasing the number of quadrature points. And, it means that there is no need to increase the number of quadrature points that will increase computing time but will not improve significantly the results. However, increasing the number of quadrature points also increases the computation time. One way for major improvement of the program is the use of multi-core (parallel) computing scheme. This scheme allows to make the computation of the contributions to the likelihood (Equation (23)) at each quadrature point separately and simultaneously on several cores. It has the advantage to save time since the contributions are computed in the same time.

Finally, our method gives reasonable computing durations with real dataset. In [12] , we make use of the full SIP data set with 10,569 individuals and 255,206 observations.

Acknowledgements

We acknowledge the Centre Maurice Halbwachs (Réseau Quetelet) for access to the SIP 2007 data set (Santé et itinéraire professionnel-2007. DARES producteur. Centre Maurice Halbwachs diffuseur).

NOTES

1For each section, specifics notations are down at the beginning of the section. Otherwise, in general ${f\left(x\right)|}_{x=a}$ denote the value of the function or the matrix f at the point a. When not specify, $|a|$ denote the integer part of the scalar a.

2Notice that even without this factor, one can use the Gauss-Hermite quadrature by using a straightforward transformation that is to multiply and divide the integrand $f\left(x\right)$ by a Gaussian factor $\mathrm{exp}\left(-{x}^{2}\right)$ .

3In this section, $\varphi \left(x\right)$ denotes the univariate normal density function, $\varphi \left(x,y,\rho \right)$ denote the bivariate normal density with correlation r, ${\Phi }_{1}\left(x\right)$ denote the univariate normal probability function, and ${\Phi }_{2}\left(x,y,\rho \right)$ denote the bivariate normal probability function with correlation ρ.

4We do the same with columns 1’, 2’ of Table 1 and Table 2 (new results are in column 3’) and with columns 4 and 5 of both tables (new results in column 6).

Conflicts of Interest

The authors declare no conflicts of interest.

Cite this paper

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 Letters, 8, 1257-1278. doi: 10.4236/tel.2018.86083.

 [1] Adams, P., Hurd, M.D., McFadden, D., Merril, A. and Ribeiro, T. (2003) Healthy, Wealthy, and Wise? Tests for Direct Causal Paths between Health and Socioeconomic Status. Journal of Econometrics, 112, 3-56. [2] Heckman, J.J. (1979) The Incidental Parameters Problem and the Problem of Initial Conditions in Estimating: A Discrete Time-Discrete Data Stochastic Process and Some Monte Carlo Evidence. Graduate School of Business and Department of Economics, University of Chicago, Chicago. [3] Alessie, R., Hochguertel, S. and Van Soest, A. (2004) Ownership of Stocks and Mutual Funds: A Panel Data Analysis. The Review of Economics and Statistics, 86, 783-796. [4] Gourieroux, C. and Monfort, A. (1996) Simulation-Based Econometric Methods. CORE Lectures, Oxford University Press, Oxford. [5] Naylor, J.C. and Smith, A.F.M. (1982) Applications of a Method for the Efficient Computation of Posterior Distributions. Applied Statistics, 31, 214-225. [6] Liu, Q. and Pierce, D. (1994) A Note on Gauss-Hermite Quadrature. Biometrika, 83, 624-629. https://doi.org/10.2307/2337136 [7] Jackel, P. (2005) A Note on Multivariate Gauss-Hermite Quadrature. ABN AMRO, London. [8] Granger, C.W.J. (1969) Investigating Causal Relations by Econometric Models and Cross-Spectral Methods. Econometrica, 37, 428-438. https://doi.org/10.2307/1912791 [9] Nair-Reichert, A. and Weinhold, D. (2000) Causality Tests for Cross-Country Panels: New Look at FDI and Economic Growth in Developing Countries. Oxford Bulletin of Economics and Statistics, 63, 153-171. [10] Gould, W., Pitblado, J. and Poi, B. (2010) Maximum Likelihood Estimation with Stata. 4th Edition, Stata Press, College Station. [11] Miranda, A. (2007) Migrant Networks, Migrant Selection and High School Graduation in Mexico, IZA dp No. 3204, Institute for the Study of Labor, Bonn. [12] Delattre, E., Moussa, R. and Sabatier, M. (2015) Health Condition and Job Status Interactions: Econometric Evidence of Causality from a French Longitudinal Survey. THEMA Working Paper, No. 19, 1-28.