*!version 1.0.1 22/12/99	     (STB-56: sg149)
/*	Mark S. Pearce, University of Newcastle upon Tyne			*/
/* 	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		*/ 

program define seast
version 6.0
local options "EDwards EXAct NOTab Generate(string) SECtor(string) LENgth(string)"
local varlist "req ex min(1) max(2)"
parse "`*'"
parse "`varlist'", parse(" ")
local obs="`1'"
if "`2'"=="" {
	tempvar pop
	local nopop=1	
	gen `pop'=1
}
else {
	local pop ="`2'"
	local nopop=0
}

if "`generate'"!="" {
	confirm new var `generate'
}

if "`sector'"=="" {
	local sector="month"
}

tempvar theta weight xmean ymean srtm cosbit sinbit sumcos sumsin nmi totb totc sumnmi ci days x ex1 ey1 ex2 ey2 sumci exp ti sumti oe2 gfit suml rate rank 

quietly {

summ `sector'
local k=r(max)

egen `totb'=sum(`pop')
summ `totb'
local M=r(max)

egen `totc'=sum(`obs')
summ `totc'
local N=r(max)

local pi=3.1415927

if "`length'"!="" {
	egen `suml'=sum(`length')
	summ `suml'
	gen `x'=(360/`suml')*`length'*(1/57.296)
	gen `theta'=`pi'/`k'
	local stheta=`pi'/`k'
	local i=2
	while `i'<=`k' {
		summ `x' if `sector'==`i'-1
		replace `theta'=`stheta'+r(max) if `sector'==`i'
		local stheta=`stheta'+r(max)
		local i=`i'+1
	}
}	

else {
	if "`exact'"=="" {
		gen `theta' =`pi'/`k'*((2*`sector')-1)
	}
	else {
		if `k'==12 {
			gen `days'=28+2/9 if `sector'==2
			for num 1 3 5 7 8 10 12: replace `days'=31 if `sector'==X
			replace `days'=30 if `days'==.
			gen `x'=(360/365.25)*`days'*(1/57.296)
			gen `theta'=`pi'/12
			local stheta=`pi'/12
			local i=2
			while `i'<=12 {
				summ `x' if `sector'==`i'-1
				replace `theta'=`stheta'+r(max) if `sector'==`i'
				local stheta=`stheta'+r(max)
				local i=`i'+1
			}
		}
		else {
			if "`length'"=="" {
		di in red "Can only use exact with monthly data or by using the length option"
			e 100
			}
		}
	}	
}
egen `weight'=sum(sqrt(`obs'))
summ `weight'
local W=r(max)

*Calcluate observed centre point
egen `xmean'=sum((sqrt(`obs'))*(cos(`theta')))	
summ `xmean'
local xbar=r(max)/`W'

egen `ymean'=sum((sqrt(`obs'))*(sin(`theta')))
summ `ymean'
local ybar=r(max)/`W'

*Getting expected centre point

if "`edwards'"=="" {
	egen `srtm' = sum(sqrt(`pop'))
	summ `srtm'
	local srmi=r(max)
	gen `ex1'= (sqrt(`pop'))*(cos(`theta'))
	egen `ex2'=sum(`ex1')
	gen `ey1' =(sqrt(`pop'))*(sin(`theta'))
	egen `ey2'=sum(`ey1')
	summ `ex2'
	local xhat=r(max)/`srmi'
	summ `ey2'
	local yhat=r(max)/`srmi'
}
else {
	local xhat=0
	local yhat =0
}

*Variances
gen `cosbit'=0.25*(cos(`theta'))*(cos(`theta'))
gen `sinbit'=0.25*(sin(`theta'))*(sin(`theta'))
egen `sumcos'=sum(`cosbit')
summ `sumcos'
local vx1=r(max)
egen `sumsin'=sum(`sinbit')
summ `sumsin'
local vy1=r(max)

if "`edwards'"=="" {
	gen `nmi'=sqrt(`N'*`pop'/`M')
	egen `sumnmi'=sum(`nmi')
	summ `sumnmi'
	local v2=r(max)*r(max)

	local varx= `vx1'/`v2'
	local vary=`vy1'/`v2'
}

