?===================================================================

? A simultaneous equations model with two binary dependent         =

? variables.                                                       =

?                                                                  =

?                                                                  =

? User must define the data matrices lines:                        =

?                                                                  =

?        NAMELIST ; X1 = RHS in first equation, may include Y2     =                              =

?                 ; X2 = RHS in second equation, may include Y1    =                              =

?                 ; X  = full set of ALL exogenous variables in    =

?                        both equations.                           =

? User's dependent variables are Y1 and Y2.                        =

?                                                                  =

?        CREATE ; Y1 = ... ; Y2 = ... $                            =

?===================================================================



? Reduced form probits, saving the coefficients



       PROBIT   ; lhs = y1 ; rhs = X $ 

       MATRIX   ; pi1 = b $ 

       PROBIT   ; lhs = y2 ; rhs = X $ 

       MATRIX   ; pi2 = b $ 

       CREATE   ; y1f = dot(x1,pi1) 

                ; y2f = dot(x2,pi2) $ 



? Structural probits



       NAMELIST ; z1 = one,y2f,x1  

                ; z2 = one,y1f,x2 $ 

       PROBIT   ; lhs = y1 ; rhs = z1 $ 

       MATRIX   ; alpha1 = b $ 

       PROBIT   ; lhs = y2 ; rhs = z2 $ 

       MATRIX   ; alpha2 = b$ 



? Calculate covariance matrices, first for alpha1 then,

? symmetrically, for alpha2



       CREATE   ; q1 =  dot(z1,alpha1)   ; q2 =  dot(z2,alpha2)

                ; a1 = -lmm(q1)/phi(q1)  ; a2 = -lmm(q2)/phi(q2)

                ; capa1 = n01(q1) * a1   ; capa2 = n01(q2) * a2

                ; u1 = y1 - phi(q1)      ; u2 = y2 - phi(q2)

                ; v = a1 * a2 * u1 * u2 $  



       MATRIX   ; w1 = xcpm(z1,z1,capa1)

                ; w2 = xcpm(x,x,capa2)

                ; w3 = alpha1(1) * xcpm(z1,x,capa1)

                ; w4 = xcpm(x,z1,v) 

                ; w324 = w3 | sinv(w2) | w4

                ; w323 = w3 | sinv(w2) | w3'

                ; va1 = w1 ~ w324 ~ w324' & w323 

                ; va1 = sinv(w1) | va1 | sinv(w1) $



? Now reverse subscripts and do it again.



       MATRIX   ; w1 = xcpm(z2,z2,capa2)

                ; w2 = xcpm(x,x,capa1)

                ; w3 = alpha2(1) * xcpm(z2,x,capa2)

                ; w4 = xcpm(x,z2,v)

                ; w324 = w3 | sinv(w2) | w4

                ; w323 = w3 | sinv(w2) | w3'

                ; va2 =  w1 ~ w324 ~ w324' & w323

                ; va2 = sinv(w1) | va2 | sinv(w1) $ 



? Display results.



       MATRIX   ; Stat(alpha1,va1)  ;  Stat(alpha2,va2)  $