*! version 1.0.0  //  Ariel Linden 29mar2019 

program define retrodesign, rclass
version 11.0

        syntax anything, Se(numlist max=1) [ alpha(real 0.05) DF(string) SEED(string) Reps(integer 10000)  ]

		numlist "`anything'", min(1)
		tokenize `anything', parse(" ")
		local kn : list sizeof anything
		
		// error checking of std err
		if `se' <0 {
			di as err "standard error cannot be negative"
			exit 411
		}
		
		// Generate the matrix shell that will hold the results
		matrix X = J(`kn',3,0)
		
		// Loop over values of effect sizes
		forvalues i = 1/`kn' { 
			local A : word `i' of `anything'

			*************************************************
			* Gelman and Carlin’s (2014) simulation approach
			*************************************************
			if "`df'" != "" {
				if `df' > 9007199254740990 {
					di as err "df cannot be greater than 9007199254740990"
					exit 198
				}
			
				// set the seed
				if "`seed'" != "" {
					`version' set seed `seed'
				}
				local z = invt(`df',(1-`alpha'/2))
				local p_hi = 1 - t(`df', `z'- `A'/`se')
				local p_lo = t(`df', -`z'- `A'/`se')
				local power = `p_hi' + `p_lo'
				local typeS = cond(`A' >= 0, `p_lo'/`power', 1-(`p_lo'/`power'))

				tempvar sim estimate significant exaggeration
				quietly set obs `reps'
				gen `sim' = rt(`df')
				gen `estimate' = `A' + `se' * `sim'
				gen `significant' = abs(`estimate') > `se'*`z'
				gen `exaggeration' = abs(abs(`estimate')/`A')
				sum `exaggeration' if `significant'==1, meanonly
				local typeM = r(mean)
			} // end if "df"
		
			***********************************************
			* The closed form solution by Lu et al. (2018)
			***********************************************
			else {
				local z = invnorm(1-`alpha'/2)
				local p_hi = 1 - normal(`z'-`A'/`se')
				local p_lo = normal(-`z'-`A'/`se')
				local power = `p_hi' + `p_lo'
				local typeS = cond(`A' >= 0, `p_lo'/`power', 1-(`p_lo'/`power'))
				local lambda = `A'/`se'
				local typeM = (normalden(`lambda' + `z') + normalden(`lambda' - `z') + ///
					`lambda' * (normal(`lambda' + `z') + normal(`lambda' - `z') - 1)) / /// 
						(`lambda' * (1 - normal(`lambda' + `z') + normal(`lambda' - `z')))
			} // end "closed solution"

			// populate matrix with results after each loop of effect sizes
			matrix X[`i',1]= `power'
			matrix X[`i',2]= `typeS'
			matrix X[`i',3]= `typeM'

		} // end forvalues

		// format matrix table
		matrix colnames X = Power S-error M-error
		matrix rownames X = `anything' 
		matlist X, tw(12) lines(eq) border(bottom) showcoleq(comb) format(%9.4f) rowtitle(Effect)

		// table footnote
		di as txt
		di as txt "Note:"	
		di as txt "- Power is the probability that the statistical test correctly rejects the null hypothesis." 
		di as txt "- Type-S (sign) error is the probability of the sign being in the opposite direction of the effect size." 
		di as txt "- Type-M (magnitude) error is the factor by which the magnitude of the effect size might be exaggerated."

		// store results 
		mata : X = st_matrix("X")
		return matrix table = X

end