subroutine bkfilter(series y, series ybp, scalar n, scalar up, scalar dn, scalar k)
' Band-Pass Filter using the Baxter and King algorithm
'
' programmed by Aubhik Khan 06.03.2005
'
' y is a data vector of length n
' ybp is hte band pass filtered series, with NA for the first and last k observations
' n is the length of the data vector y
' n = @obsrange will yield number of observations in current workfile
' up is the lower frequency
' for example, the business cycle band will have up = 6 for quarterly data
' and up = 1.5 for annual data
' dn is the higher frequency
' for example, the business cycle band will have dn = 32 for quarterly data
' and dn = 8 for annual data
' k is the number of leads and lags used by the symmetric moving-average
' representation of the band pass filter. Baxter and King (1999) find
' that k = 12 is a practical choice.
scalar pi = @acos(-1.0)
scalar omlbar = 2.0*pi/dn
scalar omubar = 2.0*pi/up
rowvector(k+1) akvec
akvec(1) = (omubar - omlbar)/pi
for !j = 1 to k
akvec(!j+1) = sin(!j*omubar) - sin(!j*omlbar)
akvec(!j+1) = akvec(!j+1)/(!j*pi)
next
scalar theta = 0.0
for !j = 2 to k+1
theta = theta + 2.0*akvec(!j)
next
theta = akvec(1) + theta
theta = theta/(2.0*k + 1.0)
theta = -1.0*theta
akvec = akvec + theta
rowvector(2*k+1) avec
avec(k+1) = akvec(1)
for !j = 1 to k
avec(k+1-!j) = akvec(!j+1)
avec(k+1+!j) = akvec(!j+1)
next
for !j = k+1 to n-k
ybp(!j) = 0.0
for !l = 1 to 2*k+1
ybp(!j) = ybp(!j) + avec(!l)*y(!j - k - 1 + !l)
next
next
endsub