********************************************************************************
*! combreg, v4, GCerulli, 12/10/2018
********************************************************************************
program combreg, eclass
version 14
#delimit ;
syntax varlist [if] [in] [fweight pweight iweight] ,
model(string)
s(numlist max=1 integer)
k(numlist max=1 integer)
seed(numlist max=1)
[
factors(varlist numeric)
vce(string)
graph
];
#delimit cr
********************************************************************************
quietly{ // open quietly
********************************************************************************
marksample touse
markout `touse' `factors'
********************************************************************************
* Keep only the "if", and save the initial dataset into "mydata"
********************************************************************************
tempfile mydata
qui save `mydata' , replace
keep if `touse'
********************************************************************************
* Generate dummies for factor variables
********************************************************************************
foreach V of local factors{
levelsof `V' , local(L)
local M: word count of `L'
local H=`M'-1
forvalues i=1/`M'{
cap drop _`V'`i'
}
qui tab `V' if `touse' , gen(_`V') mis
local sum_`V' ""
forvalues i=1/`H'{
local sum_`V' `sum_`V'' _`V'`i'
}
}
********************************************************************************
local sum_tot ""
foreach V of local factors{
local sum_tot `sum_tot' `sum_`V''
}
********************************************************************************
tokenize `varlist' `sum_tot'
local y `1'
local w `2'
macro shift
macro shift
local xvars `*'
********************************************************************************
local m=`s'
********************************************************************************
preserve
set seed `seed'
local varlist `xvars'
local N : word count `varlist'  //   N = # of covariates
********************************************************************************
* Warning 1
********************************************************************************
if `s'>=`N'{
break
di _newline(2)
di as result in red "**********************************************************"
di as result in red "Warning: 's' must be lower than the number of covariates  "
di as result in red "considered in the benchmark model.                        "
di as result in red "**********************************************************"
exit
}
********************************************************************************
clear
********************************************************************************
set obs `N'
tempvar X
gen `X'="."
********************************************************************************
local i=1
foreach x of local varlist{
replace `X'="`x'" in `i'
local i=`i'+1
}
replace `X'="" if `X'=="."
qui save `X' , replace
levelsof `X' , local(L) clean
di `"`L'"'
********************************************************************************
* K = number of re-samples
********************************************************************************
forvalues j=`s'/`m'{
tempfile D`j'
}
********************************************************************************
forvalues j=`s'/`m'{
forvalues i=1/`k'{
tempfile S_`j'_`i'
}
}
********************************************************************************
forvalues j=`s'/`m'{
use `X' , clear
sample `j' , count
cap drop id
gen id = _n
drop `X'
save `D`j'' , replace
forvalues i=1/`k'{
use `X' , clear
sample `j' , count
cap drop id
gen id = _n
rename `X' `X'`j'`i'
save `S_`j'_`i'' , replace
use `D`j'' , clear
merge 1:1 id using `S_`j'_`i''
drop _merge
save `D`j'' , replace
}
}
********************************************************************************
forvalues j=`s'/`m'{
use `D`j'' , clear
forvalues i=1/`k'{
levelsof `X'`j'`i' , local(" L`j'`i'")  clean
di `"L`j'`i'"'
}
}
cap erase `X'.dta
restore
********************************************************************************
* Model with all variables (Baseline model)
********************************************************************************
if "`model'"=="reg"{
xi: reg `y' `w' `varlist' if `touse' [`weight'`exp'] ,  vce(`vce')
local ATET_original=_b[`w']
local se_att=_se[`w']
local Tstud_original=abs(_b[`w']/_se[`w'])
local lim=1.96*_se[`w']
local upper=_b[`w']+`lim'
local lower=_b[`w']-`lim'
********************************************************************************
local B=`k'*(`m'-`s'+1)
********************************************************************************
tempname J
mat `J'=J(`B',2,.)
********************************************************************************
* Run the regressions/matchings for estimating ATET
local h=1
forvalues j=`s'/`m'{
forvalues i=1/`k'{
xi: reg `y' `w' `L`j'`i'' if `touse' [`weight'`exp'] ,  vce(`vce')
********************************************************************************
mat `J'[`h',1] = _b[`w']
mat `J'[`h',2] = abs(_b[`w']/_se[`w'])
local h=`h'+1
}
}
}
********************************************************************************
else if "`model'"=="match"{
xi: psmatch2 `w' `varlist' if `touse' , out(`y') ate
local ATET_original=r(att)
local se_att=r(seatt)
local Tstud_original=r(att)/r(seatt)
local lim=1.96*`se_att'
local upper=`ATET_original'+`lim'
local lower=`ATET_original'-`lim'
********************************************************************************
local B=`k'*(`m'-`s'+1)
********************************************************************************
tempname J
mat `J'=J(`B',2,.)
********************************************************************************
* Run the regressions/matchings for estimating ATET
local h=1
forvalues j=`s'/`m'{
forvalues i=1/`k'{
xi: psmatch2 `w' `L`j'`i'' if `touse' , out(`y') ate
*cap drop near_obs*
*teffects psmatch (`y') (`w' `L`j'`i'' , probit) , atet generate(near_obs) atet
mat `J'[`h',1] = r(att)
mat `J'[`h',2] = r(att)/r(seatt)
local h=`h'+1
}
}
}
********************************************************************************
* Simulation results
********************************************************************************
cap drop _ATET*
svmat `J', names(_ATET)
local ate_bench=round(`ATET_original',0.01)
local Tstud_bench=round(`Tstud_original',0.01)
if "`graph'"=="graph"{
kdensity _ATET1 if `touse' , xline(`ATET_original',lpattern(dash)) xtitle("ATET") note("Number of simulations: `k'" ///
"Number of covariates: `s' out of `N'" "Reference ATET: `ate_bench'" ) scheme(s1mono) title("")
}
********************************************************************************
* Testing whether there are differences between model Ho and the
* average of the simulated models
********************************************************************************
qui reg _ATET1  if `touse'
test _cons=`ATET_original'
********************************************************************************
* Returns
********************************************************************************
ereturn clear
qui sum _ATET1 if `touse' , d
local mean_b_sim=r(mean)
ereturn scalar mean_b_sim = r(mean)
ereturn scalar median_b_sim = r(p50)
ereturn scalar sd_b_sim = r(sd)
ereturn scalar b_bench = `ate_bench'
ereturn scalar t_bench = `Tstud_bench'
local delta=abs((`mean_b_sim'-`ate_bench')/`ate_bench')
ereturn scalar delta = `delta'
ereturn scalar upci = `upper'
ereturn scalar lowci = `lower'
cap drop _ATET1
********************************************************************************
qui sum _ATET2 if `touse' , d
local mean_Tstud_sim=r(mean)
ereturn scalar mean_Tstud_sim = r(mean)
cap drop _ATET2
********************************************************************************
} // end quietly
********************************************************************************
qui use `mydata' , clear
end
********************************************************************************
* END
********************************************************************************