local dsq=((`xbar'-`xhat')^2)+((`ybar'-`yhat')^2)
local d=sqrt(`dsq')

local thstar=atan((`ybar'-`yhat')/(`xbar'-`xhat'))

local alpha = 4*`d'


if "`edwards'"!="" {
	gen `ci'=(1+`alpha'*cos(`theta'-`thstar'))
	egen `sumci'=sum(`ci')
	summ `sumci'
	gen `exp'=`N'*`ci'/r(max)
	local test = 0.5*`alpha'*`alpha'*`N'
}
else {
	gen `ci'=`pop'*(1+`alpha'*cos(`theta'-`thstar'))
	egen `sumci'=sum(`ci')
	summ `sumci'
	gen `exp'=`N'*`ci'/r(max)
	local test = (((`xbar'-`xhat')/sqrt(`varx'))^2) + (((`ybar'-`yhat')/sqrt(`vary'))^2)
}

gen `oe2'=((`obs'-`exp')^2)/`exp'
egen `gfit'=sum(`oe2')
summ `gfit'
local fit=r(max)

} *end of quietly

local thstar = `thstar'*57.296

if "`notab'"=="" {
	di in green _dup(50) "-"
	if `k'==12 {
		di in green "Month" _col(20) "|" _col(22) "Observed" _col(35) "Expected" 
	}
	else {
		di in green "Sector" _col(20) "|" _col(22) "Observed" _col(35) "Expected" 
	}
	di in green _dup(19) "-" _col(20) "+" _dup(30) "-"
	local i=1
	while `i'<=`k' {
		quietly summ `obs' if `sector'==`i'
		local obsi=r(max)
		quietly summ `exp' if `sector'==`i'
		local expi=r(max)
		if `k'!=12 {
			di in gr `i' _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
		}
		else {
			if `i'==1 {
				di in gr "January" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==2 {
				di in gr "February" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==3 {
				di in gr "March" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==4 {
				di in gr "April" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==5 {
				di in gr "May" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==6 {
				di in gr "June" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==7 {
				di in gr "July" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==8 {
				di in gr "August" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==9 {
			          di in gr "September" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==10 {
				di in gr "October" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==11 {
			            di in gr "November" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
			else if `i'==12 {
			          di in gr "December" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi'
			}
		}
		local i=`i'+1
	}
	di in green _dup(50) "-"
}
di _n
if "`edwards'"!="" {
	di in gr "Edwards Test" 
	di _n
}
else {
	if "`exact'"!="" {
		di in green "Walter & Elwood Test (using exact month lengths)" 
	}
	else {
		di in gr "Walter & Elwood Test"
	}
	di _n
}

quietly {
	tempvar sumobs
	gen `sumobs'=sum(`obs')
	summ `sumobs'
	local totobs=r(max)
}

di in green "Total Number of Cases =  " in yel `totobs'
di _n
di in green "Seasonality Test" _col(50) "Goodness of fit test"

di in green "chi2(" in yel "2" in gr ")" _col(15) "= " in yel %9.4f `test' _col(50) in gr "chi2(" in yel "2" in gr ")" _col(65) "=" in yel %9.4f `fit'
di in green "Prob > chi2" _col(15) "= " in yel  %9.4f chiprob(2,`test') _col(50) in gr "Prob > chi2" _col(65) "=" in yel %9.4f chiprob(`k'-1,`fit')
di _n
di in green _dup(50) "-"
di in green "Parameter" _col(38) "|"_col(40) "Estimate"
di in green _dup(37) "-" _col(38) "+" _dup(12) "-"
di in green "Amplitude of cyclic variation" _col(38) "|" _col(40) in yellow `alpha'
di in green  "Angle of maximum rate" _col(38) "|" _col(40) in yellow `thstar'
di in green _dup(50) "-"

if "`generate'"!="" {
	gen `generate'=`exp'
}

end