/* Skewness test for normality */
/* Based on C.Minotani(1992) pp.83-84 and checked by 'dagoTest' in library 'fBasics' of R. */                   
new; cls;
x=rndn(100,1);
call skewtest(x);



proc skewtest(x);
   local n,xbar,s,sk,Y,gam,B2,delta,alpha,Z,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   if n<=8;
      errorlog "Warning: Sample size must be greater than 8.";
      retp(".");
   endif;
   xbar=meanc(x);
   s=sqrt(sumc((x-xbar)^2)/n);     /* It should be divided by n to replicate most software. */
   sk=sumc(((x-xbar)/s)^3)/n;
   Y=sk*sqrt((n+1)*(n+3)/(6*(n-2)));  
   gam=(3*(n^2+27*n-70)*(n+1)*(n+3))/((n-2)*(n+5)*(n+7)*(n+9));
   B2=-1+sqrt(2*(gam-1));
   delta=1/sqrt(ln(sqrt(B2)));
   alpha=sqrt(2/(B2-1));
   Z=delta*ln(Y/alpha+sqrt((Y/alpha)^2+1));
   pval=2*(1-cdfn(abs(Z)));
   print "stat=" Z;
   print "pval=" pval;
   retp(Z);
endp;