function [H,P,JBSTAT,CV] = jbtest(x,alpha,mctol)
%JBTEST Jarque-Bera hypothesis test of composite normality.
%   H = JBTEST(X) performs the Jarque-Bera goodness-of-fit test of composite
%   normality, i.e., that the data in the vector X came from an unspecified
%   normal distribution, and returns the result of the test in H. H=0
%   indicates that the null hypothesis ("the data are normally distributed")
%   cannot be rejected at the 5% significance level. H=1 indicates that the
%   null hypothesis can be rejected at the 5% level.
%
%   JBTEST treats NaNs in X as missing values, and ignores them.
%
%   The Jarque-Bera test is a 2-sided goodness-of-fit test suitable for
%   situations where a fully-specified null distribution is not known, and its
%   parameters must be estimated.  For large sample sizes, the test statistic
%   has a chi-square distribution with two degrees of freedom. Critical
%   values, computed using Monte-Carlo simulation, have been tabulated for
%   sample sizes N <= 2000 and significance levels 0.001 <= ALPHA <= 0.50.
%   JBTEST computes a critical value for a given test by interpolating into
%   that table, using the analytic approximation to extrapolate for larger
%   sample sizes.
%
%   The Jarque-Bera hypotheses and test statistic are:
%
%             Null Hypothesis:  X is normally distributed with unspecified
%                               mean and standard deviation.
%      Alternative Hypothesis:  X is not normally distributed.  The test is
%                               specifically designed for alternatives in the
%                               Pearson family of distributions.
%              Test Statistic:  JBSTAT = N*(SKEWNESS^2/6 + (KURTOSIS-3)^2/24),
%                               where N is the sample size and the kurtosis of
%                               the normal distribution is defined as 3.
%
%   H = JBTEST(X,ALPHA) performs the test at significance level ALPHA.  ALPHA
%   is a scalar in the range 0.001 <= ALPHA <= 0.50.  To perform the test at
%   significance levels outside that range, use the MCTOL input argument.
%
%   [H,P] = JBTEST(...) returns the p-value P, computed using inverse
%   interpolation into the look-up table of critical values. Small values of P
%   cast doubt on the validity of the null hypothesis. JBTEST warns when P is
%   not found within the limits of the table, i.e., outside the interval
%   [0.001, 0.50], and returns one or the other endpoint of that interval. In
%   this case, you can use the MCTOL input argument to compute a more
%   accurate value.
%
%   [H,P,JBSTAT] = JBTEST(...) returns the test statistic JBSTAT.
%
%   [H,P,JBSTAT,CRITVAL] = JBTEST(...) returns the critical value CRITVAL for
%   the test. When JBSTAT > CRITVAL, the null hypothesis can be rejected at a
%   significance level of ALPHA.
%
%   [H,P,...] = JBTEST(X,ALPHA,MCTOL) computes a Monte-Carlo approximation
%   for P directly, rather than using interpolation of the pre-computed
%   tabulated values.  This is useful when ALPHA or P is outside the range of
%   the look-up table.  JBTEST chooses the number of MC replications, MCREPS,
%   large enough to make the MC standard error for P, SQRT(P*(1-P)/MCREPS),
%   less than MCTOL.
%
%   See also LILLIETEST, KSTEST, KSTEST2, CDFPLOT.

%   Copyright 1993-2006 The MathWorks, Inc.
%   $Revision: 1.4.2.5 $   $ Date:  $

% References:
%   Jarque, C.M., and A.K. Bera (1987), "A test for normality of observations
%     and regression residuals," International Statistical Review.  This
%     paper proposed the original test.
%   Deb, P., and M. Sefton (1996), "The distribution of a Lagrange
%     multiplier test of normality," Economics Letters.  This paper proposed
%     a Monte Carlo simulation for determining the distribution of the
%     test statistic.  The jbtest results are based on an independent
%     Monte Carlo simulation, not the one from this paper.

if ~isvector(x) || ~isreal(x)
    error('stats:jbtest:BadInput', ...
          'Input X must be a vector of real values.');
end
% Remove missing observations indicated by NaN's.
x = x(~isnan(x));

ssize = length(x);
if ssize < 2
    error('stats:jbtest:NotEnoughData',...
          'Sample vector X must have at least 2 valid observations.');
elseif ssize == 2
    % The J-B stat is a constant when the sample size is 2
    H = 0;
    P = 1;
    JBSTAT = 1/3;
    CV = NaN;
    return;
end

