/* HP - Hodrick - Prescott Filter by M. Watson, 1991 Usage: xp = HP(x,l)
where: x = NxK matrix of K series of length N to be filtered; l =
smoothness parameter (1600 for quarterly data); xp = NxK matrix of K series
of length N representing the percentage deviations of each of the original
series in x from their trend. Note: if you want the trend, change the next
to last line to: retp(xt); if you want the smoothed series, use retp(xd).
*/ 
proc HP(x,l); 
local j,xd,xt,i,k,ty,cy,n,xhp,q,a; 
n=rows(x); 
j=1;
a=zeros(n-2,n); 
do until j>n-2; 
a[j,j]=1;
 if j+1 le n;
a[j,j+1]=-2;
endif;
if j+2 le n; 
a[j,j+2]=1; 
endif; 
j=j+1; 
endo; 
i=eye(n);
xd=zeros(rows(x),cols(x)); 
xt=xd; 
xhp=xt; 
k=1; 
do until k>cols(x);
ty=solpd(x[.,k],i+l*a'a); 
cy=x[.,k]-ty; 
xd[.,k]=cy; 
xt[.,k]=ty; 
k=k+1;
endo; 
xhp=xd./xt*100; 
retp(xhp); 
endp;









function [s]=hpfilter(y,w)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Author: Ivailo Izvorski, 
%          Department of Economics
%	   Yale University.
%          izvorski@econ.yale.edu
%  This code has been used and seems to be free of error.
%  However, it carries no explicit or implicit guarantee.
%
%  function [s]=hpfilter(y,w)
%  Hondrick Prescott filter where:
%  w - smoothing parameter; w=1600 for quarterly data
%  y - the original series that has to be smoothed
%  s - the filtered series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(y,1)<size(y,2)
   y=y';
end
t=size(y,1);
a=6*w+1;
b=-4*w;
c=w;
d=[c,b,a];
d=ones(t,1)*d;
m=diag(d(:,3))+diag(d(1:t-1,2),1)+diag(d(1:t-1,2),-1);
m=m+diag(d(1:t-2,1),2)+diag(d(1:t-2,1),-2);
%
m(1,1)=1+w;       m(1,2)=-2*w;
m(2,1)=-2*w;      m(2,2)=5*w+1;
m(t-1,t-1)=5*w+1; m(t-1,t)=-2*w;
m(t,t-1)=-2*w;    m(t,t)=1+w;
%
s=inv(m)*y;










@ This program reads and filters the citibase series @
@ Make sure that you change the directory's address @
library pgraph;
#include c:\Mydoc\pigia\courses\econ984\hp.txt;  @ execute procedure for hp filter @

nobs=136;          @ number of observations @
nser=4;
l=1600;            @ smoothing parameter @
load matser[]=c:\mydoc\pigia\courses\econ984\series;
matser2=reshape(matser,nobs,nser);
matserf=hp(matser2,l);

corr_mat = corrx(matserf);
vcov=vcx(matserf);

/*print "Correlation matrix  ";
print " Cons   Output   Invest   Hours";
print corr_mat;

stdev=vcov^(1/2);
print "standard deviations ";
print " Cons   Output   Invest   Hours";
print stdev;*/

c_raw = matser2[.,1];
c_fil = matserf[.,1];
periods = seqa(1,1,136);
title("c raw");
xy(periods,c_raw);

title("c fil");
xy(periods,c_fil);

stdevcy=sqrt(vcov[1,1])/sqrt(vcov[2,2]);
stdeviy=sqrt(vcov[3,3])/sqrt(vcov[2,2]);
stdevny=sqrt(vcov[4,4])/sqrt(vcov[2,2]);

print "sd(c)/sd(y) sd(i)/sd(y) sd(n)/sd(y)";
print stdevcy;
print stdeviy;
print stdevny;

print "corr c,i,n & y";
print corr_mat[1,2];
print corr_mat[2,2];
print corr_mat[3,2];
print corr_mat[4,2];