*
* @BKfilter( options ) series start end outser
* Baxter-King band-pass filter
*
* Parameters:
*
* SERIES = series to be filtered
*
* START, END = starting and ending dates/entry numbers specifying
* the range over which you want to run the filter.
* Defaults to the defined range of the SERIES.
*
* OUTSER = output series
*
* Options:
*
* UPPER = [2] period corresponding to upper cutoff frequency
* (nonnegative integer; UPPER=2 is a highpass filter)
*
* LOWER = [32] period corresponding to lower cutoff frequency
* (nonnegative integer less than UPPER; note that
* if LOWER > number of observations, frequency is set
* to zero)
*
* NMA = [12] controls size of symmetric moving average filter
* (nonnegative integer; number of terms is 2*NMA + 1)
*
* [DIFF]/NODIFF padding uses AR on 1st differences (NODIFF=levels)
*
* AR = [4] AR order for padding
*
* (padding adds artificial data at start
* and end of series to permit moving-average method
* of filtering using AR backcasts and forecasts)
*
*
* @BKFILTER executes a bandpass filter on a series using the method of
* Baxter and King (1999). Users can determine the selection of
* frequencies in the filter using the UPPER and LOWER cutoff variables.
* The default option passes frequencies corresponding to between 2 and
* 32 periods, a typical business-cycle frequency range with quarterly
* data.
*
* The filter uses a centered moving average method using up to NMA
* weighted leads and lags. Thus, for the filter to work it is necessary
* to pad the series at the start and end with NMA observations using AR
* backcasts and forecasts. AR controls the order of this AR.
*
* [Corollary: Missing data are a big problem! Interpolate gaps before
* using BKFILTER.]
*
* The output series is the filtered data (i.e., the cyclical component)
* and is saved in the user-specified output series.
*
* Default options are based on the parameters used for U.S. quarterly
* data by Stock and Watson (1999) and on Baxter and King's (1999)
* settings. The file BKFILTER.PRG replicates the Stock and Watson filter
* on postwar US quarterly GDP, and demonstrates the syntax.
*
* Defaults: UPPER=2, LOWER=32, NMA=12, ARPAD=4
* Suggested setting for annual data: UPPER=2, LOWER= 8, NMA= 3, ARPAD=1
*
*
* Revision History:
* 12/1998 Translation to RATS by Alan Taylor (now at UC-Davis)
* Based on James Stock and Mark Watson's GAUSS code BPASS.PRC,
* PAD.PRC. Based, in turn, on the MATLAB code of Marianne Baxter
* and Robert King. My thanks to Mark Watson for sharing his code.
* 09/1999 Modifications made to BKPAD to properly process
* explicit START/END range and/or a user-specified
* SMPL range. Also checks for NA's in the sample.
* Revisions by Tom Maycock, Estima
* 07/2005 Substantial modifications to clean up code by Tom Doan, Estima.
* Procedure renamed to BKFILTER; DIFF and AR options renamed
*
* References:
*
* Baxter and King (1999), "Measuring business cycles: Approximate
* band-pass filters for economic time series", Review of Economics and
* Statistics, vol 81, no. 4, pp. 575-593.
*
* James H. Stock and Mark W. Watson(1999), "Business Cycle Fluctuations
* in U.S. Macroeconomic Time Series," in Woodford and Taylor, eds.
* "Handbook of Macroeconomics."
*
*******************************************************************************
*******************************************************************************
procedure bkpad series start end padded
*
* -- procedure to pad the series front and back --
*
type series series
type integer start end
type series *padded
*
* Options passed through from the main procedure
*
option integer nma 12
option integer ar 4
option switch diff 1
*
local integer startl endl
local series work flipped
local equation fcast
*
* Copy in actuals and shift location up nma slots
*
compute startl=start+nma,endl=end+nma
move series start end padded startl
*
* Transform dependent variable, if needed and run the AR regression
*
set work startl endl = %if(diff,padded-padded{1},padded)
linreg(noprint,define=fcast,entries=1+6*diff) work startl endl
# constant work{1 to ar}
*
* Recast the equation in terms of padded
*
modify fcast
if diff
vreplace work with padded diff 1
else
vreplace work with padded
*
* Pad out end of series
*
uforecast(equation=fcast) padded endl+1 endl+nma
*
* Pad out beginning of series. Flip the series into work, estimate a new
* representation, forecast out, then flip the results back into padded
*
set flipped startl endl = padded(endl+startl-t)
set work startl endl = %if(diff,flipped-flipped{1},flipped)
linreg(noprint,define=fcast,entries=1+6*diff) work startl endl
# constant work{1 to ar}
modify fcast
if diff
vreplace work with flipped diff 1
else
vreplace work with flipped
uforecast(equation=fcast) flipped endl+1 endl+nma
set padded startl-nma startl-1 = flipped(endl+startl-t)
end bkpad
*******************************************************************************
*
* Procedure called from outside
*
procedure bkfilter series start end outser
*
type series series
type integer start end
type series *outser
*
option integer upper 2
option integer lower 32
option integer nma 12
option integer ar 4
option switch diff 1
*
local real omubar omlbar lam
local integer kk i startl endl total
local vect avec
local series padded
local real missing
if .not.%defined(series).or..not.%defined(outser) {
display "syntax: @bkfilter series start end cyclical-component"
return
}
inquire(series=series) startl<>missing
if missing>0 {
display "Filter cannot be applied because the sample range contains missing values."
display "Change sample range or interpolate missing values."
return
}
compute total=%nobs
*
* -- pad the series -- put result into padded, with entries shifted
* up by NMA slots to avoid negative subscripts
*
@bkpad(nma=nma,ar=ar,diff) series startl endl padded
*
* -- implied frequencies --
*
compute omubar=2*%pi/upper, omlbar=2*%pi/lower
if lower >= total
{
compute omlbar=0
}
/*
-------------------------------------------------------------------------
To construct a low pass filter, with a cutoff frequency of "ombar",
we note that the transfer function of the approximating filter
is given by:
alpha(om) = a0 + a1 cos(om) + ... aK cos(K om)
and the ak's are given by:
a0 = ombar/(pi)
ak = sin(k ombar)/(k pi)
where ombar is the cutoff frequency.
We employ the fact that a bandpass filter is the difference between two
low pass filters,
bp(L) = bu(L) - bl(L)
with bu(L) being the filter with the high cutoff point and bl(L) being
that with the low cutoff point.
-------------------------------------------------------------------------
*/
dimension avec(2*nma+1)
ewise avec(i)=(kk=(i-(nma+1))),%if(kk==0,(omubar-omlbar)/%pi,(sin(kk*omubar)-sin(kk*omlbar))/(kk*%pi))
/*
-------------------------------------------------------------------------
Impose constraint that transfer is
(i) 0 at om = 0 if oml > 0;
(ii) 1 at om = 0 if oml = 0;
This amounts to requiring that weights sum to zero.
-------------------------------------------------------------------------
*/
*
* Initial sum of weights:
*
compute lam=%sum(avec)
*
* -- lam = amount to add to each weight to get sum to add to zero --
*
if omlbar > 0.00000001
{
compute lam = -lam/(2*(nma+1))
}
else
{
compute lam = (1-lam)/(2*(nma+1))
}
ewise avec(i)=avec(i)+%if(i==nma+1,2.0*lam,lam)
*
* -- construct output series --
*
clear outser
filter padded startl+nma endl+nma outser startl
# -nma to nma
# avec
end bkfilter