/*======================================================================

  This routine estimates Maddala's Model 3, page 245. Simultaneous 

  equations:  

              y1* = c1 y2* + b1'x1 + u1

              y2* = c2 y1* + b2'x2 + u2



  y1* is directly observed, but only the sign of y2* is observed.

  (Setting up a probit model)  Variance of u1 is sigma-1 squared. 

  Since y2* is not observed, save for sign, sigma-2 squared is 

  normalized to 1.0.  (Maddala's additional material is a waste of time.

  Do this directly.) The correlation is r12.



  Steps:  Reduced forms, in which X is the union of X1 and X2 are  



              y1* = p1'X + v1

              y2* = p2'X + v2.





          By solution, the variance of v1 is 



              t1-squared = (s1sq + c1sq + 2c1s1r12)/(1-c1c2)sq



          1.  Regress y1 on X (all) to get p1 and t1, compute p1'X

          2.  Probit of y2 on X (all) to get p2.  Compute p2'X

          3.  OLS of y1 on (p2'X) and X1 to get c1 and b1, s1sq

          4.  Probit of y2 on (p1'x) and X2 to get c2 and b2.

          5.  Use estimates of t1, c1, c2, and s1sq to estimate r12

          6.  Compute c =    s1sq - 2s1c1r12

          7.  Compute d = c2sq s1sq - 2s1c2r12

          8.  Compute estimated covariance matrices.  (Correct errors

              in Maddala for the second.  His G is supposed to be

              Z2'X, which in his notation would be G'X'X.  In his

              expression for Var{alpha2}, inv(V0) should be just V0.

              Easier, in this expression, we set V0 equal to the inverse

              of the V0 that appears in Var{alpha1}

*/

?  Exogenous variables in the two equations.  X = union.

Namelist  ; X1 = ...

          ; X2 = ...

          ; X  = ...  $

?  Copies of endogenous variables in two equations.

Create    ; y1 = ...

          ; y2 = ... $

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

? This form of the procedure declaration is useable only in version 7 =

Procedure=Model3(y1,y2,x1,x2,x) $

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

? First reduced form regression

Regr      ; Lhs = y1 ; Rhs = X ; Keep = P1X $

? Retrieve reduced form variance

Calc      ; t1sq = ssqrd $

? Second reduced form

Probit    ; Lhs = y2 ; Rhs = X $

? Estimate continuous variable.

Create    ; P2X = Dot(X,b) $ 

? Structural variables in two equations

Namelist  ; Z1 = P2X,X1 $

Namelist  ; Z2 = P1X,X2 $

? Retrieve inverse of asymptotic covariance for probit estimates

Matrix    ; V0I = <Varb> $

? First structural variance

Regr      ; Lhs = y1 ; Rhs = Z1 $

? Trap coefficient estimates and disturbance variance estimators

Matrix    ; alpha1 = b $

Calc      ; gamma1 = b(1) 

          ; Sigmasq1 = ssqrd 

          ; sigma1 = s $

? Second structural estimates

Probit    ; Lhs = Y2 ; Rhs = Z2 $

? Retrieve coefficient estimates and ancillary parameters

Matrix    ; alpha2 = b $

Calc      ; gamma2 = b(1) 

          ; rho12 = (2*gamma1*sqr(sigmasq1)) *

                    ((1-gamma1*gamma2)^2*t1sq - gamma1^2 - sigmasq1)

          ; c =            sigmasq1 - 2*gamma1*rho12*sigma1

          ; d = gamma2^2 * sigmasq1 - 2*gamma2*rho12*sigma1 $

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

? Asymptotic covariance matrices.

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

Matrix ; Z1Z1I=<Z1'Z1>

       ; Va1 = c*Z1Z1I + {gamma1^2} * Z1Z1I * Z1'X * <V0I> * X'Z1 * Z1Z1I

       ; V0I   = <V0I>

       ; Z2XVI = Z2'X * V0I * X'Z2

       ; Z2XVI = <Z2XVI>

       ; Va2 = Z2XVI + d*Z2XVI * Z2'X * V0I * <X'X> * V0I * X'Z2 * Z2XVI

       ; Stat(Alpha1,Va1)

       ; Stat(Alpha2,Va2) $

Endprocedure



/* 

        To use the procedure, the command is



        Execute ; Proc = Model3(y1, y2, x1, x2, x) $



        where the arguments are the variables and namelists that you have

        defined earlier.

*/