/*--------------------------------------------------------------------------
 * longrun.ox - code to obtain dynamic analysis of dynamic systems
 *
 *       code (C) Jurgen Doornik 1998-2003
 *
 * This code computes the long-run matrix/covariance, and the
 * static long run with standard errors as in PcFiml (see
 * the chapter on System Estimation in the PcFiml book).
 *
 *
 * Syntax:
 * 	LongRun(const amInfo, const iYvar, const miXvar,
 * 		const mPi, const mVarPi, const mOmega, const fOutput, const asNames)
 *		amInfo     in:  1 x 3 array with three matrices:
 *						database indices with selected vars (m_mVarSel)
 *                      group indicator of each selected var (m_mSelGroup)
 *                      lag length of each selected var (m_mLagsel)
 *      iYvar	   in:  int, group of dependent variable
 *      miXvar	   in:  1 x g matrix of integers, groups of regressor variables
 *		mPi		   in:  n x k matrix, (reduced form) coefficient matrix
 *				        NB: the columns of mPi must first have coeffs of
 *                      lagged Y's (in selection order), then the rest in
 *                      order of miXvar (within that the selection order, i.e.
 *						determined by amInfo)
 *		mVarPi	   in:  nk x nk variance-covariance matrix of (reduced form) coefs
 *		mOmega	   in:  n x n matrix of residual covariances
 *		fOutput	   in:  TRUE: print results
 *		asNames	   in:  database names
 *
 *	Returns a 1 x 7 array of matrices:
 *      Long-run matrix \Pi(1)-I = Po (n x n)
 *      Long-run covariance (n x n)
 *      Static long run	(n x r)
 *      Standard errors of static long run (n x r)
 *      Dependent variable selection matrix (1 x n)
 *      Regressor variable selection matrix (1 x r)
 *      Variance matrix of long-run coefficients (nr x nr)
 *
 *  For example, as called from within the PcFiml class:
 * 		LongRun(GetSelInfo(), Y_VAR, X_VAR ~ U_VAR, GetPi(),
 *			GetVarRf(), GetOmega(), m_fPrint, m_asNames);
 *
 *
 *	Companion(const amInfo, const iYvar, const mPi, const mVarPi,
 *		const fOutput, const asNames, const bIsDgp)
 *		amInfo     in:  1 x 3 array with three matrices:
 *						database indices with selected vars (m_mVarSel)
 *                      group indicator of each selected var (m_mSelGroup)
 *                      lag length of each selected var (m_mLagsel)
 *      iYvar	   in:  int, group of dependent variable
 *		mPi		   in:  n x k matrix, (reduced form) coefficient matrix
 *				        NB: the columns of mPi must first have coeffs of
 *                      lagged Y's (in selection order), then the rest in
 *                      order of miXvar (within that the selection order, i.e.
 *						determined by amInfo)
 *		mVarPi	   in:  NOT USED
 *		fOutput	   in:  TRUE: print results
 *		asNames	   in:  database names (not parameter names!)
 *		bIsDgp	   in:  TRUE: is DGP (assumes coefficients are exact)
 *
	return {istationary, mc, ev, mlongpi, mmeanlag, mgamma};
 *	Returns a 1 x 6 array of matrices:
 *      istationary: -1 (eigenvalues failed), 0 (stationary), > 0 (no of e.v.s
 *					 > 0 in absolute value)
 *      Companion matrix C (nm x nm), m = highest lag of dep var
 *      eigenvalues of C (2 x n) or (1 x n)
 *      Static long run	(n x n)
 *      Mean-lag matrix	(n x n)
 *      Gamma matrix	(n x n)
 *
 *  For example, as called from within the PcFiml class:
 *		Companion(GetSelInfo(), Y_VAR, GetPi(), <>,
 *			m_fPrint, GetAllNames(), FALSE);
 *--------------------------------------------------------------------------
 * 17-06-2003: LongRun: was printing negative long-run covariances as zero.
 *--------------------------------------------------------------------------*/
#include <oxstd.h>

