
function [h chi2jb pval crit chi2jbs]=jarqueberatest(data, alpha)
if nargin == 1
    alpha=0.05;
end
   data = data(:);
   n= length(data);     
   skew = skewness(data);   kurt = kurtosis(data);
   chi2jb  =    n/6 * (skew.^2 + (kurt-3).^2/4 ) ;
   if n> 1000
           pval=1-chi2cdf(chi2jb,2);
            h = (pval < alpha);
            crit=chi2inv(1-alpha,2);
   else
    B=floor(10000000/n);
    s = RandStream('mt19937ar','Seed',0);
    RandStream.setGlobalStream(s);
    norms = normrnd(0,1, [n B]);
    skew = skewness(norms);
    kurt = kurtosis(norms);
    chi2jbs  =    n/6 * (skew.^2 + (kurt-3).^2/4 ) ;
    pval=sum(chi2jbs > chi2jb)/B;
    h = (pval < alpha);
   crit = prctile(chi2jbs, 100*(1-alpha));
end

