? *******  allowing for missing values in 1st or 2nd period  ******

?          inc1,inc2 continuous dependant variables > 0

?          d1 = 1 if inc1 > 0, 0 otherwise
?          d2 = 1 if inc2 > 0, 0 otherwise
?          d3 = 1 if inc1 missing, 0 otherwise
?          d4 = 1 if inc2 missing, 0 otherwise

options crt;
smpl 1 382;
supres @smpl;
load (file='../data/reg.dat')
id d1 d2 d3 d4 inc1 inc2 x1-x27 y1-y36;
inc1 = inc1/1000;inc2 = inc2/1000;
linc1 = log(D1*inc1+(1-D1));linc2 = log(D2*inc2+(1-D2));
v1 = 1;v2 = x3;v3 = x4/1000;v4 = x5;v5 = x6;v6 = x8;v7 = x9;v8 = x10;
v9 = x11;v10 = x12;v11 = x16/(x15*13);v12 = x19;v13 = x20;v14 = 1*x27;
v15 = x23*x27;v16 = x25*x27;v17 = x26*x27;
w1 = 1;w2 = y3;w3 = y4/1000;w4 = y5;w5 = y6;w6 = y8;w7 = y9;w8 = y10;
w9 = y11;w10 = y12;w11 = y16/(y15*13);w12 = y19;w13 = y20;w14 = 1*y27;
w15 = y23*y27;w16 = y25*y27;w17 = y26*y27;

list xvar1 v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12 v13 v14
v15 v16 v17;
list xvar2 w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12 w13 w14
w15 w16 w17;

list y1par = a1-a17;list I1par = c1-c17;
list y2par = b1-b17;list I2par = e1-e17;
form (coefpref=A) EQ1 XB1 XVAR1;form (coefpref=C) EQ3 XG1 XVAR1;
form (coefpref=B) EQ2 XB2 XVAR2;form (coefpref=E) EQ4 XG2 XVAR2;

frml L1  p1;

frml L2  norm((linc1-xb1)/s1)*p2/(s1*inc1*d1+(1-d1));

frml L3  norm((linc2-xb2)/s2)*p3/(s2*inc2*d2+(1-d2));

frml L4  (1/((1-d1)+(1-d2)+d1*d2*inc1*inc2*s1*s2*2*3.1415927*sqrt(1-rho*rho)))*
exp{(-0.5/(1-rho*rho))*[((linc1-xb1)/s1)**2
- 2*rho*((linc1-xb1)/s1)*((linc2-xb2)/s2) + ((linc2-xb2)/s2)**2 ]}*p4;

frml L5 cnorm(-XG1);

frml L6 (1/(s1*inc1*d1+(1-d1)))*norm[(linc1-xb1)/s1]
*cnorm[{xg1+(s13/(s1*s1))*(linc1-xb1)}/sqrt(1-(s13/s1)**2)];

frml L7 cnorm(-XG2);

frml L8 (1/(s2*inc2*d2+(1-d2)))*norm[(linc2-xb2)/s2]*
cnorm[{xg2+(s24/(s2*s2))*(linc2-xb2)}/sqrt(1-(s24/s2)**2)];

frml p1 cnorm(-xg1)*cnorm(-xg2) + norm(-xg1)*norm(-xg2)*
(su + su**2*(-xg1)*(-xg2)/2 + su**3*((-xg1)**2-1)*
((-xg2)**2-1)/6);

frml p2 cnorm(mu21)*cnorm(mu22) + norm(mu21)*norm(mu22)*
(rho2 + rho2**2*(mu21)*(mu22)/2 + rho2**3*((mu21)**2-1)*
((mu22)**2-1)/6);

frml p3 cnorm(mu31)*cnorm(mu32) + norm(mu31)*norm(mu32)*
(rho3 + rho3**2*(mu31)*(mu32)/2 + rho3**3*((mu31)**2-1)*
((mu32)**2-1)/6);

frml p4 cnorm(mu41)*cnorm(mu42) + norm(mu41)*norm(mu42)*
(rho4 + rho4**2*(mu41)*(mu42)/2 + rho4**3*((mu41)**2-1)*
((mu42)**2-1)/6);

