*! didplacebo 1.0.1 2023.08.19
prog def didplacebo, eclass
	version 16
    qui xtset
    if "`r(panelvar)'" == "" | "`r(timevar)'" == "" {
		di as err "panel variable or time variable missing, please use -{bf:xtset} {it:panelvar} {it:timevar}"
		exit 198
    }
	syntax anything, TREATVar(varname) [RANUnitnum(numlist min=1 max=1 >0 integer) RANTimescope(numlist min=2 max=2 >0 integer) ///
		REPeat(integer 500) seed(integer 1) PBOUnit PBOTime(numlist >0 integer) PBOMix(numlist min=1 max=3 >=1 <=3 integer) frame(string) noFIGure SAVEGraph(string) nodrop]
	loc panelvar "`r(panelvar)'"
    loc timevar "`r(timevar)'"
	local graphlist = ""
**# Check estimation
	cap qui estimates restore `anything'
	if _rc {
		di as err "{bf:`anything'} should be an estimation generated by {bf:estimates store}."
		exit 198
	}
	else {
		qui estimates restore `anything'
		loc estimatename "`anything'"
		loc cmd "`e(cmd)'"
		loc cmdline "`e(cmdline)'"
	}
**# Check treatvar()
	cap confirm numeric var `treatvar'
	if _rc {
		di as err "{bf:`treatvar'} should be a numeric variable in dataset"
		exit 198
	}
**# Check frame()
	if "`frame'" == "" tempname frame
	else {
		cap frame drop `frame'
		qui pwf
		if "`frame'" == "`r(currentframe)'" {
			di as err "invalid frame() -- current frame can not be specified"
			exit 198
		}
		loc framename "`frame'"
	}
	qui frame create `frame'
	qui frame `frame': set obs `repeat'
**# Creat treatvartip
	if "`cmd'" == "xtdidregress" | "`cmd'" == "didregress" {
		qui estimates restore `anything'
		mata: temp = st_matrixcolstripe("e(b)")[1,.]; st_local("treatvartip", temp[1,1] + ":" + temp[1,2]);
	}
	else loc treatvartip "`treatvar'"
**# Store treatment effect
	loc treff = _b["`treatvartip'"]
	tempname pbounitp pbotimeci pbomixp1 pbomixp2 pbomixp3 b V
	mata: data = st_data(., "`panelvar' `timevar' `treatvar'");
	mata: units = uniqrows(data[., 1])';
**# In-time Placebo Test
	if "`pbotime'" != "" {
		preserve
		di
	    di as txt "Implementing in-time placebo test using fake treatment time forwarded by " _continue
		mata: didplacebo_intimePBOTEST(data, "`panelvar'", "`timevar'", "`treatvar'", "`treatvartip'", "`cmdline'", strtoreal(tokens("`pbotime'")), "`drop'")
		di as txt " periods respectively."
		di
		ereturn post b V
		di as txt "Results of in-time placebo test using fake treatment times:" _continue
		di
		ereturn display
		di `"{p 0 6 2}{txt}Note: The standard errors are computed using the same method as specified by the Stata command previously used for DID estimation. For example, if "{bf:xtreg, r}" or "{bf:reghdfe, cluster(clustvar)}" is used, then cluster-robust standard errors are reported.{p_end}"'
		tempname pboframe
		qui frame create `pboframe'
		frame `pboframe'{
			matrix `pbotimeci' = (r(table))'
			qui svmat `pbotimeci', names(col)
			mata: st_store(., st_addvar("int", "forwarded"), abs(strtoreal(tokens("`pbotime'")')))
			twoway(rcap ul ll forwarded)(connected b forwarded, msymbol(smcircle_hollow)), xlabel(`pbotime') ///
			yline(0, lp(dash)) xtitle("number of periods forwarded as fake treatment time") ytitle("placebo effect") title("In-time Placebo Test") ///
			legend(order(2 "Placebo Effect" 1 "95% Confidence Interval") position(6) rows(1)) nodraw name(pbotime, replace)
			local graphlist = "`graphlist' pbotime"
			mkmat forwarded b se z pvalue ll ul, matrix(`pbotimeci')
		}
		restore
	}
**# In-space Placebo Test
	if "`pbounit'" != "" {
		preserve
		di
		di as txt "Implementing in-space placebo test using fake treatment units:"
		mata: CoefOfPboUnit = didplacebo_inspacePBOTEST(data, "`panelvar'", "`timevar'", "`treatvar'", "`treatvartip'", units,  "`cmdline'", strtoreal("`repeat'"), strtoreal("`seed'"))
		frame `frame': mata: st_store(., st_addvar("double", "CoefOfPboUnit"), CoefOfPboUnit)
		mata: didplacebo_pbomatp(strtoreal(st_local("treff")), CoefOfPboUnit, "`pbounitp'", "`treatvar'");
		di
		di as txt "Results of in-space placebo test results using fake treatment units:" _continue
		matlist `pbounitp', border(rows) showcoleq(c) cspec(& %10s | %11.6f | %10.4f o1 & %10.4f & %11.4f &)  rspec(---)
		di "{p 0 6 2}{txt}Note: (1) The two-sided p-value is the frequency that the absolute values of the placebo effects are greater than or equal to the absolute value of estimated treatment effect.{p_end}"
		di "{p 6 6 2}{txt}(2) The left-sided (right-sided) p-value is the frequency that the placebo effects are smaller (greater) than or equal to the estimated treatment effect.{p_end}"
		mata: st_local("xmax", strofreal(colmaxabs(CoefOfPboUnit \ strtoreal(st_local("treff")))))
		frame `frame': graph twoway (kdensity CoefOfPboUnit) ///
			(histogram CoefOfPboUnit, fcolor(gs8%50) lcolor(white) lalign(center) below), ///
			xline(0,lp(dash)) xline(`treff', lp(solid)) xtitle("distribution of placebo effects") ytitle("density") xsc(r(-`xmax' `xmax')) ///
			title("In-space Placebo Test") legend(order(1 "Kernel density estimate" 2 "Histogram") position(6) rows(1)) name(pbounit, replace) nodraw
		local graphlist = "`graphlist' pbounit"
		restore
	}
