Marginal effects of Hurdle Models with endogenous variable

[Yoann Morin]

Hello,

I'm currently estimating a Hurdle model (with only one hurdle) and I'm having issues computing joint marginal effects for both equations.

Here is the code i'm using to compute joint and independant marginal effect without IV:
(I used smoke dataset as an example for reproductibility)

ssc install twopm
use https://www.stata.com/data/jwooldridge/eacsap/smoke.dta

* Probit and its marginal effects
probit cigs educ restaurn lincome lcigpric, robust
margins, dydx(*)

* Poisson and its marginal effects
glm cigs educ restaurn lincome lcigpric if cigs>0, family(poisson) link(log) robust
margins, dydx(*)

* Joint marginal effects
twopm cigs educ restaurn lincome lcigpric , firstpart(probit) secondpart(glm, family(poisson) link(log)) vce(robust)
margins, dydx(*)

Then I consider for this example "lcigpric" is endogenous and instrumented with "age" and "agesq"

* IV Probit and its marginal effects
ivprobit cigs educ restaurn lincome (lcigpric=age agesq), vce(robust)
margins, dydx(*)

* IV Poisson and its marginal effects
ivpoisson gmm cigs educ restaurn lincome (lcigpric=age agesq) if cigs>0, vce(robust)
margins, dydx(*)

My issue is to find a way to compute joint marginal effect when one variable is endogenous.

I've seen two post on the stata blog dealing with the computation of marginal effects for multiple equation models, however the first method require the use of maximum likelihood estimates (which is not the case with the IV estimates I'm using), or to use GSEM but it does not allow ivprobit or ivpoisson either and doesn't allow for simple hurdle (because the second equation does not have the same number of observations as the first one).

Anyone know a method to overcome this issue please?

Best regards,
YM.

[Yoann Morin]

To make things more simple, I first try to reproduce twopm average marginal effects (AMEs)
I was able to reproduce AMEs for the probit and poisson with no issue but I'm struggling with twopm joint marginal effects.

***** To begin compute AMEs using twopm for the variable "educ" and save them in a new variable
twopm cigs educ restaurn lincome lcigpric, firstpart(probit) secondpart(glm, family(poisson) link(log)) vce(robust)
margins, dydx(educ) generate(twopm_m_educ)


*****Then compute AMEs manually

*Estimate and save probit coefficients
probit cigs educ restaurn lincome lcigpric, robust
predict xb_probit if e(sample), pr
matrix beta_prob = e(b)

*Estimate and save poisson coefficients
glm cigs educ restaurn lincome lcigpric if cigs>0, family(poisson) link(log) robust
predict xb_poisson, xb 
matrix beta_pois = e(b)

* Compute AMEs
gen manual_m_educ = normalden(invnorm(xb_probit)) * beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:educ")]*exp(xb_poisson) +
normalden(xb_probit)*exp(xb_poisson) * beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:educ")]

***** Compare results
* Summarize
summarize twopm_m_educ manual_m_educ 
Variable
Obs
Mean
Std. Dev.
Min
Max





twopm_m_ed~1
807
-.252511
.0476096
-.3538936
-.1237544

manual_m_e~c
807
-.2631823
.0515923
-.3472176
-.072222




* Compute the absolute difference between twopm and manual methods, then summarize and sum them
gen abs_dif = abs(twopm_m_educ - manual_m_educ)

summarize abs_dif 
Variable
Obs
Mean
Std. Dev.
Min
Max





abs_dif
807
.0493039
.0319773
7.69e-06
.1782453


display r(sum)
39.788252

The difference between both results is quite huge, while for the probit and poisson cases I had the exact same values.
Anyone has a idea of a mistake I could have made please?

[Yoann Morin]

Found my mistake, the right formula to compute joint margins is:

gen margin_educ = normalden(invnorm(xb_probit))*beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:educ")]*exp(xb_poisson) +
xb_probit*exp(xb_poisson)*beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:educ")]

Now I'll try to compute the variance using the delta method and I'll keep updating this post in case someone else need this in the future.

[Yoann Morin]

I'm having issues computing SEs using the delta method, so for now I've made a boostrap code to compute them instead.

Because I was having convergence issues, I've added to the code some parts to remove replications where at least one of the two estimates did not converge.
I've limited the maximum number of iterations to 100.

capture program drop margins_iv
program define margins_iv, rclass
preserve

ivprobit cigs educ restaurn lincome (lcigpric=age agesq), vce(robust) iterate(100)
predict xb_probit if e(sample), pr
matrix beta_prob = e(b)
scalar conv_probit = e(converged)

ivpoisson gmm cigs educ restaurn lincome (lcigpric=age agesq) if cigs>0, vce(robust) iterate(100)
predict xb_poisson, xb
matrix beta_pois = e(b)
scalar conv_poisson = e(converged)

scalar conv_final = conv_poisson*conv_probit

gen margin_educ = normalden(invnorm(xb_probit))*beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:educ")]*exp(xb_poisson) + xb_probit*exp(xb_poisson)*beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:educ")]
summarize margin_educ
return scalar margin_educ_manual = r(mean)

gen margin_restaurn = normalden(invnorm(xb_probit))*beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:restaurn")]*exp(xb_poisson) + xb_probit*exp(xb_poisson)*beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:restaurn")]
summarize margin_restaurn
return scalar margin_restaurn_manual = r(mean)

gen margin_lincome = normalden(invnorm(xb_probit))*beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:lincome")]*exp(xb_poisson) + xb_probit*exp(xb_poisson)*beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:lincome")]
summarize margin_lincome
return scalar margin_lincome_manual = r(mean)

gen margin_lcigpric = normalden(invnorm(xb_probit))*beta_prob[rownumb(beta_prob,"y1"),colnumb(beta_prob,"cigs:lcigpric")]*exp(xb_poisson) + xb_probit*exp(xb_poisson)*beta_pois[rownumb(beta_pois,"y1"),colnumb(beta_pois,"cigs:lcigpric")]
summarize margin_lcigpric
return scalar margin_lcigpric_manual = r(mean)

restore
    exit
end

bootstrap margin_educ=r(margin_educ_manual) margin_restaurn=r(margin_restaurn_manual) margin_lincome=r(margin_lincome_manual) margin_lcigpric=r(margin_lcigpric_manual), reps(100) seed(123) reject(conv_final!=1) : margins_iv