/**********************************************************/ /* Heckman probit d2 evoluator */ /* */ /* by Jerzy Mycielski */ /* Warsaw Univeristy */ /* 16.07.2004 */ /**********************************************************/ program define heckp_d2 version 8 args todo b lnf g negH tempvar Lj gamma1 gamma2 tempname athrho mleval `gamma1' = `b', eq(1) mleval `gamma2' = `b', eq(2) mleval `athrho' = `b', eq(3) scalar /*assuming that rho is scalar, can be generalized*/ tempname rho tempvar `Lj' scalar `rho' = tanh(`athrho') /*if $ML_y1!=0 assuming that rho is scalar, can be generalized*/ /**************Logrithm of likelihood function****************/ local prob_0 normprob(-`gamma2') local prob_11 binorm(`gamma2',`gamma1',`rho') local prob_10 (1-`prob_0'-`prob_11') quiet gen double `Lj' = `prob_0' if $ML_y2==0 quiet replace `Lj' = `prob_10' if $ML_y2!=0 & $ML_y1==0 quiet replace `Lj' = `prob_11' if $ML_y2!=0 & $ML_y1!=0 mlsum `lnf' = ln(`Lj') if (`todo'==0 | `lnf'==.) exit /**********************First order derivatives of binominal normal distribution*******************/ tempvar Fg1 Fg2 Fr quietly gen double `Fg1' = normprob((`gamma2'-`rho'*`gamma1')/sqrt(1-(`rho')^2))*normden(`gamma1') if $ML_y2!=0 quietly gen double `Fg2' = normprob((`gamma1'-`rho'*`gamma2')/sqrt(1-(`rho')^2))*normden(`gamma2') if $ML_y2!=0 quietly gen double `Fr' = normden((`gamma2'-`rho'*`gamma1')/sqrt(1-(`rho')^2))*normden(`gamma1')/sqrt(1-(`rho')^2) if $ML_y2!=0 /*First order derivatives with respect to g1, g2 i rho*/ tempvar lg1 lg2 lr /*for y2=0*/ quietly gen double `lg1' = 0 if $ML_y2==0 quietly gen double `lg2' = -normden(`gamma2')/`Lj' if $ML_y2==0 quietly gen double `lr' = 0 if $ML_y2==0 /*for y2=1 y1=0*/ quietly replace `lg1' = -`Fg1'/`Lj' if $ML_y2!=0 & $ML_y1==0 quietly replace `lg2' = (normden(`gamma2')-`Fg2')/`Lj' if $ML_y2!=0 & $ML_y1==0 quietly replace `lr' = -`Fr' /`Lj' if $ML_y2!=0 & $ML_y1==0 /*for y2=1 y1=1*/ quietly replace `lg1' = `Fg1'/`Lj' if $ML_y2!=0 & $ML_y1!=0 quietly replace `lg2' = `Fg2'/`Lj' if $ML_y2!=0 & $ML_y1!=0 quietly replace `lr' = `Fr' /`Lj' if $ML_y2!=0 & $ML_y1!=0 /*First order derivatives with respect with respect to g1,g2, tan(rho)*/ tempname dg1 dg2 dr mlvecsum `lnf' `dg1' = `lg1', eq(1) mlvecsum `lnf' `dg2' = `lg2', eq(2) mlvecsum `lnf' `dr' = `lr'*(1-(`rho')^2), eq(3) /*tanh(rho)!*/ matrix `g' = (`dg1',`dg2',`dr') if (`todo'==1 | `lnf'==.) exit /**********************Second order derivatives of binominal normal distribution***********************/ local Fg11 (-`rho'*`Fr'-`gamma1'*`Fg1') local Fg22 (-`rho'*`Fr'-`gamma2'*`Fg2') local Fg12 `Fr' local Fg1r -(`gamma1'-`rho'*`gamma2')/(1-(`rho')^2)*`Fr' local Fg2r -(`gamma2'-`rho'*`gamma1')/(1-(`rho')^2)*`Fr' local Frr 1/(1-(`rho')^2)*(1+(`gamma2'-`rho'*`gamma1')*(`gamma1'-`rho'*`gamma2')/(1-(`rho')^2))*`Fr' /*Second order derivatives with respect to g1, g2 i rho*/ tempvar lg11 lg22 lg12 lg1r lg2r lrr /*for y1=0*/ quietly gen double `lg11' = 0 if $ML_y2==0 quietly gen double `lg22' = -`gamma2'*`lg2'-`lg2'^2 if $ML_y2==0 quietly gen double `lg12' = 0 if $ML_y2==0 quietly gen double `lg1r' = 0 if $ML_y2==0 quietly gen double `lg2r' = 0 if $ML_y2==0 quietly gen double `lrr' = 0 if $ML_y2==0 /*for y2=1 y1=0*/ quietly replace `lg11' = -`Fg11'/`Lj'-`lg1'^2 if $ML_y2!=0 & $ML_y1==0 quietly replace `lg22' = (-`gamma2'*normden(`gamma2')-`Fg22')/`Lj'-`lg2'^2 if $ML_y2!=0 & $ML_y1==0 quietly replace `lg12' = -`Fg12'/`Lj'-`lg1'*`lg2' if $ML_y2!=0 & $ML_y1==0 quietly replace `lg1r' = -`Fg1r'/`Lj'-`lg1'*`lr' if $ML_y2!=0 & $ML_y1==0 quietly replace `lg2r' = -`Fg2r'/`Lj'-`lg2'*`lr' if $ML_y2!=0 & $ML_y1==0 quietly replace `lrr' = -`Frr' /`Lj'-`lr'^2 if $ML_y2!=0 & $ML_y1==0 /*for y2=1 y2=1*/ quietly replace `lg11' = `Fg11'/`Lj'-`lg1'^2 if $ML_y2!=0 & $ML_y1!=0 quietly replace `lg22' = `Fg22'/`Lj'-`lg2'^2 if $ML_y2!=0 & $ML_y1!=0 quietly replace `lg12' = `Fg12'/`Lj'-`lg1'*`lg2' if $ML_y2!=0 & $ML_y1!=0 quietly replace `lg1r' = `Fg1r'/`Lj'-`lg1'*`lr' if $ML_y2!=0 & $ML_y1!=0 quietly replace `lg2r' = `Fg2r'/`Lj'-`lg2'*`lr' if $ML_y2!=0 & $ML_y1!=0 quietly replace `lrr' = `Frr' /`Lj'-`lr'^2 if $ML_y2!=0 & $ML_y1!=0 /*Matrix of second order derivatives with respect with respect to g1,g2, tan(rho)*/ tempname dg11 dg12 dg1r dg22 dg2r drr mlmatsum `lnf' `dg11' = `lg11', eq(1) mlmatsum `lnf' `dg12' = `lg12', eq(1,2) mlmatsum `lnf' `dg1r' = `lg1r'*(1-`rho'^2), eq(1,3) /*tanh(rho)! */ mlmatsum `lnf' `dg22' = `lg22', eq(2) mlmatsum `lnf' `dg2r' = `lg2r'*(1-`rho'^2), eq(2,3) /*tanh(rho)! */ mlmatsum `lnf' `drr' = `lrr'*(1-`rho'^2)^2-2*`lr'*`rho'*(1-`rho'^2), eq(3) /*tanh(rho)! */ matrix `negH' = -(`dg11' ,`dg12' ,`dg1r' \/* */ `dg12'',`dg22' ,`dg2r' \/* */ `dg1r'',`dg2r'',`drr' ) end