**# Mixed Placebo Test
	if "`pbomix'" != "" {
		foreach i of numlist `pbomix'{
			preserve
			di
			mata: printf("{txt}Implementing" + ("`i'"!="1"? ("`i'"=="2" ? " unrestricted " : " restricted ") : " ") + "mixed placebo test for" + ("`i'"!="1"? " staggered ": " standard ") + "DID (version {res}" + "`i'" + "{txt}) using both fake treatment units and times:\n")
			if "`ranunitnum'" == "" {
				mata: units_tr = uniqrows(select(data[., 1], data[., 3]))';
				mata: ranunitnum = cols(units_tr);
			}
			else mata: ranunitnum = strtoreal(st_local("ranunitnum"));
			if "`rantimescope'" == "" {
				mata: time = uniqrows(data[., 2])';
				mata: data_tmp = select(data, data[.,3]);
				mata: _sort(data_tmp, (1,2,3));
				mata: info = panelsetup(data_tmp, 1, 2);
				mata: data_sum = data_tmp[info[.,1],1..2];
				mata: timelb = ("`i'"=="1"? sort(time', 1)[2, .]:min(data_sum[., 2]));
				mata: timeub = ("`i'"=="1"? max(time):max(data_sum[., 2]));
			}
			else mata: tmp = tokens("`rantimescope'"); timelb = strtoreal(tmp[1]); timeub = strtoreal(tmp[2]);
			mata: printf(("`i'"!="3"?"{hline 34}{c TT}{hline 38}\n": "{hline 38}\n"))
			mata: printf(("`i'"!="3"?"{txt}   The number of units randomly   {c |}     The range within which fake      \n": "{txt}     The range within which fake      \n"))
			mata: printf(("`i'"!="3"?"{txt} selected as fake treatment units {c |} treatment times are randomly selected\n":"{txt} treatment times are randomly selected\n"))
			mata: printf(("`i'"!="3"?"{hline 34}{c +}{hline 38}\n":"{hline 38}\n"))
			mata: if("`i'"!="3") printf("%18.0f{space 16}{c |}", ranunitnum); printf("%~38uds\n", ("[" + strofreal(timelb) + ", " + strofreal(timeub) + "]"));
			mata: printf(("`i'"!="3"?"{hline 34}{c BT}{hline 38}\n":"{hline 38}\n"))
			mata: times_ltd = uniqrows(select(data[., 2], (data[., 2]:>=timelb):&(data[.,2]:<=timeub)))';
			mata: CoefOfPboMix`i' = didplacebo_mixPBOTEST(data, "`panelvar'", "`timevar'", "`treatvar'", "`treatvartip'", ranunitnum, units, times_ltd, "`cmdline'", strtoreal("`repeat'"), strtoreal("`seed'"), strtoreal("`i'"))
			frame `frame': mata: st_store(., st_addvar("double", "CoefOfPboMix`i'"), CoefOfPboMix`i')
			mata: st_local("pbomixp", (("`i'" == "1") ? "`pbomixp1'" : (("`i'" == "2") ? "`pbomixp2'" : "`pbomixp3'")))
			mata: didplacebo_pbomatp(strtoreal(st_local("treff")), CoefOfPboMix`i', "`pbomixp'", "`treatvar'");
			di
			mata: printf("{txt}Results of" + ("`i'"!="1"? ("`i'"=="2" ? " unrestricted " : " restricted ") : " ") + "mixed placebo test for" + ("`i'"!="1"? " staggered ": " standard ") + "DID (version {res}" + "`i'" + "{txt}) using both fake treatment units and times:")
			matlist `pbomixp', border(rows) showcoleq(c) cspec(& %10s | %11.6f | %10.4f o1 & %10.4f & %11.4f &)  rspec(---)
			di "{p 0 6 2}{txt}Note: (1) The two-sided p-value is the frequency that the absolute values of the placebo effects are greater than or equal to the absolute value of estimated treatment effect.{p_end}"
			di "{p 6 6 2}{txt}(2) The left-sided (right-sided) p-value is the frequency that the placebo effects are smaller (greater) than or equal to the estimated treatment effect.{p_end}"
			mata: st_local("xmax", strofreal(colmaxabs(CoefOfPboMix`i'\strtoreal(st_local("treff")))))
			mata: st_local("gphtitle", ("`i'"!="1"? ("`i'"=="2" ? "Unrestricted " : "Restricted ") : "") + "Mixed Placebo Test for" + ("`i'"!="1"? " Staggered ": " Standard ") + "DID")
			frame `frame': graph twoway (kdensity CoefOfPboMix`i') ///
				(histogram CoefOfPboMix`i', fcolor(gs8%50) lcolor(white) lalign(center) below), ///
				xline(0, lp(dash)) xline(`treff', lp(solid)) xtitle("distribution of placebo effects") ytitle("density") xsc(r(-`xmax' `xmax')) ///
				legend(order(1 "Kernel density estimate" 2 "Histogram") position(6) rows(1)) name(pbomix`i', replace) nodraw title("`gphtitle'")
			local graphlist = "`graphlist' pbomix`i'"
			restore
		}
	}
	cap local graphlist = strltrim("`graphlist'")
**# Display and save graphs
	if "`figure'" != "nofigure" {
	    if "`savegraph'" != "" {
			di
			ereturn local graphlist "`graphlist'"
			didplacebo_SaveGraph `savegraph'
		}
		else foreach graph in `graphlist'{
			capture graph display `graph'
		}
    }
	ereturn clear
	cap ereturn matrix pbomix3 = `pbomixp3'
	cap ereturn matrix pbomix2 = `pbomixp2'
	cap ereturn matrix pbomix1 = `pbomixp1'
	cap ereturn matrix pbotime = `pbotimeci'
	cap ereturn matrix pbounit = `pbounitp'
	if "`pbounit'" != "" | "`pbomix'" != "" ereturn loc frame "`framename'"
	cap ereturn loc graph "`graphlist'"
	ereturn loc seed "`seed'"
	ereturn loc cmdline "didplacebo `0'"
	ereturn loc cmd "didplacebo"
	ereturn loc treatvar "`treatvar'"
	ereturn loc timevar "`timevar'"
    ereturn loc panelvar "`panelvar'"
	ereturn loc didcmdline "`cmdline'"
	ereturn loc didcmd "`cmd'"
	ereturn loc didestname "`estimatename'"
	di
	di as txt "Finished."
end

program didplacebo_SaveGraph
	version 16
	preserve
	syntax [anything], [asis replace]
	foreach graph in `e(graphlist)'{
		graph save `graph' `anything'_`graph', `asis' `replace'
	}
end

version 16
mata:
	void didplacebo_intimePBOTEST(real matrix data, string scalar panelvar, string scalar timevar, string scalar treatvar, string scalar treatvartip, string scalar cmdline, real matrix forwards, string scalar nodrop){
		real matrix data_sum_F, data_tmp, info, data_sum, data_res, key, pboeffs_b, pboeffs_V; real scalar t; string matrix names; real scalar time_max, time_tr_max, i;
		time_max = max(data[., 2]');
		data_tmp = select(data, data[.,3]);
		_sort(data_tmp, (1,2,3));
		info = panelsetup(data_tmp, 1, 2);
		data_sum = data_tmp[info[.,1],1..2];
		if(nodrop == ""){
			time_tr_max = max(data_sum[., 2]')
			stata(sprintf("qui drop if " + timevar + " >= %g", time_tr_max))
		}
		stata("qui drop if " + treatvar + " == 1")
		pboeffs_b = J(1, cols(forwards), 0)
		pboeffs_V = J(cols(forwards), cols(forwards), 0)
		names = J(cols(forwards), 2, "")
		data_res = st_data(., panelvar + " "+ timevar + " " + treatvar);
		for(t = 1; t<=cols(forwards); t++){
			data_sum_F = (data_sum[., 1], data_sum[., 2]:-forwards[., t])
			key = asarray_create("real");
			asarray_notfound(key, time_max + 1);
			for(i = 1; i <= rows(data_sum_F); i ++) asarray(key, data_sum_F[i, 1], data_sum_F[i, 2]);
			for(i = 1; i <= rows(data_res); i++) data_res[i, 3] = (data_res[i, 2] >= asarray(key, data_res[i, 1]) ? 1: 0);
			st_store(., treatvar,  data_res[, 3]);
			// stata("panelview " + treatvar + ", i(" +panelvar+ ") t(" + timevar + ")  type(treat) name(GRAPH" + strofreal(t) + ", replace)")
			stata(cmdline, 1);
			stata(sprintf(`"matrix tmpb = e(b)[., "%s"]"', treatvartip));
			stata(sprintf(`"matrix tmpV = e(V)["%s", "%s"]"', treatvartip, treatvartip));
			pboeffs_b[1, t] = st_matrix("tmpb");
			pboeffs_V[t, t] = st_matrix("tmpV");
			names[t, 2] = sprintf("F%g.%s", forwards[., t], treatvar)
			printf("{res}%g", forwards[t]);
			if(t != cols(forwards)) printf(", ");
			displayflush();
		}
		st_matrix("b", pboeffs_b)
		st_matrix("V", pboeffs_V)
		st_matrixcolstripe("b", names);
		st_matrixrowstripe("V", names);
		st_matrixcolstripe("V", names);
	}
	real matrix didplacebo_inspacePBOTEST(real matrix data, string scalar panelvar, string scalar timevar, string scalar treatvar, string scalar treatvartip, real matrix units, string scalar cmdline, real scalar reps, real scalar seed){
		real scalar time_max, seed_org, i, r; real matrix units_random, data_tmp, data_sum, info, key, data_sum_random, pboeffs, data_res;
 		seed_org = rngstate();
 		rseed(seed);
		time_max = max(data[., 2]');
		data_tmp = select(data, data[.,3]);
		_sort(data_tmp, (1,2,3));
		info = panelsetup(data_tmp, 1, 2);
		data_sum = data_tmp[info[.,1],1..2];
		data_res = data;
		pboeffs = J(reps, 1, .)
		printf("{txt}Simulations ({res}" + strofreal(reps) + "{txt}):")
		for(r = 1; r<=reps; r++){
			data_sum_random = data_sum;
			key = asarray_create("real");
			units_random = jumble(units')';
			for(i = 1; i <= cols(units); i++) asarray(key, units[1, i], units_random[1, i]);
			for(i = 1; i <= rows(data_sum_random); i ++) data_sum_random[i, 1] = asarray(key, data_sum_random[i, 1]);
			key = asarray_create("real");
			asarray_notfound(key, time_max + 1);
			for(i = 1; i <= rows(data_sum_random); i ++) asarray(key, data_sum_random[i, 1], data_sum_random[i, 2]);
			for(i = 1; i <= rows(data_res); i++) data_res[i, 3] = (data_res[i, 2] >= asarray(key, data_res[i, 1]) ? 1: 0)
			st_store(., treatvar,  data_res[,3])
			// stata("panelview " + treatvar + ", i(" +panelvar+ ") t(" + timevar + ")  type(treat) name(GRAPH" + strofreal(r) + ", replace)")
			stata(cmdline, 1)
			stata("scalar pboeff = _b["+treatvartip+"]")
			pboeffs[r, .] = st_numscalar("pboeff")
			if(mod(r, 10) == 0) printf("%g", r); else printf(".");
			displayflush();
		}
 		rngstate(seed_org);
		return(pboeffs);
	}
	real matrix didplacebo_mixPBOTEST(real matrix data, string scalar panelvar, string scalar timevar, string scalar treatvar, string scalar treatvartip, real unitnum_ltd, real matrix units, real matrix times, string scalar cmdline, real scalar reps, real scalar seed, real scalar method){
		real scalar time_max, seed_org, i, unitnum, r; real matrix data_tmp, data_sum, info, key, data_sum_random, pboeffs, units_random, times_random, times_tr, data_res;
		seed_org = rseed();
		rseed(seed);
		time_max = max(data[., 2]');
		data_tmp = select(data, data[.,3]);
		_sort(data_tmp, (1,2,3));
		info = panelsetup(data_tmp, 1);
		data_sum = data_tmp[info[.,1],1..2];
		times_tr = uniqrows(data_sum[., 2]);
		unitnum = rows(uniqrows(data_sum[., 1]));
		data_res = data;
		pboeffs = J(reps, 1, .);
		printf("{txt}Simulations ({res}" + strofreal(reps) + "{txt}):")
		for(r = 1; r<=reps; r++){
			if(method == 1){
				units_random = (jumble(units')')[.,1..unitnum_ltd];
				times_random = (jumble(times')')[.,1];
				key = asarray_create("real");
				asarray_notfound(key, time_max + 1);
				for(i = 1; i <= unitnum_ltd; i ++) asarray(key, units_random[1, i], times_random);
			}
			if(method == 2){
				units_random = (jumble(units')')[.,1..unitnum_ltd];
				times_random = times[., runiformint(1, unitnum_ltd, 1, cols(times))];
				key = asarray_create("real");
				asarray_notfound(key, time_max + 1);
				for(i = 1; i <= unitnum_ltd; i ++) asarray(key, units_random[1, i], times_random[1, i]);
			}
			if(method == 3){
				if(rows(times_tr)>cols(times)) {
					printf("{err}The range of fake treatment periods specified by the options rantimescope({it:t_min t_max}) must be greater than or equal to the number of periods in the actual treatment variable.\n")
					exit(198);
				}
				times_random = jumble(times')[1..rows(times_tr),.];
				key = asarray_create("real");
				for(i = 1; i <= rows(times_tr); i ++) asarray(key, times_tr[i, 1], times_random[i, 1]);
				data_sum = data_tmp[info[.,1], 1..2];
				for(i = 1; i <= rows(data_sum); i++) data_sum[i, 2] = asarray(key, data_sum[i, 2]);
				units_random = jumble(units')[1..unitnum, .];
				data_sum_random = (units_random, data_sum[., 2]);
				key = asarray_create("real");
				asarray_notfound(key, time_max + 1);
				for(i = 1; i <= rows(data_sum_random); i ++) asarray(key, data_sum_random[i, 1], data_sum_random[i, 2]);
			}
			for(i = 1; i <= rows(data_res); i++) data_res[i, 3] = (data_res[i, 2] >= asarray(key, data_res[i, 1]) ? 1: 0);
			st_store(., treatvar,  data_res[, 3])
			// stata("panelview " + treatvar + ", i(" +panelvar+ ") t(" + timevar + ")  type(treat) name(GRAPH" + strofreal(r) + ", replace)")
			stata(cmdline, 1)
			stata("scalar pboeff = _b["+treatvartip+"]")
			pboeffs[r, .] = st_numscalar("pboeff")
			if(mod(r, 10) == 0) printf("%g", r); else printf(".");
			displayflush();
		}
		rngstate(seed_org);
		return(pboeffs);
	}
	void didplacebo_pbomatp(real scalar treff, real matrix pboeffs, string scalar matrixname, string scalar treatvar){
		real matrix pbomatp;
		pbomatp= J(1,4,.)
		pbomatp[1,1] = treff
		pbomatp[1,2] = mean(abs(pboeffs) :>= abs(treff))
		pbomatp[1,3] = mean(pboeffs :<= treff)
		pbomatp[1,4] = mean(pboeffs :>= treff)
		st_matrix(matrixname, pbomatp)
		st_matrixrowstripe(matrixname, ("", treatvar))
		st_matrixcolstripe(matrixname, (("", "Coefficient")\("P-value", "Two-sided")\("P-value", "Left-sided")\("P-value", "right-sided")))
	}
end
* Version history
* 1.0.1 Fix the issue of excessive treated units causing the program to be unable to run
* 1.0.0 Submit the initial version