capture program drop bpass
program define bpass
* cfb
	version 6.0
/*			       STATA COMMAND FOR BAND PASS FILTER:  bpass X pl pu

		This is a Stata program that filters time series data using an approximation to the 
		band pass filter as discussed in the paper "The Band Pass Filter" by Lawrence J. 
		Christiano and Terry J. Fitzgerald (1999).

		Required Inputs:
		X     - series of data (T x 1)
		pl    - minimum period of oscillation of desired component 
		pu    - maximum period of oscillation of desired component (2<=pl<pu<infinity).

		Output:
		fX - matrix (T x 1) containing filtered data 

		Examples: 
 
		Quarterly data: pl=6, pu=32 returns component with periods between 1.5 and 8 yrs.
		Monthly data: pl=2, pu=24 returns component with all periods less than 2 yrs.

		Note:  When feasible, we recommend dropping 2 years of data from the beginning 
		       and end of the filtered data series.  These data are relatively poorly
		       estimated.

		===============================================================================
		This program contains only the default filter recommended in Christiano and 
		Fitzgerald (1999). This program was written by Eduard Pelz and any errors are 
		my own and not the authors. For those who wish to use optimal band-pass filters
		other than the default filter please use the Matlab version of code available 
		at www.clev.frb.org/economists/fitzgerald.
		===============================================================================

		Version Date: 1/24/00 (Please notify Eduard Pelz at eduard.pelz@clev.frb.org or 
		(216)579-2063 if you encounter bugs).  
*/
* cfb comment out
* set matsize 800

tempname X tl Xun a b bnot bhat	B Bprime bhatd bhatu AA
tempvar t tl B Btemp 

quietly {
	generate `X' = `1'
	local pl = `2'
	local pu = `3'
	local undrift = 1	

	if `pl' >= `pu' {di in red "Lower bound (pl) must be less than upper bound (pu)"
		exit}
	if `pl' < 2     {di in red "Lower bound (pl) must be greater than 2"
		exit}

	gen `t' = _n
	local T = _N
	mkmat `X', matrix(`X')

/* This section removes the drift from a time series using the 
   formula:  drift = (x(T) - x(1)) / (T-1).                     */

	if `undrift' == 1 {
		gen `tl' = `t'[_n-1]
		replace `tl' = 0 in 1
		mkmat `tl', matrix(`tl')

		local drift = (`X'[`T']-`X'[1])/(`T'-1)
	
		matrix `Xun' = J(`T',1,0)
		local i = 1
		while `i' <= `T' {
			matrix `Xun'[`i',1] = `X'[`i'] - (`tl'[`i']*`drift')
			local i = `i'+1    }
	}
	
	else matrix `Xun' = `X'

/*Create the ideal B's, then construct the AA matrix*/

	local T2=2*`T'
	local T1 = `T'-1		
	matrix `AA' = J(`T2',`T2',0)
	scalar `a' = 2*_pi/(`pu')
	scalar `b' = 2*_pi/(`pl')
	scalar `bnot' = (`b'-`a')/_pi
	scalar `bhat' = `bnot'/2	

	gen `B' = (sin(`t'*`b')-sin(`t'*`a'))/(_pi*`t')	
	gen `Btemp' = `B'[_n-1]
	replace `Btemp' = `bnot' in 1
	mkmat `Btemp', matrix(`B')
	matrix `Bprime' = `B''	

	local i = 1
	while `i' <= `T' {
		matrix substitute `AA'[`i',`i']=`Bprime'
		matrix substitute `AA'[`i',`i']=`B'
		local i = `i'+1 }
	
	matrix `bhatu' = J(`T',1,0)
	matrix `bhatd' = J(`T',1,0)
	matrix `bhatu'[1,1] = `bhat'
	matrix `bhatd'[`T',1] = `bhat'
	
	local i = 2
		while `i' <= `T' {
			local j = `i'-1
			local k = `T'-`j'
			matrix `bhatu'[`i',1]=(`bhatu'[`j',1]-`B'[`j',1])
			matrix `bhatd'[`k',1]=(`bhatu'[`j',1]-`B'[`j',1])
			local i = `i'+1 }
	
	matrix `AA' = `AA'[1..`T',1..`T']
	matrix substitute `AA'[1,1] = `bhatu'
	matrix substitute `AA'[1,`T'] = `bhatd'

/*Filter data using AA matrix*/
	
	matrix fX = `AA'*`Xun'}
end