/* Corrected Jarque Bera statistic */
new; cls;
x=rndn(100,1);
call CJB(x);



proc CJB(x);
   local n,xbar,s,z,sk,kur,stat,pval;
   if ismiss(x);
      errorlog "Warninig: 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=sk^2*(n+1)*(n+3)/(6*(n-2))+(n+1)^2*(n+3)*(n+5)*(kur-3*(n-1)/(n+1))^2/(24*n*(n-2)*(n-3));
   pval=cdfchic(stat,2);
/* results */
   print "Corrected Jarque-Bera Test:";
   print/rz "n   =" n;
   print    "stat=" stat;
   print    "pval=" pval;
   retp(stat);
endp;