/*
	swilk.ado
	P Royston
*/
*! version 1.0.0  
*  updated with 15 August 1991 note
program define swilk
version 2.1
	mac def _varlist "req ex min(1)"
	mac def _if "opt"
	mac def _in "opt"
	mac def _options "Quiet Lnnormal"
	/*
		Default options are: 
			display the output, 
			lnnormal (not 3-param lognormal),
			b1 (minimum skewness estimate of shift parameter, not 
				gammastar).
	*/
	parse "%_*"
	if "%_quiet"=="" { 
		if "%_lnnorma" == "" {
			di _n in gr _col(20) /* 
			*/ "Shapiro-Wilk W test for normal data"
		}
		else { 
			di _n in gr _col(10) /* 
			*/ "Shapiro-Wilk W test for 3-parameter lognormal data"
		}
		#delimit ; 
		di in gr " Variable |    Obs" _col(30) "W"
			_col(40) "V"
			_col(51) "z    Pr>z" _n 
		" ---------+" _dup(49) "-"
		;
		#delimit cr
	}
	capture {
		parse "%_varlist", parse(" ")
		while "%_1"~="" { 
			capture drop _WTEST
			capture drop _UVAR
			gen _UVAR=%_1 %_if %_in
			/*
				Calculate number of good values and SSQ of data
			*/
			summarize _UVAR
			if _result(1) < 7 { 
				di in red "insufficient observations"
				exit 2001
			}
			mac def _G=_result(1)
			mac def _SD=sqrt(_result(4))
			genrank _WTEST = _UVAR		/* calc. ranks */
			/*
				Calculate coefficients
			*/
			replace _WTEST=invnorm((_WTEST-0.375)/(%_G+0.25))
			quietly summarize _WTEST
			mac def _U=_result(4)*(%_G-1)
			mac def _R=4*(%_U-2*_WTEST[1]^2)
			mac def _U=sqrt(%_U)
			mac def _M=%_G-(%_G<=20)
			mac def _Q=exp(log(6*%_M+7)-log(6*%_M+13)+ /*
				*/ .5*(1+(%_M-2)*log(%_M+1)-(%_M-1)*log(%_M+2)))
			mac def _T=%_R/(1/%_Q-2)
			mac def _R=sqrt(%_R+2*%_T)
			replace _WTEST=2*_WTEST/%_R
			replace _WTEST=sqrt(%_T)/%_R in %_G
			replace _WTEST=-_WTEST[%_G] in 1
			/*
				Get W
			*/
			correlate _UVAR _WTEST
			mac def _W=_result(4)^2
			if %_G<=20 {		/* we already know _G>=6	*/
				mac def _X=log(%_G)-3
				mac def _L=118.898E-3+%_X*(133.414E-3+%_X*327.907E-3)
				mac def _M=exp(-(375.42E-3+%_X*(492.145E-3+ /* 
					*/ %_X*(1124.332E-3+%_X*199.422E-3))))
				mac def _S=exp(-3.15805+%_X*(729.399E-3+ /* 
					*/ %_X*(3.01855+%_X*1558.776E-3)))
			}
			else { 
				mac def _X=log(%_G)-5
				mac def _L=480.385E-3+%_X*(318.828E-3+ /*
					*/ %_X*%_X*(-24.1665E-3+ /*
					*/ %_X*(879.701E-5+%_X*298.9646E-5)))
				mac def _M=exp(-(1914.87E-3+%_X*(1378.88E-3+ /* 
					*/ %_X*(4183.209E-5+%_X*(-106.6339E-3+ /* 
					*/ %_X*(35.13666E-3+%_X*15.04614E-3))))))
				mac def _S=exp(-(3735.38E-3+%_X*(1015.807E-3+ /* 
					*/ %_X*(331.885E-3+%_X*(-177.3538E-3+ /* 
					*/ %_X*(1638.782E-5+%_X*(32.15018E-3+ /* 
					*/ %_X*(-385.2646E-5))))))))
			}
			if "%_lnnorma" != "" {
				mac def _X=log(%_G)
				if %_SD<0.75 {
					mac def _ZBAR=%_X*(-0.8625+ /*
						*/ %_X*(0.3285+ /* 
						*/ %_X*(-0.04292+0.0019*%_X)))
					mac def _ZSD=0.6615+ /* 
						*/ %_X*(0.0905-0.006395*%_X)
				}
				else {
					mac def _ZBAR=%_X*(-0.8955+ /*
						*/ %_X*(0.3552+ /* 
						*/ %_X*(-0.04820+0.002216*%_X)))
					mac def _ZSD=0.8161+ /* 
						*/ %_X*%_X*(0.01653- /*
						*/ 0.001643*%_X)
				}
				mac def _M=%_M+%_ZBAR*%_S
				mac def _S=%_S*%_ZSD
			}
			mac def _Z=((1-%_W)^%_L-%_M)/%_S
			mac def _P=normprob(-%_Z)
			/*
				V is (1-_W)/(median of 1-_W).
			*/
			mac def _V=(1-%_W)/(%_M^(1/%_L))
			if "%_quiet" == "" {
				mac def _skip = 9 - length("%_1") 
				#delimit ;
				di in gr _skip(%_skip) "%_1 |" in ye
					%7.0f %_G "    "
					%8.5f %_W "  "
					%8.3f %_V "   "
					%8.3f %_Z "  "
					%7.5f %_P ;
				;
				#delimit cr
			}
			mac shift
		}
		mac def S_1 %_G			/* # obs		*/
		mac def S_2 %_W			/* W statistic		*/
		mac def S_3 %_V			/* V			*/
		mac def S_4 %_Z			/* z			*/
		mac def S_5 %_P			/* P			*/
		drop _UVAR _WTEST
		exit
	}
	mac def _rc = _rc 
	capture drop _UVAR
	capture drop _WTEST
	exit %_rc
end