*******************************************************************************
* @BPFILTER(OPTIONS) SERIES START END OUTSER
*
* 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)
*
*     [PADDIFF]/
*       NOPADDIFF    padding uses AR on 1st differences [levels]
*
*     ARPAD  = [4]   AR order for padding
*
*                    (padding adds artifical data at start
*                    and end of series prior to permit
*                    moving-average method of filtering
*                    using AR backcasts and forecasts)
*
*
* BPFILTER execeutes a bandpass filter on a series using the method
* of Baxter and King (1995). 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 cenetered 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 observationsusing AR backcasts and
* forecasts. ARPAD controls the order of this AR.
*
* [Corollary: Missing data is a big problem! Interpolate gaps
* before using BPFILTER.]
*
* 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 (1998) and on Baxter and King's (1995) settings.
* The file DEMO.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
*
* Alan Taylor's RATS translation 12/98.
* 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.
*
* Revision History:
*  9/22/1999: Modifications made to BPPAD 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
*
* References:
*
* Marianne Baxter and Robert G. King, "Measuring Business Cycles:
* Approximate Band-Pass Filters for Economic Time Series," NBER
* Working Paper No. 5022, February 1995.
*
* James H. Stock and Mark W. Watson, "Business Cycle Fluctuations
* in U.S. Macroeconomic Time Series," NBER Working Paper No. 6528,
* April 1998.
*
* Alan M. Taylor                                            mailto:amt@nwu.edu
* Department of Economics              http://www.econ.nwu.edu/faculty/taylor/
* Northwestern University
* 2003 Sheridan Road                                          tel:847-491-8234
* Evanston IL 60208-2600                                      fax:847-491-7001
*
*******************************************************************************

*******************************************************************************
* -- procedure to pad the series front and back --
procedure  bppad  series  start  end  xvec
*
type series series
type integer start  end
type rectangular *xvec
*
option integer nma 12
option integer arpad 4
option switch paddiff 1
*
local integer startl endl kk
local series hold hold1 hold2 trend
local equation fcast bcast
*

dimension xvec(end-start+1+2*nma,1)

* -- given data --

do kk = start,end
  compute xvec(kk-start+nma+1,1)=series(kk)
end do

* -- pad forward with forecast --

clear hold hold1 hold2
* set hold = series
set hold start end = series

if paddiff == 1
{
* Note: the SET is done over full range, so that
*       first entry of HOLD is set to NA, as desired.
*diff hold / hold1; set hold = hold1
diff hold start+1 end hold1; set hold start end = hold1
linreg(noprint,define=fcast) hold start+1 end
# constant hold{1 to arpad}
}
else
{
linreg(noprint,define=fcast) hold start end
# constant hold{1 to arpad}
}

forecast 1 nma;# fcast hold2 1

if paddiff == 1
{
accumulate hold2
do kk = 1,nma
  compute xvec(kk+end-start+1+nma,1)=hold2(kk)+xvec(end-start+1+nma,1)
end do
}

if paddiff == 0
{
do kk = 1,nma
  compute xvec(kk+end-start+1+nma,1)=hold2(kk)
end do
}

* -- pad back with backcast (reverse the series) --

clear hold hold1 hold2
set hold start end = series
set trend start end = t
order(decreasing) trend start end hold

if paddiff == 1
{
* Note: the SET is done over full range, so that
*       first entry of HOLD is set to NA, as desired.
diff hold start+1 end hold1; set hold start end = hold1
linreg(noprint,define=bcast) hold start+1 end
# constant hold{1 to arpad}
forecast 1 nma;# bcast hold2 1
}
else
{
linreg(noprint,define=bcast) hold start end
# constant hold{1 to arpad}
forecast 1 nma;# bcast hold2 1
}

if paddiff == 1
{
accumulate hold2
do kk = 1,nma
  compute xvec(nma+1-kk,1)=hold2(kk)+xvec(nma+1,1)
end do

}

if paddiff == 0
{
do kk = 1,nma
  compute xvec(nma+1-kk,1)=hold2(kk)
end do
}

end bppad

*******************************************************************************
* -- procedure to do the bandpass filter --
procedure  bpfilter  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 arpad 4
option switch paddiff 1
*
local real omubar omlbar lam temp
local integer kk jj startl endl total
local rectangular akvec avec xvec

dimension akvec(nma+1,1)
dimension avec(2*nma+1,1)

if .not.(%defined(series).and.%defined(outser)) {
  display 'syntax: @bpfilter  series  start  end  cyclical-component'
  return
}

inquire(series=series) startl>>start endl>>end
compute total=endl-startl+1

if (%defined(start).and.%defined(end))
{
 stats(noprint) series start end
 if %nobs<(end-start+1) {
   display 'Filter cannot be applied because the sample range'
   display 'contains missing values. Change sample range or '
   display 'interpolate missing values.'
   return
  }
}

* -- pad the series -- put result in xvec(endl-startl+1,1)
if startl<=arpad {
 display 'Supplied start date to early given ARPAD value.'
 display 'Re-setting start date to: ' %datelabel(arpad+1)
 compute startl = arpad + 1
}

@bppad(nma=nma,arpad=arpad,paddiff) series startl endl xvec

* -- 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.
 -------------------------------------------------------------------------
*/

* -- set up the akvec weights --

compute akvec = %const(0.0)
compute akvec(1,1) = (omubar-omlbar)/%pi

do kk = 1,nma
  compute akvec(kk+1,1) = (sin(kk*omubar)-sin(kk*omlbar))/(kk*%pi)
end do

/*
 -------------------------------------------------------------------------
  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=2*%sum(akvec)-akvec(1,1)

* -- 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))
   }

do kk = 1,nma+1
  compute akvec(kk,1)=akvec(kk,1)+lam
end do

compute akvec(1,1)=akvec(1,1)+lam

* -- set up the avec weights --

compute avec = %const(0.0)
compute avec(nma+1,1) = akvec(1,1)

do kk = 1,nma
  compute avec(nma+1+kk,1) = akvec(kk+1,1)
  compute avec(nma+1-kk,1) = akvec(kk+1,1)
end do

* -- construct output series --

clear outser

do kk = startl,endl
  compute temp = 0.0
  do jj = 1,2*nma+1
    compute temp = temp + xvec(kk-startl+jj,1)*avec(jj,1)
    end do jj
  compute outser(kk) = temp
end do

end bpfilter

*******************************************************************************

