How to do simulation in Stata with a variable generated in Mata

[Cuicui]

Dear all,

I have an estimating function as y=b1+x2*b2+x3*b3+x4*b4+u. The sample size is N. u follows a multivariate distribution with a NxN correlation matrix generated from Mata.
My strategy is to first generate an NX1 vector of u from Mata and transfer it to Stata as a variable. The other variables x2, x3, x4 are generated in Stata using the simulate command.
I tried two methods, the first one is: I use the simulate command and put the Mata commands in between the Stata commands. The second one is: I first simulate the data set of x2 x3 x4 and merge u to the simulated data set.
But it seems they could not go through, can anyone help me?
Codes for first method is in attachment.
Thank you.

[Christophe Kolodziejczyk]

The end statement for the mata code is conflicting with the end statement of your Stata program. Try to put the mata code in curly brackets.
Another way is to wrap your Mata code in a function.
Numgroups should be between quotes when calling st_global and there were parentheses missing.

clear all
global numobs 400
global numgroups 100
global groupobs 4
global numsims "100"
set more off
*** Program gee ***
****************************
capture program drop gee
program define gee, rclass
    version 13
    drop _all
set obs $numgroups
gen g=_n
generate y=rnormal()
expand $groupobs
sort g
by g: gen l=_n

generate double one=1
generate double x2=1+rnormal(2)
generate double x3=0.2*x2-1.2*rnormal(1)
generate double x5=0.2*x2+0.2*x3+rnormal(0)
generate double x4=(x5>0)
scalar b1=1
scalar b2=1
scalar b3=1
scalar b4=1

*MATA TO CREATE CORRELATION MATRIX
mata {
G=strtoreal(st_global("numgroups")) // <== new
L=strtoreal(st_global("groupobs")) // <== new
N=G*L
rho=0.6
A=(J(sqrt(G),1,1))#((J(sqrt(G),1,1))#((1..sqrt(L))'#J(sqrt(L),1,1)))
B=2*(0::sqrt(G)-1)#(J(sqrt(G)*L,1,1))
I=A+B
U=J(G,1,1)#(J(sqrt(L),1,1)#(sqrt(L)::1))
W=J(sqrt(G),1,1)#((2*(sqrt(G)-1::0))#J(L,1,1))
J=U+W
I=J(N,1,1)'#I
J=J(N,1,1)'#J
D=sqrt((I-I'):^2+(J-J'):^2)
W=((1-rho):*I(N)+rho:/(D+I(N))) 
e=rnormal(N,1,0,1)
C=-0.5:+cholesky(W)*e
u=exp(C)
st_addvar("double", "u")
st_store(.,"u",u)
}

by g: replace y=b1+x2*b2+x3*b3+x4*b4+u
end

clear all
global numobs 400
global numgroups 100
global groupobs 4
global numsims "100"
set more off
*** Program gee ***
****************************
capture program drop gee
program define gee, rclass
    version 12
    drop _all
    set obs $numgroups
    gen g=_n
    generate y=rnormal()
    expand $groupobs
    sort g
    by g: gen l=_n

    generate double one=1
    generate double x2=1+rnormal(2)
    generate double x3=0.2*x2-1.2*rnormal(1)
    generate double x5=0.2*x2+0.2*x3+rnormal(0)
    generate double x4=(x5>0)
    scalar b1=1
    scalar b2=1
    scalar b3=1
    scalar b4=1

    *MATA TO CREATE CORRELATION MATRIX
    mata: cuicui()
    by g: replace y=b1+x2*b2+x3*b3+x4*b4+u

end


mata:
void cuicui() 
{
        G=strtoreal(st_global("numgroups"))
        L=strtoreal(st_global("groupobs"))
        N=G*L
        rho=0.6
        A=(J(sqrt(G),1,1))#((J(sqrt(G),1,1))#((1..sqrt(L))'#J(sqrt(L),1,1)))
        B=2*(0::sqrt(G)-1)#(J(sqrt(G)*L,1,1))
        I=A+B
        U=J(G,1,1)#(J(sqrt(L),1,1)#(sqrt(L)::1))
        W=J(sqrt(G),1,1)#((2*(sqrt(G)-1::0))#J(L,1,1))
        J=U+W
        I=J(N,1,1)'#I
        J=J(N,1,1)'#J
        D=sqrt((I-I'):^2+(J-J'):^2)
        W=((1-rho):*I(N)+rho:/(D+I(N))) 
        e=rnormal(N,1,0,1)
        C=-0.5:+cholesky(W)*e
        u=exp(C)
        st_addvar("double", "u")
        st_store(.,"u",u)
}
end


set seed 10101
simulate y x2 x3 x4, reps($numsims) saving (gee, replace) nolegend: gee

[Cuicui]

Christophe, thank you so much! It runs well this time!