program define sktestdc
version 2.1
	mac def _varlist "req ex min(1)"
	mac def _if "opt"
	mac def _in "opt"
	mac def _exp "opt"
/*
	Start of P Royston modification 1
*/
	mac def _options "noAdjust"
	parse "%_*"
	parse "%_varlist", parse(" ")
	if "%_adjust" ~= "noadjust" {
		mac def _adj "adj"
	}
	else {
		mac def _adj "   "
	}
	di _n _col(17) in gr "Skewness/Kurtosis tests for Normality" _n /*
		*/ _col(47) "------- joint -------" _n /*
		*/ " Variable |  Pr(Skewness)   Pr(Kurtosis)" /*
		*/ "  %_adj chi-sq(2)  Pr(Chi-sq)" _n /*
		*/ "----------+" _dup(56) "-"
	while ("%_1"!="") {
		quietly sum %_1 %_if %_in %_exp, detail
		if _result(1) < 8 {
			di in red "less than 8 observations"
			exit 2001
		}
/*
	End of P Royston modification 1
*/
		#delimit ;
		mac def _nobs = _result(1);
		mac def _Y = _result(14)*sqrt( 
			( (%_nobs+1)*(%_nobs+3)) /
			( 6*(%_nobs-2) )
			) ;
		mac def _Beta2= 
			(
			3*(%_nobs*%_nobs + 27*%_nobs-70)*(%_nobs+1)*(%_nobs+3)
			) / (
			(%_nobs-2)*(%_nobs+5)*(%_nobs+7)*(%_nobs+9)
			) ;
		mac def _W2 = -1 + sqrt(2*(%_Beta2-1)) ;
		mac def _delta = 1/sqrt(log(sqrt(%_W2))) ;
		mac def _alpha = sqrt(2/(%_W2-1)) ;
		mac def _Z1= %_delta*log( 
			%_Y/%_alpha + sqrt((%_Y/%_alpha)^2+1) 
			) ;

		mac def _Eb2=(3*(%_nobs-1)) / (%_nobs+1) ;
		mac def _Vb2 = 
			(
			24*%_nobs*(%_nobs-2)*(%_nobs-3)
			) / (
			(%_nobs+1)^2 * (%_nobs+3)*(%_nobs+5)
			) ;

		mac def _X = (_result(15)-%_Eb2)/sqrt(%_Vb2) ;
		mac def _RBeta1 = 
			(
			(
			6*(%_nobs*%_nobs-5*%_nobs+2)
			) / ( 
			(%_nobs+7)*(%_nobs+9)
			)
			) * sqrt(
			(
			6*(%_nobs+3)*(%_nobs+5)
			) / (
			%_nobs*(%_nobs-2)*(%_nobs-3)
			)
			) ;
		mac def _A = 6 + (8/%_RBeta1)*(
			2/%_RBeta1 + sqrt(1+4/((%_RBeta1)*(%_RBeta1)))
			) ; 
		mac def _Z2 = ( 
			(1-2/(9*%_A))
			- ( (1-2/%_A)/(1+%_X*sqrt(2/(%_A-4))) )^(1/3)
			)
			/ sqrt(2/(9*%_A)) ;
		mac def _K2 = %_Z1*%_Z1 + %_Z2*%_Z2 ;
		/*
			Start of P Royston modification 2.

			Empirical adjustment to chi-square(2) statistic.
			_ZC2 is normal deviate corresponding to putative
			 chisq(2) _K2. This is adjusted to a final normal
			deviate _Z, from which the P-value for the test
			is calculated as its upper tail area.
			The original chi-sq is adjusted so that the reported
			P value is the same as would be obtained from tables.
			Also the P value for kurtosis is adjusted.
		*/
		#delimit cr
		if "%_adjust" ~= "noadjust" {
			mac def _ZC2 = -invnorm(exp(-0.5*%_K2))
			mac def _logn = log(%_nobs)
			mac def _cut = .55*(%_nobs^.2)-.21
			mac def _a1 = (-5+3.46*%_logn)*exp(-1.37*%_logn)
			mac def _b1 = 1+(.854-.148*%_logn)*exp(-.55*%_logn)
			mac def _b2mb1 = 2.13/(1-2.37*%_logn)
			mac def _a2 = %_a1-%_b2mb1*%_cut
			mac def _b2 = %_b2mb1+%_b1
			if %_ZC2<-1 {
				mac def _Z = %_ZC2
			}
			else if %_ZC2<%_cut { 
				mac def _Z = %_a1+%_b1*%_ZC2
			}
			else {
				mac def _Z=%_a2+%_b2*%_ZC2
			}


			mac def _P = 1-normprob(%_Z)
			mac def _K2 = -2*log(%_P)
		}
		else {
			mac def _P = chiprob(2,%_K2)
		}
		#delimit ;
		/*
			End of P Royston modification 2.
		*/
		mac def _indent = 9 - length("%_1") ;
		di in gr _skip(%_indent) "%_1 |" in ye  
			_col(18) %5.3f 2-2*normprob(abs(%_Z1))
			_col(32) %5.3f 2-2*normprob(abs(%_Z2))
			_col(44) %9.2f %_K2
			_col(60) %6.4f %_P ;	/* P Royston modification 3 */
		#delimit cr
		mac shift
	}
end