/* Jarque Bera statistic with critical value by simulation */  
/* Notice: Table in Deb & Sefton(1996) might be wrong. I verify the table in Cho & Im(2002).  */             
new; cls;
x=rndn(100,1);
times=10000;                      /* more than 200000 recommended(by normal version of GAUSS) */
pp=0.05;
call JB2(x,times,pp);

/* When meantimes=1, the same as above. As it is bigger(no limits), that will take more time. */
meantimes=10; 
call JB2a(x,times,meantimes,pp);




proc JB2(x,times,pp);
   local n,xbar,s,z,sk,kur,stat,statm,i,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   xbar=meanc(x);
   s=sqrt(sumc((x-xbar)^2)/n);
   z=(x-xbar)/s;
   sk=sumc((z)^3)/n;
   kur=sumc((z)^4)/n;
   stat=n*sk^2/6+n*(kur-3)^2/24;
/* simulation */
   statm=zeros(times,1);
   i=1;
   do while i<=times;
      x=rndn(n,1);
      x=(x-meanc(x))/stdc(x);
      xbar=meanc(x);
      s=sqrt(sumc((x-xbar)^2)/n);
      z=(x-xbar)/s;
      sk=sumc((z)^3)/n;
      kur=sumc((z)^4)/n;
      statm[i]=n*sk^2/6+n*(kur-3)^2/24;
      i=i+1;
   endo;
   statm=sortc(statm,1);
   cval=statm[round(times*(1-pp))];
/* results */
   print "Jarque-Bera Test:";
   print/rz "n             =" n;
   print    "stat          =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%";
   retp(stat~cval);
endp;


/* Workable in Light version. (by using the means of each critical value many times) */

proc JB2a(x,times,meantimes,pp);
   local n,xbar,s,z,sk,kur,stat,statm,i,j,cval;
   if ismiss(x);
      errorlog "Warning: missing data found.";
      x=packr(x);
   endif;
   n=rows(x);
   xbar=meanc(x);
   s=sqrt(sumc((x-xbar)^2)/n);
   z=(x-xbar)/s;
   sk=sumc((z)^3)/n;
   kur=sumc((z)^4)/n;
   stat=n*sk^2/6+n*(kur-3)^2/24;
/* simulation */
   cval=0;
   j=1;
   do while j<=meantimes;
      statm=zeros(times,1);
      i=1;
      do while i<=times;
         x=rndn(n,1);
         x=(x-meanc(x))/stdc(x);
         xbar=meanc(x);
         s=sqrt(sumc((x-xbar)^2)/n);
         z=(x-xbar)/s;
         sk=sumc((z)^3)/n;
         kur=sumc((z)^4)/n;
         statm[i]=n*sk^2/6+n*(kur-3)^2/24;
         i=i+1;
      endo;
      statm=sortc(statm,1);
      cval=cval+statm[round(times*(1-pp))];
      j=j+1;
   endo;
   cval=cval/meantimes;
/* results */
   print "Jarque-Bera Test:";
   print/rz "n             =" n;
   print    "stat          =" stat;
   print    "critical value=" cval;; print/rz "        at" pp*100 "%  (mean)";
   retp(stat~cval);
endp;