*! 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;