/*--------------------------------------------------------------------------
 * coigamma.ox - code for asymptotic distribution of I(1) and I(2) tests
 *
 *       code (C) Jurgen Doornik 1997-99
 *
 *
 * This code implements the results reported in:
 *
 * Doornik, J.A. (1998), `Approximations to the Asymptotic	Distribution
 *   of Cointegration Tests', Journal of Economic Surveys, 12, 573--593.
 * Boswijk, H.P. and Doornik, J.A. (1999),  `Distribution Approximations
 *   for Cointegration Tests with Stationary Exogenous Regressors', mimeo.
 *
 * Syntax:
 *
 * 	GammaParTraceAsymp(const iCase, const cN, const cT);
 *		returns a 1x2 matrix with r and a for the Trace test, I(1) analysis
 *     	(r and a are the parameters for probgamma and quangamma)
 *
 *		iCase refers to treatment of Constant and Trend:
 *					   Trend         Constant      Johansen case
 * 			0 = H_z;   unrestricted  unrestricted  case 0;
 *          1 = H_c;   restricted    unrestricted  case 1*;
 *          2 = H_lc;  none          unrestricted  case 1;
 *          3 = H_l;   none          restricted    case 2*;
 *          4 = H_ql;  none          none          case 2
 *		cN is the dimension (>= 1), n-p
 *		cT argument is currently not used
 *
 * 	GammaParMaxEvAsymp(const iCase, const cN, const cT);
 *     	returns a 1x2 matrix with r and a for the MaxEv test, I(1) analysis
 *
 * 	GammaParI2Asymp(const iCase, const iTestCase, const cS, const cNminS);
 *		returns a 1x2 matrix with r and a for the Trace test, I(2) analysis
 * 		iCase      refers to treatment of Constant and Trend, as above
 *		iTestCase: 1 = S_p;  0 = Q_p
 *		cS         number of I(1) components: s
 *		cNminS     number of I(2) components: n-p-s.
 *
 * GammaParI1CondAsymp(const iCase, const cN, const cN1, const cT)
 * 		iCase      refers to treatment of Constant and Trend, as above
 *				   (H_lc and H_ql are not defined here).
 *      cN         dimension of full system
 * 		cN1        dimension in partial system
 *      cT         sample size (normally 0)
 *
 * GammaParI0CondAsymp(const iCase, const vCorr, const cT)
 * 		iCase      refers to treatment of Constant and Trend, as above
 *				   (H_lc and H_ql are not defined here).
 *      vCorr      vector with correlations; cN = sizerc(vCorr)
 *      cT         sample size (normally 0)
 *--------------------------------------------------------------------------*/

static decl s_mTrace_m_coef = < // m = 
  //  n^2     n        1    n==1     n==2  n^1/2
      2,  -1.00,    0.07,   0.07,       0,     0;
      2,   2.01,       0,   0.06,    0.05,     0;
      2,   1.05,   -1.55,  -0.50,   -0.23,     0;
      2,   4.05,    0.50,  -0.23,   -0.07,     0;
      2,   2.85,   -5.10,  -0.10,   -0.06,  1.35>;
static decl s_mTrace_v_coef = < // v = 
      3,  -0.33,  -0.55,      0,        0,     0;
      3,    3.6,   0.75,   -0.4,     -0.3,     0;
      3,    1.8,      0,   -2.8,     -1.1,     0;
      3,    5.7,    3.2,   -1.3,     -0.5,     0;
      3,    4.0,    0.8,   -5.8,    -2.66,     0>;
	  
static decl s_mTrace_m_time = < // log(m)-log(m-hat) =
//sqrt(n)/T   n/T  n^2/T^2   n==1/T     n==1     n==2     n==3
  -0.101,   0.499,   0.896,  -0.562, 0.00229, 0.00662,       0; 
       0,   0.465,   0.984,  -0.273,-0.00244,       0,       0; 
   0.134,   0.422,    1.02,    2.17,-0.00182,       0,-0.00321; 
  0.0252,   0.448,    1.09,  -0.353,       0,       0,       0; 
  -0.819,   0.615,   0.896,    2.43, 0.00149,       0,       0>;
static decl s_mTrace_v_time = < // log(v)-log(v-hat) =
  -0.204,   0.980,    3.11,   -2.14,  0.0499, -0.0103,-0.00902; 
   0.224,   0.863,    3.38,  -0.807,       0,       0, -0.0091; 
   0.422,   0.734,    3.76,    4.32,-0.00606,       0,-0.00718; 
       0,   0.836,    3.99,   -1.33,-0.00298,-0.00139,-0.00268; 
   -1.29,    1.01,    3.92,    4.67, 0.00484,-0.00127, -0.0199>;
	
