OPTIONS CRT;
name tob2 "Bivariate Tobit Estimation using ML" ;
?
? TSP 4.2A (Nov. 1992)
? Author: Clint Cummins
?
? Bivariate Tobit, using Cholesky factorization of Sigma-inverse
? (the factorization is required to keep sigma positive definite).
? This model is described in Maddala, Limited-Dependent and Qualitative
? Dependent Variables in Econometrics, pp. 205-208.
?
FRML TOB2 LOGL =
 S1*[LOG(1-GAMMA1*GAMMA2) + LOG(D1) + LOG(D2) - LOG(6.2831853)
     - (E1*E1*S11 + 2*E1*E2*S11*L + E2*E2*S22)/2] +
 S2*(-LOG(SIGMA1) + LNORM(E10/SIGMA1) +
     LCNORM[{-G2Y1-X2B2 + E10*S11*L/S22}*SQRT(S22)] ) +
 S3*( LOG(D2    ) + LNORM(E20*D2    ) +
     LCNORM[{-G1Y2-X1B1 + E20*L        }*D1       ] ) +
 S4*LOG[CNORM(-X1B1S)*CNORM(-X2B2S) + NORM(X1B1S)*NORM(X2B2S)*{
         RHO + RHO**2*X1B1S*X2B2S/2 + RHO**3*(X1B1S**2 - 1)*(X2B2S**2 - 1)/6 }
       ];
FRML S11 D1*D1          ;   ? (1,1) element of Sigma-inverse
FRML S22 S11*L*L + D2*D2;   ? (2,2) element of Sigma-inverse
FRML ERHO RHO = -L*D1/SQRT(S22);  ? original parameters in terms of Cholesky
FRML ES1 SIGMA1 = SQRT(S22)/(D1*D2);
FRML ES2 SIGMA2 = 1/D2;
FRML X1B1S X1B1*D1*D2/SQRT(S22);   ? standardized X1B1 = X1B1/sigma1
FRML X2B2S X2B2*D2;
FRML G1Y2 GAMMA1*Y2;   ? user can change the names of Y1 and Y2 in next 4 eqns
FRML G2Y1 GAMMA2*Y1;
FRML E1 Y1-G1Y2-X1B1;
FRML E2 Y2-G2Y1-X2B2;
FRML E10 Y1-X1B1;   ? E1 | Y2=0
FRML E20 Y2-X2B2;   ? E2 | Y1=0
?
?  LHS and RHS variables -- user-defined
?  Note that exclusion restrictions are required for identification
?  (if GAMMA1 and GAMMA2 are non-zero)
?
FRML Y1 MY1;   ? real names for Y1 and Y2
FRML Y2 MY2;
?  Model for Latent Variables.
FRML X1B1  B10 + B11*X1 + B12*X2 + B13*X3;   
FRML X2B2  B20 + B21*X1 + B22*X2 +        B23*X4;
LIST BS  B10-B12 B20-B22;           ? List of RHS parameters. 
?
EQSUB TOB2 E1 E2 E10 E20 G1Y2 G2Y1 X1B1S X2B2S 
           ERHO ES1 S22 S11 X1B1 X2B2 Y1 Y2;
EQSUB ERHO S22 S11;  
EQSUB ES1 S22 S11;
PARAM BS GAMMA1 GAMMA2 D1 D2 L;
?
?  Estimation example with data generated from known model
?
smpl 1,500;
?
?  Covariance parameter values
?
read(nrow=2,type=sym) covu;
.9
.5 1.5 ;
?
?  RHS parameter values
param b10 .2 b11 1.5 b12 2 b13 1.3
      b20 .3 b21 3.5 b22 4 b23 2.3
      gamma1 .3 gamma2 .5;
?
?  Draw residuals and RHS variables.
random(seedin=178948,vcov=covu) u1 u2;
random x1-x4;
?
?  Generate dependent variables from simultaneous model.
frml ey1 my1 = pos(g1y2 + x1b1 + u1);
frml ey2 my2 = pos(g2y1 + x2b2 + u2);
eqsub ey1 g1y2 x1b1 y2;
eqsub ey2 g2y1 x2b2 y1;
my1 = .1;  my2 = .1;     ? these starting values seem important
siml(endog=(my1,my2),tag=none,silent,noprnsim,maxsqz=20) ey1 ey2;
print @ifconv;  ? check for any non-converged observations.
?
?  Generate subset dummies
S1 = MY1>0 & MY2>0;
S2 = MY1>0 & MY2=0;
S3 = MY1=0 & MY2>0;
S4 = MY1=0 & MY2=0;
msd S1-S4;  ? check sample counts
?
?  Initialize covariance parameters.
?
inv covu si;   
yldfac si dsig tsig;
set D1 = sqrt(dsig(1,1));
set D2 = sqrt(dsig(2,2));
set L = tsig(1,2);
?
ML(HITER=N,HCOV=NBW) tob2;
?
? recover values + SEs for original covariance parameters
ANALYZ ERHO ES1 ES2;
?
?  Documentation for Cholesky factorization of Sigma-inverse:
?
?  instead of -.5*log(det(sigma)), make simplification:
?   det(sigma) = 1/det(sigma-inverse) = 1/D1*D1*D2*D2*...
?   log(det(sigma)) = -[2*log(D1) + 2*log(D2) + ...]
?                   = -2*ldetd     (and -.5*(-2) = 1)
?                  or could use log(D1*D1) if there are negative problems
?
?  parameterize sigma inverse in terms of its LDL' factorization,
?  so that its determinant is simplified
?  s.. = sigma inverse element (lower triangle)
?  D.  = square root of D element (.,.)
?  L   = L element (2,1) (unit triangular:  ones on diagonal,
?                                           zeros above diagonal)
?  2 x 2 case example:
?
?  sigma inverse = s11 s21 = [ 1 0 ] * [ D1*D1    0   ] * [ 1 L ]
?                  s21 s22 = [ L 1 ] * [   0    D2*D2 ] * [ 0 1 ]
?
?  Conditional densities (for LCNORM with S1 or S2):
?  e1|e2
?  mean = rho*(sigma1/sigma2)*e2 = -L
?  var  = sigma1**2*(1-rho**2) = 1/(D1*D1)
?  e2|e1
?  mean = rho*(sigma2/sigma1)*e1 = -L*D1*D1/S22
?  var  = sigma2**2*(1-rho**2) = 1/S22
stop ; end ;