/*
Hodrick-Prescott Filter (Matheny,vanNorden,Vigfusson)
(Paul Soderlin's HPMat2 program follows this,
and Lubuele's program follows that.) */

/*  HPFILTER    First Draft : December 21, 1989 by Ken Matheny (hodpres())
                Second Draft: August 3, 1994 by Simon van Norden (hpfilter())
                Third Draft : April 27, 1995 by Robert Vigfusson (hpfilter())

    This program is written and submitted for public, non-commercial use.
    The authors make no guarantees as to its performance, and request that
    use of this code be acknowledged in any published research, including
    working papers.

    This program is designed to smooth the stochastic series  y  by use
    of the Hodrick-Prescott method mentioned in footnotes of the Kydland
    & Prescott "Time to Build" paper, and Hansen's "Indivisible Labor"
    paper.

    The method involves minimizing an expression containing the squared
    differences between the actual and smoothed series, and terms involving
    the differences between successive elements of the smoothed series.
    This leads to a linear filter for the original data.

    "w" is a smoothness parameter.  The larger w, the smoother the resulting
    series, shp.  Kydland and Prescott suggest a "w" of 1600 for quarterly data.

FORMAT : {s} = hpfilter(y,w);

INPUTS : y  : The series to be filtered;
         w  : The smoothness parameter (suggest 1600).

OUTPUT : shp: The smoothed series.  The deviations can be calculated simply
              as (y-shp).

*/

proc(1) = hpfilter(y,w);

local t,a,b,c,d,m,i,j,shp;

    t       = rows(y);
    a       = 6*w+1;
    b       = -4*w;
    c       = w;
    d       = c~b~a;

d       = ones(t,1)*d;

d[1,1]=0;
d[2,1]=0;
d[1,2]=0;

d[2,2] = -2*w;
d[t,2] = -2*w;

d[1,3] = 1+w;
d[t,3] = 1+w;

d[2,3]=5*w+1;
d[t-1,3]=5*w+1;

@ The Smoothed Series @

shp = bandsolpd(eye(t),d)*y;

retp(shp);
endp;

/********************************
Gaston Luis Repetto  grepet1@icarus.cc.uic.edu offered the following comments:
For monthly data, Edward Prescott suggested (by email) to begin the analysis with
14,400  and then try to see with which value the adjusted
series better mimic the behavior of the original one.

"There is no such a thing as a "right" smoothing parameter. For purposes of
postwar business cycles, the 1600 number was the first tried and seem to
mimic well the smooth curve people typically fit throughout the data.  As a
year is four quarters, reducing it by four squared, that is by  factor of
16, resulted in the same amount of variation in trend being permitted.  This
logic suggest multiplyings 1600 by 9  (three squared) to have same amount of
variability in trend with monthly data as 1600 permits for quarterly data.

Some good people tried to develop some formal theory, that made connections
between different length time periods.  The "right" method depended upon
many very specific probabilistic assumptions.  But this is a data
representation exercise, not a measurement or estimation procedure."
***********************************/







/*-----------------------------------------------------------------------
**  HPMat2
**  Author: Paul Soderlind <nepso@HHS.SE>
**  An efficient HP-filter. Any comments are welcome.
**  Provided without guarantees for free public, non-commercial use.
**  Date: 30 May 1997
**
**
**  Purpose:   HP-filtering of a several time series, assuming no missing
**             values.
**
**  Usage:     { gy,cy } = HPMat2( y,lambda );
**
**  Input:     y            NxM matrix with M time series
**             lambda       degree of smoothness (1600 is a common value)
**
**  Output:    gy           NxM matrix of trend component
**             cy           NxM matrix of cyclical component
**
**
**
**  Note:      (a) This procedure uses Toeplitz() to set up the K matrix,
**             and the band matrix routines for inverting. This gives
**             a fast code with very small storage requirements.
**
**             (b) The trend component, gy, is calculated as
**             cy = inv(A) * y, where  A =  I + lambda*K'K. The cyclical
**             component is then cy = y - gy. The matrix K is as follows
**
**                    1 -2  1  0  0 ... 0  0  0
**                    0  1 -2  1  0 ... 0  0  0
**               K =  0  0  1 -2  1 ... 0  0  0
**                    .  .  .  .  .     .  .  .
**                    .  .  .  .  .     .  .  .
**                    0  0  0  0  0 ... 1 -2  1
**
**
**  Reference: Danthine and Girardin (1989), "Business Cycles In Switzerland:
**             A Comparative Study," European Economic Review 33, 31-50.
**
**
**  Paul Soderlind, 30 May 1997.
-----------------------------------------------------------------------*/

