Dear Statalist members, I would be very thankful for an advice. I am running a biprobit model, and use eststo/esttab commands to output marginal effects. I need to add the rho coefficient (ideally, with significe level) to the output, but so far I was not able to figure out how to do it. When I use the estadd scalar corr = e(rho), the resulting scalar is empty after I run the eststo: margins
Here is the code I have worked so far with:

biprobit y1 y2 female $controls_fainc i.date , vce(robust)
estadd scalar corr = e(rho)
eststo: margins, predict (pmarg1) dydx(female $controls_fainc) atmeans noatlegend post
sum y1 if e(sample)==1
estadd scalar mean= r(mean)
eststo: margins, predict (pmarg2) dydx(female $controls_fainc) atmeans noatlegend post
sum y2 if e(sample)==1
estadd scalar mean= r(mean)


esttab using "${results}\table1.tex", replace style(tex) keep(female $controls_fainc) order(female $controls_fainc) ///
nonumber star(* 0.10 ** 0.05 *** 0.01) varlabels(_cons Constant) ///
cells(b(star fmt(3)) se(par fmt(3))) label booktabs alignment(D{.}{.}{-1}) collabels(none) nodepvar ///
stats(mean N corr , layout("@" "@" "\multicolumn{2}{c}{@}") labels( "\text{Mean dependent variable}" "\text{N observations}" "\text{Rho}") fmt(3 0 3))


Would be very obliged for any tip! Thank you!
Olga