This a part of a simulation based on an empirical paper.Here we simulate the process about 2000 times. In each time, we get a DGP generating 2000 observations.

clear all
set matsize 5000
program drop _all
drop _all
matrix drop _all

global matrix

program  sim1
  version 10
syntax [, obs(integer 1) ]

set obs `obs'
tempvar cons x1 nj njh gap w1  de1 de2 de3 de4 u0 e1 e2 e3 e4

gen x1=invnormal(uniform())
gen u0=uniform()
gen nj=1 if u0>0.2
replace nj=0 if u0<=0.2
gen w1=4.25+1.3*uniform()
gen njh=1 if w1>=5.05
replace njh=0 if w1<5.05
gen gap=0 if u0<=0.2
replace gap=0 if w1>=5.05
replace gap=(5.05-w1)/w1 if w1<5.05
gen e1=invnormal(uniform())
gen e2=invnormal(uniform())
gen e3=invnormal(uniform())
gen e4=invnormal(uniform())
gen cons=1.99
gen de1=cons+2*x1+2.3*nj+e1
gen de2=cons+2*x1+2.01*gap+e2
gen de3=cons+2*x1+2.01*gap+2.3*nj+e3
gen de4=cons+2*x1+2.01*gap+2.3*nj+3.32*njh+e4

local mat t
quietly reg de1 x1 nj
matrix t=nullmat(t)e(b)
quietly reg de2 x1 nj
matrix t=(nullmat(t),e(b))
quietly reg de3 x1 nj
matrix t=(nullmat(t),e(b))
quietly reg de4 x1 nj
matrix t=(nullmat(t),e(b))

matrix tt=(nullmat(tt)t)
 drop _all
 matrix drop t
end

simulate, reps(2000): sim1, obs(2000)

svmat tt
matrix drop _all
quietly summarize tt1
matrix t=nullmat(t)(r(mean),r(Var))
quietly summarize tt2
matrix t=t(r(mean),r(Var))
quietly summarize tt3
matrix t=t(r(mean),r(Var))
quietly summarize tt4
matrix t=t(r(mean),r(Var))
quietly summarize tt5
matrix t=t(r(mean),r(Var))
quietly summarize tt6
matrix t=t(r(mean),r(Var))
quietly summarize tt7
matrix t=t(r(mean),r(Var))
quietly summarize tt8
matrix t=t(r(mean),r(Var))
quietly summarize tt9
matrix t=t(r(mean),r(Var))
quietly summarize tt10
matrix t=t(r(mean),r(Var))
quietly summarize tt11
matrix t=t(r(mean),r(Var))
quietly summarize tt12
matrix t=t(r(mean),r(Var))

drop _all
matrix colnames t = mean var
svmat t, names(col)
outsheet using F:Statafastfood1.xls