How to do a logit random effects in mata?

[Bing Zhang]

Dear statalists:
I want to do a logit random effects in mata.I have refered the xtlogit ,re in Cameron's book(Microeconometrics Using Stata: Revised Edition:chapter 18.4.6),and wrote the code in mata,but it always went wrong.I wonder if you can give me some suggestions and corrections to my code.
the likelihood function is the

I want to do it by simulation based method,and the code is as the following:
----------------------- copy starting from the next line -----------------------

webuse womenwk, clear
gen dwage=wage!=.
qui sum dwage
di r(sum)
program xtlogit_lf
    version 15.0
    args todo b lnf g negH 

    tempname lnL gi
    tempname xb11 xb12 xb13 xb14 xb15 xb16 xb17 xb18 xb19

    qui gen double `lnL'=.
    mleval `xb11' = `b', eq(1)    

    forvalues i=2/9 {
        qui gen double `xb1`i''=.
        qui gen double `g`i''=.
    }
    
    qui sum `touse'
    scalar nobs = r(sum)
    scalar sim = 50
    
    mata: _mtreatnb_rmat = `scale'*invnormal(halton(`nobs'*`sim',1))
    mata: xtlogit_lf("`lnL'",_mtreatnb_rmat,"nobs","sim" ///
        ,"`xb11'","`xb12'","`xb13'","`xb14'","`xb15'","`xb16'","`xb17'","`xb18'","`xb19'" ///
        ,"`g1'","`g2'","`g3'","`g4'","`g5'","`g6'","`g7'","`g8'","`g9'","H")

    mlsum `lnf' = `lnL'

    local k = colsof(`b')
    local c 1
    matrix `g' = J(1,`k',0)
    

    mlvecsum `lnf' `g1' = `g1', eq(1)
    matrix `g'[1,`c'] = `gi'
    local c = `c' + colsof(`gi')

    mat `negH' = -H

end

ml model d2 xtlogit_lf (xb11:dwage=married children educ age)
ml maximize


mata:
function xtlogit_lf(
        string scalar lnL, real matrix rmat
    , 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 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 H)

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

    y = *(findexternal("_xtlogit_y"))
    x = *(findexternal("_xtlogit_x"))
    x = (x, J(nobs,1,1))

    exb = J(nobs,sim,.)
    pmml = J(nobs,sim*neq,.)
    den = J(nobs,sim,1)
    exb[,1..sim] = exp(xb[,j]:+rmat[,1..sim])

    xtlogit = ((exb[,1..sim]/(1+exb[,1..sim)]))^y):*((1/(1+exb[,1..sim])^(1-y))

    L = (rowsum(xtlogit))/sim
    L=rowmax((L , J(nobs,1,smallestdouble())))
    vlnL = ln(L)

    gxtlogit = J(nobs,sim,.)
    vg = J(nobs,9,.)
    gxtlogit[,1..sim] = y:-(exb[,1..sim]:/(1:+exb[,1..sim]))
    vg[,1] = (1:/L):*rowsum(xtlogit:*gxtlogit[,1..sim])/sim

    nparmsmml = cols(x)
    H = J(nparmsmml,nparmsmml,.)
    hxtlogit = -(exb[,1..sim]:/(1:+exb[,1..sim]))^2
    h = ((-vg[,j]:*vg[,j] + (1:/L):*rowsum(xtlogit ///
                    :*gxtlogit[,1..sim]:*gxtlogit[,1..sim])/sim ///
                    :+ (1:/L):*rowsum(xtlogit:*hxtlogit)/sim) :* x)'*x
            H[1..rows(h),1..cols(h)] = h
            
    st_store(., lnL, vlnL)
    st_store(., (g1,g2,g3,g4,g5,g6,g7,g8,g9), vg)
    st_matrix("H", H)
    }

    mata mosave mmlogit_lf(), dir(`mydir') replace

end



end

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



Thank you for your patience and I'm looking forward to your feendback.
Ice Zhang