screen on;

@ ***************************************************** @
@     PROGRAMME BY MORTEN O. RAVN, SEPTEMBER 2001       @
@ ***************************************************** @
@ THIS PROGRAM PRODUCES THE ALTERNATIVE VERSIONS OF THE
  HP-FILTER DESCRIBED IN "THE HP-FILTER IN CROSS-COUNTRY
  COMPARISONS" BY ALBERT MARCET AND MORTEN O RAVN, 2001 @
@ ***************************************************** @

@ ***************************************************** @
@ AS DESCRIBED IN THE PAPER, WE CONSIDER TWO ALTERNATIVE
  VERSIONS OF THE HP-FILTER REFORMULATED AS A CONSTRAINED
  MINIMIZATION PROBLEM

  RULE 1: MIN SUM(Y-YTR)^2 ST.
          SUM(YTR(t+1)-2YTR(T)+YTR(t-1))^2/SUM(Y-YTR)^2 <= V
  WHERE YTR IS THE ESTIMATE OF THE TREND COMPONENT

  RULE 2: MIN SUM(Y-YTR)^2 ST.
          1/(T-2)*SUM(YTR(t+1)-2YTR(t)+YTR(t+1))^2 <= W

  BELOW WE ALLOW FOR 2 DIFFERENT WAYS OF CHOOSING V AND W

  METHOD 1: IMPOSE V AND W EXOGENOUSLY - THIS SHOULD BE CHOSEN
            IF YOU HAVE SPECIFIC VALUES IN MIND. BELOW W
            IS SPECIFIED AS 100*(1/(T-2))*var(YTR(t+1)-2YTR(T)+YTR(t-1))

  METHOD 2: IMPOSE V AND W TO MATCH THE CORRSPONDING NUMBERS
            FOR A "NUMERIAIRE" COUNTRY. THIS IS OUR CHOICE AND
            WE CHOOSE THE US AS THE "NUMERAIRE" COUNTRY BUT
            OTHER CHOICES ARE POSSIBLE @
@ ***************************************************** @

@ ***************************************************** @
@ USER DEFINITION AREA - YOU SHOULD LOAD THE DATA AND
  MAKE YOUR CHOICES BELOW                               @
@ ***************************************************** @

@ 1. input the data -  if you have T observations and
  data for N countries this data should be loaded into the
  vector data(T,N). You can shorten the sample by redifining
  data just after you have loaded it if you wish @

load yus[159,1]=c:\data\albert\usyrea.txt;
load yrea[8,159]=c:\data\risk\gdprea.txt;
yrea=yrea';

data=yus~yrea[.,1]~yrea[.,5]~yrea[.,7];
data=data[1:157,.];

