/*--------------------------------------------------------------------------
 * NormTest.ox - code for Normality test
 *
 *       code (C) Jurgen Doornik 1994-1999
 *
 * This code implements the results reported in:
 *
 * J.A. Doornik and H. Hansen (1994), `An omnibus test for
 *     univariate and multivariate normality', discussion paper,
 *     Nuffield College, Oxford.
 *
 * Syntax:
 *
 *	NormalityTest(mX, const bPrint)
 *		mX[cT][cX]    in:  matrix to test for normality
 *		bPrint		  in:  TRUE: print test, else don't print
 *	Returns a 2 x 1 matrix: the test statistic and its p-value,.
 *  This function is in testres.h/testres.ox
 *
 *	NormTest(mX)
 *		mX[cT][cX]    in:  matrix to test for normality
 *	Returns the test statistic (a double),
 * 	    which is approximately Chi^2(2*cX) distributed.
 *
 *  AsympNormTest(mX)
 *		mX[cT][cX]    in:  matrix to test for normality
 *	Returns the asymptotic form of the test statistic (a double),
 * 	    which is asymptotically Chi^2(2*cX) distributed.
 *--------------------------------------------------------------------------*/

#include <oxstd.h>
#include "normtest.h"

IN01Data(mX)                             /* transform to approximate N(0,I) */
{
    decl veval, mevec;

    eigensym(correlation(mX), &veval, &mevec);             /* corr.matrix C */

    veval = fabs(veval);
    veval = (veval .> 1e-12) .? 1 ./ sqrt(veval) .: 0;
                /* inv(sqrt(L)) = 1 / dg(sqrt(m)), sqrt(C)=E*inv(sqrt(L))E' */
return standardize(mX) * mevec * diag(veval) * mevec';  /* Xo=Zinv(sqrt(C)) */
}
SkewSU(const vSkew, const cT)
{
    decl n, n2, beta, w2, delta, alfa, vy;
                  /* transform skewness to normality, see D'Agostino [1970] */
	n = double(cT);
	if (n < 8) n = 8;
    n2 = n * n;

    beta = 3 * (n2+27*n-70)/((n-2)*(n+5)) * ((n+1)/(n+7)) * ((n+3)/(n+9));
    w2 = -1 + sqrt(2*(beta - 1));
    delta = 1 / sqrt(log(sqrt(w2)));
    alfa = sqrt(2/(w2 - 1));

    vy = vSkew * sqrt((n+1)*(n+3)/(6*(n-2))) / alfa;
    vy = delta * log(vy + sqrt(vy .^ 2 + 1));

return vy;
}
KurtGamma(const vSkew, const vKurt, const cT)
{
    decl i, n, n2, delta, a, c, k, r, p, vz, kurt;

   /* transform kurtosis to gamma
                     see Shenton and Bowman [1977], Johnson and Kotz [1970] */
	n = double(cT);
	if (n < 8) n = 8;
    n2 = n * n;

	delta = ((n+5)/(n-3)) * ((n+7)/(n+1)) / (6*(n2+15*n-4));
    a = (n-2) * (n2+27*n-70) * delta;
    c = (n-7) * (n2+2*n-5) * delta;
    k = (n*n2+37*n2+11*n-313) * delta / 2;
    r = 1 + c / k;
    p = 3 * (n-1)/(n+1) - r *6*(n-2)/((n+1)*(n+3));
    vz = a + (vSkew .^ 2) * c;                                      /* alfa */

    kurt = (vKurt - 1 - vSkew .^ 2) * k * 2;   /*  gamma standard. to chisq */
	kurt = matrix(kurt);

    for (i = 0; i < columns(kurt); ++i)       /* individual Wilson-Hilferty */
        kurt[0][i] = ( ((kurt[0][i] / (2 * vz[0][i])) ^ (1/3)) - 1 +
            1/(9*vz[0][i])) * sqrt(9*vz[0][i]);

return kurt;
}
NormTest(mX)
{
    decl ct, vskew, vkurt, newskew, newkurt;

    ct = rows(mX);
    mX = IN01Data(mX);

    vskew = sumc(mX .^ 3) / ct;
    vkurt = sumc(mX .^ 4) / ct;

    newskew = SkewSU(vskew, ct);         /* transform skewness to normality */
    newkurt = KurtGamma(vskew, vkurt, ct);  /* transform kurtosis via gamma */

return double(sumr(newskew .^ 2) + sumr(newkurt .^ 2));
}
AsympNormTest(mX)
{
    decl ct, vskew, vkurt;

    ct = rows(mX);
    mX = IN01Data(mX);

    vskew = sumc(mX .^ 3) / ct;
    vkurt = sumc(mX .^ 4) / ct - 3;

return double(ct * sumr(vskew .^ 2) / 6 + ct * sumr(vkurt .^ 2) / 24);
}
/*----------------------------- END NormTest -------------------------------*/
