/*--------------------------------------------------------------------------
 * pcnaive.ox - code for PcNaiveDgp class
 *
 *       (C) Jurgen Doornik 1994-98
 *
 * 25 -2-1998 (JAD): small adjustments to asymptotic analysis
 * 22 -9-1997 (JAD): storing distribution arguments in rows
 * ..- 5-1997 (JAD): added support for breaks; I(2) analysis in Asymp()
 *--------------------------------------------------------------------------*/

#include <oxstd.h>
#include <oxprob.h>
#include <pcnaive.h>

const decl sPcNaiveVersion = "2.00";

/*
       y_t = A0 y_t + A1 y_t-1 + A2 z_t + a3 + u_t   y_t,a3 (n x 1), z_t (q x 1)
	         + A5 y_t-2
       u_t = B0 u_t-1 + e_t + B1 e_t-1               u_t (n x 1), e_t (n x 1)
       z_t = C0 z_t-1 + c1 + c2 t + v_t  C0 (q x q), c1 (q x 1), c2 (q x 1)
matrix form:
       Y = Y A0' + Y_1 A1' + Z A2' + O a3' + U       Y (T x n), Z (T x q)
       U = U_1 B0' + E + E_1 B1'                     O (T x 1) vector of ones
       Z = Z_1 C0' + O c1' + R c2' + V               R (T x 1) trend
                                                     
reduced form:
       Y  = Y_1 A1'inv(D) + Z A2'inv(D) + O a3'*inv(D) + U*inv(D)   D=(I-A0)'
	   		+ Y_2 A5'inv(D)


       Dy_t = AB'(y_t-1:z_t) + A2 z_t + a3 + u_t   y_t,a3 (n x 1), z_t (q x 1)
	         + A5 Dy_t-1						   A (n x p), B (n x p+q)
	   translated to:
       y_t = (AB'+I+A5) y_t-1 + (AB'+A2) z_t + a3 + u_t + (-A5) y_t-2

Initial conditions (disregarding A0):, write e for static part (A2 z_t + a3 + u_t)
1) y_t = 0, t < 0, can be done using cumulate()
       y_0 =                   e_0
       y_1 = A1 y_0 +          e_1
       y_2 = A1 y_1 + A5 y_0 + e_2

1) y[-2] | y[-1], t < 0, cannot be done using cumulate()
       y_0 = A1 y[-1] + A5 y[-2] + e_0
       y_1 = A1 y_0   + A5 y[-1] + e_1
       y_2 = A1 y_1   + A5 y_0   + e_2
   construct
	   u_0 = A1 y[-1] + A5 y[-2] + e_0
	   u_1 = 	        A5 y[-1] + e_1
	   u_t =			         e_t, t >= 2
   then in cumulate():
   		y_0 =                   u_0 = A1 y[-1] + A5 y[-2] + e_0
   		y_1 = A1 y_0 +          u_1 = A1 y_0   + A5 y[-1] + e_1
   		y_2 = A1 y_1 + A5 y_0 + u_2
*/

