how to do a maximum simulated likelihood estimation?

[Bing Zhang]

Dear statlist:
I want to do a MSLE for the following equation systems:
the choice function: D_it = 1 {z_itγ + ηa_i + v_it}
the wage function:y_it = x_itβ + λa_i + u_it
the a_i is an unobservable heterogeneity,the choice equation is estimated by probit and the wage equation is estimated by finite mixture model which is consisted of two ols functions.
then we will have
f (y_it, D_it = 1|a_i) = f (y_it|a_i) × Pr (D_it =1|a_i)
= [p/σ₁ (x_it)ϕ ((y_it−x_itβ₁ -λa_i ) /σ₁ (x_it))+ (1-p)/σ₂ (x_it)ϕ ((y_it−x_itβ₁ -λa_i )/σ₂(x_it))]× Φ (z_itγ + ηa_i)

so,we can have a unconditional likelihood function like that:
ln L =∑ln [∫ (∑D_itf (y_it, D_it =1|a_i) + (1- D_it)Pr (D_{it =} 0|a_i))f(a_i) da_i]

After reading Maximum simulated likelihood estimation of a negative binomial regression model with multinomial endogenous treatment(The Stata Journal (2006)
6, Number 2, pp. 246–255),I have coded following：

----------------------- copy starting from the next line -----------------------

mata:
function mmprobit_lf(
        string scalar lnL, real matrix rmat, string scalar neq
    , string scalar neqall, string scalar nobs, string scalar sim
    , string scalar xb11, string scalar xb12, string scalar xb13
    , string scalar xb14, string scalar xb15, string scalar xb16
    , string scalar xb17, string scalar xb18, string scalar xb19
    , string scalar xb2
    , string scalar eta1, string scalar eta2, string scalar eta3
    , string scalar eta4, string scalar eta5, string scalar eta6
    , string scalar eta7, string scalar eta8, string scalar eta9
    , string scalar lam1, string scalar lam2, string scalar lam3
    , string scalar lam4, string scalar lam5, string scalar lam6
    , string scalar lam7, string scalar lam8, string scalar lam9
    , string scalar g1, string scalar g2, string scalar g3
    , string scalar g4, string scalar g5, string scalar g6
    , string scalar g7, string scalar g8, string scalar g9
    , string scalar g10, string scalar g11, string scalar g12
    , string scalar g13, string scalar g14, string scalar g15
    , string scalar g16, string scalar g17, string scalar g18
    , string scalar g19, string scalar g20
    , string scalar H)

    {
    neq=st_numscalar(neq)
    nobs=st_numscalar(nobs)
    sim=st_numscalar(sim)

    yl = *(findexternal("_mtreatfmm_yl"))
    y2  = *(findexternal("_mtreatfmm_y2"))
    
    xbl = *(findexternal("_mtreatnb_xmml"))
    xbl = (xmml, J(nobs,1,1))
    xb2 = *(findexternal("_mtreatnb_xb2"))
    xb2 = (xnb, J(nobs,1,1))
    xb3 = *(findexternal("_mtreatnb_xb2"))
    xb3 = (xnb, J(nobs,1,1))    
    st_view(xbl, ., (xb11,xb12,xb13,xb14,xb15,xb16,xb17,xb18,xb19))
    st_view(xb2, ., xb2)
    st_view(lam, ., (lam1,lam2,lam3,lam4,lam5,lam6,lam7,lam8,lam9))
    
    
    xb1mat = J(nobs,sim,0)
    for (j=1; j<=sim; j++) {
        xb1mat[,j] = xb1
    }
    for (j=1; j<=neq; j++) {
        xb1mat = xb1mat + eta[,j]:*rmat[,((j-1)*sim+1)..(j*sim)]
    }

prob = exp(lnnormal(xb1mat))+(1-y1)ln(1-normal(xb1mat))
prob1=    exp(lnnormal(xb1))+(1-y1)ln(1-normal(xb1mat))
    
    xb2mat = J(nobs,sim,0)
    for (j=1; j<=sim; j++) {
    xb2mat[,j] = xb2
    }
    for (j=1; j<=neq; j++) {
        xb2mat = xb2mat + lam[,j]:*rmat[,((j-1)*sim+1)..(j*sim)]
    }

    xb3mat = J(nobs,sim,0)
    for (j=1; j<=sim; j++) {
    xb3mat[,j] = xb2
    }
    for (j=1; j<=neq; j++) {
    xb3mat = xb3mat + lam[,j]:*rmat[,((j-1)*sim+1)..(j*sim)]
    }

fmm = p*normalden(y2,xb2mat,sigma1)+(1-p)*normalden(y2,xb3mat,sigma1)
p = exp(xb2mat)/(exp(xb2mat)+exp(xb3mat))

    fllk=y1*(p*normalden(y2,xb2mat,sigma2)+(1-p)*normalden(y2,xb3mat,sigma3))+(1-y1)*(1-normal(y1,xb1mat,sigma1))
    
    L = (rowsum(fllk))/sim
    L = rowmax((L , J(nobs,1,smallestdouble())))
    vlnL = ln(L)
end

------------------ copy up to and including the previous line ------------------

but I don't know how to do the following work.Mata code provided by professor Partha Deb in the refered paper has no technical notes.I wonder if you can tell me more details about the following work,or another mata code with notes.
Thank you for your patience and I am looking forward to your feedback