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