*! version 1.1.0  25feb2022
capture program drop srtfreq
program define srtfreq, rclass
version 15.1

syntax varlist(numeric fv) [if] [in], INTervention(varlist numeric fv max=1) [, NPerm(integer 0) NBoot(integer 0) SEED(integer 1020252) noDOT PERCentile BASIC noIsily PASTE *]
quietly {
        preserve
		if "`paste'"!="" {
        if "`nperm'" != "0" {
        cap drop PermC_I* PermUnc_I*
        }
        if "`nboot'" != "0" {
        cap drop BootC_I* BootUnc_I*
        }
		}
        if "`nboot'" == "0" & "`basic'"!="" {
		noi disp as error "Please specify number of bootstraps"
        error 198
		}
		if "`percentile'"!="" & "`basic'"!="" | "`nperm'" != "0" & "`nboot'" != "0" {
		noi disp as error "you have included resampling options that cannot be specified at the same time"
        error 198
		}

		local ci percentile
		if "`basic'"!="" local ci `basic'


        tempfile orig
        save `orig'

        local intervraw: copy local intervention

        fvrevar `intervention', list
        local intervention `r(varlist)'

        cap {
                local newvarlist
                foreach var of local varlist {
                        if "`var'" != "`intervention'" & !regexm("`var'","i[^\.]*\.`intervention'$") & !regexm("`var'","\((#[0-9])\)*\.`intervention'$") | regexm("`var'","#")   {
                                local newvarlist `newvarlist' `var'
                                }
                        else {
                        noi disp as txt "Note: Inclusion of the intervention variable in the variable list is redundant."
                        }
                        }
                local varlist: copy local newvarlist
                }

        fvrevar `varlist', list
        local varlist_clean `r(varlist)'

        marksample touse
        markout `touse' `intervention'
        keep if `touse'


        tempvar r brkn_fctor
        tempname test1 test0 Beta Beta1 Beta2 b Sigma1 Sigma2 Sigma x max nt nc N_total colnumber coef

        gettoken depvar indepvar: varlist

        tab `intervention'
        scalar `max' = r(r)

        baseset, max(`max') intervention(`intervraw')
        local refcat `r(refcat)'

        levelsof `intervention', local(levels)
        tokenize `levels'

        tempfile srt
        save `srt'

        tempfile srt2
        save `srt2'

        reg `depvar', `options'
        scalar `Sigma2' = e(rmse)

        `isily' reg `depvar' i.`intervention' `indepvar' , `options'

        scalar `Sigma1' = e(rmse)
        matrix `Beta' = e(b)
        matrix rownames `Beta' = "Beta"
        matrix `coef' = `Beta'[1,1..`=`max'']

        mata funczero("`coef'")

        matrix `coef'=`coef''

        matrix `Sigma' = round(`Sigma1'^2,.01),round(`Sigma2'^2,.01)
        matrix rownames `Sigma' = "Sigma2"
        matrix colnames `Sigma' = "Conditional" "Unconditional"

        scalar `colnumber' = colnumb(r(table),"`=`refcat'+0'b.`intervention'")

        if "`=`colnumber''"=="1" {
                matrix `test0' = r(table)[1...,2...]
                }
        else {
                matrix `test0' = r(table)[1...,1..`=`colnumber'-1'],r(table)[1...,`=`colnumber'+1'...]
                }
        matrix `test1' = `test0'
        forvalues i = 1/`=rowsof(`test0')' {
                forvalues j = 1/`=colsof(`test0')' {
                        matrix `test1'[`i',`j']= round(`test0'[`i',`j'],.01)
                }
        }

        matrix `Beta1' = (`test1'["b", .] \ `test1'["ll".."ul", . ])
        matrix `Beta2'=`Beta1''
        matrix colnames `Beta2' = "Estimate" "95% LB" "95% UB"


        forvalues i = 1/`=`max''{
                if "`=`refcat'+0'" != "``i''" {
                        local rowname `rowname' "`intervention'``i''"
                        local fnc `fnc' `i'
                        }
                else {
                local m `i'
                }
                }
        tokenize `fnc'
        qui tab `intervention', generate(`brkn_fctor') matcell(`x') /* extract matrix x of broken Intervention to use in calculations below*/
        drop `brkn_fctor'1
        clear


        forvalues s = 1/2 {
                forvalues j = 1/`=`max'-1' {
                tempname cd`j'`s' varcd`j'`s' secd`j'`s' cdlb`j'`s' cdub`j'`s' jdf g`j'`s' varg`j'`s' seg`j'`s'         glb`j'`s' gub`j'`s' A`j'`s' g`j'`s'
                        scalar `cd`j'`s'' = (`coef'[`j',1]/ `Sigma`s'' )
                        scalar `nt' = `x'[`m',1]
                        scalar `nc' = `x'[``j'',1]
                        scalar `varcd`j'`s''= ((`nt'+`nc')/(`nt'*`nc')+`cd`j'`s''^2/(2*(`nt'+`nc')))
                        scalar `secd`j'`s'' = sqrt(`varcd`j'`s'')
                        scalar `cdlb`j'`s'' = (`cd`j'`s'' - 1.96*`secd`j'`s'')
                        scalar `cdub`j'`s'' = (`cd`j'`s'' + 1.96*`secd`j'`s'')
                        scalar `jdf'     = (1 - (3/(4*(`nt'+`nc'-2)-1)))
                        scalar `g`j'`s''    = (`jdf'*`cd`j'`s'')
                        scalar `varg`j'`s'' = (`jdf'^2 * `varcd`j'`s'')
                        scalar `seg`j'`s''  = sqrt(`varg`j'`s'')
                        scalar `glb`j'`s''  = (`g`j'`s'' - 1.96*`seg`j'`s'')
                        scalar `gub`j'`s''  = (`g`j'`s'' + 1.96*`seg`j'`s'')
                        matrix `A`j'`s'' = round(`g`j'`s'',.01),round(`glb`j'`s'',.01),round(`gub`j'`s'',.01)
                        }
                }
                forvalues s = 1/2 {
                                tempname G`s' P`s'
                                matrix `G`s'' = `A1`s''
                                matrix `P`s'' = round(`g1`s'',.01)
                                }

                forvalues s = 1/2 {
                        if `=`max'-1'>1 {
                                matrix `G`s'' = `A1`s''
                                foreach i of numlist 2/`=`max'-1'{ /*extract g's (coefficients of the long calculation to be used later)*/
                                        matrix `G`s'' = `G`s'' \ `A`i'`s''
                                        }
                                matrix `P`s'' = round(`g1`s'',.01)
                                foreach i of numlist 2/`=`max'-1'{
                                        matrix `P`s'' = (`P`s'' \ round(`g`i'`s'',.01))
                                        }
                                }
                        }
                matrix CondES = `G1'
                matrix rownames CondES = `rowname'
                matrix colnames CondES = "Estimate" "95% LB" "95% UB"
                matrix UncondES = `G2'
                matrix rownames UncondES = `rowname'
                matrix colnames UncondES = "Estimate" "95% LB" "95% UB"
                use `srt'

   //====================================================//
  //===================                =================//
 //==================  PERMUTATIONS  ==================//
//=================                ===================//
//====================================================//

        if "`nperm'" != "0" {
				tempvar shuffle
                noisily di as txt "  Running Permutations..."

                set seed `seed'

                if `nperm'<1000         {
                        display as error "error: nPerm must be greater than 1000"
                        error 7
                        }
                forvalues k = 1/`nperm' {
                if "`dot'" == "" {
                                if !mod(`k', 100) {
                                noi di _c "`k'"
                                }
                        else {
                                if !mod(`k', 10) {
                                        noi di _c "."
                                        }
                                }
                        }
                        keep `intervention'
                        gen `shuffle' = runiform()
                        mata funcsh("`intervention'","`shuffle'")

                        drop `shuffle'

                        tempfile permute
                        save `permute'

                        use `srt2'
                        merge 1:1 _n using `permute', update replace nogenerate
                        tempfile srt2
                        save `srt2'
                        gettoken depvar indepvar: varlist

                        reg `depvar', `options'
                        scalar `Sigma2' = e(rmse)

                        `isily' reg `depvar' i.`intervention' `indepvar' , `options'
                        scalar `Sigma1' = e(rmse)
                        matrix `b' = e(b)
                        capture drop `coef'
                        matrix `coef' = `b'[1,1..`=`max'']
                        mata funczero("`coef'")
                        matrix `coef'=`coef''

                        qui tab `intervention', generate(`brkn_fctor') matcell(`x')
                        drop `brkn_fctor'1
                        clear

                        forvalues s = 1/2 {
                                forvalues j = 1/`=`max'-1' {
                                tempname c`k'd`j'`s' var`k'cd`j'`s' g`k'`j'`s'
                                        scalar `c`k'd`j'`s'' = (`coef'[`j',1]/ `Sigma`s'' )
                                        scalar `nt' = `x'[`m',1]
                                        scalar `nc' = `x'[``j'',1]
                                        scalar `var`k'cd`j'`s''= ((`nt'+`nc')/(`nt'*`nc')+`c`k'd`j'`s''^2/(2*(`nt'+`nc')))
                                        scalar `jdf'     = (1 - (3/(4*(`nt'+`nc'-2)-1)))
                                        scalar `g`k'`j'`s''    = (`jdf'*`c`k'd`j'`s'')
                                        }
                                }
                        use `srt2'
                        }
                        clear
                        use `srt'
                        count
                        scalar `N_total' = `r(N)'
                        if `nperm'>`=`N_total'' {
                                set obs `nperm'
                                }

                forvalues k = 1/`nperm' {
                        forvalues s = 1/2 {
                                forvalues i = 1/`=`max'-1' {
                                        capture gen Treat_PEst`s'_`i' = .
                                        replace Treat_PEst`s'_`i' = `g`k'`i'`s'' in `k'
                                        }
                                }
                        }

				/* Permutation Test*/
					tempname pval

					matrix `pval'=J(2,`=`max'-1',.)

					matrix rownames `pval' = "Conditional ES" "Unconditional ES"
					matrix colnames `pval' = `rowname'

				forvalues s=1/2 {
					forvalues i=1/`=`max'-1' {
							tempvar p`s'_`i'

							 gen `p`s'_`i'' = abs(Treat_PEst`s'_`i')>=abs(`cd`i'`s'')
							 summarize `p`s'_`i'', meanonly
							 matrix `pval'[`s',`i'] =round(r(mean),.01)
						}
					}


			return matrix Pv = `pval'

                tokenize `levels'
				if "`paste'"!="" {
                local m
                forvalues i = 1/`=`max'' {
                        if "`=`refcat'+0'" != "``i''" {
                        local m = `m' + 1
                        rename (Treat_PEst1_`m' Treat_PEst2_`m') (PermC_I``i'' PermUnc_I``i'')
                        }
                }
                keep PermC_I* PermUnc_I*

                tempfile perMES
                save `perMES'

                use `orig'
                merge 1:1 _n using `perMES', nogenerate
                tempfile orig
                save `orig'
				}

                if "`dot'" == "" {
                                noi di as txt ""
                                }
                noi di as txt "  Permutations completed."
        } /*if nperm*/



   //====================================================//
  //====================            ====================//
 //==================== BOOTSTRAPS ====================//
