/* Jarque Bera statistic */
/* Based on C.Minotani(1992) p.34 and checked by 'jarqueberaTest' in library 'fBasics' of R. */   
new; cls;
x=rndn(100,1);
call JB(x);



proc JB(x);
   local n,xbar,s,z,sk,kur,stat,pval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
/* skewness & kurtosis */
   xbar=meanc(x);
   s=sqrt(sumc((x-xbar)^2)/n);     /* It should be divided by n to replicate most software. */
   z=(x-xbar)/s;
   sk=sumc((z)^3)/n;
   kur=sumc((z)^4)/n;
/* main */
   stat=n*sk^2/6+n*(kur-3)^2/24;
   pval=cdfchic(stat,2);
/* results */
   print "Jarque-Bera Test:";
   print/rz "n   =" n;
   print    "stat=" stat;
   print    "pval=" pval;
   retp(stat);
endp;