#!python

# bkfilter.py
# Ryan Banerjee
# 17th April 2009
# Purpose: Implement the Baxter King 1999 symmetric bandpass filter
# From matlab code: 
# http://www.clevelandfed.org/Research/Models/bandpass/bpassm.txt
# pl=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
# pu=dn is the higher frequency
#	for example, the business cycle band will have dn = 32 for quarterly data
#	and dn = 8 for annual data
# nfix 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 = 3 years is a practical choice.

from __future__ import division
import scipy as sp

def bkfilter(x,pl,pu,nfix):
    root=1
    drift=1
    ifilt=3
    theta=[1]

    nq=len(theta)-1
    b=2*sp.pi/pl
    a=2*sp.pi/pu

    T=len(x)
    assert (T>5)
    assert (pu>pl)
    assert (pl>1.1)

    cc = 1
    # Compute ideal Bs
    j = sp.linspace(1,2*T,2*T)
    B = sp.append([(b-a)/sp.pi], (sp.sin(j*b)-sp.sin(j*a))/(j*sp.pi))

    # Compute R using closed form integral solution
    R=sp.zeros(T)
    R0 = B[0]*cc
    R[0] = sp.pi*R0
    for j in range(1,T):
        dj = 2*sp.pi*B[j-1]*cc
        R[j] = R[j-1] - dj

    AA = sp.zeros((T,T))

    # Baxter King filter

    bb = sp.zeros(2*nfix+1)
    bb[nfix:2*nfix]=B[0:nfix]
    reverse = B[1:nfix]
    reverse = reverse[::-1]
    bb[1:nfix] = reverse 
    bb1 = bb
    bb1 = bb1 - sp.sum(bb1)/(2*nfix+1)
    for i in range(nfix,T-nfix):
        AA[i,i-nfix:i+nfix+1] = bb1

    # Compute filtered time series using selected filter matrix AA
    fX = sp.dot(AA,x)
    return fX

# test function
#import matplotlib.pyplot as plt
#a = sp.random.random(100)
#rho=0.95
#b=[0]
#for i in xrange(0,100):
#    b.append(b[i]*rho + a[i]-.5)
#
#c = bkfilter(b,1.5,8,3)
#
#plt.figure()
#plt.plot(b,'r')
#plt.plot(c,'b')
#plt.show()