//====================            ====================//
//===================================================//

                if "`nboot'" != "0"  {
                clear
                use `srt'
                count
                cap scalar `N_total' = `r(N)'

                noisily di as txt "  Running Bootstraps..."
                if `nboot'<1000 {
                        display as error "error: nBoot must be greater than 1000"
                        error 7
                        }

                keep `varlist_clean' `intervention'


                tempfile touseit
                save `touseit'

                set seed `seed'

                forvalues i = 1/`nboot' {
                if "`dot'" == "" {
                                if !mod(`i', 100) {
                                noi di _c "`i'"
                                }
                        else {
                                if !mod(`i', 10) {
                                        noi di _c "."
                                        }
                                }
                        }

                        bsample
                        keep `varlist_clean' `intervention'
                        gettoken depvar indepvar: varlist

                        reg `depvar', `options'
                        scalar `Sigma2' = e(rmse)

                        `isily' reg `depvar' i.`intervention' `indepvar', `options'
                        scalar `Sigma1' = e(rmse)

                        matrix `b' = e(b)
                        capture drop `coef'
                        matrix `coef' = `b'[1,1..`=`max'']
                        mata funczero("`coef'")
                        matrix `coef'=`coef''

                forvalues s = 1/2 {
                        forvalues j = 1/`=`max'-1' {
                        tempvar Treat`j'_RM`i'SE`s'
                                scalar `Treat`j'_RM`i'SE`s'' = (`coef'[`j',1]/ `Sigma`s'' )
                                }
                        }
                        clear
                        use `touseit'
                        }
                        if `nboot'>`=`N_total'' {
                                set obs `nboot'
                                }
                forvalues i = 1/`nboot' {
                        forvalues s = 1/2 {
                                forvalues j = 1/`=`max'-1' {
                                        cap gen Treat`j'_RMSE`s' = .
                                        replace Treat`j'_RMSE`s' = `Treat`j'_RM`i'SE`s'' in `i'
                                        }
                                }
                        }

                forvalues s=1/2 {
				if "`ci'"=="percentile" {
                        forvalues i=1/`=`max'-1' {
                                centile Treat`i'_RMSE`s', centile(2.5)
                                tempname x_`i'`s'  y_`i'`s'
                                scalar `x_`i'`s''=r(c_1)
                                centile Treat`i'_RMSE`s', centile(97.5)
                                scalar `y_`i'`s''=r(c_1)
                                }
								}
				if "`ci'"=="basic" {
                        forvalues i=1/`=`max'-1' {
                                centile Treat`i'_RMSE`s', centile(2.5)
                                tempname x_`i'`s'  y_`i'`s'
                                scalar `y_`i'`s''=2*`cd`i'`s''-r(c_1)
                                centile Treat`i'_RMSE`s', centile(97.5)
                                scalar `x_`i'`s''=2*`cd`i'`s''-r(c_1)
                                }
								}
                                tempname L1`s' L2`s' L`s'
                        matrix `L1`s''=round(`x_1`s'',.01)
                        matrix `L2`s''=round(`y_1`s'',.01)

                        if `=`max'-1'>1 {
                                foreach i of numlist 2/`=`max'-1' {
                                        matrix `L1`s''=`L1`s'' \ round(`x_`i'`s'',.01)
                                        }
                                foreach i of numlist 2/`=`max'-1' {
                                        matrix `L2`s''=`L2`s'' \ round(`y_`i'`s'',.01)
                                        }
                                }
                        matrix `L`s'' = (`P`s'',`L1`s'',`L2`s'')
                        matrix rownames `L`s'' = `rowname'
                        matrix colnames `L`s'' = "Estimate" "95% LB (BT)" "95% UB (BT)"
                        }
                tokenize `levels'
				if "`paste'"!="" {
                local m
                forvalues i = 1/`=`max'' {
                        if "`=`refcat'+0'" != "``i''" {
                        local m = `m' + 1
                        rename (Treat`m'_RMSE1 Treat`m'_RMSE2) (BootC_I``i'' BootUnc_I``i'')
                        }
                }

                keep BootC_I* BootUnc_I*
                tempfile origBoot
                save `origBoot'
                clear
                use `orig'
                merge 1:1 _n using `origBoot', nogenerate
                tempfile orig
                save `orig'
				}

                if "`dot'" == "" {
                                noi di as txt ""
                                }
                noi di as txt "  Bootstraps completed."
         matrix CondES = `L1'
         matrix UncondES = `L2'
                }/*if nboot*/
        cap noi {
                noisily {
                        return matrix Beta = `Beta2'
                        return matrix Sigma2 = `Sigma'

                        matrix list CondES
                        return matrix CondES = CondES

                        matrix list UncondES
                        return matrix UncondES = UncondES
                        }
                }
        clear
