*! version 1.0.0  12apr2004 M. Lokshin, Z. Sajaia              (SJ4-3: st0000)
#delim ;
capture program drop movestay_d2;
program define movestay_d2, eclass;
	args todo b lnf g negH g1 g2 g3 g4 g5 g6 g7;

	tempname lns1 lns2 theta1 theta2 xb1 xb2 xb3;

	mleval `xb1' 	=`b', eq(1);
	mleval `xb2' 	=`b', eq(2);
	mleval `xb3' 	=`b', eq(3);
	mleval `lns1'	=`b', eq(4) scalar;
	mleval `lns2'	=`b', eq(5) scalar;
	mleval `theta1'	=`b', eq(6) scalar;
	mleval `theta2'	=`b', eq(7) scalar;

	quietly {;

	tempname sig1 sig2 rho1 rho2 rr1 rr2;

	scalar `sig1' 	= exp(`lns1');
	scalar `sig2' 	= exp(`lns2');
	scalar `rho1' 	= tanh(`theta1');
	scalar `rho2' 	= tanh(`theta2');
	scalar `rr1'	= 1/sqrt(1-`rho1'^2);
	scalar `rr2'	= 1/sqrt(1-`rho2'^2);

	tempvar eta1 eta2 lf1 eps1 eps2;

	generate double `eta1' = (`xb3' +  ($ML_y1 - `xb1') * `rho1'/`sig1')*`rr1';
	generate double `eta2' = (`xb3' +  ($ML_y2 - `xb2') * `rho2'/`sig2')*`rr2';
	generate double `eps1' = $ML_y1 - `xb1';
	generate double `eps2' = $ML_y2 - `xb2';

	tempname const1 const2;

	scalar const1=0.5*ln(2*_pi*`sig1'^2);
	scalar const2=0.5*ln(2*_pi*`sig2'^2) ;

	generate double `lf1' =	cond($ML_y3==1, ln(norm(  `eta1'))-const1-0.5*(`eps1'/`sig1')^2,
	                                        ln(norm(- `eta2'))-const2-0.5*(`eps2'/`sig2')^2);

	mlsum `lnf' =`lf1';

	if (`todo'==0 | `lnf'==.)  exit;  // calculate first derivative

	// derivative v.r. to x's in first regression equation
	tempname  deta_dx deta_dsig deta_drho M;

	generate double `deta_dx'   = -cond($ML_y3==1,`rho1'/`sig1'*`rr1',`rho2'/`sig2'*`rr2');
	generate double `deta_dsig' = -cond($ML_y3==1,`eps1'*`rho1'/`sig1'^2*`rr1',`eps2'*`rho2'/`sig2'^2*`rr2');

	generate double `deta_drho' =  cond($ML_y3==1,`eps1'*`rr1'/`sig1'+(`rho1'*(`xb3'+`rho1'*`eps1'/`sig1'))*`rr1'^3,
	                                              `eps2'*`rr2'/`sig2'+(`rho2'*(`xb3'+`rho2'*`eps2'/`sig2'))*`rr2'^3);

	generate double `M'         =  cond($ML_y3==1,normden(`eta1')/norm(`eta1'),-normden(`eta2')/norm(-`eta2'));

    // derivatives w.r. to x1 regression 1
	replace `g1' 	 =  cond($ML_y3==1,`M'*`deta_dx'+`eps1'/(`sig1'^2),0);

	// derivatives w.r. to x2 regression 2
    replace `g2' 	 =  cond($ML_y3==0,`M'*`deta_dx'+`eps2'/(`sig2'^2),0);

	// derivatives w.r. to x3 probit
	replace `g3'	 =  cond($ML_y3==1,`M'*`rr1',`M'*`rr2');

	// derivatives w.r. to sigma1
	replace `g4' 	 =  cond($ML_y3==1,`M'*`deta_dsig'-(1/`sig1')+(`eps1'^2)/(`sig1'^3),0);

	// derivatives w.r. to sigma2
	replace `g5' 	 =  cond($ML_y3==0,`M'*`deta_dsig'-(1/`sig2')+(`eps2'^2)/(`sig2'^3),0);

	// derivatives w.r to rho1
	replace `g6'	 =  cond($ML_y3==1,`M'*`deta_drho',0);

	// derivatives w.r to rho2
	replace `g7'	 =  cond($ML_y3==0,`M'*`deta_drho',0);


	tempname d1 d2 d3 d4 d5 d6 d7 drho1_dtheta1 drho2_dtheta2;

  	scalar `drho1_dtheta1' 	= 4*exp(2*`theta1')/((1+exp(2*`theta1'))^2);	  // derivative of rho w.r.t theta
	scalar `drho2_dtheta2' 	= 4*exp(2*`theta2')/((1+exp(2*`theta2'))^2);	  // derivative of rho w.r.t theta

    mlvecsum `lnf' `d1' = `g1' 			    	 , eq(1);
	mlvecsum `lnf' `d2' = `g2' 		    		 , eq(2);
	mlvecsum `lnf' `d3' = `g3' 		     		 , eq(3);
	mlvecsum `lnf' `d4' = `g4' * `sig1'          , eq(4);
	mlvecsum `lnf' `d5' = `g5' * `sig2'          , eq(5);
	mlvecsum `lnf' `d6' = `g6' * `drho1_dtheta1' , eq(6);
	mlvecsum `lnf' `d7' = `g7' * `drho2_dtheta2' , eq(7);

	matrix `g' = (`d1' , `d2' , `d3' , `d4' , `d5' , `d6', `d7' );

	if (`todo'==1 | `lnf'==. )  exit;   // calculate second derivative

	tempvar  DM d2eta_drho2 d2rho1_dtheta12 d2rho2_dtheta22;

	generate double `DM'          = -cond($ML_y3==1,`M'*(`eta1'+`M'),`M'*(`eta2'+`M'));

    generate double `d2eta_drho2' =  cond($ML_y3==1,`rr1'^3*(2*`rho1'*`eps1'/`sig1'+(`xb3'+`rho1'*`eps1'/`sig1')*(1+3*(`rho1'*`rr1')^2)),
                                                    `rr2'^3*(2*`rho2'*`eps2'/`sig2'+(`xb3'+`rho2'*`eps2'/`sig2')*(1+3*(`rho2'*`rr2')^2)));

	scalar `d2rho1_dtheta12'      = 8*exp(2*`theta1')*(1-exp(2*`theta1'))/(1+exp(2*`theta1'))^3;	  // 2nd derivative of rho w.r.t theta
	scalar `d2rho2_dtheta22'      = 8*exp(2*`theta2')*(1-exp(2*`theta2'))/(1+exp(2*`theta2'))^3;

	tempvar	 	g11 g12 g13 g14 g15 g16 g17
			    	g22 g23 g24 g25 g26 g27
			        	g33 g34 g35 g36 g37
				    		g44 g45 g46 g47
								g55 g56 g57
			 		    			g66 g67
										g77
				;

	tempvar	 	d11 d12 d13 d14 d15 d16 d17
			    	d22 d23 d24 d25 d26 d27
						d33 d34 d35 d36 d37
		                    d44 d45 d46 d47
								d55 d56 d57
					    			d66 d67
										d77
				;

	generate double `g11' = - cond($ML_y3==1,`DM'*`deta_dx'^2-1/`sig1'^2,0);
    generate double `g13' = - cond($ML_y3==1,`DM'*`deta_dx'*`rr1',0);
    generate double `g14' = - cond($ML_y3==1,(`DM'*`deta_dsig'-`M'/`sig1')*`deta_dx'-2*`eps1'/`sig1'^3,0);
    generate double `g16' = - cond($ML_y3==1,(`DM'*`deta_drho'+`M'*(1/`rho1'+`rho1'*`rr1'^2))*`deta_dx',0);

	generate double `g22' = - cond($ML_y3==0,`DM'*`deta_dx'^2-1/`sig2'^2,0);
    generate double `g23' = - cond($ML_y3==0,`DM'*`deta_dx'*`rr2',0);
    generate double `g25' = - cond($ML_y3==0,(`DM'*`deta_dsig'-`M'/`sig2')*`deta_dx'-2*`eps2'/`sig2'^3,0);
    generate double `g27' = - cond($ML_y3==0,(`DM'*`deta_drho'+`M'*(1/`rho2'+`rho2'*`rr2'^2))*`deta_dx',0);

    generate double `g33' = - cond($ML_y3==1,`DM'*`rr1'^2,`DM'*`rr2'^2);
    generate double `g34' = - cond($ML_y3==1,`DM'*`deta_dsig'*`rr1',0);
    generate double `g35' = - cond($ML_y3==0,`DM'*`deta_dsig'*`rr2',0);;
    generate double `g36' = - cond($ML_y3==1,`DM'*`deta_drho'*`rr1'+`M'*`rho1'*`rr1'^3,0);
    generate double `g37' = - cond($ML_y3==0,`DM'*`deta_drho'*`rr2'+`M'*`rho2'*`rr2'^3,0);

	generate double `g44' = - cond($ML_y3==1,(`DM'*`deta_dsig'^2-2*`M'*`deta_dsig'/`sig1'+1/`sig1'^2-3*`eps1'^2/`sig1'^4)*`sig1',0)-`g4';
    generate double `g46' = - cond($ML_y3==1,`DM'*`deta_drho'*`deta_dsig'+`M'*`deta_dsig'*(1/`rho1'+`rho1'*`rr1'^2),0);

    generate double `g55' = - cond($ML_y3==0,(`DM'*`deta_dsig'^2-2*`M'*`deta_dsig'/`sig2'+1/`sig2'^2-3*`eps2'^2/`sig2'^4)*`sig2',0)-`g5';
    generate double `g57' = - cond($ML_y3==0,`DM'*`deta_drho'*`deta_dsig'+`M'*`deta_dsig'*(1/`rho2'+`rho2'*`rr2'^2),0);

    generate double `g66' = - cond($ML_y3==1,(`DM'*`deta_drho'^2+`M'*`d2eta_drho2')*`drho1_dtheta1'^2,0)-`g6'*`d2rho1_dtheta12';

    generate double `g77' = - cond($ML_y3==0,(`DM'*`deta_drho'^2+`M'*`d2eta_drho2')*`drho2_dtheta2'^2,0)-`g7'*`d2rho2_dtheta22';

	mlmatsum `lnf' `d11'=`g11',  			  			eq(1);
	mlmatsum `lnf' `d12'= 0, 				  			eq(1,2);
	mlmatsum `lnf' `d13'=`g13', 						eq(1,3);
	mlmatsum `lnf' `d14'=`g14'*`sig1',          		eq(1,4);
	mlmatsum `lnf' `d15'= 0, 				  			eq(1,5);
	mlmatsum `lnf' `d16'=`g16'*`drho1_dtheta1',			eq(1,6);
	mlmatsum `lnf' `d17'= 0, 				  			eq(1,7);

	mlmatsum `lnf' `d22'=`g22', 						eq(2);
	mlmatsum `lnf' `d23'=`g23', 				 		eq(2,3);
	mlmatsum `lnf' `d24'= 0, 				  			eq(2,4);
	mlmatsum `lnf' `d25'=`g25'*`sig2',       			eq(2,5);
	mlmatsum `lnf' `d26'= 0, 				  			eq(2,6);
	mlmatsum `lnf' `d27'=`g27'*`drho2_dtheta2', 		eq(2,7);

	mlmatsum `lnf' `d33'=`g33', 				  		eq(3);
	mlmatsum `lnf' `d34'=`g34'*`sig1',       			eq(3,4);
	mlmatsum `lnf' `d35'=`g35'*`sig2',	        		eq(3,5);
	mlmatsum `lnf' `d36'=`g36'*`drho1_dtheta1',			eq(3,6);
	mlmatsum `lnf' `d37'=`g37'*`drho2_dtheta2', 		eq(3,7);

	mlmatsum `lnf' `d44'=`g44'*`sig1',              	eq(4);
	mlmatsum `lnf' `d45'= 0, 				  			eq(4,5);
	mlmatsum `lnf' `d46'=`g46'*`sig1'*`drho1_dtheta1',	eq(4,6);
	mlmatsum `lnf' `d47'= 0, 				  			eq(4,7);

	mlmatsum `lnf' `d55'=`g55'*`sig2',                  eq(5);
	mlmatsum `lnf' `d56'= 0, 				  			eq(5,6);
	mlmatsum `lnf' `d57'=`g57'*`sig2'*`drho2_dtheta2',  eq(5,7);

	mlmatsum `lnf' `d66'=`g66',	                	  	eq(6);
	mlmatsum `lnf' `d67'= 0, 		                  	eq(6,7);

	mlmatsum `lnf' `d77'=`g77',              		  	eq(7);

	matrix `negH'=(	`d11' , `d12' ,`d13' ,`d14' ,`d15' ,`d16' ,`d17'   \
					`d12'', `d22' ,`d23' ,`d24' ,`d25' ,`d26' ,`d27'   \
					`d13'', `d23'',`d33' ,`d34' ,`d35' ,`d36' ,`d37'   \
					`d14'', `d24'',`d34'',`d44' ,`d45' ,`d46' ,`d47'   \
					`d15'', `d25'',`d35'',`d45'',`d55' ,`d56' ,`d57'   \
					`d16'', `d26'',`d36'',`d46'',`d56'',`d66' ,`d67'   \
					`d17'', `d27'',`d37'',`d47'',`d57'',`d67'',`d77'   );
}; // end quietly
end;