*! 2.0.0 17 May 2005     		(SJ 5-3: sg149_1; STB-56: sg149)

/*	
	tests for seasonality with a variable population at risk
	
 	follows the tests given in Walter & Elwood (1975), Edwards (1961):
 	default is Walter & Elwood test using equal month lengths
	when a population variable is specified and Edwards test otherwise
*/ 		 

program seast, rclass  
	version 8   
	syntax varlist(min=1 max=2) /// 
	[ , EXact NOTab Generate(string) SECtor(varname) LENgth(varname) ]

	if "`sector'" == "" { 
		capture confirm variable month 
		if _rc { 
			di as err "default month variable not found" 
			exit 498 
		}	
		local sector "month"
	} 	
	
	summ `sector', meanonly 
	local k = r(max)

	if "`exact'" != "" & "`length'" == ""  & `k' != 12 {
		di as err "can only use exact with monthly data or length()"
		exit 198 
	} 
	
	if "`generate'" != "" confirm new var `generate'

	tokenize `varlist'
	args obs pop 
	if "`pop'" == "" {
		local edwards 1 
		tempvar pop
		gen byte `pop' = 1
	}
	else local edwards 0 
	
	tempvar theta work days x exp     

	quietly {
		summ `pop', meanonly 
		local M = r(sum)

		summ `obs', meanonly 
		local N = r(sum)

		if "`length'" != "" {
			summ `length', meanonly 
			gen `x' = 2 * _pi / r(sum) * `length' 
			gen `theta' = _pi / `k'
			local stheta = _pi / `k'
			
			forval i = 2 / `k' {
				summ `x' if `sector' == `i' - 1, meanonly 
				replace `theta' = ///
				`stheta' + r(max) if `sector' == `i'
				local stheta = `stheta' + r(max)
			}
		}	
		else {
			if "`exact'" == ""  { 
				gen `theta' = _pi / `k' * ((2 * `sector') - 1)
			}
			else {
				// `k' == 12 
				local exacttxt " (using exact month lengths)"
				gen `days' = 31 if inlist(`sector',1,3,5,7,8,10,12)  
				replace `days' = 28.25 if `sector' == 2
				replace `days' = 30 if `days' == .
				gen `x' = 2 * _pi / 365.25 * `days'
				gen `theta' = _pi / 12
				local stheta = _pi / 12 
				
				forval i = 2 / 12 {
					summ `x' if `sector' == `i' - 1, meanonly 
					replace `theta' = `stheta' + r(max) if `sector' == `i'
					local stheta = `stheta' + r(max)
				}
			}
		}
		
		gen `work' = sum(sqrt(`obs'))
		local W = `work'[_N] 

		// Calculate observed centre point
		replace `work' = sum((sqrt(`obs')) * (cos(`theta')))	
		local xbar = `work'[_N] / `W'
		replace `work' = sum((sqrt(`obs')) * (sin(`theta')))
		local ybar = `work'[_N] / `W'

		// Getting expected centre point
		if !`edwards' {
			replace `work' = sum(sqrt(`pop'))
			local srmi = `work'[_N] 
			replace `work' = sum(sqrt(`pop') * cos(`theta')) 
			local xhat = `work'[_N] / `srmi'
			replace `work' = sum(sqrt(`pop') * sin(`theta')) 
			local yhat = `work'[_N] / `srmi'
		}
		else {
			local xhat = 0
			local yhat = 0
		}

		// Variances
		replace `work' = 0.25 * cos(`theta') * cos(`theta')
		sum `work', meanonly 
		local vx1 = r(sum)
		replace `work' = 0.25 * sin(`theta') * sin(`theta')
		summ `work', meanonly 
		local vy1 = r(sum)
		
		if !`edwards' {
			replace `work' = sqrt(`N' * `pop' / `M')
			summ `work', meanonly 
			local v2 = r(sum) * r(sum)
			local varx = `vx1' / `v2'
			local vary = `vy1' / `v2'
		}
		
		local alpha = 4 * sqrt((`xbar' - `xhat')^2 + (`ybar' - `yhat')^2)
		local thstar = atan((`ybar' - `yhat') / (`xbar' - `xhat'))
		if `xbar' - `xhat' < 0 local thstar = _pi + `thstar'
		
		if `edwards' {
			replace `work' = (1 + `alpha' * cos(`theta' - `thstar'))
			summ `work', meanonly
			gen `exp' = `N' * `work' / r(sum)
			local test = 0.5 * `alpha' * `alpha' * `N'
		}
		else {
			replace `work' = `pop' * (1 + `alpha' * cos(`theta' - `thstar'))
			summ `work', meanonly 
			gen `exp' = `N' * `work' / r(sum)
			local test = ///
			(((`xbar' - `xhat') / sqrt(`varx'))^2) + (((`ybar' - `yhat') / sqrt(`vary'))^2)
		}

		replace `work' = ((`obs' - `exp')^2)/`exp'
		summ `work', meanonly 
		local fit = r(sum)

		local thstar = `thstar' * 180 / _pi 
		local foy = cond(`thstar' < 0, 1 + `thstar'/360, `thstar'/360) 
	} // end quietly 
	
	if "`notab'" == "" {
		di _n as txt "{hline 14}{c TT}{hline 23}"
		local which = cond(`k' == 12, "Month", "Sector") 
		di as txt "`which'" "{col 15}{c |} Observed" "{col 31}Expected" 
		di as txt "{hline 14}{c +}{hline 23}"

		tokenize "`c(Months)'" 
		forval i = 1 / `k' {
			summ `obs' if `sector' == `i', meanonly 
			local obsi = r(max)
			summ `exp' if `sector' == `i', meanonly 
			local expi = r(max)
			if `k' != 12 di as txt `i' "{col 15}{c |}"     ///
			                as res %9.0f `obsi'  ///
					as res "{col 30}" %9.3f `expi'
			else         di as txt "``i''" "{col 15}{c |}" ///
					as res %9.0f `obsi'    ///
					as res "{col 30}" %9.3f `expi'
		}
		di as txt "{hline 14}{c BT}{hline 23}"
	}
	
	if `edwards' di _n as txt "Edwards test" 
	else         di as txt _n "Walter & Elwood test`exacttxt'" 

	di _n as txt "Total number of cases =  " as res `N'
	di _n as txt "Seasonality test" "{col 35}Goodness of fit test"
	local chitxt1 `"as txt "chi2(" as res "2" as txt ")""' 
	local chitxt2 `"as txt "chi2(" as res "`= `k' -1'" as txt ")""' 
	di `chitxt1' "{col 15}= " as res %9.4f `test' ///
	   "{col 35}" `chitxt2' "{col 50}=" as res %9.4f `fit'	
	di as txt "Prob > chi2" "{col 15}= " as res %9.4f chiprob(2,`test') ///
	   as txt "{col 35}Prob > chi2" "{col 50}=" as res %9.4f chiprob(`k' - 1,`fit')
	   
	di _n as txt "{hline 30}{c TT}{hline 9}"
	di as txt "Parameter" "{col 31}{c |} Estimate"
	di as txt "{hline 30}{c +}{hline 9}" 
	di as txt "Amplitude of cyclic variation" "{col 31}{c |}" ///
	as res %9.3f `alpha'
	di as txt "Angle of maximum rate" "{col 31}{c |}" ///
	as res %9.1f `thstar'
	di as txt "Fraction of year" "{col 31}{c |}" ///
	as res %9.3f `foy'

	di as txt "{hline 30}{c BT}{hline 9}"

	return scalar foy = `foy' 
	return scalar angmax = `thstar' 
	return scalar amplitude = `alpha' 
	return scalar chi2_seasonal = `fit' 
	return scalar chi2_fit = `test' 
	return scalar n_cases = `N' 

	if "`generate'" != "" gen `generate' = `exp'
end