Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • How to do a logit random effects in mata?

    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
    Click image for larger version

Name:	1.jpg
Views:	1
Size:	35.2 KB
ID:	1494254


    I want to do it by simulation based method,and the code is as the following:
    ----------------------- copy starting from the next line -----------------------
    Code:
    
    
    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


Working...
X