*! Ver 1.71 bug with weights
* Ver 1.7 adds version for easier update
* Ver 1.68 Bug with drimp
* Ver 1.67 Bug with IPW
* Ver 1.66 Bug with binit, skip fixed
* Ver 1.65 Bug with binit. Fixed
* Ver 1.64 added binit for faster csdid (if slighly)
* Ver 1.63 added Dryrun
* Ver 1.63 Bug with option all and RC estimator
* Ver 1.62 Added level
* Ver 1.61 Change sample for drimp
* Ver 1.6 Change in check for 2x2 data. And updated site!
* added correction to teffects for ALL
* Ver 1.5 New output, with panel GMM estimators. Also WB standard errors with cluster
* Ver 1.38 Adding Cluster Stantandard errors
* version 1.37 2jun2021 Add extra messages of 2x2 balance. Checks that you indeed have panel data
* version 1.36 17may2021 Changes with EP display options
** Goal. Make the estimator modular. That way other options can be added
* v1.35 DRDID for Stata
* Added by AN and adapted by FRA
* Added Version control (v14), Variables are no longer left after DRDID. But option stub is added
* for later, as CSDID will use those
* v1.31 DRDID for Stata by FRA. Adding program for Bootstrap multiplier
* v1.3 DRDID for Stata by FRA Changing ATT RIF creation
* v1.2 DRDID for Stata by FRA All Estimators are ready For panel and RC
* Next help!
* v1.0 DRDID for Stata by FRA All Estimators but IPWRA have are available for panel
* Need to add weights.
* v0.8 DRDID for Stata by FRA Other estimators are available. Onlyone missing iwp panel
* v0.7 DRDID for Stata by FRA IPW estimator for panel (Asjad original Rep)
* v0.5 DRDID for Stata by FRA Incorporates RC1 and RC2 estimators
* v0.2 DRDID for Stata by FRA Fixes typo with tag
* v0.2 DRDID for Stata by FRA Allows for Factor notation
* v0.1 DRDID for Stata by FRA Typo with ID TIME
/* !! epg: added method gmm . Also added parsing routing _Vce_Parse
gmm and vce(...)
*/
/* !! epg: what to do with bootstrap csboot(csboot_opts)
csboot_opts -- reps()
-- rseed()
-- wbtype()
-- ...
Perfecto ya todo fue integrado
!! Crear la tabla Funciona
!! Necesitamos regresar e(V) o solamente e(b)?
Creo que ya esta.
!! Deberia cambiar el t estadistico?
*/
/*
csdid loop over
*/
program define drdid, eclass byable(onecall)
version 14
if _by() {
local BY `"by `_byvars'`_byrc0':"'
}
`BY' _vce_parserun drdid, noeqlist jkopts(eclass): `0'
if "`s(exit)'" != "" {
ereturn local cmdline `"drdid `0'"'
exit
}
syntax [anything(everything)] [iw pw aw], [* version]
/**Version**/
if "`version'"!="" {
display "version: 1.71"
addr scalar version = 1.71
exit
}
if replay() {
if `"`e(cmd)'"' != "drdid" {
error 301
}
else if _by() {
error 190
}
else {
Display `0'
}
exit
}
if runiform()<.001 {
easter_egg
}
`vv' ///
`BY' drdid_wh `0'
ereturn local cmdline `"drdid `0'"'
end
mata
void is_2x2(string scalar time, treat, touse, isok){
real matrix dta
dta=st_data(.,time,touse),st_data(.,treat,touse)
dta=uniqrows(dta)
if ((rows(dta)==4) & (rows(uniqrows(dta[,1]))==2) & (rows(uniqrows(dta[,2]))==2)) st_numscalar(isok,1)
else st_numscalar(isok,0)
}
void is_balp(string scalar ivar, touse, balp){
real matrix dta
dta=st_data(.,ivar,touse)
st_numscalar(balp,mean(uniqrows(dta,2)[,2]))
}
end
program define drdid_wh, eclass sortpreserve byable(recall)
syntax varlist(fv numeric) [if] [in] [iw], ///
[Ivar(varname)] ///
Time(varname) ///
TReatment(varname) ///
[noisily ///
drimp ///
dripw ///
reg ///
stdipw ///
ipw ///
ipwra ///
all ///
rc1 ///
WBOOT(string) ///
WBOOT1 ///
*reps(int 999) ///
*wbtype(int 1) /// Hidden option
rseed(str) /// set seed
Level(int 95) /// CI level
stub(name) replace /// to avoid overwritting
cluster(varname) /// For Cluster
vce(string) ///
gmm ///
pscore(string) ///
csdid ///
binit(string) ///
dryrun ///
* ///
]
* so it returns Nothing
if "`dryrun'"!="" error 1111
_get_diopts diopts other, `options'
quietly capture Display, `diopts' `other' level(`level')
if _rc==198 {
Display, `diopts' `other' level(`level')
}
marksample touse
markout `touse' `ivar' `time' `treatment' `cluster'
_Vce_Parse if `touse', `gmm' `wboot1' wboot(`wboot') vce(`vce')
** Move this into VCE_parse
*if "`binit'"!="" local binit `binit', skip
local semethod "`r(semethod)'"
if ("`semethod'"=="wildboot") {
local wboot "wboot"
local reps = r(reps)
local wbtype = r(wbtype)
}
** Verify 2x2 data
tempname isok
mata:is_2x2("`time'","`treatment'","`touse'","`isok'")
if scalar(`isok')==0 {
display in red "You do not have a 2x2 design."
error 555
}
if "`cluster'"=="" local cluster "`r(cluster)'"
if "`cluster'"!="" {
tempvar clvar
qui:egen double `clvar' = group(`cluster') if `touse'
local ocluster `cluster'
local cluster `clvar'
}
if "`rseed'"=="" local rseed "`r(seed)'"
capture:set seed `rseed'
**# Verifies if data is panel data when "ivar" is declared
if "`ivar'"!="" {
capture xtset
if _rc!=0 {
qui:xtset `ivar' `time'
qui:xtset , clear
}
if "`cluster'"!="" {
_xtreg_chk_cl2 `cluster' `ivar'
}
**# verify ALL data is locally balanced
tempname balp
*mata:is_balp("`ivar'", "`touse'", "`balp'")
qui:bysort `touse' `ivar':gen `balp'=_N if `touse'
sum `balp', meanonly
if r(mean)<2 {
display in red "{p}Some panel units are observed only once.{p_end}" _n ///
"{p}Those observations will be excluded from the sample. " ///
"If you want to keep them, do not use `ivar' as the panel id -ivar-, " ///
"which will request repeated crossection estimators {p_end}"
quietly replace `touse' = 0 if `balp'==1
}
}
**# Verifies If variable exists. For CSDID
cap unab allnew: `stub'att*
if "`stub'"!="" & "`replace'"=="" & "`allnew'"!="" {
foreach v of var `allnew' {
conf new var `v'
}
}
*** Add Char to variables. They will be ours. Perhaps for CSDID
**# Verifies REPS is a number
** Which option
if "`drimp'`dripw'`reg'`stdipw'`ipw'`ipwra'`all'"=="" {
display "No estimator selected. Using default {bf:drimp}"
local drimp drimp
}
else {
if `:word count `drimp' `dripw' `reg' `stdipw' `ipw' `ipwra' `all' '!=1 {
display "Only one option allowed, more than 1 selected."
error 1
}
}
_S__est, `drimp' `dripw' `reg' `stdipw' `ipw' `ipwra' `all'
local estimator "`s(estimator)'"
** First determine outcome and xvars
gettoken y xvar:varlist
** Sanity Checks for Time. Only 2 values
** just in case for xvar not be empty
*local xvar `xvar'
tempvar vals vals2 vals3
tempvar tmt
qui:egen byte `tmt'=group(`time') if `touse'
qui:replace `tmt'=`tmt'-1
tempvar trt
qui:egen byte `trt'=group(`treatment') if `touse'
qui:replace `trt'=`trt'-1
** Sanity Check Weights
** Weights
if "`exp'"=="" {
tempname wgt
gen byte `wgt'=1
}
else {
tempvar wgt
qui:gen double `wgt'`exp'
qui:sum `wgt' if `touse' & `tmt'==0, meanonly
qui:replace `wgt'=`wgt'/r(mean)
}
**# Here we collect all options
** This are all the options send to DRDID
local 01 touse(`touse') tmt(`tmt') trt(`trt') y(`y') ///
xvar(`xvar') `isily' ivar(`ivar') ///
weight(`wgt') stub(`stub') ///
treatvar(`treatment') `rc1' cluster(`cluster') ///
`wboot' reps(`reps') wbtype(`wbtype') level(`level') binit(`binit')
*seed(`seed')
** Default will be IPT
if "`estimator'"!="all" {
if ("`semethod'"!="gmm") {
drdid_`estimator', `01' `diopts'
ereturn local semethod `semethod'
ereturn local clustvar `ocluster'
ereturn local seed `rseed'
if "`ivar'"!="" ereturn local datatype "panel"
else ereturn local datatype "rcs"
Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' level(`level')
exit
}
else {
tempvar touse2
tempname b V
if ("`ivar'"=="") {
local repeated repeated
if ("`estimator'"=="sipwra") {
display as error ///
"{bf:ipwra} not allowed with {bf:gmm} and repeated" ///
" cross-sectional data"
exit 198
}
}
if "`weight'" != "" {
local wgtgmm [`weight' `exp']
}
quietly generate byte `touse2' = . if `touse'
_het_did_gmm `y' `xvar' if `touse' `wgtgmm', ///
estimator(`estimator') groupvar(`ivar') ///
psvars(`xvar') treatvar(`trt') ///
timevar(`time') vce(`vce') touse2(`touse2') ///
t0(`tmt') pscore(`pscore') `repeated' ///
treatname(`treatment')
matrix `b' = r(b)
matrix `V' = r(V)
local vce "`r(vce)'"
local vcetype "`r(vcetype)'"
local N = e(N)
if ("`vce'"=="cluster") {
local N_clust = r(N_clust)
local clustvar "`r(clustvar)'"
}
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
ereturn local cmd drdid
ereturn local method `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
ereturn local semethod `semethod'
if ("`vce'"=="cluster") {
ereturn scalar N_clust = `N_clust'
ereturn local clustvar "`clustvar'"
}
ereturn local vce "`vce'"
ereturn local vcetype "`vcetype'"
ereturn local policy "`treatment'"
ereturn hidden local method2 ///
`drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
exit
}
}
if "`estimator'"=="all" {
** DR
tempvar ttouse
qui:clonevar `ttouse'=`touse'
tempname bb VV
qui:drdid_dripw , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname ATET:dripw
**capture drop `stub'att_dripw
*ren `stub'att `stub'att_dripw
if "`ivar'"=="" {
qui:capture gen byte `touse'=`ttouse'
qui:drdid_dripw , `01' rc1
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:dripw_rc1
**capture drop `stub'att_dripwrc
*clonevar `stub'att_dripwrc=`stub'att
}
** DRIMP
qui:capture gen byte `touse'=`ttouse'
qui:drdid_imp , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:drimp
**capture drop `stub'att_drimp
*ren `stub'att `stub'att_drimp
if "`ivar'"=="" {
qui:capture gen byte `touse'=`ttouse'
qui:drdid_imp , `01' rc1
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:drimp_rc1
**capture drop `stub'att_drimprc
*ren `stub'att `stub'att_drimprc
}
** REG
qui:capture gen byte `touse'=`ttouse'
qui:drdid_reg , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:reg
**capture drop `stub'att_reg
*ren `stub'att `stub'att_reg
** TRAD_IPW
qui:capture gen byte `touse'=`ttouse'
qui:drdid_aipw , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:ipw
**capture drop `stub'att_ipw
*ren `stub'att `stub'att_ipw
** STD IPW
qui:capture gen byte `touse'=`ttouse'
qui:drdid_stdipw , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
**capture drop `stub'att_stdipw
*ren `stub'att `stub'att_stdipw
local clname `clname' ATET:stdipw
if "`ivar'"!="" {
qui:capture gen byte `touse'=`ttouse'
qui:drdid_sipwra , `01'
matrix `bb' = nullmat(`bb')\e(b)
matrix `VV' = nullmat(`VV')\e(V)
local clname `clname' ATET:sipwra
}
matrix `bb'=`bb''
matrix colname `bb' = `clname'
matrix `VV'=diag(`VV')
matrix colname `VV' = `clname'
matrix rowname `VV' = `clname'
ereturn post `bb' `VV'
*ereturn display
}
ereturn local cmd drdid
ereturn local policy "`treatment'"
ereturn local method `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
ereturn hidden local method2 `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' level(`level')
end
program define _Vce_Parse, rclass
syntax [anything] [if] [in], [gmm WBOOT1 vce(string) WBOOT(string)]
marksample touse
local semethod standard
if (("`wboot'"!=""|"`wboot1'"!="") & "`vce'"!="") {
opts_exclusive "vce wboot"
}
if ("`gmm'"!="" & ("`wboot'"!=""|"`wboot1'"!="")) {
opts_exclusive "gmm wboot"
local semethod gmm
}
else if ("`gmm'"!="") {
local semethod "gmm"
}
else if (("`wboot'"!=""|"`wboot1'"!="")) {
if ("`wboot'"!="" & "`wboot1'"!="") {
display as error "incorrect wildbootstrap specification"
di as txt "{p 4 4 2}"
di as smcl as err ///
"You may specify {bf:wboot} or " ///
"{bf:wboot()} with arguments but not both."
di as smcl as err "{p_end}"
exit 198
}
local semethod "wildboot"
_Parse_Wildboot if `touse', `wboot1' `wboot'
local seed "`r(seed)'"
local reps = r(reps)
local wbtype = r(wbtype)
local cluster "`r(cluster)'"
if ("`cluster'"!="") {
local wncl = r(wncl)
return scalar wncl = r(wncl)
}
return local cluster "`cluster'"
return local seed "`seed'"
return local reps = `reps'
return local wbtype = `wbtype'
}
if ("`vce'"!="") {
_Vce_Parse_Clust, vce(`vce') `gmm'
return local cluster "`s(cluster)'"
}
return local semethod "`semethod'"
end
program _Parse_Wildboot, rclass
syntax [anything] [if] [in], [ ///
WBOOT1 ///
reps(integer 999) ///
rseed(string) ///
wbtype(string) ///
cluster(string) ///
]
marksample touse
if ("`wboot1'"=="") {
return local reps = `reps'
return local seed "`rseed'"
if ("`wbtype'"=="") {
local wbtypen = 1
}
else if ("`wbtype'"=="mammen") {
local wbtypen = 1
}
else if ("`wbtype'"=="rademacher") {
local wbtypen = 2
}
else if ("`wbtype'"!="rademacher" & "`wbtype'"!="mammen") {
display as error "invalid {bf:wbtype()}"
di as txt "{p 4 4 2}"
di as smcl as err ///
"{bf:wbtype()} should be one of {bf:mammen} or " ///
"{bf:rademacher}."
di as smcl as err "{p_end}"
exit 198
}
return local wbtype = `wbtypen'
if ("`cluster'"!="") {
tempvar nclust wncl0
capture confirm numeric variable `cluster'
local rc = _rc
if (`rc') {
capture destring `rest', generate(`nclust')
local rc = _rc
if (`rc') {
display in red "option {bf:cluster()} incorrectly specified"
exit 198
}
capture confirm numeric variable `nclust'
local rc = _rc
if (`rc') {
display in red "option {bf:cluster()} incorrectly specified"
exit 198
}
}
*quietly egen `wncl0' = group(`cluster') if `touse'
*summarize `wncl0', meanonly
*local wncl = r(max)
*return scalar wncl = `wncl'
}
return local cluster "`cluster'"
}
else {
return local reps = `reps'
return local seed "`rseed'"
return local wbtype = 1
}
end
program define _Vce_Parse_Clust, sclass
syntax [anything], [vce(string) gmm *]
gettoken key rest : vce, parse(", ")
local lkey = length(`"`key'"')
local nvce: list sizeof vce
local iscluster = 0
if (`"`key'"' == bsubstr("cluster",1,max(2,`lkey'))) {
local iscluster = 1
}
if (`nvce'>1 & `iscluster'==0 & "`gmm'"=="") {
display as error "{bf:vce()} option {bf:`key'} not allowed"
exit 198
}
if (`nvce'==1 & "`vce'"!="if" & "`gmm'"=="") {
display as error "{bf:vce()} option {bf:`key'} not allowed"
exit 198
}
if ("`vce'"=="if" & "`gmm'"!="") {
display as error "{bf:vce()} option {bf:if} not allowed"
exit 198
}
if (`nvce'>1) {
gettoken key rest : vce, parse(", ")
local lkey = length(`"`key'"')
local voy = 0
if `"`key'"' == bsubstr("cluster",1,max(2,`lkey')) {
capture confirm numeric variable `rest'
local rc = _rc
if (`rc') {
tempname nclust
capture destring `rest', generate(`nclust')
local rc = _rc
if (`rc') {
display in red "option {bf:vce()} incorrectly specified"
exit 198
}
capture confirm numeric variable `nclust'
local rc = _rc
if (`rc') {
display in red "option {bf:vce()} incorrectly specified"
exit 198
}
}
local voy = 1
}
if ("`key'"=="hac" & "`gmm'"!="") {
local nvce: list sizeof rest
local voy = 1
}
if ("`key'"=="hac" & "`gmm'"=="") {
display as error ///
"{bf:vce()} option {bf:`key'} not allowed with" ///
" estimator {bf:gmm}"
exit 198
}
if (`voy'==0) {
display in red "option {bf:vce()} incorrectly specified"
exit 198
}
}
if (`iscluster'==1) {
local cluster: word 2 of `vce'
sreturn local cluster "`cluster'"
}
end
program define _S_Me_thod, sclass
if ("`e(method)'"=="drimp") {
local tmodel "inverse probability tilting"
local omodel "weighted least squares"
}
if ("`e(method)'"=="aipw") {
local tmodel "inverse probability"
local omodel "weighted mean"
}
if ("`e(method)'"=="dripw") {
local tmodel "inverse probability"
local omodel "least squares"
}
if ("`e(method)'"=="reg") {
local tmodel "none"
local omodel "regression adjustment"
}
if ("`e(method)'"=="stdipw") {
local tmodel "stabilized inverse probability"
local omodel "weighted mean"
}
if ("`e(method)'"=="sipwra") {
local tmodel "inverse probability"
local omodel "regression adjustment"
}
sreturn local omodel "`omodel'"
sreturn local tmodel "`tmodel'"
end
program define Display
syntax [, bmatrix(passthru) vmatrix(passthru) COEFLegend *]
_get_diopts diopts rest, `options'
local myopts `bmatrix' `vmatrix'
if ("`rest'"!="") {
display in red "option {bf:`rest'} not allowed"
exit 198
}
_S_Me_thod
local omodel "`s(omodel)'"
local tmodel "`s(tmodel)'"
if ("`e(method)'"!="all") {
_coef_table_header, title(Doubly robust difference-in-differences)
noi display as text "Outcome model : {res:`omodel'}"
noi display as text "Treatment model: {res:`tmodel'}"
}
else {
_coef_table_header, ///
title(Doubly robust difference-in-differences estimator summary)
}
if ("`e(policy)'"!="" & "`e(semethod)'"=="wildboot") {
if "`e(clustvar)'"!="" {
display as text "(Std. err. adjusted for" ///
as result %9.0gc e(N_clust) ///
as text " clusters in " as result e(clustvar) as text ")"
}
_my_tab_drdid, `diopts'
}
if ("`e(semethod)'"!="wildboot") {
_coef_table, `diopts' `myopts' `coeflegend' neq(1)
}
if ("`e(method)'"=="all") {
local reg Outcome regression or Regression augmented estimator
display "{p}Note: This table is provided for comparison" ///
" across estimations only. You cannot use it to compare" ///
" estimates across different estimators{p_end}"
display "{cmd:dripw} :Doubly Robust IPW"
display "{cmd:drimp} :Doubly Robust Improved estimator"
display "{cmd:reg} :`reg'"
display "{cmd:ipw} :Abadie(2005) IPW estimator"
display "{cmd:stdipw}:Standardized IPW estimator"
display "{cmd:sipwra}:IPW and Regression adjustment estimator."
}
end
program _my_tab_drdid, rclass
syntax [, level(int `c(level)') noci cformat(string) sformat(string) *]
_get_diopts diopts rest, `options'
local cf %9.0g
local pf %5.3f
local sf %7.2f
if ("`cformat'"!="") {
local cf `cformat'
}
if ("`sformat'"!="") {
local sf `sformat'
}
tempname mytab z t ll ul cimat rtab
local policy "`e(policy)'"
local labn: value label `policy'
local largo = strlen("`policy'")
local newvs = 0
local novstab = 0
if ("`labn'"!="") {
local uno: label `labn' 1
local zero: label `labn' 0
local vs "(`uno' vs `zero')"
local widet = strlen("`vs'")
local widet0 = `widet'
local widet = max(`widet', `largo')
if (`largo'>`widet0' & `widet'>16) {
local policy = abbrev("`policy'", 20)
local widet = 21
}
if (`widet0'>24) {
local kn1: list sizeof uno
local kn0: list sizeof zero
local kvs0 = max(max(max(`kn1', `kn0'), `largo'), 13)
local policy = abbrev("`policy'", `kvs0')
local uno0 = abbrev("`uno'", `kvs0')
local zero0 = abbrev("`zero'", `kvs0')
local vs1 "(`uno0' "
local vs2 "vs"
local vs3 "`zero0')"
local newvs = 1
local widet = `kvs0' + 2
}
}
else {
local vs "(1 vs 0)"
local widet = 13
local widet0 = `widet'
local widet = max(`widet', `largo')
if (`largo'>`widet0' & `widet'>21) {
local policy = abbrev("`policy'", 20)
local widet = 21
}
}
.`mytab' = ._tab.new, col(6) lmargin(0)
.`mytab'.width `widet' |12 12 8 12 12
.`mytab'.titlefmt . . . %6s %24s .
.`mytab'.pad . 2 1 0 3 3
.`mytab'.numfmt . `cf' `cf' `sf' `cf' `cf'
local stat t
/*if "`e(df_r)'" != "" {
local stat t
scalar `z' = invttail(e(df_r),(100-`level')/200)
}
else {
local stat z
scalar `z' = invnorm((100+`level')/200)
}*/
local namelist : colname e(b)
local eqlist : coleq e(b)
local k : word count `namelist'
local name : word 1 of `namelist'
matrix `cimat'= e(ciband)
scalar `ll' = `cimat'[1,1]
scalar `ul' = `cimat'[1,2]
scalar `t' = `beq'_b[`name']/`beq'_se[`name']
matrix `rtab' = _b[`name'] \ ///
_se[`name'] \ ///
`t' \ ///
. \ ///
`ll' \ ///
`ul' \ ///
. \ ///
. \ ///
0
matrix colnames `rtab' = ATET:r1vs0.`policy'
matrix rownames `rtab' = b se t pvalue ll ul df crit eform
.`mytab'.sep, top
if `:word count `e(depvar)'' == 1 {
local depvar "`e(depvar)'"
local depvar = abbrev("`depvar'", 12)
}
.`mytab'.titles "`depvar'" /// 1
" Coefficient" /// 2
"Std. err." /// 3
"`stat'" /// 4
"[`level'% conf. interval]" "" // 5 6
local eq : word 1 of `eqlist'
forvalues i = 1/3 {
if "`eq'" != "_" {
if "`eq'" != "`eq0'" {
.`mytab'.strcolor result . . . . .
.`mytab'.strfmt %-12s . . . . .
if (`i'==1) {
.`mytab'.sep
.`mytab'.row "`eq'" "" "" "" "" ""
}
.`mytab'.strcolor text . . . . .
.`mytab'.strfmt %12s . . . . .
}
local beq "[`eq']"
}
if (`i'==3) {
.`mytab'.row "`vs'" ///
`beq'_b[`name'] ///
`beq'_se[`name'] ///
`t' ///
`ll' `ul'
}
if (`i'==2) {
.`mytab'.row "`policy'" "" "" "" "" ""
}
}
.`mytab'.sep, bottom
return matrix table = `rtab'
end
program define easter_egg
display "{p}This is just for fun. Its my attempt to an Easter Egg within my program. {p_end}" _n ///
"{p} Also, if you are reading this, it means you are lucky," ///
"only 0.1% of people using this program will see this message. {p_end}" _n ///
"{p} This program was inspired by challenge post by Scott Cunningham. Unforunately, I arrived late due to " ///
"lack of understanding. I tried to do CSDID (RUN) before learning DRDID (walk). {p_end} " _n ///
"{p} Asjad, was the first person who started to properly implement the code in Stata. I got inspired on his work that, " ///
"Reread the paper, and voala. Everything fit in the form of the first version of this code. {p_end} " _n ///
"{p} Many thanks go to Pedro, who spend a lot of time explaining details that would otherwise would remain confusing (he is the author fo the original paper after all; {p_end} " _n ///
"{p} to Miklos, who pushed us to mantain a Github repository for this program. He is taking the side of the user {p_end}" ///
"{p} To Enrique Pinzon, who helped from the shadows, for a smooth transition from R to Stata (Yay Stata) {p_end}" _n ///
"{p} and Austin, who was working in parallel when Asjad and I took the lead on this. {p_end}"
end
////////////////////////////////////////////////////////////////////////////////
*** FIRST
**# Abadies
program define drdid_aipw, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
stub(name) ///
treatvar(string) ///
wboot ///
reps(int 999) ///
level(int 95) ///
wbtype(int 1) ///
seed(string) ///
cluster(str) ///
binit(str) ///
* ///
]
** PS
*set trace on
_get_diopts diopts other, `options'
quietly capture Display, `diopts' `other'
if _rc==198 {
Display, `diopts' `other'
}
tempvar att psxb __dy__ xb touse2
tempname psb psV
qui:gen double `att'=.
if "`ivar'"!="" {
*display "Estimating IPW logit"
qui {
`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
predict double `psxb', xb
matrix `psb'=e(b)
matrix `psV'=e(V)
** _delta
bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
gen byte `touse2'=`touse'*(`tmt'==0)
replace `touse'=0 if `__dy__'==.
tempname b V ciband ncl
mata:ipw_abadie_panel("`__dy__'","`xvar' ","`xb'","`psb' ","`psV' ","`psxb'","`trt'","`tmt'","`touse2'","`att'","`weight'")
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b'= ATET:r1vs0.`treatvar'
matrix colname `V'= ATET:r1vs0.`treatvar'
matrix rowname `V'= ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
}
*display "DiD with IPW"
*display "Abadie (2005) inverse probability weighting DiD estimator"
*ereturn display
}
else if "`ivar'"=="" {
qui {
`isily' logit `trt' `xvar' if `touse' [iw = `weight']
tempvar psxb
predict double `psxb', xb
tempname psb psV
matrix `psb'=e(b)
matrix `psV'=e(V)
tempname b V ciband ncl
mata:ipw_abadie_rc("`y'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
** Wbootstrap Multipler
local ci = `level'/100
local touse2 `touse'
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
}
quietly count if `touse'
local N = r(N)
tempvar touse2
clonevar `touse2'=`touse'
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
*ereturn display
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
}
** if STUB is used, then RIF is saved
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
ereturn local vcetype "Wboot"
}
ereturn local cmd drdid
ereturn local method aipw
ereturn hidden local method2 aipw
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
** can be more efficient
**#drdid_dripw
program define drdid_dripw, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
rc1 ///
stub(name) ///
treatvar(string) ///
wboot ///
level(int 95) ///
reps(int 999) ///
wbtype(int 1) ///
seed(string) ///
cluster(str) ///
binit(str) ///
* ///
]
** PS
tempvar att psxb __dy__ xb touse2
tempname psb psV regb regV b V ciband ncl
qui:gen double `att'=.
if "`ivar'"!="" {
*display "Estimating IPW logit"
qui {
`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
predict double `psxb', xb
matrix `psb'=e(b)
matrix `psV'=e(V)
** _delta
bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
** Reg for outcome
`isily' reg `__dy__' `xvar' if `touse' & `trt'==0 & `tmt'==0 [iw = `weight']
matrix `regb'=e(b)
matrix `regV'=e(V)
predict double `xb'
*capture drop `stub'att
*gen double `stub'att=.
gen byte `touse2'=`touse'*(`tmt'==0)
replace `touse'=0 if `__dy__'==.
mata:drdid_panel("`__dy__'","`xvar' ","`xb'","`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'","`att'","`weight'")
**replace `stub'att=. if `tmt'==1
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b'= ATET:r1vs0.`treatvar'
matrix colname `V'= ATET:r1vs0.`treatvar'
matrix rowname `V'= ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
ereturn matrix regb `regb'
ereturn matrix regV `regV'
}
/* display "DR DiD with IPW and OLS"
display "Sant'Anna and Zhao (2020)" _n "{p}Doubly robust DiD estimator based on stabilized inverse probability weighting and ordinary least squares{p_end}"
ereturn display*/
}
else if "`ivar'"=="" {
qui {
`isily' logit `trt' `xvar' if `touse' [iw = `weight']
tempvar psxb
predict double `psxb', xb
tempname psb psV
matrix `psb'=e(b)
matrix `psV'=e(V)
tempname b V ciband ncl
*capture drop `stub'att
*gen double `stub'att=.
**ols
tempvar y00 y01 y10 y11
tempname regb00 regb01 regb10 regb11
tempname regV00 regV01 regV10 regV11
`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==0 [iw = `weight']
predict double `y00'
matrix `regb00'=e(b)
matrix `regV00'=e(V)
`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==1 [iw = `weight']
predict double `y01'
matrix `regb01'=e(b)
matrix `regV01'=e(V)
`isily' reg `y' `xvar' if `trt'==1 & `tmt' ==0 [iw = `weight']
predict double `y10'
matrix `regb10'=e(b)
matrix `regV10'=e(V)
`isily' reg `y' `xvar' if `trt'==1 & `tmt' ==1 [iw = `weight']
predict double `y11'
matrix `regb11'=e(b)
matrix `regV11'=e(V)
if "`rc1'"=="" {
mata:drdid_rc("`y'","`y00' `y01' `y10' `y11'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
}
else {
mata:drdid_rc1("`y'","`y00' `y01' `y10' `y11'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
local nle "Not Locally efficient"
}
////
local touse2 `touse'
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
}
*display "DR DiD with IPW and OLS for Repeated Crossection: `nle'"
*display "Sant'Anna and Zhao (2020)" _n "{p}Doubly robust DiD estimator based on stabilized inverse probability weighting and ordinary least squares{p_end}"
quietly count if `touse'
local N = r(N)
tempvar touse2
clonevar `touse2'=`touse'
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
*ereturn post `b' `V'
*ereturn display
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
ereturn matrix regb00 `regb00'
ereturn matrix regV00 `regV00'
ereturn matrix regb01 `regb01'
ereturn matrix regV01 `regV01'
ereturn matrix regb10 `regb10'
ereturn matrix regV10 `regV10'
ereturn matrix regb11 `regb11'
ereturn matrix regV11 `regV11'
}
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
}
ereturn local cmd drdid
ereturn local method dripw
ereturn hidden local method2 dripw
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
**#drdid_reg
program define drdid_reg, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
stub(name) ///
treatvar(string) ///
wboot ///
level(int 95) ///
reps(int 999) ///
wbtype(int 1) ///
seed(string) ///
cluster(str) ///
binit(str) ///
* ///
]
** Simple application. But right now without RIF
tempvar att
qui:gen double `att'=.
if "`ivar'"!="" {
qui {
tempvar __dy__ xb touse2
tempname regb regV b V ciband ncl
bysort `touse' `ivar' (`tmt'):gen double ///
`__dy__'=`y'[2]-`y'[1] if `touse'
`isily' reg `__dy__' `xvar' if `touse' & `trt'==0 ///
& `tmt'==0 [iw = `weight']
predict double `xb'
//////////////////////
matrix `regb'=e(b)
matrix `regV'=e(V)
gen byte `touse2'=`touse'*(`tmt'==0)
replace `touse'=0 if `__dy__'==.
mata:reg_panel("`__dy__'", "`xvar' ", "`xb' " , "`trt'", ///
"`tmt'" , "`touse2'","`att'","`weight'")
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix regb `regb'
ereturn matrix regV `regV'
}
*display "DiD with OR" _n "Outcome regression DiD estimator based on ordinary least squares"
*ereturn display
}
else if "`ivar'"=="" {
qui {
tempvar y00 y01
tempname regb00 regb01 regV00 regV01
`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==0 [iw = `weight']
predict double `y00'
matrix `regb00' = e(b)
matrix `regV00' = e(V)
`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==1 [iw = `weight']
predict double `y01'
matrix `regb01' = e(b)
matrix `regV01' = e(V)
tempname b V ciband ncl
*capture drop `stub'att
*gen double `stub'att=.
mata:reg_rc("`y'","`y00' `y01'","`xvar' ", ///
"`tmt'","`trt'","`weight'","`touse'","`att'")
local ci = `level'/100
local touse2 `touse'
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
tempvar touse2
clonevar `touse2'=`touse'
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
ereturn matrix regb00 `regb00'
ereturn matrix regV00 `regV00'
ereturn matrix regb01 `regb01'
ereturn matrix regV01 `regV01'
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
}
*display "DiD with OR for RC" _n "Outcome regression DiD estimator based on ordinary least squares"
*ereturn display
}
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
}
ereturn local cmd drdid
ereturn local method reg
ereturn hidden local method2 reg
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
**#drdid_sipw
program define drdid_stdipw, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
stub(name) ///
treatvar(string) ///
wboot ///
level(int 95) ///
reps(int 999) ///
wbtype(int 1) ///
seed(string) ///
cluster(str) ///
binit(str) ///
* ///
]
tempvar att
qui:gen double `att'=.
** Simple application. But right now without RIF
if "`ivar'"!="" {
*display "Estimating IPW logit"
qui {
`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
tempvar psxb __dy__ xb touse2
tempname psb psV b V ciband ncl
predict double `psxb', xb
matrix `psb'=e(b)
matrix `psV'=e(V)
** _delta
bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] ///
if `touse'
** Reg for outcome
}
*display "Estimating Counterfactual Outcome"
qui {
*`isily' reg `__dy__' `xvar' if `trt'==0
*gen double `stub'att=.
gen byte `touse2'=`touse'*(`tmt'==0)
replace `touse'=0 if `__dy__'==.
mata:std_ipw_panel("`__dy__'","`xvar' ","`xb'", ///
"`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'", ///
"`att'","`weight'")
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b'= ATET:r1vs0.`treatvar'
matrix colname `V'= ATET:r1vs0.`treatvar'
matrix rowname `V'= ATET:r1vs0.`treatvar'
}
*display "DiD with stabilized IPW" _n "{p}Abadie (2005) inverse probability weighting DiD estimator with stabilized weights{p_end}"
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
*ereturn display
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
}
else if "`ivar'"=="" {
qui {
`isily' logit `trt' `xvar' if `touse' [iw = `weight']
tempvar psxb
predict double `psxb', xb
tempname psb psV
matrix `psb'=e(b)
matrix `psV'=e(V)
tempname b V ciband ncl
*capture drop `stub'att
*gen `stub'att=.
mata:std_ipw_rc("`y'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
local ci = `level'/100
local touse2 `touse'
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
}
quietly count if `touse'
local N = r(N)
tempvar touse2
clonevar `touse2'=`touse'
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix ipwb `psb'
ereturn matrix ipwV `psV'
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
}
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
}
ereturn local cmd drdid
ereturn local method stdipw
ereturn hidden local method2 stdipw
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
// only one without Mata writting. Consider working on it
**#drdid_sipwra
program define drdid_sipwra, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
stub(name) ///
treatvar(string) ///
reps(int 999) /// Notused here
wbtype(int 1) /// not used here
seed(string) ///
level(int 95) /// not used here
cluster(str) ///
* ///
]
tempvar att
qui:gen double `att'=.
** Simple application. But right now without RIF
if "`cluster'"!="" {
local clopt vce(cluster `cluster')
}
if "`ivar'"=="" {
display "Estimator not implemented for RC"
exit 198
}
else {
qui {
tempvar __dy__ sy
bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
sum `__dy__' if `touse'
local scl = r(mean)
gen double `sy'= `__dy__'/`scl'
qui:teffects ipwra (`sy' `xvar') (`trt' `xvar', logit) ///
if `touse' & `tmt'==0 [iw = `weight'] , atet `clopt' iter(5)
tempname b V ciband ncl aux
matrix `aux'=e(b)*`scl'
matrix `b'=`aux'[1,1]
matrix `aux'=e(V)*`scl'^2
matrix `V'=`aux'[1,1]
matrix colname `b' = ATET:r1vs0.`treatvar'
matrix colname `V' = ATET:r1vs0.`treatvar'
matrix rowname `V' = ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
}
*ereturn display
}
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
ereturn local cmd drdid
ereturn local method sipwra
ereturn hidden local method2 sipwra
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
**#drdid_dript
program define drdid_imp, eclass
syntax, [ ///
touse(str) ///
trt(str) ///
y(str) ///
xvar(str) ///
noisily ///
ivar(str) ///
tmt(str) ///
weight(str) ///
rc1 stub(name) ///
treatvar(string) ///
wboot ///
level(int 95) ///
reps(int 999) ///
wbtype(int 1) ///
seed(string) ///
cluster(str) ///
binit(str) ///
* ///
]
_get_diopts diopts other, `options'
quietly capture Display, `diopts' `other'
if _rc==198 {
Display, `diopts' `other'
}
tempvar att psxb __dy__ w0 xb
tempname iptb iptV regb regV b V ciband ncl
qui:gen double `att'=.
if "`ivar'"!="" {
qui {
`isily' mlexp (`trt'*{xb:`xvar' _cons}-(`trt'==0)*exp({xb:})) ///
if `touse' & `tmt'==0 [iw = `weight'], from(`binit') ///
derivative(/xb=`trt'-(`trt'==0)*exp({xb:}))
*vce(robust)
matrix `iptb'=e(b)
matrix `iptV'=e(V)
predict double `psxb',xb
** Determine dy and dyhat
bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
** determine weights
gen double `w0' = ((logistic(`psxb')*(1-`trt')))/(1-logistic(`psxb'))*`weight'
sum `w0' if `touse' , meanonly
replace `w0'=`w0'/r(mean)
** estimating dy_hat for a counterfactual
`isily' reg `__dy__' `xvar' [w=`w0'] if `trt'==0 & `tmt'==0,
matrix `regb' =e(b)
matrix `regV' =e(V)
qui:predict double `xb'
tempvar touse2
gen byte `touse2'=`touse'*(`tmt'==0)
replace `touse'=0 if `__dy__'==.
mata:drdid_imp_panel("`__dy__'","`xvar' ","`xb'", ///
"`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'", ///
"`att'","`weight'")
********************************************************************************
********************************************************************************
local ci = `level'/100
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b'= ATET:r1vs0.`treatvar'
matrix colname `V'= ATET:r1vs0.`treatvar'
matrix rowname `V'= ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix iptb `iptb'
ereturn matrix iptV `iptV'
ereturn matrix regb `regb'
ereturn matrix regV `regV'
ereturn local cmd drdid
ereturn local method drimp
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
}
}
*display "DR DiD with IPT and WLS" _n "{p}Sant'Anna and Zhao (2020) Improved doubly robust DiD estimator based on inverse probability of tilting and weighted least squares{p_end}"
*ereturn display
}
else {
**# for Crossection estimator
qui {
/*`isily' gmm ((`trt'==1)-(`trt'==0)*exp({b:`xvar' _cons})) if `touse' [iw = `weight'], ///
instrument(`xvar' ) derivative(/b=-(`trt'==0)*exp({b:})) ///
onestep winit(identity) */
`isily' mlexp (`trt'*{xb:`xvar' _cons}-(`trt'==0)*exp({xb:})) ///
if `touse' [iw = `weight'], vce(robust) from(`binit') ///
derivative(/xb=`trt'-(`trt'==0)*exp({xb:}))
//& `tmt'==0
tempname iptb iptV regb00 regV00 regb01 regV01 regb10 ///
regV10 regb11 regV11
tempvar psxb w1 w0 y01 y00 y10 y11
matrix `iptb'=e(b)
matrix `iptV'=e(V)
predict double `psxb', xb
** outcomes
gen double `w0' = (1-`trt')*logistic(`psxb')/(1-logistic(`psxb'))
`isily' reg `y' `xvar' [w=`w0'] if `trt'==0 & `tmt'==0,
predict double `y00'
matrix `regb00' =e(b)
matrix `regV00' =e(V)
`isily' reg `y' `xvar' [w=`w0'] if `trt'==0 & `tmt'==1,
predict double `y01'
matrix `regb01' =e(b)
matrix `regV01' =e(V)
`isily' reg `y' `xvar' if `trt'==1 & `tmt'==0,
predict double `y10'
matrix `regb10' =e(b)
matrix `regV10' =e(V)
`isily' reg `y' `xvar' if `trt'==1 & `tmt'==1,
predict double `y11'
matrix `regb11' =e(b)
matrix `regV11' =e(V)
tempname b V ciband ncl
if "`rc1'"=="" {
mata:drdid_imp_rc("`y'","`y00' `y01' `y10' `y11'", ///
"`xvar' ","`tmt'","`trt'","`iptV'","`psxb'", ///
"`weight'","`touse'","`att'")
}
else {
mata:drdid_imp_rc1("`y'","`y00' `y01' `y10' `y11'", ///
"`xvar' ","`tmt'","`trt'","`iptV'","`psxb'", ///
"`weight'","`touse'","`att'")
local nle "Not Locally efficient"
}
local ci = `level'/100
local touse2 `touse'
mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
matrix colname `b'=ATET:r1vs0.`treatvar'
matrix colname `V'=ATET:r1vs0.`treatvar'
matrix rowname `V'=ATET:r1vs0.`treatvar'
quietly count if `touse'
local N = r(N)
tempvar touse2
clonevar `touse2'=`touse'
ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
ereturn local cmd drdid
ereturn local method drimp
if ("`wboot'"=="wboot") {
ereturn matrix ciband = `ciband'
ereturn local semethod "wildboot"
}
}
*display "{p}DR DiD with IPT and WLS for OLS `nle'{p_end}" _n "{p}Sant'Anna and Zhao (2020) Improved doubly robust DiD estimator based on inverse probability of tilting and weighted least squares{p_end}"
*ereturn display
local att1 =`=_b[r1vs0.`treatvar']'
local attvar1 =`=_se[r1vs0.`treatvar']'^2
ereturn scalar att1 =`att1'
ereturn scalar attvar1 =`attvar1'
ereturn matrix iptb `iptb'
ereturn matrix iptV `iptV'
ereturn matrix regb00 `regb00'
ereturn matrix regV00 `regV00'
ereturn matrix regb01 `regb01'
ereturn matrix regV01 `regV01'
ereturn matrix regb10 `regb10'
ereturn matrix regV10 `regV10'
ereturn matrix regb11 `regb11'
ereturn matrix regV11 `regV11'
}
if "`stub'"!="" {
qui:capture drop `stub'att
qui:gen double `stub'att=`att'
}
if "`cluster'"!="" {
ereturn scalar N_clust =`=scalar(`ncl')'
ereturn local clustvar `cluster'
}
ereturn local policy "`treatvar'"
*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts'
end
program define _S__est, sclass
syntax [anything], [drimp dripw reg stdipw ipw ipwra all * ]
if ("`drimp'"!="") {
local drimp imp
}
if ("`ipw'"!="") {
local ipw "aipw"
}
if ("`ipwra'"!="") {
local ipwra "sipwra"
}
sreturn local estimator ///
"`drimp'`dripw'`reg'`stdipw'`ipw'`ipwra'`all'"
end
mata
////////////////////////////////////////////////////////////////////////////////////////////////
// dript
void drdid_imp_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
// This code is based on Asjad Replication
// Gather all data
real scalar nn
dy =st_data(.,dy_ ,touse)
nn=rows(dy)
//st_view(xvar=.,xvar_,touse)
xb=st_data(.,xb_,touse)
psxb=st_data(.,psxb_,touse)
real matrix psc
psc=logistic(psxb)
trt =st_data(.,trt_ ,touse)
real matrix w
w=st_data(.,ww,touse)
real matrix w_1 , w_0, att, att_inf_func
w_1 = w :* trt
w_0 = w :* psc :* (1:-trt):/(1:-psc)
w_1 = w_1:/mean(w_1)
w_0 = w_0:/mean(w_0)
att=(dy:-xb):*(w_1:-w_0)
att_inf_func = mean(att) :+ att :- w_1:*mean(att)
st_store(.,rif,touse, att_inf_func)
// Variance should be divided by n not n-1
}
// standard IPW
void std_ipw_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
// Gather all data
real scalar nn
dy =st_data(.,dy_ ,touse)
nn=rows(dy)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(dy),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
//xb=st_data(.,xb_,touse)
psxb=st_data(.,psxb_,touse)
real matrix psc
psc=logistic(psxb)
trt =st_data(.,trt_ ,touse)
tmt =st_data(.,tmt_ ,touse)
// and matrices
//psb =st_matrix(psb_ )
psv =st_matrix(psV_ )
// for now assume weights = 1
real matrix w
w=st_data(.,ww,touse)
real matrix w_1, w_0, att_cont, att_treat,
eta_treat, eta_cont,
lin_ps, att_inf_func
w_1= w :* trt
w_0= w :* psc :* (1 :- trt):/(1 :- psc)
att_treat = w_1:* dy
att_cont = w_0:* dy
eta_treat = mean(att_treat)/mean(w_1)
eta_cont = mean(att_cont)/mean(w_0)
ipw_att = eta_treat :- eta_cont
inf_treat = (att_treat :- w_1 :* eta_treat)/mean(w_1)
inf_cont_1 = (att_cont :- (w_0 :* eta_cont))
lin_ps = (w:* (trt :- psc) :* xvar)*(psv * nn)
//M2 =
inf_cont_2 = lin_ps * mean(w_0 :* (dy :- eta_cont) :* xvar)'
inf_control = (inf_cont_1 :+ inf_cont_2)/mean(w_0)
att_inf_func = mean(ipw_att):+inf_treat :- inf_control
st_store(.,rif,touse, att_inf_func)
}
void reg_panel(string scalar dy_, xvar_, xb_ , trt_,tmt_,touse,rif,ww) {
real matrix dy, xvar, xb, trt, tmt
// This code is based on Asjad Replication
// Gather all data
real scalar nn
dy =st_data(.,dy_ ,touse)
nn=rows(dy)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(dy),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
xb=st_data(.,xb_,touse)
// psxb=st_data(.,psxb_,touse)
// real matrix psc
// psc=logistic(psxb)
trt =st_data(.,trt_ ,touse)
tmt =st_data(.,tmt_ ,touse)
// and matrices
// psb =st_matrix(psb_ )
// psv =st_matrix(psV_ )
// for now assume weights = 1
real matrix w
w=st_data(.,ww,touse)
real matrix w_1, w_0, att_cont, att_treat,
eta_treat, eta_cont, wols_x, wols_eX,
lin_ols, att_inf_func
w_1 = w :* trt
w_0 = w :* trt
att_treat = w_1:* dy
att_cont = w_0:* xb
eta_treat = mean(att_treat):/mean(w_1)
eta_cont = mean(att_cont) :/mean(w_0)
reg_att = eta_treat :- eta_cont
wols = w :* (1 :- trt)
wols_x = wols :* xvar
wols_eX = wols :* (dy:-xb) :* xvar
XpX_inv = invsym(quadcross(wols_x, xvar))*nn
lin_ols = wols_eX * XpX_inv
inf_treat = (att_treat :- w_1 * eta_treat):/mean(w_1)
inf_cont_1 = (att_cont :- w_0 * eta_cont)
inf_cont_2 = lin_ols * mean(w_0 :* xvar )'
inf_control = (inf_cont_1 :+ inf_cont_2):/mean(w_0)
att_inf_func = mean(reg_att):+(inf_treat :- inf_control)
st_store(.,rif,touse, att_inf_func)
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// TIPW drdid_tipw
void ipw_abadie_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
// Gather all data
real scalar nn
dy =st_data(.,dy_ ,touse)
nn=rows(dy)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(dy),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
//xb=st_data(.,xb_,touse)
psxb=st_data(.,psxb_,touse)
real matrix psc
psc=logistic(psxb)
trt =st_data(.,trt_ ,touse)
tmt =st_data(.,tmt_ ,touse)
// and matrices
psb =st_matrix(psb_ )
psv =st_matrix(psV_ )
// for now assume weights = 1
real matrix w
w=st_data(.,ww,touse)
real matrix w_1, w_0, att_cont, att_treat,
eta_treat, eta_cont, ipw_att,
lin_ps, att_lin1, mom_logit, att_lin2, att_inf_func
w_1= w :* trt
w_0= w :* psc :* (1 :- trt):/(1 :- psc)
att_treat = w_1:* dy
att_cont = w_0:* dy
eta_treat = mean(att_treat)/mean(w_1)
eta_cont = mean(att_cont)/mean(w_1)
ipw_att = eta_treat - eta_cont
lin_ps = (w:* (trt :- psc) :* xvar)*(psv * nn)
att_lin1 = att_treat :- att_cont
mom_logit = mean(att_cont :* xvar)
att_lin2 = lin_ps * mom_logit'
att_inf_func = mean(ipw_att):+(att_lin1 :- att_lin2 :- w_1 :* ipw_att)/mean(w_1)
st_store(.,rif,touse, att_inf_func)
// Variance should be divided by n not n-1
}
/// drdid_ipw
void drdid_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
// This code is based on Asjad Replication
// Gather all data
real scalar nn
dy =st_data(.,dy_ ,touse)
nn=rows(dy)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(dy),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
//xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
xb=st_data(.,xb_,touse)
psxb=st_data(.,psxb_,touse)
real matrix psc
psc=logistic(psxb)
trt =st_data(.,trt_ ,touse)
tmt =st_data(.,tmt_ ,touse)
// and matrices
psb =st_matrix(psb_ )
psv =st_matrix(psV_ )
// for now assume weights = 1
real matrix w
w=st_data(.,ww,touse)
// TRAD DRDID
real matrix w_1, w_0
w_1 = (w:*trt)
w_1 = w_1 /mean(w_1)
w_0 = (w:*psc:*(-trt:+1):/(-psc:+1))
w_0 = w_0 /mean(w_0 )
// ATT
real matrix dy_xb, att, w_ols, wols_eX,
XpX_inv, lin_wols, lin_ps,
n1, n0, nest, a, att_inf_func
dy_xb=dy:-xb
att = mean((w_1:-w_0):*(dy_xb))
// influence functions OLS
w_ols = w :* (1 :- trt)
wols_eX = w_ols:* (dy_xb):* xvar
XpX_inv = invsym(quadcross(xvar,w_ols,xvar))*nn
lin_wols = wols_eX * XpX_inv
// IF for logit
lin_ps = (w :* (trt:-psc) :* xvar) * (psv * nn)
// Components for RIF
n1 = w_1:*((dy_xb):-mean(dy_xb,w_1))
n0 = w_0:*((dy_xb):-mean(dy_xb,w_0))
a = ((1:-trt):/(1:-psc):^2)/ mean(psc:*(1:-trt):/(1:-psc))
// This only works because w_1 and w_0 are mutually exclusive
nest = lin_wols * (mean(xvar,w_1):-mean(xvar,w_0))' :+
lin_ps * mean( a :* (dy_xb :- mean(dy_xb,w_0)) :* exp(psxb):/(1:+exp(psxb)):^2:*xvar)'
// RIF att_inf_func = inf_treat' :- inf_control
att_inf_func = att:+n1:-n0:-nest
st_store(.,rif,touse, att_inf_func)
// Variance should be divided by n not n-1
}
//// standard IPW
//# RC
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void std_ipw_rc(string scalar y_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real scalar nn
nn = rows(y)
real matrix w10, w11, w00, w01
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
w00 = w00:/mean(w00 )
w01 = w01:/mean(w01 )
w10 = w10:/mean(w10 )
w11 = w11:/mean(w11 )
real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post,
eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post
att_treat_pre = w10 :* y
att_treat_post = w11 :* y
att_cont_pre = w00 :* y
att_cont_post = w01 :* y
eta_treat_pre = mean(att_treat_pre)
eta_treat_post = mean(att_treat_post)
eta_cont_pre = mean(att_cont_pre)
eta_cont_post = mean(att_cont_post)
real matrix ipw_att, lin_rep_ps
ipw_att = (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre)
//score_ps = wgt :* (trt :- psc) :* xvar
//Hessian_ps = psv :* nn
lin_rep_ps = (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
real matrix inf_treat, inf_cont, M2_pre, M2_post, inf_cont_ps, att_inf_func
inf_treat = (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre :- w10 :* eta_treat_pre)
inf_cont = (att_cont_post :- w01 :* eta_cont_post) :- (att_cont_pre :- w00 :* eta_cont_pre)
M2_pre = mean(w00 :* (y :- eta_cont_pre) :* xvar)
M2_post = mean(w01 :* (y :- eta_cont_post) :* xvar)
inf_cont_ps = lin_rep_ps * (M2_post :- M2_pre)'
inf_cont = inf_cont :+ inf_cont_ps
att_inf_func = ipw_att :+ inf_treat :- inf_cont
st_store(.,rif,touse,att_inf_func)
// -15.80330618
// 9.087929526
}
// Regression approach
void reg_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, wgt_ , touse, rif) {
// main Loading variables
real matrix y, y00, y01,
xvar, tmt, trt, wgt
y = st_data(.,y_, touse)
yy = st_data(.,yy_, touse)
y00 = yy[,1]
y01 = yy[,2]
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
//psc = logistic(st_data(.,pxb_, touse))
//psv = st_matrix(psv_)
wgt = st_data(.,wgt_, touse)
real scalar nn
nn = rows(y)
real matrix w10, w11, w0
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w0 = wgt :* trt
w10 = w10:/mean(w10 )
w11 = w11:/mean(w11 )
w0 = w0:/mean(w0 )
real matrix att_treat_pre, att_treat_post, att_cont,
eta_treat_pre, eta_treat_post, eta_cont, reg_att
att_treat_pre = w10 :* y
att_treat_post = w11 :* y
att_cont = w0 :* (y01 :- y00)
eta_treat_pre = mean(att_treat_pre)
eta_treat_post = mean(att_treat_post)
eta_cont = mean(att_cont)
reg_att = (eta_treat_post :- eta_treat_pre) :- eta_cont
real matrix w_ols_pre, wols_eX_pre, XpX_inv_pre, lin_rep_ols_pre
w_ols_pre = wgt :* (1 :- trt) :* (1 :- tmt)
//wols_x_pre = w_ols_pre :* xvar
wols_eX_pre = w_ols_pre :* (y :- y00) :* xvar
XpX_inv_pre = invsym(quadcross(xvar,w_ols_pre, xvar)):*nn
lin_rep_ols_pre = wols_eX_pre * XpX_inv_pre
real matrix w_ols_post, wols_eX_post, XpX_inv_post, lin_rep_ols_post
w_ols_post = wgt :* (1 :- trt) :* tmt
//wols_x_post = w_ols_post :* xvar
wols_eX_post = w_ols_post :* (y :- y01) :* xvar
XpX_inv_post = invsym(quadcross(xvar, w_ols_post, xvar)):*nn
lin_rep_ols_post = wols_eX_post * XpX_inv_post
real matrix inf_treat, inf_cont_1, inf_cont_2_post, inf_cont_2_pre, inf_control
//inf_treat_pre = (att_treat_pre :- w10 :* eta_treat_pre)
//inf_treat_post = (att_treat_post :- w11 :* eta_treat_post)
inf_treat = (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre :- w10 :* eta_treat_pre)
inf_cont_1 = (att_cont :- w0 :* eta_cont)
//M1 = mean(w0 :* xvar)
inf_cont_2_post = lin_rep_ols_post * mean(w0 :* xvar)'
inf_cont_2_pre = lin_rep_ols_pre * mean(w0 :* xvar)'
inf_control = (inf_cont_1 :+ inf_cont_2_post :- inf_cont_2_pre)
real matrix att_inf_func
att_inf_func = reg_att :+ (inf_treat :- inf_control)
st_store(.,rif,touse,att_inf_func)
// Variance should be divided by n not n-1
}
/// Abadies IPW
void ipw_abadie_rc(string scalar y_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real scalar nn
nn = rows(y)
real matrix w10, w11, w00, w01
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
real matrix pi_hat, lambda_hat, one_lambda_hat
pi_hat = mean(wgt :* trt)
lambda_hat = mean(wgt :* tmt)
one_lambda_hat = mean(wgt :* (1 :- tmt))
real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post,
eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post
att_treat_pre = w10 :* y:/(pi_hat :* one_lambda_hat)
att_treat_post = w11 :* y:/(pi_hat :* lambda_hat)
att_cont_pre = w00 :* y:/(pi_hat :* one_lambda_hat)
att_cont_post = w01 :* y:/(pi_hat :* lambda_hat)
eta_treat_pre = mean(att_treat_pre)
eta_treat_post = mean(att_treat_post)
eta_cont_pre = mean(att_cont_pre)
eta_cont_post = mean(att_cont_post)
real matrix ipw_att
ipw_att = (eta_treat_post - eta_treat_pre) - (eta_cont_post - eta_cont_pre)
real matrix lin_rep_ps, inf_treat_post, inf_treat_pret, inf_cont_post, inf_cont_pret
//score_ps = wgt :* (trt :- psc) :* xvar
//Hessian_ps = psv :* nn
lin_rep_ps = (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
inf_treat_post =(att_treat_post :- eta_treat_post) :-
((wgt :* trt :- pi_hat) :* eta_treat_post:/pi_hat) :-
((wgt :* tmt :- lambda_hat) :* eta_treat_post:/lambda_hat)
///inf_treat_post = inf_treat_post1 :+ inf_treat_post2 :+ inf_treat_post3
inf_treat_pret = att_treat_pre :- eta_treat_pre :-
(wgt :* trt :- pi_hat) :* eta_treat_pre:/pi_hat :-
(wgt :* (1 :- tmt) :- one_lambda_hat) :* eta_treat_pre:/one_lambda_hat
// inf_treat_pret = inf_treat_pre1 :+ inf_treat_pre2 :+ inf_treat_pre3
inf_cont_post = att_cont_post :- eta_cont_post :-
(wgt :* trt :- pi_hat) :* eta_cont_post:/pi_hat :-
(wgt :* tmt :- lambda_hat) :* eta_cont_post:/lambda_hat
///inf_cont_post = inf_cont_post1 :+ inf_cont_post2 :+ inf_cont_post3
inf_cont_pret = att_cont_pre :- eta_cont_pre :-
(wgt :* trt :- pi_hat) :* eta_cont_pre:/pi_hat :-
(wgt :* (1 :- tmt) :- one_lambda_hat) :* eta_cont_pre:/one_lambda_hat
//inf_cont_pret= inf_cont_pre1 :+ inf_cont_pre2 :+ inf_cont_pre3
real matrix inf_logit, att_inf_func
//mom_logit_pre = mean(-att_cont_pre :* xvar)
//mom_logit_pre = mean(mom_logit_pre)
//mom_logit_post = mean(-att_cont_post :* xvar)
//mom_logit_post = mean(mom_logit_post)
inf_logit = lin_rep_ps * (mean(-att_cont_post :* xvar) :-
mean(-att_cont_pre :* xvar))'
att_inf_func = ipw_att :+ (inf_treat_post :- inf_treat_pret) :- (inf_cont_post :- inf_cont_pret) :+ inf_logit
st_store(.,rif,touse,att_inf_func)
// -19.89330192
// 53.86822411
}
/// DRDID_RC
void drdid_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, y00, y01, y10, y11,
xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
yy = st_data(.,yy_, touse)
y00 = yy[,1]
y01 = yy[,2]
y10 = yy[,3]
y11 = yy[,4]
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else {
xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
}
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real matrix y0
real scalar nn
y0 = y00:*(-tmt:+1) + y01:*tmt
nn = rows(y)
real matrix w10, w11, w00, w01, w1
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
w1 = wgt :* trt
w00 = w00:/mean(w00 )
w01 = w01:/mean(w01 )
w10 = w10:/mean(w10 )
w11 = w11:/mean(w11 )
w1 = w1:/mean(w1 )
real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post,
att_trt_post, att_trtt1_post, att_trt_pre, att_trtt0_pre,
eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post,
eta_trt_post, eta_trtt1_post, eta_trt_pre, eta_trtt0_pre
att_treat_pre = w10 :* (y :- y0)
att_treat_post = w11 :* (y :- y0)
att_cont_pre = w00 :* (y :- y0)
att_cont_post = w01 :* (y :- y0)
att_trt_post = w1 :* (y11 :- y01)
att_trtt1_post = w11 :* (y11 :- y01)
att_trt_pre = w1 :* (y10 :- y00)
att_trtt0_pre = w10 :* (y10 :- y00)
eta_treat_pre = mean(att_treat_pre)
eta_treat_post = mean(att_treat_post)
eta_cont_pre = mean(att_cont_pre)
eta_cont_post = mean(att_cont_post)
eta_trt_post = mean(att_trt_post)
eta_trtt1_post = mean(att_trtt1_post)
eta_trt_pre = mean(att_trt_pre)
eta_trtt0_pre = mean(att_trtt0_pre)
real matrix trtr_att
trtr_att = (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre) :+ (eta_trt_post :- eta_trtt1_post) :- (eta_trt_pre :- eta_trtt0_pre)
real matrix wgt_ols_pre, XpX_inv_pre, lin_ols_pre,
wgt_ols_post, XpX_inv_post, lin_ols_post,
XpX_inv_pre_treat, lin_ols_pre_treat,
XpX_inv_post_treat, lin_ols_post_treat
wgt_ols_pre = wgt :* (1 :- trt) :* (1 :- tmt)
//wols_x_pre = wgt_ols_pre :* xvar
//wols_eX_pre = wgt_ols_pre :* (y :- y00) :* xvar
XpX_inv_pre = invsym(quadcross(xvar,wgt_ols_pre, xvar)):*nn
lin_ols_pre = ( wgt_ols_pre :* (y :- y00) :* xvar) * XpX_inv_pre
wgt_ols_post = wgt :* (1 :- trt) :* tmt
//wols_x_post = wgt_ols_post :* xvar
//wols_eX_post = wgt_ols_post :* (y :- y01) :* xvar
XpX_inv_post = invsym(quadcross(xvar,wgt_ols_post, xvar)):*nn
lin_ols_post = (wgt_ols_post :* (y :- y01) :* xvar) * XpX_inv_post
//wols_x_pre_treat = w10 :* xvar
//wols_eX_pre_treat = w10 :* (y :- y10) :* xvar
XpX_inv_pre_treat = invsym(quadcross(xvar, w10, xvar)):*nn
lin_ols_pre_treat = ( w10 :* (y :- y10) :* xvar) * XpX_inv_pre_treat
//wols_x_post_treat = w11 :* xvar
//wols_eX_post_treat = w11 :* (y :- y11) :* xvar
XpX_inv_post_treat = invsym(quadcross(xvar, w11, xvar)):*nn
lin_ols_post_treat = (w11 :* (y :- y11) :* xvar) * XpX_inv_post_treat
real matrix lin_rep_ps, inf_treat_pre, inf_treat_post
// check psv for probit
//score_ps = wgt :* (trt :- psc) :* xvar
//Hessian_ps = psv :* nn
lin_rep_ps = (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
inf_treat_pre = att_treat_pre :- w10 :* eta_treat_pre
inf_treat_post = att_treat_post :- w11 :* eta_treat_post
real matrix M1_post, M1_pre, inf_treat_or_post, inf_treat_or_pre
M1_post = -mean(w11 :* tmt :* xvar)
M1_pre = -mean(w10 :* (1 :- tmt) :* xvar)
inf_treat_or_post = lin_ols_post * M1_post'
inf_treat_or_pre = lin_ols_pre * M1_pre'
real matrix inf_treat_or, inf_treat, inf_cont_pre, inf_cont_post
inf_treat_or = inf_treat_or_post :+ inf_treat_or_pre
inf_treat = inf_treat_post :- inf_treat_pre :+ inf_treat_or
inf_cont_pre = att_cont_pre :- w00 :* eta_cont_pre
inf_cont_post = att_cont_post :- w01 :* eta_cont_post
real matrix M2_pre, M2_post, inf_cont_ps, M3_post, M3_pre, inf_cont_or_post, inf_cont_or_pre
M2_pre = mean(w00 :* (y :- y0 :- eta_cont_pre) :* xvar)
M2_post = mean(w01 :* (y :- y0 :- eta_cont_post) :* xvar)
inf_cont_ps = lin_rep_ps * (M2_post :- M2_pre)'
M3_post = -mean(w01 :* tmt :* xvar)
M3_pre = -mean(w00 :* (1 :- tmt) :* xvar)
inf_cont_or_post = lin_ols_post * M3_post'
inf_cont_or_pre = lin_ols_pre * M3_pre'
real matrix inf_cont_or, inf_cont, trtr_eta_inf_func1
inf_cont_or = inf_cont_or_post :+ inf_cont_or_pre
inf_cont = inf_cont_post :- inf_cont_pre :+ inf_cont_ps :+ inf_cont_or
trtr_eta_inf_func1 = inf_treat :- inf_cont
real matrix inf_eff, mom_post, mom_pre, inf_or
//inf_eff1 = att_trt_post :- w1 :* eta_trt_post
//inf_eff2 = att_trtt1_post :- w11 :* eta_trtt1_post
//inf_eff3 = att_trt_pre :- w1 :* eta_trt_pre
//inf_eff4 = att_trtt0_pre :- w10 :* eta_trtt0_pre
inf_eff = ((att_trt_post :- w1 :* eta_trt_post) :-
(att_trtt1_post :- w11 :* eta_trtt1_post)) :-
((att_trt_pre :- w1 :* eta_trt_pre) :-
(att_trtt0_pre :- w10 :* eta_trtt0_pre))
mom_post = mean((w1 :- w11) :* xvar)
mom_pre = mean((w1 :- w10) :* xvar)
// check this
//inf_or_post = ((lin_ols_post_treat :- lin_ols_post) * mom_post')
//inf_or_pre =
inf_or = ((lin_ols_post_treat :- lin_ols_post) * mom_post') :- ((lin_ols_pre_treat :- lin_ols_pre) * mom_pre')
att_inf_func = trtr_att :+ trtr_eta_inf_func1 :+ inf_eff :+ inf_or
st_store(.,rif,touse,att_inf_func)
// Variance should be divided by n not n-1
// -.1677954483
// .2008991705
}
/// DRDID_RC1
void drdid_rc1(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, y00, y01, y10, y11,
xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
yy = st_data(.,yy_, touse)
y00 = yy[,1]
y01 = yy[,2]
y10 = yy[,3]
y11 = yy[,4]
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real matrix y0
real scalar nn
y0 = y00:*(-tmt:+1) + y01:*tmt
nn = rows(y)
real matrix w10, w11, w00, w01
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
w00 = w00:/mean(w00 )
w01 = w01:/mean(w01 )
w10 = w10:/mean(w10 )
w11 = w11:/mean(w11 )
real matrix eta_trt_pre, eta_trt_post, eta_cont_pre, eta_cont_post,
att_trt_pre, att_trt_post, att_cont_pre, att_cont_post
eta_trt_pre = w10 :* (y :- y0)
eta_trt_post = w11 :* (y :- y0)
eta_cont_pre = w00 :* (y :- y0)
eta_cont_post = w01 :* (y :- y0)
att_trt_pre = mean(eta_trt_pre)
att_trt_post = mean(eta_trt_post)
att_cont_pre = mean(eta_cont_pre)
att_cont_post = mean(eta_cont_post)
real matrix trtr_att, w_ols_pre, wols_eX_pre, lin_rep_ols_pre,
w_ols_post, wols_eX_post, lin_rep_ols_post
trtr_att = (att_trt_post :- att_trt_pre) :- (att_cont_post :- att_cont_pre)
w_ols_pre = wgt :* (1 :- trt) :* (1 :- tmt)
wols_eX_pre = w_ols_pre :* (y :- y00) :* xvar
lin_rep_ols_pre = wols_eX_pre * invsym(quadcross(xvar, w_ols_pre, xvar)):*nn
w_ols_post = wgt :* (1 :- trt) :* tmt
wols_eX_post = w_ols_post :* (y :- y01) :* xvar
lin_rep_ols_post = wols_eX_post * invsym(quadcross( xvar,w_ols_post, xvar)):*nn
real matrix lin_rep_ps, inf_trt_pre, inf_trt_post, M1_post, M1_pre
//score_ps = wgt :* (trt :- psc) :* xvar
//Hessian_ps = psv :* nn
lin_rep_ps = (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
inf_trt_pre = eta_trt_pre :- w10 :* att_trt_pre
inf_trt_post = eta_trt_post :- w11 :* att_trt_post
M1_post = -mean(w11 :* tmt :* xvar)
M1_pre = -mean(w10 :* (1 :- tmt) :* xvar)
real matrix inf_trt_or, inf_trt
inf_trt_or = (lin_rep_ols_post * M1_post') :+ (lin_rep_ols_pre * M1_pre')
inf_trt = inf_trt_post :- inf_trt_pre :+ inf_trt_or
real matrix inf_cont_pre , inf_cont_post, M2_pre, M2_post, inf_cont_ps, M3_post, M3_pre
inf_cont_pre = eta_cont_pre :- w00 :* att_cont_pre
inf_cont_post = eta_cont_post :- w01 :* att_cont_post
M2_pre = mean(w00 :* (y :- y0 :- att_cont_pre) :* xvar)
M2_post = mean(w01 :* (y :- y0 :- att_cont_post) :* xvar)
inf_cont_ps = lin_rep_ps * (M2_post :- M2_pre)'
M3_post = -mean(w01 :* tmt :* xvar)
M3_pre = -mean(w00 :* (1 :- tmt) :* xvar)
real matrix inf_cont_or, inf_cont, att_inf_func
inf_cont_or = (lin_rep_ols_post * M3_post') :+ (lin_rep_ols_pre * M3_pre')
inf_cont = inf_cont_post :- inf_cont_pre :+ inf_cont_ps :+ inf_cont_or
att_inf_func = trtr_att :+ inf_trt :- inf_cont
st_store(.,rif,touse,att_inf_func)
// Variance should be divided by n not n-1
// -3.633433441
//3.107123089
}
/// drdid_imp_rc
void drdid_imp_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, y00, y01, y10, y11,
xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
yy = st_data(.,yy_, touse)
y00 = yy[,1]
y01 = yy[,2]
y10 = yy[,3]
y11 = yy[,4]
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
// xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real matrix y0
real scalar nn
y0 = y00:*(-tmt:+1) + y01:*tmt
nn = rows(y)
real matrix w10, w11, w00, w01, w1
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
w1 = wgt :* trt
w00 = w00:/mean(w00 )
w01 = w01:/mean(w01 )
w10 = w10:/mean(w10 )
w11 = w11:/mean(w11 )
w1 = w1:/mean(w1 )
real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post, att_trt_post, att_trtt1_post,
att_trt_pre, att_trtt0_pre
att_treat_pre = w10 :* (y :- y0)
att_treat_post = w11 :* (y :- y0)
att_cont_pre = w00 :* (y :- y0)
att_cont_post = w01 :* (y :- y0)
att_trt_post = w1 :* (y11 :- y01)
att_trtt1_post = w11 :* (y11 :- y01)
att_trt_pre = w1 :* (y10 :- y00)
att_trtt0_pre = w10 :* (y10 :- y00)
real matrix eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post, eta_trt_post, eta_trtt1_post,
eta_trt_pre, eta_trtt0_pre
eta_treat_pre = mean(att_treat_pre)
eta_treat_post = mean(att_treat_post)
eta_cont_pre = mean(att_cont_pre)
eta_cont_post = mean(att_cont_post)
eta_trt_post = mean(att_trt_post)
eta_trtt1_post = mean(att_trtt1_post)
eta_trt_pre = mean(att_trt_pre)
eta_trtt0_pre = mean(att_trtt0_pre)
real matrix trtr_att
trtr_att = (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre) :+ (eta_trt_post :- eta_trtt1_post) :- (eta_trt_pre :- eta_trtt0_pre)
real matrix inf_treat, inf_cont, att_inf_func1, inf_eff, att_inf_func
inf_treat = (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre :- w10 :* eta_treat_pre)
inf_cont = (att_cont_post :- w01 :* eta_cont_post) :- (att_cont_pre :- w00 :* eta_cont_pre)
att_inf_func1 = inf_treat :- inf_cont
inf_eff = ((att_trt_post :- w1 :* eta_trt_post) :-
(att_trtt1_post :- w11 :* eta_trtt1_post)) :-
((att_trt_pre :- w1 :* eta_trt_pre) :-
(att_trtt0_pre :- w10 :* eta_trtt0_pre))
att_inf_func = trtr_att :+ att_inf_func1 :+ inf_eff
st_store(.,rif,touse,att_inf_func)
// Variance should be divided by n not n-1
//-.2088586355
//.2003375215
}
/// drdid_imp_rc1
void drdid_imp_rc1(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ , touse, rif){
// main Loading variables
real matrix y, y00, y01, y10, y11,
xvar, tmt, trt, psv, psc, wgt
y = st_data(.,y_, touse)
yy = st_data(.,yy_, touse)
y00 = yy[,1]
y01 = yy[,2]
y10 = yy[,3]
y11 = yy[,4]
// verify xvar
if (xvar_==" ") {
xvar=J(rows(y),1,1)
}
else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
tmt = st_data(.,tmt_, touse)
trt = st_data(.,trt_, touse)
psc = logistic(st_data(.,pxb_, touse))
wgt = st_data(.,wgt_, touse)
psv = st_matrix(psv_)
real matrix y0
real scalar nn
y0 = y00:*(-tmt:+1) + y01:*tmt
nn = rows(y)
real matrix w10, w11, w00, w01,
eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post,
att_treat_pre, att_treat_post, att_cont_pre, att_cont_post
w10 = wgt :* trt :* (1 :- tmt)
w11 = wgt :* trt :* tmt
w00 = wgt :* psc:* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
w01 = wgt :* psc:* (1 :- trt) :* tmt:/(1 :- psc)
eta_treat_pre = w10 :* (y :- y0):/mean(w10)
eta_treat_post = w11 :* (y :- y0):/mean(w11)
eta_cont_pre = w00 :* (y :- y0):/mean(w00)
eta_cont_post = w01 :* (y :- y0):/mean(w01)
att_treat_pre = mean(eta_treat_pre)
att_treat_post = mean(eta_treat_post)
att_cont_pre = mean(eta_cont_pre)
att_cont_post = mean(eta_cont_post)
real matrix trtr_att, inf_treat, inf_cont, att_inf_func
trtr_att = (att_treat_post :- att_treat_pre) :- (att_cont_post :- att_cont_pre)
inf_treat = (eta_treat_post :- w11 :* att_treat_post:/mean(w11) ):- (eta_treat_pre :- w10 :* att_treat_pre:/mean(w10))
inf_cont = (eta_cont_post :- w01 :* att_cont_post:/mean(w01)) :- ( eta_cont_pre :- w00 :* att_cont_pre:/mean(w00))
att_inf_func = trtr_att :+ inf_treat :- inf_cont
// Wrapping up
st_store(.,rif,touse,att_inf_func)
// Variance should be divided by n not n-1
// -3.683728719
// 3.114495585
}
// Clustered Standard errors
void clusterse(real matrix iiff, cl, V, real scalar cln){
/// estimates Clustered Standard errors
real matrix ord, xcros, ifp, info, vv
ord = order(cl,1)
iiff = iiff[ord,]
cl = cl[ord,]
// check how I cleaned data!
info = panelsetup(cl,1)
// faster Cluster? Need to do this for mmqreg
ifp = panelsum(iiff,info)
xcros = quadcross(ifp,ifp)
real scalar nt, nc
nt=rows(iiff)
nc=rows(info)
V = xcros/(nt^2)
cln=nc
// Esto es para ver como hacer clusters.
//*nc/(nc-1)
//st_matrix(V, vv)
//st_numscalar(ncl, nc)
// ^ ^
// | |
// stata mata
}
//** Simple Bootstrap.
real matrix mboot_did(real matrix rif, mean_rif, real scalar reps, wbtype) {
real matrix yy, bsmean
yy=rif:-mean_rif
bsmean=J(reps,cols(yy),0)
real scalar i,n, k1, k2
n=rows(yy)
k1=((1+sqrt(5))/(2*sqrt(5)))
k2=0.5*(1+sqrt(5))
// WBootstrap:Mammen
if (wbtype==1) {
for(i=1;i<=reps;i++){
bsmean[i,]=mean(yy:*(k2:-sqrt(5)*(rbinomial(n,1,1,k1))) )
}
}
else if (wbtype==2) {
for(i=1;i<=reps;i++){
bsmean[i,]=mean(yy:*(1:-2*rbinomial(n,1,1,0.5) ) )
}
}
return(bsmean)
}
real matrix mboot_didc(real matrix rif, mean_rif, clv , real scalar reps, wbtype, nc) {
real matrix yy, bsmean
real matrix sclv, wmult
real scalar i,n, k1, k2, nn
yy=rif:-mean_rif
bsmean=J(reps,cols(yy),0)
nn=rows(uniqrows(clv))
nc=nn
k1=((1+sqrt(5))/(2*sqrt(5)))
k2=0.5*(1+sqrt(5))
if (wbtype==1) {
for(i=1;i<=reps;i++){
wmult=rbinomial(nn,1,1,k1)
bsmean[i,]=mean(yy:*(k2:-sqrt(5)*wmult[clv] ) )
}
}
else if (wbtype==2) {
for(i=1;i<=reps;i++){
wmult=(rbinomial(nn,1,1,0.5))
bsmean[i,]=mean(yy:*(1:-2* wmult[clv] ) )
}
}
return(bsmean)
}
///clust("`att'","`cluster'","`touse'", "`V'",
///mboot("`att'", "`touse'", "`V'","`ciband'","`cluster'", `reps', `wbtype', `ci')
void make_tbl(string scalar iiff, clv, touse, // RIF, Cluster variable, and sample
bb_ , VV_, ccbb, cln, // WHere to save bb and VV. Also #clv
wboot, // options for Wbootstrap. Reps, wbtype and CI
real scalar reps, wbtype, ci ){
real matrix rif, clvar , bb, VV, ord, info, ifp
real matrix cband
real scalar nobs, nc
// nc number of clusters
rif = st_data(.,iiff,touse)
// real scalar cln
bb = mean(rif)
nobs = rows(rif)
// simple
if ((clv==" ") & (wboot==" ")) {
VV=quadcrossdev(rif,bb,rif,bb):/ (nobs^2)
}
// cluster std
if ((clv!=" ") & (wboot==" ")) {
clvar=st_data(.,clv,touse)
ord = order(clvar,1)
rif = rif[ord,]
clvar = clvar[ord,]
info = panelsetup(clvar,1)
ifp = panelsum((rif:-bb),info)
VV = quadcross(ifp,ifp)/nobs:^2
nc=rows(info)
}
// wboot no cluster
if (wboot!=" ") {
mboot(rif,bb, VV, cband, clv, touse ,reps, wbtype, ci, nc )
}
st_matrix(bb_,bb)
st_matrix(VV_,VV)
st_matrix(ccbb,cband)
st_numscalar(cln, nc)
}
real matrix iqrse(real matrix y) {
real scalar q25,q75
q25=ceil((rows(y)+1)*.25)
q75=ceil((rows(y)+1)*.75)
real scalar j
real matrix iqrs
iqrs=J(1,cols(y),0)
for(j=1;j<=cols(y);j++){
y=sort(y,j)
iqrs[,j]=(y[q75,j]-y[q25,j]):/(invnormal(.75)-invnormal(.25) )
}
return(iqrs)
}
real vector qtp(real matrix y, real scalar p) {
real scalar k, i, q
real matrix yy, qq
qq=J(1,0,.)
k = cols(y)
for(i=1;i<=k;i++){
yy=sort(y[,i],1)
q=ceil((rows(yy)+1)*p)
qq=qq,yy[q,]
}
return(qq)
}
void mboot(real matrix rif,mean_rif, vv, cband, string scalar clv, touse, real scalar reps, wbtype, ci , nc ) {
//, real scalar reps, wbtype, ci
real matrix fr, qt
//real scalar reps, wbtype
real matrix ifse
// this gets the Bootstraped values
if (clv ==" ") {
fr=mboot_did(rif,mean_rif , reps, wbtype)
ifse = iqrse(fr)
qt = qtp(abs(fr :/ ifse),ci)
cband=( mean_rif':-qt':* ifse' ,
mean_rif':+qt':* ifse' , mean_rif, ifse, qt' )
}
else {
clvar=st_data(.,clv, touse)
fr=mboot_didc(rif,mean_rif, clvar, reps, wbtype, nc )
ifse = iqrse(fr)
qt = qtp(abs(fr :/ ifse),ci)
cband=( mean_rif':-qt':* ifse' ,
mean_rif':+qt':* ifse' , mean_rif, ifse, qt' )
}
vv=quadcross(ifse,ifse):*I(cols(ifse))
}
/// qtp(abs(xx/ iqrse(xx)),.95)
end
program addr, rclass
return `0'
end