/*

     Computing the selection variables from the bivariate probit model.

     The univariate probit model has the ;Hold(IMR = name) option avail-

     able for retaining the inverse mills ratio for corrections for sample

     selection.  There is no direct counterpart in the bivariate probit

     model, but it is possible to obtain the variables by other means.



     The model is as follows:  The bivariate probit model is



                z1  =  sign ( a1'w1  +  e1 )

                z2  =  sign ( a2'w2  +  e2 )



                e1,e2 ~ Normal[0,0,1,1,rho]



                y   =  b'x  +  u, u,e1,e2 ~ trivariate normal

                                  Cov[u,ej]  =  sigma * rj, j=1,2.

                                             =  dj



                [y,x] are observed when z1=1,z2=1 (or some other pairing).



      The sample selection corrected regression is



                E[y|z1=1,z2=1] = b'x +  d1*L1 + d2*L2



      where L1 and L2 are the counterparts to the inverse mills ratio in

      the univariate probit model.  These two variables are computed as 

      follows:

                Let  q1  =  2*z1 - 1,   q2 = 2*z2 - 1

                     c1  =  q1 * a1'w1, c2 = q2 * a2'w2

                     v1  =  (c2 - rho * v1)/sqr(1 - rho*rho)

                     v2  =  (c1 - rho * v2)/sqr(1 - rho*rho)



      Then,

                        L1  =  q1 * f(c1) * F(v1) / B(c1,c2,q1*q2*rho)

                        L2  =  q2 * f(c2) * F(v2) / B(c1,c2,q1*q2*rho)



      (These account for all pairings of z1 and z2, not just 1,1.)



      A set of LIMDEP commands that can be used to obtain these two

      variables is given below.  To use this program, the user needs only

      to insert the appropriate definitions of X1 and X2.  The

      procedure is self contained.  To use it, the command is



                Execute ; Proc = Lambda2(X1,X2,z1,z2,l1,l2) $



      where z1 and z2 are the two Lhs variables in your bivariate

      probit model (you use your names), and in place of l1 and l2, you

      insert the names you want to use for the two "Lambda" variables.



*/

        Namelist  ;  X1  =  ... the Rh1 of your bivariate probit

                  ;  X2  =  ... the Rh2 of your bivariate probit $



        Procedure = Lambda2(W1,W2,Z1,Z2,L1,L2) $

          Bivariate ; Lhs = Z1,Z2 

                    ; Rh1 = W1 

                    ; Rh2 = W2 $

          Calculate ; K1 = Col(W1) 

                    ; J1 = K1+1 

                    ; J2 = K1+Col(W2) $

          Matrix    ; B1 = Part(B,1,K1) ; B2 = Part(B,J1,J2) $

          Create    ; q1 = 2*z1-1       ; q2 = 2*z2-1

                    ; c1 = q1 * b1'W1   ; c2 = q2 * b2'w2

                    ; v1  =  (c2 - rho * c1)/sqr(1 - rho*rho)

                    ; v2  =  (c1 - rho * c2)/sqr(1 - rho*rho) $

          Namelist  ; V = v1,v2 $

          Create    ; L1 = q1 * N01(c1) * Phi(v1) / Bvn(V,(q1*q2*rho))

                    ; L2 = q2 * N01(c2) * Phi(v2) / Bvn(V,(q1*q2*rho)) $

        Endprocedure

 

        Execute ; Proc = Lambda2(X1,X2,z1,z2,Lambda1,Lambda2) $