# HP-Filter function
# by Christian Zimmermann, 3 Jan 2005
# Based on Fortran subroutine by Edward Prescott
# Takes as input the penalty value and a series (both in a single array)
# Output is the filtered trend.  
sub hpfilter{
    $S=$_[0];
    @Y=@_;
    $N=$#Y;
#
#     compute sequences of covariance matrix for f[x(t),x(t-1) | y(<t)]
#
    $SS=$S;
    $NN=$N;
    $V11=1;
    $V22=1;
    $V12=0;
for($I=3;$I<=$N;$I++){
    $X=$V11;
    $Z=$V12;
    $V11=1/$S + 4*($X-$Z) + $V22;
    $V12=2*$X - $Z;
    $V22=$X;
    $DET=$V11*$V22-$V12*$V12;
    $V[$I][1]=$V22/$DET;
    $V[$I][3]=$V11/$DET;
    $V[$I][2]=-$V12/$DET;
    $X=$V11+1;
    $Z=$V11;
    $V11=$V11-$V11*$V11/$X;
    $V22=$V22-$V12*$V12/$X;
    $V12=$V12-$Z*$V12/$X;
}
#
#     this is the forward pass
#
$M1=$Y[2];
$M2=$Y[1];
for($I=3;$I<$N+1;$I++){
    $X=$M1;
    $M1=2*$M1-$M2;
    $M2=$X;
    $T[$I-1]=$V[$I][1]*$M1+$V[$I][2]*$M2;
    $D[$I-1]=$V[$I][2]*$M1+$V[$I][3]*$M2;
    $DET=$V[$I][1]*$V[$I][3]-$V[$I][2]*$V[$I][2];
    $V11=$V[$I][3]/$DET;
    $V12=-$V[$I][2]/$DET;
    $Z=($Y[$I]-$M1)/($V11+1);
    $M1=$M1+$V11*$Z;
    $M2=$M2+$V12*$Z;
}
$T[$N]=$M1;
$T[$N-1]=$M2;
#
#       this is the backward pass
#
$M1=$Y[$N-1];
$M2=$Y[$N];
#        DO 15 I=N-2,1,-1
for($I=$N-2;$I>0;$I--){
    $I1=$I+1;
    $IB=$N-$I+1;
    $X=$M1;
    $M1=2*$M1 - $M2;
    $M2=$X;
#
#           combine info for y(.lt.i) with info for y(.ge.i)
#
    if($I>2){
	$E1=$V[$IB][3]*$M2 + $V[$IB][2]*$M1 + $T[$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];
	$DET=$B11*$B22-$B12*$B12;
	$T[$I]=(-$B12*$E1+$B11*$E2)/$DET;
    }
#
#           end of combining
#
    $DET=$V[$IB][1]*$V[$IB][3]-$V[$IB][2]*$V[$IB][2];
    $V11=$V[$IB][3]/$DET;
    $V12=-$V[$IB][2]/$DET;
    $Z=($Y[$I]-$M1)/($V11+1);
    $M1=$M1+$V11*$Z;
    $M2=$M2+$V12*$Z;
}
$T[1]=$M1;
$T[2]=$M2;
shift(@T);
return @T;
}
1;
