diff --git a/csdid2 - Copy/README.md b/csdid2 - Copy/README.md new file mode 100644 index 0000000..91a6f4e --- /dev/null +++ b/csdid2 - Copy/README.md @@ -0,0 +1,94 @@ +# csdid2 +New Version of CSDID. All in Mata + +Hi, this is a new version of csdid that is now fully integrated into Mata. +Thus is faster a more veratile that before! + +You can see for yourself how it works: + +```stata +* loads data from a repository +ssc install frause +frause mpdta, clear +* This will generate everything, but show nothing! unless you request it. +* this can be done using the options agg(attgt) or agg(group) etc +. csdid2 lemp, ivar(countyreal) tvar(year) gvar(first) +Always Treated units have been excluded +----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 +............ +* However after that is done, you can just use estat to produce outcomes you want + +. estat event +------------------------------------------------------------------------------ + | Coefficient Std. err. z P>|z| [95% conf. interval] +-------------+---------------------------------------------------------------- + Pre_avg | .0018283 .007657 0.24 0.811 -.0131791 .0168357 + Post_avg | -.0772398 .019965 -3.87 0.000 -.1163705 -.0381092 + tm3 | .0305067 .0150336 2.03 0.042 .0010414 .0599719 + tm2 | -.0005631 .0132916 -0.04 0.966 -.0266142 .0254881 + tm1 | -.0244587 .0142364 -1.72 0.086 -.0523616 .0034441 + tp0 | -.0199318 .0118264 -1.69 0.092 -.0431111 .0032474 + tp1 | -.0509574 .0168935 -3.02 0.003 -.084068 -.0178468 + tp2 | -.1372587 .0364357 -3.77 0.000 -.2086713 -.0658461 + tp3 | -.1008114 .0343592 -2.93 0.003 -.1681542 -.0334685 +------------------------------------------------------------------------------ +* This produces Asymptotic Standard errors by default. But you can also reqyest bootstrap (no saverif anymore) + +. estat event, wboot +--------------------------------------------------------------------- + | Coefficient Std. err. t [95% conf. interval] +------------+-------------------------------------------------------- + Pre_avg | .0018283 .0076744 0.24 -.0172037 .0208603 + Post_avg | -.0772398 .020526 -3.76 -.1281428 -.0263368 + tm3 | .0305067 .0156614 1.95 -.0083324 .0693458 + tm2 | -.0005631 .013331 -0.04 -.0336231 .0324969 + tm1 | -.0244587 .0147201 -1.66 -.0609636 .0120462 + tp0 | -.0199318 .0116118 -1.72 -.0487284 .0088648 + tp1 | -.0509574 .0162214 -3.14 -.0911853 -.0107294 + tp2 | -.1372587 .0367764 -3.73 -.2284616 -.0460559 + tp3 | -.1008114 .0358325 -2.81 -.1896734 -.0119493 +--------------------------------------------------------------------- +WildBootstrap Standard errors +with 999 Repetitions + +** But the fun doesnt end there. You can also plot! +. estat event, wboot plot +[Plot not included but you can check it out] +``` + +Now the biggest difference (if you notice) is that all the IF information is kept in memory using Mata. +So, if you do something else, you may want to clean the created objects: + +``` +csdid2 , clear +``` + +Or, probably better, save the object in disk, so you can come back to the analysis If needed. + +``` +csdid2 save_ex1, save +clear all +csdid2 save_ex1, load +``` + +Of course, the only caveat you may have to use csdid2 to do the deed, instead of estat + +``` +csdid2 event, estat +``` + +Couple of notes. When using very large samples (many groups/periods/observations) you will need a lot of memory to keep it all in memory. +In general, you should plan to have atleast: + +(nobs x (# periods) x ( # groups) x 16 /1024^3) GB +of memory to store all the info needed. Which will stay in memory until you use `csdid2, clear`. + +For most applications that would be fine, but I recall some people using 10gb datasets, which may find problems. + +Last note! + +you need to copy all files in this repository on your ado/personal folder. And start from a new Stata session for it to work. +Fernando + + + diff --git a/csdid2 - Copy/compiler.do b/csdid2 - Copy/compiler.do new file mode 100644 index 0000000..349e900 --- /dev/null +++ b/csdid2 - Copy/compiler.do @@ -0,0 +1,9 @@ +cd "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2" +clear all + +run drdid.mata +run csdid.mata +run csdid_stats.mata +mata:mata mlib create lcsdid , replace +mata:mata mlib add lcsdid *() , complete + diff --git a/csdid2 - Copy/csdid.mata b/csdid2 - Copy/csdid.mata new file mode 100644 index 0000000..2352932 --- /dev/null +++ b/csdid2 - Copy/csdid.mata @@ -0,0 +1,530 @@ +*! v1.11 Bug with Cluster and Panel +* v1.1 Allows for new Estimators. +* v1 Allows for anticipation + +mata: +// This ones are to get the names of events and Other Stuff + void event_p( string scalar newvars, string scalar tblx){ + real matrix tbl, ntbl2 + string matrix ntbl + tbl = st_matrix(tblx) + + ntbl = st_matrixcolstripe(tblx) + ntbl = usubinstr(ntbl,"tp","+",.) + ntbl = usubinstr(ntbl,"tm","-",.) + ntbl2= strtoreal(ntbl) + tbl = tbl[(1,5,6),]' + tbl = select(tbl,(ntbl2[,2]:!=.)) + ntbl2= select(ntbl2[,2],(ntbl2[,2]:!=.)) + real matrix ss + ss= _st_addvar("double",tokens(newvars)) + st_store((1::rows(tbl)) ,tokens(newvars),(ntbl2,tbl)) + } + + void other_p(string scalar newvars, string scalar tblx){ + real matrix tbl + string matrix ntbl + tbl = st_matrix(tblx) + ntbl = st_matrixcolstripe(tblx) + //ntbl = usubinstr(ntbl,"g","",.) + //ntbl = usubinstr(ntbl,"t","",.) + ntbl = ntbl [,2] + tbl = tbl[(1,5,6),]' + string matrix tnv + tnv = tokens(newvars) + real matrix ss + ss= _st_addvar(sprintf("str%f",max(strlen(ntbl))),tnv[1]) + ss= _st_addvar("double",tnv[2..4]) + st_sstore((1::rows(tbl)) ,tnv[1],ntbl) + st_store((1::rows(tbl)) ,tnv[2..4],tbl) + } +end + + +mata +// Mata Class Gets all Relevant info for CSDID. +// Version 3 will use Sparse matrix +class csdid { + // input data + real matrix yvar + real matrix xvar + real matrix cvar + real matrix tvar,stvar + real matrix gvar,sgvar + real scalar delta + real matrix wvar + real matrix ivar + real matrix oid + real matrix foid + real matrix ord + // final RIF + real matrix frif + // weights + real matrix frwt + /// type of model, options + real scalar notyet + real scalar shortx + real scalar asinr + real scalar type_data + real scalar type_est + real scalar antici + /// Special Option for JW Rolling Regressopm + real scalar rolljw + + // output for selection + real matrix fgtvar + real matrix ogtvar + real matrix eventvar, sevent + // Did it converged? + real matrix convar + // useful functions + real matrix sample_select() + real matrix nvalid() + real matrix nsample_select() + + void gtvar() + void makeid() + void sevent() + void fixrif() + // model setup + void csdid() + void csdid_type() + void csdid_setup() + void setup_yvar() + void setup_xvar() + void setup_wvar() + void setup_cvar() + void setup_ivar() + void setup_tvar() + void setup_gvar() + //////////////////////////////////////////////////////////////////////////// + void easter_egg() +} + + +// This Module FIXES the missing records. +void csdid::easter_egg() { +if (runiform(1,1,0,1)<.5) { +" .-^-." +" .'=^=^='." +" /=^=^=^=^=\" +" .-~-. :^= HAPPY =^;" +" .'~~*~~'.|^ EASTER! ^|" +" /~~*~~~*~~\^=^=^=^=^=^:" +" :~*~~~*~~~*~\.-*))`*-,/" +" |~~~*~~~*~~|/* ((* *'." +" :~*~~~*~~~*| *)) * *\" +" \~~*~~~*~~| * ((* * /" +" '.~~*~~.' \ *)) * .'" +" '~~~' '-.((*_.-' " +} +else { +" .==." +" ()''()-." +" .---. ;--; /" +" .'_:___'. _..'. __'." +" |__ --==|'-''' \'...;" +" [ ] :[| |---\" +" |__| I=[| .' '." +" / / ____| : '._" +" |-/.____.' | : :" +" /___\ /___\ '-'._----'" +" ObiWan Kenobi-CSDID is our only hope" +} +} + +void csdid::fixrif(){ + real matrix mn_rif, rif2 + real scalar cnmiss + cnmiss = colnonmissing(frif) + mn_rif= colsum(frif):/cnmiss + frif = frif:-mn_rif + frif = editmissing(frif,0) + frif = mn_rif:+frif:*(rows(frif):/cnmiss) +} +/// Everything else will be used + +// Load ALL data +/* +1. get data +2. set method and data type +3. final tuning +*/ +/// SETS Data from Stata +/// yvar +void csdid::setup_yvar(string scalar ys, touse) st_view(yvar =.,.,ys,touse) +/// xvar, may be null +void csdid::setup_xvar(string scalar xs, touse) st_view(xvar =.,.,xs,touse) +/// ivar May be null +void csdid::setup_ivar(string scalar is, touse) st_view(ivar =.,.,is,touse) +/// tvar +void csdid::setup_tvar(string scalar ts, touse) st_view(tvar =.,.,ts,touse) +void csdid::setup_gvar(string scalar gs, touse) st_view(gvar =.,.,gs,touse) +// wvar may be null +void csdid::setup_wvar(string scalar ws, touse) st_view(wvar =.,.,ws,touse) +void csdid::setup_cvar(string scalar cs, touse) st_view(cvar =.,.,cs,touse) + +//Setting up model +void csdid::csdid_type(type_est, notyet, shrt, asinr){ + // notyet + this.notyet = notyet + // isshort + this.shortx = shrt + // as in r + this.asinr = asinr + // panel or RC + // this.type_data = type_data + // DRIPW DRIMP SIPW REG + this.type_est = type_est +} + + +void csdid::csdid_setup(){ + /// final Tunning to model + // oid. Original ID. Useful for RC. DROP + // IF RC + // antici for Anticipation + antici = 0 + rolljw = 0 + if (length(ivar)==0) { + // Repeated CrosSection + type_data = 2 + oid = 1::rows(yvar) + ivar = oid + // Ivar here does not matter, its just a sorting variable + if (length(cvar)>0) ord = order((cvar,gvar,tvar,ivar),(1,2,3,4)) + else ord = order((gvar,tvar,ivar),(1,2,3)) + // If data was already loaded sorted. No need to Resort + if (ord!=oid) { + yvar=yvar[ord,] + if (length(xvar)>0) xvar=xvar[ord,] + tvar=tvar[ord,] + gvar=gvar[ord,] + + if (length(wvar)>0) wvar=wvar[ord,] + + ivar=ivar[ord,] + + if (length(cvar)>0) cvar=cvar[ord,] + } + + } + else { + // if panel, first sort + type_data = 1 + oid = 1::rows(yvar) + if (length(cvar)>0) ord = order((cvar,ivar,tvar),(1,2,3)) + else ord = order((ivar,tvar),(1,2)) + + if (ord!=oid) { + // Order all variables in case order changed + yvar=yvar[ord,] + tvar=tvar[ord,] + gvar=gvar[ord,] + ivar=ivar[ord,] + // This runs only if they exist + if (length(xvar)>0) xvar=xvar[ord,] + if (length(wvar)>0) wvar=wvar[ord,] + if (length(cvar)>0) cvar=cvar[ord,] + + } + // then recode + makeid() + // if panel + foid = uniqrows(oid) + } + + if (length(wvar)==0) wvar=J(rows(yvar),1,1) +} +// Justs Puts all into Running Data +void csdid::makeid(){ + real scalar i,j, in + real matrix aux + in = rows(ivar) + oid=runningsum(0\(ivar[2..in,]-ivar[1..in-1,]):!=0):+1 + // recoded + //oid=uniqrows(id2,1)[id2,] +} + +// gtvar creates 2 things +// fgtvar <- Full set of Cohort year for regression +// ogtvar <- All combinations observed in data. To detect what is missing +// This IDS data and Creates Gvar + +void csdid::gtvar(){ + real matrix aux + // ids all cohorts + stvar=uniqrows(tvar) + sgvar=uniqrows(gvar) + // Check if there is overlap + if (length(uniqrows(stvar\sgvar))>=(length(stvar)+length(sgvar))) { + "There is no overlapping between Tvar and Gvar" + "Check to verify Gvar is correctly Specified" + exit(1) + } + + // If Never treated in Time window + /*if (max(sgvar):>max(stvar)) { + stata(`"display "Some Records were treated After the last period in data ""') + stata(`"display "They will be treated as Never treated""') + // Recode as Never Treated + gvar=gvar:*(gvar:<=max(stvar)) + sgvar=uniqrows(gvar) + }*/ + // if Treated After Max T. Never treated + //notyet -> Verify if there is any Not treated + if (notyet==0) notyet=min(sgvar:>0) + + sgvar=select(sgvar, ( (sgvar:>0) :* (sgvar:>min(stvar)) :* (sgvar:<=max(stvar)) ) ) + // delta time. Period Change + aux = uniqrows(sgvar\stvar) + + delta=min(aux[2..length(aux)]:- + aux[1..length(aux)-1]) + + // rescale antici to folow Delta. + antici=antici*delta + sgvar=sgvar:-antici + gvar=gvar:-antici:*(gvar:>0) + // recreate stvar + stvar=range(min(stvar),max(stvar),delta) + real matrix stvar2 + stvar2=stvar[2..length(stvar)] + // fullgtvar + // may not be necessary. And may be waste of space + // except for names! + fgtvar=sgvar#J(rows(stvar2),1,1), J(rows(sgvar),1,stvar2) + // fgtvar = gvar tvar : ATT(G,T) + ogtvar=uniqrows((gvar,tvar),1) + + // Next ID all good fgtvar + real matrix sel_gtvar + sel_gtvar=nsample_select() + + // Contents. + // gvar tvar t0 t1 te y00 y01 y10 y11 sel + // fgtvar ,sel_gtvar N Select + fgtvar = select((fgtvar,sel_gtvar),sel_gtvar[,cols(sel_gtvar)]) + convar = J(rows(fgtvar),1,0) + // IDs all events into eventvar + // uses short and long differences + sevent() + sevent=uniqrows(eventvar) + sgvar=uniqrows(fgtvar[,1]) + stvar=uniqrows(fgtvar[,2]) +} +// Define event for all FGTVAR +// IDS Event For Dynamic Effects +void csdid::sevent(){ + eventvar=J(rows(fgtvar),1,1) + eventvar = fgtvar[,2]-fgtvar[,1] + if (shortx==0) eventvar = fgtvar[,5]-fgtvar[,1] + +} +//!asd + +real matrix csdid::sample_select(real matrix gvtv) { + // Cohort Never + // As in R uses all Nontreated units up to point T. + // THat has a problem. Units are used as treatment and control all the time. + // Not the right way in my opinion + real matrix tsel, gsel + real scalar gv, tv0, tv1,tv + gv = gvtv[1] + tv = gvtv[2] + tv0 = gvtv[3] + tv1 = gvtv[4] + if (notyet==0) gsel = (gvar:==0 :| gvar:==gv) + else { + gsel = (gvar:==0 :| gvar:>max((gv,tv1)) :| gvar:==gv) + if ((asinr==1) & (tvtv0 :| gvar:==gv) + else gsel = (gvar:==0 :| gvar:>tv1 :| gvar:==gv) + } + } + + /// time Selection. Rolljw uses All pre-treatment. May need to do something about it. + if (rolljw==0) tsel = (tvar:==tv0 :| tvar:==tv1) + else if (rolljw==1) tsel = (tvar:<=tv0 :| tvar:==tv1) + + return(tsel:*gsel) +} + +// Based on Summary, check if a particular combination is valid +real matrix csdid::nvalid(real matrix sgtvar, + real scalar gv ){ + real scalar tv0, tv1 + tv0 = min(sgtvar[,2]) + tv1 = max(sgtvar[,2]) + + return( sum(sgtvar[,3]:*(sgtvar[,1]:!=gv):*(sgtvar[,2]:==tv0)), + sum(sgtvar[,3]:*(sgtvar[,1]:!=gv):*(sgtvar[,2]:==tv1)), + sum(sgtvar[,3]:*(sgtvar[,1]:==gv):*(sgtvar[,2]:==tv0)), + sum(sgtvar[,3]:*(sgtvar[,1]:==gv):*(sgtvar[,2]:==tv1)) ) +} + +real matrix csdid::nsample_select() { + real scalar i + real matrix tsel, gsel, sgtvar + real scalar gv, tv + real matrix toreturn + real matrix eftime + toreturn = J(rows(fgtvar),7,0) + // event + eftime=J(1,2,0) + + // This is to determine Obs + for(i=1;i<=rows(fgtvar);i++){ + gv = fgtvar[i,1];tv = fgtvar[i,2] + + if (notyet==0) gsel = (ogtvar[,1]:==0 :| ogtvar[,1]:==gv) + else { + gsel = (ogtvar[,1]:==0 :| ogtvar[,1]:>max((gv,tv)) :| ogtvar[,1]:==gv) + // Check asin R above. This is the same but with N (counting Obs) + if ((asinr==1) & (tvtv0 :| ogvar:==gv) + else gsel = (ogvar:==0 :| ogvar:>tv1 :| ogvar:==gv) + } + + } + /// time Selection + if (tv>=gv) { + tsel = (ogtvar[,2]:==gv-delta :| ogtvar[,2]:==tv) + /// T0 T1 T + eftime=( gv-delta ,tv , tv) + } + else { + if (shortx==0 ) { + tsel = (ogtvar[,2]:==gv-delta :| ogtvar[,2]:==tv-delta) + /// T0 T1 T + eftime=( tv-delta ,gv-delta, tv-delta ) + } + else { + tsel = (ogtvar[,2]:==tv-delta :| ogtvar[,2]:==tv) + /// T0 T1 T + eftime=(tv-delta ,tv , tv-delta) + } + } + // Up to here we "select" sample + // next, we need to see if sample is feasible based on numbers + + sgtvar = select(ogtvar ,tsel:*gsel) + + toreturn[i,] = eftime, nvalid(sgtvar,gv) + } + + toreturn=toreturn,rowmin(toreturn:>0) + + return(toreturn) + // T is used for Event . Capture relative numbers. + // Otherwise use t0 t1 fgtvar34 +} + + +void csdid::csdid(){ + //frif=oid,gvar,wvar + class drdid scalar drdid + real matrix smsel + real scalar gv, tv, dots, i + // To make it faster + real matrix index1 + + frwt=frif=J( ((type_data==1) ? max(oid) : rows(oid)) , + rows(fgtvar),.) + index1 = 1::rows(oid) + + stata("_dots 0") + drdid.rolljw=rolljw + + for(i=1;i<=rows(fgtvar);i++) { + dots = 1 + drdid.init() + + //smsel = sample_select(fgtvar[i,]) + smsel = select(index1, sample_select(fgtvar[i,])) + + gv = fgtvar[i,1]; tv = fgtvar[i,4] + + ///minn = min(fgtvar[i,6..9]) + drdid.yvar=yvar[smsel,] + //if (cols(xvar)>0) drdid.xvar=select(xvar,smsel) + if (cols(xvar)>0) drdid.xvar=xvar[smsel,] + //drdid.tmt =select(tvar,smsel):==tv + //drdid.trt =select(gvar,smsel):==gv + drdid.tmt =tvar[smsel,]:==tv + drdid.trt =gvar[smsel,]:==gv + //if (cols(wvar)>0) drdid.wvar=select(wvar,smsel) + if (cols(wvar)>0) drdid.wvar=wvar[smsel,] + //drdid.id =select(ivar,smsel) + //drdid.oid =select(oid ,smsel) + drdid.id =ivar[smsel,] + drdid.oid =oid[smsel,] + drdid.data_type = type_data + drdid.method_type = type_est + drdid.conv=1 + drdid.drdid() + + /// Stores RIFS + if ((drdid.conv==1) & sum(abs(drdid.rif))>0 ) { + convar[i,]=1 + if (shortx==0) frif[drdid.oid,i]=drdid.rif:*sign(eventvar[i]+.01) + else frif[drdid.oid,i]=drdid.rif + frwt[drdid.oid,i]=drdid.wvar + //:*drdid.trt + + } + dots = dots-convar[i,] + stata(sprintf("_dots %f %f",i,dots)) + } + + /// fixes missing in rifs + /// Next version will try to avoid this + fixrif() + /// and will avoid having a copyu of weights. + /// This weights are simple to use. Not affected by RW + _editmissing(frwt,0) + /// Extra clean up + frwt = select(frwt , convar') + frif = select(frif , convar') + eventvar = select(eventvar, convar) + fgtvar = select(fgtvar , convar) + convar = select(convar , convar) + + //ConVar Converged Variable + + /// cleaning all else + // TRT not needed + ord=tvar=yvar=J(0,0,.) + /// One Risk. Missing data after drdid or else. + // if panel. we keep one per + if (type_data==1) { + real matrix aux + + if (length(cvar)>0) { + aux = oid, gvar, wvar, cvar + aux=uniqrows(aux) + ord = order(aux,(4,1,2,3)) + aux = aux[ord,] + oid = aux[,1] + gvar= aux[,2] + wvar= aux[,3] + cvar= aux[,4] + frwt = frwt[ord,] + frif = frif[ord,] + } + else { + aux = oid, gvar, wvar + aux=uniqrows(aux) + oid = aux[,1] + gvar= aux[,2] + wvar= aux[,3] + } + + } + + /// Very last step. Sort important variables by Cvar? + + ord=aux = J(0,0,.) +} + +end diff --git a/csdid2 - Copy/csdid2.ado b/csdid2 - Copy/csdid2.ado new file mode 100644 index 0000000..3f8aa78 --- /dev/null +++ b/csdid2 - Copy/csdid2.ado @@ -0,0 +1,369 @@ +*! v1.3 Allows for Reg2 and drimp2: hdidregress compatible +*! v1.21 Allows for treatvar +*! v1.2 Allows for Anticipation +*! v1.13 Cluster Corrections for RCS and RC +*! v1.12 Adds Option for CSname +*! v1.11 Corrects For missing in cluster +*! v1.1 adds Rolling +*! v1 Wrapper for CSDID2-Mata version + +program csdid2, sortpreserve eclass + version 14 + ///version checker + syntax [anything(everything)] [iw aw pw], [* version estat /// + save clear load replace plot csname(name)] + if "`save'`clear'`load'`replace'"!="" { + + exit + } + + if "`estat'"!="" { + csdid2_estat `anything', `options' + exit + } + + if "`plot'"!="" { + csdid2_plot, `options' + exit + } + + if "`version'"!="" { + display "version: 1" + addr scalar version = 1 + exit + } + + ///Replay + + if replay() { + if `"`e(cmd)'"' != "csdid2" { + error 301 + } + else { + Display `0' + } + exit + } + if runiform()<.001 { + easter_egg + } + + csdid_r `0' + +end + +program define easter_egg + display "This Easter Egg is Broken, but try not to get lost" +end + + + mata + + void balance(string scalar ivar, + string scalar tvar, + string scalar touse){ + real matrix i, t + real scalar ntt, ni, nt, max + st_view(i=.,.,ivar,touse) + st_view(t=.,.,tvar,touse) + ntt=rows(i) + ni=rows(uniqrows(i)) + nt=rows(uniqrows(t)) + max = max(uniqrows((i,t),1)[,3]) + if (max>1) { + _error(451, "repeated time values within panel") + } + if ( (ni*nt) > ntt) { + stata(`"display in red "Panel is not balanced""') + stata(`"display in red "Will use observations with Pair balanced (observed at t0 and t1)""') + } + + } +end + +program _gcsgvar, sortpreserve + syntax newvarname =/exp [if] [in], tvar(varname) ivar(varname) + local exp = subinstr("`exp'","(","",.) + local exp = subinstr("`exp'",")","",.) + tempvar touse + qui:gen byte `touse'=0 + qui:replace `touse'=1 `if' `in' + qui:replace `touse'=0 if `tvar'==. | `ivar'==. | `exp'==. + tempvar vals + sum `exp' if `touse' , meanonly + local lmin=r(min) + local lmax=r(max) + + if `lmin'<0 | `lmax'>1 { + display in r "`exp' can only have values between 0 and 1" + error 4444 + } + qui: { + tempvar aux auxd + qui: gen byte `auxd'=`exp'>0 if `exp'!=. + bysort `touse' `ivar' `auxd':egen `aux'=min(`tvar') + replace `aux'=0 if `exp'==0 + by `touse' `ivar':egen `varlist'=max(`aux') + replace `varlist'=. if `exp'==. | !`touse' + } + label var `varlist' "Group Variable based on `exp'" +end + program csdid_r, sortpreserve eclass + syntax varlist(fv ) [if] [in] [aw pw iw/], /// Basic syntax allows for weights + [Ivar(varname numeric) group(varname)] /// + [TIme(varname numeric)] /// + [Tvar(varname numeric)] /// + [Gvar(varname numeric) trtvar(varname)] /// + [cluster(varname numeric) /// + notyet /// att_gt basic option. May prepare others as Post estimation + method(str) /// + long long2 /// to allow for "long gaps" + asinr /// For pretreatment + agg(string) /// type of aggregation + rolljw /// Rolling Regression Estimator + csname(name) /// To change name to csobj + short /// for compatability + antici(int 0) /// anticipation + ] + ** Checking if Setup + local issetup:char _dta[csdid2] + if "`issetup'"=="csdid2" { + local tvar:char _dta[tvar] + local gvar:char _dta[gvar] + local ivar:char _dta[ivar] + } + ** Marking Sample + if "`csname'"=="" local csname csdid + marksample touse + ** First determine outcome and xvars + gettoken y xvar:varlist + ** Is time set? + if "`time'`tvar'"=="" { + display in red "Option tvar() or time() is required" + error 198 + } + else { + // Tvar will superseed time + if "`tvar'"=="" local tvar `time' + } + ** is gvar or trtvar set? + if "`gvar'`trtvar'"=="" { + display in red "option gvar/trtvar() required" + error 198 + } + else if "`trtvar'"!="" & "`gvar'"!="" { + display as error "You can only specify gvar or trtvar. Not both" + error 198 + } + else if "`trtvar'"!="" { + capture drop __gvar + if "`ivar'`group'"=="" display "You need to specifcy {it: ivar()} for panel or {it: group()} for Repeated Crossection" + else { + if "`ivar'"!="" local iivar "`ivar'" + else local iivar "`group'" + + } + qui:_gcsgvar __gvar=`trtvar', tvar(`tvar') ivar(`iivar') + local gvar __gvar + } + *** OLD for SHORT + if "`short'"=="" { + local long long + display "Producing Long Gaps by default" + } + else { + local long + } + + + markout `touse' `ivar' `tvar' `gvar' `y' `xvar' `cluster' + local wvar `exp' + + ** Aggregation + if "`agg'"!="" & inlist("`agg'","simple","attgt","event","calendar","group") { + if "`agg'"=="simple" local aggtype 2 + if "`agg'"=="attgt" local aggtype 1 + if "`agg'"=="event" local aggtype 5 + if "`agg'"=="calendar" local aggtype 4 + if "`agg'"=="group" local aggtype 3 + } + + ** Verifying basics + if "`method'"=="" { + local method dripw + if "`xvar'"=="" { + local method reg + } + } + else if !inlist("`method'","drimp","dripw","reg","reg2","drimp2","stdipw") { + display in red "Method `method' not allowed" + error 1 + } + ** What method is being used + display "Using method `method'" + + ** Always Treated Excluded + qui: sum `touse', meanonly + local pre_mean `r(mean)' + sum `tvar' if `touse', meanonly + qui:replace `touse'=0 if `gvar'<`r(min)' & `gvar'>0 + + qui: sum `touse', meanonly + local post_mean `r(mean)' + + if `pre_mean'!=`post_mean' display "Always Treated units have been excluded" + + ** is gvar nested iwthing county + if "`ivar'"!="" { + _xtreg_chk_cl2 `gvar' `ivar' `touse' + } + ** is cluster correct? + if "`cluster'"!="" { + if "`ivar'"!="" _xtreg_chk_cl2 `cluster' `ivar' `touse' + ** else _xtreg_chk_cl2 `gvar' `cluster' `touse' + } + ** Is a balanced panel? + if "`ivar'"!="" { + mata:balance("`ivar'","`tvar'","`touse'") + } + *** Here we SETUP CSDID + // Create the Mata Object + local cvar `cluster' + mata:`csname'=csdid() + mata:`csname'.setup_yvar("`y'" ,"`touse'") + mata:`csname'.setup_tvar("`tvar'" ,"`touse'") + mata:`csname'.setup_gvar("`gvar'" ,"`touse'") + if "`xvar'"!="" mata:`csname'.setup_xvar("`xvar'" ,"`touse'") + if "`ivar'"!="" mata:`csname'.setup_ivar("`ivar'" ,"`touse'") + if "`wvar'"!="" mata:`csname'.setup_wvar("`wvar'" ,"`touse'") + if "`cvar'"!="" mata:`csname'.setup_cvar("`cvar'" ,"`touse'") + ** Setup as panel or rc + mata:`csname'.csdid_setup() + ** + ** For Anticipation + mata:`csname'.antici=`antici' + + ** Rolling Regression + if "`rolljw'"!="" mata:`csname'.rolljw = 1 + + // Type_est 1 dripw 2 drimp 3 stipw 4 reg + // not_yet 0 Never 1 Notyet + // shrt 1 short 0 long + // asinr 0 stata 1 R + + if "`method'"=="dripw" local type_est=1 + else if "`method'"=="drimp" local type_est=2 + else if "`method'"=="stdipw" local type_est=3 + else if "`method'"=="reg" local type_est=4 + else if "`method'"=="reg2" local type_est=5 + else if "`method'"=="drimp2" local type_est=6 + + + + local ntyet = 0 + if "`tyet'"!="" local ntyet = 1 + local shrt = 1 + + if "`long'"!="" local shrt = 0 + local asr= 0 + if "`asinr'"!="" local asr = 1 + // Setup CSDID Type of estimation. Allowing Reg2 and Drimp2 + mata:`csname'.csdid_type(`type_est',`ntyet',`shrt',`asr') + // Get all GTs + mata:`csname'.gtvar() + // Estimate Model + mata:`csname'.csdid() + // Done + ereturn clear + ereturn local cmd csdid2 + ereturn local cmdline csdid2 `0' + ereturn local estat_cmd csdid2_estat + ereturn local method `method' + ereturn local asinr `asinr' + + if "`aggtype'"!="" { + mata:csdidstat=csdid_estat() + mata:csdidstat.init() + mata:csdidstat.test_type=`aggtype' + mata:csdidstat.atts_asym(`csname') + ereturn post _bb _vv + ereturn local vcetype Robust + ereturn local cmd csdid2 + ereturn local cmdline csdid2 `0' + ereturn local estat_cmd csdid2_estat + ereturn local method `method' + ereturn local asinr `asinr' + if "`cvar'"!="" ereturn local cluster `cvar' + Display + } + + + if `ntyet'==1 ereturn local cntrl "Not yet treated" + else ereturn local cntrl "Never treated" + if `shrt'==1 ereturn local base "Varying Base" + else ereturn local base "Base Universal" +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)'"=="reg") { + local tmodel "none" + local omodel "regression adjustment" + } + if ("`e(method)'"=="reg2") { + local tmodel "none" + local omodel "regression adjustment: RC2" + } + if ("`e(method)'"=="drimp2") { + local tmodel "inverse probability tilting" + local omodel "WLS: From RC2 reg2" + } + if ("`e(method)'"=="stdipw") { + local tmodel "stabilized inverse probability" + local omodel "weighted mean" + } + sreturn local omodel "`omodel'" + sreturn local tmodel "`tmodel'" +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' + } + + +end + diff --git a/csdid2 - Copy/csdid2_clean.ado b/csdid2 - Copy/csdid2_clean.ado new file mode 100644 index 0000000..aaca977 --- /dev/null +++ b/csdid2 - Copy/csdid2_clean.ado @@ -0,0 +1,23 @@ +*! v 1.1 adds csname +* v 1 Save and clear +program csdid2_clean + syntax [namelist(max=1 id="Names")], [save clear load replace dir(string asis) /// + csname(name) ] + if "`csname'"=="" local csname csdid + if "`clear'"!="" { + mata:mata drop `csname' + capture: mata:mata drop csdidstat + } + else if "`save'"!="" { + local old `"`c(pwd)'"' + if `"`dir'"' !="" qui: cd `dir' + mata:mata matsave `namelist' `csname', `replace' + qui: cd `"`old'"' + } + else if "`load'"!="" { + local old `"`c(pwd)'"' + if `"`dir'"' !="" qui: cd `dir' + mata:mata matuse `namelist', `replace' + qui: cd `"`old'"' + } +end \ No newline at end of file diff --git a/csdid2 - Copy/csdid2_estat.ado b/csdid2 - Copy/csdid2_estat.ado new file mode 100644 index 0000000..438e72f --- /dev/null +++ b/csdid2 - Copy/csdid2_estat.ado @@ -0,0 +1,336 @@ +*! v1.3 Correct with group +*! v1.2 for anticipation +*! v1.1 adds option for CSname +*! v1. FRA +* Allows plot with extra stuff + +program define clrreturn, rclass + exit +end + +program define adde, eclass + ereturn `0' +end + +program define addr, rclass + return add + return `0' +end + +program define adds, sclass + sreturn `0' +end + + program define csdid2_estat, sortpreserve + version 14 + syntax anything, [* plot csname(name)] + if "`csname'"=="" local csname csdid + capture mata:`csname' + if _rc!=0 error 301 + gettoken key rest : 0, parse(", ") + if `"`rest'"'=="" local rest ", csname(`csname')" + else local rest `rest' csname(`csname') + + if inlist("`key'","attgt","simple","pretrend","group","calendar","event","cevent") { + csdid_do `key' `rest' + addr local cmd estat + addr local cmd2 csdid2 + if "`plot'"!="" csdid2_plot , ktype(`ktype') + } + else if inlist("`key'","plot") { + csdid2_plot, `options' + } + else { + display in red "Option `key' not recognized" + error 199 + } + +end + + program define csdid_do, rclass + syntax namelist(min=1 max=1 ), [ estore(name) /// + esave(name) replace /// + post level(int 95) /// + WBOOT /// + WBOOT1(str) /// + reps(integer 999) /// + rseed(string) /// + wbtype(string) /// + rgroup(numlist) /// + rcalendar(numlist) /// + revent(numlist) /// + window(numlist min=2 max=2) /// + REBALance(numlist) /// <-- restricts groups and event, unless event is used too + max_mem(real 1) /// + noavg /// + csname(name) /// + plot * ] + + // confirm csdid exists and if csdidstat=csdid_estat() + local key `namelist' + capture mata:csdidstat + if _rc!=0 mata:csdidstat=csdid_estat() + + // check Keys + + if "`key'"=="pretrend" local ktype = 0 + if "`key'"=="attgt" local ktype = 1 + if "`key'"=="simple" local ktype = 2 + if "`key'"=="group" local ktype = 3 + if "`key'"=="calendar" local ktype = 4 + if "`key'"=="event" local ktype = 5 + if "`key'"=="cevent" local ktype = 6 + + /// initialize + if "`rseed'"!="" set seed `rseed' + local nov 0 + if "`avg'"!="" local nov 1 + + mata: csdidstat.cilevel = `level'/100 + mata: csdidstat.bwtype = 1 + mata: csdidstat.noavg = `nov' + mata: csdidstat.reps = `reps' + mata: csdidstat.max_mem = `max_mem' + mata: csdidstat.range.selgvar = J(0,0,.) + mata: csdidstat.range.seltvar = J(0,0,.) + mata: csdidstat.range.selevent= J(0,0,.) + mata: csdidstat.range.selbal = J(0,0,.) + + /// Check if we have to make sample selection + mata:st_local("adj",strofreal(`csname'.antici)) + + + if "`rgroup'"!="" { + numlist "`rgroup'", int + mata:csdidstat.range.selgvar=csdidstat.rtokens("`r(numlist)'"):-`adj' + } + if "`rcalendar'"!="" { + numlist "`rcalendar'", int + mata:csdidstat.range.seltvar=csdidstat.rtokens("`r(numlist)'") + } + if "`revent'"!="" { + numlist "`revent'", int + mata:csdidstat.range.selevent=csdidstat.rtokens("`r(numlist)'"):-`adj' + } + if "`window'"!="" { + numlist "`window'", min(2) max(2) sort integer + local window `r(numlist)' + local n1: word 1 of `window' + local n2: word 2 of `window' + numlist "`n1'/`n2'", int + mata:csdidstat.range.selevent=csdidstat.rtokens("`r(numlist)'"):-`adj' + } + if "`rebalance'"!="" { + numlist "`rebalance'", int + mata:csdidstat.range.selbal=csdidstat.rtokens("`r(numlist)'"):-`adj' + } + + + if `ktype'>0 mata: csdidstat.test_type = `ktype' + else { + mata: csdidstat.pretrend(`csname') + display "Pre-trend test" + display "H0: All ATTGT=g for all T`wdt' local wdt = length("`i'")+3 + } + if `wdt'<15 local wdt = 12 +*** + tempname mytab b se z t ll ul cimat rtab + .`mytab' = ._tab.new, col(6) lmargin(0) + .`mytab'.width `wdt' |12 12 8 12 12 + .`mytab'.titlefmt . . . %6s %24s . + .`mytab'.pad . 2 1 0 3 3 + .`mytab'.numfmt . %9.0g %9.0g %7.2f %9.0g %9.0g + + + local stat t + + local namelist : rowname `tablex' + local eqlist : roweq `tablex' + local k : word count `namelist' + local knew = `k' + matrix `rtab' = J(9, `k', .) + matrix `cimat'= `tablex' + * pvalue + matrix rownames `rtab' = b se t p ll ul df crit eform + matrix colnames `rtab' = `namelist' + forvalues i = 1/`k' { + local kxc: word `i' of `eqlist' + if ("`kxc'"=="wgt") { + local knew = `knew' -1 + } + matrix `rtab'[1,`i'] = `cimat'[`i',1] + matrix `rtab'[2,`i'] = `cimat'[`i',2] + matrix `rtab'[3,`i'] = `cimat'[`i',3] + matrix `rtab'[5,`i'] = `cimat'[`i',4] + matrix `rtab'[6,`i'] = `cimat'[`i',5] + matrix `rtab'[8,`i'] = `cimat'[`i',6] + } + .`mytab'.sep, top + if `:word count `e(depvar)'' == 1 { + local depvar "`e(depvar)'" + } + .`mytab'.titles "`depvar'" /// 1 + " Coefficient" /// 2 + "Std. err." /// 3 + "`stat'" /// 4 "P>|`stat'|" /// 5 + "[`level'% conf. interval]" "" // 6 7 + + forvalues i = 1/`knew' { + local name : word `i' of `namelist' + local eq : word `i' of `eqlist' + if ("`eq'" != "_") { + if "`eq'" != "`eq0'" { + .`mytab'.sep + local eq0 `"`eq'"' + .`mytab'.strcolor result . . . . . + .`mytab'.strfmt %-12s . . . . . + .`mytab'.row "`eq'" "" "" "" "" "" + .`mytab'.strcolor text . . . . . + .`mytab'.strfmt %12s . . . . . + } + local beq "[`eq']" + } + else if `i' == 1 { + local eq + .`mytab'.sep + } + + scalar `b' = `cimat'[`i',1] + scalar `se' = `cimat'[`i',2] + scalar `t' = `cimat'[`i',3] + scalar `ll' = `cimat'[`i',4] + scalar `ul' = `cimat'[`i',5] + .`mytab'.row "`name'" /// + `b' /// + `se' /// + `t' /// `p' /// + `ll' `ul' + } + .`mytab'.sep, bottom + return matrix table = `rtab' +end diff --git a/csdid2 - Copy/csdid2_plot.ado b/csdid2 - Copy/csdid2_plot.ado new file mode 100644 index 0000000..0d1e330 --- /dev/null +++ b/csdid2 - Copy/csdid2_plot.ado @@ -0,0 +1,323 @@ +*! v1.1 csdid2_plot Better names when Date has formats +* v1.03 csdid2_plot Allows for Anticpiation +* v1.02 csdid2_plot Fixes `' +* v1.01 csdid2_plot for csdid2 only + +program csdid2_plot, rclass + + syntax, [* ktype(string)] + local cmd `r(cmd)' + local agg `r(agg)' + tempname b V table bb vv + matrix `b' = r(b) + matrix `V' = r(V) + matrix `table' = r(table) + + if "`agg'"=="group" local ktype = 3 + if "`agg'"=="calendar" local ktype = 4 + if "`agg'"=="event" local ktype = 5 + + + capture noisily csdid2_plot_wh, `options' ktype(`ktype') table(`table') + return local cmd `cmd' + return local agg `agg' + + return matrix b =`b' + return matrix V =`V' + return matrix table =`table' + +end + +program csdid2_plot_wh + syntax, [ * ktype(int 5) table(str) asy level(int 95)] + + tempname tbl + matrix `tbl'=`table' + capture: confirm matrix `tbl' + + *if det(`tbl')==. matrix `tbl'=rtb + if "`asy'"=="" { + if `ktype'==5 { + tempvar t b ll uu + + mata:event_p("`t' `b' `ll' `uu'","`tbl'") + + csdid_plot_eventx `t' `b' `ll' `uu', `options' + } + else if `ktype'==3 | `ktype'==4 { + // Group Calendar + tempvar t b ll uu + mata:other_p("`t' `b' `ll' `uu'","`tbl'") + + csdid_plot_other `t' `b' `ll' `uu', `options' ktype(`ktype') + } + else { + display in red "Plot option not allowed" + } + } + else { + if `ktype'==5 { + tempvar t b ll uu se all auu + + mata:event_p2("`t' `b' `ll' `uu' `se' `all' `auu'","`tbl'", `=`level'/100') + csdid_plot_eventx `t' `b' `ll' `uu' `all' `auu', `options' asy + } + else if `ktype'==3 | `ktype'==4 { + // Group Calendar + tempvar t b ll uu se all auu + mata:other_p2("`t' `b' `ll' `uu' `se' `all' `auu'","`tbl'", `=`level'/100') + csdid_plot_other `t' `b' `ll' `uu' `all' `auu', `options' ktype(`ktype') asy + } + else { + display in red "Plot option not allowed" + } + } + +end + +mata: + void event_p2( string scalar newvars, string scalar tblx, real scalar level){ + real matrix tbl, ntbl2 + string matrix ntbl + real scalar alpha + tbl = st_matrix(tblx) + //asume always here + + alpha = 1-level + ntbl = st_matrixcolstripe(tblx) + ntbl = usubinstr(ntbl,"tp","+",.) + ntbl = usubinstr(ntbl,"tm","-",.) + ntbl2= strtoreal(ntbl) + tbl = tbl[(1,5,6,2),]' + tbl = tbl,tbl[,1]-tbl[,4]*invnormal(1-alpha/2),tbl[,1]+tbl[,4]*invnormal(1-alpha/2) + + tbl = select(tbl,(ntbl2[,2]:!=.)) + ntbl2= select(ntbl2[,2],(ntbl2[,2]:!=.)) + real matrix ss + ss= _st_addvar("double",tokens(newvars)) + st_store((1::rows(tbl)) ,tokens(newvars),(ntbl2,tbl)) + } + + void other_p2(string scalar newvars, string scalar tblx, real scalar level){ + real matrix tbl + string matrix ntbl + real scalar alpha + alpha = 1-level + tbl = st_matrix(tblx) + ntbl = st_matrixcolstripe(tblx) + //ntbl = usubinstr(ntbl,"g","",.) + //ntbl = usubinstr(ntbl,"t","",.) + ntbl = ntbl [,2] + tbl = tbl[(1,5,6,2),]' + tbl = tbl,tbl[,1]-tbl[,4]*invnormal(1-alpha/2),tbl[,1]+tbl[,4]*invnormal(1-alpha/2) + + string matrix tnv + tnv = tokens(newvars) + real matrix ss + ss= _st_addvar(sprintf("str%f",max(strlen(ntbl))),tnv[1]) + ss= _st_addvar("double",tnv[2..7]) + + st_sstore((1::rows(tbl)) ,tnv[1],ntbl) + st_store((1::rows(tbl)) ,tnv[2..7],tbl) + } +end + +program csdid2_default, sclass + syntax, [style(str) PSTYle1(str) color1(str) /// + PSTYLE2(str) color2(str) /// + LWidth1(str) lwidth2(str) /// + BARWidth1(str) barwidth2(str) * asy] + + if "`style'"=="" local style rspike + + + if "`pstyle1'"=="" local pstyle1 pstyle(p1) + else local pstyle1 pstyle(`pstyle1') + if "`pstyle2'"=="" local pstyle2 pstyle(p2) + else local pstyle2 pstyle(`pstyle2') + + if `"`color1'"'=="" local color1 color(%40) + else local color1 color(`"`color1'"') + if `"`color2'"'=="" local color2 color(%40) + else local color2 color(`"`color2'"') + + if "`style'"=="rspike" { + if "`lwidth1'"=="" local lwidth1 lwidth(3) + else local lwidth1 lwidth(`lwidth1') + if "`lwidth2'"=="" local lwidth2 lwidth(3) + else local lwidth2 lwidth(`lwidth2') + if "`asy'"!="" { + local lwidth1 lwidth(1) + local lwidth2 lwidth(1) + } + } + + if "`style'"=="rarea" { + if "`lwidth1'"=="" local lwidth1 lwidth(0) + else local lwidth1 lwidth(`lwidth1') + if "`lwidth2'"=="" local lwidth2 lwidth(0) + else local lwidth2 lwidth(`lwidth2') + local conn connect(l) + } + + if "`style'"=="rcap" { + if "`lwidth1'"=="" local lwidth1 lwidth(1) + else local lwidth1 lwidth(`lwidth1') + if "`lwidth2'"=="" local lwidth2 lwidth(1) + else local lwidth2 lwidth(`lwidth2') + + local conn connect(l) + } + + if "`style'"=="rbar" { + if "`lwidth1'"=="" local lwidth1 lwidth(0) + else local lwidth1 lwidth(`lwidth1') + if "`lwidth2'"=="" local lwidth2 lwidth(0) + else local lwidth2 lwidth(`lwidth2') + if "`barwidth1'"=="" local barwidth1 barwidth(0.5) + else local barwidth1 barwidth(`barwidth1') + if "`barwidth2'"=="" local barwidth2 barwidth(0.5) + else local barwidth2 barwidth(`barwidth2') + local conn connect(l) + } + + + sreturn local style `style' + sreturn local df11 `pstyle1' `color1' `lwidth1' `barwidth1' + sreturn local df12 `pstyle1' `conn' + sreturn local df21 `pstyle2' `color2' `lwidth2' `barwidth2' + sreturn local df22 `pstyle2' `conn' + sreturn local delse `options' +end + +program csdid_plot_eventx + syntax varlist, [ xtitle(passthru) ytitle(passthru) /// + legend(passthru) asy * ] + gettoken t rest:varlist + gettoken b rest:rest + gettoken ll rest:rest + gettoken uu rest:rest + + + ** defaults + + if `"`xtitle'"'=="" local xtitle xtitle("Periods to Treatment") + if `"`ytitle'"'=="" local ytitle ytitle("ATT") + + if `"`legend'"'=="" local legend legend(order(1 "Pre-treatment" 3 "Post-treatment")) + csdid2_default , `options' `asy' + + local gf11 `s(df11)' + local gf12 `s(df12)' + local gf21 `s(df21)' + local gf22 `s(df22)' + local style `s(style)' + local dels `s(delse)' + + mata:st_local("adj",strofreal(csdid.antici)) + if "`asy'"=="" { + two (`style' `ll' `uu' `t' if `t'<=(-1- `adj'), `gf11') /// + (scatter `b' `t' if `t'<=(-1- `adj'), `gf12') /// + (`style' `ll' `uu' `t' if `t'> (-1- `adj'), `gf21') /// + (scatter `b' `t' if `t'> (-1- `adj'), `gf22') , /// + `legend' `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `dels' + } + else { + //gettoken se rest:rest + gettoken all rest:rest + gettoken auu rest:rest + + two (`style' `ll' `uu' `t' if `t'<=(-1- `adj'), `gf11') /// + (scatter `b' `t' if `t'<=(-1- `adj'), `gf12') /// + (`style' `ll' `uu' `t' if `t'> (-1- `adj'), `gf21') /// + (scatter `b' `t' if `t'> (-1- `adj'), `gf22') /// + (`style' `all' `auu' `t' if `t'<=(-1- `adj'), `gf11' lwidth(3)) /// + (`style' `all' `auu' `t' if `t'> (-1- `adj'), `gf21' lwidth(3)), /// + `legend' `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `dels' + + } +end + + +program csdid_plot_other + syntax varlist, [ktype(int 3) * /// + xtitle(passthru) ytitle(passthru) /// + asy legend(passthru) format(passthru)] + gettoken t rest:varlist + gettoken b rest:rest + gettoken ll rest:rest + gettoken uu rest:rest + + + local xlab `xlab' `=`j'+1' " " + + csdid2_default , `options' + local gf11 `s(df11)' + local gf12 `s(df12)' + local style `s(style)' + local dels `s(delse)' + + if `"`xtitle'"'=="" & `ktype' ==3 local xtitle xtitle("Group") + else if `"`xtitle'"'=="" & `ktype' ==4 local xtitle xtitle("Calendar") + if `"`ytitle'"'=="" local ytitle ytitle("ATT") + + mata:st_local("adj",strofreal(csdid.antici)) + + tempvar t2 + qui:myencode `t', gen(`t2') `format' + if "`asy'"=="" { + two (`style' `ll' `uu' `t2' , `gf11' ) /// + (scatter `b' `t2' , `gf12' ) , /// + legend(off) `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `title' /// + xlab( ,val) `name' `dels' + } + else { + //gettoken se rest:rest + gettoken all rest:rest + gettoken auu rest:rest + + two (`style' `ll' `uu' `t2' , `gf11' lwidth(1)) /// + (`style' `all' `auu' `t2' , `gf11' ) /// + (scatter `b' `t2' , `gf12' ) , /// + legend(off) `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `title' /// + xlab( ,val) `name' `dels' + } + + +end + +*program drop myencode +program define myencode + syntax varname, gen(name) [format(string asis)] + if "`=`varlist''"=="GAverage" { + local torep "GAverage" + *replace `varlist' = `varlist'[2]-(`varlist'[3]-`varlist'[2]) + } + if "`=`varlist''"=="TAverage" { + local torep "TAverage" + *replace `varlist' = `varlist'[2]-(`varlist'[3]-`varlist'[2]) + } + + if "`torep'"!="" { + replace `varlist'=subinstr(`varlist',"g","",1) if _n>1 + replace `varlist'=subinstr(`varlist',"t","",1) if _n>1 + } + else { + replace `varlist'=subinstr(`varlist',"g","",1) + replace `varlist'=subinstr(`varlist',"t","",1) + } + qui:destring `varlist', force gen(`gen') + if "`torep'"!="" { + replace `gen' = `gen'[2]-(`gen'[3]-`gen'[2]) in 1 + tempname aux + label define _aux_ `=`gen'[1]' "`torep'", modify + label values `gen' _aux_ + } +qui format `format' `gen' +end + + diff --git a/csdid2 - Copy/csdid2_plotx.ado b/csdid2 - Copy/csdid2_plotx.ado new file mode 100644 index 0000000..b0f9ce3 --- /dev/null +++ b/csdid2 - Copy/csdid2_plotx.ado @@ -0,0 +1,281 @@ +*! v1.03 csdid2_plot Allows for Anticpiation +* v1.02 csdid2_plot Fixes `' +* v1.01 csdid2_plot for csdid2 only + +program csdid2_plot, rclass + + syntax, [* ktype(string)] + local cmd `r(cmd)' + local agg `r(agg)' + tempname b V table bb vv + matrix `b' = r(b) + matrix `V' = r(V) + matrix `table' = r(table) + + if "`agg'"=="group" local ktype = 3 + if "`agg'"=="calendar" local ktype = 4 + if "`agg'"=="event" local ktype = 5 + + + capture noisily csdid2_plot_wh, `options' ktype(`ktype') table(`table') + return local cmd `cmd' + return local agg `agg' + + return matrix b =`b' + return matrix V =`V' + return matrix table =`table' + +end + +program csdid2_plot_wh + syntax, [ * ktype(int 5) table(str) asy level(int 95)] + + tempname tbl + matrix `tbl'=`table' + capture: confirm matrix `tbl' + + *if det(`tbl')==. matrix `tbl'=rtb + if "`asy'"=="" { + if `ktype'==5 { + tempvar t b ll uu + + mata:event_p("`t' `b' `ll' `uu'","`tbl'") + + csdid_plot_eventx `t' `b' `ll' `uu', `options' + } + else if `ktype'==3 | `ktype'==4 { + // Group Calendar + tempvar t b ll uu + mata:other_p("`t' `b' `ll' `uu'","`tbl'") + + csdid_plot_other `t' `b' `ll' `uu', `options' ktype(`ktype') + } + else { + display in red "Plot option not allowed" + } + } + else { + if `ktype'==5 { + tempvar t b ll uu se all auu + + mata:event_p2("`t' `b' `ll' `uu' `se' `all' `auu'","`tbl'", `=`level'/100') + csdid_plot_eventx `t' `b' `ll' `uu' `all' `auu', `options' asy + } + else if `ktype'==3 | `ktype'==4 { + // Group Calendar + tempvar t b ll uu se all auu + mata:other_p2("`t' `b' `ll' `uu' `se' `all' `auu'","`tbl'", `=`level'/100') + `t' `b' `ll' `uu' `se' `all' `auu' + csdid_plot_other `t' `b' `ll' `uu' `all' `auu', `options' ktype(`ktype') asy + } + else { + display in red "Plot option not allowed" + } + } + +end + +mata: + void event_p2( string scalar newvars, string scalar tblx, real scalar level){ + real matrix tbl, ntbl2 + string matrix ntbl + real scalar alpha + tbl = st_matrix(tblx) + //asume always here + + alpha = 1-level + ntbl = st_matrixcolstripe(tblx) + ntbl = usubinstr(ntbl,"tp","+",.) + ntbl = usubinstr(ntbl,"tm","-",.) + ntbl2= strtoreal(ntbl) + tbl = tbl[(1,5,6,2),]' + tbl = tbl,tbl[,1]-tbl[,4]*invnormal(1-alpha/2),tbl[,1]+tbl[,4]*invnormal(1-alpha/2) + + tbl = select(tbl,(ntbl2[,2]:!=.)) + ntbl2= select(ntbl2[,2],(ntbl2[,2]:!=.)) + real matrix ss + ss= _st_addvar("double",tokens(newvars)) + st_store((1::rows(tbl)) ,tokens(newvars),(ntbl2,tbl)) + } + + void other_p2(string scalar newvars, string scalar tblx, real scalar level){ + real matrix tbl + string matrix ntbl + real scalar alpha + alpha = 1-level + tbl = st_matrix(tblx) + ntbl = st_matrixcolstripe(tblx) + //ntbl = usubinstr(ntbl,"g","",.) + //ntbl = usubinstr(ntbl,"t","",.) + ntbl = ntbl [,2] + tbl = tbl[(1,5,6,2),]' + tbl = tbl,tbl[,1]-tbl[,4]*invnormal(1-alpha/2),tbl[,1]+tbl[,4]*invnormal(1-alpha/2) + + string matrix tnv + tnv = tokens(newvars) + real matrix ss + ss= _st_addvar(sprintf("str%f",max(strlen(ntbl))),tnv[1]) + ss= _st_addvar("double",tnv[2..4]) + st_sstore((1::rows(tbl)) ,tnv[1],ntbl) + st_store((1::rows(tbl)) ,tnv[2..4],tbl) + } +end + +program csdid2_default, sclass + syntax, [style(str) PSTYle1(str) color1(str) /// + PSTYLE2(str) color2(str) /// + LWidth1(str) lwidth2(str) /// + BARWidth1(str) barwidth2(str) * asy] + + if "`style'"=="" local style rspike + + + if "`pstyle1'"=="" local pstyle1 pstyle(p1) + else local pstyle1 pstyle(`pstyle1') + if "`pstyle2'"=="" local pstyle2 pstyle(p2) + else local pstyle2 pstyle(`pstyle2') + + if `"`color1'"'=="" local color1 color(%40) + else local color1 color(`"`color1'"') + if `"`color2'"'=="" local color2 color(%40) + else local color2 color(`"`color2'"') + + if "`style'"=="rspike" { + if "`lwidth1'"=="" local lwidth1 lwidth(3) + if "`lwidth2'"=="" local lwidth2 lwidth(3) + if "`asy'"!="" { + local lwidth1 lwidth(1) + local lwidth2 lwidth(1) + } + } + + if "`style'"=="rarea" { + if "`lwidth1'"=="" local lwidth1 lwidth(0) + if "`lwidth2'"=="" local lwidth2 lwidth(0) + local conn connect(l) + } + + if "`style'"=="rcap" { + if "`lwidth1'"=="" local lwidth1 lwidth(1) + if "`lwidth2'"=="" local lwidth2 lwidth(1) + local conn connect(l) + } + + if "`style'"=="rbar" { + if "`lwidth1'"=="" local lwidth1 lwidth(0) + if "`lwidth2'"=="" local lwidth2 lwidth(0) + if "`barwidth1'"=="" local barwidth1 barwidth(0.5) + if "`barwidth2'"=="" local barwidth2 barwidth(0.5) + local conn connect(l) + } + + + sreturn local style `style' + sreturn local df11 `pstyle1' `color1' `lwidth1' `barwidth1' + sreturn local df12 `pstyle1' `conn' + sreturn local df21 `pstyle2' `color2' `lwidth2' `barwidth2' + sreturn local df22 `pstyle2' `conn' + sreturn local delse `options' +end + +program csdid_plot_eventx + syntax varlist, [ xtitle(passthru) ytitle(passthru) /// + legend(passthru) asy * ] + gettoken t rest:varlist + gettoken b rest:rest + gettoken ll rest:rest + gettoken uu rest:rest + + + ** defaults + + ** defaults + if `"`xtitle'"'=="" local xtitle xtitle("Periods to Treatment") + if `"`ytitle'"'=="" local ytitle ytitle("ATT") + + if `"`legend'"'=="" local legend legend(order(1 "Pre-treatment" 3 "Post-treatment")) + csdid2_default , `options' `asy' + + local gf11 `s(df11)' + local gf12 `s(df12)' + local gf21 `s(df21)' + local gf22 `s(df22)' + local style `s(style)' + local dels `s(delse)' + + mata:st_local("adj",strofreal(csdid.antici)) + if "`asy'"=="" { + two (`style' `ll' `uu' `t' if `t'<=(-1- `adj'), `gf11') /// + (scatter `b' `t' if `t'<=(-1- `adj'), `gf12') /// + (`style' `ll' `uu' `t' if `t'> (-1- `adj'), `gf21') /// + (scatter `b' `t' if `t'> (-1- `adj'), `gf22') , /// + `legend' `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `dels' + } + else { + //gettoken se rest:rest + gettoken all rest:rest + gettoken auu rest:rest + + two (`style' `ll' `uu' `t' if `t'<=(-1- `adj'), `gf11') /// + (scatter `b' `t' if `t'<=(-1- `adj'), `gf12') /// + (`style' `ll' `uu' `t' if `t'> (-1- `adj'), `gf21') /// + (scatter `b' `t' if `t'> (-1- `adj'), `gf22') /// + (`style' `all' `auu' `t' if `t'<=(-1- `adj'), `gf11' lwidth(3)) /// + (`style' `all' `auu' `t' if `t'> (-1- `adj'), `gf21' lwidth(3)), /// + `legend' `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `dels' + + } +end + + +program csdid_plot_other + syntax varlist, [ktype(int 3) * /// + xtitle(passthru) ytitle(passthru) /// + legend(passthru)] + gettoken t rest:varlist + gettoken b rest:rest + gettoken ll rest:rest + gettoken uu rest:rest + + + local xlab `xlab' `=`j'+1' " " + + csdid2_default , `options' + local gf11 `s(df11)' + local gf12 `s(df12)' + local style `s(style)' + local dels `s(delse)' + + if `"`xtitle'"'=="" & `ktype' ==3 local xtitle xtitle("Group") + else if `"`xtitle'"'=="" & `ktype' ==4 local xtitle xtitle("Calendar") + if `"`ytitle'"'=="" local ytitle ytitle("ATT") + + mata:st_local("adj",strofreal(csdid.antici)) + + tempvar t2 + + qui:encode `t', gen(`t2') + if "`asy'"=="" { + two (`style' `ll' `uu' `t2' , `gf11' ) /// + (scatter `b' `t2' , `gf12' ) , /// + legend(off) `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `title' /// + xlab( ,val) `name' `dels' + } + else { + //gettoken se rest:rest + gettoken all rest:rest + gettoken auu rest:rest + two (`style' `ll' `uu' `t2' , `gf11' ) /// + (`style' `all' `auu' `t2' , lwidth(3) `gf11' ) /// + (scatter `b' `t2' , `gf12' ) , /// + legend(off) `xtitle' `ytitle' /// + yline(0 , lp(dash) lcolor(black)) `title' /// + xlab( ,val) `name' `dels' + } + + +end + diff --git a/csdid2 - Copy/csdid2_setup.ado b/csdid2 - Copy/csdid2_setup.ado new file mode 100644 index 0000000..12a893e --- /dev/null +++ b/csdid2 - Copy/csdid2_setup.ado @@ -0,0 +1,16 @@ +*! v0 Simple Setup. Adds info so we do not need to type all the time +program csdid2_setup + syntax, [ivar(varname) gvar(varname) tvar(varname) clear] + if "`clear'"!="" { + char _dta[ivar] + char _dta[gvar] + char _dta[tvar] + char _dta[csdid2] + } + else { + char _dta[ivar] `ivar' + char _dta[gvar] `gvar' + char _dta[tvar] `tvar' + char _dta[csdid2] csdid2 + } +end \ No newline at end of file diff --git a/csdid2 - Copy/csdid_stats.mata b/csdid2 - Copy/csdid_stats.mata new file mode 100644 index 0000000..e4320b7 --- /dev/null +++ b/csdid2 - Copy/csdid_stats.mata @@ -0,0 +1,701 @@ +*! v1.2 Balance event +*! v1.1 Corrects Group SE +*! v1 Allows for anticipation +mata +class select_range { + real matrix selgvar + real matrix seltvar + real matrix selevent + real matrix selbal +} +// mata drop estat +// mata drop csdid_estat() + +class csdid_estat { + // functions to create tables + void bvcv_asym() + void bvcv_clus() + // Aggregatprs + // void attgt() + void atts_asym() + void atts_wboot() + void group_att() + void calendar_att() + void simple_att() + void cevent_att() + void event_att() + void pretrend() + // makes WB and the table for WB + void mboot_any() + void make_table() + // void init only for testing + void init() + string matrix attgt_names() + + real matrix aggte() + real matrix rtokens() + real matrix select_data() + real matrix wmult() + real matrix iqrse() + real scalar qtc() + + + real matrix erif, table + real matrix bb, vv, sderr, bsmean + string matrix onames + // to be initialized + real scalar t_stat + // Required Info + real scalar cilevel, bwtype, reps, max_mem, test_type, noavg + // info created + real scalar nclust, nobs, ggroups, ccalendar, eevent, error + // to transfer Range Info. + class select_range scalar range +} + +// Creates ASYM VCV +void csdid_estat::init() { + cilevel = 0.95 + bwtype = 1 + reps = 999 + max_mem = 1 + test_type = 1 +} +void csdid_estat::bvcv_asym(real matrix rif) { + bb = mean(rif) + nobs= rows(rif) + vv = quadcrossdev(rif,bb, rif,bb) :/ (nobs^2) +} +// Creates ASYM Cluster VCV +// Need to think how to Compress data when clustered + +void csdid_estat::bvcv_clus(real matrix rif, + real matrix cvar) { + real matrix ord, info + bb = mean(rif) + nobs= rows(rif) + // sort //ord = order(cvar,1) //rif = rif[ord,] //cvar= cvar[ord,] + // Standard Errors + info = panelsetup(cvar,1) + nclust= rows(info) + real matrix sumrif + sumrif= panelsum(rif:-bb,info) + vv = quadcross(sumrif,sumrif):/(nobs^2) + // unsort //rif = rif[invorder(ord),] //cvar= cvar[invorder(ord),] +} + + +real matrix csdid_estat::aggte(real matrix rif,| real matrix wgt ) { + + real matrix mn_all, mn_rif, mn_wgt + if (args()==1) { + wgt = J(1,cols(rif),1) + } + // Avg Effect + mn_rif = mean(rif) + mn_wgt = mean(wgt) + mn_all = sum(mn_rif:*mn_wgt):/sum(mn_wgt) + // gets agg rif + real matrix wgtw, attw + wgtw = (mn_wgt ) :/sum(mn_wgt) + attw = (mn_rif ) :/sum(mn_wgt) + // r1 r2 r3 + real matrix r1 , r2 , r3 + r1 = (wgtw:*(rif :-mn_rif)) + r2 = (attw:*(wgt :-mn_wgt )) + r3 = (wgt :- mn_wgt) :* (mn_all :/ sum(mn_wgt) ) + // Aggregates into 1 + return(rowsum(r1):+rowsum(r2):-rowsum(r3):+mn_all) +} + + + +// Will use Separate function for WB bc it process data differently +void csdid_estat::atts_asym(class csdid scalar csdid){ + // Estimate effects + error = 0 + if (test_type==1) { + // ATTGT + erif=select(csdid.frif,select_data(csdid)') + if (length(csdid.cvar)==0) { + bvcv_asym(erif) + } + else { + bvcv_clus(erif,csdid.cvar) + } + // names + onames=attgt_names(csdid)' + } + else if (test_type==2) { + //simple att + simple_att(csdid) + if (length(csdid.cvar)==0) bvcv_asym(erif) + else bvcv_clus(erif,csdid.cvar) + onames = J(rows(onames),1,""),onames + } + else if (test_type==3) { + //group att + group_att(csdid) + if (length(csdid.cvar)==0) bvcv_asym(erif) + else bvcv_clus(erif,csdid.cvar) + onames = J(rows(onames),1,""),onames + } + else if (test_type==4) { + //calendar att + calendar_att(csdid) + if (length(csdid.cvar)==0) bvcv_asym(erif) + else bvcv_clus(erif,csdid.cvar) + onames = J(rows(onames),1,""),onames + } + else if (test_type==5) { + //event att + event_att(csdid) + if (length(csdid.cvar)==0) bvcv_asym(erif) + else bvcv_clus(erif,csdid.cvar) + onames = J(rows(onames),1,""),onames + } + else if (test_type==6) { + //cevent att + cevent_att(csdid) + if (length(csdid.cvar)==0) bvcv_asym(erif) + else bvcv_clus(erif,csdid.cvar) + onames = J(rows(onames),1,""),onames + } + + if (error == 0) { + st_matrix("_bb",bb) + st_matrix("_vv",vv) + st_matrixcolstripe("_bb", onames) + st_matrixrowstripe("_vv", onames) + st_matrixcolstripe("_vv", onames) + } + else { + stata(`"display in red "There was an error estimating aggregation" "') + } + // drops vv, the largest matrix + vv=0 + erif = 0 +} + +void csdid_estat::make_table(){ + real matrix serr + real matrix ci + // Standard error + serr = iqrse(bsmean) + // Critical value + + t_stat =qtc( abs(bsmean:/serr) , cilevel ) + + // bb are point estimates + ci = (bb:-serr:*t_stat) \ (bb:+serr:*t_stat) + + table = bb \ serr \ (bb:/serr) \ ci \ J(1,cols(bb),t_stat) + bsmean=J(0,0,.) + +} + +void csdid_estat::atts_wboot(class csdid scalar csdid){ + + // Estimate effects + if (test_type==1) { + // ATTGT + error=0 + erif=select(csdid.frif,select_data(csdid)') + onames=attgt_names(csdid)' + mboot_any(csdid) + make_table() + // names + } + else if (test_type==2) { + //simple att + simple_att(csdid) + onames = J(rows(onames),1,""),onames + mboot_any(csdid) + make_table() + + } + else if (test_type==3) { + //group att + group_att(csdid) + onames = J(rows(onames),1,""),onames + mboot_any(csdid) + make_table() + } + else if (test_type==4) { + //calendar att + calendar_att(csdid) + onames = J(rows(onames),1,""),onames + mboot_any(csdid) + make_table() + } + else if (test_type==5) { + //event att + event_att(csdid) + onames = J(rows(onames),1,""),onames + mboot_any(csdid) + make_table() + } + else if (test_type==6) { + //cevent att + cevent_att(csdid) + onames = J(rows(onames),1,""),onames + mboot_any(csdid) + make_table() + } + + if (error == 0) { + string matrix xnames + xnames =J(6,1,""),("b"\"se"\"t"\"ll"\"ul"\"crit") + + st_matrix("_table",table) + st_matrixrowstripe("_table", xnames) + st_matrixcolstripe("_table", onames) + } + else { + stata(`"display in red "There was an error estimating aggregation" "') + } + // drops vv, the largest matrix + } + + +string csdid_estat::attgt_names(class csdid scalar csdid){ + real scalar i + string matrix toreturn + real matrix sfgtvar + + sfgtvar=select(csdid.fgtvar,select_data(csdid)) + toreturn=J(2,rows(sfgtvar),"") + for(i=1;i<=cols(toreturn);i++){ + toreturn[1,i]=sprintf("g%f",sfgtvar[i,1]:+csdid.antici) + toreturn[2,i]=sprintf("t%f_%f",sfgtvar[i,3],sfgtvar[i,4]) + } + return(toreturn) +} + +real matrix csdid_estat::rtokens(string scalar totok){ + return(uniqrows(strtoreal(tokens(totok))' )') +} +//////////////////////////////////////////////////////////////////////////////// +/// Group Aggregations +//////////////////////////////////////////////////////////////////////////////// + +void csdid_estat::group_att(class csdid scalar csdid ){ + // Counter i + // kgroups (max) + real scalar i , iic + real matrix toselect0,toselect, aux_rif, sumwgt + real matrix aux_wgt, aux + toselect0 =select_data(csdid)' + //:*csdid.convar' + error=0 + if (sum(toselect0)>0) { + ggroups = rows(csdid.sgvar) + nobs = rows(csdid.frif) + + onames=J(ggroups+1,1,"") + onames[1,]="GAverage" + + aux =J(nobs,ggroups,.) + sumwgt =J(nobs,ggroups,.) + iic=0 + for(i=1;i<=ggroups;i++){ + // select + toselect=toselect0:*(csdid.fgtvar[,1]:==csdid.sgvar[i]:& csdid.eventvar:>=0)' + if (sum(toselect)>0) { + iic++ + // if any selected -> Estimate + // Anticip only affects how things look + onames[iic+1,] = sprintf("g%f", csdid.sgvar[i]:+csdid.antici) + aux_wgt = select(csdid.frwt,toselect) + aux [,iic] = aggte(select(csdid.frif,toselect),aux_wgt ) + sumwgt[,iic] = rowsum(aux_wgt):/cols(aux_wgt) + //sumwgt[,iic] = aggte(aux_wgt) + } + } + // Drop Zeroes + sumwgt = sumwgt[,1..iic] + aux = aux[,1..iic] + onames = onames[1..iic+1,] + //sumwgt = colsum(sumwgt) + erif= aggte(aux,sumwgt ), aux + // If request no AVG + + if (noavg==1) { + erif = erif[.,2..cols(erif)] + onames =onames[2..rows(onames),] + } + + } + else { + error=1 + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Calendaar Aggregations +//////////////////////////////////////////////////////////////////////////////// + +void csdid_estat::calendar_att(class csdid scalar csdid ){ + // Counter i + // kgroups (max) + real scalar i , iic + real matrix toselect0,toselect, aux_rif, sumwgt + real matrix aux_wgt, aux + + toselect0 =select_data(csdid)' + //:*csdid.convar' + error=0 + if (sum(toselect0)>0) { + ccalendar = rows(csdid.stvar) + nobs = rows(csdid.frif) + + onames=J(ccalendar+1,1,"") + onames[1,]="TAverage" + + aux =J(nobs,ccalendar,.) + //sumwgt =J(nobs,ccalendar,.) + iic=0 + for(i=1;i<=ccalendar;i++){ + // select + toselect=toselect0:*(csdid.fgtvar[,2]:==csdid.stvar[i] :& csdid.eventvar:>=0)' + + if (sum(toselect)>0) { + iic++ + // if any selected -> Estimate + onames[iic+1,] = sprintf("t%f", csdid.stvar[i]) + aux_wgt = select(csdid.frwt,toselect) + aux [,iic] = aggte(select(csdid.frif,toselect),aux_wgt ) + //sumwgt[i,] = rowsum(aux_wgt):/cols(aux_wgt) + } + } + + //sumwgt = sumwgt[,1..iic] + aux = aux[,1..iic] + onames = onames[1..iic+1,] + erif = aggte(aux, J(1,cols(aux),1) ), aux + // If request no AVG + if (noavg==1) { + erif = erif[.,2..cols(erif)] + onames =onames[2..rows(onames),] + } + + } + else { + error=1 + } +} + +void csdid_estat::pretrend(class csdid scalar csdid ){ + // should be always drop v? + real scalar df + real matrix toselect,toselect0 + toselect0=select_data(csdid)' + //:*csdid.convar' + toselect=toselect0:*(csdid.eventvar :< 0)' + error=0 + if (sum(toselect)>0) { + if (length(csdid.cvar)==0) bvcv_asym(select(csdid.frif,toselect)) + else bvcv_clus(select(csdid.frif,toselect),csdid.cvar) + + real scalar chi2 + chi2=bb*invsym(vv)*bb' + df = cols(bb) + // Drops V matrix + vv=0 + st_numscalar("chi2_",chi2) + st_numscalar("df_",df) + st_numscalar("pchi2_",chi2tail(df,chi2)) + } + else { + error = 1 + } +} + +void csdid_estat::simple_att(class csdid scalar csdid ){ + real matrix toselect,toselect0 + // nobs = rows(csdid.frif) + // Select based on some criteria + toselect0 = select_data(csdid)' + //:*csdid.convar' + toselect = toselect0:*(csdid.eventvar:>=0)' + error = 0 + if (sum(toselect)>0) { + onames = "SimpleATT" + erif = aggte(select(csdid.frif,toselect),select(csdid.frwt,toselect) ) + } + else { + error = 1 + } +} + +real matrix csdid_estat::select_data(class csdid scalar csdid){ + real matrix toselect1,toselect2,toselect3,toselect4 + real scalar i, i1, i2, i3, i4 + real scalar rws + // Can we adapt this to other?. Yes we should be able to! + rws = rows(csdid.eventvar) + toselect1 =toselect2=toselect3=toselect4= J(rws,1,1) + i1=length(range.selgvar) + i2=length(range.seltvar) + i3=length(range.selevent) + i4=length(range.selbal) + + if (i1>0) { + toselect1 =J(rws,1,0) + for(i=1;i<=i1;i++){ + toselect1=toselect1:+(csdid.fgtvar[,1]:==range.selgvar[i]) + } + } + if (i2>0){ + toselect2 =J(rws,1,0) + for(i=1;i<=i2;i++){ + toselect2=toselect2:+(csdid.fgtvar[,2]:==range.seltvar[i]) + } + } + if (i3>0){ + toselect3 =J(rws,1,0) + for(i=1;i<=i3;i++){ + toselect3=toselect3:+(csdid.eventvar[,1]:==range.selevent[i]) + } + } + if (i4>0){ + real matrix gg, selbal2 + // Get smaller event: This is to verify we have all "reasonable" events in the list + gg = uniqrows(csdid.eventvar) + selbal2 = J(1,0,0) + for(i=1;i<=i4;i++){ + if ( sum(gg[,1]:==range.selbal[i]) == 1) { + selbal2=selbal2,range.selbal[i] + + } + } + i4 = length(selbal2) + toselect4 =J(rws,1,0) + for(i=1;i<=i4;i++){ + toselect4=toselect4:+(csdid.eventvar[,1]:==selbal2[i]) + } + + // select Ggroups + + gg = select(csdid.fgtvar[,1],toselect4:>0) + gg = uniqrows(gg,1) ; gg=select(gg[,1],gg[,2]:==i4) + + // do gvar again + toselect4 =J(rws,1,0) + for(i=1;i<=length(gg);i++){ + toselect4=toselect4:+(csdid.fgtvar[,1]:==gg[i]) + } + if (sum(toselect4)==0) { + _error(123,"No observations after imposing balance") + + } + } + + // tosel1 gvar + // tosel2 tvar + // tosel3 evar + // tosel4 evar + // if rbalance ---> tsel4 event Balance & tsel1 group balance + // \-> But if tsel3 exist It trumps tsel4 + + return( csdid.convar:* (toselect1:*toselect2:*toselect3:*toselect4):>0 ) +} + +void csdid_estat::cevent_att(class csdid scalar csdid ){ + real matrix toselect + toselect = select_data(csdid)' + error = 0 + if (sum(toselect)>0) { + onames = "ATTC" + erif = aggte(select(csdid.frif,toselect),select(csdid.frwt,toselect) ) + } + else { + error = 1 + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// event Aggregations +//////////////////////////////////////////////////////////////////////////////// + +void csdid_estat::event_att(class csdid scalar csdid){ + // Counter i + // kgroups (max) + real scalar i , iic , ievent + real matrix toselect,toselect0, aux_rif, sumwgt + real matrix aux_wgt, iim + real matrix aux_event, aux + + // + toselect0=select_data(csdid)' + error = 0 + if (sum(toselect0)>0) { + ievent = 0 + aux_event=select(csdid.eventvar,toselect0') + + eevent = rows(csdid.sevent) + nobs = rows(csdid.frif) + onames=J(eevent+sum(csdid.sevent:>=0)+sum(csdid.sevent:< 0),1,"") + + // Is there a Pre or post + iic = 0 + if (sum(aux_event:<0 )) { + iic++ + ievent++ + onames[iic,]="Pre_avg" + } + if (sum(aux_event:>=0)) { + iic++ + ievent++ + onames[iic,]="Post_avg" + } + aux =J(nobs,eevent,.) + + + iim =J(1,0,.) + for(i=1;i<=eevent;i++){ + // select + toselect=toselect0:*(csdid.eventvar:==csdid.sevent[i])' + + if (sum(toselect)>0) { + iic++ + // if any selected -> Estimate + // For Two. One is for names. The other one for iim + + if ((csdid.sevent[i])<0) { + iim = iim , 0 + } + else { + iim = iim , 1 + + } + + if ((csdid.sevent[i]:-csdid.antici)<0) { + onames[iic,] = sprintf("tm%f", abs(csdid.sevent[i]:-csdid.antici)) + } + else { + onames[iic,] = sprintf("tp%f", abs(csdid.sevent[i]:-csdid.antici)) + } + + aux_wgt = select(csdid.frwt,toselect) + aux[,iic-ievent] = aggte(select(csdid.frif,toselect),aux_wgt ) + //sumwgt[i,] = rowsum(aux_wgt):/cols(aux_wgt) + + } + } + // drop zeroes + + aux = aux[,1..iic-ievent] + onames = onames[1..iic,] + // iim ids pre and post effects + erif =J(nobs,0,.) + if (sum(iim:==0)) { + erif =erif, aggte(select(aux,iim:==0) ) + } + if (sum(iim:==1)) { + erif =erif, aggte(select(aux,iim:==1) ) + } + real scalar col_erif + col_erif=cols(erif)+1 + // erif + erif = erif, aux + + // If request no AVG + if (noavg==1) { + erif = erif[.,col_erif..cols(erif)] + onames =onames[col_erif..rows(onames),] + } + } + else { + error=1 + } +} + +/// Auxiliary programs +// Gets the pth value of the rowmas matrix sent +// used to get t-critical for uniform matrix +real scalar csdid_estat::qtc(real matrix y, real scalar p){ + // idea. maximizar + y =rowmax(y) + _sort(y,1) + if ( p>0 & p<1) return(y[ ceil( (rows(y)+1)*p ) ]) + else if (p==0) return(y[ 1 ]) + else if (p==1) return(y[ rows(y) ]) +} + +real matrix csdid_estat::iqrse(real matrix y) { + real scalar q25,q75 + // saves q25 and q75 + q25=ceil(rows(y)*.25);q75=ceil(rows(y)*.75) + real scalar j + real matrix iqrs, sy + iqrs=J(1,cols(y),0) + + for(j=1;j<=cols(y);j++){ + sy=sort(y[,j],1) + iqrs[,j]=(sy[q75,]-sy[q25,]):/(invnormal(.75)-invnormal(.25) ) + } + return(iqrs) +} + +real matrix csdid_estat::wmult(real scalar mdsize_eff) { + real scalar k1, k2 + k1=((1+sqrt(5))/(2*sqrt(5))) + k2=0.5*(1+sqrt(5)) + + if (bwtype==1) return( k2:-sqrt(5)*rbinomial(nobs,mdsize_eff,1, k1) ) + else if (bwtype==2) return( 1 :-2* rbinomial(nobs,mdsize_eff,1,0.5) ) + +} + + +void csdid_estat::mboot_any(class csdid scalar csdid ) { + // RIF is FED from out. + real matrix mean_rif + real scalar i, ncols, xnobs + real scalar coord1 + real matrix ccrd + //real matrix bsmean + // First Means + bb=mean_rif=mean(erif) + // Re-estimate RIF + erif = erif:-mean_rif + // contains all iterations, for reps parameters + ncols=cols(erif) + xnobs=rows(erif) + // + + bsmean=J(reps,ncols,0) + // Options for cluster. + real matrix info + if (length(csdid.cvar)>0) { + info=panelsetup(csdid.cvar,1) + erif= panelsum(erif,info) + } + + nobs =rows(erif) + + // check Repetitions and parameters + // This is to use BLOCKS of Stuff. But not sure about how large it can be + real scalar mdsize_eff, mdsize, mmax_mem + // 134217728 <-- Total number of observations in 1gb of memory. We can select More + // Need to initialize max men<- Global. MMaxmem local + if (max_mem==0) mmax_mem=134217728 + else mmax_mem=max_mem*134217728 + + mdsize = min( (reps, max( ( 1 , floor(mmax_mem/nobs/ncols) ) ) ) ) + + coord1=1 + mdsize_eff = mdsize + + for(i=1;i<=reps;i=i+mdsize){ + ccrd = (coord1,1) \ ( coord1+mdsize_eff-1 ,ncols) + coord1 = coord1+mdsize_eff + bsmean[|ccrd|]= cross(erif, wmult(mdsize_eff))':/xnobs + mdsize_eff = min( (mdsize, reps-(coord1-1)) ) + } + erif = J(0,0,.) + //return(bsmean) +} + +end diff --git a/csdid2 - Copy/ddd_csdid.ipynb b/csdid2 - Copy/ddd_csdid.ipynb new file mode 100644 index 0000000..d5cb937 --- /dev/null +++ b/csdid2 - Copy/ddd_csdid.ipynb @@ -0,0 +1,448 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Untitled\n", + "\n", + "## How to run DDD with CSDID in Stata\n", + "\n", + "Here I will show a small example, using the mpdta datafile.\n", + "\n", + "First the setup" + ], + "id": "e5b05ba8-b265-4716-af14-0d741ec20b61" + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "(Written by R. )" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "File f1.dta will be replaced" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + ".." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + ".\n", + "file f1.dta saved" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "\n", + "Difference-in-difference with Multiple Time Periods\n", + "\n", + " Number of obs = 1,295\n", + "Outcome model : regression adjustment\n", + "Treatment model: none" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "------------------------------------------------------------------------------" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + " | Coefficient Std. err. z P>|z| [95% conf. interval]\n", + "-------------+----------------------------------------------------------------\n", + "g2004 |" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + " t_2003_2004 | -.0211008 .041165 -0.51 0.608 -.1017827 .0595812\n", + " t_2003_2005 | -.1292921 .0454923 -2.84 0.004 -.2184553 -.0401288\n", + " t_2003_2006 | -.2074393 .0571321 -3.63 0.000 -.319416 -.0954625\n", + " t_2003_2007 | -.1355455 .0595128 -2.28 0.023 -.2521884 -.0189026\n", + "-------------+----------------------------------------------------------------\n", + "g2006 |\n", + " t_2003_2004 | .0054858 .0350106 0.16 0.875 -.0631337 .0741054\n", + " t_2004_2005 | -.0217113 .0381993 -0.57 0.570 -.0965805 .0531579\n", + " t_2005_2006 | -.024916 .0290084 -0.86 0.390 -.0817715 .0319395\n", + " t_2005_2007 | -.0544724 .0379769 -1.43 0.151 -.1289058 .0199611\n", + "-------------+----------------------------------------------------------------\n", + "g2007 |\n", + " t_2003_2004 | .0273227 .0273815 1.00 0.318 -.0263441 .0809894\n", + " t_2004_2005 | -.0117455 .0323164 -0.36 0.716 -.0750844 .0515935\n", + " t_2005_2006 | -.0344658 .031155 -1.11 0.269 -.0955285 .0265969\n", + " t_2006_2007 | -.0117583 .0322454 -0.36 0.715 -.0749581 .0514416\n", + "------------------------------------------------------------------------------" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "Control: Never Treated\n", + "\n", + "See Callaway and Sant'Anna (2021) for details\n", + "File f2.dta will be replaced" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + ".." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "." + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + ".\n", + "file f2.dta saved\n", + "\n", + "Difference-in-difference with Multiple Time Periods\n", + "\n", + " Number of obs = 1,205\n", + "Outcome model : regression adjustment\n", + "Treatment model: none\n", + "------------------------------------------------------------------------------\n", + " | Coefficient Std. err. z P>|z| [95% conf. interval]\n", + "-------------+----------------------------------------------------------------\n", + "g2004 |\n", + " t_2003_2004 | -.0036454 .017134 -0.21 0.832 -.0372275 .0299367\n", + " t_2003_2005 | -.016655 .0232565 -0.72 0.474 -.0622369 .0289269\n", + " t_2003_2006 | -.0705129 .0266075 -2.65 0.008 -.1226626 -.0183631\n", + " t_2003_2007 | -.0725062 .0217586 -3.33 0.001 -.1151523 -.0298601\n", + "-------------+----------------------------------------------------------------\n", + "g2006 |\n", + " t_2003_2004 | -.000845 .0293484 -0.03 0.977 -.0583669 .0566768\n", + " t_2004_2005 | .0069579 .0178941 0.39 0.697 -.0281139 .0420296\n", + " t_2005_2006 | .0125415 .024053 0.52 0.602 -.0346015 .0596845\n", + " t_2005_2007 | -.0352532 .0212842 -1.66 0.098 -.0769695 .006463\n", + "-------------+----------------------------------------------------------------\n", + "g2007 |\n", + " t_2003_2004 | .0270388 .01238 2.18 0.029 .0027745 .0513031\n", + " t_2004_2005 | .0024533 .0111352 0.22 0.826 -.0193712 .0242778\n", + " t_2005_2006 | -.025611 .0216096 -1.19 0.236 -.067965 .016743\n", + " t_2006_2007 | -.042635 .0129499 -3.29 0.001 -.0680163 -.0172538\n", + "------------------------------------------------------------------------------\n", + "Control: Never Treated\n", + "\n", + "See Callaway and Sant'Anna (2021) for details" + ] + } + ], + "source": [ + "qui:ssc install frause, replace\n", + "frause mpdta, clear\n", + "** Assume two groups, based on lpop\n", + "gen g1=lpop>3.3\n", + "** Now use csdid to estimate the model, without covariates\n", + "csdid lemp if g1==0, ivar(countyreal) time(year) gvar(first_treat) saverif(f1) replace\n", + "csdid lemp if g1==1, ivar(countyreal) time(year) gvar(first_treat) saverif(f2) replace" + ], + "id": "07b7ca43-5e9f-4a04-964a-ad6eb5f95668" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now there are two files with the corresponding RIFs, so I can get\n", + "aggregations for each separately" + ], + "id": "c5941cea-0f8a-4991-ace5-72af19358579" + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "(Written by R. )\n", + "\n", + "------------------------------------------------------------------------------\n", + " | Coefficient Std. err. z P>|z| [95% conf. interval]\n", + "-------------+----------------------------------------------------------------\n", + " ATT | -.0526544 .0223489 -2.36 0.018 -.0964575 -.0088513\n", + "------------------------------------------------------------------------------\n", + "file f1.dta saved\n", + "(Written by R. )\n" + ] + }, + { + "output_type": "stream", + "name": "stdout", + "text": [ + "------------------------------------------------------------------------------\n", + " | Coefficient Std. err. z P>|z| [95% conf. interval]\n", + "-------------+----------------------------------------------------------------\n", + " ATT | -.0328 .0094731 -3.46 0.001 -.051367 -.014233\n", + "------------------------------------------------------------------------------\n", + "file f2.dta saved" + ] + } + ], + "source": [ + "use f1, clear\n", + "csdid_stats simple, save\n", + "ren ATT ATT_g1\n", + "save f1, replace\n", + "\n", + "use f2, clear\n", + "csdid_stats simple, save\n", + "ren ATT ATT_g2\n", + "save f2, replace" + ], + "id": "836c5aa7-126e-4e3b-becf-b7a3ad1ec877" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next I need to merge them together" + ], + "id": "64ab3147-657a-467d-84c9-415c546ea194" + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "(Written by R. )\n", + "\n", + " Result Number of obs\n", + " -----------------------------------------\n", + " Not matched 500\n", + " from master 259 (_merge==1)\n", + " from using 241 (_merge==2)\n", + "\n", + " Matched 0 (_merge==3)\n", + " -----------------------------------------\n", + "\n", + " Variable | Obs Mean Std. dev. Min Max\n", + "-------------+---------------------------------------------------------\n", + " ATT_g1 | 259 -.0526544 .3603682 -3.026657 2.10358\n", + " ATT_g2 | 241 -.0328 .1473684 -.5103483 .7988653" + ] + } + ], + "source": [ + "use f1, clear\n", + "merge 1:1 countyreal using f2\n", + "sum ATT_g1 ATT_g2" + ], + "id": "ad7af8b5-d765-4637-84eb-e7341171a84d" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now you have both files together, with different treatment effects. To\n", + "create a table you can use the command csdid_rif, to compare both\n", + "results, and make test for differences" + ], + "id": "0ed16375-b221-4b4b-a49a-a26e37f0ca69" + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "------------------------------------------------------------------------------\n", + " | Robust\n", + " | Coefficient std. err. z P>|z| [95% conf. interval]\n", + "-------------+----------------------------------------------------------------\n", + " ATT_g1 | -.0526544 .0223489 -2.36 0.018 -.0964575 -.0088513\n", + " ATT_g2 | -.0328 .0094731 -3.46 0.001 -.051367 -.014233\n", + "------------------------------------------------------------------------------\n", + "\n", + " ( 1) ATT_g1 - ATT_g2 = 0\n", + "\n", + " chi2( 1) = 0.67\n", + " Prob > chi2 = 0.4134" + ] + } + ], + "source": [ + "csdid_rif ATT_g1 ATT_g2\n", + "test ATT_g1= ATT_g2" + ], + "id": "dc674da3-6687-4d94-adeb-4b8a11f719f2" + } + ], + "nbformat": 4, + "nbformat_minor": 5, + "metadata": { + "kernelspec": { + "name": "nbstata", + "display_name": "Stata (nbstata)", + "language": "stata" + }, + "language_info": { + "name": "stata", + "file_extension": ".do", + "mimetype": "text/x-stata", + "version": "17" + } + } +} \ No newline at end of file diff --git a/csdid2 - Copy/drdid.mata b/csdid2 - Copy/drdid.mata new file mode 100644 index 0000000..a3e2750 --- /dev/null +++ b/csdid2 - Copy/drdid.mata @@ -0,0 +1,1127 @@ +*! v1.1 Two new Methods RC2 and DRIMP2 +*! Also Addressed Bug with Clustered ST errors with RC +* v1 Allows for anticipation + +mata +// IPT +void fipt(transmorphic M, + real scalar todo, + real rowvector b, + real scalar lnf, + real matrix gg, + real matrix hh){ + + real colvector y1, xb + y1 = moptimize_util_depvar(M,1) + + xb = moptimize_util_xb(M,b,1) + + lnf = moptimize_util_sum(M, y1:*xb :- (y1:==0):*exp(xb) ) + if (todo>=1){ + gg = moptimize_util_vecsum(M, 1, y1:-(y1:==0):*exp(xb) , lnf) + if (todo==2){ + hh = moptimize_util_matsum(M, 1,1, -(y1:==0):*exp(xb) , lnf) + } + } +} +// logit +void flogit(transmorphic M, + real scalar todo, + real rowvector b, + real scalar lnf, + real matrix gg, + real matrix hh){ + real colvector y1, pxb, xb, exb + + y1 = moptimize_util_depvar(M,1) + xb = moptimize_util_xb(M,b,1) + pxb = logistic(xb) + exb = exp(xb) + + lnf = moptimize_util_sum(M, y1:*xb :- ln(exb:+1)) + if (todo>=1){ + gg = moptimize_util_vecsum(M, 1, y1:-exb:/(exb:+1) , lnf) + if (todo==2){ + hh = moptimize_util_matsum(M, 1,1, -exb:/(exb:+1):^2 , lnf) + } + } +} + + +class drdid { + // For later. + real scalar data_type + real scalar method_type + + // Data needed + real matrix yvar + + real matrix xvar + real matrix wvar + real matrix id + real matrix oid + real matrix trt + //real matrix wtrt + real matrix tmt + real matrix tmt_trt + + // also minn. if minn=1 -> reg regardless + //real scalar minn + + // data created + real matrix xb + real matrix yhat + real matrix b + //real matrix ixx + real matrix psv + real scalar nn + real scalar conv + real scalar kx + real scalar minn + real scalar rolljw + real scalar err + real matrix select4() + // regressions + ///void fipt() + void ipt() + ///void flogit() + void ilogit() + void ols() + void ols_ipw_rc() + void ols_ipt_rc() + // out + real matrix rif + + // Other functions + void reg_panel() + void drimp_panel() + void stdipw_panel() + void dripw_panel() + + void reg_rc() + // Reg2 is the Double reg + void reg2_rc() + void drimp_rc() + // This is Double Reg based on Drimp/ + void drimp2_rc() + + void stdipw_rc() + void dripw_rc() + + void makeid() + //void makeid2() + // setting up data + void msetup_panel() + void msetup_panel2() + void msetup_rc() + // this one "fixes stuff" + // void setup() + // void setup2() + void csdid_setup() + // master + void drdid() + // To "export, results" + // void drdid_outpt() + void init() +} + +//# Regression models + +void drdid::init(){ + yvar=xvar=wvar=id=oid=trt=tmt=tmt_trt=J(0,0,.) + xb=yhat=b= psv=J(0,0,.) + nn=conv=kx=minn=. +} + +void drdid::ipt(){ + ///real matrix sy,sx,sw + transmorphic M + M = moptimize_init() + moptimize_init_evaluator(M, &fipt()) + moptimize_init_depvar(M,1, trt) + + moptimize_init_weight(M, wvar) + moptimize_init_eq_indepvars(M,1, xvar) + + b=J(1,cols(xvar),0),logit(mean(trt)) + moptimize_init_eq_coefs(M, 1, b) + moptimize_init_evaluatortype(M, "d2") + moptimize_init_conv_maxiter(M, 100) + moptimize_init_tracelevel(M, "none") + moptimize_init_conv_warning(M, "off") + moptimize_init_verbose(M, "off") + moptimize(M) + b= moptimize_result_coefs(M) + conv=moptimize_result_converged(M) + + xb =(xvar,J(nn,1,1))*b' + +} + +void drdid::ilogit(){ + transmorphic M + M = moptimize_init() + + moptimize_init_evaluator(M, &flogit()) + moptimize_init_depvar(M,1, trt) + moptimize_init_weight(M, wvar) + moptimize_init_eq_indepvars(M,1, xvar) + b=J(1,cols(xvar),0),logit(mean(trt,wvar)) + moptimize_init_eq_coefs(M, 1, b) + moptimize_init_evaluatortype(M, "d2") + moptimize_init_conv_maxiter(M, 50) + moptimize_init_tracelevel(M, "none") + moptimize_init_conv_warning(M, "off") + moptimize_init_verbose(M, "off") + moptimize(M) + b =moptimize_result_coefs(M) + psv =moptimize_result_V(M) + conv=moptimize_result_converged(M) + xb =(xvar,J(nn,1,1))*b' + +} + + +void drdid::ols(real matrix sw, ixx ){ + real matrix xy + if (kx>0) { + ixx = invsym(quadcross(xvar,1,sw,xvar,1)) + xy = quadcross(xvar,1,sw,yvar,0) + b = ixx*xy + } + else { + b=mean(yvar,sw) + ixx=1/sum(sw:!=0) + } +} + + +/// Setting Data UP + +void drdid::msetup_panel(){ + // assume data is sorted + // makeid to ID balance panel + // expands current ID + err = 0 + makeid() + // keeps only those with 2 observations + tmt = tmt :==max(tmt) ; trt = trt :==max(trt) + + yvar = select(yvar,(id[,2]:==2)) + kx = cols(xvar) + if (kx>0) { + + // Keep Data from T0 (earlier) + xvar = select(xvar,((id[,2]:==2):*(tmt:==0))) + + // Keep only Xs with variation + xvar = select(xvar,diagonal(variance(xvar))':!=0) + } + + kx = cols(xvar) + + // if Method is not OLS then check. + // Basically If using other methods and there are no covariates, then + // OLS1 is the best option + if (method_type != 4) { + if (kx ==0) { + method_type=4 + } + else if (kx >=minn) { + // if More covariates than obs Keep minn-1 variables + // Very Inneficient way to adress the problem. But may still be better than the alternative + method_type=4 + stata(`"display in red "More X's than Observations. Dropping Variables""') + xvar = xvar[,1..minn-1] + + } + } + + //wtrt = + wvar = select(wvar,((id[,2]:==2):*(tmt:==0))) + wvar = wvar:/mean(wvar) + + // Original Copy of selected cases + // I can probably get rid of ID and keep only OID for matching + oid = select(oid ,(id[,2]:==2):*(tmt:==0)) + + trt = select(trt,(id[,2]:==2):*(tmt:==0)) + //wtrt = wtrt:*trt + tmt = select(tmt,(id[,2]:==2)) + + id = select(id ,(id[,2]:==2)) + id = select(id ,(tmt:==0)) + + + yvar = select(yvar,(tmt:==1)):-select(yvar,(tmt:==0)) + nn = rows(yvar) + + if (rows(yvar)==0) err=1 +} + +// For Rolling Regressions + +void drdid::msetup_panel2(){ + // GOAL. get the Diff. Getting Mean + // + // assume data is sorted + // makeid to ID balance panel + // expands current ID + err = 0 + + makeid() + // keeps only those with 2 observations + + tmt = tmt :==max(tmt) ; trt = trt :==max(trt) + // Read Data in two blocks. + // One for data at T= + // Data for G-1 or earlier + + real matrix yvar_post, yvar_pre + yvar_post = select(yvar,(tmt:==1)) + yvar_pre = select(yvar,(tmt:==0)) + real matrix id2 + id2 = select(id,(tmt:==0)) + real matrix info + info = panelsetup(id2,1) + + yvar_pre = panelsum(yvar_pre,info):/(info[,2]:-info[,1]:+1) + + kx = cols(xvar) + // DOES NOT HANDLE time varying data + if (kx>0) { + // Keep Data from T0 (earlier): average Pre treatment + xvar = select(xvar,(tmt:==0)) + xvar = panelsum(xvar,info):/(info[,2]:-info[,1]:+1) + // Keep only Xs with variation + xvar = select(xvar,diagonal(variance(xvar))':!=0) + } + + kx = cols(xvar) + + // if Method is not OLS then check + + if (method_type != 4) { + if (kx ==0) { + method_type=4 + + } + else if (kx >=minn) { + // if More covariates than obs Keep minn-1 variables + + method_type=4 + stata(`"display in red "More X's than Observations. Dropping Variables""') + xvar = xvar[,1..minn-1] + + } + } + + //wtrt = + wvar = select(wvar,(tmt:==1)) + wvar = wvar:/mean(wvar) + + // Original Copy of selected cases + + oid = select(oid ,(tmt:==1)) + trt = select(trt,(tmt:==1)) + //wtrt = wtrt:*trt + id = select(id ,(tmt:==1)) + + yvar = yvar_post :-yvar_pre + nn = rows(yvar) + + if (rows(yvar)==0) err=1 + +} + +// Data setup +void drdid::msetup_rc(){ + // assume data is sorted + // makeid to ID balance panel + // expands current ID + // makeid2() + err=0 + + //Doesnt need to be Redone. It comes from the data + //oid + id=1::rows(id) + // keeps only those with 2 observations + tmt = tmt :==max(tmt) + trt = trt :==max(trt) + // How many variables in the data + kx = cols(xvar) + + tmt_trt = (tmt:+ 2*trt) + + if (kx>0) { + // If any variable. Check which variables have variation. + // Does Gets rid of Constants. + xvar = select(xvar, select4() ) + kx = cols(xvar) + } + // If no variables left, using SIMPLE OLS + if (kx ==0) { + method_type=4 + } + + kx = cols(xvar) + + //wtrt = wvar:*trt + // Standardizes Weights to avoid overshooting + wvar = wvar:/mean(wvar) + // Sample size + nn = rows(yvar) + //if min( uniqrows(tmt_trt,1))[,2] ==0 ) err =1 + +} + +// Selects covariates +real matrix drdid::select4(){ + real matrix s4way + //tmt_trt = (tmt:+ 2*trt) + s4way=(diagonal(variance(xvar,tmt_trt :==0))':!=0):* + (diagonal(variance(xvar,tmt_trt :==1))':!=0):* + (diagonal(variance(xvar,tmt_trt :==2))':!=0):* + (diagonal(variance(xvar,tmt_trt :==3))':!=0) + + return(s4way) +} + +// makes IDS for panel +// This makes it easy to transfer things from one side to the other. + +void drdid::makeid(){ + real scalar i,j + real matrix id2 + //makes a copy + //oid=id + id2=J(rows(id),1,0) + // Recode ID. Assumes ID are ordered + // + j = 1 + for(i=1;i<=rows(id);i++){ + if (i>1) { + if (id[i]:!=id[i-1]) { + j++ + } + } + id2[i]=j + } + // recoded + id=uniqrows(id2,1)[id2,] + +} + +/*void drdid::makeid2(){ + oid=id + id=1::rows(id) +}*/ + +//// REG Panel + +void drdid::reg_panel(){ + + real matrix wols, w_1 + real scalar mw_1 + real matrix ixx + + wols = wvar :* (1 :- trt) + w_1 = wvar :* trt + mw_1 = mean(w_1) + + // OLS Simple. no checks. + ols(wols,ixx) + + if (kx>0) xb = xvar*b[1..kx]:+b[kx+1] + else xb = b + + // adds constant + if (kx>0) xvar = xvar, J(nn,1,1) + else xvar = J(nn,1,1) + + /// ATTs + real matrix att_treat, att_cont + att_treat = w_1:* yvar + att_cont = w_1:* xb + + real scalar eta_treat, eta_cont + eta_treat = mean(att_treat):/mw_1 + eta_cont = mean(att_cont) :/mw_1 + + //real matrix XpX_inv + //XpX_inv = invsym(quadcross(xvar,wols,xvar))*nn + //wols_eX = wols :* (dy:-xb) :* xvar + real matrix lin_ols + lin_ols = ( wols :* (yvar:-xb) :* xvar ) * + ( ixx * nn ) + + //real matrix inf_treat, inf_cont_1, inf_cont_2, inf_control + + rif = (eta_treat :- eta_cont):+ + (att_treat :- w_1 * eta_treat)/mw_1 :- + (att_cont :- w_1 * eta_cont)/mw_1 :- + lin_ols * (mean(xvar, w_1))' + +} + + +void drdid::drimp_panel(){ + // estimate psxb + ipt() + if (conv==1) { + //xb=quadcross(xvar',b) + real matrix psc, ixx + psc=logistic(xb) + + real matrix w_1 , w_0, att + w_1 = wvar :* trt + w_0 = wvar :* psc :* (1:-trt):/(1:-psc) + w_1 = w_1:/mean(w_1) + w_0 = w_0:/mean(w_0) + + ols(w_0,ixx) + + if (kx>0) xb = xvar*b[1..kx]:+b[kx+1] + else xb = b + + att=(yvar:-xb):*(w_1:-w_0) + + rif = mean(att) :+ att :- w_1:*mean(att) + } +} + +void drdid::stdipw_panel() { + ilogit() + if (conv==1) { + real matrix psc, inf_cont_1 + psc=logistic(xb) + // and matrices + //psb =st_matrix(psb_ ) + + real matrix w_1, w_0, att_cont, att_treat, + eta_treat, eta_cont, + lin_ps + + w_1= wvar :* trt + w_0= wvar :* psc :* (1 :- trt):/(1 :- psc) + + real scalar mw_1 , mw_0 + mw_1=mean(w_1) + mw_0=mean(w_0) + + att_treat = w_1:* yvar + att_cont = w_0:* yvar + + eta_treat = mean(att_treat)/mw_1 + eta_cont = mean(att_cont) /mw_0 + //ipw_att = eta_treat :- eta_cont + //inf_treat = (att_treat :- (w_1 :* eta_treat))/mw_1 + inf_cont_1 = (att_cont :- (w_0 :* eta_cont )) + + xvar=xvar,J(nn,1,1) + //lin_ps = (wvar:* (trt :- psc) :* xvar)*(psv * nn) + //M2 = + ///inf_cont_2 = ( (wvar:* (trt :- psc) :* xvar) ) * ( psv * mean(inf_cont_1 :* xvar)' ) * nn + // inf_control = (inf_cont_1 :+ inf_cont_2)/mw_0 + + rif = ( eta_treat :- eta_cont ) :+ + ( (att_treat :- (w_1 :* eta_treat)) /mw_1 ) :- + ( inf_cont_1 :+ (wvar:* (trt :- psc) :* xvar ) * + ( psv * mean(inf_cont_1 :* xvar)') * nn )/mw_0 + } + } + + +void drdid::dripw_panel() { + + real matrix exb, psc + real matrix w_1, w_0, wols + real matrix dy_xb, ixx + real matrix lin_ols, lin_ps + + ilogit() + if (conv==1) { + real matrix nest + exb=exp(xb); psc=logistic(xb) + + w_1 = (wvar:*trt) + w_0 = (wvar:*psc:*(1:-trt):/(1:-psc)) + w_1 = w_1 / mean(w_1) + w_0 = w_0 / mean(w_0) + + // ols part + wols = wvar :* (1 :- trt) + + ols(wols,ixx) + + if (kx>0) xb = xvar*b[1..kx]:+b[kx+1] + else xb = b + + // adds a constant for the rest of the analysis + + xvar=xvar,J(nn,1,1) + dy_xb = yvar:-xb + //att = mean((w_1:-w_0):*(dy_xb)) + + lin_ols = (wols :* (dy_xb) :* xvar) * (ixx * nn) + lin_ps = (wvar :* (trt:-psc) :* xvar) * (psv * nn) + + // Components for RIF + real matrix aa + + //n1 = w_1:*(dy_xb:-mean(dy_xb,w_1)) + //n0 = w_0:*(dy_xb:-mean(dy_xb,w_0)) + aa = ((1:-trt):/(1:-psc):^2)/ mean(psc:*(1:-trt):/(1:-psc)) + + nest = lin_ols * (mean(xvar,w_1) :- mean(xvar ,w_0))' :+ + lin_ps * mean( aa :* (dy_xb :- mean(dy_xb,w_0)) :* + exb:/(1:+exb):^2:*xvar )' + + // RIF att_inf_func = inf_treat' :- inf_control + rif = mean((w_1:-w_0):*(dy_xb)):+ + w_1:*(dy_xb:-mean(dy_xb,w_1)):- + w_0:*(dy_xb:-mean(dy_xb,w_0)):- + nest + } + } + +void drdid::ols_ipw_rc(real scalar ii, real matrix yy, real matrix ixx ){ + ///tmt_trt = (tmt:+ 2*trt) + /// 0 0 = 0 + /// 1 0 = 1 + /// 0 1 = 2 + /// 1 1 = 3 + real matrix xy, sw + sw = wvar:*(tmt_trt:==ii) + sw = sw :/mean(sw) + if (kx>0) { + ixx = invsym(quadcross(xvar,1, sw ,xvar,1)) + xy = quadcross(xvar,1, sw ,yvar,0) + b = ixx*xy + yy = (xvar,J(nn,1,1))*b + } + else { + b=mean(yvar,sw) + ixx=1/sum(sw:!=0) + yy = b + } +} + +void drdid::ols_ipt_rc( real scalar ii, + real matrix ww, + real matrix yy, + real matrix ixx ){ + ///tmt_trt = (tmt:+ 2*trt) + /// 0 0 = 0 + /// 1 0 = 1 + /// 0 1 = 2 + /// 1 1 = 3 + + real matrix xy, sw + + sw = ww:*(tmt_trt:==ii) + + if (kx>0) { + + ixx = invsym(quadcross(xvar,1,sw,xvar,1)) + xy = quadcross(xvar,1,sw,yvar,0) + b = ixx*xy + yy = (xvar,J(nn,1,1))*b + } + else { + + b=mean(yvar,sw) + ixx=1/sum(sw:!=0) + yy = b + } +} + +void drdid::dripw_rc(){ + // main Loading variables + ilogit() + if (conv==1) { + real matrix psc + psc=logistic(xb) + /// tmt trt + real matrix y00, y01, y10, y11 + real matrix ixx00, ixx01, ixx10, ixx11 + real matrix w00 , w01, w10, w11, w1 + real matrix y0 + + ols_ipw_rc(0,y00,ixx00) + ols_ipw_rc(1,y01,ixx01) + ols_ipw_rc(2,y10,ixx10) + ols_ipw_rc(3,y11,ixx11) + + y0 = y00:*(-tmt:+1) + y01:*tmt + + w00 = wvar :* (tmt_trt:==0) :* psc :/(1 :- psc) ; w00 = w00:/mean(w00 ) + w01 = wvar :* (tmt_trt:==1) :* psc :/(1 :- psc) ; w01 = w01:/mean(w01 ) + w10 = wvar :* (tmt_trt:==2) ; w10 = w10:/mean(w10 ) + w11 = wvar :* (tmt_trt:==3) ; w11 = w11:/mean(w11 ) + w1 = wvar :* trt ; 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 + + // adds constant + real matrix y_y0 + xvar = xvar,J(nn,1,1) + y_y0 = yvar:-y0 + att_treat_pre = w10 :* (y_y0) ; eta_treat_pre = mean(att_treat_pre) + att_treat_post = w11 :* (y_y0) ; eta_treat_post = mean(att_treat_post) + att_cont_pre = w00 :* (y_y0) ; eta_cont_pre = mean(att_cont_pre) + att_cont_post = w01 :* (y_y0) ; eta_cont_post = mean(att_cont_post) + att_trt_post = w1 :* (y11 :- y01); eta_trt_post = mean(att_trt_post) + att_trtt1_post = w11 :* (y11 :- y01); eta_trtt1_post = mean(att_trtt1_post) + att_trt_pre = w1 :* (y10 :- y00); eta_trt_pre = mean(att_trt_pre) + att_trtt0_pre = w10 :* (y10 :- y00); 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 wgt00, XpX_inv_pre, lin_ols_pre, + wgt01, XpX_inv_post, lin_ols_post, + XpX_inv_pre_treat, lin_ols_pre_treat, + XpX_inv_post_treat, lin_ols_post_treat + + // cannot be simplified + wgt00 = wvar :* (tmt_trt:==0) + lin_ols_pre = ( wgt00 :* (yvar :- y00) :* xvar) * invsym(quadcross(xvar,wgt00, xvar)):*nn + + wgt01 = wvar :* (tmt_trt:==1) + lin_ols_post = (wgt01 :* (yvar :- y01) :* xvar) * invsym(quadcross(xvar,wgt01, xvar)):*nn + + //wols_x_pre_treat = w10 :* xvar + //wols_eX_pre_treat = w10 :* (y :- y10) :* xvar + + // These ones CAN + lin_ols_pre_treat = ( w10 :* (yvar :- y10) :* xvar) * ixx10 *nn + lin_ols_post_treat = (w11 :* (yvar :- y11) :* xvar) * ixx11 *nn + + 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 = (wvar :* (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 + + inf_treat_or_post = -lin_ols_post * mean(w11 :* tmt :* xvar)' + inf_treat_or_pre = -lin_ols_pre * mean(w10 :* (1 :- tmt) :* xvar)' + + real matrix inf_treat_or, inf_treat, inf_cont_post_pre + + //inf_treat_or = inf_treat_or_post :+ inf_treat_or_pre + inf_treat = inf_treat_post :- inf_treat_pre :+ + (inf_treat_or_post :+ inf_treat_or_pre) + + inf_cont_post_pre = (att_cont_post :- w01 :* eta_cont_post) :- + (att_cont_pre :- w00 :* eta_cont_pre) + + real matrix M2_pre, M2_post, inf_cont_ps, M3_post, M3_pre, inf_cont_or_post, inf_cont_or_pre + + inf_cont_ps = lin_rep_ps * + (mean(w01 :* (y_y0 :- eta_cont_post) :* xvar):- + mean(w00 :* (y_y0 :- eta_cont_pre) :* xvar))' + + real matrix inf_cont_or, inf_cont, trtr_eta_inf_func1 + + inf_cont_or = -lin_ols_post * mean(w01 :* tmt :* xvar)' :- + lin_ols_pre * mean(w00 :* (1 :- tmt) :* xvar)' + + inf_cont = inf_cont_post_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_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)) + + inf_or = (lin_ols_post_treat :- lin_ols_post) * mean((w1 :- w11) :* xvar)' :- + (lin_ols_pre_treat :- lin_ols_pre) * mean((w1 :- w10) :* xvar)' + + + rif = trtr_att :+ trtr_eta_inf_func1 :+ inf_eff :+ inf_or + } + } + +void drdid::drimp_rc(){ + // main Loading variables + real matrix psc, ipw + ipt() + if (conv==1) { + psc=logistic(xb) + ipw=psc:/(1:-psc) + /// tmt trt + real matrix y00, y01, y10, y11, + ixx00, ixx01, ixx10, ixx11, + w00 , w01, w10, w11, w1 + + ols_ipt_rc(0,ipw,y00,ixx00) + ols_ipt_rc(1,ipw,y01,ixx01) + ols_ipt_rc(2,1 ,y10,ixx10) + ols_ipt_rc(3,1 ,y11,ixx11) + + real matrix y0 + + y0 = y00:*(1:-tmt) :+ y01:*tmt + + + w00 = wvar :* (tmt_trt:==0) :* ipw; w00 = w00:/mean(w00 ) + w01 = wvar :* (tmt_trt:==1) :* ipw; w01 = w01:/mean(w01 ) + w10 = wvar :* (tmt_trt:==2) ; w10 = w10:/mean(w10 ) + w11 = wvar :* (tmt_trt:==3) ; w11 = w11:/mean(w11 ) + w1 = wvar :* trt ; 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 + 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 + real matrix y_y0 + y_y0 = yvar :- y0 + + att_treat_pre = w10 :* (y_y0) ; eta_treat_pre = mean(att_treat_pre) + att_treat_post = w11 :* (y_y0) ; eta_treat_post = mean(att_treat_post) + att_cont_pre = w00 :* (y_y0) ; eta_cont_pre = mean(att_cont_pre) + att_cont_post = w01 :* (y_y0) ; eta_cont_post = mean(att_cont_post) + att_trt_post = w1 :* (y11 :- y01); eta_trt_post = mean(att_trt_post) + att_trtt1_post = w11 :* (y11 :- y01); eta_trtt1_post = mean(att_trtt1_post) + att_trt_pre = w1 :* (y10 :- y00); eta_trt_pre = mean(att_trt_pre) + att_trtt0_pre = w10 :* (y10 :- y00); 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)) + + rif = trtr_att :+ att_inf_func1 :+ inf_eff + } +} + +void drdid::drimp2_rc(){ + // main Loading variables + real matrix psc, ipw + ipt() + xvar=xvar,J(nn,1,1) + if (conv==1) { + ipw=exp(xb) + /// tmt trt + + real matrix w10, w11, w00, w01, w1 + // Avoid Weights that are Too large + // Weights Scale is irrelevant + w00 = wvar :* (tmt_trt:==0) :*ipw + w01 = wvar :* (tmt_trt:==1) :*ipw + w10 = wvar :* (tmt_trt:==2) + w11 = wvar :* (tmt_trt:==3) + w1 = wvar :* trt + + xvar=xvar,J(nn,1,1) + + real matrix y00, y01, y10, y11, + ixx00, ixx01, ixx10, ixx11 + + // Need to get ixx for all cases + // as well as Predicted + ixx00 = invsym(cross(xvar,w00,xvar));y00=xvar*(ixx00*cross(xvar,w00,yvar)) + ixx01 = invsym(cross(xvar,w01,xvar));y01=xvar*(ixx01*cross(xvar,w01,yvar)) + ixx10 = invsym(cross(xvar,w10,xvar));y10=xvar*(ixx10*cross(xvar,w10,yvar)) + ixx11 = invsym(cross(xvar,w11,xvar));y11=xvar*(ixx11*cross(xvar,w11,yvar)) + + real matrix iff00, iff01, iff10 , iff11 + iff00 = (nn * ixx00*(xvar:*(yvar:-y00):*w00)')' + iff01 = (nn * ixx01*(xvar:*(yvar:-y01):*w01)')' + iff10 = (nn * ixx10*(xvar:*(yvar:-y10):*w10)')' + iff11 = (nn * ixx11*(xvar:*(yvar:-y11):*w11)')' + // IFs for all cases + real scalar nw1 + nw1 = nn/sum(w1) + + real matrix ifd00, ifd01, ifd10 , ifd11 + real scalar mn_y00, mn_y01,mn_y10,mn_y11, csumx + + mn_y00=mean(y00,w1);mn_y01=mean(y01,w1) + mn_y10=mean(y10,w1);mn_y11=mean(y11,w1) + csumx = colsum(xvar:*w1)' + + ifd00 = nw1*(w1:*(y00:-mn_y00):+ iff00 *csumx/nn) + ifd01 = nw1*(w1:*(y01:-mn_y01):+ iff01 *csumx/nn) + ifd10 = nw1*(w1:*(y10:-mn_y10):+ iff10 *csumx/nn) + ifd11 = nw1*(w1:*(y11:-mn_y11):+ iff11 *csumx/nn) + + real scalar att_gt + att_gt = mn_y11-mn_y10 - mn_y01+mn_y00 + + rif = (ifd11-ifd10)-(ifd01-ifd00):+att_gt + + } +} + // This is similar to Statas hdidregress +void drdid::reg2_rc(){ + // main Loading variables + real matrix w10, w11, w00, w01, w1 + // Avoid Weights that are Too large + // Weights Scale is irrelevant + w00 = wvar :* (tmt_trt:==0) + w01 = wvar :* (tmt_trt:==1) + w10 = wvar :* (tmt_trt:==2) + w11 = wvar :* (tmt_trt:==3) + w1 = wvar :* trt + + xvar=xvar,J(nn,1,1) + + real matrix y00, y01, y10, y11, + ixx00, ixx01, ixx10, ixx11 + + // Need to get ixx for all cases + // as well as Predicted + ixx00 = invsym(cross(xvar,w00,xvar));y00=xvar*(ixx00*cross(xvar,w00,yvar)) + ixx01 = invsym(cross(xvar,w01,xvar));y01=xvar*(ixx01*cross(xvar,w01,yvar)) + ixx10 = invsym(cross(xvar,w10,xvar));y10=xvar*(ixx10*cross(xvar,w10,yvar)) + ixx11 = invsym(cross(xvar,w11,xvar));y11=xvar*(ixx11*cross(xvar,w11,yvar)) + + real matrix iff00, iff01, iff10 , iff11 + iff00 = (nn * ixx00*(xvar:*(yvar:-y00):*w00)')' + iff01 = (nn * ixx01*(xvar:*(yvar:-y01):*w01)')' + iff10 = (nn * ixx10*(xvar:*(yvar:-y10):*w10)')' + iff11 = (nn * ixx11*(xvar:*(yvar:-y11):*w11)')' + // IFs for all cases + real scalar nw1 + nw1 = nn/sum(w1) + + real matrix ifd00, ifd01, ifd10 , ifd11 + real scalar mn_y00, mn_y01,mn_y10,mn_y11, csumx + + mn_y00=mean(y00,w1);mn_y01=mean(y01,w1) + mn_y10=mean(y10,w1);mn_y11=mean(y11,w1) + csumx = colsum(xvar:*w1)' + + ifd00 = nw1*(w1:*(y00:-mn_y00):+ iff00 *csumx/nn) + ifd01 = nw1*(w1:*(y01:-mn_y01):+ iff01 *csumx/nn) + ifd10 = nw1*(w1:*(y10:-mn_y10):+ iff10 *csumx/nn) + ifd11 = nw1*(w1:*(y11:-mn_y11):+ iff11 *csumx/nn) + + real scalar att_gt + att_gt = mn_y11-mn_y10 - mn_y01+mn_y00 + + rif = (ifd11-ifd10)-(ifd01-ifd00):+att_gt +} + + +void drdid::reg_rc() { + // main Loading variables + real matrix y00, y01, ixx00, ixx01 + + ols_ipt_rc(0,wvar,y00,ixx00) + ols_ipt_rc(1,wvar,y01,ixx01) + + // add constant + if (kx > 0) xvar=xvar,J(nn,1,1) + else xvar= J(nn,1,1) + + real matrix w10, w11, w1 + w10 = wvar :* trt :* (1 :- tmt);w10 = w10:/mean(w10 ) + w11 = wvar :* trt :* tmt ;w11 = w11:/mean(w11 ) + w1 = wvar :* trt ;w1 = w1 :/mean(w1 ) + + real matrix att_treat_pre, att_treat_post, att_cont, + eta_treat_pre, eta_treat_post, eta_cont, reg_att + + att_treat_pre = w10 :* yvar ; eta_treat_pre = mean(att_treat_pre) + att_treat_post = w11 :* yvar ; eta_treat_post = mean(att_treat_post) + att_cont = w1 :* (y01 :- y00); eta_cont = mean(att_cont) + + reg_att = (eta_treat_post :- eta_treat_pre) :- eta_cont + + real matrix w_ols_pre, wols_eX_pre, lin_rep_ols_pre + + w_ols_pre = wvar :* (1 :- trt) :* (1 :- tmt) + + wols_eX_pre = w_ols_pre :* (yvar :- y00) :* xvar + lin_rep_ols_pre = wols_eX_pre * ixx00 * nn + real matrix w_ols_post, wols_eX_post, lin_rep_ols_post + w_ols_post = wvar :* (1 :- trt) :* tmt + wols_eX_post = w_ols_post :* (yvar :- y01) :* xvar + lin_rep_ols_post= wols_eX_post * ixx01 :* nn + + real matrix inf_treat, inf_cont_1, inf_cont_2_post, inf_cont_2_pre, inf_control + + inf_treat = (att_treat_post :- w11 :* eta_treat_post) :- + (att_treat_pre :- w10 :* eta_treat_pre) + + inf_cont_1 = (att_cont :- w1 :* eta_cont) + + //M1 = mean(w0 :* xvar) + inf_cont_2_post = lin_rep_ols_post * mean(w1 :* xvar)' + inf_cont_2_pre = lin_rep_ols_pre * mean(w1 :* xvar)' + inf_control = (inf_cont_1 :+ inf_cont_2_post :- inf_cont_2_pre) + + rif = reg_att :+ (inf_treat :- inf_control) + + } + + +void drdid::stdipw_rc(){ + // main Loading variables + ilogit() + if (conv==1) { + real matrix psc, w00, w01, w10, w11 + psc=logistic(xb) + w00 = wvar :* (tmt_trt:==0) :* psc:/(1 :- psc) + w01 = wvar :* (tmt_trt:==1) :* psc:/(1 :- psc) + w10 = wvar :* (tmt_trt:==2) + w11 = wvar :* (tmt_trt:==3) + + 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 :* yvar; eta_treat_pre = mean(att_treat_pre) + att_treat_post = w11 :* yvar; eta_treat_post = mean(att_treat_post) + att_cont_pre = w00 :* yvar; eta_cont_pre = mean(att_cont_pre) + att_cont_post = w01 :* yvar; eta_cont_post = mean(att_cont_post) + // add constant + xvar=xvar,J(nn,1,1) + + 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 = (wvar :* (trt :- psc) :* xvar) * (psv :* nn) + + real matrix inf_treat, inf_cont, 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 :* (yvar :- eta_cont_pre) :* xvar) + //M2_post = mean(w01 :* (yvar :- eta_cont_post):* xvar) + inf_cont_ps = lin_rep_ps * ( mean(w01 :* (yvar :- eta_cont_post):* xvar) :- + mean(w00 :* (yvar :- eta_cont_pre) :* xvar))' + //inf_cont = inf_cont :+ inf_cont_ps + // ipw_att = (eta_treat_post :- eta_treat_pre) :- + // (eta_cont_post :- eta_cont_pre) + + rif = (eta_treat_post :- eta_treat_pre) :- + (eta_cont_post :- eta_cont_pre) :+ + inf_treat :- ( inf_cont :+ inf_cont_ps ) + + + } + } + +/// SETUP Interactive Version + +/*void drdid::setup(){ + yvar = st_data(.,"re") + xvar = st_data(.,"age educ black married nodegree hisp re74") + tmt = st_data(.,"year") + trt = st_data(.,"experimental") + wvar = J(rows(yvar),1,1) + id = st_data(.,"id") + oid = st_data(.,"id") +} + +void drdid::setup2(real scalar dt, real scalar mt){ + data_type = dt + method_type = mt + yvar = st_data(.,"re") + xvar = st_data(.,"age educ black married nodegree hisp re74") + tmt = st_data(.,"year") + trt = st_data(.,"experimental") + wvar = J(rows(yvar),1,1) + id = st_data(.,"id") + oid = st_data(.,"id") +}*/ + +// Non Interactive. +// This IMPLEMENTS DRDID +void drdid::drdid(){ + // setup2(dt,mt) + // IF N=0 NOT CONVERGED. eRROR + if (data_type == 1) { + + if (rolljw==0) msetup_panel() + else if (rolljw==1) msetup_panel2() + + if (err ==0) { + if (method_type ==1) dripw_panel() + else if (method_type ==2) drimp_panel() + else if (method_type ==3) stdipw_panel() + else if (method_type ==4) reg_panel() + } + else conv=0 + } + else { + // IF RC 6 methods. Method 6 may be redundant + msetup_rc() + if (method_type ==1) dripw_rc() + else if (method_type ==2) drimp_rc() + else if (method_type ==3) stdipw_rc() + else if (method_type ==4) reg_rc() + else if (method_type ==5) reg2_rc() + else if (method_type ==6) drimp2_rc() + } + + if (conv==1) minn=rows(rif)*(1+(data_type:==1)) + else minn=0 + +} + +void drdid::csdid_setup(real matrix syvar, sxvar , + stvar, sgvar, swvar, sivar, soid, + real scalar dt, mt){ + yvar = syvar + xvar = sxvar + tmt = stvar + trt = sgvar + wvar = swvar + id = sivar + oid = soid + data_type = dt + method_type = mt + conv=1 + +} + + + +/*void drdid::drdid_outpt(real matrix xout, minx, cnv){ + //xout = oid, id[,1], tmt, trt, wvar, rif + //cnv = conv + //minx=rows(rif)*(1+data_type:=1) +}*/ + + + +end diff --git a/csdid2 - Copy/lcsdid.mlib b/csdid2 - Copy/lcsdid.mlib new file mode 100644 index 0000000..8f8aec0 Binary files /dev/null and b/csdid2 - Copy/lcsdid.mlib differ diff --git a/csdid2 - Copy/saux.dta b/csdid2 - Copy/saux.dta new file mode 100644 index 0000000..415ed77 Binary files /dev/null and b/csdid2 - Copy/saux.dta differ diff --git a/csdid2 - Copy/spcsdid.do b/csdid2 - Copy/spcsdid.do new file mode 100644 index 0000000..85e658e --- /dev/null +++ b/csdid2 - Copy/spcsdid.do @@ -0,0 +1,13 @@ +mata + class spcsdid { + real matrix index + real matrix attgt + real scalar gvar + real scalar tvar + } + pmn = spcsdid(2) + pmn[1].attgt = runiform(6,1,0,1) + pmn[1].index = (1,4,5,6,10,12)' + pmn[2].attgt = runiform(6,1,0,1) + pmn[2].index = (2,3,7,10,11,12)' +end \ No newline at end of file diff --git a/csdid2 - Copy/unit_test.do b/csdid2 - Copy/unit_test.do new file mode 100644 index 0000000..5454bdc --- /dev/null +++ b/csdid2 - Copy/unit_test.do @@ -0,0 +1,489 @@ +frause mpdta, clear +gen st = floor(county/1000) +mata mata clear +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\drdid.mata" +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\csdid.mata" +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\csdid_stats.mata" + + +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +assert `"`e(cmdline)'"' == `"csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .017595643753745 +mat T_b[1,2] = -.0772398214573161 +mat T_b[1,3] = .003306356692512 +mat T_b[1,4] = .0250218295975542 +mat T_b[1,5] = .024458744971169 +mat T_b[1,6] = -.0199318167892598 +mat T_b[1,7] = -.050957367065195 +mat T_b[1,8] = -.1372587388894044 +mat T_b[1,9] = -.1008113630854053 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .000286290277302 +mat T_V[1,2] = -.0000699197196847 +mat T_V[1,3] = .0003815389226184 +mat T_V[1,4] = .0002804513165675 +mat T_V[1,5] = .0001968805927202 +mat T_V[1,6] = .0000333051709288 +mat T_V[1,7] = -.0000551650679429 +mat T_V[1,8] = -.0001468122320985 +mat T_V[1,9] = -.0001110067496262 +mat T_V[2,1] = -.0000699197196847 +mat T_V[2,2] = .0003986007882398 +mat T_V[2,3] = -.000109822195412 +mat T_V[2,4] = -.0000581427116154 +mat T_V[2,5] = -.0000417942520266 +mat T_V[2,6] = .0000776843638391 +mat T_V[2,7] = .0002581489637783 +mat T_V[2,8] = .0006493134095847 +mat T_V[2,9] = .0006092564157568 +mat T_V[3,1] = .0003815389226184 +mat T_V[3,2] = -.000109822195412 +mat T_V[3,3] = .0005978940904662 +mat T_V[3,4] = .0003359073547961 +mat T_V[3,5] = .0002108153225929 +mat T_V[3,6] = .0000284213133248 +mat T_V[3,7] = -.0000824643289709 +mat T_V[3,8] = -.0002135535062099 +mat T_V[3,9] = -.0001716922597921 +mat T_V[4,1] = .0002804513165675 +mat T_V[4,2] = -.0000581427116154 +mat T_V[4,3] = .0003359073547961 +mat T_V[4,4] = .0003282952872387 +mat T_V[4,5] = .0001771513076679 +mat T_V[4,6] = .0000382262676265 +mat T_V[4,7] = -.000042674305463 +mat T_V[4,8] = -.0001311357514049 +mat T_V[4,9] = -.0000969870572202 +mat T_V[5,1] = .0001968805927202 +mat T_V[5,2] = -.0000417942520266 +mat T_V[5,3] = .0002108153225929 +mat T_V[5,4] = .0001771513076679 +mat T_V[5,5] = .0002026751478997 +mat T_V[5,6] = .0000332679318352 +mat T_V[5,7] = -.0000403565693949 +mat T_V[5,8] = -.0000957474386806 +mat T_V[5,9] = -.0000643409318662 +mat T_V[6,1] = .0000333051709288 +mat T_V[6,2] = .0000776843638391 +mat T_V[6,3] = .0000284213133248 +mat T_V[6,4] = .0000382262676265 +mat T_V[6,5] = .0000332679318352 +mat T_V[6,6] = .0001398628868337 +mat T_V[6,7] = .0000615316221245 +mat T_V[6,8] = .0000437211453232 +mat T_V[6,9] = .0000656218010752 +mat T_V[7,1] = -.0000551650679429 +mat T_V[7,2] = .0002581489637783 +mat T_V[7,3] = -.0000824643289709 +mat T_V[7,4] = -.000042674305463 +mat T_V[7,5] = -.0000403565693949 +mat T_V[7,6] = .0000615316221245 +mat T_V[7,7] = .0002853895404404 +mat T_V[7,8] = .0003604010457427 +mat T_V[7,9] = .0003252736468058 +mat T_V[8,1] = -.0001468122320985 +mat T_V[8,2] = .0006493134095847 +mat T_V[8,3] = -.0002135535062099 +mat T_V[8,4] = -.0001311357514049 +mat T_V[8,5] = -.0000957474386806 +mat T_V[8,6] = .0000437211453232 +mat T_V[8,7] = .0003604010457427 +mat T_V[8,8] = .001327557632085 +mat T_V[8,9] = .0008655738151882 +mat T_V[9,1] = -.0001110067496262 +mat T_V[9,2] = .0006092564157568 +mat T_V[9,3] = -.0001716922597921 +mat T_V[9,4] = -.0000969870572202 +mat T_V[9,5] = -.0000643409318662 +mat T_V[9,6] = .0000656218010752 +mat T_V[9,7] = .0003252736468058 +mat T_V[9,8] = .0008655738151882 +mat T_V[9,9] = .0011805563999581 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + + +csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +assert `"`e(cmdline)'"' == `"csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .0169262437014923 +mat T_b[1,2] = -.069393526595745 +mat T_b[1,3] = .007932802034353 +mat T_b[1,4] = .0214173492404631 +mat T_b[1,5] = .0214285798296607 +mat T_b[1,6] = -.0237705768478331 +mat T_b[1,7] = -.0407911825853478 +mat T_b[1,8] = -.1153380429226659 +mat T_b[1,9] = -.0976743040271333 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0002311442272515 +mat T_V[1,2] = -.0000913516255354 +mat T_V[1,3] = .0002942556613254 +mat T_V[1,4] = .0002300102113643 +mat T_V[1,5] = .0001691668090648 +mat T_V[1,6] = 1.97724945317e-06 +mat T_V[1,7] = -.0000719224656945 +mat T_V[1,8] = -.0001704115735685 +mat T_V[1,9] = -.0001250497123318 +mat T_V[2,1] = -.0000913516255354 +mat T_V[2,2] = .0002956989691059 +mat T_V[2,3] = -.0001229611368996 +mat T_V[2,4] = -.0000827218925021 +mat T_V[2,5] = -.0000683718472045 +mat T_V[2,6] = .0000497459253988 +mat T_V[2,7] = .0001986659277167 +mat T_V[2,8] = .0005048986425726 +mat T_V[2,9] = .0004294853807356 +mat T_V[3,1] = .0002942556613254 +mat T_V[3,2] = -.0001229611368996 +mat T_V[3,3] = .0004321568114057 +mat T_V[3,4] = .0002650905107159 +mat T_V[3,5] = .0001855196618545 +mat T_V[3,6] = 4.29689321764e-06 +mat T_V[3,7] = -.0000922749994412 +mat T_V[3,8] = -.0002301133832139 +mat T_V[3,9] = -.0001737530581611 +mat T_V[4,1] = .0002300102113643 +mat T_V[4,2] = -.0000827218925021 +mat T_V[4,3] = .0002650905107159 +mat T_V[4,4] = .0002623370288506 +mat T_V[4,5] = .0001626030945263 +mat T_V[4,6] = 1.91446099020e-06 +mat T_V[4,7] = -.0000640406288387 +mat T_V[4,8] = -.0001547330524719 +mat T_V[4,9] = -.0001140283496882 +mat T_V[5,1] = .0001691668090648 +mat T_V[5,2] = -.0000683718472045 +mat T_V[5,3] = .0001855196618545 +mat T_V[5,4] = .0001626030945263 +mat T_V[5,5] = .0001593776708136 +mat T_V[5,6] = -2.79605848321e-07 +mat T_V[5,7] = -.0000594517688036 +mat T_V[5,8] = -.0001263882850197 +mat T_V[5,9] = -.0000873677291463 +mat T_V[6,1] = 1.97724945317e-06 +mat T_V[6,2] = .0000497459253988 +mat T_V[6,3] = 4.29689321764e-06 +mat T_V[6,4] = 1.91446099020e-06 +mat T_V[6,5] = -2.79605848321e-07 +mat T_V[6,6] = .0000730459300969 +mat T_V[6,7] = .0000481226655546 +mat T_V[6,8] = .000031455063973 +mat T_V[6,9] = .0000463600419708 +mat T_V[7,1] = -.0000719224656945 +mat T_V[7,2] = .0001986659277167 +mat T_V[7,3] = -.0000922749994412 +mat T_V[7,4] = -.0000640406288387 +mat T_V[7,5] = -.0000594517688036 +mat T_V[7,6] = .0000481226655546 +mat T_V[7,7] = .0002266976249852 +mat T_V[7,8] = .0002812880537518 +mat T_V[7,9] = .0002385553665754 +mat T_V[8,1] = -.0001704115735685 +mat T_V[8,2] = .0005048986425726 +mat T_V[8,3] = -.0002301133832139 +mat T_V[8,4] = -.0001547330524719 +mat T_V[8,5] = -.0001263882850197 +mat T_V[8,6] = .000031455063973 +mat T_V[8,7] = .0002812880537518 +mat T_V[8,8] = .0010145967593686 +mat T_V[8,9] = .0006922546931968 +mat T_V[9,1] = -.0001250497123318 +mat T_V[9,2] = .0004294853807356 +mat T_V[9,3] = -.0001737530581611 +mat T_V[9,4] = -.0001140283496882 +mat T_V[9,5] = -.0000873677291463 +mat T_V[9,6] = .0000463600419708 +mat T_V[9,7] = .0002385553665754 +mat T_V[9,8] = .0006922546931968 +mat T_V[9,9] = .0007407714211994 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + + +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st) + + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(cluster)'"' == `"st"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +_assert_streq `"`e(cmdline)'"' `"csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .017595643753745 +mat T_b[1,2] = -.0772398214573161 +mat T_b[1,3] = .003306356692512 +mat T_b[1,4] = .0250218295975542 +mat T_b[1,5] = .024458744971169 +mat T_b[1,6] = -.0199318167892598 +mat T_b[1,7] = -.050957367065195 +mat T_b[1,8] = -.1372587388894044 +mat T_b[1,9] = -.1008113630854053 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0007829071217662 +mat T_V[1,2] = -.0002717218539991 +mat T_V[1,3] = .0010315906382616 +mat T_V[1,4] = .000760194622637 +mat T_V[1,5] = .0005569361044 +mat T_V[1,6] = -.0001115421517223 +mat T_V[1,7] = -.0002578712924792 +mat T_V[1,8] = -.0003893609598156 +mat T_V[1,9] = -.0003281130119792 +mat T_V[2,1] = -.0002717218539991 +mat T_V[2,2] = .0002183504515259 +mat T_V[2,3] = -.0003294819691427 +mat T_V[2,4] = -.0002994560390609 +mat T_V[2,5] = -.0001862275537937 +mat T_V[2,6] = .0000568837146026 +mat T_V[2,7] = .0002265694813186 +mat T_V[2,8] = .0003051592816874 +mat T_V[2,9] = .0002847893284949 +mat T_V[3,1] = .0010315906382616 +mat T_V[3,2] = -.0003294819691427 +mat T_V[3,3] = .0015201009874528 +mat T_V[3,4] = .0009289254297727 +mat T_V[3,5] = .0006457454975592 +mat T_V[3,6] = -.000110003095552 +mat T_V[3,7] = -.0002160765790923 +mat T_V[3,8] = -.0005383288323113 +mat T_V[3,9] = -.0004535193696153 +mat T_V[4,1] = .000760194622637 +mat T_V[4,2] = -.0002994560390609 +mat T_V[4,3] = .0009289254297727 +mat T_V[4,4] = .0008118381067114 +mat T_V[4,5] = .0005398203314268 +mat T_V[4,6] = -.0001164242778234 +mat T_V[4,7] = -.0003698006047936 +mat T_V[4,8] = -.0003826496670879 +mat T_V[4,9] = -.0003289496065386 +mat T_V[5,1] = .0005569361044 +mat T_V[5,2] = -.0001862275537937 +mat T_V[5,3] = .0006457454975592 +mat T_V[5,4] = .0005398203314268 +mat T_V[5,5] = .0004852424842139 +mat T_V[5,6] = -.0001081990817914 +mat T_V[5,7] = -.0001877366935518 +mat T_V[5,8] = -.0002471043800477 +mat T_V[5,9] = -.0002018700597838 +mat T_V[6,1] = -.0001115421517223 +mat T_V[6,2] = .0000568837146026 +mat T_V[6,3] = -.000110003095552 +mat T_V[6,4] = -.0001164242778234 +mat T_V[6,5] = -.0001081990817914 +mat T_V[6,6] = .0001057850796644 +mat T_V[6,7] = .0000670490666338 +mat T_V[6,8] = .0000127123457306 +mat T_V[6,9] = .0000419883663815 +mat T_V[7,1] = -.0002578712924792 +mat T_V[7,2] = .0002265694813186 +mat T_V[7,3] = -.0002160765790923 +mat T_V[7,4] = -.0003698006047936 +mat T_V[7,5] = -.0001877366935518 +mat T_V[7,6] = .0000670490666338 +mat T_V[7,7] = .0004120548254679 +mat T_V[7,8] = .0002160765790923 +mat T_V[7,9] = .0002110974540805 +mat T_V[8,1] = -.0003893609598156 +mat T_V[8,2] = .0003051592816874 +mat T_V[8,3] = -.0005383288323113 +mat T_V[8,4] = -.0003826496670879 +mat T_V[8,5] = -.0002471043800477 +mat T_V[8,6] = .0000127123457306 +mat T_V[8,7] = .0002160765790923 +mat T_V[8,8] = .0005383288323113 +mat T_V[8,9] = .0004535193696153 +mat T_V[9,1] = -.0003281130119792 +mat T_V[9,2] = .0002847893284949 +mat T_V[9,3] = -.0004535193696153 +mat T_V[9,4] = -.0003289496065386 +mat T_V[9,5] = -.0002018700597838 +mat T_V[9,6] = .0000419883663815 +mat T_V[9,7] = .0002110974540805 +mat T_V[9,8] = .0004535193696153 +mat T_V[9,9] = .0004325521239021 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + +csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st) + + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(cluster)'"' == `"st"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +_assert_streq `"`e(cmdline)'"' `"csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .0169262437014923 +mat T_b[1,2] = -.069393526595745 +mat T_b[1,3] = .007932802034353 +mat T_b[1,4] = .0214173492404631 +mat T_b[1,5] = .0214285798296607 +mat T_b[1,6] = -.0237705768478331 +mat T_b[1,7] = -.0407911825853478 +mat T_b[1,8] = -.1153380429226659 +mat T_b[1,9] = -.0976743040271333 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0008553362098125 +mat T_V[1,2] = -.0003729056186157 +mat T_V[1,3] = .0011344934983186 +mat T_V[1,4] = .0008256904623471 +mat T_V[1,5] = .0006058246687719 +mat T_V[1,6] = -.0001762809606567 +mat T_V[1,7] = -.0003472355983124 +mat T_V[1,8] = -.0005319225374714 +mat T_V[1,9] = -.0004361833780222 +mat T_V[2,1] = -.0003729056186157 +mat T_V[2,2] = .0002955924619583 +mat T_V[2,3] = -.0004470501707058 +mat T_V[2,4] = -.0004022834464612 +mat T_V[2,5] = -.0002693832386801 +mat T_V[2,6] = .0000797094956907 +mat T_V[2,7] = .0003083500186895 +mat T_V[2,8] = .0004160569766317 +mat T_V[2,9] = .0003782533568214 +mat T_V[3,1] = .0011344934983186 +mat T_V[3,2] = -.0004470501707058 +mat T_V[3,3] = .0016447127057925 +mat T_V[3,4] = .0010073691441048 +mat T_V[3,5] = .0007513986450585 +mat T_V[3,6] = -.0001769648458234 +mat T_V[3,7] = -.0002916471911943 +mat T_V[3,8] = -.0007284565174212 +mat T_V[3,9] = -.0005911321283841 +mat T_V[4,1] = .0008256904623471 +mat T_V[4,2] = -.0004022834464612 +mat T_V[4,3] = .0010073691441048 +mat T_V[4,4] = .0008715078927645 +mat T_V[4,5] = .0005981943501718 +mat T_V[4,6] = -.0001981435967092 +mat T_V[4,7] = -.0004765814799506 +mat T_V[4,8] = -.0005077733730396 +mat T_V[4,9] = -.0004266353361453 +mat T_V[5,1] = .0006058246687719 +mat T_V[5,2] = -.0002693832386801 +mat T_V[5,3] = .0007513986450585 +mat T_V[5,4] = .0005981943501718 +mat T_V[5,5] = .0004678810110854 +mat T_V[5,6] = -.0001537344394376 +mat T_V[5,7] = -.0002734781237922 +mat T_V[5,8] = -.0003595377219535 +mat T_V[5,9] = -.0002907826695372 +mat T_V[6,1] = -.0001762809606567 +mat T_V[6,2] = .0000797094956907 +mat T_V[6,3] = -.0001769648458234 +mat T_V[6,4] = -.0001981435967092 +mat T_V[6,5] = -.0001537344394376 +mat T_V[6,6] = .0001225240012724 +mat T_V[6,7] = .0001146947512824 +mat T_V[6,8] = .0000207157532791 +mat T_V[6,9] = .0000609034769289 +mat T_V[7,1] = -.0003472355983124 +mat T_V[7,2] = .0003083500186895 +mat T_V[7,3] = -.0002916471911943 +mat T_V[7,4] = -.0004765814799506 +mat T_V[7,5] = -.0002734781237922 +mat T_V[7,6] = .0001146947512824 +mat T_V[7,7] = .0005344335781815 +mat T_V[7,8] = .0002974894749616 +mat T_V[7,9] = .0002867822703325 +mat T_V[8,1] = -.0005319225374714 +mat T_V[8,2] = .0004160569766317 +mat T_V[8,3] = -.0007284565174212 +mat T_V[8,4] = -.0005077733730396 +mat T_V[8,5] = -.0003595377219535 +mat T_V[8,6] = .0000207157532791 +mat T_V[8,7] = .0002974894749616 +mat T_V[8,8] = .000743048976445 +mat T_V[8,9] = .0006029737018409 +mat T_V[9,1] = -.0004361833780222 +mat T_V[9,2] = .0003782533568214 +mat T_V[9,3] = -.0005911321283841 +mat T_V[9,4] = -.0004266353361453 +mat T_V[9,5] = -.0002907826695372 +mat T_V[9,6] = .0000609034769289 +mat T_V[9,7] = .0002867822703325 +mat T_V[9,8] = .0006029737018409 +mat T_V[9,9] = .0005623539781832 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + diff --git a/csdid2 - Copy/wboot_sbs.do b/csdid2 - Copy/wboot_sbs.do new file mode 100644 index 0000000..5006975 --- /dev/null +++ b/csdid2 - Copy/wboot_sbs.do @@ -0,0 +1,66 @@ +frause oaxaca, clear +drop if lnwage==. + +mata +y = st_data(.,"lnwage"); n=rows(y) +x = st_data(.,"educ exper tenure female"),J(n,1,1); k =cols(x) + +//beta + +xx = cross(x,x); ixx = invsym(xx); xy = cross(x,y) +beta = ixx * xy + +// errors + +e = y:-x*beta +end + +** This is what you do with asymptotic theory: +** First lets get the IF +mata + iff = (x:*e)*ixx*n + vcv = cross(iff,iff)/n^2 + //sd + sd = sqrt(diagonal(vcv)) +end + * you can compare this with Reg, robust +mata:sd + +// Now the **Wboot done in CSDID: + +mata +// need a place to store ALL betas +beta = J(999,k,.) +for(i=1;i<=999;i++){ + beta[i,]=mean(iff:*rnormal(n,k,0,1)) +} +// This is the Equivalent to the bootstrapped coefficients +// Also i could have used a different noise "v". But normal is just easiest to program + +// Get SE using the Interquartile difference +iqrs=J(1,k,.) +for(i=1;i<=k;i++){ + ///999 is the Nrows in Beta + aux=sort(beta,i) + q25=ceil((999+1)*.25) + q75=ceil((999+1)*.75) + // Thse are the SE reported in CSDID when wboot is used + iqrs[,i]=(aux[q75,i]-aux[q25,i]):/(invnormal(.75)-invnormal(.25) ) +} + +// Get the t-stats based on Beta and IQRS + +tstat = beta:/iqrs + +// This is important. To get the Uniform CI, you need to first the the distribution of the "largest absolute tstat". +// the Absolute part makes this symetric + +largest_tstat = rowmax(abs(tstat)) + +// Finally the Critical is based by looking at the t-stat in this vector that is in the 95th possition +_sort(largest_tstat ,1) +largest_tstat[ceil((999+1)*.95)] +// CIs are constructed using this largest tstat +end + + diff --git a/csdid2/csdid.mata b/csdid2/csdid.mata index 38ebfbd..024ae3a 100644 --- a/csdid2/csdid.mata +++ b/csdid2/csdid.mata @@ -186,11 +186,13 @@ void csdid::csdid_setup(){ // antici for Anticipation antici = 0 rolljw = 0 + if (length(ivar)==0) { // Repeated CrosSection type_data = 2 oid = 1::rows(yvar) ivar = oid + // Ivar here does not matter, its just a sorting variable if (length(cvar)>0) ord = order((cvar,gvar,tvar,ivar),(1,2,3,4)) else ord = order((gvar,tvar,ivar),(1,2,3)) // If data was already loaded sorted. No need to Resort @@ -212,15 +214,19 @@ void csdid::csdid_setup(){ // if panel, first sort type_data = 1 oid = 1::rows(yvar) - ord = order((ivar,tvar),(1,2)) + if (length(cvar)>0) ord = order((cvar,ivar,tvar),(1,2,3)) + else ord = order((ivar,tvar),(1,2)) + if (ord!=oid) { + // Order all variables in case order changed yvar=yvar[ord,] - if (length(xvar)>0) xvar=xvar[ord,] tvar=tvar[ord,] gvar=gvar[ord,] + ivar=ivar[ord,] + // This runs only if they exist + if (length(xvar)>0) xvar=xvar[ord,] if (length(wvar)>0) wvar=wvar[ord,] - ivar=ivar[ord,] - if (length(cvar)>0) cvar=cvar[ord,] + if (length(cvar)>0) cvar=cvar[ord,] } // then recode @@ -236,7 +242,7 @@ void csdid::makeid(){ real scalar i,j, in real matrix aux in = rows(ivar) - oid=runningsum(0\(ivar[2..in,]-ivar[1..in-1,]):>0):+1 + oid=runningsum(0\(ivar[2..in,]-ivar[1..in-1,]):!=0):+1 // recoded //oid=uniqrows(id2,1)[id2,] } @@ -361,24 +367,35 @@ real matrix csdid::nvalid(real matrix sgtvar, real matrix csdid::nsample_select() { real scalar i real matrix tsel, gsel, sgtvar - real scalar gv, tv + real scalar gv, tv, tv0, tv1 real matrix toreturn real matrix eftime toreturn = J(rows(fgtvar),7,0) // event eftime=J(1,2,0) + + // This is to determine Obs for(i=1;i<=rows(fgtvar);i++){ - gv = fgtvar[i,1];tv = fgtvar[i,2] - + + gv = fgtvar[i,1]; tv = fgtvar[i,2] + + /// Group Selection if (notyet==0) gsel = (ogtvar[,1]:==0 :| ogtvar[,1]:==gv) else { - gsel = (ogtvar[,1]:==0 :| ogtvar[,1]:>max((gv,tv)) :| ogtvar[,1]:==gv) - // Check asin R above. This is the same but with N (counting Obs) + gsel = (ogtvar[,1]:==0 :| ogtvar[,1]:>max((gv,tv)) :| ogtvar[,1]:==gv) + // Check asin R: Include other units not treated at tv if ((asinr==1) & (tvtv0 :| ogvar:==gv) - else gsel = (ogvar:==0 :| ogvar:>tv1 :| ogvar:==gv) + if (shortx==0) { + // If long pre-att + // gsel = ( (gvar:==0) :| (gvar:>tv0) :| (gvar:==gv)) + gsel = ( (gvar:==0) :| (gvar :> tv0) :| (gvar:==gv)) + } + else { + // If short pre-att + gsel = ( (gvar:==0) :| (gvar :> tv1) :| (gvar:==gv)) + } } } @@ -390,11 +407,13 @@ real matrix csdid::nsample_select() { } else { if (shortx==0 ) { + // Long Pre ATT tsel = (ogtvar[,2]:==gv-delta :| ogtvar[,2]:==tv-delta) /// T0 T1 T eftime=( tv-delta ,gv-delta, tv-delta ) } else { + // Short Pre ATT tsel = (ogtvar[,2]:==tv-delta :| ogtvar[,2]:==tv) /// T0 T1 T eftime=(tv-delta ,tv , tv-delta) @@ -450,6 +469,7 @@ void csdid::csdid(){ drdid.trt =gvar[smsel,]:==gv //if (cols(wvar)>0) drdid.wvar=select(wvar,smsel) if (cols(wvar)>0) drdid.wvar=wvar[smsel,] + //drdid.id =select(ivar,smsel) //drdid.oid =select(oid ,smsel) drdid.id =ivar[smsel,] @@ -465,8 +485,11 @@ void csdid::csdid(){ if (shortx==0) frif[drdid.oid,i]=drdid.rif:*sign(eventvar[i]+.01) else frif[drdid.oid,i]=drdid.rif frwt[drdid.oid,i]=drdid.wvar:*drdid.trt + + //:*drdid.trt } + dots = dots-convar[i,] stata(sprintf("_dots %f %f",i,dots)) } @@ -495,14 +518,15 @@ void csdid::csdid(){ real matrix aux if (length(cvar)>0) { - aux = oid, gvar, wvar, cvar + aux = cvar, oid, gvar, wvar aux=uniqrows(aux) - ord = order(aux,(4,1,2,3)) + ord = order(aux,(1,2,3,4)) aux = aux[ord,] - oid = aux[,1] - gvar= aux[,2] - wvar= aux[,3] - cvar= aux[,4] + cvar= aux[,1] + oid = aux[,2] + gvar= aux[,3] + wvar= aux[,4] + frwt = frwt[ord,] frif = frif[ord,] } diff --git a/csdid2/csdid_stats.mata b/csdid2/csdid_stats.mata index 6e675ac..58d3a56 100644 --- a/csdid2/csdid_stats.mata +++ b/csdid2/csdid_stats.mata @@ -108,6 +108,7 @@ real matrix csdid_estat::aggte(real matrix rif,| real matrix wgt ) { return(rowsum(r1):+rowsum(r2):-rowsum(r3):+mn_all) } + // Will use Separate function for WB bc it process data differently void csdid_estat::atts_asym(class csdid scalar csdid){ @@ -516,7 +517,7 @@ void csdid_estat::cevent_att(class csdid scalar csdid ){ //////////////////////////////////////////////////////////////////////////////// /// event Aggregations //////////////////////////////////////////////////////////////////////////////// - + void csdid_estat::event_att(class csdid scalar csdid){ // Counter i // kgroups (max) @@ -576,7 +577,8 @@ void csdid_estat::event_att(class csdid scalar csdid){ onames[iic,] = sprintf("tp%f", abs(csdid.sevent[i]:-csdid.antici)) } - aux_wgt = select(csdid.frwt,toselect) + aux_wgt = select(csdid.frwt,toselect) + aux[,iic-ievent] = aggte(select(csdid.frif,toselect),aux_wgt ) //sumwgt[i,] = rowsum(aux_wgt):/cols(aux_wgt) diff --git a/csdid2/drdid.mata b/csdid2/drdid.mata index a63dd73..a3e2750 100644 --- a/csdid2/drdid.mata +++ b/csdid2/drdid.mata @@ -401,7 +401,7 @@ void drdid::makeid(){ j = 1 for(i=1;i<=rows(id);i++){ if (i>1) { - if (id[i]>id[i-1]) { + if (id[i]:!=id[i-1]) { j++ } } diff --git a/csdid2/saux.dta b/csdid2/saux.dta new file mode 100644 index 0000000..415ed77 Binary files /dev/null and b/csdid2/saux.dta differ diff --git a/csdid2/unit_test.do b/csdid2/unit_test.do new file mode 100644 index 0000000..252991e --- /dev/null +++ b/csdid2/unit_test.do @@ -0,0 +1,503 @@ +frause mpdta, clear +gen st = floor(county/1000) +gen w = runiformint(1,7) +bysort county:replace w=w[1] +mata mata clear +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\drdid.mata" +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\csdid.mata" +run "C:\Users\Fernando\Documents\GitHub\stpackages\csdid2\csdid_stats.mata" +csdid lemp , ivar( countyreal) time(year) gvar(first) agg(attgt) method(reg) long2 +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(attgt) method(reg) + + +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) +csdid lemp , ivar( countyreal) time(year) gvar(first) agg(event) long2 + +csdid2 lemp lpop , ivar( countyreal) tvar(year) gvar(first) agg(event) method(reg) +csdid lemp lpop , ivar( countyreal) time(year) gvar(first) agg(event) method(reg) long2 + +csdid2 lemp [w=w] , ivar( countyreal) tvar(year) gvar(first) agg(attgt) method(reg) + +csdid lemp [w=w], ivar( countyreal) time(year) gvar(first) agg(attgt) method(reg) long2 + +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +assert `"`e(cmdline)'"' == `"csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .017595643753745 +mat T_b[1,2] = -.0772398214573161 +mat T_b[1,3] = .003306356692512 +mat T_b[1,4] = .0250218295975542 +mat T_b[1,5] = .024458744971169 +mat T_b[1,6] = -.0199318167892598 +mat T_b[1,7] = -.050957367065195 +mat T_b[1,8] = -.1372587388894044 +mat T_b[1,9] = -.1008113630854053 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .000286290277302 +mat T_V[1,2] = -.0000699197196847 +mat T_V[1,3] = .0003815389226184 +mat T_V[1,4] = .0002804513165675 +mat T_V[1,5] = .0001968805927202 +mat T_V[1,6] = .0000333051709288 +mat T_V[1,7] = -.0000551650679429 +mat T_V[1,8] = -.0001468122320985 +mat T_V[1,9] = -.0001110067496262 +mat T_V[2,1] = -.0000699197196847 +mat T_V[2,2] = .0003986007882398 +mat T_V[2,3] = -.000109822195412 +mat T_V[2,4] = -.0000581427116154 +mat T_V[2,5] = -.0000417942520266 +mat T_V[2,6] = .0000776843638391 +mat T_V[2,7] = .0002581489637783 +mat T_V[2,8] = .0006493134095847 +mat T_V[2,9] = .0006092564157568 +mat T_V[3,1] = .0003815389226184 +mat T_V[3,2] = -.000109822195412 +mat T_V[3,3] = .0005978940904662 +mat T_V[3,4] = .0003359073547961 +mat T_V[3,5] = .0002108153225929 +mat T_V[3,6] = .0000284213133248 +mat T_V[3,7] = -.0000824643289709 +mat T_V[3,8] = -.0002135535062099 +mat T_V[3,9] = -.0001716922597921 +mat T_V[4,1] = .0002804513165675 +mat T_V[4,2] = -.0000581427116154 +mat T_V[4,3] = .0003359073547961 +mat T_V[4,4] = .0003282952872387 +mat T_V[4,5] = .0001771513076679 +mat T_V[4,6] = .0000382262676265 +mat T_V[4,7] = -.000042674305463 +mat T_V[4,8] = -.0001311357514049 +mat T_V[4,9] = -.0000969870572202 +mat T_V[5,1] = .0001968805927202 +mat T_V[5,2] = -.0000417942520266 +mat T_V[5,3] = .0002108153225929 +mat T_V[5,4] = .0001771513076679 +mat T_V[5,5] = .0002026751478997 +mat T_V[5,6] = .0000332679318352 +mat T_V[5,7] = -.0000403565693949 +mat T_V[5,8] = -.0000957474386806 +mat T_V[5,9] = -.0000643409318662 +mat T_V[6,1] = .0000333051709288 +mat T_V[6,2] = .0000776843638391 +mat T_V[6,3] = .0000284213133248 +mat T_V[6,4] = .0000382262676265 +mat T_V[6,5] = .0000332679318352 +mat T_V[6,6] = .0001398628868337 +mat T_V[6,7] = .0000615316221245 +mat T_V[6,8] = .0000437211453232 +mat T_V[6,9] = .0000656218010752 +mat T_V[7,1] = -.0000551650679429 +mat T_V[7,2] = .0002581489637783 +mat T_V[7,3] = -.0000824643289709 +mat T_V[7,4] = -.000042674305463 +mat T_V[7,5] = -.0000403565693949 +mat T_V[7,6] = .0000615316221245 +mat T_V[7,7] = .0002853895404404 +mat T_V[7,8] = .0003604010457427 +mat T_V[7,9] = .0003252736468058 +mat T_V[8,1] = -.0001468122320985 +mat T_V[8,2] = .0006493134095847 +mat T_V[8,3] = -.0002135535062099 +mat T_V[8,4] = -.0001311357514049 +mat T_V[8,5] = -.0000957474386806 +mat T_V[8,6] = .0000437211453232 +mat T_V[8,7] = .0003604010457427 +mat T_V[8,8] = .001327557632085 +mat T_V[8,9] = .0008655738151882 +mat T_V[9,1] = -.0001110067496262 +mat T_V[9,2] = .0006092564157568 +mat T_V[9,3] = -.0001716922597921 +mat T_V[9,4] = -.0000969870572202 +mat T_V[9,5] = -.0000643409318662 +mat T_V[9,6] = .0000656218010752 +mat T_V[9,7] = .0003252736468058 +mat T_V[9,8] = .0008655738151882 +mat T_V[9,9] = .0011805563999581 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + + +csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +assert `"`e(cmdline)'"' == `"csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .0169262437014923 +mat T_b[1,2] = -.069393526595745 +mat T_b[1,3] = .007932802034353 +mat T_b[1,4] = .0214173492404631 +mat T_b[1,5] = .0214285798296607 +mat T_b[1,6] = -.0237705768478331 +mat T_b[1,7] = -.0407911825853478 +mat T_b[1,8] = -.1153380429226659 +mat T_b[1,9] = -.0976743040271333 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0002311442272515 +mat T_V[1,2] = -.0000913516255354 +mat T_V[1,3] = .0002942556613254 +mat T_V[1,4] = .0002300102113643 +mat T_V[1,5] = .0001691668090648 +mat T_V[1,6] = 1.97724945317e-06 +mat T_V[1,7] = -.0000719224656945 +mat T_V[1,8] = -.0001704115735685 +mat T_V[1,9] = -.0001250497123318 +mat T_V[2,1] = -.0000913516255354 +mat T_V[2,2] = .0002956989691059 +mat T_V[2,3] = -.0001229611368996 +mat T_V[2,4] = -.0000827218925021 +mat T_V[2,5] = -.0000683718472045 +mat T_V[2,6] = .0000497459253988 +mat T_V[2,7] = .0001986659277167 +mat T_V[2,8] = .0005048986425726 +mat T_V[2,9] = .0004294853807356 +mat T_V[3,1] = .0002942556613254 +mat T_V[3,2] = -.0001229611368996 +mat T_V[3,3] = .0004321568114057 +mat T_V[3,4] = .0002650905107159 +mat T_V[3,5] = .0001855196618545 +mat T_V[3,6] = 4.29689321764e-06 +mat T_V[3,7] = -.0000922749994412 +mat T_V[3,8] = -.0002301133832139 +mat T_V[3,9] = -.0001737530581611 +mat T_V[4,1] = .0002300102113643 +mat T_V[4,2] = -.0000827218925021 +mat T_V[4,3] = .0002650905107159 +mat T_V[4,4] = .0002623370288506 +mat T_V[4,5] = .0001626030945263 +mat T_V[4,6] = 1.91446099020e-06 +mat T_V[4,7] = -.0000640406288387 +mat T_V[4,8] = -.0001547330524719 +mat T_V[4,9] = -.0001140283496882 +mat T_V[5,1] = .0001691668090648 +mat T_V[5,2] = -.0000683718472045 +mat T_V[5,3] = .0001855196618545 +mat T_V[5,4] = .0001626030945263 +mat T_V[5,5] = .0001593776708136 +mat T_V[5,6] = -2.79605848321e-07 +mat T_V[5,7] = -.0000594517688036 +mat T_V[5,8] = -.0001263882850197 +mat T_V[5,9] = -.0000873677291463 +mat T_V[6,1] = 1.97724945317e-06 +mat T_V[6,2] = .0000497459253988 +mat T_V[6,3] = 4.29689321764e-06 +mat T_V[6,4] = 1.91446099020e-06 +mat T_V[6,5] = -2.79605848321e-07 +mat T_V[6,6] = .0000730459300969 +mat T_V[6,7] = .0000481226655546 +mat T_V[6,8] = .000031455063973 +mat T_V[6,9] = .0000463600419708 +mat T_V[7,1] = -.0000719224656945 +mat T_V[7,2] = .0001986659277167 +mat T_V[7,3] = -.0000922749994412 +mat T_V[7,4] = -.0000640406288387 +mat T_V[7,5] = -.0000594517688036 +mat T_V[7,6] = .0000481226655546 +mat T_V[7,7] = .0002266976249852 +mat T_V[7,8] = .0002812880537518 +mat T_V[7,9] = .0002385553665754 +mat T_V[8,1] = -.0001704115735685 +mat T_V[8,2] = .0005048986425726 +mat T_V[8,3] = -.0002301133832139 +mat T_V[8,4] = -.0001547330524719 +mat T_V[8,5] = -.0001263882850197 +mat T_V[8,6] = .000031455063973 +mat T_V[8,7] = .0002812880537518 +mat T_V[8,8] = .0010145967593686 +mat T_V[8,9] = .0006922546931968 +mat T_V[9,1] = -.0001250497123318 +mat T_V[9,2] = .0004294853807356 +mat T_V[9,3] = -.0001737530581611 +mat T_V[9,4] = -.0001140283496882 +mat T_V[9,5] = -.0000873677291463 +mat T_V[9,6] = .0000463600419708 +mat T_V[9,7] = .0002385553665754 +mat T_V[9,8] = .0006922546931968 +mat T_V[9,9] = .0007407714211994 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + + +csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st) + + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(cluster)'"' == `"st"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +_assert_streq `"`e(cmdline)'"' `"csdid2 lemp , ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .017595643753745 +mat T_b[1,2] = -.0772398214573161 +mat T_b[1,3] = .003306356692512 +mat T_b[1,4] = .0250218295975542 +mat T_b[1,5] = .024458744971169 +mat T_b[1,6] = -.0199318167892598 +mat T_b[1,7] = -.050957367065195 +mat T_b[1,8] = -.1372587388894044 +mat T_b[1,9] = -.1008113630854053 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0007829071217662 +mat T_V[1,2] = -.0002717218539991 +mat T_V[1,3] = .0010315906382616 +mat T_V[1,4] = .000760194622637 +mat T_V[1,5] = .0005569361044 +mat T_V[1,6] = -.0001115421517223 +mat T_V[1,7] = -.0002578712924792 +mat T_V[1,8] = -.0003893609598156 +mat T_V[1,9] = -.0003281130119792 +mat T_V[2,1] = -.0002717218539991 +mat T_V[2,2] = .0002183504515259 +mat T_V[2,3] = -.0003294819691427 +mat T_V[2,4] = -.0002994560390609 +mat T_V[2,5] = -.0001862275537937 +mat T_V[2,6] = .0000568837146026 +mat T_V[2,7] = .0002265694813186 +mat T_V[2,8] = .0003051592816874 +mat T_V[2,9] = .0002847893284949 +mat T_V[3,1] = .0010315906382616 +mat T_V[3,2] = -.0003294819691427 +mat T_V[3,3] = .0015201009874528 +mat T_V[3,4] = .0009289254297727 +mat T_V[3,5] = .0006457454975592 +mat T_V[3,6] = -.000110003095552 +mat T_V[3,7] = -.0002160765790923 +mat T_V[3,8] = -.0005383288323113 +mat T_V[3,9] = -.0004535193696153 +mat T_V[4,1] = .000760194622637 +mat T_V[4,2] = -.0002994560390609 +mat T_V[4,3] = .0009289254297727 +mat T_V[4,4] = .0008118381067114 +mat T_V[4,5] = .0005398203314268 +mat T_V[4,6] = -.0001164242778234 +mat T_V[4,7] = -.0003698006047936 +mat T_V[4,8] = -.0003826496670879 +mat T_V[4,9] = -.0003289496065386 +mat T_V[5,1] = .0005569361044 +mat T_V[5,2] = -.0001862275537937 +mat T_V[5,3] = .0006457454975592 +mat T_V[5,4] = .0005398203314268 +mat T_V[5,5] = .0004852424842139 +mat T_V[5,6] = -.0001081990817914 +mat T_V[5,7] = -.0001877366935518 +mat T_V[5,8] = -.0002471043800477 +mat T_V[5,9] = -.0002018700597838 +mat T_V[6,1] = -.0001115421517223 +mat T_V[6,2] = .0000568837146026 +mat T_V[6,3] = -.000110003095552 +mat T_V[6,4] = -.0001164242778234 +mat T_V[6,5] = -.0001081990817914 +mat T_V[6,6] = .0001057850796644 +mat T_V[6,7] = .0000670490666338 +mat T_V[6,8] = .0000127123457306 +mat T_V[6,9] = .0000419883663815 +mat T_V[7,1] = -.0002578712924792 +mat T_V[7,2] = .0002265694813186 +mat T_V[7,3] = -.0002160765790923 +mat T_V[7,4] = -.0003698006047936 +mat T_V[7,5] = -.0001877366935518 +mat T_V[7,6] = .0000670490666338 +mat T_V[7,7] = .0004120548254679 +mat T_V[7,8] = .0002160765790923 +mat T_V[7,9] = .0002110974540805 +mat T_V[8,1] = -.0003893609598156 +mat T_V[8,2] = .0003051592816874 +mat T_V[8,3] = -.0005383288323113 +mat T_V[8,4] = -.0003826496670879 +mat T_V[8,5] = -.0002471043800477 +mat T_V[8,6] = .0000127123457306 +mat T_V[8,7] = .0002160765790923 +mat T_V[8,8] = .0005383288323113 +mat T_V[8,9] = .0004535193696153 +mat T_V[9,1] = -.0003281130119792 +mat T_V[9,2] = .0002847893284949 +mat T_V[9,3] = -.0004535193696153 +mat T_V[9,4] = -.0003289496065386 +mat T_V[9,5] = -.0002018700597838 +mat T_V[9,6] = .0000419883663815 +mat T_V[9,7] = .0002110974540805 +mat T_V[9,8] = .0004535193696153 +mat T_V[9,9] = .0004325521239021 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V + +csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st) + + +assert `"`e(base)'"' == `"Base Universal"' +assert `"`e(cntrl)'"' == `"Never treated"' +assert `"`e(cluster)'"' == `"st"' +assert `"`e(method)'"' == `"reg"' +assert `"`e(estat_cmd)'"' == `"csdid2_estat"' +_assert_streq `"`e(cmdline)'"' `"csdid2 lemp [w=lpop], ivar( countyreal) tvar(year) gvar(first) agg(event) cluster(st)"' +assert `"`e(cmd)'"' == `"csdid2"' +assert `"`e(vcetype)'"' == `"Robust"' +assert `"`e(properties)'"' == `"b V"' + +qui { +mat T_b = J(1,9,0) +mat T_b[1,1] = .0169262437014923 +mat T_b[1,2] = -.069393526595745 +mat T_b[1,3] = .007932802034353 +mat T_b[1,4] = .0214173492404631 +mat T_b[1,5] = .0214285798296607 +mat T_b[1,6] = -.0237705768478331 +mat T_b[1,7] = -.0407911825853478 +mat T_b[1,8] = -.1153380429226659 +mat T_b[1,9] = -.0976743040271333 +} +matrix C_b = e(b) +assert mreldif( C_b , T_b ) < 1E-8 +_assert_streq `"`: rowfullnames C_b'"' `"y1"' +_assert_streq `"`: colfullnames C_b'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_b T_b + +qui { +mat T_V = J(9,9,0) +mat T_V[1,1] = .0008553362098125 +mat T_V[1,2] = -.0003729056186157 +mat T_V[1,3] = .0011344934983186 +mat T_V[1,4] = .0008256904623471 +mat T_V[1,5] = .0006058246687719 +mat T_V[1,6] = -.0001762809606567 +mat T_V[1,7] = -.0003472355983124 +mat T_V[1,8] = -.0005319225374714 +mat T_V[1,9] = -.0004361833780222 +mat T_V[2,1] = -.0003729056186157 +mat T_V[2,2] = .0002955924619583 +mat T_V[2,3] = -.0004470501707058 +mat T_V[2,4] = -.0004022834464612 +mat T_V[2,5] = -.0002693832386801 +mat T_V[2,6] = .0000797094956907 +mat T_V[2,7] = .0003083500186895 +mat T_V[2,8] = .0004160569766317 +mat T_V[2,9] = .0003782533568214 +mat T_V[3,1] = .0011344934983186 +mat T_V[3,2] = -.0004470501707058 +mat T_V[3,3] = .0016447127057925 +mat T_V[3,4] = .0010073691441048 +mat T_V[3,5] = .0007513986450585 +mat T_V[3,6] = -.0001769648458234 +mat T_V[3,7] = -.0002916471911943 +mat T_V[3,8] = -.0007284565174212 +mat T_V[3,9] = -.0005911321283841 +mat T_V[4,1] = .0008256904623471 +mat T_V[4,2] = -.0004022834464612 +mat T_V[4,3] = .0010073691441048 +mat T_V[4,4] = .0008715078927645 +mat T_V[4,5] = .0005981943501718 +mat T_V[4,6] = -.0001981435967092 +mat T_V[4,7] = -.0004765814799506 +mat T_V[4,8] = -.0005077733730396 +mat T_V[4,9] = -.0004266353361453 +mat T_V[5,1] = .0006058246687719 +mat T_V[5,2] = -.0002693832386801 +mat T_V[5,3] = .0007513986450585 +mat T_V[5,4] = .0005981943501718 +mat T_V[5,5] = .0004678810110854 +mat T_V[5,6] = -.0001537344394376 +mat T_V[5,7] = -.0002734781237922 +mat T_V[5,8] = -.0003595377219535 +mat T_V[5,9] = -.0002907826695372 +mat T_V[6,1] = -.0001762809606567 +mat T_V[6,2] = .0000797094956907 +mat T_V[6,3] = -.0001769648458234 +mat T_V[6,4] = -.0001981435967092 +mat T_V[6,5] = -.0001537344394376 +mat T_V[6,6] = .0001225240012724 +mat T_V[6,7] = .0001146947512824 +mat T_V[6,8] = .0000207157532791 +mat T_V[6,9] = .0000609034769289 +mat T_V[7,1] = -.0003472355983124 +mat T_V[7,2] = .0003083500186895 +mat T_V[7,3] = -.0002916471911943 +mat T_V[7,4] = -.0004765814799506 +mat T_V[7,5] = -.0002734781237922 +mat T_V[7,6] = .0001146947512824 +mat T_V[7,7] = .0005344335781815 +mat T_V[7,8] = .0002974894749616 +mat T_V[7,9] = .0002867822703325 +mat T_V[8,1] = -.0005319225374714 +mat T_V[8,2] = .0004160569766317 +mat T_V[8,3] = -.0007284565174212 +mat T_V[8,4] = -.0005077733730396 +mat T_V[8,5] = -.0003595377219535 +mat T_V[8,6] = .0000207157532791 +mat T_V[8,7] = .0002974894749616 +mat T_V[8,8] = .000743048976445 +mat T_V[8,9] = .0006029737018409 +mat T_V[9,1] = -.0004361833780222 +mat T_V[9,2] = .0003782533568214 +mat T_V[9,3] = -.0005911321283841 +mat T_V[9,4] = -.0004266353361453 +mat T_V[9,5] = -.0002907826695372 +mat T_V[9,6] = .0000609034769289 +mat T_V[9,7] = .0002867822703325 +mat T_V[9,8] = .0006029737018409 +mat T_V[9,9] = .0005623539781832 +} +matrix C_V = e(V) +assert mreldif( C_V , T_V ) < 1E-8 +_assert_streq `"`: rowfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +_assert_streq `"`: colfullnames C_V'"' `"Pre_avg Post_avg tm4 tm3 tm2 tp0 tp1 tp2 tp3"' +mat drop C_V T_V +