what's going wrong with my msm for logit fixed effects

[Bing Zhang]

Dear statalists:
I want to use maximum simulation method to estimate a logit with fixed effects,but it always going wrong and shows "too few variables specified".I wonder if you can tell me the solution.I will be very greatful if other mistakes in this program are pointed out.
Thank you for your patience and I'm looking forward to your reply.
Ice Zhang

the code is here:

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

* Example generated by -dataex-. To install: ssc install dataex

webuse union
bys idcode:gen nt=_N
keep if nt==8

sort year
bys idcode (year):gen cid=_n
global draws "20"
mdraws if cid==1,neq(1) dr($draws) prefix(z) burn(15)

forvalues i=1/20 {
    bys idcode (year): replace z1_`i' = z1_`i'[_n-1] if z1_`i'==.
    }

forvalues r=1/$draws{
gen random_1`r'=invnormal(z1_`r')
 }
 
 program define xtlogit_lf
args todo b lnf
tempvar xb random1 lj p1 p2 sum lnpi 
tempname lnsig sigma 
mleval `xb' = `b', eq(1)

qui {

    gen double `random1' = 0
    gen double `lnpi'=0
    gen double `sum'=0
    gen double `L1'=0
    gen double `L2'=0
    by idcode: gen byte `last'=(_n==_N)
    gen double `p1'=0
    gen double `p2'=0
}

forvalues r=1/$draws{
qui {
    replace `random1' = random_1`r'
    replace `p1'= 1/(1 + exp(`xb'+`random1')
    replace `p2'= exp(`x'+`random1')*`p1'
    
    replace `lnpi'=ln(`p1')*(1-y)+ln(`p2')*y
    by idcode: replace `sum'=sum(`lnpi')
    by idcode: replace `L1' =exp(`sum'[_N]) if _n==_N
    by idcode: replace `L2'=`L2'+`L1' if _n==_N
}
}
qui gen `lj'=cond(!`last',0, ln(`L2'/`$draws'))
qui mlsum `lnf'=`lj'

end

ml model d0 xtlogit_lf ( union=age grade i.not_smsa south##c.year)
ml check


end

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