/*--------------------------------------------------------------------------
 * pcfimldgp.ox - definitions/declarations for PcFimlDgp class
 *
 *       (C) Jurgen Doornik 1998
 *
 *--------------------------------------------------------------------------*/

#include <oxstd.h>
#include <oxfloat.h>
#include <arma.h>
#import <lib/longrun>
#include <pcfimldgp.h>

static decl sPcFimlDgpVersion = "2.00";

PcFimlDgp::PcFimlDgp(const cY, const cZ)
{
    this.Database();
	m_cY = cY;
	m_cZ = cZ;
	m_cYU = 2 * cY;
	m_cYUZ = 2 * cY + cZ;
	m_cYUZV = 2 * cY + 2 * cZ;

    m_mPit = 0;	 m_mC0t = 0.5 * unit(cZ);
    m_iEdist = NORMAL;
    m_mEdistM = zeros(1, cY);  m_mEdistS = ones(1, cY);
    m_iZdist = NORMAL;
    m_mZdistM = zeros(1, cZ);  m_mZdistS = ones(1, cZ);

	m_bFixedZ = TRUE;
	m_bHaveZ = FALSE;
	m_bUseObsLoop = FALSE;
	m_mxDgpLag = 4;
	
	SetInit(Y_DGP, 0);
	SetInit(Z_DGP, 0);
}
PcFimlDgp::Create(const iFreq, const iYear1, const iPeriod1,
	const cTdiscard, const mxDgpLag, const mxT)
{
	decl db = new Database();
	db.Create(iFreq, iYear1, iPeriod1, iYear1 + 2, iPeriod1);

	m_mxDgpLag = max(1, mxDgpLag);
	
	Database::Create(iFreq,
		db.ObsYear(- cTdiscard - m_mxDgpLag), db.ObsPeriod(- cTdiscard - m_mxDgpLag),
		db.ObsYear(mxT - 1), db.ObsPeriod(mxT - 1));
	delete db;

	decl i, j, as;

	for (i = j = 0, as = new array[m_cYUZV]; i < m_cY; ++i, ++j)
		as[j] = sprint("Y", string('a' + i));
	for (i = 0; i < m_cY; ++i, ++j)
		as[j] = sprint("U", string('a' + i));
	for (i = 0; i < m_cZ; ++i, ++j)
		as[j] = sprint("Z", string('a' + i));
	for (i = 0; i < m_cZ; ++i, ++j)
		as[j] = sprint("V", string('a' + i));

	Append(zeros(GetSize(), m_cYUZV), as);
	Deterministic(0);
}

PcFimlDgp::SetYParameter(const mPit)
{
	m_mPit = mPit;
}
PcFimlDgp::SetZParameter(const mC0)
{
	m_mC0t = mC0';
}

PcFimlDgp::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);
    }
	if (iDist > EXPONENTIAL)
		m_bUseObsLoop = TRUE;
}

PcFimlDgp::SetInit(const iDgp, const mInit)
{
	if (iDgp == Y_DGP)
	{
		m_mInitY = mInit;
		if (columns(m_mInitY) == 1)					   /* want a row vector */
			m_mInitY = m_mInitY';
		if (columns(m_mInitY) < m_cY)
			m_mInitY = zeros(1, m_cY);
	}
	else if (iDgp == Z_DGP)
	{
		m_mInitZ = mInit;
		if (columns(m_mInitZ) == 1)					   /* want a row vector */
			m_mInitZ = m_mInitZ';
		if (columns(m_mInitZ) < m_cZ)
			m_mInitZ = zeros(1, m_cZ);
	}
}

