* Next step Integrate all into csdid.mata <- so Do not need to play with many files *! v1.72 by FRA. Drops always treated * v1.71 by FRA. adds weights back * v1.7 by FRA. Changes on Datacheks. This should avoid time gaps problems. * Adds asinr for pretreatmetn * adds version * Need to add Checks for Covariates and variables. * v1.6 by FRA. Changes how agg looks * v1.58 by FRA. Dropping binit. Problems with drimp * v1.58 by FRA. adding base info * v1.57 by FRA. change default to dripw from drimp. More stable * Also corrects agg(event) * Does not create Unecessary RIFS * v1.57 by FRA. Takes Pretrend from AGG NOT ALLOWED * v1.56 by FRA. bug for NOT YET, 2021 ref * v1.55 by FRA. Allows for more periods * v1.53 by FRA. changes e(gtt). THis give sdetailed sample * v1.52 by FRA. Adds option "from" for simple aggregation * v1.51 by FRA. Added seed and method. for bootstrap * v1.3 by FRA. Overhaul Change in how Weights are estimated. All effects are now in Mata. * All effects in mata. Quick! * v1.03 by FRA. Adding Balance checks * v1.02 by FRA. added touse * v1.01 by FRA. Make sure we have at least 1 not treated period. // and check that if no Never treated exist, we add extra variable for dummies. // exclude if There are no treated observations * v1 by FRA csdid Added all tests! and ready to play * v0.2 by FRA csdid or Multiple periods did. Adds Notyet * v0.1 by FRA mpdid or Multiple periods did. * This version implements did::attgp * Logic. Estimate the ATT of the base (first) year against all subsequent years * using data BY groups ** assumes all years are available. For now /*program csdid_check, syntax varlist(fv ) [if] [in] [iw], /// Basic syntax allows for weights [Ivar(varname)] /// Time(varname) /// [gvar(varname)] /// Ppl either declare gvar [trvar(varname)] /// or declare treat var, and we create gvar. Only works if Panel data [att_gt] /// [notyet] /// att_gt basic option. May prepare others as Post estimation [method(str) ] // This allows other estimators end*/ *capture program drop sdots *capture program drop csdid *program drop mynote *This should break Notes into smaller pieces. Still unsure about how to pull them all together program mynote, syntax anything gettoken lc rest:anything local k =0 if length("`rest'")>60000 { local cnt = 0 foreach i in `rest' { local ll `ll' `i' local cnt = `cnt'+length("`i'") if `cnt'>60000 { local cnt=0 local k = `k'+1 local ll`k' `ll' local ll } } } *return local knote `k' *return local note_nm `lc' if `k'==0 { note : `lc' `rest' } else { note : `lc' "see_char `k'" forvalues kk = 1/`k' { char _dta[`lc'`kk'] "`ll`kk''" } } end ** Parsing notes: program sdots syntax , [mydots(integer 10) maxdots(int 50) bad] if mod(`mydots',`maxdots')==0 { if "`bad'"=="" display "." else display "x" } else { if "`bad'"=="" display "." _c else display "x" _c } end mata: // This should check all combinations for Tvar and Gvar // Brute force approach to fix missing data void data_check(string scalar tvar, string scalar gvar, string scalar touse){ real matrix tgvar , ttvar, ggvar, dtvar st_view(ttvar=.,.,tvar,touse) st_view(ggvar=.,.,gvar,touse) // Get full set of tvar gvar ttvar = uniqrows(ttvar)' ggvar = uniqrows(ggvar)' st_local("no_notyet",strofreal(sum(ggvar:==0))) ggvar = select(ggvar,ggvar:>0) // get MIN time delta dtvar = min(ttvar[2..cols(ttvar)]:-ttvar[1..cols(ttvar)-1]) // Pre treatment tgvar = ggvar:-dtvar // get all possible Times ttvar = range(min(ttvar),max(ttvar),dtvar)' // st_local("mintime",strofreal(min(ttvar))) st_local("time0" ,strofreal(min(ttvar))) st_local("gdelta" ,strofreal(min(dtvar))) st_local("mingvar",strofreal(min(ggvar))) st_local("maxgvar",strofreal(max(ggvar))) st_local("glev" ,invtokens(strofreal(ggvar))) st_local("tlev" ,invtokens(strofreal(ttvar[2..cols(ttvar)]))) st_local("tlev1",invtokens(strofreal(ttvar))) st_local("precheck",strofreal(max(ggvar))) } end program csdid, sortpreserve eclass version 14 ///version checker syntax [anything(everything)] [iw aw pw], [* version] if "`version'"!="" { display "version: 1.72" addr scalar version = 1.72 exit } capture drdid, version if _rc!=0 { display in red "Program DRDID is outdated, or not installed. " display in red "Please install ssc install drdid" *exit 101 } if _rc==0 { if r(version)<1.71 display in red "Program DRDID is outdated. Please update" as text *exit 101 } //// if replay() { if `"`e(cmd)'"' != "csdid" { error 301 } else { Display `0' } exit } if runiform()<.001 { easter_egg } csdid_r `0' *ereturn local cmdline `"drdid `0'"' end program define Display syntax [, bmatrix(passthru) vmatrix(passthru) *] _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(Difference-in-difference with Multiple Time Periods) noi display as text "Outcome model : {res:`omodel'}" noi display as text "Treatment model: {res:`tmodel'}" } if ("`e(vcetype)'"=="WBoot") { 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 ")" } csdid_table, `diopts' } else { _coef_table, `diopts' `myopts' neq(`e(neqr)') } end program define _S_Me_thod, sclass if ("`e(method)'"=="drimp") { local tmodel "inverse probability tilting" local omodel "weighted least squares" } if ("`e(method)'"=="dripw") { local tmodel "inverse probability" local omodel "least squares" } if ("`e(method)'"=="aipw") { local tmodel "inverse probability" local omodel "weighted mean" } if ("`e(method)'"=="reg") { local tmodel "none" local omodel "regression adjustment" } if ("`e(method)'"=="stdipw") { local tmodel "stabilized inverse probability" local omodel "weighted mean" } sreturn local omodel "`omodel'" sreturn local tmodel "`tmodel'" end program _Parse2_wboot, rclass syntax, [reps(integer 999) wbtype(str) rseed(str) cluster(str)] return scalar reps = `reps' return local seed `rseed' return local cluster `cluster' if ("`wbtype'"=="") { return scalar wbtype = 1 } else if ("`wbtype'"=="mammen") { return scalar wbtype = 1 } else if ("`wbtype'"=="rademacher") { return scalar wbtype = 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 } end program _Parse_Wildboot, rclass syntax , [ /// WBOOT1 /// WBOOT(str) /// reps(integer 999) /// rseed(string) /// wbtype(string) /// cluster(string) /// ] marksample touse if ("`wboot1'"=="" & "`wboot'"!="") { _Parse2_wboot, `wboot' return scalar reps= `r(reps)' return local seed `r(seed)' return scalar wbtype= `r(wbtype)' return local cluster `r(cluster)' } else if ("`wboot1'"!="") { return scalar 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 scalar 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 } } } return local cluster `cluster' } end program issaverif, syntax, [saverif(str) replace] if "`saverif'"!="" { capture confirm new file "`saverif'.dta" if _rc!=0 { if "`replace'"=="" { display in red "File `saverif'.dta already exists" display in red "Please use a different file name or use -replace-" exit 602 } else { display "File `saverif'.dta will be replaced" } } else { display "File `saverif'.dta will be used to save all RIFs" } } end program csdid_r, sortpreserve eclass syntax varlist(fv ) [if] [in] [iw], /// Basic syntax allows for weights [Ivar(varname numeric)] /// Time(varname numeric) /// Gvar(varname numeric) /// Ppl either declare gvar [cluster(varname numeric)] /// [notyet] /// att_gt basic option. May prepare others as Post estimation [saverif(name) replace ] /// [method(str) /// agg(str) /// WBOOT(str) /// WBOOT1 /// *reps(int 999) /// *wbtype(str) /// Hidden option rseed(str) /// set seed Level(int 95) /// CI level pointwise /// with Wboot from(int 0) /// for aggregations long long2 /// to allow for "long gaps" dryrun /// for testing asinr /// For pretreatment ] // This allows other estimators ** Weight correction *if "`weight'"!="" local weight iw marksample touse ** First determine outcome and xvars gettoken y xvar:varlist markout `touse' `ivar' `time' `gvar' `y' `xvar' `cluster' tempname cband ** Parsing WBOOT _Parse_Wildboot, wboot(`wboot') `wboot1' rseed(`rseed') cluster(`cluster') if "`rseed'"!="" set seed `rseed' if "`wboot'`wboot1'"!="" { if "`rseed'"!="" set seed `rseed' local cluster `r(cluster)' local ocluster `r(cluster)' local seed `r(seed)' if "`wbtype'"=="" local owbtype mammen local owbtype `wbtype' local wbtype `r(wbtype)' local reps `r(reps)' local vcetype WBoot if "`pointwise'"=="" local citype "uniform" else local citype "pointwise" } else { local wbtype 1 } *********************************************** if "`method'"=="" { local method dripw if "`xvar'"=="" { local method reg } } else if !inlist("`method'","drimp","dripw","reg","ipw","stdipw") { display as error "Method `method' not allowed" } ** Default. If no controls, use reg if "`seed'"!="" set seed `seed' ** Original cluster local ocluster `cluster' ** Confirm if new file exists issaverif, saverif(`saverif') `replace' if "`agg'"=="" { local agg attgt } if !inlist("`agg'","attgt","simple","group","calendar","event") { display in red "Aggregation not Allowed" exit 10 } /// added data check qui:mata:data_check("`time'","`gvar'","`touse'") // verify not yet if "`tyet'"=="" { // if never treated if `no_notyet'==0 { local tyet notyet display "No never treated observations found. Using Not yet treated data" } else { sum `time' if `touse' & `gvar'==0, meanonly qui:replace `touse'=0 if `time'>`r(max)' } } /// Added Restrictions to sample *!! Stop here. How to restrict from going over a group that isnt going to estimate anything? if "`tyet'"!="" { // if never treated sum `time' if `touse' & `gvar'==0, meanonly local tt_max = r(max) sum `gvar' if `touse' , meanonly local gg_max = r(max) qui:replace `touse'=0 if `time'>`gg_max' & `time'>`tt_max' } qui:mata:data_check("`time'","`gvar'","`touse'") ** checking for Min gvar being larger than min time /* Now done with data_check qui:sum `time' if `touse', meanonly local mintime=r(min) qui:sum `gvar' if `touse' & `gvar'>0, meanonly local mingvar=r(min) local maxgvar=r(max) */ if `maxgvar'<=`mintime' { display "Gvar max value is `maxgvar' and there are no periods available before that treatment" display "Verify that Gvar is correctly defined" error 190 } if `mintime'>=`mingvar' { display "Units always treated found. These will be ignored" } qui:replace `touse'=0 if (`gvar'<=`mintime') & (`gvar'>0) & `touse' qui:mata:data_check("`time'","`gvar'","`touse'") ** determine time0 /*if "`time0'"=="" { qui:sum `time' if `touse' local time0 `r(min)' } */ ** Check if panel is Full Balance qui { // Only if panel if "`ivar'"!="" { ** First check if data is Panel tempname ispan capture xtset if _rc!=0 { qui:xtset `ivar' `time' qui:xtset , clear } if "`cluster'"!="" { _xtreg_chk_cl2 `cluster' `ivar' } tempname isbal mata:isbalanced("`ivar'","`touse'","`isbal'") if scalar(`isbal')!=1 { display in red "Panel is not balanced" if "`bal'"=="full" { display in red "Will use Only observations fully balanced (Max periods observed)" replace `touse'=0 if `x1'<`mxcol' } if "`bal'"=="" { display in red "Will use observations with Pair balanced (observed at t0 and t1)" } if "`bal'"=="unbal" { display in red "Will use all observations (as if CS)" local cluster `ivar' local ivar } } } } ** prepare loops over gvar. *local att_gt att_gt tempvar tr tempname gtt * qui:gen byte `tr'=`gvar'!=0 if `gvar'!=. *qui:levelsof `gvar' if `gvar'>0 & `touse', local(glev) *qui:levelsof `time' if `time'>`time0' & `touse', local(tlev) ** Checks If Gvar is ever untreated if `precheck'==1 { display in red "All observations require at least 1 not treated period" display in red "See cross table below, and verify All Gvar have at least 1 not treated period" tab `time' `gvar' if `touse' error 1 } tempname b v tempvar gsel tempvar cs_sample capture qui:gen byte _blank_=. qui:gen byte `cs_sample'=0 //////////////////////////////////////////////////////////////////////////////// if "`long'`long2'"=="" { foreach i of local glev { local time1 = `time0' foreach j of local tlev { local mydots=`mydots'+1 * Stata quirk: `notyet' is called `tyet' because "no" is removed ** This implements the Never treated capture drop `gsel' qui:gen byte `gsel'=inlist(`gvar',0,`i') if `touse' if "`tyet'"!="" { if "`asinr'"=="" qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | (`gvar'> `i' & `gvar' >`j' )) & `touse' if "`asinr'"!="" { if `j'< `i' qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | `gvar' >`j' ) & `touse' if `j'>=`i' qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | (`gvar'> `i' & `gvar' >`j' )) & `touse' } } capture drop `tr' qui:gen byte `tr'=(`gvar'==`i') if `touse' qui:capture:drdid `varlist' if `gsel' /// & inlist(`time',`time1',`j') /// & `touse' [`weight'`exp'], /// ivar(`ivar') time(`time') treatment(`tr') /// `method' stub(__) replace `dryrun' *binit(`bii' ) *tempname bii if _rc == 1 { display in red "Stopped by user" capture drop `vlabrif' capture drop __att error 11882 } if _rc ==0 { qui:replace `cs_sample'=e(sample) if `cs_sample'==0 qui:count if e(sample)==1 & `tr'==1 matrix `gtt'=nullmat(`gtt')\[`i',min(`time1',`j'),max(`time1',`j'),_rc,e(N),`=e(N)-r(N)',r(N)] *matrix `bii'=e(iptb) *matrix list `bii' } if _rc!=0 { local bad bad matrix `gtt'=nullmat(`gtt')\[`i',min(`time1',`j'),max(`time1',`j'),_rc,.,.,.] } sdots, mydots(`mydots') `bad' local bad local eqname `eqname' g`i' local colname `colname' t_`time1'_`j' capture drop _g`i'_`time1'_`j' capture confirm variable __att if _rc==111 { local vlabrif `vlabrif' _blank_ local rifvar `rifvar' _blank_ if `j'<`i' local time1 = `j' } else { ren __att _g`i'_`time1'_`j' local vlabrif `vlabrif' _g`i'_`time1'_`j' if `j'<`i' local time1 = `j' local rifvar `rifvar' _g`i'_`time1'_`j' } } } } //////////////////////////////////////////////////////////////////////////////// if "`long'`long2'"!="" { ** times goes from t0 *qui:levelsof `time' if `time'>=`time0' & `touse', local(tlev) foreach i of local glev { *qui:sum `time' if `time'<`i' & `touse', meanonly local time1 = `i'-`gdelta' foreach j of local tlev1 { ****************************************** if `time1'!=`j' { local mydots=`mydots'+1 * Stata quirk: `notyet' is called `tyet' because "no" is removed ** This implements the Never treated capture drop `gsel' qui:gen byte `gsel'=inlist(`gvar',0,`i') if `touse' if "`tyet'"!="" { if "`asinr'"=="" qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | (`gvar'> `i' & `gvar' >`j' )) & `touse' if "`asinr'"!="" { if `j'< `i' qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | `gvar' >`j' ) & `touse' if `j'>=`i' qui:replace `gsel'=1 if (`gvar'==0 | `gvar'==`i' | (`gvar'> `i' & `gvar' >`j' )) & `touse' } } capture drop `tr' qui:gen byte `tr'=(`gvar'==`i') if `touse' qui:capture:drdid `varlist' if `gsel' /// & inlist(`time',`time1',`j') /// & `touse' [`weight'`exp'], /// ivar(`ivar') time(`time') treatment(`tr') /// `method' stub(__) replace `dryrun' noisily *binit(`bii',skip) *tempname bii if _rc == 1 { display in red "Stopped by user" capture drop `vlabrif' capture drop __att error 11882 } if _rc ==0 { qui:replace `cs_sample'=e(sample) if `cs_sample'==0 qui:count if e(sample)==1 & `tr'==1 matrix `gtt'=nullmat(`gtt')\[`i',min(`time1',`j'),max(`time1',`j'),_rc,e(N),`=e(N)-r(N)',r(N)] *matrix `bii'=e(iptb) } if _rc!=0 { local bad bad matrix `gtt'=nullmat(`gtt')\[`i',min(`time1',`j'),max(`time1',`j'),_rc,.,.,.] } sdots, mydots(`mydots') `bad' local bad local eqname `eqname' g`i' if `time1'<`j' { local colname `colname' t_`time1'_`j' capture drop _g`i'_`time1'_`j' capture confirm variable __att if _rc==111 { *qui: gen _blank_=. local vlabrif `vlabrif' _blank_ local rifvar `rifvar' _blank_ } else { ren __att _g`i'_`time1'_`j' local vlabrif `vlabrif' _g`i'_`time1'_`j' local rifvar `rifvar' _g`i'_`time1'_`j' } } else { local colname `colname' t_`j'_`time1' capture drop _g`i'_`j'_`time1' capture confirm variable __att if _rc==111 { *qui: gen _blank_=. local vlabrif `vlabrif' _blank_ local rifvar `rifvar' _blank_ } else { ren __att _g`i'_`j'_`time1' if "`long2'"!="" qui:replace _g`i'_`j'_`time1'=-_g`i'_`j'_`time1' local vlabrif `vlabrif' _g`i'_`j'_`time1' local rifvar `rifvar' _g`i'_`j'_`time1' } } } ****************************************** } } } ** names for Weights To be changed foreach i of local glev { local neqr = `neqr'+1 foreach j of local tlev { local eqname `eqname' wgt local colname `colname' w`i'_`j' } } //////////////////////////////////////////////////////////////////////////////// preserve qui:notes drop _all tempvar wgtt if "`exp'"!="" clonevar `wgtt'`exp' else gen byte `wgtt'=1 qui:keep if `touse' qui:keep `ivar' `time' `gvar' `vlabrif' `wgtt' `cluster' _blank_ local oivar `ivar' if "`ivar'"=="" { qui:gen double ivar=_n label var ivar "indicator" local ivar ivar } local svlist `time' `gvar' `vlabrif' `wgtt' `cluster' local svlist = subinstr("`svlist'","_blank_","",.) collapse `svlist' _blank_ , by(`ivar' ) fast ren `wgtt' __wgt__ label var __wgt__ "Weight Variable" *qui:levelsof `gvar' if `gvar'!=0, local(gglev) local kc = 0 foreach i of local glev { foreach j of local tlev { local kc=`kc'+1 local rvl:word `kc' of `vlabrif' if "`rvl'"!="_blank_" { qui:gen double w`i'_`j'=0 qui:replace w`i'_`j'=1 if (`rvl'!=.) & (`gvar'==`i') local lvl_gvar `lvl_gvar' w`i'_`j' } else { local lvl_gvar `lvl_gvar' _blank_ } } } *** Organizing all data. *** first Gvars if "`cluster'"!="" { qui:egen long _cl_var=group(`cluster') *qui:sort _cl_var label var _cl_var "Effective cluster" local cluster _cl_var } /// saving RIF note: Data created with -csdid-. Contains all -RIFs- associated with model estimation. note: cmdline csdid `0' note: time0 `time0' note: tlvls `tlev' note: glvls `glev' note: clvar `cluster' mynote rifgt `vlabrif' *note: rifgt `vlabrif' mynote rifwt `lvl_gvar' *note: rifwt `lvl_gvar' mynote colname `colname' *note: colname `colname' mynote eqname `eqname' *note: eqname `eqname' note: cmd csdid note: typebase `long2'`long' /// parsing notes. forvalues i = 3/8 { local `:char _dta[note`i']' } ** checking Notes 7 and 8 if "`:word 1 of `rifgt''"=="see_char" { local mk "`:word 2 of `rifgt''" local rifgt forvalues i = 1/`mk'{ local rifgt `rifgt' `:char _dta[rifgt`i']' } } if "`:word 1 of `rifwt''"=="see_char" { local mk "`:word 2 of `rifwt''" local rifwt forvalues i = 1/`mk'{ local rifwt `rifwt' `:char _dta[rifwt`i']' } } /// "Option will define the following // 1 simple // 2 dynamic // 3 calendar // 4 group // 5 attgt // 6 pretrend // 7 all tempname b1 b2 b3 b4 b5 tempname s1 s2 s3 s4 s5 ** New idea. Hacerlo todo desde makerif *mata:makerif("`rifgt'","`rifwt'","__wgt__","`b'","`v'","`cluster' ") *save extra, replace local ci `level'/100 local citp=1 if "`wboot'`wboot1'"!="" { local wboot wboot if "`citype'"=="uniform" local citp = 1 if "`citype'"=="pointwise" local citp = 2 } local l2 = 0 if "`long2'"!="" local l2 = 1 noisily mata: makerif2("`rifgt'" , "`rifwt'","__wgt__","`agg'", /// "`time0'","`glvls'","`tlvls'", /// "`b1'", /// `b2' `b3' `b4' `b5' `b6' "`s1'", /// `s2' `s3' `s4' `s5' `s6' "`clvar' ","`wboot' ", "`cband'", /// `ci', `reps', `wbtype', `citp',`from', `l2') ** Time0 `time0' if "`saverif'"!="" { save `saverif', replace } restore capture drop `vlabrif' capture drop _blank_ //////////////////////////////////////////////////////////////////////////////// if "`agg'"=="attgt" { matrix colname `b1'=`colname' matrix coleq `b1'=`eqname' matrix colname `s1'=`colname' matrix coleq `s1'=`eqname' matrix rowname `s1'=`colname' matrix roweq `s1'=`eqname' } matrix colname b_attgt=`colname' matrix coleq b_attgt=`eqname' matrix colname V_attgt=`colname' matrix coleq V_attgt=`eqname' matrix rowname V_attgt=`colname' matrix roweq V_attgt=`eqname' *** Mata to put all together. matrix colname `gtt' = cohort t0 t1 error N N_cntr N_trt quietly count if `cs_sample' local N = r(N) ereturn post `b1' `s1', esample(`cs_sample') obs(`N') ereturn local cmd csdid ereturn local cmdline csdid `0' ereturn local estat_cmd csdid_estat ereturn matrix b_attgt b_attgt ereturn matrix V_attgt V_attgt ereturn matrix gtt `gtt' ereturn local agg `agg' ereturn local glev `glev' ereturn local tlev `tlev' ereturn local time0 `time0' ereturn local rif `rifvar' ereturn local ggroup `gvar' ereturn local id `ivar' ereturn scalar neqr = `neqr' ereturn local riffile `saverif' ereturn local method `method' if "`long'`long2'"!="" ereturn local typebase `long'`long2' else ereturn local typebase default if "`ocluster'"!="" { ereturn local clustvar `ocluster' ereturn scalar N_clust = `=scalar(clnm)' } /// Add here Cluster var and number of Clusters. if "`tyet'"=="" ereturn local control_group "Never Treated" if "`tyet'"!="" ereturn local control_group "Not yet Treated" if "`vcetype'"=="WBoot" { ereturn local seed `seed' ereturn local wbtype `wbtype' ereturn scalar reps = `reps' matrix colname `cband'= b se t ll ul ereturn matrix cband `cband' ereturn local vcetype WBoot ereturn local citype `citype' } Display, level(`level') display "Control: `e(control_group)'" display _n "See Callaway and Sant'Anna (2021) for details" end /// This can be used for aggregation. Creates the matrixes we need. 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. " /// "It is the second part of Pedro, Brantly and Juns's contribution to the DID world{p_end} " _n /// "{p} Remember One Difference is good, and 2x2 DiD is twice as good!. " /// " Just dont confuse it with DnD (Dungeons and Dragons){p_end} " end mata: vector event_list(real matrix glvl, tlvl,window){ real matrix toreturn, toreturn2 real scalar i,j toreturn=J(1,0,.) toreturn2=J(1,0,.) for(i=1;i<=cols(glvl);i++) { for(j=1;j<=cols(tlvl);j++) { toreturn=toreturn,(tlvl[j] -glvl[i]) } } toreturn=uniqrows(toreturn')' if (cols(window)==0) return(toreturn) else { for(i=1;i<=cols(toreturn);i++){ if ( (toreturn[i]>=window[1]) & (toreturn[i]<=window[2]) ) toreturn2=toreturn2,toreturn[i] } return(toreturn2) } } // Next task. // amek all elements separete RIF_siple RIF event, etc // Think how to save all elements. void makerif2(string scalar rifgt_ , rifwt_ , wgt_, agg, time0, glvl_, tlvl_, bb_, ss_, clvar_, wboot_ , cband, real scalar ci, reps, wbtype, citype, from, real scalar lng ) { real matrix rifgt , rifwt, wgt, t0, glvl, tlvl real scalar i,j,k,h rifgt = st_data(.,rifgt_) rifwt = st_data(.,rifwt_) // Because of blanks we need to make it zero. Otherwise there is an error _editmissing(rifwt,0) wgt = st_data(.,wgt_) wgt = wgt:/mean(wgt) rifwt = rifwt:*wgt rifwt = rifwt:/rowsum(mean(rifwt)) /// pg here is just a dummy // stp1 all together?? No //all=att_gt,pg // stp2 get Mean(RIF) // This just rescales the IFs RIF's to make the statistics later. real matrix mean_y, exp_factor mean_y = colsum(rifgt):/colnonmissing(rifgt) _editmissing(mean_y,0) exp_factor = (rows(rifgt):/colnonmissing(rifgt)) rifgt = rifgt:-mean_y _editmissing(rifgt,0) _editmissing(exp_factor,0) rifgt =rifgt:*exp_factor:+mean_y mean_y =colsum(rifwt):/colnonmissing(rifwt) _editmissing(mean_y,0) exp_factor = (rows(rifwt):/colnonmissing(rifwt)) rifwt =rifwt:-mean_y _editmissing(rifwt,0) _editmissing(exp_factor,0) rifwt =rifwt:*exp_factor:+mean_y st_store(.,tokens(rifgt_+" "+rifwt_),(rifgt,rifwt)) /// Extra stuff t0 = strtoreal(tokens(time0)) glvl = strtoreal(tokens(glvl_)) tlvl = strtoreal(tokens(tlvl_)) real matrix ag_rif, ag_wt real matrix bb, VV, VV1, aux real vector ind_gt, ind_wt string matrix coleqnm ///////////////////////////////////////// // Always make attgt, even if not shown. //VV-VV1 if (wboot_==" ") { make_tbl( (rifgt,rifwt) ,bb,VV,clvar_, wboot_ ,VV1 , cband, ci, reps, wbtype, citype) st_matrix("b_attgt",bb) st_matrix("V_attgt",VV) } else { make_tbl( (rifgt,rifwt) ,bb,VV,clvar_, " " ,VV1, cband, ci, reps, wbtype,citype ) st_matrix("b_attgt",bb) st_matrix("V_attgt",VV) make_tbl( (rifgt,rifwt) ,bb,VV,clvar_, wboot_ ,VV1 , cband, ci, reps, wbtype,citype) } //if (wboot_!=" ") make_tbl( (rifgt,rifwt) ,bb,VV,clvar_, wboot_ , VV1) ///////////////////////////////////////// if (agg=="simple") { k=0 ind_gt=J(1,0,.) ind_wt=colsum(abs(rifgt)) for(i=1;i<=cols(glvl);i++) { for(j=1;j<=cols(tlvl);j++) { k++ // G <= T if ( (tlvl[j]-glvl[i]>=from) & (ind_wt[k]!=0) ) { ind_gt=ind_gt,k } } } // Above gets the Right elements Below, aggregates them ag_rif = rifgt[.,ind_gt] ag_wt = rifwt[.,ind_gt] aux = aggte(ag_rif, ag_wt) make_tbl(aux ,bb,VV,clvar_, wboot_ , VV1, cband, ci, reps, wbtype, citype) coleqnm = "ATT" } /////////////////////////////////////////simple pretrend real scalar dfchi, flag if (agg=="pretrend") { k=0 ind_gt=J(1,0,.) ind_wt=colsum(abs(rifgt)) for(i=1;i<=cols(glvl);i++) { for(j=1;j<=cols(tlvl);j++) { k++ if ( (glvl[i]>tlvl[j]) & (ind_wt[k]!=0) ) { //ag_rif=ag_rif, rifgt[.,k] ind_gt=ind_gt,k } } } real scalar nn2 ag_rif = rifgt[.,ind_gt] nn2 = rows(ag_rif)^2 bb = nn2 * mean(ag_rif) * invsym(crossdev(ag_rif,mean(ag_rif),ag_rif,mean(ag_rif))) * mean(ag_rif)' dfchi = cols(ag_rif) } /////////////////////////////////////////simple pretrend group if (agg=="group") { // i groups j time k=0 aux =J(rows(rifwt),0,.) sumwgt =J(rows(rifwt),0,.) coleqnm="GAverage" ind_wt=colsum(abs(rifgt)) /// ag_wt=J(rows(rifwt),0,.) for(i=1;i<=cols(glvl);i++) { ind_gt=J(1,0,.) flag=0 ag_rif=J(rows(rifwt),0,.) for(j=1;j<=cols(tlvl);j++) { k++ if ((tlvl[j]-glvl[i]>=from) & (ind_wt[k]!=0)) { //ag_rif=ag_rif, rifgt[.,k] flag=1 ind_gt=ind_gt,k } } if (flag==1) { coleqnm=coleqnm+sprintf(" G%s",strofreal(glvl[i])) ag_rif = rifgt[.,ind_gt] ag_wt = rifwt[.,ind_gt] sumwgt = sumwgt, rowsum(ag_wt):/cols(ag_wt) aux = aux, aggte(ag_rif, ag_wt) } } _editmissing(sumwgt,0) //# Bookmark #1 Uncertain if this is the best way to fix this, but right now gives best results sumwgt=J(rows(sumwgt),1,colsum(sumwgt)) aux = aggte(aux,sumwgt ), aux // get table elements make_tbl(aux ,bb,VV,clvar_, wboot_ ,VV1, cband, ci, reps, wbtype, citype) } ///////////////////////////////////////// if (agg=="calendar") { // i groups j time sumwgt=aux =J(rows(rifwt),0,.) coleqnm="CAverage " ind_wt=colsum(abs(rifgt)) for(h=1;h<=cols(tlvl);h++){ k=0 flag=0 ind_gt=J(1,0,.) /// ag_wt=J(rows(rifwt),0,.) for(i=1;i<=cols(glvl);i++) { for(j=1;j<=cols(tlvl);j++) { k++ if ( (tlvl[j]-glvl[i]>=from) & (tlvl[h]==tlvl[j]) & (ind_wt[k]!=0) ){ //ag_rif=ag_rif, rifgt[.,k] //ag_wt =ag_wt , rifwt[.,i] ind_gt=ind_gt,k //ind_wt=ind_wt,i if (flag==0) coleqnm=coleqnm+sprintf(" T%s",strofreal(tlvl[h])) flag=1 } } } if (flag==1) { ag_rif = rifgt[.,ind_gt] ag_wt = rifwt[.,ind_gt] ///sumwgt = sumwgt, rowsum(ag_wt):/cols(ag_wt) aux = aux, aggte(ag_rif, ag_wt) } } ///_editmissing(sumwgt,0) aux = aggte(aux, J(rows(aux),cols(aux),1) ), aux // get table elements make_tbl(aux ,bb,VV,clvar_, wboot_ ,VV1, cband, ci, reps, wbtype, citype ) } if (agg=="event") { // i groups j time real matrix evnt_lst, iit evnt_lst=event_list(glvl,tlvl,wndw) coleqnm="Pre_avg Post_avg " ind_wt=colsum(abs(rifgt)) sumwgt = aux =J(rows(rifwt),0,.) iit = J(1,0,.) if (lng==1) { lng=max( evnt_lst[1,2..length(evnt_lst)]-evnt_lst[1,1..(length(evnt_lst)-1)] ) } for(h=1;h<=cols(evnt_lst);h++){ k=0 flag=0 ind_gt=J(1,0,.) /// ag_wt=J(rows(rifwt),0,.) for(i=1;i<=cols(glvl);i++) { for(j=1;j<=cols(tlvl);j++) { k++ if ( ((glvl[i]+evnt_lst[h])==tlvl[j]) & (ind_wt[k]!=0) ) { //ag_rif=ag_rif, rifgt[.,k] //ag_wt =ag_wt , rifwt[.,i] ind_gt=ind_gt,k //ind_wt=ind_wt,i if (flag==0) { if (evnt_lst[h]< 0) coleqnm=coleqnm+sprintf(" Tm%s" ,strofreal(abs(evnt_lst[h]-lng))) if (evnt_lst[h]==0) coleqnm=coleqnm+" Tp0" if (evnt_lst[h]> 0) coleqnm=coleqnm+sprintf(" Tp%s",strofreal(abs(evnt_lst[h]))) } flag=1 } } } if (flag==1) { ag_rif = rifgt[.,ind_gt] ag_wt = rifwt[.,ind_gt] //sumwgt = sumwgt, rowsum(ag_wt):/rowsum(ag_wt:!=0) aux = aux, aggte(ag_rif, ag_wt ) iit = iit , evnt_lst[h]>=0 } } //_editmissing(sumwgt,0) //sumwgt = J(rows(aux),cols(aux),1) aux = aggte(select(aux,iit:==0), J(rows(aux),cols(aux)-sum(iit),1) ), aggte(select(aux,iit) , J(rows(aux),sum(iit),1) ), aux /// NEW line for Missing _editmissing(aux,0) // get table elements make_tbl(aux ,bb,VV,clvar_, wboot_ ,VV1, cband, ci, reps, wbtype, citype) } st_matrix(bb_,bb) st_matrix(ss_,VV) if (agg!="attgt") { stata("matrix colname "+bb_+" ="+coleqnm) stata("matrix colname "+ss_+" ="+coleqnm) stata("matrix rowname "+ss_+" ="+coleqnm) } } void make_tbl(real matrix rif,bb,VV, clv , wboot ,VV1 , string scalar cband_, real scalar ci, reps, wbtype, citype){ real matrix aux, nobs, clvar 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) clusterse((rif:-bb),clvar,VV) } real matrix cband // wboot no cluster if ( wboot!=" ") { mboot(rif,bb, VV, cband, clv, VV1, ci, reps, wbtype, citype) st_matrix(cband_, cband) } if (clv!=" ") { cln = rows(uniqrows(st_data(.,clv))) } st_numscalar("clnm",cln) } void clusterse(real matrix iiff, cl, V){ /// estimates Clustered Standard errors real matrix ord, xcros, ifp, info, vv //1st get the IFS and CL variable. //iiff = st_data(.,rif,touse) //cl = st_data(.,clvar,touse) // order and sort them, Make sure E(IF) is zero. ord = order(cl,1) //iiff = iiff:-mean(iiff) 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) // Esto es para ver como hacer clusters. //*nc/(nc-1) //st_matrix(V, vv) //st_numscalar(ncl, nc) // ^ ^ // | | // stata mata } void isbalanced(string ivar, touse, isbal) { real matrix xx, info xx= st_data(.,ivar,touse) xx=sort(xx,1) info=panelsetup(xx,1) info=info[,2]:-info[,1]:+1 st_numscalar(isbal,max(info)==min(info)) } void ispanel(string itvar, touse, ispan) { real matrix xx, info xx= st_data(.,itvar,touse) info=uniqrows(xx,1)[,3] st_numscalar(ispan,max(info)) } real colvector aggte(real matrix attg, wgt){ real scalar atte, mn_attg, mn_wgt real vector wgtw, attw real matrix r1, r2, r3 mn_attg = mean(attg) mn_wgt = mean(wgt) atte = sum(mn_attg:*mn_wgt):/sum(mn_wgt) wgtw = (mn_wgt ) :/sum(mn_wgt) attw = (mn_attg) :/sum(mn_wgt) r1 = (wgtw:*(attg:-mn_attg)) r2 = (attw:*(wgt :-mn_wgt )) r3 = (wgt :- mn_wgt) :* (atte :/ sum(mn_wgt) ) return(rowsum(r1):+rowsum(r2):-rowsum(r3):+atte) } ///////////////////////////////////////////////////////// real matrix mboot_did(real matrix rif, mean_rif, real scalar reps, bwtype) { 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 (bwtype==1) { for(i=1;i<=reps;i++){ bsmean[i,]=mean(yy:*(k2:-sqrt(5)*(rbinomial(n,1,1,k1))) ) } } else if (bwtype==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, real scalar reps, bwtype, clv) { real matrix yy, bsmean yy=rif:-mean_rif bsmean=J(reps,cols(yy),0) real scalar i,n, k1, k2, nn n=rows(yy) k1=((1+sqrt(5))/(2*sqrt(5))) k2=0.5*(1+sqrt(5)) real matrix sclv, wmult sclv=uniqrows(clv) nn=rows(sclv) if (bwtype==1) { for(i=1;i<=reps;i++){ wmult=(rbinomial(nn,1,1,k1)) bsmean[i,]=mean(yy:*(k2:-sqrt(5)*wmult[clv] ) ) } } else if (bwtype==2) { for(i=1;i<=reps;i++){ wmult=(rbinomial(nn,1,1,0.5)) bsmean[i,]=mean(yy:*(1:-2* wmult[clv] ) ) } } return(bsmean) } 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 qtp2(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) } 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) y=rowmax(y) for(i=1;i<=k;i++){ yy=sort(y,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, real matrix vv1, real scalar ci, reps, wbtype, citype) { real matrix fr, tt //real scalar reps, wbtype //reps = 999 //wbtype = 1 //ci = 0.95 real matrix ifse , ccb // this gets the Bootstraped values if (clv ==" ") { fr=mboot_did(rif,mean_rif, reps, wbtype) ifse = iqrse(fr) // this gets Tvalue //qtp(abs(fr :/ ifse),ci) //qtp2(abs(fr :/ ifse),ci) if (citype ==1) tt = qtp(abs(fr :/ ifse),ci) else if (citype ==2) tt = qtp2(abs(fr :/ ifse),ci) cband=( mean_rif', ifse', mean_rif':/ifse', mean_rif':-tt':* ifse' , mean_rif':+tt':* ifse' ) } else { clvar=st_data(.,clv) fr=mboot_didc(rif,mean_rif, reps, wbtype, clvar) ifse = iqrse(fr) if (citype ==1) tt = qtp(abs(fr :/ ifse),ci) if (citype ==2) tt = qtp2(abs(fr :/ ifse),ci) cband=( mean_rif', ifse', mean_rif':/ifse', mean_rif':-tt':* ifse' , mean_rif':+tt':* ifse' ) } vv=quadcross(ifse,ifse):*I(cols(ifse)) vv1=quadcross(fr,fr)/rows(fr) //sqrt(variance(fr)) //st_matrix(vv,iqrse(fr)^2) //st_matrix(cband,ccb) } end program addr, rclass return `0' end