Implements the IdealBandPass filter as described in Corbae and Ouliaris (2006), "Extracting Cycles From Non Stationary Data", Econometric Theory and Practice, pp. 167-177, Cambridge University Press.
Involves detrending the series in the frequency domain. Series can be stationary, or it can have a unit root.
It tends to produce output similar to the B-K filter, though without the end-point issue (estimates end-points directlty).
Subroutine DFT() computes the discrete fourier transform a series.
Subroutine INVDFT() computed the inverse of the discrete fourier transform. DFT() needs to be called first.
Subroutine DFTE() annihilates appropriate frequencies according to user requirements
Subroutine IdealBandPass() implements the routine (at the end of the code section).
The arguments to this routine are:
'x: series name
's: starting period (can be a fraction of pi, or the number of periods)
'e: ending period (can be a fraction of pi, or the number of periods)
'%s_name = name of the filtered series
'%s_date -- start from this date
'%e_date -- end at this date.
Last two arguments are for series that having missing values at the beginning and end of the workfile range.
'Example call: call IdealBandPass(r_apsp,0.25,0.5,"i_bp_oil","1980m05", "2009m05")
in which: r_apsp is the actual name of the series in the active workfile. "i_bp_oil" will be the name of the filtered series in the same workfile.
Sam Ouliaris, October 29, 2009
souliaris@imf.org
Code: Select all
subroutine dft(series x, series y, string %s_date, string %e_date)
' Subroutine to compute the DFT of a (possibly complex) variable.
' Y is mostly set to 0.
' Output is saved in two temporary series: dft_r (real_part), dft_m (imaginary part)
' Sam Ouliaris, October 29, 2009, souliaris@imf.org
series dft_r = 0.0
series dft_m = 0.0
!TwoPi = 2 * 3.141592653589793238462
!n = @dtoo(%e_date)-@dtoo(%s_date)+1
!p = @dtoo(%s_date)-1
for !i = 0 to !n-1
!ae = -!TwoPi*(!i/!n)
!rp=0.0
!cp = 0.0
for !k = 1 to !n
!a = Cos((!k-1)*!ae)
!b = Sin((!k-1)*!ae)
!rp = !rp + ((x(!k+!p)*!a)-(y(!k+!p)*!b))
!cp = !cp + ((x(!k+!p)*!b)+(y(!k+!p)*!a))
next
dft_r(!i+!p+1) = !rp
dft_m(!i+!p+1) = !cp
next
series dft_r = dft_r/!n
series dft_m = dft_m/!n
endsub
subroutine invdft(string %s_date, string %e_date)
'Subroutine to compute the inverse of a DFT.
'Result is saved in two series dft_inv_r (real part)
'dft_inv_m (complex part).
!TwoPi = 2*3.141592653589793238462
series dft_inv_r=0.0
series dft_inv_m=0.0
!n = @dtoo(%e_date)-@dtoo(%s_date)+1
!p = @dtoo(%s_date)-1
for !i = 0 to !n-1
!ae = !TwoPi*(!i/!n)
!rp=0.0
!cp = 0.0
for !k = 1 to !n
!a = Cos((!k-1)*!ae)
!b = Sin((!k-1)*!ae)
!rp = !rp + ((dft_r(!k+!p)*!a)-(dft_m(!k+!p)*!b))
!cp = !cp + ((dft_r(!k+!p)*!b)+(dft_m(!k+!p)*!a))
next
dft_inv_r(!i+!p+1) = !rp
dft_inv_m(!i+!p+1) = !cp
next
endsub
subroutine dfte(series x, scalar s, scalar e, string %s_name,string %s_date, string %e_date)
' Extracts the relevant frequencies between s, e [inclusive], assuming s, e are between 0 and 1.
!n = @dtoo(%e_date)-@dtoo(%s_date)+1
series ibp_zero=0.0
call dft(x,ibp_zero,%s_date,%e_date)
call rm("ibp_zero")
!p = @dtoo(%s_date)-1
!k = 0
if (@mod(!n,2) = 0.0) then
!m = !n / 2
else
!m = (!n+1)/2
endif
for !j = 1 to !m
!k = !j/!m
if ((!k < s) or (!k > e)) then
if !j = 1 then
dft_r(!p+1) = 0.0
dft_m(!p+1) = 0.0
else
dft_r(!j+!p)=0.0
dft_r(!n-(!j-2)+!p)=0.0
dft_m(!j+!p)=0.0
dft_m(!n-(!j-2)+!p)=0.0
endif
endif
next
if (s > 0.0) then 'This is to ensure that the 0 frequency is always eliminated for small N
dft_r(!p+1) = 0.0
dft_m(!p+1) = 0.0
endif
if (@mod(!n,2)=0.0) then
if (e < 1.0) then
dft_r(!m+1+!p) = 0.0
dft_m(!m+1+p!) = 0.0
endif
endif
call invdft(%s_date,%e_date) 'Inverse; result goes back to dft_inv_r, dft_inv_e
series {%s_name} = dft_inv_r
'Clean up...
call rm("dft_inv_r")
call rm("dft_inv_m")
call rm("dft_r")
call rm("dft_m")
endsub
subroutine rm(string %s_name)
if (@isObject(%s_name) = 1) then
delete {%s_name}
endif
endsub
' Examplel: call IdealBandPass(r_gdp,6,32,"Ibp_r_gdp", "start_date", "end_date")
'r_gdp = name of the series
'6,32: extract cycle between 6 and 32 periods. Can also be called with fractions" 1/16=2/32, 1/3=2/6
subroutine IdealBandPass(series x, scalar s, scalar e, string %s_name,string %s_date, string %e_date)
'Implements the IdealBandPass filter, as in Corbae and Ouliaris (2006), "Extracting Cycles From Non Stationary Data", Econometric Theory and Practice, pp. 167-177, Cambridge University Press.
'Series can be stationary, or it can have a unit root.
'Calling arguments:
'x: series name
's: starting period (can be a fraction of pi, or the number of periods)
'e: ending period (can be a fraction of pi, or the number of periods)
'%s_name = name of the filtered series
'%s_date -- start from this date
'%e_date -- end at this date.
'Last two arguments are for series that having missing values at the beginning and end of the workfile range.
if (s > 1) then 'Allow the user to use periods when calling the routine.
s1 = 2/e
e = 2/s
s = 2/e
endif
if s <=0 or (s>e) then
stop
endif
if e <=0 or (e>1) then
stop
endif
smpl %s_date %e_date
series dfft_gg = (@trend+1)/@obs(x) 'Scaled time trend.
call dfte(dfft_gg,s,e,"dfft_st",%s_date, %e_date) 'Compute the DFT of gg, and extract the relative frequencies.
call dfte(x,s,e,"dfft_sx",%s_date, %e_date) 'Compute the DFT of x, and extract the relative frequencies.
equation dft_yy
freeze(mode=overwrite,dfft_ibp_reg) dft_yy.ls dfft_sx dfft_st 'Run the regression in the time domain.
series {%s_name} = resid 'Residuals contain the ideal band pass filter.
call rm("dfft_sx")
call rm("dfft_st")
call rm("dfft_gg")
call rm("dfft_ibp_reg")
call rm("dft_yy")
endsub