I have written the following do file to conduct some Monte Carlo simulation. However, for some reason the mata function, controlf(), is not being read while executing the program "cfsimu." When I open up the cfsimu and put the mata function just before it is called, then I don't get the problem. Any help would be most welcomed.
Best wishes
Amaresh
cd "C:\DATA\MCExp"
set more off
clear all
program drop _all
program define cfsimu, rclass
version 15
drop _all
set obs 500 // n= 500 individuals (change the number of Individuals here)
gen id = _n // idendification number of the individuals
//sz1 = 5; sz2 = 2; st = 4; sa1 = 6; sa2= 1
//r12 = 0.25;
//r1t =0.1 ; r2t =0.15 ;
//r1a1 =0.2; r2a1= 0.3; r1a2 =0.25 ; r2a2 =0.3 ; ra1a2= 0.5
//rta1 = 0.5; rta2 = 0.25
mat ZTA = (25, 2.5, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6, 1.25 \ ///
2.5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 25, 2.5, 0, 0, 0, 0, 0, 0, 2, 6, 1.25 \ ///
0, 0, 2.5, 4, 0, 0, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 25, 2.5, 0, 0, 0, 0, 2, 6, 1.25 \ ///
0, 0, 0, 0, 2.5, 4, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 0, 0, 25, 2.5, 0, 0, 2, 6, 1.25 \ ///
0, 0, 0, 0, 0, 0, 2.5, 4, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 0, 0, 0, 0, 25, 2.5, 2, 6, 1.25 \ ///
0, 0, 0, 0, 0, 0, 0, 0, 2.5, 4, 1.2, 3.6, 0.6 \ ///
2, 1.2, 2, 1.2, 2, 1.2, 2, 1.2, 2, 1.2, 16, 12, 1 \ ///
6, 3.6, 6, 3.6, 6, 3.6, 6, 3.6, 6, 3.6, 12, 36, 3 \ ///
1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1, 3, 1 )
mat ZTA = (1/36)*ZTA
drawnorm z11 z21 z12 z22 z13 z23 z14 z24 z15 z25 th1 al11 al21, cov(ZTA)
gen al12 = al11
gen al13 = al11
gen al14 = al11
gen al15 = al11
gen al22 = al21
gen al23 = al21
gen al24 = al21
gen al25 = al21
gen th2 = th1
gen th3 = th1
gen th4 = th1
gen th5 = th1
reshape long z1 z2 al1 al2 th, i(id) j(time)
/*Discretizing z1_it and z2_it*/
gen dz1 = (z1>0)
gen dz2 = (z2>.25)
bysort id: egen mdz1 = mean(dz1) // computing the mean of the z1_it for every individial
bysort id: egen mdz2 = mean(dz2) // computing the mean of the z2_it for every individial
/*EE is the covariance matrix, Σ_ζϵ, in section 3 of the paper.*/
matrix EE =( 1,.75,.25 \ .75, 1, .5 \ .25, .5, 1)
/*Generating the idiosyncratic errors, ζ_it (et) and ϵ_it (ep).*/
drawnorm ep1 ep2 et, cov(EE)
/*Generating the endogenous explanatory variable, x_it.*/
gen x1 = -1*dz1 + .05*dz2 + al1 + ep1
gen x2 = .025*dz1 + .75*dz2 + al2 + ep2
/*Generating the endogenous binary outcome, y_it.*/
gen ys = -1*x1 + 0.5*x2 + th + et
gen y = (ys >0)
matrix drop _all
/*Estimating the reduced form equation.*/
xtset id time
xtsur (x1 dz1 dz2 mdz1 mdz2 )(x2 dz1 dz2 mdz1 mdz2 )
mat b = e(b)
mat B1 = b[1, "x1:"]
mat B2 = b[1, "x2:"]
mat VAL = e(sigma_u)
mat VEP = e(sigma_e)
gen mzb1 = 0 //Chamberlain-Mundlak specification for conditional expectation of $\alpha_{x1}$ given exogenous regressors
gen mzb2 = 0 //Chamberlain-Mundlak specification for conditional expectation of $\alpha_{x2}$ given exogenous regressors
foreach var of local mzvar {
mat b1`var' = b[1, "x1:`var'"]
sca sb1`var' = b1`var'[1,1]
replace mzb1 = mzb1+`var'*sb1`var'
mat b2`var' = b[1, "x2:`var'"]
sca sb2`var' = b2`var'[1,1]
replace mzb2 = mzb2+`var'*sb2`var'
}
local idvar "id"
local xvar "x1 x2"
local zvar "dz1 dz2"
local mzvar "mdz1 mdz2"
local mzbvar "mzb1 mzb2" // Conditional expectation of α_1i and α_2i given exogenous variables Z_i
sort id time
/*Calling the mata function, controlf, to compute the control functions, \hat{α}_i= ex_x and \hat{ϵ}_it=ep_x*/
marksample touse
mata: controlf("B1", "B2", "VAL", "VEP", "`idvar'", "`xvar'", "`zvar'", "`mzvar'", "`mzbvar'", "`touse'")
/*Estimating the structural equation augmented with the control functions, \hat{α}_i and \hat{ϵ}_it.*/
xtset id time
probit y x1 x2 ex_x1 ex_x2 ep_x1 ep_x2, nocon
mat BS = e(b) //Saving the coefficients of the structural equation to subsequently use them for computing the average partial effect.
replace x1 = 0.5 // To compute the average partial effect at x1=0.5
replace x2 = 1 // To compute the average partial effect at x2=1
local svar "x1 x2 ex_x1 ex_x2 ep_x1 ep_x2"
gen xsb = 0
foreach var of local svar {
mat bs`var' = BS[1, "y:`var'"]
sca sbs`var' = bs`var'[1,1]
replace xsb = xsb+`var'*sbs`var'
}
gen pe1_smp = sbsx1*normalden(xsb) // average partial effect of x1 at x1=0.5 & x2=1 from the CRECF estimator for every individual i and time t
gen pe2_smp = sbsx2*normalden(xsb) // average partial effect of x2 at x1=0.5 & x2=1 from the CRECF estimator for every individual i and time t
gen pe1_pop = -1*normalden( -1*0.5 + 0.5 + th ) // true average partial effect of x1 evaluated at x1= 0.5 and x2 = 1 for every individual i and time t
gen pe2_pop = 0.5*normalden( -1*0.5 + 0.5 + th ) // true average partial effect of x2 evaluated at x1= 0.5 and x2 = 1 for every individual i and time t
sca drop _all
su pe1_smp
return scalar ape1_smp = r(mean) // average partial effect (APE) of x1 from the CRECF estimator
su pe1_pop
return scalar ape1_pop = r(mean) // true average partial effect of x1
su pe2_smp
return scalar ape2_smp = r(mean) // average partial effect (APE) of x2 from the CRECF estimator
su pe2_pop
return scalar ape2_pop = r(mean) // true average partial effect of x2
end
/* Mata function to compute control functions */
clear mata
mata:
function controlf(string scalar B1, string scalar B2, string scalar VAL, string scalar VEP, string scalar idvar, string scalar xvar, string scalar zvar, string scalar mzvar, string scalar mzbvar, touse)
{
st_view(ID=., ., tokens(idvar), touse ) // Variable identifying individuals
st_view(X=., ., tokens(xvar), touse) // Exogenous regressors
st_view(Z=., ., tokens(zvar), touse) // Exogenous regressors
st_view(MZ=., ., tokens(mzvar), touse) // Means of Exogenous regressors
st_view(MZB=., ., tokens(mzbvar), touse) // Conditional expectation of α_1i and α_2i given exogenous variables Z_i
BC1 = st_matrix(B1)
BC2 = st_matrix(B2)
V_AL = st_matrix(VAL) //Estimate of the variance of α_i
V_EP = st_matrix(VEP) //Estimate of the variance of ϵ_it
ZV = Z,MZ
IAL = invsym(V_AL)
IEP = invsym(V_EP)
SG = (invsym(5*IEP + IAL))*IEP // 5 is the number of time periods
R1 = X[.,1] - ZV*BC1' //Residuals
R2 = X[.,2] - ZV*BC2' //Residuals
R = R1,R2
EX_POS1 =J(1,1,0); EX_POS2 =J(1,1,0)
info = panelsetup(ID, 1, 5, 5)
for (i=1; i<=rows(info); i++) {
/*i for each individual*/
RI= panelsubmatrix(R, i, info)
SRI=J(1,cols(RI),0)
j=1
while (j<=rows(RI)){
RIJ = RI[j,.]
SRI=SRI + RIJ
j=j+1
}
EX_POSI = SG*SRI'
EX_POS1= EX_POS1 \ J(rows(RI),1,1)*EX_POSI[1,.]
EX_POS2= EX_POS2 \ J(rows(RI),1,1)*EX_POSI[2,.]
}
EX_POS1 = EX_POS1[2..rows(EX_POS1),.]; EX_POS2 = EX_POS2[2..rows(EX_POS2),.]
EX_POSM1 = MZB[.,1] + EX_POS1 //$\hat{\alpha}_{i}$ for income
EX_POSM2 = MZB[.,2] + EX_POS2 //$\hat{\alpha}_{i}$ for index of productive assets
ER1 = R1 - EX_POS1 //$\hat{\epsilon}_{it}$ for x1
ER2 = R2 - EX_POS2 //$\hat{\epsilon}_{it}$ for x2
EX_POSM = EX_POSM1, EX_POSM2
ER = ER1, ER2
NVAR= EX_POSM, ER
NVAR= EX_POSM, ER
st_addvar("float", ( "ex_x1", "ex_x2", "ep_x1", "ep_x2"))
st_store(., ( "ex_x1", "ex_x2", "ep_x1", "ep_x2"), NVAR) //the control functions, \hat{α}_i= ex_x and \hat{ϵ}_it=ep_x
}
end
set seed 339487731
simulate ape1_smp=r(ape1_smp) ape1_pop=r(ape1_pop) ape2_smp=r(ape2_smp) ape2_pop=r(ape2_pop) _b, reps(100) : cfsimu //simulating 10000 averrage partial effects for 500 individuals
gen dape1_cf = ape1_smp- ape1_pop //Estimated APE of x1 - True APE of x1
gen ser1 = (dape1_cf)^2
gen sser1 = sum(ser1)
gen rmse1 = sqrt(sser1[_N]/_N) //computing root mean square
gen dape2_cf = ape2_smp- ape2_pop //Estimated APE of x1 - True APE of x1
gen ser2 = (dape2_cf)^2
gen sser2 = sum(ser2)
gen rmse2 = sqrt(sser2[_N]/_N) //computing root mean square
save "data_crecf.dta", replace // saving 10000 simulated averrage partial effects
Best wishes
Amaresh
cd "C:\DATA\MCExp"
set more off
clear all
program drop _all
program define cfsimu, rclass
version 15
drop _all
set obs 500 // n= 500 individuals (change the number of Individuals here)
gen id = _n // idendification number of the individuals
//sz1 = 5; sz2 = 2; st = 4; sa1 = 6; sa2= 1
//r12 = 0.25;
//r1t =0.1 ; r2t =0.15 ;
//r1a1 =0.2; r2a1= 0.3; r1a2 =0.25 ; r2a2 =0.3 ; ra1a2= 0.5
//rta1 = 0.5; rta2 = 0.25
mat ZTA = (25, 2.5, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6, 1.25 \ ///
2.5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 25, 2.5, 0, 0, 0, 0, 0, 0, 2, 6, 1.25 \ ///
0, 0, 2.5, 4, 0, 0, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 25, 2.5, 0, 0, 0, 0, 2, 6, 1.25 \ ///
0, 0, 0, 0, 2.5, 4, 0, 0, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 0, 0, 25, 2.5, 0, 0, 2, 6, 1.25 \ ///
0, 0, 0, 0, 0, 0, 2.5, 4, 0, 0, 1.2, 3.6, 0.6 \ ///
0, 0, 0, 0, 0, 0, 0, 0, 25, 2.5, 2, 6, 1.25 \ ///
0, 0, 0, 0, 0, 0, 0, 0, 2.5, 4, 1.2, 3.6, 0.6 \ ///
2, 1.2, 2, 1.2, 2, 1.2, 2, 1.2, 2, 1.2, 16, 12, 1 \ ///
6, 3.6, 6, 3.6, 6, 3.6, 6, 3.6, 6, 3.6, 12, 36, 3 \ ///
1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1.25, 0.6, 1, 3, 1 )
mat ZTA = (1/36)*ZTA
drawnorm z11 z21 z12 z22 z13 z23 z14 z24 z15 z25 th1 al11 al21, cov(ZTA)
gen al12 = al11
gen al13 = al11
gen al14 = al11
gen al15 = al11
gen al22 = al21
gen al23 = al21
gen al24 = al21
gen al25 = al21
gen th2 = th1
gen th3 = th1
gen th4 = th1
gen th5 = th1
reshape long z1 z2 al1 al2 th, i(id) j(time)
/*Discretizing z1_it and z2_it*/
gen dz1 = (z1>0)
gen dz2 = (z2>.25)
bysort id: egen mdz1 = mean(dz1) // computing the mean of the z1_it for every individial
bysort id: egen mdz2 = mean(dz2) // computing the mean of the z2_it for every individial
/*EE is the covariance matrix, Σ_ζϵ, in section 3 of the paper.*/
matrix EE =( 1,.75,.25 \ .75, 1, .5 \ .25, .5, 1)
/*Generating the idiosyncratic errors, ζ_it (et) and ϵ_it (ep).*/
drawnorm ep1 ep2 et, cov(EE)
/*Generating the endogenous explanatory variable, x_it.*/
gen x1 = -1*dz1 + .05*dz2 + al1 + ep1
gen x2 = .025*dz1 + .75*dz2 + al2 + ep2
/*Generating the endogenous binary outcome, y_it.*/
gen ys = -1*x1 + 0.5*x2 + th + et
gen y = (ys >0)
matrix drop _all
/*Estimating the reduced form equation.*/
xtset id time
xtsur (x1 dz1 dz2 mdz1 mdz2 )(x2 dz1 dz2 mdz1 mdz2 )
mat b = e(b)
mat B1 = b[1, "x1:"]
mat B2 = b[1, "x2:"]
mat VAL = e(sigma_u)
mat VEP = e(sigma_e)
gen mzb1 = 0 //Chamberlain-Mundlak specification for conditional expectation of $\alpha_{x1}$ given exogenous regressors
gen mzb2 = 0 //Chamberlain-Mundlak specification for conditional expectation of $\alpha_{x2}$ given exogenous regressors
foreach var of local mzvar {
mat b1`var' = b[1, "x1:`var'"]
sca sb1`var' = b1`var'[1,1]
replace mzb1 = mzb1+`var'*sb1`var'
mat b2`var' = b[1, "x2:`var'"]
sca sb2`var' = b2`var'[1,1]
replace mzb2 = mzb2+`var'*sb2`var'
}
local idvar "id"
local xvar "x1 x2"
local zvar "dz1 dz2"
local mzvar "mdz1 mdz2"
local mzbvar "mzb1 mzb2" // Conditional expectation of α_1i and α_2i given exogenous variables Z_i
sort id time
/*Calling the mata function, controlf, to compute the control functions, \hat{α}_i= ex_x and \hat{ϵ}_it=ep_x*/
marksample touse
mata: controlf("B1", "B2", "VAL", "VEP", "`idvar'", "`xvar'", "`zvar'", "`mzvar'", "`mzbvar'", "`touse'")
/*Estimating the structural equation augmented with the control functions, \hat{α}_i and \hat{ϵ}_it.*/
xtset id time
probit y x1 x2 ex_x1 ex_x2 ep_x1 ep_x2, nocon
mat BS = e(b) //Saving the coefficients of the structural equation to subsequently use them for computing the average partial effect.
replace x1 = 0.5 // To compute the average partial effect at x1=0.5
replace x2 = 1 // To compute the average partial effect at x2=1
local svar "x1 x2 ex_x1 ex_x2 ep_x1 ep_x2"
gen xsb = 0
foreach var of local svar {
mat bs`var' = BS[1, "y:`var'"]
sca sbs`var' = bs`var'[1,1]
replace xsb = xsb+`var'*sbs`var'
}
gen pe1_smp = sbsx1*normalden(xsb) // average partial effect of x1 at x1=0.5 & x2=1 from the CRECF estimator for every individual i and time t
gen pe2_smp = sbsx2*normalden(xsb) // average partial effect of x2 at x1=0.5 & x2=1 from the CRECF estimator for every individual i and time t
gen pe1_pop = -1*normalden( -1*0.5 + 0.5 + th ) // true average partial effect of x1 evaluated at x1= 0.5 and x2 = 1 for every individual i and time t
gen pe2_pop = 0.5*normalden( -1*0.5 + 0.5 + th ) // true average partial effect of x2 evaluated at x1= 0.5 and x2 = 1 for every individual i and time t
sca drop _all
su pe1_smp
return scalar ape1_smp = r(mean) // average partial effect (APE) of x1 from the CRECF estimator
su pe1_pop
return scalar ape1_pop = r(mean) // true average partial effect of x1
su pe2_smp
return scalar ape2_smp = r(mean) // average partial effect (APE) of x2 from the CRECF estimator
su pe2_pop
return scalar ape2_pop = r(mean) // true average partial effect of x2
end
/* Mata function to compute control functions */
clear mata
mata:
function controlf(string scalar B1, string scalar B2, string scalar VAL, string scalar VEP, string scalar idvar, string scalar xvar, string scalar zvar, string scalar mzvar, string scalar mzbvar, touse)
{
st_view(ID=., ., tokens(idvar), touse ) // Variable identifying individuals
st_view(X=., ., tokens(xvar), touse) // Exogenous regressors
st_view(Z=., ., tokens(zvar), touse) // Exogenous regressors
st_view(MZ=., ., tokens(mzvar), touse) // Means of Exogenous regressors
st_view(MZB=., ., tokens(mzbvar), touse) // Conditional expectation of α_1i and α_2i given exogenous variables Z_i
BC1 = st_matrix(B1)
BC2 = st_matrix(B2)
V_AL = st_matrix(VAL) //Estimate of the variance of α_i
V_EP = st_matrix(VEP) //Estimate of the variance of ϵ_it
ZV = Z,MZ
IAL = invsym(V_AL)
IEP = invsym(V_EP)
SG = (invsym(5*IEP + IAL))*IEP // 5 is the number of time periods
R1 = X[.,1] - ZV*BC1' //Residuals
R2 = X[.,2] - ZV*BC2' //Residuals
R = R1,R2
EX_POS1 =J(1,1,0); EX_POS2 =J(1,1,0)
info = panelsetup(ID, 1, 5, 5)
for (i=1; i<=rows(info); i++) {
/*i for each individual*/
RI= panelsubmatrix(R, i, info)
SRI=J(1,cols(RI),0)
j=1
while (j<=rows(RI)){
RIJ = RI[j,.]
SRI=SRI + RIJ
j=j+1
}
EX_POSI = SG*SRI'
EX_POS1= EX_POS1 \ J(rows(RI),1,1)*EX_POSI[1,.]
EX_POS2= EX_POS2 \ J(rows(RI),1,1)*EX_POSI[2,.]
}
EX_POS1 = EX_POS1[2..rows(EX_POS1),.]; EX_POS2 = EX_POS2[2..rows(EX_POS2),.]
EX_POSM1 = MZB[.,1] + EX_POS1 //$\hat{\alpha}_{i}$ for income
EX_POSM2 = MZB[.,2] + EX_POS2 //$\hat{\alpha}_{i}$ for index of productive assets
ER1 = R1 - EX_POS1 //$\hat{\epsilon}_{it}$ for x1
ER2 = R2 - EX_POS2 //$\hat{\epsilon}_{it}$ for x2
EX_POSM = EX_POSM1, EX_POSM2
ER = ER1, ER2
NVAR= EX_POSM, ER
NVAR= EX_POSM, ER
st_addvar("float", ( "ex_x1", "ex_x2", "ep_x1", "ep_x2"))
st_store(., ( "ex_x1", "ex_x2", "ep_x1", "ep_x2"), NVAR) //the control functions, \hat{α}_i= ex_x and \hat{ϵ}_it=ep_x
}
end
set seed 339487731
simulate ape1_smp=r(ape1_smp) ape1_pop=r(ape1_pop) ape2_smp=r(ape2_smp) ape2_pop=r(ape2_pop) _b, reps(100) : cfsimu //simulating 10000 averrage partial effects for 500 individuals
gen dape1_cf = ape1_smp- ape1_pop //Estimated APE of x1 - True APE of x1
gen ser1 = (dape1_cf)^2
gen sser1 = sum(ser1)
gen rmse1 = sqrt(sser1[_N]/_N) //computing root mean square
gen dape2_cf = ape2_smp- ape2_pop //Estimated APE of x1 - True APE of x1
gen ser2 = (dape2_cf)^2
gen sser2 = sum(ser2)
gen rmse2 = sqrt(sser2[_N]/_N) //computing root mean square
save "data_crecf.dta", replace // saving 10000 simulated averrage partial effects
Comment