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

/*-------------------------------- GetFapprox ------------------------------*/
GetFapproximation(const cW, const cT, const cYa, const cWa,
	const adRsq, const adF1, const adF2)
{
    decl k = cW, r = cYa, p = cWa, s, q, n, u;

    s = (r .<= 2 .&& p .<= 1) .|| (r .<= 1 .&& p .<= 2)
		.? 1 .:	sqrt( (sqr(r) .* sqr(p) - 4) ./ (sqr(r) + sqr(p) - 5) );
    q = r .* p / 2 - 1;
    n = cT - k - p - (r - p + 1) / 2;

    u = (1 - adRsq[0]) .^ (1 ./ s);
    adRsq[0] = ((1 - u) ./ u) .* ((n .* s - q) ./ (r .* p));
	r = r .* p;
	s = n .* s - q;
    adF1[0] = sizerc(r) <= 1 ? int(r) : floor(r);
    adF2[0] = sizerc(s) <= 1 ? int(s) : floor(s);
}
/*---------------------------- END GetFapprox ------------------------------*/

/*-------------------------- PrintTestValFixed -----------------------------*/
PrintTestValFixed(const sPre, const dTest, const cR, const cTdf, const sPost)
{
    decl pval, s;

    if (cR <= 0 || dTest < 0)
		return;
    if (cTdf > 0)
    {   
        pval = tailf(dTest, cR, cTdf);
		s = sprint(" F(", cR, ",", cTdf, ")");
        print(sPre, "%-10s", s, "=", "%#9.5g", dTest, " [", "%.4f", pval, "]");
    }
    else
    {   pval = tailchi(dTest, cR);
		s = sprint(" Chi^2(", cR, ")");
        print(sPre, "%-10s", s, "=", "%#9.5g", dTest, " [", "%.4f", pval, "]");
    }
    print(pval < 0.01 ? "**" : pval < 0.05 ? "* " : "  ", sPost);
}
/*------------------------- END PrintTestValFixed --------------------------*/

