*! 1.0.0 Ariel Linden 22Sep2023 


capture program drop power_cmd_concord
program define power_cmd_concord, rclass
version 12.0


			/* obtain settings */
			syntax anything(id="numlist"),	///
				[ Alpha(real 0.05)        	/// significance level
				n(string)					/// total sample size
				Power(string)				///
				ONESIDed					///
				]

				numlist "`anything'", min(2) max(2)
				
				gettoken concord0 rest : anything
				gettoken concord1 rest : rest
				
				if `concord0' < 0 | `concord0' > 1 {
					noi di in red "concord0 must be a number between 0 and 1"
					exit 198
				} 
				if `concord1' < 0 | `concord1' > 1 {
					noi di in red "concord1 must be a number between 0 and 1"
					exit 198
				} 
				
				if `concord0' == `concord1' {
					noi di in red "concord0 must be different than concord1"
					exit 198					
				}
				if `alpha' < 0 | `alpha' > 1 {
					di as err "option {bf:alpha()} must contain numbers between 0 and 1"
					exit 121
				}
				

				*******************
				// Sample size  //
				*******************
				if ("`n'" == "") {
					
					if `power' < 0 | `power' > 1 {
						di as err "option {bf:power()} must contain numbers between 0 and 1"
						exit 121
					} 
					tempname zalpha zbeta delta n z0 z1
				
					// set default test type to "twosided"  
					if "`onesided'" == "" {
						local test = `alpha'/ 2
					}
					else local test = `alpha'

					scalar `zalpha' = invnorm(1-`test')
					scalar `zbeta' = invnorm(`power')
					scalar `delta' = abs(`concord0' - `concord1')
					
					scalar `z0' = log((1 + `concord0') / (1 - `concord0')) / 2
					scalar `z1' = log((1 + `concord1') / (1 - `concord1')) / 2					

 					scalar `n' =  ceil((((`zalpha' + `zbeta') / (`z1' - `z0'))^2)+2) 
				
				} // end sample size

				
				*******************
				// Power  //
				*******************				
				else if ("`n'" !="") {	
					
					tempname zalpha delta power z0 z1
				
					// set default test type to "twosided"  
					if "`onesided'" == "" {
						local test = `alpha'/ 2
					}
					else local test = `alpha'

					scalar `zalpha' = invnorm(1-`test')
					scalar `delta' = abs(`concord0' - `concord1')
					scalar `z0' = log((1 + `concord0') / (1 - `concord0')) / 2
					scalar `z1' = log((1 + `concord1') / (1 - `concord1')) / 2	

					scalar `power' = normal(((`z1' - `z0') * sqrt(`n' - 2) - (`zalpha')))
							
				} // end power
			
				// saved results
				return scalar N = `n'
				return scalar alpha = `alpha'
				return scalar power = `power'
				return scalar concord0 = `concord0'
				return scalar concord1 = `concord1'
				return scalar delta = `delta'



end