Dear Statalisters,
I would like to calculate the average marginal effects from a recursive bivariate probit. Austin Nichols suggests three methods in his presentation "causal inference for binary regression". The first method computes slightly different AME compared to the other two procedures.

Code:
biprobit (y=x R) (R=x A)
margins, dydx(R) predict(pmarg1) force
loc ATEm=el(r(b),1,1)
predict double xb2, xb2
preserve
ren R TR
g R=0
predict double p0, pmarg1
predict double xb0, xb1
replace R=1
predict double p1, pmarg1
predict double xb1, xb1
g double dp=p1-p0
su dp, mean
loc ATE1=r(mean)
su dp if TR==1, mean
loc TOT1=r(mean)
loc r=e(rho)
gen double pdx=(binormal(xb1,xb2,`r')-binormal(xb0,xb2,`r'))/normal(xb2) if TR==1
su pdx, mean
loc TOT2=r(mean)
qui replace pdx=normal(xb1)-normal(xb0)
su pdx, mean
loc ATE2=r(mean)


A previous thread on this topic suggested to use the following method to calculate the AME, which computes the predicted conditional probabilities of success using the bivariate predicted probabilities and the univariate predicted marginal probability. The results in this case, however, are completely different from the ones obtained using one of the three methods suggested above.

Code:
biprobit (y=x R) (R=x A)
gen wasr=R
replace R=1
predict p1a if e(sample), p11
predict p1b if e(sample), p10
predict p1c if e(sample), p01
predict p1d if e(sample), p00
gen p1=p1a/(p1a+p1c)
replace R=0
predict p0a if e(sample), p11
predict p0b if e(sample), p10
predict p0c if e(sample), p01
predict p0d if e(sample), p00
gen p0=p0b/(p0b+p0d)
replace R=wasr
gen dp=p1-p0
sum dp
Isn't this second piece of code equivalent to ATE1 above? Any suggestion on what is the best way to proceed?

Thank you in advance