% Ensure the significance level ALPHA is a scalar.  Set default if necessary.
% Monte Carlo results are available only for 0.001 <= ALPHA <= 0.25.
if (nargin < 2) || isempty(alpha)
    alpha  =  0.05;
elseif ~isscalar(alpha) || ~(0 < alpha && alpha < 1)
    error('stats:jbtest:BadAlpha',...
          'Significance level ALPHA must be a scalar between 0 and 1.');
end

if nargin < 3 || isempty(mctol)
    mctol = [];
else
   if ~isscalar(mctol) || mctol<=0
      error('stats:jbtest:BadStdErrTol',...
            'Monte-Carlo standard error tolerance MCTOL must be a positive scalar value.');
   end
end

% Compute moments and test statistic.
ssize = length(x);
zscores = (x-mean(x))/std(x,1);
skew = sum(zscores.^3)/ssize;
kurt = sum(zscores.^4)/ssize - 3;

JBSTAT = ssize*(skew.*skew/6 + kurt.*kurt/24);

if isempty(mctol)
    % Get a row of the critical value table for the current sample size.
    [alphas,CVs] = CVtbl(ssize);

    % Make sure requested alpha is within the look-up table.
    if (alpha < alphas(1)) || (alphas(end) < alpha)
        error('stats:jbtest:BadAlpha',...
              ['Significance level ALPHA is outside the range of the tabulated values.\n' ...
               'Use a value in the interval [%g, %g], or use the MCTOL input argument.'], ...
               alphas(1),alphas(end));
    end

    %  1-D interpolation into the tabulated quantiles.
    CV = spline(alphas, CVs, alpha);

    % Compute the P-value. Warn if the P-value is not found within the
    % available 'alphas' of the table and return one of the extremes.
    if (JBSTAT < CVs(end)) % smallest critval at end
        warning('stats:jbtest:OutOfRangeP',...
                'P is greater than the largest tabulated value, returning %.3g.',alphas(end));
        P = alphas(end);
    elseif (CVs(1) < JBSTAT) % largest critval at beginning
        warning('stats:jbtest:OutOfRangeP',...
                'P is less than the smallest tabulated value, returning %.3g.',alphas(1));
        P = alphas(1);
    else
        P = spline(CVs, alphas, JBSTAT);
    end


% Compute the critical value and p-value on the fly using Monte-Carlo simulation.
else
    [CV,P] = jbMC(JBSTAT,ssize,alpha,mctol);
end

% Returning "H = 0" implies that we "Do not reject the null hypothesis at the
% significance level of alpha" and "H = 1" implies that we "Reject the null
% hypothesis at significance level of alpha."
H = double(JBSTAT > CV);


% ----------------------------------------------------------------

function [alphas,CVs] = CVtbl(n)
% Tabulated critical values for JB test

% Significance levels used in the simulation -- logspace(-3,-1,10), plus the
% other common levels.
alphas = [0.001 0.0016681 0.0027826 0.0046416 0.0077426 0.01 0.012915 ...
    0.021544 0.025 0.035938 0.05 0.059948 0.1 0.15 0.2 0.25 0.50];

% Sample sizes used in the simulation.
sampleSizes = [3 4 5 10 15 20 25 30 35 40 45 50 60 70 80 90 100 ...
    125 150 175 200 250 300 400 500 800 1000 1200 1400 1600 1800 2000 inf];

