/* D'Agostino's kurtosis test for normality (by Anscombe&Glynn's approximation(1983)) */
/* Based on C.Minotani(1992) pp.84-85 and checked by 'dagoTest' in library 'fBasics' of R. */                   
new; cls;
x=rndn(100,1);
call kurtosistest(x);



proc kurtosistest(x);
   local n,xbar,s,kur,varb2,Eb2,Y,k3,A,Z,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   if n<=20;
      errorlog "Warning: Sample size must be greater than 20.";
      retp(".");
   endif;
   xbar=meanc(x);
   s=sqrt(sumc((x-xbar)^2)/n);     /* It should be divided by n to replicate most software. */
   kur=sumc(((x-xbar)/s)^4)/n;
   varb2=24*(n-2)*(n-3)*n/((n+1)^2*(n+3)*(n+5));
   Eb2=3*(n-1)/(n+1);
   Y=(kur-Eb2)/sqrt(varb2);
   k3=(6*(n^2-5*n+2)/((n+7)*(n+9)))*sqrt((6*(n+3)*(n+5))/(n*(n-2)*(n-3)));
   A=6+(8/k3)*((2/k3)+sqrt(1+(2/k3)^2));
   Z=((1-2/(9*A))-((1-2/A)/(1+Y*sqrt(2/(A-4))))^(1/3))/sqrt(2/(9*A));
   /* Z4 = (1-2/(9*A)-pos)/jm */
   pval=2*(1-cdfn(abs(Z)));
   print "stat=" Z;
   print "pval=" pval;
   retp(Z);
endp;