PcNaiveDgp::PcNaiveDgp(const cY, const cZ)
{
    m_mB0t = m_mA5 = m_mA5tRf = m_mAlpha = m_mBeta = m_mZC = 0;
	m_bCustomZ = FALSE;

    m_cY = cY;
    if (m_cY <= 0)
        eprint("PcNaiveDgp(): cY must be > 0");

    SetDistribution(U_DGP, MVNORMAL, zeros(cY,1), unit(cY));
    SetUParameter(zeros(cY, cY), zeros(cY, cY));

    m_cZ = cZ;
    if (cZ > 0)
    {   SetDistribution(Z_DGP, MVNORMAL, zeros(cZ,1), unit(cZ));
        SetZParameter(zeros(cZ,cZ), zeros(cZ,1), zeros(cZ,1));
    }
    else
        m_mC0t = m_mC1t = m_mC2t = m_iZdist = m_mZdistM =
            m_mZdistS = 0;
    SetFixedZ(TRUE);

    SetYParameter(zeros(cY,cY), zeros(cY,cY), zeros(cY,cZ), zeros(cY,1));

	SetInit(Y_DGP, 0);
	SetInit(Z_DGP, 0);
}
PcNaiveDgp::SetUParameter(const mLagAr, const mLagMa)
{
    if (rows(mLagAr) != columns(mLagAr) || rows(mLagAr) != m_cY
     || rows(mLagMa) != columns(mLagMa) || rows(mLagMa) != m_cY)
        eprint("U parameter error\n");

    m_bArmaErrors = FALSE;
	
    m_mB0t = mLagAr';
    m_mB1t = mLagMa';

	if (!(m_mB1t == 0))
	    m_bArmaErrors = TRUE;
    if (!(m_mB0t == 0))
	    m_bArmaErrors = TRUE;
}
PcNaiveDgp::SetZParameter(const mLag, const vConst, const vTrend)
{
    if (rows(mLag) != columns(mLag) || rows(mLag) != m_cZ
        || rows(vConst) != m_cZ || rows(vTrend) != m_cZ)
        eprint("Z parameter error\n");

    m_mC0t   = mLag';
    m_mC1t = vConst';
    m_mC2t = vTrend';
}
PcNaiveDgp::SetYParameter(const mLevel, const mLag, const mZpar,
	const vConst, ...)
{
    if (rows(mLevel) != columns(mLevel) || (!(mLevel == 0) && rows(mLevel) != m_cY)
        || rows(mLag) != columns(mLag) || rows(mLag) != m_cY
        || (m_cZ && (rows(mZpar) != m_cY)) || (m_cZ && (columns(mZpar) != m_cZ))
        || rows(vConst) != m_cY || columns(vConst) != 1)
        eprint("Y parameter error\n");

    m_mA0 = setdiagonal(mLevel, 0);
    m_mA1 = mLag;
    m_mA2 = mZpar;
    m_mA3 = vConst;

	m_mA5 = m_mA5tRf = 0;
	decl va = va_arglist();
	if (sizeof(va) >= 1)
	{
		m_mA5 = va[0];
		if (m_mA5 == 0)
			m_mA5 = 0;
	    else if (rows(m_mA5) != columns(m_mA5) || rows(m_mA5) != m_cY)
        {	eprint("Y parameter error: A5\n");
			m_mA5 = 0;
		}
	}

	if (!(m_mA0 == 0))
	{	m_mDinv = 1 / (unit(m_cY) - m_mA0');
		m_bSimultaneous = TRUE;
	}
	else
	{	m_mDinv = 1;
		m_bSimultaneous = FALSE;
		m_mA0 = 0;
	}
    m_mA1tRf = m_mA1' * m_mDinv;
    m_mA2tRf = (m_cZ) ? m_mA2' * m_mDinv : 0;
    m_mA3tRf = m_mA3' * m_mDinv;
	m_mA5tRf =  m_mA5' * m_mDinv;

	m_cYlag = 0;
	if (!(m_mA1tRf == 0))
		m_cYlag = 1;
	if (!(m_mA5tRf == 0))
		m_cYlag = 2;
	m_bUsedEcm = FALSE;
}
PcNaiveDgp::SetYParameterEcm(const mAlpha, const mBeta, const mZpar,
	const vConst, ...)
{
	decl mpi1, mpi2 = 0, ma5 = 0;

	mpi1 = mAlpha * mBeta[:m_cY - 1][]';
	if (rows(mBeta) > m_cY)
		mpi2 = mAlpha * mBeta[m_cY:][]';

	decl va = va_arglist();
	if (sizeof(va) >= 1)
		ma5 = va[0];
		  /* y_t = (AB'+I+A5) y_t-1 + (AB'+A2) z_t + a3 + u_t + (-A5) y_t-2 */
	SetYParameter(0, mpi1 + unit(m_cY) + ma5, mpi2 + mZpar, vConst, -ma5);
	m_bUsedEcm = TRUE;
	m_mAlpha = mAlpha;
	m_mBeta = mBeta;
}
PcNaiveDgp::SetDistribution(const iEqn, const iDist, mPar1, mPar2)
{
    if (iEqn == U_DGP)
    {
		[m_iEdist, m_mEdistM, m_mEdistS] = RanMC::CheckDist("U_DGP", iDist, mPar1, mPar2);
    }
    else if (iEqn == Z_DGP)
    {
		[m_iZdist, m_mZdistM, m_mZdistS] = RanMC::CheckDist("Z_DGP", iDist, mPar1, mPar2);
    }
}
PcNaiveDgp::SetInit(const iDgp, const mInit)
{
	if (iDgp == Y_DGP)
	{
		m_mInitY = mInit;
		if (columns(m_mInitY) == 1)					    /* want a T in rows */
			m_mInitY = m_mInitY';
		m_bInitY = !(m_mInitY == 0);

		if (m_bInitY)
		{	if (rows(m_mInitY) == 1)				    /* allow for 2 lags */
				m_mInitY |= m_mInitY;
		}
		else
			m_mInitY = zeros(2, m_cY);
	}
	else if (iDgp == Z_DGP)
	{
		m_mInitZ = mInit;
		if (columns(m_mInitZ) == 1)					   /* want a row vector */
			m_mInitZ = m_mInitZ';
	}
}

PcNaiveDgp::SetFixedZ(const fSetting)
{
    m_fFixedZ = fSetting;  m_mZ = 0;
}
PcNaiveDgp::DiscardZ()
{
    m_mZ = 0;
}
PcNaiveDgp::GetZ()
{
	return m_mZ;
}
PcNaiveDgp::GetY()
{
    return m_mY;
}
PcNaiveDgp::GetU()
{
    return m_mU;
}
PcNaiveDgp::SetZCustom(mZC)
{
	m_mZC = mZC;
	m_bCustomZ = rows(m_mZC) > 1;
}

PcNaiveDgp::generateZ(const cT)
{
    if (m_cZ && (rows(m_mZ) == 0 || !m_fFixedZ))   /* Z, but only if needed */
    {
		m_mZ = RanMC::RanDist(m_iZdist, cT, m_cZ, m_mZdistM, m_mZdistS);

		if (m_bCustomZ)
			m_mZ += m_mZC;

		m_mZ +=	m_mC1t + range(1,cT)' .* m_mC2t;
		m_mZ[0][] += m_mInitZ * m_mC0t;	 /* + C0 z: initial value for lagged z */
		
	    m_mZ = cumulate(m_mZ, m_mC0t);
	}
}
PcNaiveDgp::generateU(const cT)
{
	if (m_iEdist == MVNHETERO)
        return generateHeteroObs(cT);

	m_mU = RanMC::RanDist(m_iEdist, cT, m_cY, m_mEdistM, m_mEdistS);

	if (m_bArmaErrors)
	{	if (!(m_mB1t == 0))
	        m_mU = cumulate(m_mU + (zeros(1,m_cY) | m_mU[0:cT-2][]) * m_mB1t, m_mB0t);
    	else if (!(m_mB0t == 0))
        	m_mU = cumulate(m_mU, m_mB0t);
	}
}

PcNaiveDgp::Generate(const cT)
{
	PcNaiveDgp::generateZ(cT);								 /* Z if needed */
	PcNaiveDgp::generateU(cT);				    /* generate U, errors for Y */

	m_mY = m_mZ * m_mA2tRf + m_mA3tRf + m_mU * m_mDinv;		/* 'baseline' Y */

	if (m_bInitY)						     /* start up, up to 2 lags in y */
	{										  /* see comment at top of file */
		m_mY[0:1][] += m_mInitY * m_mA5tRf;
		m_mY[0][] += m_mInitY[1][] * m_mA1tRf;
	}
													  
	m_mY = cumulate(m_mY, m_mA1tRf,	m_mA5tRf);				  /* generate Y */
	return m_mY;
}

PcNaiveDgp::GenerateTo(const cT)
{
	PcNaiveDgp::Generate(cT);
	return m_mY ~ m_mU;
}

PcNaiveDgp::generateY(const iT1, const iT2, const mU, const mInitY,
	const mA1tRf, const mA2tRf, const mA3tRf, const mA5tRf, const mDinv)
{
	decl u, ymin2, ymin1;
	
	u = mA3tRf + mU[iT1 : iT2][] * mDinv;
	if (m_cZ)
		u += m_mZ[iT1 : iT2][] * mA2tRf;

	ymin2 = mInitY[0][];
	ymin1 = mInitY[1][];
	u[0][] += ymin1 * mA1tRf + ymin2 * mA5tRf;
	u[1][] += ymin1 * mA5tRf;
	/* Similar as in comment at top of file. In cumulate():
       		y_0 = u_0*          = A1 y[-1] + A5 y[-2] + u_0
   	   		y_1 = A1 y_0 + u_1* = A1 y_0   + A5 y[-1] + u_1
       		y_2 = A1 y_1 + A5 y_0 + u_2
	*/
    return cumulate(u, mA1tRf, mA5tRf);
}

PcNaiveDgp::GenerateBreakTo(const cT, const iTbreak,
	const iTreset, const mA0, const mA1, const mA2, const mA3, const mA4)
{
	decl y0;
								 				    /* generate all Z and U */
	if (m_iEdist == MVNHETERO)
		eprint("Breaks with MVNHETERO not supported");

	PcNaiveDgp::generateZ(cT);								 /* Z if needed */
	PcNaiveDgp::generateU(cT);				    /* generate U, errors for Y */
	m_mY = zeros(cT, m_cY);		           /* create matrix of correct size */

	y0 = m_mInitY;				 /* initial values allowing for 2 lags of y */

	if (iTbreak > 0)					  /* generate over pre-break sample */
	{	m_mY[][] = PcNaiveDgp::generateY(0, min(cT, iTbreak) - 1,
			m_mU, y0, m_mA1tRf, m_mA2tRf, m_mA3tRf, m_mA5tRf, m_mDinv);

		if (iTbreak == 1)
			y0[1][] = m_mY[iTbreak-1 : iTbreak-1][];   /* y0[0][] unchanged */
		else
			y0 = m_mY[iTbreak-2 : iTbreak-1][];
	}

	if (iTbreak < cT && iTbreak < iTreset)   /* generate over broken sample */
	{
		decl mdinv, ma1trf, ma2trf, ma3trf, ma5trf;
		
		mdinv  = 1 / (unit(m_cY) - mA0');
    	ma1trf = mA1' * mdinv;
	    ma2trf = (m_cZ) ? mA2' * mdinv : 0;
    	ma3trf = mA3' * mdinv;
	    ma5trf = m_mA5' * mdinv;
		mdinv = mA4 * mdinv;      /* this will multiply u in broken process */
		
		m_mY[iTbreak:][] = PcNaiveDgp::generateY(iTbreak, min(cT, iTreset) - 1,
			m_mU, y0, ma1trf, ma2trf, ma3trf, ma5trf, mdinv);
	}
	if (iTreset < cT)
	{	y0 = m_mY[iTreset-2 : iTreset-1][];
		m_mY[iTreset:][] = PcNaiveDgp::generateY(iTreset, cT - 1,
			m_mU, y0, m_mA1tRf, m_mA2tRf, m_mA3tRf, m_mA5tRf, m_mDinv);
	}
		
	return m_mY ~ m_mU;
}
PcNaiveDgp::generateHeteroObs(const cT)
{
    decl  i, trend, m0 = zeros(1, m_cY);
                                      /* generate one observation at a time */
    m_mU = new matrix[cT][m_cY];
    m_mY = new matrix[cT][m_cY];

    m_mU[0][] = RanMC::RanDist1(MVNHETERO, m_cY, m0, m_mEdistS, 0, m0);
    m_mY[0][] = m_mZ[0][] * m_mA2tRf + m_mA3tRf + m_mU[0][] * m_mDinv;

    for (i = 1; i < cT; i++)
    {
        m_mU[i][] = RanMC::RanDist1(MVNHETERO, m_cY, m_mEdistM, m_mEdistS,
			0, m_mY[i-1][]);
			
        if (!(m_mB1t == 0))
            m_mU[i][] += m_mU[i-1][] * m_mB1t;
        else if (!(m_mB0t == 0))
            m_mU[i][] += m_mU[i-1][] * m_mB0t;

		if (m_cZ)
	        m_mY[i][] = m_mZ[i][] * m_mA2tRf + m_mA3tRf + m_mU[i][] * m_mDinv
    	         + m_mY[i-1][] * m_mA1tRf + i > 1 ? m_mY[i-2][] * m_mA5tRf : 0;
		else
	        m_mY[i][] = m_mA3tRf + m_mU[i][] * m_mDinv
				+ m_mY[i-1][] * m_mA1tRf + i > 1 ? m_mY[i-2][] * m_mA5tRf : 0;
    }
return m_mY;
}
PcNaiveDgp::Print()
{
    decl i, fyparz = FALSE;

    print("\n---- PcNaive (", sPcNaiveVersion, ") DGP ----\ny is (",
		m_cY, " x 1), z is (", m_cZ, " x 1) and",
		(m_fFixedZ) ? " " : " NOT ", "fixed.\n");

    if (!(m_mB0t == 0) || !(m_mB1t == 0))
        print("\ny[t] = u[t]");
    else
        print("\ny[t] = e[t]");
    if (!(m_mA0 == 0))
        print(" + A0 y[t]");
    if (!(m_mA1 == 0))
        print(" + A1 y[t-1]");
    if (!(m_mA2 == 0))
        print(" + A2 z[t]"), fyparz = TRUE;
    if (!(m_mA3 == 0))
        print(" + a3");
    if (!(m_mA5 == 0))
        print(" + A5 y[t-2]");
    print("\n");

    if (!(m_mA0 == 0))
        print("A0 = ", m_mA0);
    if (!(m_mA1 == 0))
        print("A1 = ", m_mA1);
    if (fyparz != 0)
        print("A2 = ", m_mA2);
    if (!(m_mA3 == 0))
        print("a3 = ", m_mA3);
    if (!(m_mA5 == 0))
        print("A5 = ", m_mA5);
	if (m_bInitY)
		print("Y initial values:", "%r", {"y[-2]", "y[-1]"}, m_mInitY);
    print("\n");

    if (!(m_mB0t == 0) || !(m_mB1t == 0))
    {
        print("\nu[t] = B0 u[t-1] + e[t] + B1 e[t-1]\n");
        print("B0 = ", !(m_mB0t == 0) ? m_mB0t' : "0 ");
        print("B1 = ", !(m_mB1t == 0) ? m_mB1t' : "0 ", "\n");
    }
    RanMC::WriteDist("e", m_iEdist, m_mEdistM, m_mEdistS);

    if (!(m_mA0 == 0))
    {   print("\nreduced form y[t] = P0 y[t-1] + P1 z[t] + p2 + P3 y[t-2] + w[t]\n");
        print("P0 = ", m_mA1tRf', "P1 = ", !(m_mA2tRf == 0) ? m_mA2tRf' : "0 ",
			"p2 = ", m_mA3tRf', "P3 = ", !(m_mA5tRf == 0) ? m_mA5tRf' : "0\n");
	    if (m_iEdist == MVNORMAL)
   	        print("V[w]=", m_mDinv'm_mEdistS'm_mEdistS*m_mDinv);
	}

    if (m_cZ == 0 || !fyparz)
        return;

    print("\nz[t] = v[t]");
    if (!(m_mC0t == 0))
        print(" + C0 z[t-1]");
    if (!(m_mC1t == 0))
        print(" + c1");
    if (!(m_mC2t == 0))
        print(" + c2 t");
	if (m_bCustomZ)
        print(" + custom part");
    print("\n");

    if (!(m_mC0t == 0))
        print("C0 = ", m_mC0t');
    if (!(m_mC1t == 0))
        print("c1 = ", m_mC1t');
    if (!(m_mC2t == 0))
        print("c2 = ", m_mC2t');
	if (!(m_mInitZ == 0))
		print("Z initial values:", m_mInitZ);
    print("\n");
    RanMC::WriteDist("v", m_iZdist, m_mZdistM, m_mZdistS);
}

PcNaiveDgp::AsympEx(const asModelY, const asModelX, const cT)
{
	decl i, clag, asy, md, ms, ev, bstationary, mm, mmy, mmz, mpi, mmeanlag, mgamma;
	decl bdoy_2 = !(m_mA5 == 0 && m_mB0t == 0);

	[bstationary, md, ms, ev, mm, mmy, mmz, mpi, mmeanlag, mgamma] = companion(bdoy_2);

	asy = new array[0];
	for (i = 0; i < m_cY; ++i)
		asy ~= sprint("Y", "%c", 'a' + i);
	for (i = 0, clag = 1; i < m_cY; ++i)
		asy ~= sprint("Y", "%c", 'a' + i, "_1");
	if (bdoy_2)
	{	for (i = 0; i < m_cY; ++i)
			asy ~= sprint("Y", "%c", 'a' + i, "_2");
		++clag;
	}
	for (i = 0; i < m_cZ; ++i)
		asy ~= sprint("Z", "%c", 'a' + i);

    print("\nCompanion matrix:", "%10.3g", "%r", asy, md);
    if ( (m_iEdist == MVNORMAL || m_iEdist == NORMAL)
		&& (m_iZdist == MVNORMAL || m_iZdist == NORMAL || m_iZdist == NO_DIST) )
		print("\nVariance matrix of companion form:", "%10.3g", "%r", asy, ms);
	
	if (!(m_mB1t == 0) || !(m_mC2t == 0))
	{
   		print("(ignores ");
		if (!(m_mB1t == 0))
    		print("MA errors");
	    if (!(m_mC2t == 0))
    	    print("trend");
   	    print(")\n");
	}
	if (m_bCustomZ)
	    print("Ignoring custom Z's in asymptotic analysis\n");
	

	if (bstationary >= 0)
    {	print("\neigenvalues of companion matrix:",	"%10.3g", 
			"%c", {"real", "imag", "modulus"}, ev' ~ cabs(ev)');
	}
	else
    {   print("\nFailed to compute eigenvalues of companion matrix\n");
		return;
	}

	print("The DGP is", bstationary ? " " : " NOT ", "stationary.\n\n");

	if (bstationary && m_mC2t == 0)
	{
    	print("Long-run mean of y_t:", "%c", asy, mmy');
		print("M matrix", "%c", asy, mm);
		if (isarray(asModelY) && isarray(asModelX))
		{
			decl sely = strfind(asy, asModelY);
			decl selx = strfind(asy, asModelX);
			decl idxc = strfind(asModelX, "Constant");
			if (idxc >= 0) selx = dropc(selx, idxc);
			if (sely >= 0 && sizerc(selx) && selx >= 0)
			{
				print("Plims for OLS estimates of linear model",
					mmy == 0 ? "" : " with constant", ", T=", cT, ":");
				decl mxx = invertgen(mm[selx][selx]);
				decl mxy = mm[selx][sely];
				decl beta = mxx * mxy;
				decl ese = mm[sely][sely] - mxy'beta;
				decl bse = sqrt(diagonal(mxx)' .* diagonal(ese)) / sqrt(cT);
				if (sizerc(sely) == 1)
					println("%r", {"beta", "SE[beta]"}, "%c", array(asy[selx]), beta' | bse',
							"sigma^2 = ", double(ese));
				else
					println("\nbeta", "%c", asy[sely], "%r", array(asy[selx]), beta,
							"SE[beta]", "%c", asy[sely], "%r", array(asy[selx]), bse,
							"sigma^2", ese);
			}
		}
		print("\n");
	}

	decl bstationary_z = TRUE;
	if (m_cZ > 0 && !(m_mC0t == 0))
	{
		decl evz;
		if (eigen(m_mC0t, &evz) == 0)
			bstationary_z = !any(cabs(evz) .>= 1);
	}
	if (m_bUsedEcm)
	{
		print("Specified alpha:", m_mAlpha);
		print("Specified beta:", m_mBeta);
	}
	
	if (bstationary_z)	
	{	decl rankpi, vsv, maorth, mborth, ranki2, m, blag2;

    	print("Long-run matrix Pi(1)-I = Po:", mpi);
		rankpi = rank(mpi);
    	print("rank of long-run matrix: ", rankpi, "\n");

	    print("Mean-lag matrix sum pi_i:", mmeanlag);
		blag2 = !(m_mA5 == 0);
		if (blag2)
    		print("I(2) matrix Gamma:", mgamma);

		ranki2 = m_cY - rankpi;
		if (rankpi < m_cY)
		{
		    decsvd(mpi, &maorth, &vsv, &mborth);
			if (rankpi > 0)
    			print("Cointegrating space:", mborth[][:rankpi - 1]');
			m = maorth[][rankpi:]'mgamma*mborth[][rankpi:];
			if (rows(m) == 1 && columns(m) == 1)
				ranki2 = fabs(m[0][0]) < 1e-12 ? 0 : 1;
			else
				ranki2 = rank(m);
			if (blag2)
    			print("rank of I(2) condition: ", ranki2, "\n");
		}
		print("This DGP is ");
		if (blag2 && ranki2 < m_cY - rankpi)
   			print("I(2).\n");
		else if (rankpi < m_cY)
   			print("I(1).\n");
		else if (bstationary)
   			print("I(0).\n");
		else if (max(cabs(ev)) > 1)
   			print("explosive?\n");
		else
   			print("not stationary.\n");
	}
}
PcNaiveDgp::Asymp()
{
	AsympEx(0, 0, 0);
}
PcNaiveDgp::companion(const bdoY_2)
{
    decl cs, mlongmeany, mlongmeanz, mpi, mm, ms, mlongpi, ev, cy2 = 0, mp1,
		mp2, mp3, mh, istationary, mmeanlag, mgamma;

    /* using H = inv(I-A0), companion matrix D =                            */
    /*     H(A1+B0(I-A0)) :  H(A5-B0.A1) : -H.B0.A5 : H(A2.C0-B0.A2)        */
    /*          I                 0              0                0         */
    /*          0                 I              0                0         */
    /*          0                 0              0               C0         */
    /* variance matrix S =                                                  */
    /*     H(E1+A2.E2.A2')H' :  0  :  0  : H.A2.E2      */
    /*          0               0     0       0         */
    /*          0               0     0       0         */
    /*      E2.A2'.H'           0     0      E2         */
    /* long-run mean of y_t =                           */
    /*     inv(D-A1-A5).[A2.inv(I-C0)c1 + a3]           */
	mh = m_mDinv';

	mp1 = mh*(m_mA1 + m_mB0t' * (unit(m_cY) - m_mA0));
	mp2 = mh*(m_mA5 - m_mB0t'm_mA1);
	mp3 = 0;
	
	mpi = mp1 ~ mp2;		 // start building companion matrix
	
	if (bdoY_2)
	{	mp3 = -mh*(m_mB0t'm_mA5);
		mpi ~= mp3;
		cy2 = m_cY;
	}
	if (m_cZ)
    	mpi ~= mh*(m_mA2 * m_mC0t' - m_mB0t' * m_mA2);

	mpi |= (unit(m_cY) ~ zeros(m_cY, m_cY + cy2 + m_cZ));
	if (cy2 != 0)
    	mpi |= (zeros(m_cY, m_cY) ~ unit(m_cY) ~ zeros(m_cY, m_cY + m_cZ));
	
	if (m_cZ)
    	mpi |= (zeros(m_cZ, 2 * m_cY + cy2) ~ m_mC0t');

	decl mzsig = 0, mesig = 0;
    if (m_iEdist == MVNORMAL)	 mesig = m_mEdistS'm_mEdistS;
    else if (m_iEdist == NORMAL) mesig = diag(sqr(m_mEdistS));
    if (m_iZdist == MVNORMAL)	 mzsig = m_mZdistS'm_mZdistS;
    else if (m_iZdist == NORMAL) mzsig = diag(sqr(m_mZdistS));
		
	cs = rows(mpi);
	ms = zeros(cs, cs);	 // construct variance of companion form
	ms[][] = mh*(mesig + m_mA2 * mzsig * m_mA2')*mh';

	if (m_cZ)			 
    {	ms[][cs - m_cZ:] = mh * m_mA2 * mzsig;
		ms[cs - m_cZ:][] = (mh * m_mA2 * mzsig)';
		ms[cs - m_cZ:][cs - m_cZ:] = mzsig;
	}
	
	istationary = -1;
	if (eigen(mpi, &ev) == 0)
	{	istationary = !any(cabs(ev) .>= 1 - 1e-12);
		if (rows(ev) == 1)
			ev |= zeros(1, columns(ev));
	}

	mm = 0;
	mlongmeany = zeros(m_cY,1);
	mlongmeanz = zeros(m_cZ,1);
	
	if (istationary > 0 && m_mC2t == 0)
	{
		decl m1 = m_mA3, m2 = m_mC1t';
	    if (m_iEdist == MVNORMAL || m_iEdist == NORMAL)
			m1 += m_mEdistM';
	    if (m_iZdist == MVNORMAL || m_iZdist == NORMAL)
			m2 += m_mZdistM';

		if (m_cZ)
		{	mlongmeanz = invert(unit(m_cZ) - m_mC0t') * m2;
			m1 += m_mA2 * mlongmeanz;
		}
		mlongmeany = invert(unit(m_cY) - m_mA0 - m_mA1 - m_mA5) * m1;
		mm = shape(invert(unit(sqr(rows(mpi))) - mpi ** mpi) * vec(ms), rows(mpi), rows(mpi));
	}
	mlongpi = mp1 + mp2 + mp3 - unit(m_cY);
	mmeanlag = mp1 + 2*mp2 + 3*mp3;
	mgamma = unit(m_cY) + mp2 + 2*mp3;
	
	return {istationary, mpi, ms, ev, mm, mlongmeany, mlongmeanz, mlongpi, mmeanlag, mgamma};
}