PcFimlDgp::SetFixedZ(const fSetting)
{
    m_bFixedZ = fSetting;
	m_bHaveZ = FALSE;
}
PcFimlDgp::DiscardZ()
{
	m_bHaveZ = FALSE;
}
PcFimlDgp::UseObsLoop(const bUseObsLoop)
{
	m_bUseObsLoop = bUseObsLoop;
}
PcFimlDgp::Prepare()
{
	decl msel, mselw, csel, i, j, cy;
	
	m_mxSelLag = GetMaxSelLag();

	if (m_mxSelLag > m_mxDgpLag)
	{
		println("Not enough lags allowed for in Create()");
		exit(1);
	}
								 /* set Y and Z to specified initial values */
	m_mData[][: m_cY - 1] = m_mInitY;
	if (m_cYUZ > m_cYU)
		m_mData[][m_cYU : m_cYUZ - 1] = m_mInitZ;

	m_iT1 = m_mxDgpLag;
	m_iT2 = Database::GetSize() - 1;
			
					    /* translate all groups not Y_VAR or X_VAR as X_VAR */
	m_mSelgroup = m_mSelgroup .> X_VAR .? X_VAR .: m_mSelgroup;
	
					  /* make selection in order in order Y, Y lagged, rest */
	csel = columns(m_mVarsel);
	msel = zeros(1, csel);
	for (i = j = 0; i < csel; ++i)
	{	if (m_mSelgroup[i] == Y_VAR && m_mLagsel[i] == 0)
			msel[j++] = i;
	}
	cy = j;
	mselw = zeros(1, csel - cy);
	for (i = 0; i < csel; ++i)
	{	if (m_mSelgroup[i] == Y_VAR && m_mLagsel[i] != 0)
			msel[j++] = i;
	}
	for (i = 0; i < csel; ++i)
	{	if (m_mSelgroup[i] != Y_VAR)
			msel[j++] = i;
	}
	decl msely = sortr(msel[][:cy-1]);
	for (i = cy - 1, mselw = msel[cy:]; i >= 0; --i)
		mselw = mselw .> msely[i] .? mselw - 1 .: mselw;

	m_mVarsel = m_mVarsel[msel];
	m_mLagsel = m_mLagsel[msel];
	m_mSelgroup = m_mSelgroup[msel];
	m_mPit = m_mPit[mselw][];
	
						  /* delete columns corresponding to dependent vars */
	m_mRegY = deleteifc(m_mVarsel, (m_mSelgroup .== Y_VAR) .&& (m_mLagsel .== 0));
	m_mLagY = deleteifc(m_mLagsel, (m_mSelgroup .== Y_VAR) .&& (m_mLagsel .== 0));
	m_mDepY = range(0, m_cY - 1);

	m_mLagZ = ones(1, m_cZ);
	m_mDepZ = range(m_cYU, m_cYUZ - 1);

	decl berror = FALSE;
	if (cy != m_cY)
	{	println("Error: all Ys must be a dependent variable");
		berror = TRUE;
	}
	if (columns(m_mRegY) != rows(m_mPit))
	{	println("Error: expected ", columns(m_mRegY)," rows in Pi\'");
		berror = TRUE;
	}
	if (columns(m_mPit) != m_cY)
	{	println("Error: expected ", m_cY," columns in Pi\'");
		berror = TRUE;
	}
	if (berror)
	{	print("Pi=", m_mPit);
		m_mPit = 0;
	}
	return !berror;
}

