/*
	Shapiro-Francia W' test
	P Royston
*/
program define sfrancia
version 2.1
	mac def _varlist "req ex"
	mac def _if "opt"
	mac def _in "opt"
	mac def _options "quiet"
	parse "%_*"
	parse "%_varlist", parse(" ")
	if "%_quiet"=="" { 
		#delimit ; 
		di in gr _n 
			_col(16) "Shapiro-Francia W' test for normal data" _n
			" Variable |    Obs" _col(29) "W'"
			_col(39) "V'"
			_col(51) "z    Pr>z" _n 
		" ---------+" _dup(49) "-"
		;
		#delimit cr
	}
	capture {
		while "%_1"~="" { 
			/*
				Calculate ranks
			*/
			capture drop _WD1
			genrank _WD1 = %_1 %_if %_in
			/*
				Calculate number of good values
			*/
			summarize _WD1
			if _result(1)<5 { 
				di in red "insufficient observations"
				exit 2001
			}
			mac def _G=_result(1)
			replace _WD1 = invnorm((_WD1-0.375)/(%_G+0.25))
			corr %_1 _WD1
			mac def _X=log(%_G)-5
			mac def _L=-0.0480157+%_X*(0.01971964-0.0119065*%_X*%_X)
			mac def _M=-exp(1.6930674+%_X*(0.1441647+%_X*(-0.01849276  /*
			*/        +%_X*(0.031074485+%_X*0.0055717663))))
			mac def _S=exp(-0.510725+%_X*(-0.1160364+%_X*(-0.006702098 /*
			*/        +%_X*(0.054465944+%_X*0.0087397329))))
			mac def _Y=(((1-(_result(4))^2)^%_L)-1)/%_L
			mac def _Z=(%_Y-%_M)/%_S
			mac def _R=(_result(4))^2
			mac def _w=(1-%_R)/( (%_L*%_M+1)^(1/%_L) )
			mac def S_5 = normprob(-%_Z)
			if "%_quiet"=="" {
				mac def _skip = 9 - length("%_1") 
				#delimit ;
				di in gr _skip(%_skip) "%_1 |" in ye
					%7.0f %_G "    "
					%8.5f %_R "  "
					%8.3f %_w "   "
					%8.3f %_Z "  "
					%7.5f max(%S_5,.00001)
				;
				#delimit cr
			}
			mac shift 
		}
		drop _WD1
		mac def S_1 %_G		/* # observations		*/
		mac def S_2 %_R		/* W'				*/
		mac def S_3 %_w		/* V'				*/
		mac def S_4 %_Z		/* z				*/
		/* mac def S_5 */	/* P				*/
		exit
	}
	mac def _rc = _rc 
	capture drop _WD1
	exit %_rc
end