LongRunInfo(const mVarsel, const mSelgroup, const mLagsel, const iYvar,
	const miXvar, const cP)
{
	decl i, ck, cy, cw, ctot, vsely, vselx, dbgroup, dbindex, dblag, sel;

	dbindex = mVarsel;						  /* variable index in database */
	dbgroup = mSelgroup;								  /* variable group */
	dblag   = mLagsel;										  /* lag length */
	
	sel = dbgroup .== vec(iYvar)|vec(miXvar);  /* drop all other group vars */
	dbindex = selectifc(dbindex, sel);
	dbgroup = selectifc(dbgroup, sel);
	dblag   = selectifc(dblag, sel);

	/* reassign new group numbers: -1:Y; 0: lagged Y; 1:miXvar[0], etc.  */
	decl lrgroup = constant(.NaN, sizer(dbgroup), sizec(dbgroup));
	decl cyg = sizerc(iYvar);
	for (i = 0; i < cyg; ++i)
	{	lrgroup = (dbgroup .== iYvar[i] .&& dblag .== 0) .? i - cyg .: lrgroup;
		lrgroup = (dbgroup .== iYvar[i] .&& dblag .!= 0) .? i       .: lrgroup;
	}
	for (i = 0; i < sizerc(miXvar); ++i)
		lrgroup = (dbgroup .== miXvar[i]) .? i + cyg .: lrgroup;

	/* sort selection according to groups; keeping order within group */
	ctot = columns(lrgroup);
	sel = sortbyc(lrgroup' ~ range(0, ctot-1)' ~ dbindex' ~ dblag', <0,1>);
	lrgroup = sel[][0]';
	dbindex = sel[][2]';
	dblag   = sel[][3]';

	/* build long-run selection matrices according to variable index */
	decl dbtolr, lrtodb, idx;/*use dbtolr to assign unique long-run var nos */
	dbtolr = constant(-1, 1, max(dbindex) + 1);/* lookup table for lr index */

	/* now the selection is in same order as the coeff estimates (we assume)*/	
	for (cy = 0; cy < ctot && lrgroup[cy] < 0; ++cy)	/* get no of y vars */
		dbtolr[dbindex[cy]] = cy;			/* first assign dep var indices */
	lrtodb = dbindex[0:cy - 1];			   /* dbindex of long-run variables */
	cw = idiv(cP, cy);						  		 /* no of rhs variables */

	for (i = ck = cy, sel = zeros(cy, cw); i < ctot; ++i)
	{
		idx = dbtolr[dbindex[i]];					 /* find long-run index */
		if (idx < 0)							   /* no index assigned yet */
		{	dbtolr[dbindex[i]] = idx = ck++;/* none: create new long-run var*/
			sel |= zeros(1, cw);
			lrtodb ~= dbindex[i];
		}
		sel[idx][i - cy] = 1;
	}
	ctot -= cy;							/* subtract rhs for dimension check */
	if (ctot < cw) /* yes: remaining coefs were not in the specified groups */			 
		sel |= zeros(cw - ctot, ctot) ~ unit(cw - ctot);   /* add remainder */

	if (sizerc(lrtodb) <= cy)					 /* no right-hand variables */
		return {lrtodb, <>, sel};
		
	vsely = lrtodb[ : cy - 1];
	vselx = lrtodb[cy : ];
	return { vsely, vselx, sel };
}