/*--------------------------- PortmanteauTest ------------------------------*/
PortmanteauTest(const mY, iLag, const cDfLoss, const iPrint)
{
    decl i, lbtest, v_0i, df, cy, ct, ytyinv, pval;

    cy = columns(mY);  ct = rows(mY);
	if (iLag >= ct - 1) iLag = ct - 1;

    for (i = 1, lbtest = 0, ytyinv = invertgen(mY, 4); i <= iLag; ++i)
    {
        v_0i = mY'lag0(mY, i);                 /* missing ones are set to 0 */
        lbtest += 1 / (ct - i) * trace(v_0i' ytyinv * v_0i * ytyinv);
    }
    lbtest *= sqr(ct);
    df = cy * cy * iLag - cDfLoss;
	pval = tailchi(lbtest, df);

    if (iPrint == 2)
    {	
		if (cy == 1)
		{	decl vacf = acf(mY, iLag);
			println("autocorrelation function (ACF) from lag 1 to ", iLag, ":",
				vacf[1:]');
			vacf = pacf(vacf);
			if (rows(vacf))
				println("Partial autocorrelation function (PACF):",	vacf[1:]');
		}
	}
    if (iPrint >= 1)
    {
	 	decl s = sprint(cy > 1 ? "Vector " : "", "Portmanteau(", "%2d", iLag, "):");
		decl s1 = sprint(cy > 1 ? "%-24s" : "%-17s", s);
		PrintTestValFixed(s1, lbtest, df, 0, "\n");
		if (df <= 0 && cDfLoss)
			println("Portmanteau: degrees of freedom correction >= n^2*lag");
	}

return lbtest | pval;
}
/*-------------------------- END PortmanteauTest ---------------------------*/

/*----------------------------- NormalityTest ------------------------------*/
NormalityTest(mX, const iPrint)
{
    decl cy = columns(mX), test, df, pval;

    df = 2 * cy;
    test = NormTest(mX);
	pval = tailchi(test, df);

	if (iPrint == 2)
	{
    	decl ms = IN01Data(mX), ct = rows(ms);

	    decl vskew = sumc(ms .^ 3) / ct;
    	decl vkurt = sumc(ms .^ 4) / ct;
	    decl newskew = SkewSU(vskew, ct);
    	decl newkurt = KurtGamma(vskew, vkurt, ct);

		if (cy == 1)
		{	vskew = double(vskew);
			vkurt = double(vkurt);
			decl mean = double(meanc(mX));
			println("Observations    ", "%12d", ct);
			println("Mean            ", "%#12.5g", fabs(mean) < 1e-10 ? 0.0 : mean);
			println("Std.Devn.       ", "%#12.5g", double(sqrt(varc(mX))));
			println("Skewness        ", "%#12.5g", vskew);
			println("Excess Kurtosis ", "%#12.5g", vkurt - 3);
			println("Minimum         ", "%#12.5g", min(mX));
			println("Maximum         ", "%#12.5g", max(mX));
			PrintTestValFixed("Asymptotic test: ", ct * (sqr(vskew)/6 + sqr(vkurt-3)/24), 2, 0, "\n");
		}
		else
		{	println("Skewness", vskew,
				"Excess kurtosis", vkurt);
			println("Skewness (transformed)", newskew,
				"Excess kurtosis (transformed)", newkurt);
		}
	}
	if (iPrint)
	{	decl s = sprint(cy > 1 ? "Vector " : "", "Normality test:");
		decl s1 = sprint(cy > 1 ? "%-24s" : "%-17s", s);
		PrintTestValFixed(s1, test, df, 0, "\n");
	}
	
return test | pval;
}
/*--------------------------- END NormalityTest ----------------------------*/

/*------------------------------- ArchTest ---------------------------------*/
ArchTest(const mY, const iAr1, const iAr2, const cDfLoss, const iPrint)
{
    decl i, my2 = sqr(mY), cy, car, ct, w, y, mpi, ftest, rss, urss,
		df1, df2, sign, usign, pval, wtwi;

	ct = rows(mY) - iAr2;	cy = columns(mY);
    car = iAr2 - iAr1 + 1;
    for (w = 1, i = iAr1; i <= iAr2; i++)  		 /* append lagged residuals */
        w ~= lag0(my2, i);          		   /* missing ones are set to 0 */

	w = w[iAr2:][];  my2 = my2[iAr2:][];
	olsc(my2, w, &mpi, &wtwi);
	w = my2 - w * mpi; 					                       /* residuals */
    rss = w'w;                                                       /* V'V */
	my2 -= meanc(my2);
	urss = my2'my2;

	ftest = 1 - exp(logdet(rss / ct, &sign) - logdet(urss / ct, &usign));
    if (sign != 1)
        println("ArchTest(): problems with det|V\'V|");
    GetFapproximation(1 + cDfLoss, ct, cy, car * cy, &ftest, &df1, &df2);
	pval = tailf(ftest, df1, df2);

	if (iPrint == 2 && cy == 1)
	{	
		decl cfrom = 1;
		decl se = sqrt(diagonal(wtwi) * rss / df2)';
		print("\nARCH coefficients:");
		println("%c", {"Lag", "Coefficient", "Std.Error"},
			"%cf", {"%5.0f","%12.5g","%11.4g"},
			range(iAr1, iAr2)' ~ mpi[cfrom : ][] ~ se[cfrom : ][]);
		println("RSS = ", double(rss), "  sigma = ", double(sqrt(rss / df2)));

		println("\nTesting for ", cy > 1 ? "Vector " : "",
			"error ARCH from lags ", iAr1, " to ", iAr2);
	}
    if (iPrint >= 1)
    {
	 	decl s = sprint(cy > 1 ? "Vector " : "", "ARCH ", iAr1, "-", iAr2, " test:");
		decl s1 = sprint(cy > 1 ? "%-24s" : "%-17s", s);
		PrintTestValFixed(s1, ftest, df1, df2, "\n");
	}
return ftest | pval;
}
/*----------------------------- END ArchTest -------------------------------*/