% Critical value table for the above sample sizes (rows) and significance
% levels (cols). Last row is chi2inv(1-alpha,2), an asymptotic approximation.
criticalValues = [ ...
 0.5312  0.5312  0.5312  0.5312  0.5312  0.5312  0.5311 0.5310 0.5309 0.5305 0.5297 0.5290 0.5251 0.5176 0.5074 0.4946 0.4063; ...
 0.9606  0.9590  0.9564  0.9520  0.9448  0.9396  0.9329 0.9133 0.9056 0.8817 0.8519 0.8316 0.7553 0.6721 0.6303 0.5947 0.4739; ...
 1.8289  1.8053  1.7727  1.7276  1.6661  1.6275  1.5828 1.4723 1.4343 1.3294 1.2185 1.1516 0.9442 0.7945 0.7302 0.6878 0.5285; ...
10.9719  9.8430  8.6751  7.4815  6.2927  5.7077  5.1350 4.0456 3.7474 3.0691 2.5239 2.2555 1.6231 1.2821 1.1235 1.0198 0.6951; ...
19.5425 16.7207 14.0454 11.5527  9.2814  8.2365  7.2613 5.5238 5.0729 4.0773 3.2985 2.9215 2.0533 1.5965 1.3779 1.2336 0.7916; ...
25.0722 21.0001 17.2981 13.9762 11.0602  9.7531  8.5489 6.4401 5.8998 4.7176 3.8011 3.3596 2.3505 1.8239 1.5631 1.3885 0.8577; ...
28.4885 23.6332 19.2818 15.4635 12.1658 10.7058  9.3658 7.0389 6.4463 5.1501 4.1494 3.6684 2.5707 1.9986 1.7063 1.5079 0.9071; ...
30.6106 25.2501 20.5242 16.4088 12.8805 11.3263  9.9063 7.4472 6.8222 5.4578 4.4039 3.8968 2.7431 2.1390 1.8216 1.6040 0.9457; ...
31.9343 26.2695 21.3105 17.0176 13.3572 11.7488 10.2818 7.7395 7.0942 5.6856 4.5973 4.0734 2.8827 2.2547 1.9172 1.6835 0.9771; ...
32.7514 26.9025 21.8016 17.4145 13.6749 12.0355 10.5399 7.9520 7.2949 5.8598 4.7481 4.2126 2.9987 2.3524 1.9980 1.7508 1.0033; ...
33.2273 27.2686 22.1030 17.6594 13.8848 12.2302 10.7201 8.1098 7.4469 5.9945 4.8689 4.3258 3.0973 2.4361 2.0674 1.8085 1.0256; ...
33.4825 27.4765 22.2863 17.8172 14.0345 12.3739 10.8586 8.2340 7.5661 6.1045 4.9697 4.4214 3.1834 2.5097 2.1283 1.8592 1.0449; ...
33.5009 27.5458 22.3790 17.9431 14.1788 12.5255 11.0182 8.3968 7.7287 6.2620 5.1203 4.5688 3.3246 2.6316 2.2297 1.9434 1.0768; ...
33.2610 27.3738 22.2848 17.9209 14.2092 12.5801 11.0884 8.4921 7.8286 6.3696 5.2305 4.6803 3.4374 2.7298 2.3115 2.0114 1.1023; ...
32.8742 27.0975 22.1135 17.8305 14.1911 12.5841 11.1155 8.5520 7.8938 6.4462 5.3134 4.7662 3.5292 2.8104 2.3789 2.0674 1.1231; ...
32.3418 26.7278 21.8770 17.6977 14.1283 12.5542 11.1112 8.5850 7.9361 6.5028 5.3796 4.8378 3.6071 2.8788 2.4361 2.1151 1.1408; ...
31.8884 26.3990 21.6407 17.5510 14.0491 12.5067 11.0882 8.6004 7.9589 6.5429 5.4314 4.8957 3.6730 2.9372 2.4851 2.1557 1.1557; ...
30.6089 25.4748 21.0100 17.1460 13.8271 12.3565 10.9976 8.6020 7.9816 6.6071 5.5277 5.0071 3.8025 3.0523 2.5819 2.2363 1.1852; ...
29.4290 24.6043 20.3979 16.7425 13.5878 12.1804 10.8779 8.5719 7.9718 6.6387 5.5919 5.0863 3.8987 3.1384 2.6543 2.2964 1.2072; ...
28.4225 23.8582 19.8658 16.3846 13.3664 12.0139 10.7589 8.5307 7.9496 6.6571 5.6408 5.1474 3.9732 3.2050 2.7106 2.3436 1.2243; ...
27.5140 23.1914 19.3922 16.0635 13.1643 11.8602 10.6494 8.4904 7.9259 6.6688 5.6783 5.1957 4.0327 3.2584 2.7559 2.3812 1.2382; ...
26.0239 22.0764 18.5877 15.5123 12.8124 11.5923 10.4542 8.4114 7.8751 6.6803 5.7343 5.2681 4.1224 3.3402 2.8250 2.4388 1.2591; ...
24.8353 21.1758 17.9384 15.0633 12.5206 11.3653 10.2839 8.3384 7.8276 6.6842 5.7728 5.3194 4.1873 3.3988 2.8749 2.4808 1.2744; ...
23.0771 19.8542 16.9669 14.3843 12.0743 11.0157 10.0226 8.2266 7.7526 6.6877 5.8248 5.3889 4.2748 3.4790 2.9434 2.5382 1.2956; ...
21.8170 18.9068 16.2666 13.8903 11.7469 10.7604  9.8318 8.1454 7.6984 6.6882 5.8581 5.4336 4.3320 3.5318 2.9886 2.5760 1.3097; ...
19.5539 17.1748 14.9865 12.9765 11.1441 10.2914  9.4852 8.0027 7.6040 6.6858 5.9096 5.5043 4.4246 3.6180 3.0631 2.6388 1.3334; ...
18.6707 16.4879 14.4745 12.6143 10.9072 10.1106  9.3521 7.9501 7.5691 6.6847 5.9282 5.5296 4.4580 3.6500 3.0907 2.6624 1.3421; ...
18.0126 15.9830 14.0990 12.3493 10.7362  9.9803  9.2587 7.9127 7.5441 6.6832 5.9408 5.5471 4.4814 3.6718 3.1099 2.6787 1.3484; ...
17.5190 15.6036 13.8149 12.1501 10.6072  9.8811  9.1869 7.8828 7.5239 6.6800 5.9482 5.5587 4.4979 3.6876 3.1237 2.6904 1.3530; ...
17.1305 15.3029 13.5883 11.9883 10.5059  9.8072  9.1333 7.8611 7.5089 6.6789 5.9546 5.5677 4.5105 3.6997 3.1345 2.6996 1.3564; ...
16.8158 15.0600 13.4102 11.8682 10.4296  9.7478  9.0911 7.8446 7.4979 6.6777 5.9600 5.5755 4.5214 3.7101 3.1433 2.7072 1.3595; ...
16.5585 14.8633 13.2657 11.7666 10.3636  9.6981  9.0549 7.8296 7.4872 6.6763 5.9635 5.5806 4.5289 3.7173 3.1499 2.7129 1.3619; ...
13.8155 12.7921 11.7687 10.7454 9.7220   9.2103  8.6987 7.6753 7.3778 6.6519 5.9915 5.6286 4.6052 3.7942 3.2189 2.7726 1.3863]; 

