program define amstat 
version 2.1
	mac def _varlist "req ex max(1)"
	mac def _options "Quiet"
	parse "%_*"
	quietly sum %_varlist, detail
	#delimit ;
	mac def _nobs = _result(1);
	quietly {
		gen double Y = _result(14)*sqrt( 
			( (%_nobs+1)*(%_nobs+3)) /
			( 6*(%_nobs-2) )
			) in 1 ;
		gen double Beta2= 
			(
			3*(%_nobs*%_nobs + 27*%_nobs-70)*(%_nobs+1)*(%_nobs+3)
			) / (
			(%_nobs-2)*(%_nobs+5)*(%_nobs+7)*(%_nobs+9)
			) in 1;
		gen double W2 = -1 + sqrt(2*(Beta2-1)) in 1 ;
		gen double delta = 1/sqrt(log(sqrt(W))) in 1;
		gen double alpha = sqrt(2/(W2-1)) in 1;
		gen double Z1= delta*log( Y/alpha + sqrt((Y/alpha)^2+1) ) 
			in 1 ; 

		gen double Eb2=(3*(%_nobs-1)) / (%_nobs+1) in 1 ;
		gen double Vb2 = 
			(
			24*%_nobs*(%_nobs-2)*(%_nobs-3)
			) / (
			(%_nobs+1)^2 * (%_nobs+3)*(%_nobs+5)
			) in 1 ;

		gen double X = (_result(15)-Eb2)/sqrt(Vb2) in 1 ;
		gen double RBeta1 = 
			(
			(
			6*(%_nobs*%_nobs-5*%_nobs+2)
			) / ( 
			(%_nobs+7)*(%_nobs+9)
			)
			) * sqrt(
			(
			6*(%_nobs+3)*(%_nobs+5)
			) / (
			%_nobs*(%_nobs-2)*(%_nobs-3)
			)
			) in 1 ;
		gen double A = 6 + (8/RBeta1)*(
			2/RBeta1 + sqrt(1+4/(RBeta1*RBeta1))
			) in 1 ; 
		gen double Z2 = ( 
			(1-2/(9*A))
			- ( (1-2/A)/(1+X*sqrt(2/(A-4))) )^(1/3)
			)
			/ sqrt(2/(9*A)) in 1 ;
		#delimit cr
	}
	mac def _K2 = Z1*Z1 + Z2*Z2
	if "%_quiet"=="" {
		di "Chi2(2) = %_K2   p = " chiprob(2,%_K2) 
	}
	mac def amstat "%_K2"
	drop Y-Z2
		
/*
	mac def _s = _result(1)*((_result(14)^2/6) + (_result(15)-3)^2/24)
	if "%_quiet"=="" {
		di "Chi2(2) = %_s   p = " chiprob(2,%_s)
	}
	mac def simple "%_s"
*/
end