* . type rifle.ado *! rifle v1.0 CBerry AFowler 13July2017 version 14 program rifle, rclass version 14 set more off set matsize 5000 *! varlist contains outcome leader period unit syntax varlist(min=4 max=4) [if] [, Controls(varlist numeric fv) Permnum(integer 100) Report(integer 10) SAVING(string) REPLACE noSTITCH noGRAPH GROWTH noTREND] tokenize `varlist' local permutations = `permnum' local reportevery = `report' local outcome `1' local leader `2' local period `3' local unit `4' preserve * Convert outcome to percentage change if requested if "`growth'" == "growth" { sort `unit' `period' tempvar growth quietly gen `growth' = (`outcome' - `outcome'[_n-1])/`outcome'[_n-1] if `unit' == `unit'[_n-1] & `period' == `period'[_n-1] + 1 local outcome `growth' } if "`trend'" != "notrend" { * Residualize outcome by period quietly reg `outcome' i.`period' `controls' quietly predict resgrowth, residual local outcome resgrowth } * Create unit and leader IDs egen unitid = group(`unit') egen leaderid = group(`leader') local unit unitid local leader leaderid * Drop if missing on any of the 4 required variables quietly drop if `outcome'==. | `unit'==. | `leader'==. | `period'==. if "`stitch'" != "nostitch" { *stitch together blocks of continuguous years *and generate new time variable that starts at 1 in each unit egen unit_time = group(`unit' `period') egen min = min(unit_time), by(`unit') gen stitchtime = unit_time - min + 1 local period stitchtime } if "`stitch'" == "nostitch" { *drop cases where we only have 1 leader for a country-period quietly { egen leader_id = group(`unit' `leader') sort `unit' `period' gen unitperiod = _n replace unitperiod = unitperiod[_n-1] if `unit' == `unit'[_n-1] & `period' == `period'[_n-1] + 1 egen sdleaderid = sd(leader_id), by(unitperiod) drop if sdleaderid == 0 drop sdleaderid leader_id egen upid = group(unitperiod) drop unitperiod local unit upid egen min_period = min(`period'), by(`unit') g period2 = `period' - min_period + 1 local period period2 drop min_period } } * Run initial regression on real data quietly { tempname memhold tempfile Placebos postfile `memhold' fstat r2 adjr2 real using `Placebos', replace areg `outcome', absorb(`leader') post `memhold' (e(F_absorb)) (e(r2)) (e(r2_a)) (1) } *generate identifier for leader stints sort `unit' `period' gen stint = _n quietly replace stint = stint[_n-1] if `leader' == `leader'[_n-1] & `unit' == `unit'[_n-1] egen stintid = group(stint) drop stint rename stintid stint * Do the permutations quietly sum `unit' local nunits = r(max) disp "Permutations Completed:" forvalues k = 1/`permutations' { quietly { *randomly permute leader stints as blocks gen random = uniform() sort `unit' `period' replace random = random[_n-1] if stint == stint[_n-1] egen stint2 = group(`unit' random) egen minstint = min(stint2), by(`unit') replace stint2 = stint2 - minstint + 1 drop minstint gen leader2 = . forvalues j = 1/`nunits' { scalar counter = 0 sum stint2 if `unit' == `j' local nstints = r(max) forvalues i = 1/`nstints' { sum `leader' if stint2 == `i' & `unit' == `j' replace leader2 = r(mean) if `period' > counter & `period' <= counter + r(N) & `unit' == `j' scalar counter = counter + r(N) } } areg `outcome', absorb(leader2) post `memhold' (e(F_absorb)) (e(r2)) (e(r2_a)) (0) drop leader2 stint2 random } if `k'/`reportevery' == round(`k'/`reportevery') { disp `k' * disp "$S_TIME $S_DATE" } } postclose `memhold' clear use `Placebos' if ("`saving'"!="") { save `saving', `replace' } *--------------------- Compute and display p-values quietly { foreach i in fstat r2 adjr2 { gen `i'_temp = `i' if real == 1 egen `i'_real = max(`i'_temp) /* Why is this step necessary */ gen `i'_bigger = `i' > `i'_real replace `i'_bigger = . if real==1 /* Ensure real data ommitted when pvalue compted below */ } } quietly { sum r2_real scalar realr2 = r(mean) local realr2 = r(mean) sum r2 scalar avg = r(mean) scalar diff = realr2 - avg sum r2_bigger scalar pval = r(mean) local pval = r(mean) } display as text "--------------------------------------------------" display as text "R2 in real data: " as result %-5.4g realr2 display as text "Avg R2 in permuted data: " as result %-5.4g avg display as text "Difference: " as result %-5.4f diff display as text "p-value: " as result %-5.4g pval display as text "--------------------------------------------------" return scalar real_r2 = realr2 return scalar avg_r2 = avg return scalar diff = diff return scalar pval = pval *--------------------- Make graph if "`graph'" != "nograph" { quietly { sum r2 local min = floor(r(min)*100)/100 local max = ceil(r(max)*100)/100 hist r2 if real == 0, xline(`realr2', lpattern(dash) lwidth(thick)) bin(100) fraction title("RIFLE: Distribution Under the Null") /// xtitle("R-squared from Permuted Data") /// caption("Dashed line shows R-squared from the real data." "p-value = `pval'", tstyle(smbody) box just(center)) gr_edit .xaxis1.reset_rule `min' `max' .01 , tickset(major) ruletype(range) } } end