@ 2. specify the names of the countries for which you have data
  in the Nx1 vector dname (it should be dname="UK"|"GBR" etc. @

dname="us"|"aus"|"jap"|"uk";

@ 3. specify the starting date of the data (sdate) and the frequency
  of the observations (freqd, 1=annual, 4 = quarterly etc)    @

sdate=1960;
freqd=4;

@ 4. choose your method for selecting the values of V and W
  you do this by setting the parameter method=1 (for method
  1) or method=2 (for method 2). If you choose method 1 you will
  be asked to set V and W when you run the program. If you
  choose method 2 you will be asked to select the base
  country. Note that for method 2 the values of V and W depend
  on your choice of lstan specified below in step 5            @

method=2;
 

@ 5. choose the name of your output file
  this file will report standard deviations and authocorrelations
  of the data and will also compare the results to those for
  standard HP-methods. To that end you also need to specify
  the "standard" value of the smoothing parameter (lstan) - this
  will usually correspond to 1600 for quarterly data and
  100 for annual data (or 6.25 as argued by Ravn and Uhlig, 2001 @

output file = morten.out reset;
output off;
lstan=1600;

@ ********************************************************** @
@ END OF USER DEFINITION AREA - USUALLY YOU NEED NOT CHANGE
  ANYTHING BELOW                                             @
@ ********************************************************** @

DATS=SEQA(SDATE,1/FREQD,rows(DATA));
DATA=LN(DATA);
uu=0; #include c:\gauss1\albert\prog\hpnew;

nlags=5;
auci=zeros(nlags,cols(data));
auc1=zeros(nlags,cols(data));
auc2=zeros(nlags,cols(data));
sti=zeros(cols(data),1);
st1=zeros(cols(data),1);
st2=zeros(cols(data),1);
lamnew1=zeros(cols(data),1);
lamnew2=zeros(cols(data),1);

cls;
"The Marcet-Ravn Filter. Programme by M.O. Ravn ";

LOCATE 2,10;
IF METHOD==1;
"YOU HAVE CHOSEN TO SPECIFY MANUALLY THE VALUES OF V AND W";
"PLEASE TYPE YOUR PREFERRED VALUED OF V:";; V0=CON(1,1);
"PLEASE TYPE YOUR PREFERRED VALUED OF W:";; W0=CON(1,1);
ELSEIF METHOD==2;
"YOU HAVE CHOSEN TO SPECIFY A `NUMERAIRE' COUNTRY";
"PLEASE TYPE IN THE NUMBER OF THE COUNTRY FROM THE
FOLLOWING LIST";
I=1;
DO UNTIL I>COLS(DATA);
I;;$DNAME[I];
I=I+1;
ENDO;
NUM=CON(1,1);
ENDIF;
"DO YOU WANT TO SAVE THE FILTERED DATA? (Y=1, N=0)";
ANS=CON(1,1);
IF ANS==1;
"DATA AND DETRENDED DATA WILL BE SAVED IN HPRAVN.OUT";
DATARES = ZEROS(ROWS(DATA),COLS(DATA)*4); DATARES[.,1:COLS(DATA)]=DATA;
ENDIF;
"";
"PRESS ENTER TO START";
WAIT;
 
 

if method==2;
datnum = data[.,num];
endif;

if method==1;
lstan=1600/((4/freqd)^4);  @ this is based on ravn and uhlig, (2001) @

x2c=hpnew(data,lstan);
x2t = data-x2c;
dx2t = x2t[2:rows(x2t)]-x2t[1:rows(x2t)-1];
adx2 = dx2t[2:rows(dx2t)]-dx2t[1:rows(dx2t)-1];

W1 = 100*((stdc(adx2))^2);
V1 = (Stdc(adx2)^2)./(stdc(x2c)^2);
 

elseif method==2;
x1c=hpnew(datnum,lstan);
x1t = datnum-x1c;
dx1t = x1t[2:rows(x1t)]-x1t[1:rows(x1t)-1];
adx1 = dx1t[2:rows(dx1t)]-dx1t[1:rows(dx1t)-1];
x2c=hpnew(data,lstan);
x2t = data-x2c;
dx2t = x2t[2:rows(x2t),.]-x2t[1:rows(x2t)-1,.];
adx2 = dx2t[2:rows(dx2t),.]-dx2t[1:rows(dx2t)-1,.];
if ans>0; DATARES[.,COLS(dATA)+1:2*COLS(DATA)]=X2C; endif;

W1 = 100*((stdc(adx2))^2);
V1 = (Stdc(adx2)^2)./(stdc(x2c)^2);
W0 = 100*((stdc(adx1))^2);
V0 = ((Stdc(adx1))^2)/((stdc(x1c))^2);

endif;
dat1=uu;
sti = stdc(uu);
i=1; do until i>cols(data);
dat2=zeros(rows(uu)-nlags,nlags);
j=1;
do until j>nlags;
dat2[.,j]=trimr(dat1[.,i],nlags+1-j,j-1);
j=j+1;
endo;
corr=corrx(dat2);
auci[.,i]=corr[.,1];
i=i+1;
endo;

@ now iterate on the hp-filter so that K1 is matched for the
  series @

zh=1;
do until zh>cols(data);

data1=data[.,zh];

ZZ=1;
DO UNTIL ZZ>2;
cls;

if zz==1;     kold=W1[zh]; k1=w0;
elseif zz==2; kold=v1[zh]; k1=v0;
endif;

lamold = LSTAN;
lamorg=lstan;
lam=lstan;

dk=1;
i=1;

do until dk<.00000001;
uu=hpnew(data1,lam);
x2ta = data1-uu;
dx2ta = x2ta[2:rows(x2ta)]-x2ta[1:rows(x2t)-1];
adx2a = dx2ta[2:rows(dx2ta)]-dx2ta[1:rows(dx2ta)-1];
x2ca = uu;

if zz == 1;
kn = 100*(stdc(adx2a))^2;
elseif zz==2;
Kn = ((stdc(adx2a))^2)/((stdc(x2ca))^2);
endif;
 

locate 5,5;
"iteration number: ";;i;
locate 6,5;
"value of smoothing parameter: ";;lam;
locate 7,5;"previous value: ";;lamorg;
locate 8,5;
if zz==1;
"value of W: ";;kn;;" target: ";;w0;;
elseif zz==2;
"value of V: ";;kn;;" target: ";;v0;;
endif;

lamold = lam;

if zz == 1;
lam = lam*(1+.1*(kn-k1)/k1);
elseif zz==2;
lam = lam*(1+.1*(kn-k1)/k1);
endif;
 

kold = kn;
dk = abs((kn-k1)/k1);
i=i+1;
endo;

if zz==1;
lamnew1[zh]=lamold;
elseif zz==2;
lamnew2[zh]=lamold;

endif;

dat1 = uu;

@ compute moments @

if zz==1;
st1[zh,1] = stdc(dat1);
dat2=zeros(rows(dat1)-nlags,nlags);
j=1;
do until j>nlags;
dat2[.,j]=trimr(dat1[.,1],nlags+1-j,j-1);
j=j+1;
endo;
corr=corrx(dat2);
auc1[.,zh]=corr[.,1];
elseif zz==2;
st2[zh,1] = stdc(dat1);
dat2=zeros(rows(dat1)-nlags,nlags);
j=1;
do until j>nlags;
dat2[.,j]=trimr(dat1[.,1],nlags+1-j,j-1);
j=j+1;
endo;
corr=corrx(dat2);
auc2[.,zh]=corr[.,1];
endif;

if ans>0; DATARES[.,COLS(DATA)*(ZZ+1)+ZH]=UU; endif;

zz=zz+1;
endo;
 
 

zh=zh+1;
endo;
 
 
 

format /lds 8,8;
output on;
"results of the Marcet and Ravn (2001) filter";
"You chose the following method for seeting W and V:";
if method==1;"You have set V and W manually";
"Values: W:";;w0;;" and V:";;v0;
elseif method==2;
"You chose the following numeraire country:";;$dname[num];
"This produced the values: W:";;w0;;" and V:";;v0;
endif;
"your sample period: ";;dats[1];;" - ";;dats[rows(data)];
" ============================================================ ";
"A. Results for standard HP-method";
" Value of smoothing parameter: ";;
if method==1;
1600/((4/freqd)^4);  @ this is based on ravn and uhlig, (2001) @
elseif method==2; lstan;
endif;
"1. standard deviations of the cycle component * 100";
if method==1;
" country    std.";
i=1; do until i>cols(data);
$dname[i];; sti[i]*100;
i=i+1;
endo;
elseif method==2;
"country     std.      std/std. of num.country";
i=1; do until i>cols(data);
$dname[i];; sti[i]*100;; sti[i]/sti[num];
i=i+1;
endo;
endif;
"2. autocorrelations";
"country          order of autocorrelation";
"         ";;seqa(1,1,nlags)';
i=1; do until i>cols(data);
$dname[i];; auci[.,i]';
i=i+1;
endo;
" ============================================================ ";

"B. Results for Marcet and Ravn's method 1";
"1. implied value of lambda";
i=1; do until i>cols(data);
$dname[i];; lamnew2[i];
i=i+1;
endo;

"2. standard deviations of the cycle component * 100";
if method==1;
"country    std.";
i=1; do until i>cols(data);
$dname[i];; st2[i]*100;
i=i+1;
endo;
elseif method==2;
"country     std.     std/std. of num.country";
i=1; do until i>cols(data);
$dname[i];; st2[i]*100;; st2[i]/st2[num];
i=i+1;
endo;
endif;
"3. autocorrelations";
"country        order of autocorrelation";
"         ";;seqa(1,1,nlags)';
i=1; do until i>cols(data);
$dname[i];; auc2[.,i]';
i=i+1;
endo;

" ============================================================ ";

"B. Results for Marcet and Ravn's method 2";
"1. implied value of lambda";
i=1; do until i>cols(data);
$dname[i];; lamnew1[i];
i=i+1;
endo;

"2. standard deviations of the cycle component * 100";
if method==1;
" country    std.";
i=1; do until i>cols(data);
$dname[i];; st1[i]*100;
i=i+1;
endo;
elseif method==2;
"country     std.    std/std. of num.country";
i=1; do until i>cols(data);
$dname[i];; st1[i]*100;; st1[i]/st1[num];
i=i+1;
endo;
endif;
"3. autocorrelations";
"country              order of autocorrelation";
"         ";;seqa(0,1,nlags)';
i=1; do until i>cols(data);
$dname[i];; auc1[.,i]';
i=i+1;
endo;

output off;

screen off;

IF ANS>0;
OUTPUT FILE=HPRAVN.OUT RESET;
OUTPUT ON;
DATS=SEQA(SDATE,1/FREQD,ROWS(DATA));
ZZ=1;
DO UNTIL ZZ>4;
IF ZZ==1;
"THE DATA YOU PROVIDED:";
"DATE";;$DNAME';
X=DATS~DATA;
X;
ELSEIF ZZ==2;
"RESULTS OF HP FOR L=";;LSTAN;
"DATE";;$DNAME';
X=DATS~DATARES[.,COLS(DATA)+1:2*COLS(DATA)];
X;
ELSEIF ZZ==3;
"RESULTS FOR ADJUSTMENT RULE 1 WHERE LAMBDA WAS CHOSEN AS:";
$DNAME';
LAMNEW2';
"DATE";;$DNAME';
X=DATS~DATARES[.,COLS(DATA)*2+1:3*COLS(DATA)];
X;
ELSEIF ZZ==4;
"RESULTS FOR ADJUSTMENT RULE 2 WHERE LAMBDA WAS CHOSEN AS:";
$DNAME';
LAMNEW1';
"DATE";;$DNAME';
X=DATS~DATARES[.,COLS(DATA)*3+1:4*COLS(DATA)];
X;
ENDIF;
ZZ=ZZ+1;
ENDO;
OUTPUT OFF;
ENDIF;

PROC(1)=HPNEW(DAT,lambda);

local u1, u, l, ff, t, v11, v22, v12, x, z, v, tt, d, i, m1, m2, de,
      i1, ib, e1, e2, b11, b12, b22, jj, hu;
l=lambda;
u1 = dat;
ff=rows(dat');
t = rows(dat);
 

v11 = 1;
v22 = 1;
v12 = 0;
i=3; v=zeros(t,3); tt =zeros(t,ff); d = zeros(t,ff);
do until i>t;
 x=v11;
 z=v12;
 v11=(1/L)+4*(x-z)+v22;
 v12 = 2*x -z;
 v22=x;
 de = v11*v22-v12*v12;
 v[i,1] = v22/de;
 v[i,3] = v11/de;
 v[i,2] = -v12/de;
 x = v11+1;
 z=v11;
 v11 = v11-(v11*v11)/x;
 v22 = v22-(v12*v12)/x;
 v12 = v12-(z*v12)/x;
 i=i+1;
endo;

u=u1;
m1=u[2,.]; m2 = u[1,.];
i=3;
do until i>t;
 x = m1;
 m1=2*m1-m2;
 m2=x;
 tt[i-1,.] = v[i,1]*m1+v[i,2]*m2;
 d[i-1,.] = v[i,2]*m1+v[i,3]*m2;
 de = v[i,1]*v[i,3]-v[i,2]*v[i,2];
 v11 = v[i,3]/de;
 v12 = -v[i,2]/de;
 z = (u[i,.]-m1)/(v11+1);
 m1=m1+v11*z;
 m2=m2+v12*z;
 i=i+1;
endo;
tt[t,.] = m1;
tt[t-1,.] = m2;
m1 = u[t-1,.];
m2 = u[t,.];
i=t-2;
do until i<1;
 i1=i+1;
 ib=t-i+1;
 x=m1;
 m1=2*m1-m2;
 m2=x;
  if i>2;
   e1 = v[ib,3]*m2+v[ib,2]*m1+tt[i,.];
   e2 = v[ib,2]*m2+v[ib,1]*m1+d[i,.];
   b11 = v[ib,3]+v[i1,1];
   b12 = v[ib,2]+v[i1,2];
   b22 = v[ib,1]+v[i1,3];
   de = b11*b22-b12*b12;
   tt[i,.] = (-b12*e1+b11*e2)/de;
  endif;
 de = v[ib,1]*v[ib,3]-v[ib,2]*v[ib,2];
 v11 = v[ib,3]/de;
 v12 = -v[ib,2]/de;
 z = (u[i,.]-m1)/(v11+1);
 m1 = m1 +v11*z;
 m2 = m2 +v12*z;
 i=i-1;
endo;
tt[1,.] = m1;
tt[2,.] = m2;
i=1;
do until i>t;
 d[i,.] = u[i,.]-tt[i,.];
 i=i+1;
endo;
hu=tt;
uU = d-ones(t,ff)*MEANC(D);
retp(uu);
endp;