frml rho S12/(S1*S2);
frml mu21  [XG1 + (S13/(S1*S1))*(LINC1-XB1)]/SQRT(1-(S13/S1)**2);
frml mu22 -[XG2 + (S14/(S1*S1))*(LINC1-XB1)]/SQRT(1-(S14/S1)**2);
frml mu31 -[XG1 + (S23/(S2*S2))*(LINC2-XB2)]/SQRT(1-(S23/S2)**2);
frml mu32  [XG2 + (S24/(S2*S2))*(LINC2-XB2)]/SQRT(1-(S24/S2)**2);
frml rho2 -(SU - S13*S14/(S1*S1))/SQRT[(1-(S13/S1)**2)*(1-(S14/S1)**2)];
frml rho3 -(SU - S23*S24/(S2*S2))/SQRT[(1-(S23/S2)**2)*(1-(S24/S2)**2)];
frml mu41 (XG1 + MU1PRIME)/S1PRIME;
frml mu42 (XG2 + MU2PRIME)/S2PRIME;
frml rho4  S12PRIME/(S1PRIME*S2PRIME);
frml del S1*S1*S2*S2-S12*S12;
frml mu1prime (S13*S2*S2-S23*S12)*(LINC1-XB1)/DEL + 
              (S23*S1*S1-S13*S12)*(LINC2-XB2)/DEL ;
frml mu2prime (S14*S2*S2-S24*S12)*(LINC1-XB1)/DEL + 
              (S24*S1*S1-S14*S12)*(LINC2-XB2)/DEL ;
frml s1prime sqrt{1 - (S13*S13*S2*S2 - 2*S13*S23*S12 + S23*S23*S1*S1)/DEL};
frml s2prime sqrt{1 - (S14*S14*S2*S2 - 2*S14*S24*S12 + S24*S24*S1*S1)/DEL};
frml s12prime SU-(S13*S2*S2*S14- S23*S12*S14- S13*S12*S24+ S23*S1*S1*S24)/DEL;

frml logl 
     (1-d3)*(1-d4)*[
     (1-d1)*(1-d2)*log(L1+1-(1-d3)*(1-d4))+
         d1*(1-d2)*log(L2+1-(1-d3)*(1-d4))+
         d2*(1-d1)*log(L3+1-(1-d3)*(1-d4))+
         d1*d2    *log(L4+1-(1-d3)*(1-d4))]+
     (1-d3)*d4*[(1-d1)*log(L5 + 1-(1-d3)*d4) + d1*log(L6 + 1-(1-d3)*d4)] +
     (1-d4)*d3*[(1-d2)*log(L7 + 1-(1-d4)*d3) + d2*log(L8 + 1-(1-d4)*d3)];

eqsub S1PRIME DEL;eqsub S2PRIME DEL; eqsub S12PRIME DEL;
eqsub MU1PRIME DEL EQ1 EQ2;eqsub MU2PRIME DEL EQ1 EQ2;
eqsub MU41 EQ3 MU1PRIME S1PRIME;eqsub MU42 EQ4 MU2PRIME S2PRIME;
eqsub RHO4 S12PRIME S1PRIME S2PRIME;
eqsub MU21 EQ1 EQ3;eqsub MU22 EQ1 EQ4;eqsub MU31 EQ2 EQ3;eqsub MU32 EQ2 EQ4;
eqsub P1 EQ3 EQ4;eqsub P2 MU21 MU22 RHO2;
eqsub P3 MU31 MU32 RHO3;eqsub P4 MU41 MU42 RHO4;
eqsub L1 P1;eqsub L2 EQ1 P2;eqsub L3 EQ2 P3;eqsub L4 EQ1 EQ2 P4 RHO;
eqsub L5 EQ3;eqsub L6 EQ1 EQ3;eqsub L7 EQ4;eqsub L8 EQ2 EQ4;
eqsub LOGL L1 L2 L3 L4 L5 L6 L7 L8;

param s1 s2 s12 s13 s14 s23 s24 su y1par y2par I1par I2par;

?**** probits and ols to get starting values ****;

select d3 = 0;probit d1 xvar1;
smpl 1 17;unmake @coef i1par;smpl 1 382;
select d3 = 0 & inc1 > 0;olsq linc1 xvar1;
smpl 1 17;unmake @coef y1par;set s1 = sqrt(@s2);smpl 1 382;
select d4 = 0;probit d2 xvar2;
smpl 1 17;unmake @coef i2par;smpl 1 382;
select d4 = 0 & inc2 > 0;olsq linc2 xvar2;
smpl 1 17;unmake @coef y2par;set s2 = sqrt(@s2);smpl 1 382;

set s12=.1;set s13=.1;set s14=.1;set su=.1;set s23=.1;set s24=.1;
ml (MAXIT=1000) LOGL;

smpl 1 76;
write (file='pol.p1') @coef;
write (file='pol.v1') @vcov;
write (file='pol.t1') @t;