/*
Power calculations with a moving part.

Computes power for different sample size and effect size, and stores them in a
matrix whose columns are the effect size and the rows the sample size.
Each cell in each matrix stores the statistical power.

A random sample of n casefiles are to be divided into 'status-quo' and 
'shift casefiles'. During the span of a year this casefiles are going to be 
notified in 60 zones (by 60 notifiers in each arm). After a casefile is notified
it is archived, otherwise it is reasigned to the same notifier in the first arm
and to another notifier in the second arm. Each zone and notifier has its own
fixed effect of probability of notified/notifying. 




The followig model is considered
	
	Y = \sum_{i}\alpha_i + \beta Z + \epsilon				(1)

	
The variables denote the following
	
	\alpha_i - Zone FE
	Z - exogenous randomization (treatment variable)
	Y - dummy dependent variable
	

Within each cell `m' replications are run to estimate the power.

*/



clear all
set more off
global directorio C:\Users\xps-seira\Downloads
timer on 1

*Sample size
local nn "500(250)6000"
local length_nn=23
*Number of replications
local m=1000
*Baseline dependent variable
local p_baseline=0.30
*Standard deviation of notification
local sd=.125
*Standard deviation of fixed effect notifier
local sd_not=.175
*Standard deviation of fixed effect zone
local sd_zone=.25
*Probability of treatment
local p_treat=0.5
*Significance level
local signif=0.010
*Grid for effect 
local beta_eff "0.02(0.02)0.20"
local length_fs=10
local sd_eff=0.1

********************************************************************************
********************************************************************************
/*
local j=0
	forvalues mu_x = 0.02(0.02)0.20 {
		local j=`j'+1
		di `mu_x'
	}
di `j'	
*/

local np=`length_nn'*`length_fs'*`m'
di "--------"
di "Número de procesos:"
di `np'
di "--------"
		
		
local k=0
local pb=0


	*Define matrices in which power is stored 
	matrix pwr= J(`length_nn',`length_fs',0)
qui{	
	local j=0
	forvalues mu_x = `beta_eff' {
		local j=`j'+1

		local i=0
		forvalues n = `nn' {
			local i=`i'+1

			forvalues t=1/`m' {
				local pb=`pb'+1
				
				clear 
				 set obs `n'
				
				*Simulate random noise
				 gen epsilon=rnormal(0,`sd')
				
				*Simulate a normal dist from where effect is going to be drawn
				 gen yy=rnormal(`mu_x',`sd_eff')
				
				*Zone-notifier dummies
				 gen unif=uniform()
				 gen zone=.
				 gen zone_fe=.
				 gen not_1_fe=.
				 gen not_2_fe=.
				 forvalues zn=1/60 {
					replace zone=`zn' if inrange(unif, (`zn'-1)/60 , `zn'/60)
					*FE + 'uniform' noise
					local fe=rnormal(0,`sd_zone')
					replace zone_fe=`fe'+(-.05+(.1/60)*`zn') if zone==`zn'
					local fe=rnormal(0,`sd_not')
					replace not_1_fe=`fe'+(-.05+(.1/60)*`zn') if zone==`zn'
					local fe=rnormal(0,`sd_not')
					replace not_2_fe=`fe'+(-.05+(.1/60)*`zn') if zone==`zn'
					}		
				
				*Exogenous randomization
				 gen zz=(uniform()<=`p_treat')
				
				*Simulate dep var
				 gen y=.
					*Baseline dep var probability
				 replace y=(uniform()<=`p_baseline'+zone_fe+not_1_fe+epsilon) if zz==0
					*ATT
				 replace y=(uniform()<=`p_baseline'+zone_fe+not_2_fe+yy+epsilon) if zz==1

				*REG
				 reg y zz i.zone, r cluster(zone)
				 test _b[zz]=0
				matrix pwr[`i',`j']=pwr[`i',`j']+(`r(p)'<=`signif')
				
				
								*Progress bar
			noi {
				if `pb'==1 {
					di "Progress"
					di "--------"
					}
				if `pb'==floor(`np'/10) {
					di "10%"
					}
				if `pb'==floor(`np'*2/10) {
					di "20%"
					}
				if `pb'==floor(`np'*3/10) {
					di "30%"
					}
				if `pb'==floor(`np'*4/10) {
					di "40%"
					}
				if `pb'==floor(`np'*5/10) {
					di "50%"
					}
				if `pb'==floor(`np'*6/10) {
					di "60%"
					}
				if `pb'==floor(`np'*7/10) {
					di "70%"
					}
				if `pb'==floor(`np'*8/10) {
					di "80%"
					}
				if `pb'==floor(`np'*9/10) {
					di "90%"
					}
				if `pb'==floor(`np'-10) {
					di "100%"
					di "--------"
					di "        "
					}
				}
				}
				
			matrix pwr[`i',`j']=pwr[`i',`j']/`m'
			}
		}
	
	
	clear
	 svmat pwr
	 export delimited using "$directorio\pwr1.csv", replace novarnames
}
	


timer off 1
timer list


*Power calculation graph
gen sample_size=_n*250+250

twoway (line pwr1 sample_size, lwidth(medthick)) ///
	(line pwr2 sample_size, lwidth(medthick)) ///
	(line pwr3 sample_size, lwidth(medthick)) ///
	(line pwr4 sample_size, lwidth(medthick)) ///
	(line pwr5 sample_size, lwidth(medthick)) ///
	(line pwr6 sample_size, lwidth(medthick)) ///
	(line pwr7 sample_size, lwidth(medthick)) ///
	(line pwr8 sample_size, lwidth(medthick)) ///
	(line pwr9 sample_size, lwidth(medthick)) ///
	(line pwr10 sample_size, lwidth(medthick)) ///
		, graphregion(color(white)) xtitle("Sample size") ytitle("Statistical power") ///
		legend(order(1 "2%" 2 "4%" 3 "6%" 4 "8%" 5 "10%" ///
					6 "12%" 7 "14%" 8 "16%" 9 "18%" 10 "20%") rows(2))
graph export "$directorio\stat_pwr.pdf", replace