% Interpolate a row of critical values for the given sample size.
CVs = spline(1./sampleSizes,criticalValues',1./n);

% ----------------------------------------------------------------

function [crit,p] = jbMC(JBstat,n,alpha,mctol)
%JBMC Simulated critical values and p-values for the Jarque-Bera test.
%   [CRIT,P] = JBMC(JBSTAT,N,ALPHA,MCTOL) returns the critical value CRIT
%   and p-value P for the Jarque-Bera test of the null hypothesis that data
%   were drawn from an unspecified normal distribution, for a sample size N
%   and confidence level 100*(1-ALPHA)%.  P is the p-value for the observed
%   value JBSTAT of the Jarque-Bera statistic.  JBMC uses Monte-Carlo
%   simulation to approximate CRIT and P, and chooses the number of MC
%   replications, MCREPS, large enough to make the standard error for P,
%   SQRT(P*(1-P)/MCREPS), less than MCTOL.
%
%   This function can used to recreate table 1 from Deb and Sefton (1996).

if n < 500 % it's faster to block the trials for small sample sizes
    k = 20;
    rep = ones(1,n);
else
    k = 1;
    rep = 1;
end

vartol = mctol^2;

crit = 0;
p = 0;
mcRepsTot = 0;
mcRepsMin = 1000;
while true
    mcRepsOld = mcRepsTot;
    mcReps = k*ceil((mcRepsMin - mcRepsOld) / k);
    JBstatMC = zeros(mcReps,1);
    for i = 1:k:length(JBstatMC)
        x = randn(k,n);
        xbar = sum(x,2)/n;      % sample mean
        devs = x - xbar(:,rep); % deviations from sample mean
        sqdevs = devs.*devs;
        s2 = sum(sqdevs,2)/n;   % sample variance

        skew = sum(devs.*sqdevs,2)./(s2.^1.5 * n);
        kurt = sum(sqdevs.*sqdevs,2)./(s2.^2 * n) - 3;
        JBstatMC(i:i+k-1) = n*(skew.*skew/6 + kurt.*kurt/24);
    end
    critMC = prctile(JBstatMC,100*(1-alpha));
    pMC = sum(JBstatMC > JBstat) ./ mcReps;
    
    mcRepsTot = mcRepsOld + mcReps;
    crit = (mcRepsOld*crit + mcReps*critMC) / mcRepsTot;
    p = (mcRepsOld*p + mcReps*pMC) / mcRepsTot;
    
    % Compute a std err for p, with lower bound (1/N)*(1-1/N)/N when p==0.
    sepsq = max(p*(1-p)/mcRepsTot, 1/mcRepsTot^2);
    if sepsq < vartol
        break
    end
    
    % Based on the current estimate, find the number of trials needed to
    % make the MC std err less than the specified tolerance.
    mcRepsMin = 1.2 * (mcRepsTot*sepsq)/vartol;
end
