@=================================================================@
@:::::::::::::: HPFILT.SET :::::::::::::::::::::::::::::::::::::::@
@ Hodrick-Prescott Filter-----------------------------------------@
@-----------------------------------------------------------------@
@ USAGE: {yc,yt}=hpfilt(y,lam);                                   @
@ Inputs:   y: (T x 1) vector of observations to filter           @
@         lam: (scalar). HP's lambda                              @
@              1600 for quarterly obs.                            @
@              129600 for monthly obs.                            @
@ Output:  yc: The cyclical part                                  @
@          yt: The trend part                                     @
@-----------------------------------------------------------------@
proc(2) = hpfilt(y,lam);
local nc,t0,t1,G,e0,e,j,k0,k1,sto,tau,yfilt,sto1;

nc=cols(y); t0=1; t1=rows(y);
@===============================================================@
@:::::::::::::::: Begin Filtering ::::::::::::::::::::::::::::::@
@---------------------------------------------------------------@

G=zeros(t1,t1);
let e0=1 -2 1; G[1,1:3]=e0';
let e0=-2 5 -4 1; G[2,1:4]=e0';
let e0=1 -4 5 -2; G[t1-1,t1-3:t1]=e0';
let e0=1 -2 1; G[t1,t1-2:t1]=e0';

let e=1 -4 6 -4 1;
j=3; do while j le t1-2;
 k0=j-2; k1=j+2;
 G[j,k0:k1]=e';
j=j+1; endo;                    @ The TxT G matrix @

sto=lam*G+eye(t1); 

j=1; do while j le nc;
 if j eq 1; tau=inv(sto)*y[.,j]; endif; 
 if j gt 1; sto1=inv(sto)*y[.,j]; tau=tau~sto1; endif;
j=j+1; endo;

yfilt=y-tau;         @ The cyclical (filtered) data @

retp(yfilt,tau);
endp;