#include <oxstd.h>

#import <lib/hpfilter>

bad_HPfilter(const vY, dLambda)
{
    decl i, mk, ma, ct = rows(vY);

    if (dLambda <= 1e-20)
        dLambda = 1600;

	for (i = 0, mk = new matrix[ct-2][ct]; i < ct - 2; i++)
        mk[i][i] = 1, mk[i][i + 1] = -2, mk[i][i + 2] = 1;
    ma = unit(ct) + dLambda*mk'mk; 

return (1/ma) * vY;	              /* inverts large (T x T) matrix */
}

main()
{
    decl my, ct = 100;

    my = rann(ct,1);	                     /* generate the data */
              /* spline is nearly the same; uses < 0 for bandwith */
    print("y,g =", my ~ HPfilter(my) ~ bad_HPfilter(my, 0)
		~ spline(my, 0, -1600));	
}

