*! version 3.0.1 10Mar2008 M. Lokshin, Z. Sajaia *! adapted for mtefe by Martin E. Andresen, 19feb2018 #delim ; program define mtefe_movestay_lf2, eclass; args todo b lnf g1 g2 g3 g4 g5 g6 g7 negH; tempname lns1 lns0 theta0 theta1 xb0 xb1 xb3; mleval `xb0' =`b', eq(1); mleval `xb1' =`b', eq(2); mleval `xb3' =`b', eq(3); mleval `lns0' =`b', eq(4) scalar; mleval `lns1' =`b', eq(5) scalar; mleval `theta0' =`b', eq(6) scalar; mleval `theta1' =`b', eq(7) scalar; quietly {; loc sig0 = exp(`lns0'); loc sig1 = exp(`lns1'); loc rho0 = tanh(`theta0'); loc rho1 = tanh(`theta1'); loc rr0 = 1/sqrt(1-`rho0'^2); loc rr1 = 1/sqrt(1-`rho1'^2); tempvar eta0 eta1 eps0 eps1; generate double `eta0' = (`xb3' + ($ML_y1 - `xb0') * `rho0'/`sig0')*`rr0'; replace `eta0'= 37 if (`eta0'> 37); generate double `eta1' = (`xb3' + ($ML_y2 - `xb1') * `rho1'/`sig1')*`rr1'; replace `eta1'=-37 if (`eta1'<-37); generate double `eps0' = $ML_y1 - `xb0'; generate double `eps1' = $ML_y2 - `xb1'; loc const0=0.5*ln(2*_pi*`sig0'^2); loc const1=0.5*ln(2*_pi*`sig1'^2); replace `lnf' = cond($ML_y3==1, ln(normal( `eta1'))-`const1'-0.5*(`eps1'/`sig1')^2, ln(normal(- `eta0'))-`const0'-0.5*(`eps0'/`sig0')^2); if (`todo'==0 | `lnf'==.) exit; // calculate first derivative tempname deta_dx deta_dsig deta_drho M; generate double `deta_dx' = -cond($ML_y3==1,`rho1'/`sig1'*`rr1',`rho0'/`sig0'*`rr0'); generate double `deta_dsig' = -cond($ML_y3==1,`eps1'*`rho1'/`sig1'^2*`rr1',`eps0'*`rho0'/`sig0'^2*`rr0'); generate double `deta_drho' = cond($ML_y3==1,`eps1'*`rr1'/`sig1'+(`rho1'*(`xb3'+`rho1'*`eps1'/`sig1'))*`rr1'^3, `eps0'*`rr0'/`sig0'+(`rho0'*(`xb3'+`rho0'*`eps0'/`sig0'))*`rr0'^3); generate double `M' = cond($ML_y3==1,normalden(`eta1')/normal(`eta1'),-normalden(`eta0')/normal(-`eta0')); loc drho0_dtheta0 = 4*exp(2*`theta0')/((1+exp(2*`theta0'))^2); // derivative of rho w.r.t theta loc drho1_dtheta1 = 4*exp(2*`theta1')/((1+exp(2*`theta1'))^2); // derivative of rho w.r.t theta replace `g1' = cond($ML_y3==0,`M'*`deta_dx'+`eps0'/(`sig0'^2),0); replace `g2' = cond($ML_y3==1,`M'*`deta_dx'+`eps1'/(`sig1'^2),0); replace `g3' = cond($ML_y3==1,`M'*`rr1',`M'*`rr0'); replace `g4' = cond($ML_y3==0,(`M'*`deta_dsig'-(1/`sig0')+(`eps0'^2)/(`sig0'^3))*`sig0',0); replace `g5' = cond($ML_y3==1,(`M'*`deta_dsig'-(1/`sig1')+(`eps1'^2)/(`sig1'^3))*`sig1',0); replace `g6' = cond($ML_y3==0,`M'*`deta_drho'*`drho0_dtheta0',0); replace `g7' = cond($ML_y3==1,`M'*`deta_drho'*`drho1_dtheta1',0); if (`todo'==1 | `lnf'==. ) exit; // calculate second derivative tempvar DM d2eta_drho2 lf; generate double `DM' = -cond($ML_y3==1,`M'*(`eta1'+`M'),`M'*(`eta0'+`M')); generate double `d2eta_drho2' = cond($ML_y3==1,`rr1'^3*(2*`rho1'*`eps1'/`sig1'+(`xb3'+`rho1'*`eps1'/`sig1')*(1+3*(`rho1'*`rr1')^2)), `rr0'^3*(2*`rho0'*`eps0'/`sig0'+(`xb3'+`rho0'*`eps0'/`sig0')*(1+3*(`rho0'*`rr0')^2))); loc d2rho1_dtheta12 = 8*exp(2*`theta1')*(1-exp(2*`theta1'))/(1+exp(2*`theta1'))^3; // 2nd derivative of rho w.r.t theta loc d2rho0_dtheta02 = 8*exp(2*`theta0')*(1-exp(2*`theta0'))/(1+exp(2*`theta0'))^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 `g22' = - cond($ML_y3==1,`DM'*`deta_dx'^2-1/`sig1'^2,0); generate double `g23' = - cond($ML_y3==1,`DM'*`deta_dx'*`rr1',0); generate double `g25' = - cond($ML_y3==1,(`DM'*`deta_dsig'-`M'/`sig1')*`deta_dx'-2*`eps1'/`sig1'^3,0); generate double `g27' = - cond($ML_y3==1,(`DM'*`deta_drho'+`M'*(1/`rho1'+`rho1'*`rr1'^2))*`deta_dx',0); generate double `g33' = - cond($ML_y3==1,`DM'*`rr1'^2,`DM'*`rr0'^2); generate double `g35' = - cond($ML_y3==1,`DM'*`deta_dsig'*`rr1',0); generate double `g34' = - cond($ML_y3==0,`DM'*`deta_dsig'*`rr0',0); generate double `g36' = - cond($ML_y3==0,`DM'*`deta_drho'*`rr0'+`M'*`rho0'*`rr0'^3,0); generate double `g37' = - cond($ML_y3==1,`DM'*`deta_drho'*`rr1'+`M'*`rho1'*`rr1'^3,0); generate double `g55' = - cond($ML_y3==1,(`DM'*`deta_dsig'^2-2*`M'*`deta_dsig'/`sig1'+1/`sig1'^2-3*`eps1'^2/`sig1'^4)*`sig1',0)-`g5'; generate double `g57' = - cond($ML_y3==1,`DM'*`deta_drho'*`deta_dsig'+`M'*`deta_dsig'*(1/`rho1'+`rho1'*`rr1'^2),0); generate double `g44' = - cond($ML_y3==0,(`DM'*`deta_dsig'^2-2*`M'*`deta_dsig'/`sig0'+1/`sig0'^2-3*`eps0'^2/`sig0'^4)*`sig0',0)-`g4'; generate double `g46' = - cond($ML_y3==0,`DM'*`deta_drho'*`deta_dsig'+`M'*`deta_dsig'*(1/`rho0'+`rho0'*`rr0'^2),0); generate double `g66' = - cond($ML_y3==0,(`DM'*`deta_drho'^2+`M'*`d2eta_drho2')*`drho0_dtheta0'^2,0)-`g6'*`d2rho0_dtheta02'; generate double `g77' = - cond($ML_y3==1 ,(`DM'*`deta_drho'^2+`M'*`d2eta_drho2')*`drho1_dtheta1'^2,0)-`g7'*`d2rho1_dtheta12'; mlmatsum `lnf' `d11'=-(`DM'*`deta_dx'^2-1/`sig0'^2) if ($ML_y3==0), eq(1); mlmatsum `lnf' `d12'= 0, eq(1,2); mlmatsum `lnf' `d13'=-`DM'*`deta_dx'*`rr0' if ($ML_y3==0), eq(1,3); mlmatsum `lnf' `d14'=-((`DM'*`deta_dsig'-`M'/`sig0')*`deta_dx'-2*`eps0'/`sig0'^3)*`sig0' if ($ML_y3==0), eq(1,4); mlmatsum `lnf' `d15'= 0, eq(1,5); mlmatsum `lnf' `d16'=-(`DM'*`deta_drho'+`M'*(1/`rho0'+`rho0'*`rr0'^2))*`deta_dx'*`drho0_dtheta0' if ($ML_y3==0), 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'*`sig1', eq(2,5); mlmatsum `lnf' `d26'= 0, eq(2,6); mlmatsum `lnf' `d27'=`g27'*`drho1_dtheta1', eq(2,7); mlmatsum `lnf' `d33'=`g33', eq(3); mlmatsum `lnf' `d34'=`g34'*`sig0', eq(3,4); mlmatsum `lnf' `d35'=`g35'*`sig1', eq(3,5); mlmatsum `lnf' `d36'=`g36'*`drho0_dtheta0', eq(3,6); mlmatsum `lnf' `d37'=`g37'*`drho1_dtheta1', eq(3,7); mlmatsum `lnf' `d44'=`g44'*`sig0', eq(4); mlmatsum `lnf' `d45'= 0, eq(4,5); mlmatsum `lnf' `d46'=`g46'*`sig0'*`drho0_dtheta0', eq(4,6); mlmatsum `lnf' `d47'=0, eq(4,7); mlmatsum `lnf' `d55'=`g55'*`sig1', eq(5); mlmatsum `lnf' `d56'=0, eq(5,6); mlmatsum `lnf' `d57'=`g57'*`sig1'*`drho1_dtheta1', eq(5,7); mlmatsum `lnf' `d66'=`g66', eq(6); mlmatsum `lnf' `d67'= 0, eq(6,7); mlmatsum `lnf' `d77'=`g77', eq(7); capture matrix drop `negH'; forvalues k = 1/7 {; tempname H; forvalues l = 1/7 {; if (`l'>`k') matrix `H' =nullmat(`H'),`d`k'`l''; else matrix `H' =nullmat(`H'),`d`l'`k'''; }; // l matrix `negH' = nullmat(`negH') \ -`H'; }; // k }; // end quietly end;