use `orig'
restore, not
}

end


capture program drop baseset
program define baseset, rclass
syntax, max(name) INTervention(varlist fv)

        if regexm("`intervention'", "bn\.") | regexm("`intervention'", "^i\(?([0-9] ?)+\)?\.") {
        noi disp as error "i(numlist) not allowed; you must specify a base level"
        error 198
        }

        local refcat
        if regexm("`intervention'", "([0-9]?)[ ]*\.") local refcat = regexs(1)
        local allow opt1

        if "`refcat'" == "" {
                if regexm("`intervention'", "\(\#*([0-9]?)\)[ ]*\.") local refcat = regexs(1)
                if "`refcat'"!="" local allow opt2
        }
        if "`refcat'" == "" {
                if regexm("`intervention'", "\(([a-zA-Z]+)*\)\.") local refcat = regexs(1)
                if "`refcat'"!="" local allow opt3
                }

        fvrevar `intervention', list
        local intervention `r(varlist)'

        levelsof `intervention', local(levels)
        tokenize `levels'

        tempname min Max
        scalar `min'=`1'
        scalar `Max' = ``=`max'''

        if "`allow'" == "opt1" {
                forvalues i=1/`=`max'' {
                        if "`refcat'" != "``i''" local s = `s'+1 /*checking cases were intervention is irregular (i.e. 1 4 9) and user has specified a number of baseline that is not 1,4 or 9 and not below 1 a
> nd not above 9*/
                        }
                if "`refcat'" != "" {
                        if "`refcat'">"`=`Max''" | "`refcat'"<"`=`min''" | "`s'" == "`=`max''" {
                        noi disp as error "{bf:Warning:} selected baseline level `refcat' is out of bounds; level `=`Max'' chosen instead."
                                }
                        }
                else {
                        local refcat = `=`min''
                }
                if "`s'" == "`=`max''" & "`refcat'" != "`=`min''" {
                        local refcat = `=`Max''
                        }

                if "`refcat'" != "" {
                fvset base `refcat' `intervention'
                }
        }
        else if "`allow'"=="opt2" {
                fvset base ``refcat'' `intervention'
                local refcat = ``refcat''
                }
        else if "`allow'"=="opt3" {
                fvset base `refcat' `intervention'
                if strpos("`refcat'","first") >0 {
                local refcat = `=`min''
                }
                if strpos("`refcat'","last") >0 {
                local refcat = `=`Max''
                }
                if strpos("`refcat'","freq")>0 {
                tempname maximum z
                tab `intervention', matcell(`maximum')
                mata: st_local("matr", strofreal(max(st_matrix("`maximum'"))))
                forvalues i = 1/`=`max''{
                scalar `z' = `maximum'[`i',1]
                if "`matr'"== "`=`z''" local refcat = ``i''
                                }
                        }
                }
                return local refcat = `refcat'
                end