try this one, not sure if it's what you want:
probit ys z
predict double zg, xb
predict double p
* PI, put your values hear separated by space
local pis ".2 .2 .2 .2 .2"
* RHO, put your values hear separated by space
local rhos ".2 .2 .2 .2 .2"
gen double f = 0
local J=5
forv j=1/`J' {
local rho : word `j' of `rhos'
local pi : word `j' of `pis'
reg lny`j' x
predict double xb`j'
local sig`j'=e(rmse)
gen double cf`j' = normalden((lny`j'-xb`j')/`sig`j'')/`sig`j''/p*normprob((zg+(`rho'/`sig`j'')*(lny`j'-xb`j'))/sqrt(1-(`rho')^2))
replace f=f+`pi'*cf`j'
}