*
* @CFFilter( options )  x  start end  fx
* Christiano-Fitzgerald filter
*
* Parameters:
*   x            Input series
*   start end    Range of data to filter (by default, full range of x)
*   fx           Output series
*
* Options:
*   pl = minimum period of oscillation of desired component [6]
*   pu = maximum period of oscillation of desired component [32]
*      (2<=pl<pu<infinity).
*   [detrend]/nodetrend. With DETREND, removes a drift term before
*
* Examples:
*   Quarterly data: pl=6, pu=32 returns component with periods between 1.5 and 8 yrs.
*   Monthly data: pl=2, pu=24 returns component with all periods less than 2 yrs.
*
*  Note:  When feasible, we recommend dropping 2 years of data from the beginning
*         and end of the filtered data series.  These data are relatively poorly estimated.
*
* Reference: Christiano and Fitzgerald, "The Band Pass Filter",
*  International Economic Review, 2003, vol 44, no. 2, pp. 435-465
*=========================================================================
* This program contains only the default filter recommended in CF. This program was written
* by Eduard Pelz and modified by Fabio Canova
*=========================================================================
* Revision Schedule:
*   08/2003 Written by Eduard Pelz, eduard.pelz@clev.frb.org or (216)579-2063
*   07/2005 Syntax changed to allow sample range and to convert pl and pu to options, and
*     to add an option to control detrending. General code clean up. By Tom Doan, Estima.
*   05/2009 Correct problems with lead missing values in series
************************************************************************************
*
procedure cffilter x start end fx
type series x *fx
type integer start end
*
option real   pl      6
option real   pu      32
option switch detrend 1
*
local integer  startl endl m i j kk
local real     b a bnot bhat
local vect     bb2 fvec
local rect     aa  xvec
local real     mu

if .not.(%defined(x).and.%defined(fx))
   {
   display "Syntax: @cffilter( options ) x start end fx"
   return
}
*
if  pu <= pl
   {
   display "@cffilter: pu must be larger than pl"
   return
}
*
if  pl < 2
   {
   display "Warning: @cffilter, pl less than 2, reset to 2"
   return
   }

inquire(series=x) startl<<start endl<<end
compute m=endl-startl+1

*******************************************************************
*  This section is designed to remove the drift from a time series using
*  the formula:  drift = [x(endl) - x(startl)] / (m-1).
*******************************************************************

make xvec startl endl
# x
if detrend {
   compute mu=(x(endl)-x(startl))/(m-1)
   ewise xvec(i,j)=xvec(i,j)-(i-1)*mu
}

*******************************************************************
*Create the ideal B's then construct the AA matrix.
*******************************************************************
compute b    = 2*%pi/pl
compute a    = 2*%pi/pu
compute bnot = (b-a)/%pi
compute bhat = bnot/2

dimension bb2(m)
ewise bb2(i)=kk=(i-1),%if(kk==0,bnot,(sin(kk*b) - sin(kk*a))/(kk*%pi))
*
* Do the bulk of aa as a Toeplitz matrix from bb2
*
dimension aa(m,m)
ewise aa(i,j)=bb2(%iabs(i-j)+1)
*
* Reset the rim entries
*
compute aa(1,1) = bhat
compute aa(m,m) = bhat
do i = 1,m-1
   compute aa(i+1,1) = aa(i,1) - bb2(i)
   compute aa(m-i,m) = aa(i,1) - bb2(i)
end do
*
* Filter data using AA matrix
*
compute fvec=aa*xvec
*
* And copy over to the output series
*
set fx startl endl = fvec(t-startl+1)
*
end cffilter

