*==============================================================================
> ===
* Check if the ideas presented in
* http://blog.stata.com/2011/08/22/use-poisson-rather-than-regress-tell-a-frien
> d/
* also apply to -propcnsreg-
* not working as well as I would have liked/hoped, especially the standard erro
> rs
*
* also requires -simsum- and -simpplot-, both can be installed using -ssc-
*==============================================================================
> ===


clear all
program define sim, rclass
        drop _all
        set obs 5000
        gen z1 = rnormal()
        gen z2 = rnormal()

        gen lat = .25*z1 + .5*z2
        gen x1 = rnormal()
        gen x2 = rnormal()
        gen lam = 1 + .5*x1
        gen y = exp(.25*x1 + .25*x2 + lam*lat + .2*rnormal())

        propcnsreg y x1 x2 , constr(z1 z2) lambda(x1) lcons poisson robust
        return scalar w     = e(w_p)
        return scalar bz1   =   [constrained]_b[z1]
        return scalar bz2   =   [constrained]_b[z2]
        return scalar blx1  =        [lambda]_b[x1]
        return scalar bx1   = [unconstrained]_b[x1]
        return scalar sez1  =   [constrained]_se[z1]
        return scalar sez2  =   [constrained]_se[z2]
        return scalar selx1 =        [lambda]_se[x1]
        return scalar sex1  = [unconstrained]_se[x1]
end

set seed 123456

simulate w=r(w)                  ///
bz1  = r(bz1)  sez1  = r(sez1)   ///
bz2  = r(bz2)  sez2  = r(sez2)   ///
blx1 = r(blx1) selx1 = r(selx1)  ///
bx1  = r(bx1)  sex1  = r(sex1),  ///
reps(5000):sim

simsum bz1 , se(sez1)  true(.25) mcse
simsum bz2 , se(sez2)  true(.50) mcse
simsum blx1, se(selx1) true(.50) mcse 
simsum bx1 , se(sex1)  true(.25) mcse

gen pz1  = 2*normal(-abs(( bz1-.25)/sez1 ))
gen pz2  = 2*normal(-abs(( bz2-.50)/sez2 ))
gen plx1 = 2*normal(-abs((blx1-.50)/selx1))
gen px1  = 2*normal(-abs(( bx1-.25)/sex1 ))

simpplot w pz1 pz2 plx1 px1, scheme(s2color) name(pois_coef, replace)

exit