#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>
#include <pcfiml.h>

#import <modelbase>
#import <maximize>
#import <lib/testres>
#include <lib/normtest.h>

#include <../lib/coigamma.ox>


PcFiml::PcFiml()
{
    this->Modelbase();
    m_fCointegration = FALSE;  m_iMethod = M_NONE;
    m_fPrint = TRUE;  m_fPrintUrf = FALSE;
	m_asA = 0;
    m_mModsel = 0;
}
PcFiml::printVariables(iGroup, sGroupname)
{
    decl i, csel;
	decl status = sGroupname ? " status " ~ sGroupname : "";

    csel = columns(m_mVarsel);              /* number of selected variables */

    for (i = 0; i < csel; ++i)                         /* walk the database */
        if (m_mSelgroup[0][i] == iGroup)
            println("%-12s", m_asNames[m_mVarsel[0][i]],
                  " lag ", m_mLagsel[0][i], status);
}
PcFiml::Estimate()
{
    decl mu, logdetvtv, sign, iret, ma;

    m_iT1est = m_iT1sel;  m_iT2est = m_iT2sel;
    m_cT = m_iT2est - m_iT1est + 1;

    m_mY = GetGroupLag(Y_VAR, 0, 0);
    m_cY = columns(m_mY);
    m_mW = GetGroupLag(Y_VAR, 1, 10000);
	m_cYlag = columns(m_mW);
	m_mW ~= GetGroup(X_VAR);                               /* append X vars */
    mu = GetGroup(U_VAR);					       /* get unrestricted vars */
    m_cU = columns(mu);
    if (m_cU) m_mW ~= mu; /* note: no vars concentrated out at system level */
    m_cW = columns(m_mW);

    ma = GetGroup(A_VAR);				      /* get additional instruments */
	m_cA = columns(ma);
	m_asA = 0;
	
	if (m_cY == 0)
	{	eprint("Estimate(): no Y's\n");
		return;
	}
	if (m_cT < m_cW)
	{	eprint("Estimate(): not enough observations\n");
		return;
	}
	if (m_cA && m_cA < m_cY - 1)
	{	eprint("Estimate(): IVE needs more instruments than non-modelled endogenous variables\n");
		return;
	}

	if (m_cA)
	{
	    GetGroupNames(A_VAR, &m_asA);/* get names of additional instruments */
 		// data: (y:Y:W:A) = (y:X:A)
		decl mwa = m_mW ~ ma, mpirf;
		if (olsc(m_mY, mwa, &mpirf) < 0)
			eprint("Estimate(): Reduced form W\'W is singular\n");
		decl yhatrf = mwa * mpirf;
		mwa = ma = 0;
		olsc(m_mY[][0], yhatrf[][1:] ~ m_mW, &m_mPi, &m_mWtWinv);
		// adjust Y and W
		m_mW = m_mY[][1:] ~ m_mW;
		m_mY = m_mY[][0];
		m_cW += m_cY - 1;
		m_cY = 1;
	}
	else if (m_cW)
	{
		if (olsc(m_mY, m_mW, &m_mPi, &m_mWtWinv) < 0)
			eprint("Estimate(): W\'W is singular\n");
	}
	if (m_cW)
		m_mResidual = m_mY - m_mW * m_mPi;
	else
		m_mResidual = m_mY, m_mPi = new matrix[0][m_cY], m_mWtWinv = <>;

	SetParCount(sizerc(m_mPi));
	SetPar(vec(m_mPi));
	SetModelStatus(MS_ESTIMATED);

    m_mPiUrf = m_mPi = m_mPi';		 // store equations in rows
    m_mOmega = m_mResidual'm_mResidual;                              /* V'V */
    m_mVtVinvUrf = m_mVtVinv = invert(m_mOmega, &logdetvtv,&sign);/*inv(V'V)*/
    if (sign != 1)
        eprint("Estimate(): problems with det|V\'V|\n");
    if (columns(m_mVtVinv) == 0)
        return;
    m_dLogLikUrf = m_dLogLik = -m_cT / 2 * (logdetvtv - m_cY * log(m_cT));
    m_mOmega /= (m_cT - m_cW);                                 /* V'V/(T-k) */
    m_mVarPi = ( diagonal(m_mOmega) ** diagonal(m_mWtWinv)' )';

    m_fCointegration = FALSE;
	m_iMethod = m_cA ? M_GIVE : M_OLS;

    if (m_fPrint)
        Output(TRUE, FALSE);

	return TRUE;
}
PcFiml::Covar()
{
    if (m_iMethod >= M_FIML)
        m_mCovar = m_mVar;
    else if (m_iMethod != M_NONE)
        m_mCovar = m_mOmega ** m_mWtWinv;

return 0;
}
PcFiml::GetStatus(const aiConst, const aiTrend)
{
    decl ai, i, iconststat, itrendstat, csel;

	aiConst[0] = aiTrend[0] = iconststat = itrendstat = -1;
	
	ai = strfind(m_asNames[m_mVarsel], {"Trend","Constant"});
	i = ai[0];
	if (i >= 0)
    {   aiTrend[0] = i;
		itrendstat = m_mSelgroup[0][i];
	}
	i = ai[1];
	if (i >= 0)
    {   aiConst[0] = i;
		iconststat = m_mSelgroup[0][i];
	}
	if (itrendstat == -1)	 // return case of test: 0,..,4
	{
		if (iconststat == -1)
			i = 0;  // no constant;
		else if (iconststat == U_VAR)
			i = 2;  // unrestricted constant;
		else
			i = 1;  // restricted constant;
	}
	else if (itrendstat == U_VAR)
		i = 4;  // unrestricted trend;
	else
		i = 3;  // restricted trend;

	return i;
}

PcFiml::fCointData(const amDy, const amY_1)
{
    decl  i, w, y_1, y_i, y_i_1, dy, maxylag, pi;

    maxylag = GetMaxGroupLag(Y_VAR);

    y_1 = GetGroupLag(Y_VAR, 1, 1);

    if (m_cY != columns(y_1))
    {   eprint("Variables missing at lag 1\n");
        return FALSE;
    }
    dy  = m_mY - y_1;

    w = GetGroup(U_VAR);
    for (i = 1, y_i = y_1; i < maxylag; i++, y_i = y_i_1)
    {
        y_i_1 = GetGroupLag(Y_VAR, i + 1, i + 1);

        if (m_cY != columns(y_i_1))
        {   eprint("Variables missing at lag ", i + 1, "\n");
            return FALSE;
        }
        w ~= y_i - y_i_1;
    }
	y_1 ~= GetGroup(X_VAR);

    if (columns(w))
    {
		olsc(y_1, w, &pi);
        amY_1[0] = y_1 - w * pi;		                    /* residual y_1 */

		olsc(dy, w, &pi);
        amDy[0] = dy - w * pi;			                     /* residual dy */
    }
    else                                      /* nothing to concentrate out */
    {   amY_1[0] = y_1;                                              /* y_1 */
        amDy[0] = dy;                                                 /* dy */
    }

return TRUE;
}
PcFiml::olsRank(mDy, mY_1, const mW, const amEval, const amBeta, const amAlpha)
{
    decl  pi, s01, eval, beta, ct = rows(mDy);

	if (columns(mW) > 0)
	{	olsc(mDy,  mW, &pi);  mDy  -= mW * pi;
		olsc(mY_1, mW, &pi);  mY_1 -= mW * pi;
	}
	s01 = mDy'mY_1 / ct;
	if (eigensymgen(s01' / ((mDy'mDy)/ct) * s01, mY_1'mY_1/ct, &eval, &beta) != 0)
		eval[] = M_NAN;
	eval = eval .< 0 .|| eval .> 1 .? M_NAN .: eval;
	
	beta = beta[][ : columns(mDy) - 1]; /* remove vectors with zero eigvals */
    eval = eval[][ : columns(mDy) - 1];         /* also remove zero eigvals */
    amAlpha[0] = s01 * beta;
	amBeta[0] = beta;
	amEval[0] = eval;

	return -ct * reversec(cumulate(log(1-reversec(eval'))))';
}
PcFiml::Cointegration()
{
    decl  y_1, dy, eval, test, icase, gammapar, pval, i, j;

    m_fCointegration = FALSE;
    if (m_iMethod != M_OLS || !fCointData(&dy, &y_1))
	{	m_mEval = zeros(5, m_cY);
        return 0;
	}
						  /* NB: w already concentrated out in fCointData() */
	test = olsRank(dy, y_1, 0, &m_mEval, &m_mBeta, &m_mAlpha);
	eval = m_mEval;
	icase = GetStatus(&i, &j);

	gammapar = GammaParTraceAsymp(icase, range(m_cY,1)', 0);
	pval = 1 - probgamma(test', gammapar[][0], gammapar[][1]);
	m_mEval |= test | pval';

	test = -m_cT * log(1-eval);
	gammapar = GammaParMaxEvAsymp(icase, range(m_cY,1)', 0);
	pval = 1 - probgamma(test', gammapar[][0], gammapar[][1]);
	m_mEval |= test | pval';

    m_fCointegration = TRUE;

    if (m_fPrint)
        Output(FALSE, TRUE);
    return 1;
}

PcFiml::iCointDataI2(const amDdy, const amY_1, const amDy_1, const amDdy_1,
	const amX, const amDx, const amU)
{
    decl  i, w, dw, y_1, y_2, dy_1, dy_i, dy_i_1, dy, maxylag;

    maxylag = GetMaxGroupLag(Y_VAR);

    y_1 = GetGroupLag(Y_VAR, 1, 1);
    y_2 = GetGroupLag(Y_VAR, 2, 2);

    if (m_cY != columns(y_1) || m_cY != columns(y_2))
    {   eprint("Variables missing at lag 1 or 2\n");
        return 0;
    }
    dy  = m_mY - y_1;
	dy_1 = y_1 - y_2;

	amDdy[0] = dy - dy_1;
    amY_1[0] = y_1;
    amDy_1[0] = dy_1;

	w = new matrix[0][0];
    for (i = 1, dy_i = dy_1; i < maxylag - 1; i++, dy_i = dy_i_1)
    {
        y_1 = GetGroupLag(Y_VAR, i + 1, i + 1);
        y_2 = GetGroupLag(Y_VAR, i + 2, i + 2);

        if (m_cY != columns(y_1) || m_cY != columns(y_2))
        {   eprint("Variables missing at lag ", i + 1, " or ", i + 2,"\n");
            return 0;
        }
		dy_i_1 = y_1 - y_2;
        w ~= dy_i - dy_i_1;
    }
    amDdy_1[0] = w;

	decl icase, iconst, itrend, ixu, ixr;
	icase = GetStatus(&iconst, &itrend);
	if (icase < 1)
    {   eprint("Need Constant for I(2) analysis\n");
        return 0;
    }
	if (icase > 3)
    {   eprint("Need restricted Trend for I(2) analysis\n");
        return 0;
    }
	for (i = 0, ixu = -1; i < columns(m_mSelgroup) && i <= iconst; ++i)
		if (m_mSelgroup[0][i] == U_VAR)
			++ixu;
	if (icase == 1)
		itrend = iconst;
	for (i = 0, ixr = -1; i < columns(m_mSelgroup) && i <= itrend; ++i)
		if (m_mSelgroup[0][i] == X_VAR)
			++ixr;
							    /* get restricted vars and their difference */
	amX[0] = amDx[0] = new matrix[0][0];
	--m_iT1sel;				  /* get one extra observation for differencing */
    w = GetGroup(X_VAR);
	++m_iT1sel;
	if (ixr >= 0)			/* drop the trend (constant for case 1), if any */
	{	w[0][ixr] = M_NAN;
		w = deletec(w);
		if (w == 0)
			w = new matrix[0][0];
	}
	if (columns(w))
	{	amX[0] = w[1:][]; 
		dw = diff0(w, 1)[1:][];
		dw[0][] = (w[0][] .== m_dMisVal .|| isdotnan(w[0][])) .? 0 .: dw[0][];
		amDx[0] = dw;
	}
	
	w = GetGroup(U_VAR);
	if (ixu >= 0)							   /* drop the constant, if any */
	{	w[0][ixu] = M_NAN;
		w = deletec(w);
		if (w == 0)
			w = new matrix[0][0];
	}
    amU[0] = w;

return icase;
}
PcFiml::CointegrationI2()
{
    decl  i, ddy, y_1, dy_1, ddy_1, x, dx, u, dy, eval, beta, alpha,
		b, b_orth, a_orth, b_orth_bar, xi, eta, test0, test, tests,
		t, dt, oner, icase, gammapar, pval;

	icase = iCointDataI2(&ddy, &y_1, &dy_1, &ddy_1, &x, &dx, &u);
    if (m_iMethod != M_OLS || icase == 0)
        return 0;
	// case		step 1: DDy = DDy_1 + 	step 2: DDy*a_o = DDy_1 + 		 example
	//  3: 2* 	Pi [y_1 t] + [Dy_1 1]	G [Dy_1*b_ob 1]  + [Dy_1 1]*b	 Jorgensen et.al (96): ukm1
	//  2: 1 	Pi [y_1]   + [Dy_1 1]	G [Dy_1*b_ob 1]  + [Dy_1]*b		 Johansen (95): Australian data
	//  1: 1* 	Pi [y_1 1] + [Dy_1]		G [Dy_1*b_ob]    + [Dy_1]*b	     Paruolo (96): ukm1
	// not done:
	//  4: 2 	Pi [y_1] + [Dy_1 1 t]	G [Dy_1*b_ob 1 t]+ [Dy_1]*b
	//  0: 0 	Pi [y_1] + [Dy_1]		G [Dy_1*b_ob]    + [Dy_1]*b
	// case		formulation: constant  trend   example
	//  3: 2* 	             U_VAR     X_VAR   Jorgensen et.al (96): ukm1
	//  2: 1 	             U_VAR     -----   Johansen (95): Australian data, ukm1:_1Q
	//  1: 1* 	             X_VAR     -----   Paruolo (96): ukm1:_0Q

	t = icase == 3 ? range(1, m_cT)' : icase == 1 ? 1 : new matrix[0][0];
	dt = icase == 1 ? new matrix[0][0] : 1;
	oner = icase == 3 ? 1 : new matrix[0][0];
	
	test0 = olsRank(ddy, y_1 ~ t ~ x, dy_1 ~ ddy_1 ~ dt ~ u, &eval, &beta, &alpha);
	if (icase == 1)
		beta = beta[:m_cY-1][];// restricted constant drops out after differencing

	gammapar = GammaParI2Asymp(icase, 1, range(m_cY,1)', 0);
	pval = 1 - probgamma(test0', gammapar[][0], gammapar[][1])';
	
    if (m_fPrint)
	{	println("\nI(2) analysis");
		println("s is # I(1), n-p-s # I(2) components, p is rank(Pi)");
		println("Q_p:   H0 = rank(Pi)<=p | rank(Pi)<=n");
		println("S_p,s: H0 = rank(Pi)<=p and n-p-s I(2) components | rank(Pi)<=n");
		println("Listed are eigenvalues, tests and p-values.");
		print("%r", {"      p="}, "%10.0f", range(0, m_cY - 1));
		print("%r", {"n-p-s=0 ","Q_p"}, "%10.5g", eval | test0);
		print("%r", {"[pval]  "}, "%10.4f", pval);
//		print("%r", {"5% crit "}, "%10.4f", quangamma(0.95, gammapar[][0]', gammapar[][1]'));
	}

	for (i = 0; i < m_cY; ++i)
	{
		if (i == 0)
			b = 0, b_orth_bar = a_orth = unit(m_cY);
		else
		{	b = beta[][:i - 1];
			b_orth = nullspace(b[:m_cY-1][]);
			a_orth = nullspace(alpha[][:i - 1]);
			b_orth_bar = b_orth * invert(b_orth'b_orth);
		}
		test = olsRank(ddy * a_orth, dy_1 * b_orth_bar ~ dt ~ dx,
			ddy_1 ~ (dy_1 ~ oner ~ dx) * b ~ u, &eval, &eta, &xi);
		tests = test + test0[i];

		gammapar = GammaParI2Asymp(icase, 1, range(0, m_cY-i-1)', range(m_cY-i, 1)');
		pval = 1 - probgamma(tests', gammapar[][0], gammapar[][1])';
		if (i > 0)
		{
			pval = zeros(1, i) ~ pval; 
			eval = zeros(1, i) ~ eval; 
			test = zeros(1, i) ~ test; 
			tests = zeros(1, i) ~ tests;
		}		
    	if (m_fPrint)
		{	if (i == 0)
				print("%r", {"  n-p-s="}, "%10.0f", range(m_cY,1));
			print("%r", {sprint("p = ",i),"S_p,s   "},
				"%10.5g", eval | tests);
			print("%r", {"[pval]  "}, "%10.4f", pval);
//			print("%r", {"5% crit "}, "%10.4f", zeros(1, i) ~ quangamma(0.95, gammapar[][0]', gammapar[][1]'));

		}
	}
//	still need to return or store everything for Monte Carlo;
}

PcFiml::SetPrintUrf(fPrintUrf)
{
    m_fPrintUrf = fPrintUrf;
}

PcFiml::GetVarNames(const aasY, const aasW)
{
    decl csel, cy, cyl, asy, asx, asw;

	GetGroupLagNames(Y_VAR, 0, 0, &asy);	              /* dependent vars */
	cy = sizeof(asy);
	GetGroupLagNames(Y_VAR, 1, 1000, &asw);				 /* lagged dep.vars */
	cyl = sizeof(asw);
	GetGroupNames(X_VAR, &asx);
	asw ~= asx;
	GetGroupNames(U_VAR, &asx);						   /* unrestricted vars */
	asw ~= asx;
	csel = sizeof(asy) + sizeof(asw);

	aasY[0] = asy;
	aasW[0] = asw;

	return cy ~ cyl ~ csel;
}
PcFiml::Print()
{
	if (rows(m_mModsel) != 0)
		Output(FALSE, FALSE);
	else
		Output(TRUE, FALSE);
}

PcFiml::IsUnivariate()
{
	return FALSE;
}
PcFiml::GetPackageName()
{
	return "PcFiml";
}
PcFiml::GetPackageVersion()
{
	return "1.00";
}
PcFiml::GetParNames()
{
    decl cy, cw, asy, asw, as, i, j, k;

	GetVarNames(&asy, &asw);
	cy = sizeof(asy);
	cw = sizeof(asw);

	as = new array[cy * cw];
	if (cy > 1)
	{	for (i = k = 0; i < cy; ++i)
			for (j = 0; j < cw; ++j, ++k)
				as[k] = sprint(asw[j], "@", asy[i]);
	}
	else
		as = asw;

	if (m_iMethod >= M_FIML)
		return as[m_mVarsel];
	else
		return as;
}

PcFiml::Output(const fSys, const fCoint)
{
    decl i, j, k, csel, cy, cyl, asy, asw;

	cy = GetVarNames(&asy, &asw);
	csel = int(cy[0][2]);
	cyl  = int(cy[0][1]);
	cy   = int(cy[0][0]);
	
    if (fSys)
    {
		decl bisgive = m_iMethod == M_GIVE ||
			(m_iMethod == M_NONE && any(m_mSelgroup .== A_VAR));
			
		print("\n\n---- System estimation by ", bisgive ? "IVE" : "OLS", " ----\n");
        if (m_iModelStatus == MS_ESTIMATED)
            println("The estimation sample is ",
                ObsYear(m_iT1est), " (", ObsPeriod(m_iT1est), ") ",
                ObsYear(m_iT2est), " (", ObsPeriod(m_iT2est), ")");
        else if (m_iT1sel >= 0)
            println("The available estimation sample is ", GetSelSample(),
                " (T=", m_iT2sel - m_iT1sel + 1, ")");
        else
            print("No sample selected yet.\n");

        printVariables(Y_VAR, "Y");
        printVariables(X_VAR, 0);
        printVariables(U_VAR, "U");
		if (m_iMethod == M_GIVE)
			print("Additional instruments:", m_asA);
		else
	        printVariables(A_VAR, "A");

		if (m_iMethod == M_GIVE)
			asw = asy[1:] ~ asw;
        if (m_iMethod == M_OLS || m_iMethod == M_GIVE)
        {   print("\ncoefficients", "%r", asw, "%c", asy, m_mPiUrf');
            print("\ncoefficient standard errors", "%r", asw, "%c", asy, sqrt(m_mVarPi'),
                  "\nequation standard errors", "%c", asy, sqrt(diagonal(m_mOmega)),
                  "\nresidual covariance",  "%r", asy,  "%c", asy, m_mOmega);
            print("\nlog-likelihood=", "%.9g", m_dLogLikUrf,
                  " det-omega=", exp(-2 / m_cT * m_dLogLikUrf), " T=", m_cT,
                  "\n");
        }
    }
    else if (rows(m_mModsel) != 0)
    {
        decl asyw = asy ~ asw, asmodels = {"?", "OLS", "IVE", "FIML", "3SLS", "2SLS" };
		
        print("\n\n---- Model estimation by ", asmodels[m_iMethod], " ----\n");
        if (m_iT1sel >= 0)
            print("The estimation sample is ",
                ObsYear(m_iT1est), " (", ObsPeriod(m_iT1est), ") ",
                ObsYear(m_iT2est), " (", ObsPeriod(m_iT2est), ")\n");
        else
            print("No sample selected yet.\n");

        if (m_iMethod <= M_OLS)
        {   setModelIndices();
            print("\nmodel formulation, -1: dependent variable, 1: coefficient ",
                "%r", asyw, "%c", asy, m_mIndex');
        }
        else
        {   decl va, mv;
            va = zeros(1, csel * cy);
            va[0][m_vSelect] = diagonal(m_mVar);
            mv = shape(va, csel, cy);

            print("\ncoefficients", "%r", asyw, "%c", asy, -m_mA',
                  "\ncoefficient standard errors", "%r", asyw, "%c", asy, sqrt(mv),
                  "\nequation standard errors", "%c", asy, sqrt(diagonal(m_mSigma)),
                  "\nresidual covariance",  "%r", asy,  "%c", asy, m_mSigma,
                  "\nlog-likelihood=", "%.9g", m_dLogLik,
                  " det-omega=", exp(-2 / m_cT * m_dLogLik), " T=", m_cT,
                  "\n");

            if (m_iMethod == M_FIML)
                print("FIML estimation: ", MaxConvergenceMsg(m_iResult), "\n");

            if (m_fPrintUrf)
                print("\nRRF coefficients", "%r", asw, "%c", asy, m_mPi',
                      "\nRRF coefficient standard errors", "%r", asw, "%c", asy, sqrt(m_mVarPi'),
                      "\nRRF equation standard errors", "%c", asy, sqrt(diagonal(m_mOmega)),
                      "\n");
        }
    }
    if (fCoint)
    {
		if (csel-cy-1 >= cyl)
			asy ~= asw[cyl:csel-cy-1];	   /* append restricted vars to asy */

		print("\nCointegration analysis\n",
              "%c", {"eigenvalues", "trace", "[pval]", "max-eval", "[pval]"},
			  "%cf", {"%#12.5g", "%#12.5g", "%10.4f", "%#12.5g", "%10.4f"}, m_mEval',
              "\nbeta", "%r", asy, m_mBeta,
              "\nalpha", "%r", asy, m_mAlpha,
              "\nstandardized beta", "%r", asy, m_mBeta ./ diagonal(m_mBeta),
              "\nstandardized alpha", "%r", asy, m_mAlpha .* diagonal(m_mBeta),
              "\nlong run matrix", "%r", asy, "%c", asy, m_mAlpha * m_mBeta'
              );
		decl iconst, itrend;
		decl icase = GetStatus(&i, &j);
		decl ascases = 
			{ "No constant", "Restricted constant", "Unrestricted constant", 
			  "Restricted trend, unrestricted constant","Unrestricted trend and constant"};
		print(ascases[icase], "\n");
    }
}
PcFiml::Chow(const iYear, const iPeriod)
{
    decl t2, ct1, ct, mw, my, mpi, df1, df2, chowtst, ftest, pval, chpval;
	decl rss, rssres, ldet, ldetres, sign, year1, period1, year2, period2;

    if (m_iMethod == M_NONE || m_iMethod == M_GIVE)
        return 0;

    t2 = GetIndex(iYear, iPeriod);               /* database index of break */
    ct1 = t2 - m_iT1est + 1;             /* no of observations in subsample */
           /* ct1 is index in m_mY and m_mW of first post-break observation */
	ct = m_cT;
	if (ct1 <= m_cW)                             /* not enough observations */
    {   eprint("Chow(): invalid break point\n");
        return 0;
    }
	if (ct1 < ct)					  /* break before end of current sample */
    {
		my = m_mY[:ct1-1][];
	    mw = m_mW[:ct1-1][];
		olsc(my, mw, &mpi);
		my -= mw * mpi;						                   /* residuals */
    	rssres = my'my;                                              /* V'V */
		ldetres = logdet(rssres / ct, &sign);

		rss = m_mOmega * (ct - m_cW);
		ldet = -2 * m_dLogLik / ct;

		year1 = iYear;	period1 = iPeriod;
		year2 = ObsYear(m_iT2sel);	period2 = ObsPeriod(m_iT2sel);
	}
	else						 /* break on or after end of current sample */
	{	decl t1sel_save, t2sel_save;
		t1sel_save = m_iT1sel,  t2sel_save = m_iT2sel;
    	ct1 = m_cT;      /* no of observations in subsample (current model) */

		year1 = ObsYear(m_iT1sel);	period1 = ObsPeriod(m_iT1sel); 								                
		/* select extended period: don't use Modelbase version */
		Database::SetSelSample(year1, period1, iYear, iPeriod);
		ct = m_iT2sel - m_iT1sel + 1;	         /* size of extended period */

		if (ct != GetIndex(iYear, iPeriod) - GetIndex(year1, period1) + 1) 
    	{
			eprint("Chow(): invalid forecast point\n");
			m_iT1sel = t1sel_save;  m_iT2sel = t2sel_save;  /* reset sample */
	        return 0;
    	}
    	my = GetGroupLag(Y_VAR, 0, 0);
	    mw = GetGroupLag(Y_VAR, 1, 10000) ~ GetGroup(X_VAR) ~ GetGroup(U_VAR);

		ols2c(my, mw, &mpi);
		my -= mw * mpi;						                   /* residuals */
		rss = my'my;
		ldet = logdet(rss / ct, &sign);					      /* log|V'V|/T */

		rssres = m_mOmega * (ct1 - m_cW);
		//ldetres = logdet(rssres / ct, &sign);
		ldetres = (-2 * m_dLogLik / ct1 + m_cY * log(ct1)) - m_cY * log(ct);

		m_iT1sel = t1sel_save;  m_iT2sel = t2sel_save;      /* reset sample */
		year1 = ObsYear(m_iT2sel);	period1 = ObsPeriod(m_iT2sel);
		year2 = iYear;	period2 = iPeriod;
	}
    if (sign != 1)
        eprint("Chow(): problems with det|V\'V|\n");

	if (ct == ct1)
	{	chowtst = zeros(1, m_cY); ftest = 0;
		GetFapproximation(m_cW, ct, m_cY, ct - ct1, &ftest, &df1, &df2);
		ftest = 0.0;  pval = 1;
		chpval = ones(1, m_cY);
	}
	else
	{	chowtst = diagonal(rss - rssres) ./ diagonal(rssres)
        	   * (ct1 - m_cW) / (ct - ct1);
	    ftest = 1 - exp(ldetres - ldet);
		GetFapproximation(m_cW, ct, m_cY, ct - ct1, &ftest, &df1, &df2);
		pval = tailf(ftest, df1, df2);
		chpval = tailf(chowtst, m_cT - ct1, ct1 - m_cW);
	}

    if (m_fPrint)
    {
        print("\nChow test for break after ",	
			year1, " (", period1, ") in sample up to ",
			year2, " (", period2, "):\n");
		if (m_cY > 1)
	        print("Scalar Chow tests:  F(", (ct - ct1), ",", (ct1 - m_cW), ")=",
				chowtst);
        println(m_cY > 1 ? "Vector " : "",
			"Chow test:   F(", df1, ",", df2, ")=", ftest,
              " [", "%.4f", pval, "]");
    }

return ftest ~ chowtst | pval ~ chpval;
}
PcFiml::ArTest(const iAr1, const iAr2)
{
    decl i, car, w, y, mpi, ftest, test, rss, df1, df2, df, sign;

    if (m_iMethod == M_NONE || m_iMethod == M_GIVE)
        return 0;

    car = iAr2 - iAr1 + 1;
    for (w = m_mW, i = iAr1; i <= iAr2; i++)     /* append lagged residuals */
        w ~= lag0(m_mResidual, i);             /* missing ones are set to 0 */

	olsc(m_mY, w, &mpi);
	w = m_mY - w * mpi; 					                   /* residuals */
    rss = w'w;                                                       /* V'V */

	ftest = 1 - exp(logdet(rss / m_cT, &sign) + 2 * m_dLogLikUrf / m_cT);
    if (sign != 1)
        eprint("ArTest(): problems with det|V\'V|\n");
    GetFapproximation(m_cW, m_cT, m_cY, car * m_cY, &ftest, &df1, &df2);

    test = m_cT * m_cY * (1 - trace(rss * m_mVtVinvUrf) / m_cY );
    df = car * m_cY * m_cY;

    if (m_fPrint)
    {
        println(m_cY > 1 ? "Vector " : "",
		      "AR ", iAr1, "-", iAr2,
              " test: F(", df1, ",", df2, ")=", ftest,
              " [", "%.4f", tailf(ftest, df1, df2), "]",
              " Chi(", df, ")=", test,
              " [", "%.4f", tailchi(test, df), "]"
              );
    }
return ftest ~ test | tailf(ftest, df1, df2) ~ tailchi(test, df);
}
PcFiml::EgeArTest(const iAr1, const iAr2)
{
    decl r, u, mu, mvtvinv0, mvtv, dloglik0, mbinv;
    decl maold, msigmaold;
    decl i, cp, car, cx, ftest, test, rss, df1, df2, df;

    if (m_iMethod != M_FIML)                /* not yet working for 2SLS/3SLS */
        return 0;

    car = iAr2 - iAr1 + 1;
    cp = car * m_cY;
    cx = m_cY + m_cW;

    u = (m_mY ~ m_mW) * m_mA';              /* compute structural residuals */
    r = new matrix[1][0];
    for (i = iAr1; i <= iAr2; i++)           /* form lagged residual matrix */
        r ~= lag0(u, i);                       /* missing ones are set to 0 */

    mvtvinv0 = m_mVtVinv;                                 /* restricted V'V */
    dloglik0 = m_dLogLik;                              /* restricted loglik */
    maold = m_mA;  msigmaold = m_mSigma;                /* save old A,Sigma */

    m_mW ~= r;                           /* append lagged residuals to m_mW */
    m_cW += cp;  m_cU += cp;                      /* these are unrestricted */

    estimateEGE(m_iMethod, FALSE, FALSE);                 /* re-estimate EGE */

    ftest = 1 - exp(2 / m_cT * (dloglik0 - m_dLogLik));
    GetFapproximation(m_cKbar, m_cT, m_cY, cp, &ftest, &df1, &df2);

    mbinv = 1 / (m_mA[][ :m_cY-1]);                               /* inv(B) */
    mvtv = m_cT * mbinv * m_mSigma * mbinv';                       /* V'V/T */
    test = m_cT * m_cY * (1 - trace(mvtv * mvtvinv0) / m_cY );
    df = car * m_cY * m_cY;

    m_cW -= cp;  m_cU -= cp;
    m_mW = m_mW[][ :m_cW-1];
    u = m_mY ~ m_mW;  m_mXtX = u'u;
    m_mA = maold;
    m_mSigma = msigmaold;  m_dLogLik = dloglik0;      /* restore A,Sigma,ll */

    if (m_fPrint)
    {
        println("Vector EGE-AR ", iAr1, "-", iAr2,
              " test: F(", df1, ",", df2, ")=", ftest,
              " [", "%.4f", tailf(ftest, df1, df2), "]",
              " Chi(", df, ")=", test,
              " [", "%.4f", tailchi(test, df), "]"
              );
    }
return ftest ~ test | tailf(ftest, df1, df2) ~ tailchi(test, df);
}
PcFiml::Portmanteau(const iLag)
{
    if (m_iMethod == M_NONE)
        return 0;

return ::PortmanteauTest(m_mResidual, iLag, (m_iMethod >= M_FIML)
    ? columns(m_vSelect) - m_cU * m_cY : (m_cW - m_cU) * m_cY, m_fPrint);
}
PcFiml::HeteroTest(const fStand, const fCross)
{
    decl i, cp, w, cw, r, res, cr, mu, mw, mrtr, sign;
    decl ftest, test, rss, df, df1, df2, logdety;

    if (m_iMethod == M_NONE || m_cW-m_cU <= 0)
        return 0;

    if (fStand)
        r = IN01Data(m_mResidual) .^ 2;    /* squared transformed residuals */
    else
    {   r = m_mResidual .^ 2;                          /* squared residuals */
        for (i = 1; i < m_cY; i++)                     /* add crossproducts */
            r ~= submat(m_mResidual, 0, m_cT, 0, m_cY - 1 - i) .*
                 submat(m_mResidual, 0, m_cT, i, m_cY - 1);
    }
    cr = columns(r);
    r -= sumc(r) / m_cT;                          /* in deviation from mean */
    mrtr = r'r;
    logdety = logdet(mrtr / m_cT, &sign);
    if (sign != 1)
        eprint("HeteroTest(): problems with det|R\'R|\n");

	w = m_mW[][ : m_cW - m_cU - 1];              /* exclude U_VAR variables */
    cw = columns(w);

    w ~= w .* w;                                      /* levels and squares */
    if (fCross)
        for (i = 1; i < cw; i++)                       /* add crossproducts */
            w ~= submat(w, 0, m_cT, 0, cw - 1 - i) .*
                 submat(w, 0, m_cT, i, cw - 1);

    w -= meanc(w);			                      /* in deviation from mean */

    decsvd(w, &mu, &mw);                               /* concentrate out w */
    cp = int(sumr(mw .>= 1e-9*max(sumr(fabs(w)))));/*compute numerical rank */
	mu = mu[][ : cp - 1];						 /* JAD: changed 16-05-2000 */

    r -= mu * (mu'r);                                          /* residuals */
    rss = r'r;                                   /* residual sum of squares */

    ftest = 1 - exp(logdet(rss / m_cT, &sign) - logdety);
          /* constant is only regressors in restricted auxiliary regression */
    if (sign != 1)
        eprint("HeteroTest(): problems with det|V\'V|\n");
//  GetFapproximation(1, m_cT, cr, cp, &ftest, &df1, &df2);
    GetFapproximation(1+m_cW, m_cT, cr, cp, &ftest, &df1, &df2);

    test = m_cT * cr * (1 - trace(rss / mrtr) / cr );
    df = cp * cr;

    if (m_fPrint)
    {
        if (m_cY > 1)
            print("Vector ");
        if (fCross)
            print("hetero-X");
        else
            print("hetero");
        if (fStand)
            print(" (stand)");
        println(" test: F(", df1, ",", df2, ")=", ftest,
              " [", "%.4f", tailf(ftest, df1, df2), "]",
              " Chi(", df, ")=", test,
              " [", "%.4f", tailchi(test, df), "]"
              );
    }
return ftest ~ test | tailf(ftest, df1, df2) ~ tailchi(test, df);
}
PcFiml::NormalityTest()
{
    if (m_iMethod == M_NONE)
        return 0;

return ::NormalityTest(m_mResidual, m_fPrint);
}

PcFiml::UnitRootTest()
{
    decl idx, alpha1, se, test, k, d, i, j, crit1, crit5,
		theta_a1, theta_11, theta_a5, theta_15, sig, star = {""," *"," **"};

    if (m_iMethod != M_OLS || m_cY != 1 || m_cYlag == 0)
        return 0;

	idx = zeros(1, m_cW);
	idx[ : m_cYlag - 1] = 1;
	alpha1 = sumr(m_mPi[ : m_cYlag - 1]) - 1;
	se = sqrt(m_mOmega * (idx * m_mWtWinv * idx'));

    test = double(alpha1 / se);
	d = GetStatus(&i, &j);
	d = (d == 1 || d == 2) ? 1 : (d == 3 || d == 4) ? 2 : 0;
	// d = 0: none, 1: constant, 2: trend
	k = sumr( (m_mSelgroup .== X_VAR .|| m_mSelgroup .== U_VAR) .&& m_mLagsel .== 0);
	k = 1 + int(k) - d;		// k: no of regressors including y

	if (k == 1)
	{	theta_a5 = <-1.95,-2.86,-3.41>[0][d];
		theta_a1 = <-2.57,-3.43,-3.96>[0][d];
		theta_15 = <-0.35,-2.81,-4.38>[0][d];
		theta_11 = <-2.19,-6.52,-8.99>[0][d];
	}
	else
	{	theta_a5 = -2.056 - 0.535 * d * (1-k / 20) - 0.343 * k * (1 - k / 40);
		theta_a1 = theta_a5 - 0.607;
		theta_15 = 7.09 - 0.483 * k + 3.64 * theta_a5 * (1 - k / 10);
		theta_11 = 13.8 - 0.447 * k + 6.06 * theta_a1 * (1 - k / 20);
	}
	crit5 = theta_a5 + theta_15 / (m_cT - (2*k - 1) - d);
	crit1 = theta_a1 + theta_11 / (m_cT - (2*k - 1) - d);

	if (test < crit1)
		sig = 2;
	else if (test < crit5)
		sig = 1;
	else
		sig = 0;
		
	if (m_fPrint)
    {
		if (k == 1)
			println("ADF(", m_cYlag - 1, ") = ", test, star[sig]);
		else
			println("ECM(", k, ") = ", test, star[sig]);
		println("Critical values used in test: 5%=", "%.4g", crit5, ", 1%=", "%.4g", crit1);
		if (d > 0) println("Test based on presence of Constant", d > 1 ? " and Trend" : "");
    }
return test | (sig == 0);
}

PcFiml::GetPi()
{
    return m_mPi;
}
PcFiml::GetVarPi()
{
    return m_mVarPi;
}
PcFiml::GetVarRf()
{
    if (m_iMethod >= M_FIML)	
	{
		decl mbi, mj;
	    mbi = 1 / (m_mA[][ :m_cY-1]);                             /* inv(B) */
    	mj = (-mbi ** (m_mPi'~unit(m_cW)))[][m_vSelect];/*(inv(B)**[Pi':I])u*/
    	return mj * m_mVar * mj';
	}
    else if (m_iMethod != M_NONE)
        return m_mOmega ** m_mWtWinv;

return 0;
}
PcFiml::GetVarTheta()
{
    if (m_iMethod >= M_FIML)
        return m_mVar;
    else if (m_iMethod != M_NONE)
        return m_mOmega ** m_mWtWinv;

return 0;
}
PcFiml::GetOmega()
{
    return m_mOmega;
}
PcFiml::GetResiduals()
{
    return m_mResidual;
}
PcFiml::SetEquation(const sEquation, const aModel)
{
    decl i, j, k, ieq, ivar, csel, carg, cp, cy;

    csel = columns(m_mVarsel);              /* number of selected variables */

    for (i = j = 0, ieq = -1; i < csel; ++i)               /* find equation */
        if (m_mSelgroup[0][i] == Y_VAR && m_mLagsel[0][i] == 0)
        {
            if (sEquation == m_asNames[m_mVarsel[0][i]])
                ieq = j, ivar = i;
            j++;                                         /* no of equations */
        }
    cy = j;

    if (rows(m_mModsel) == 0)
        m_mModsel = zeros(cy, csel);

    if (ieq < 0)
    {   eprint(sEquation, " is not an equation\n");
        return 0;
    }
    for (i = cp = 0; i < csel; ++i)                  /* initialize equation */
        if (m_mSelgroup[0][i] == U_VAR)
            m_mModsel[ieq][i] = 1, ++cp;    /* select unrestricted variable */
        else
            m_mModsel[ieq][i] = 0;                 /* don't select variable */
    m_mModsel[ieq][ivar] = -1;                        /* dependent variable */

    carg = sizeof(aModel);
    for (i = 0; i < carg; i += 3)            /* "name", first lag, last lag */
    {
        for (k = aModel[i + 1]; k <= aModel[i + 2]; k++)
        {   ivar = FindSelection(aModel[i], k);
            if (ivar < 0)
                return 0;
            m_mModsel[ieq][ivar] = 1;
            cp++;
        }
    }
    if (cp > csel - cy)
        print("Order condition for ", sEquation, " not satisfied!");

                       /* m_mModsel: -1 = dep.var 0=not selected 1=selected */
}
PcFiml::fCheckRank()
{
    decl i, j, k, cx = m_cY + m_cW, cex, czero, fsubset, fzero, vp;

    /* m_mIndex has -1 for dep vars,1 for included vars,0 for excluded vars */
                    /* ma[m_cY+cI][cx] = A, we'll only look at ma[m_cY][cx] */

    for (i = 0; i < m_cY; i++)
    {
        for (k = 0; k < m_cY; k++) /* check if equation k is subset of eq i */
        {                          
            if (k == i)                      /* equation is itself --> skip */
                continue;                    
            for (j = 0, fsubset = TRUE; fsubset && j < cx; j++)
                if (m_mIndex[i][j] == 0 && m_mIndex[k][j] != 0)
                    fsubset = FALSE;
            if (fsubset)
            {   print(
                  "\nEq. ", i + 1, " is not identified on the rank condition: "
                    "eq. ", k + 1, " is a subset!\n");
                return FALSE;
            }
        }
        for (j = cex = czero = 0; j < cx; j++)
        {                                  /* check rank of excluded system */
            if (m_mIndex[i][j] != 0)
                continue;                          /* do excluded vars only */
            cex++;                             /* count no of excluded vars */
            for (k = 0, fzero = 1; fzero && k < m_cY; k++)
                if (m_mIndex[k][j] != 0)
                    fzero = 0;
            czero += fzero;          /* count no of columns with all zero's */
        }
        if (cex - czero < m_cY - 1)
        {   print(
                "\nEq.", i + 1, " is not identified on the rank condition: "
                "rank(excluded system)=", cex - czero, "\n");
            return FALSE;
        }
    }

return TRUE;
}
PcFiml::setModelIndices()
{
    decl i, j, cx, va;

    cx = columns(m_mSelgroup);
    va = new matrix[1][cx];

    for (i = j = 0; i < cx; ++i)     /* order the selection, first dep.vars */
        if (m_mSelgroup[0][i] == Y_VAR && m_mLagsel[0][i] == 0)
            va[0][j++] = i;

    for (i = 0; i < cx; ++i)                        /* then lagged dep.vars */
        if (m_mSelgroup[0][i] == Y_VAR && m_mLagsel[0][i] > 0)
            va[0][j++] = i;

    for (i = 0; i < cx; ++i)                                 /* then X vars */
        if (m_mSelgroup[0][i] == X_VAR)
            va[0][j++] = i;

    for (i = 0; i < cx; ++i)                      /* then unrestricted vars */
        if (m_mSelgroup[0][i] == U_VAR)
            va[0][j++] = i;

    m_mIndex = m_mModsel[][va];

    m_vSelect = vecindex( (m_mIndex .> 0)' )';
}
PcFiml::initializeEGE()
{
    decl  va, mx;

    if (m_iMethod == M_NONE)
    {   eprint("No system estimated yet\n");
        return FALSE;
    }
    setModelIndices();
    if (!fCheckRank())
        return FALSE;

    mx = m_mY ~ m_mW;
    m_mXtX = mx'mx;
    m_mA = m_mIndex .< 0;

return TRUE;
}
PcFiml::solveEGE(m_mPi, m_mSigma)
{
    decl m_1, sigmai, vp, va, cx = m_cY + m_cW;

    m_1 = (m_mPi | unit(m_cW)) * m_mXtX[m_cY: ][];                  /* QW'X */
    sigmai = (1 / m_mSigma);
                                                    /* (inv(Sigma)Pi W'X)vu */
    vp = (vec((sigmai * m_1[ :m_cY-1][])'))[m_vSelect][0];
                                              /* inv((inv(Sigma) ** QW'X)u) */
    m_mVar = 1 / ((sigmai ** m_1)[m_vSelect][m_vSelect]);
    va = vec(m_mA');
    va[m_vSelect][0] = -vec(m_mVar' * vp);
    m_mA = shape(va, cx, m_cY)';
    this->m_mSigma = m_mA * (m_mXtX) * m_mA' / m_cT;   /* ML Sigma=AX'XA'/T */
}
PcFiml::accumulateEGE(const iMethod)
{
    decl mbinv, mj, vp, va, cx = m_cY + m_cW, logdetvtv, sign;

    mbinv = 1 / (m_mA[][ :m_cY-1]);                               /* inv(B) */
    m_mOmega = mbinv * m_mSigma * mbinv';                          /* V'V/T */
    m_mVtVinv = invert(m_mOmega, &logdetvtv, &sign) / m_cT;     /* inv(V'V) */
    if (sign != 1)
        eprint("accumulateEGE(): problems with det|V\'V|\n");
    m_dLogLik = -m_cT / 2 * logdetvtv;

    m_mPi = -mbinv * m_mA[][m_cY: ];
    m_cKbar = int(columns(m_vSelect) / m_cY);
    m_mOmega *= m_cT / (m_cT - m_cKbar);                    /* V'V/(T-kbar) */
    this->m_mSigma *= m_cT / (m_cT - m_cKbar);

    m_mResidual = (m_mY ~ m_mW) * m_mA' * mbinv';

    if (iMethod == M_2SLS)              /* 2SLS: (Sigma .* EGE Var(theta))u */
        m_mVar = (this->m_mSigma ** ones(cx,cx))[m_vSelect][m_vSelect] .* m_mVar;
    else if (iMethod == M_FIML)
    {
        decl mq, m_1;
        mq = m_mPi | unit(m_cW);                              /* Q'=[Pi':I] */
        m_1 = mq * m_mXtX[m_cY: ][m_cY: ] * mq';                  /* QW'WQ' */
                                             /* inv((inv(Sigma) ** QW'WQ)u) */
        m_mVar = 1 / ( ((1 / m_mSigma) ** m_1)[m_vSelect][m_vSelect] );
    }
    else
        m_mVar *= m_cT / (m_cT - m_cKbar);

    mj = (-mbinv ** (m_mPi' ~ unit(m_cW)))[][m_vSelect];/*(inv(B)**[Pi':I])u*/
    m_mVarPi = shape(diagonal(mj * m_mVar * mj'), m_cW, m_cY)';

    m_iMethod = iMethod;
}
PcFiml::fFimlLogLik(const vP, const adFunc, const avScore, const amHessian)
{
    decl va, cx = m_cY + m_cW, mbinv, logdetsig, sign;

    va = vec(m_mA');
    va[m_vSelect][0] = -vP;
    m_mA = shape(va, cx-m_cU, m_cY)';
    m_mSigma = m_mA * (m_mXtX) * m_mA' / m_cT;            /* Sigma=AX'XA'/T */
    mbinv = 1 / (m_mA[][ :m_cY-1]);
    logdetsig = logdet(mbinv * m_mSigma * mbinv', &sign);
    if (sign != 1)
        eprint("accumulateEGE(): problems with det|Sigma|\n");
    adFunc[0] = -logdetsig / 2;

    if (avScore)                  /* q = (inv(Sigma) * [-AX'ZPi' : -AX'Z])vu */
    {   avScore[0] = vec(
            ( (1 / m_mSigma) * (m_mA * m_mXtX)[][m_cY: ] *
            ( (-mbinv * m_mA[][m_cY: ])' ~ unit(m_cW-m_cU)) ) ')
            [m_vSelect][0] / m_cT;
    }
return 1;
}
PcFiml::estimateEGE(iMethod, const fInit, const fAccumulate)
{
    m_iResult = 0;

    if (fInit && !initializeEGE())
    {   m_iResult = -1;
        return FALSE;
    }

    if (iMethod == M_2SLS || iMethod == M_3SLS || (iMethod == M_FIML && m_iMethod == M_OLS) )
        solveEGE(m_mPiUrf, unit(m_cY));                    /* do 2SLS first */

    if (iMethod == M_3SLS || (iMethod == M_FIML && m_iMethod == M_OLS) )
        solveEGE(m_mPiUrf, m_mSigma);        /* 3SLS uses Sigma_2sls Pi_urf */

    if (iMethod == M_FIML)
    {
        decl vp, f, dg, mh, mu, mw, u, x, cx = m_cY + m_cW;

        m_vSelect = vecindex( (m_mIndex[][ :cx-m_cU-1] .> 0)' )';
        m_mA = m_mA[][ :cx-m_cU-1];

        if (m_cU)
        {
            x = m_mY ~ m_mW[][ :m_cW-m_cU-1];
            u = m_mW[][m_cW-m_cU: ];
            decsvd(u, &mu, &mw);     /* concentrate out u, assume full rank */
            m_mXtX = x'x - (x'mu)*(mu'x);        /* residual sum of squares */
        }
        vp = -(vec(m_mA'))[m_vSelect][0];
        fFimlLogLik(vp, &f, 0, 0);
        mh = unit(columns(m_vSelect));

        m_iResult = MaxBFGS(fFimlLogLik, &vp, &f, &mh, FALSE);
        if (m_fPrint && m_iResult > 1)
            print("FIML estimation: ", MaxConvergenceMsg(m_iResult), "\n");

        m_dLogLik = f * m_cT;

        if (m_cU)
        {   m_mA ~= -m_mA * ((1 / u) * x)';
            x = m_mY ~ m_mW;  m_mXtX = x'x;
        }
        m_vSelect = vecindex( (m_mIndex .> 0)' )';
    }
    if (fAccumulate && m_iResult <= 1)
        accumulateEGE(iMethod);

return (m_iResult == 0 || m_iResult == 1); /* 0: MAX_CONV, 1: MAX_WEAK_CONV */
}
PcFiml::TwoSLS()
{
    if (!estimateEGE(M_2SLS, TRUE, TRUE))
	    m_iMethod = M_NONE;

    if (m_fPrint)
       Output(FALSE, FALSE);

return TRUE;
}
PcFiml::ThreeSLS()
{
    if (!estimateEGE(M_3SLS, TRUE, TRUE))
	    m_iMethod = M_NONE;

    if (m_fPrint)
       Output(FALSE, FALSE);

return TRUE;
}
PcFiml::Fiml()
{
    decl result = estimateEGE(M_FIML, TRUE, TRUE);
	if (!result)
	    m_iMethod = M_NONE;

    if (m_fPrint)
       Output(FALSE, FALSE);

return result;
}

