Dear Statalist,

I am trying to implement a routine based on the following recent article in the Stata Journal (i've attached it for your convenience):

Terza, J. V. (2017). Two-stage residual inclusion estimation: A practitioners guide to Stata implementation. The Stata Journal, 17(4), 916-938.


I am following the codes on p15, reproduced here:

Step c:
/*step c*/ probit Xe Xo Wplus Step d:
/*step d*/ predict phiWalpha, p gen Xuhat=Y-phiWalpha Step e:
/*step e*/ mata: alphahat=st_matrix("e(b)") ́ mata: Valphahat=st_matrix("e(V)") Step f:
/*step f*/ glm Y Xe Xo Xuhat,family(gaussian) link(probit) vce(robust) Step g:
/*step g*/ mata: betahat=st_matrix("e(b)") ́ mata: Vbetahat=st_matrix("e(V)") mata: Bu=betahat[3] Step h:
/*step h*/ putmata Y Xe Xo Wplus Xuhat mata: X=Xe, Xo, Xuhat, J(rows(Xo),1,1) mata: W=Xo, Wplus, J(rows(Xo),1,1) Step i:
/*step i*/ mata: gradbeta=normalden(X*betahat):*X mata: gradalpha=-Bu:*normalden(X*betahat):*/*
*/normalden(W*alphahat):*W Step j:
/*step j*/ mata: B1 = gradbeta ́*gradbeta mata: B2 = gradbeta ́*gradalpha Step k:
/*step k*/ mata: AVARBeta=invsym(B1)*B2*Valphahat*B2 ́invsym(B1)/* */+ Vbetahat


Step l:
/*step l*/ mata: ACSE = sqrt(diagonal(AVARBeta))
Step m:
/*step m*/ mata: ACtstats=betahat:/ACSE
The goal of these steps is to be able to obtain the asymptotically correct standard errors (ACSE) when implementing the two-stage residual inclusion. I seem to get stuck on step J onwards. The problem is that when I checked the contents of B1 and B2, it seemed to return a full matrix of missing values. As a result, other variables (e.g. AVARBeta, ACSE, ACtstats) from step J onwards (i.e. steps K, L, M) all return missing values.

Not sure if it helps, but when I looked into gradbeta and gradalpha they seemed to have a mix of both missing and non-missing values.

I was wondering if you have any thoughts on this. Any help would be much appreciated.

Thank you very much.

Best,
GuiDeng