******************************************************************************** *! "cub14s", v.29, Cerulli, 14dic2021 (cub with shelter) ******************************************************************************** program cub14s , eclass version 14.1 args todo b lnf tempvar theta1 theta2 theta3 mleval `theta1' = `b', eq(1) mleval `theta2' = `b', eq(2) mleval `theta3' = `b', eq(3) // new for shelter local y "\$ML_y1" // this is just for readability local m=e(M) tempvar p M R S D delta * Calculate p quietly generate double `p' = 1/(1+exp(-`theta1')) * Calculate M local c = exp(lnfactorial(`m'-1)) tempname cmb mat `cmb' = J(`m',1,.) levelsof `y' , local(LEV_Y) ******************************************************************************** forvalues i=1/`m'{ foreach j of local LEV_Y { if `j'==`i'{ scalar d = (exp(lnfactorial(`j'-1))*exp(lnfactorial(`m'-`j'))) mat `cmb'[`i',1] = `c'/d } } } ******************************************************************************** qui gen double `M' = `cmb'[`y',1] ******************************************************************************** * Calculate R quietly generate double `R' = ((exp(-`theta2'))^(`y'-1))/((1+exp(-`theta2'))^(`m'-1)) * Calculate S quietly generate double `S' = 1/`m' * Calculate D quietly generate double `D'=(`y'==`e(SHELTER)') // new for shelter * Calculate delta quietly generate double `delta'= 1/(1+exp(-`theta3')) // new for shelter mlsum `lnf' = ln(`delta'*`D' + (1-`delta')*(`p'*(`M'*`R'-`S')+`S')) // new for shelter ereturn scalar M=`m' end ******************************************************************************** * END ********************************************************************************