proc (2) = HPMat2( y,lambda );

  local N,gk,K,A,bandA,CholA,gy,cy;


  N = rows(y);

  if N < 5;                         /*pattern of A different if N < 5*/
    retp( error(0),error(0) );
  endif;

  gk = 1|(-2)|1|zeros(2,1);       /*prototype 5x5 K matrix*/
  K = upmat(toeplitz(gk));
  K = trimr(k,0,2);

  A = eye(5)+lambda*k'*k;         /*prototype 5x5 A matrix*/
  bandA = band(A,2);            /*band from A, extending to TxT matrix*/
  bandA = bandA[1:2,.] | (ones(N-4,1).*bandA[3,.]) | bandA[4:5,.];

  CholA = bandchol( bandA );    /*Cholesky decomp of band representation*/

  gy = bandcholsol( y,CholA );      /*solving for "trend" component*/
  cy = y - gy;                      /*"cyclical" component*/

  retp( gy,cy );

endp;
/*-----------------------------------------------------------------------*/









/*
ANOTHER HPFILTER:
Author: Luan Lubuele
		Senior Economist
		City of New York
		Independent Budget Office
		lubuele@nwu.edu
Date:   11 May 1997
Provided without guarantees for free public, non-commercial use.


The HP filter code posted by M. Heiko uses the inverse method. It should do
the job and is the shortest and fastest method. It would be simpler to write
the code as a procedure yielding directly the trend so as to avoid holding a
huge inverse in memory. There is however the issue of interdmediate swelling
and loss of accuracy. Indeed the matrix to be inverted is of the size n^2 (n
being the sample size), which can get big for say weekly data for a long
period of time. Plus, there will be a point where Gauss will need to hold
temporarely another matrix of the same size. As a result, one may run out of
memory. On the other end, the use of the straight inverse function on a huge
matrix may yield less accurate results. I have seen occasions when the
function inv(.) failed to pick up a clear case of a matrix with incomplete
rank. I am enclosing the HP-filter procedure which (1) does not need an
inverse, (2) uses just storing space of the order of 5*n instead of n^2. Its
shortfalls are (1) it is slow since it relies on loops, and (2) it does not
allow the computation of the hptrends of many series at once. I am also
enclosing the MATLAB version of the method based on the inverse, just in
case there is now a way of handling sparse matrices and someone is
interested in translating in Gauss.
The problems discussed above are addressed by using sparse matrix methods
and LU decomposition.


/****************** Lubuele GAUSS CODE STARTS HERE **********************/
declare __lambda != 1600;

proc (2) = hpfilter(y);
/* Computes the Hodrck-Prescott filtered trend and cycles.
** Based on Prescott algorithm using the two pass kalman filter approach
** This code is a translation of Prescott,s fortran code.
** The procedure needs the global variable __lambda (default=1600).
*/
local s,n,ierr,t,d,v,m1,m2,v11,v12,v22,x,z,b11,b12,b22,dett,e1,e2;
local istep,i1,i2,i;
s=__lambda;
n=rows(y);
t=zeros(n,1);
d=t;
v=zeros(n,3);
if s<= 0;
  s=1600;
endif;
if n<=3;
  retp(0,0);
endif;
istep=1;
m1=y[2];
m2=y[1];
i1=3;
i2=n;

v11=1;
v22=1;
v12=0;
i=i1;
do while i<=i2; /* first pass */
  gosub pass;
endo;
t[n]=m1;
t[n-1]=m2;
istep=-1;
m1=y[n-1];
m2=y[n];
i1=n-2;
i2=1;
v11=1;
v22=1;
v12=0;
i=i1;
do while i>=i2; /* second backward pass */
  gosub pass;
endo;
t[1]=m1;
t[2]=m2;
i=1;
do while i<=n;
  d[i]=y[i]-t[i];
  i=i+1;
endo;
retp(t,d);

pass:
  x=m1;
  m1=2*m1-m2;
  m2=x;
  x=v11;
  z=v12;
  v11=1/s+4*(x-z)+v22;
  v12=2*x-z;
  v22=x;
  dett=v11*v22-v12*v12;
  if istep==1;
    v[i-1,1]=v22/dett;
    v[i-1,3]=v11/dett;
    v[i-1,2]=-v12/dett;
    t[i-1]=v[i-1,1]*m1+v[i-1,2]*m2;
    d[i-1]=v[i-1,2]*m1+v[i-1,3]*m2;
  else;
    if i>=2;
      b11=v11/dett;
      b12=-v12/dett;
      b22=v22/dett;
      e1=b11*m2+b12*m1+t[i];
      e2=b12*m2+b22*m1+d[i];
      b12=b12+v[i,2];
      b22=b22+v[i,3];
      b11=b11+v[i,1];
      dett=b11*b22-b12*b12;
      t[i]=(-b12*e1+b11*e2)/dett;
    endif;
  endif;
  x=v11+1;
  z=(y[i]-m1)/x;
  m1=m1+v11*z;
  m2=m2+v12*z;
  z=v11;
  v11=v11-v11*v11/x;
  v22=v22-v12*v12/x;
  v12=v12-z*v12/x;
  i=i+istep;
return;
endp;

/*********************** Lubuele GAUSS CODE ENDS HERE ***************************/

