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