/*--------------------------------------------------------------------------
 * hpfilter.ox - code to compute the Hodrick-Prescott filter
 *
 *       code (C) Jurgen Doornik 1995
 *
 * This code uses the structure of the projection matrix
 * to compute the HP filter. A bad implementation would
 * try to invert a T x T matrix. This is avoided here.
 *
 *
 * Syntax:
 *	HPfilter(const vY);
 *  HPfilter(const vY, const dLambda);
 *		vY[cT][1]     in:  T x 1 data matrix
 *      dLambda	      in:  bandwidth (default is 1600)
 *
 *	Returns a T x 1 matrix of filtered values.
 *--------------------------------------------------------------------------*/

#include <oxstd.h>

HPfilter(const vY, ...)
{
    decl i, vx, mh, ct = rows(vY), dlambda, args = va_arglist();

	if (sizeof(args) > 0)
	{	dlambda = args[0];
	    if (dlambda <= 1e-20)
    	    dlambda = 1600;
	}
	else
		dlambda = 1600;

	vx = vY - 2*lag0(vY, -1) + lag0(vY, -2);

    mh = new matrix[ct+2][1];
    mh[2:][] = solvetoeplitz( 6*dlambda+1 ~ -4*dlambda ~ dlambda,
		ct - 2, vx[ : ct - 3]);

	mh = dlambda * (mh - 2*lag0(mh, -1) + lag0(mh, -2));

return vY - mh[ : ct - 1];
}