PcFimlDgp::GenerateV(const cT)
{
	return RanMC::RanDist(m_iZdist, cT, m_cZ, m_mZdistM, m_mZdistS);
}
PcFimlDgp::GenerateZ(const cT, const mC0t, const mV)
{
	return modelforc(mV, m_mData[ : m_mxSelLag + cT - 1][m_cYU : m_cYUZ - 1],
		m_mDepZ, m_mDepZ, m_mLagZ, m_mC0t', m_mxSelLag);
}
PcFimlDgp::GenerateU(const cT)
{
	return RanMC::RanDist(m_iEdist, cT, m_cY, m_mEdistM, m_mEdistS);
}
PcFimlDgp::GenerateY(const cT, const mPit, const mU)
{
	return modelforc(mU, m_mData[ : m_mxSelLag + cT - 1][],
			m_mDepY, m_mRegY, m_mLagY, m_mPit', m_mxSelLag);
}

PcFimlDgp::GenerateV_t(const iT)
{
	if (m_iZdist <= EXPONENTIAL)
	{
		if (iT == m_iT1)
			SetV(RanMC::RanDist(m_iZdist, m_iT2-m_iT1+1, m_cZ, m_mZdistM, m_mZdistS));
		return m_mData[iT][m_cYUZ : m_cYUZV - 1];
	}
	return RanMC::RanDist1(m_iZdist, m_cZ, m_mZdistM, m_mZdistS,
		m_mData[iT - 1][m_cYUZ : m_cYUZV - 1], m_mData[iT - 1][m_cYU : m_cYUZ - 1]);
}
PcFimlDgp::GenerateZ_t(const iT, const mC0t)
{
	return m_mData[iT - 1][m_cYU : m_cYUZ - 1] * mC0t;
}
PcFimlDgp::GenerateU_t(const iT)
{
	if (m_iEdist <= EXPONENTIAL)
	{
		if (iT == m_iT1)
			SetU(RanMC::RanDist(m_iEdist, m_iT2-m_iT1+1, m_cY, m_mEdistM, m_mEdistS));
		return m_mData[iT][m_cY : m_cYU - 1];
	}
	return RanMC::RanDist1(m_iEdist, m_cY, m_mEdistM, m_mEdistS,
		m_mData[iT - 1][m_cY : m_cYU - 1], m_mData[iT - 1][ : m_cY - 1]);
}
PcFimlDgp::GenerateY_t(const iT, const mPit)
{
	return selectrc(m_mData, iT - m_mLagY, m_mRegY) * mPit;
}

PcFimlDgp::GenerateTo(const cT)
{
	// m_mData[][0 : m_cY-1]          Y: dependent variables at lag 0
	// m_mData[][m_cY : m_cYU-1]      U: error terms for Y
	// m_mData[][m_cYU : m_cYUZ-1]    Z: regressors at lag 0
	// m_mData[][m_cYUZ : m_cYUZV-1]  V: error terms for Z

	m_iT1 = m_mxDgpLag;
	m_iT2 = m_mxDgpLag + cT - 1;

	if (m_bUseObsLoop)									 /* do with looping */
	{
		decl i, e;
		
		for (i = m_iT1; i <= m_iT2; ++i)
		{
			if (!m_bHaveZ && m_cZ)					/* step 1: generate V,Z */
			{	
				e = GenerateV_t(i);
				m_mData[i][m_cYU : ] = (GenerateZ_t(i, m_mC0t) + e) ~ e;
			}
			e = GenerateU_t(i);						/* step 2: generate U,Y */
			m_mData[i][] = (GenerateY_t(i, m_mPit) + e) ~ e;
		}
	}
	else											  /* do without looping */
	{
		if (!m_bHaveZ && m_cZ) 						/* step 1: generate V,Z */
		{	SetV(GenerateV(cT));
			SetZ(GenerateZ(cT, m_mC0t, GetV()));
		}
		SetU(GenerateU(cT));			   			/* step 2: generate U,Y */
		SetY(GenerateY(cT, m_mPit, GetU()));
	}

	m_bHaveZ = TRUE;

	return m_mData[m_iT1 : m_iT2][ : m_cYU - 1];
}

PcFimlDgp::GetY()
{
	return m_mData[m_iT1 : m_iT2][ : m_cY - 1];
}
PcFimlDgp::GetU()
{
	return m_mData[m_iT1 : m_iT2][m_cY : m_cYU - 1];
}
PcFimlDgp::GetZ()
{
	return m_mData[m_iT1 : m_iT2][m_cYU : m_cYUZ - 1];
}
PcFimlDgp::GetV()
{
	return m_mData[m_iT1 : m_iT2][m_cYUZ : m_cYUZV - 1];
}
PcFimlDgp::SetY(const m)
{
	m_mData[m_iT1 : ][ : m_cY - 1] = m;
}
PcFimlDgp::SetU(const m)
{
	m_mData[m_iT1 : ][m_cY : m_cYU - 1] = m;
}
PcFimlDgp::SetZ(const m)
{
	m_mData[m_iT1 : ][m_cYU : m_cYUZ - 1] = m;
}
PcFimlDgp::SetV(const m)
{
	m_mData[m_iT1 : ][m_cYUZ : m_cYUZV - 1] = m;
}

PcFimlDgp::Print()
{
	decl asy, asylag, asw, i;
	
    print("\n---- PcFiml (", sPcFimlDgpVersion, ") DGP ----\ny is (",
		m_cY, " x 1), z is (", m_cZ, " x 1) and",
		(m_bFixedZ) ? " " : " NOT ", "fixed.\n");
		
	print("\nDGP: y[t] = e[t] + Pi w\n");
	GetGroupLagNames(Y_VAR,	0, 0, &asy);
    for (i = 0; i < sizeof(asy); ++i)
        println("%-12s", asy[i], " Y_VAR: dependent variable");
	GetGroupLagNames(Y_VAR,	1, 1000, &asylag);
    for (i = 0; i < sizeof(asylag); ++i)
        println("%-12s", asylag[i], " Y_VAR: lagged dependent variable");
	GetGroupNames(X_VAR, &asw);
    for (i = 0; i < sizeof(asw); ++i)
        println("%-12s", asw[i], " X_VAR: regressor");
	println("Database sample: ", GetSample());

	print("\nCoefficients, Pi\'=", "%r", asylag ~ asw, "%c", m_asNames, m_mPit);
	if (!(m_mInitY == 0))
		print("Y initial values:", m_mInitY);
    RanMC::WriteDist("e", m_iEdist, m_mEdistM, m_mEdistS);
    print("\n");

    if (m_cZ == 0)
        return;

    print("z[t] = v[t]");
    if (!(m_mC0t == 0))
        print(" + C0 z[t-1]");
    print("\n");

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

PcFimlDgp::Asymp()
{
	::Companion(GetSelInfo(), Y_VAR, m_mPit', <>, TRUE, m_asNames, TRUE);
}

PcFimlDgp::AsympEx(const asModelY, const asModelX, const cT)
{
	Asymp();
}

