subroutine bkfilter(n, y, up, dn, k, ybp)

! Aubhik 20.02.2005
!
! Baxter and King (1999) Band-Pass filter
!
! y is a data vector of length n
! ybp is the band pass filtered series, with 0.0 for the first and last k observations
! n is the length of the data vector y
! 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.
! 
! Algorithm provided by Bob King.  This subroutine is the Fortran 90 version of 
! bpf.m and filtk.m written by Baxter and King (1999).  The band pass filter is 
! a moving average involving 2k+1 terms.  The weights are symmetric and calculated
! in akvec.  

implicit none

integer, parameter:: rk = selected_real_kind(15,307), ik = selected_int_kind(9)

integer(ik):: n, up, dn, k
real(rk):: y(n), ybp(n)

integer(ik):: j
real(rk):: pi, omlbar, omubar, kr, kj, theta
real(rk):: akvec(k+1), avec(2*k+1)

intent(in):: n, up, dn, k, y
intent(out):: ybp

pi = 2.0_rk*acos(0.0_rk)

omlbar = 2.0_rk*pi/dble(dn)
omubar = 2.0_rk*pi/dble(up)

akvec = 0.0_rk
avec = 0.0_rk

kr = dble(k)

akvec(1) = (omubar - omlbar)/pi

do j = 1, k
    kj = dble(j)
    akvec(j+1) = (dsin(kj*omubar) - dsin(kj*omlbar))/(kj*pi)
end do

theta = akvec(1) + 2.0_rk*sum(akvec(2:k+1))
theta = -1.0_rk*(theta/(2.0_rk*kr + 1.0_rk))

akvec = akvec + theta

avec(k+1) = akvec(1)

do j = 1, k
    avec(k+1 - j) = akvec(j+1)
    avec(k+1 + j) = akvec(j+1)
end do

ybp = 0.0_rk

do j = k + 1, n - k
    ybp(j) = dot_product(avec, y(j - k: j + k))
end do

end subroutine bkfilter