*! Version 3.0 *! Author: Alfonso Miranda *! Date: 17.12.2009 program define petpoisson_ll version 11 args todo b lnf g NegH tempvar xb1 xb2 xb3 last vartype edv ssv ordv #delimit ; tempname lambda_1 lambda_2 sigma_u lambdaM Vu lf hvec; #delimit cr local Vs `Vu' local xb `xb1' `xb2' `xb3' local nex1_il =$S_snex1_il local nex2_il =$S_snex2_il local nex3_il =$S_snex3_il local nscores = 6 local lambda `lambda_1' `lambda_2' local sigma_u `sigma_u' local id "$S_id" local rep "$S_rep" local endogvs "$S_endogvs" local method "$S_method" local samp "$ML_samp" local vartype "$S_vartype" local resp "$S_resp" local k = 0 scalar `lf' = . foreach v in `xb' { mleval `v' = `b', eq(`++k') } foreach v in `lambda' { mleval `v' = `b', scalar eq(`++k') } foreach v in `sigma_u' { mleval `v' = `b', scalar eq(`++k') } qui sort $S_id $S_vartype qui gen `edv'=cond(`vartype'==1,1,0) qui gen `ssv'=cond(`vartype'==2,1,0) qui gen `ordv'=cond(`vartype'==3,1,0) qui sort $S_id qui by $S_id : gen `last' = cond(_n==_N,1,0) qui { by $S_id: replace `xb1'=. if `edv'==0 & `samp' by $S_id: replace `xb2'=. if `ssv'==0 & `samp' by $S_id: replace `xb3'=. if `ordv'==0 & `samp' } /* ensure that variance of y is within admisible range */ if (`sigma_u' < -20) scalar `sigma_u' = -20 if (`sigma_u' > 20) scalar `sigma_u' = 20 /* generate some stuff for scores */ tempname c_sigma_u scalar `c_sigma_u' = exp(`sigma_u') /* ensure that SEs of y is positive */ scalar `sigma_u' = exp(`sigma_u') /* prepare input matrices */ mat `lambdaM' =(`lambda_1',`lambda_2') mat `Vs' = (`sigma_u') /* prepare variables for scores */ local D "" forval i=1/`nscores' { tempvar d`i' local D "`D' `d`i''" qui gen double `d`i'' = . } /* matrix hvec */ mat `hvec' = ($S_hv_1, $S_hv_2, $S_hv_3) /* calculate log-likelihood */ tempvar pr qui gen double `pr' = . mata: petpoisson_v2d1("xb","resp","endogvs","vartype","`lambdaM'","id","Vs",`method',`rep',"`hvec'","`pr'","`D'","`samp'") qui replace `pr' = log(`pr') if `samp' mlsum `lnf' = `pr' if `last' & `samp' if (`todo'==0 | `lnf'>=.) { matrix `g' = J(1,colsof(`b'),.) matrix `NegH' = J(colsof(`b'),colsof(`b'),.) exit } /* do final adjustment to D */ qui replace `d`nscores'' = `d`nscores''*`c_sigma_u' /* Calculate the scores at equation level and feed them to Stata */ tempvar g1 g2 g3 local j = `nscores' - 3 *local j = `nscores' - 2 local k_1 = `j' + 1 local k_2 = `k_1' + 1 local h = `k_2' + 1 tempvar g`k_1' g`k_2' g`h' mlvecsum `lnf' `g1' = `d1' if `edv'==1, eq(1) /* equations */ mlvecsum `lnf' `g2' = `d2' if `ssv'==1, eq(2) mlvecsum `lnf' `g3' = `d3' if `ordv'==1, eq(3) mlvecsum `lnf' `g`k_1'' = `d`k_1'' if `edv'==1, eq(`k_1') /* lambdas */ mlvecsum `lnf' `g`k_2'' = `d`k_2'' if `ssv'==1, eq(`k_2') mlvecsum `lnf' `g`h'' = `d`h'' if `last'==1, eq(`h') /* sigma_u */ /* deliver matrix of first derivatives */ matrix `g' = (`g1',`g2',`g3') matrix `g' = (`g',`g`k_1'',`g`k_2'',`g`h'') /* Obtain the scores at variable level to calculate OPG approx to Variance */ local scores "" /* calculate scores eqn 1 */ local scores_x_1 "" local i = 1 foreach var of global S_x_1 { tempvar scores_x_1_`i' by $S_id: gen double `scores_x_1_`i'' = sum(`d1'*`var') if `samp' local scores_x_1 "`scores_x_1' `scores_x_1_`i''" local i = `i' + 1 } if $S_sncedv == 1 { local i = `i' + 1 tempvar scores_x_1_`i' by $S_id: gen double `scores_x_1_`i'' = sum(`d1') if `samp' local scores_x_1 "`scores_x_1' `scores_x_1_`i''" } local scores "`scores' `scores_x_1'" /* calculate scores eqn 2 */ local scores_x_2 "" local i = 1 foreach var of global S_x_2 { tempvar scores_x_2_`i' by $S_id: gen double `scores_x_2_`i'' = sum(`d2'*`var') if `samp' local scores_x_2 "`scores_x_2' `scores_x_2_`i''" local i = `i' + 1 } if $S_sncssv==1 { local i = `i' + 1 tempvar scores_x_2_`i' by $S_id: gen double `scores_x_2_`i'' = sum(`d2') if `samp' local scores_x_2 "`scores_x_2' `scores_x_2_`i''" } local scores "`scores' `scores_x_2'" /* calculate scores eqn 3 */ local scores_x_3 "" local i = 1 foreach var of global S_x_3 { tempvar scores_x_3_`i' by $S_id: gen double `scores_x_3_`i'' = sum(`d3'*`var') if `samp' local scores_x_3 "`scores_x_3' `scores_x_3_`i''" local i = `i' + 1 } if $S_sncord==1 { local i = `i' + 1 tempvar scores_x_3_`i' by $S_id: gen double `scores_x_3_`i'' = sum(`d3') if `samp' local scores_x_3 "`scores_x_3' `scores_x_3_`i''" } local scores "`scores' `scores_x_3'" /* lambdas */ tempvar score_lambda_1 score_lambda_2 by $S_id: gen double `score_lambda_1' = sum(`d`k_1'') if `samp' by $S_id: gen double `score_lambda_2' = sum(`d`k_2'') if `samp' local scores "`scores' `score_lambda_1' `score_lambda_2'" /* sigma_y */ tempvar score_sigma_u by $S_id: gen double `score_sigma_u' = sum(`d`h'') if `samp' local scores "`scores' `score_sigma_u'" /* Keep the scores for future use */ local i = 1 tokenize $S_snames foreach var of local scores { qui replace `1' = `var' if `samp' mac shift } if (`todo'==1 | `lnf'>=.) exit /* Calculate OPG covariance matrix */ //qui matrix accum `NegH' = `scores' if `last' & `samp', noconst tempvar tus gen `tus' = (`last'==1 & `samp'==1) mata: OPG("scores","tus","`NegH'") end