GammaParTraceAsymp(const iCase, const cN, const cT)
{
	decl x = sqr(cN) ~ cN ~ 1 ~ (cN .== 1) ~ (cN .== 2) ~ sqrt(cN);
	decl m = x * s_mTrace_m_coef[iCase][]';
	decl v = x * s_mTrace_v_coef[iCase][]';

	if (cT > 0)
    {	
		x = sqrt(cN)./cT ~ cN./cT ~ sqr(cN./cT) ~
			(cN.==1)./cT ~ (cN.==1) ~ (cN.==2) ~ (cN.==3);
	    m = m .* exp(x * s_mTrace_m_time[iCase][]');
    	v = v .* exp(x * s_mTrace_v_time[iCase][]');
	}
	return sqr(m)./v ~ m./v;
}

static decl s_mMaxev_m_coef = < // m = 
       //   n        n^1/2         n==1         n==2            1
       6.0019,     -2.7764,     0.67185,     0.11490,     -2.7558;  
       5.9498,     -2.3669,    0.048360,    0.018198,     0.43402;  
       5.8271,     -1.5666,     -1.6118,    -0.25949,     -1.6487;  
       5.8658,     -1.7552,    -0.34443,   -0.077991,      2.5595;  
       5.6364,    -0.21447,     -3.5166,    -0.47966,    -0.90531>; 
static decl s_mMaxev_v_coef = < // v = 
       1.8806,      14.714,      1.1136,    0.070508,     -15.499;  
       2.2231,      12.058,     0.58592,   -0.034324,     -7.9064;  
       2.0785,      13.074,     -3.3680,    -0.24528,     -9.7846;  
       1.9955,      12.841,      1.2425,     0.41949,     -5.5428;  
       2.0899,      12.393,     -7.1523,    -0.25260,     -5.3303>; 
	
GammaParMaxEvAsymp(const iCase, const cN, const cT)
{
	decl x = cN ~ sqrt(cN) ~ (cN .== 1) ~ (cN .== 2) ~ 1;
	decl m = x * s_mMaxev_m_coef[iCase][]';
    decl v = x * s_mMaxev_v_coef[iCase][]';
	return sqr(m)./v ~ m./v;
}

static decl s_mQtrace_m_coef = < // Qm =  
   // n^2            n            1       n(n-s)      (n-s)^2          n-s         s==0 
   1.9906,      2.1248,    -0.42846,      1.9998,    0.017957,    -0.20907,     0.00000;
   1.9966,      1.0758,     -1.7843,      1.9980,    0.014504,     -2.1649,      1.4170;
   1.9954,      4.1140,     0.19778,      2.0043,    0.010704,    -0.21254,     0.00000>;
static decl s_mQtrace_v_coef = < // Qv =  
   3.1448,      2.2065,      4.9710,      3.1351,     0.00000,     0.76645,     0.00000;
   2.8999,      4.1618,     -5.6183,      3.2175,     0.00000,     -2.9496,      12.173;
   2.9350,      6.5371,      1.9913,      3.4897,   -0.068124,     0.00000,     -1.0696>;
static decl s_mStrace_m_coef = < // Sm =  
   1.9906,      2.1248,    -0.42846,      1.9946,      2.0112,     -1.1079,     0.00000;
   1.9966,      1.0758,     -1.7843,      1.9929,      2.0073,   -0.060440,      1.3392;
   1.9954,      4.1140,     0.19778,      2.0042,      2.0091,      1.8118,     0.00000>;
static decl s_mStrace_v_coef = < // Sv =  
   3.1448,      2.2065,      4.9710,      2.3807,      2.0770,     0.00000,     0.00000;
   2.8999,      4.1618,     -5.6183,      2.2688,      2.2175,     0.00000,      9.8887;
   2.9350,      6.5371,      1.9913,      2.6055,      1.9891,      2.6848,     0.00000>;
  
GammaParI2Asymp(const iCase, const iTestCase, const cS, const cNminS)
{
    decl x, m, v, n = cNminS + cS;

	if (cNminS == 0)
		return GammaParTraceAsymp(iCase, n, 0);

    x = sqr(n) ~ n ~ 1 ~ n.*cNminS ~ sqr(cNminS) ~ cNminS ~ (cS.==0);

	if (iTestCase == 0)
	{
		m = x * s_mQtrace_m_coef[iCase-1][]';
   		v = x * s_mQtrace_v_coef[iCase-1][]';
	}
	else
	{
		m = x * s_mStrace_m_coef[iCase-1][]';
   		v = x * s_mStrace_v_coef[iCase-1][]';
	}
	return sqr(m)./v ~ m./v;
}	   

static decl s_mTrace_cov_coef = < // cov = 
	-1.270, -1.066, ., -1.35, . >;

GammaParI1CondAsymp(const iCase, const cN, const cN1, const cT)
{
	decl parn = GammaParTraceAsymp(iCase, cN, cT);
	decl cov = s_mTrace_cov_coef[iCase];// covariance term
	decl mean = parn[][0] ./ parn[][1];
	decl var = parn[][0] ./ sqr(parn[][1]);
	if (any(cN1 .> cN))
	{	println("\ninput error");
		return;
	}
	mean = cN1 ./ cN .* mean;
	var = cN1 ./ cN .* var - (cN - cN1) .* cN1 * cov;

	return sqr(mean)./var ~ mean./var;
}

GammaParI0CondAsymp(const iCase, const vCorr, const cT)
{
	decl rho = vCorr[], rho2 = sqr(rho);// force into row
	decl cn = columns(rho);
	decl q = (iCase == 1 || iCase == 3) ? cn + 1 : cn;
	decl sumrho2 = double(sumr(rho2));	// sum of squared rho's
	decl ra = GammaParTraceAsymp(iCase, cn, cT);
	decl m = ra[][0] ./ ra[][1];		// mean of uncorrelated case
	decl v = ra[][0] ./ sqr(ra[][1]);	// var of uncorrelated case
	decl cov = s_mTrace_cov_coef[iCase];// covariance term
	decl mc, vc, cc;
	cc = setdiagonal(rho2'rho2, zeros(1,cn));
	mc = sumrho2 / cn * m + (1 - sumrho2 / cn) * cn * q;
	vc = sumsqrr(rho2) * (v / cn - (cn - 1) * cov) + sumc(sumr(cc)) * cov
		+ 4 * sumr(rho2 .* (1 - rho2)) * m / cn
		+ 2 * q * sumr(sqr(1 - rho2));

	return sqr(mc)./vc ~ mc./vc;
}