LongRun(const amInfo, const iYvar, const miXvar,
	const mPi, const mVarPi, const mOmega, const fOutput, const asNames, ...)
{
	decl cq, mp, mp0, mp0inv, mlrun, mj, mlrcov, mlrse, asy, ask, mlrvar;
	decl vsely, vselx, sel, cy;

	[vsely, vselx, sel] = LongRunInfo(amInfo[0], amInfo[1], amInfo[2],
		iYvar, miXvar, sizerc(mPi));
	cy = sizerc(vsely);
	if (cy == 0 || sizerc(vselx) == 0)
		return {<>, <>, <>, <>, vsely, vselx, <>};
	
	mp = mPi * sel[ : cy - 1][]';
	mp0 = sizerc(mp) ? mp - unit(cy) :  -unit(cy);
//	mp0inv = -(1 / mp0);
	mp0inv = -invert(mp0);
	if (mp0inv == 0)
	{
		println("\n*** Warning: long-run matrix is singular: long-run covariance unavailable");
		mlrcov = mlrun = <>;
	}
	else
	{	mlrcov = sizerc(mOmega) ? mp0inv * mOmega * mp0inv' : <>;
		mlrun = mp0inv * (mPi * sel[cy : ][]');
	}
	mlrse = mlrvar = <>;
	cq = columns(mlrun);
		
	if (sizerc(mVarPi) && cq)
	{
		mj = mp0inv ** ((mlrun' ~ unit(cq)) * sel);

		decl args = va_arglist();
		if (sizeof(args) && sizec(mj) == sizer(args[0]))
			mj *= args[0];
		if (! (sizeof(args) > 1 && args[1] == 1) )
			mlrvar = outer(mj, mVarPi);							  /* JVarJ' */

		mlrse = reshape(outer(mj, mVarPi, 'd'), cy, cq);/* diagonal(JVarJ') */
		mlrse = sqrt(mlrse);
	}
	
	if (fOutput)
	{
		if (sizerc(vsely))
		{
			asy = asNames[vsely];  if (isstring(asy))  asy = {asy};
			print("\nLong-run matrix Pi(1)-I = Po", "%r", asy, "%c", asy, mp0);
			if (sizerc(mlrcov))
			{
				decl m = fabs(mlrcov) .< 1e-12 .? 0 .: mlrcov;
				print("\nLong-run covariance", "%r", asy, "%c", asy, m);
			}
			if (sizerc(vselx) && sizerc(mlrun))
			{	ask = asNames[vselx];  if (isstring(ask))  ask = {ask};
				decl asva = va_arglist();
				if (sizeof(asva))
					ask ~= asva[0];
				print("\nStatic long run", "%r", asy, "%c", ask, mlrun);
				if (sizerc(mlrse))
			  		print("\nStandard errors of static long run",
						  "%r", asy, "%c", ask, mlrse);
			}
		}	
	}
	return {mp0, mlrcov, mlrun, mlrse, vsely, vselx, mlrvar};
}

EigenPrint(const mEv, const sDesc)
{
	decl mev = mEv';
	if (columns(mev) == 1) mev ~= 0;	// always print imag
	mev ~= cabs(mEv)';
	// sort in decreasing modulus
	mev = mev[sortcindex(-mev[][2])][];
	print("\nEigenvalues of ", sDesc, ":", "%#11.4g", 
		"%c", {"real", "imag", "modulus"}, mev);
}

CompanionPrint(const bIsDgp, const bDynamicAnalysis, const mxYlag,
	const asY, const iStationary, const mComp, const mEv, const mPi,
	const mMeanLag, const mGamma)
{
	decl cy = columns(mPi), i, j, k;

	if (rows(mComp))
	{
		decl asy = new array[mxYlag * cy];

		for (i = 0; i < cy; ++i)
			asy[i] = asY[i];
		for (j = 1, k = cy; j < mxYlag; ++j)
			for (i = 0; i < cy; ++i, ++k)
				asy[k] = sprint(asY[i], "_", "%d", j);
		
	    print("\nCompanion matrix:", "%9.3g", "%r", asy, mComp);
	}
	if (iStationary >= 0)
    {
//		print("\nEigenvalues of companion matrix:",	"%#11.4g", 
//			"%c", {"real", "imag", "modulus"}, mEv' ~ cabs(mEv)');
		EigenPrint(mEv, "companion matrix");
	}
	else
    {   print("\nFailed to compute eigenvalues of companion matrix\n");
		return;
	}

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

	
	if (bDynamicAnalysis)
	{
		decl rankpi, vsv, maorth, mborth, ranki2, m;

		if (rows(mPi))
			print("\nLong-run matrix Pi(1)-I:", mPi);
		if (bIsDgp)
		{	rankpi = rank(mPi);
    		print("rank of long-run matrix: ", rankpi, "\n");
		}

		print("\nMean-lag matrix sum pi_i:", "%r", asY, "%c", asY, mMeanLag);
		if (mxYlag > 1)
    		print("\nI(2) matrix Gamma:", "%r", asY, "%c", asY, mGamma);


		if (bIsDgp)
		{	ranki2 = cy - rankpi;
			if (rankpi < cy)
			{
			    decsvd(mPi, &maorth, &vsv, &mborth);
				if (rankpi > 0)
    				print("\nCointegrating 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 (mxYlag > 1)
    				print("rank of I(2) condition: ", ranki2, "\n");
			}
			print("\nThis DGP is ");
			if (mxYlag > 1 && ranki2 < cy - rankpi)
   				print("I(2).\n");
			else if (rankpi < cy)
   				print("I(1).\n");
			else if (iStationary)
   				print("I(0).\n");
			else if (max(cabs(mEv)) > 1)
   				print("explosive?\n");
			else
   				print("not stationary.\n");
		}
	}
}

Companion(const amInfo, const iYvar, const mPi, const mVarPi,
	const fOutput, const asNames, const bIsDgp)
{
	decl i, j, cy, cyl, msel, mreg, myreg, mlag, midxr, midxc, mxylag, mpi_i;

						/* NB: assumes columns of mPi first have lagged Y's */
				   /* select columns corresponding to lagged dependent vars */
			 /* msel has one row for each iYvar, but mreg/mlag have one row */
	msel = (amInfo[1] .== vec(iYvar)) .&& (amInfo[2] .> 0);
			    /* do the selection row by row: Y vars become before I vars */
	for (i = 0, mreg = mlag = <>; i < sizer(msel); ++i)
	{	mreg ~= selectifc(amInfo[0], msel[i][]);	  /* the lagged depvars */
		mlag ~= selectifc(amInfo[2], msel[i][]);		  /* and their lags */
	}
	cyl = columns(mreg);								  /* number of lags */
	midxc = range(0, cyl - 1);			           /* column indices in mPi */
	mxylag = int(max(mlag));							 /* the highest lag */
	
				/* select columns corresponding to (current) dependent vars */
	msel  = (amInfo[1] .== vec(iYvar)) .&& (amInfo[2] .== 0);
	for (i = 0, myreg = <>; i < sizer(msel); ++i)
		myreg ~= selectifc(amInfo[0], msel[i][]);			 /* the depvars */
	cy = columns(myreg);								   /* no of depvars */
	for (i = 0, midxr = mreg; i < cy; ++i)			  /* row indices in mPi */
		midxr = mreg .== myreg[i] .? i .: midxr;

	decl mc, msum, mlongpi, mmeanlag, mgamma, istationary, ev;

	mc = zeros(mxylag * cy, mxylag * cy);
	mlongpi = -unit(cy);
	mmeanlag = zeros(cy,cy);
	mgamma = unit(cy);

	for (i = 1, j = 0; i <= mxylag; ++i, j += cy) /* run over lagged Ys */
	{
		msel = selectifc(midxc, mlag .== i);
		mpi_i = zeros(cy, cy);
		if (sizerc(msel) != 0)					/* yes: depvars at this lag */
			mpi_i[][selectifc(midxr, mlag .== i)] = mPi[][msel];
		mc[0 : ][j : ] = mpi_i;
		if (i <	mxylag)
			mc[j + cy : ][j : ] = unit(cy);

		mlongpi += mpi_i;
		mmeanlag += i * mpi_i;
		mgamma += (i - 1) * mpi_i;
	}

	istationary = -1;
	ev = <>;

	if (sizerc(mc))
	{
		if (eigen(mc, &ev) == 0)
		{	istationary = !any(cabs(ev) .>= 1 - 1e-12);
			if (rows(ev) == 1)
				ev |= zeros(1, columns(ev));
		}	
		if (fOutput)
		{
			decl asy;
			asy = asNames[myreg];  if (isstring(asy)) asy = { asy };
			CompanionPrint(bIsDgp, TRUE, mxylag, asy, istationary, mc,
				ev, mlongpi, mmeanlag, mgamma);
		}
	}

	return {istationary, mc, ev, mlongpi, mmeanlag, mgamma};
}
