Ideal Band Pass Filter For Stationary/Non-Stationary Series

For posting your own programs to share with others

Moderators: EViews Gareth, EViews Moderator

souliaris
Posts: 7
Joined: Thu Oct 29, 2009 5:17 am

Ideal Band Pass Filter For Stationary/Non-Stationary Series

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

terrya
Posts: 107
Joined: Wed Aug 26, 2009 2:37 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

souliaris wrote:
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

I just tried this and found that it didn't save the filtered series. I'm using EV 7 but that shouldn't make a difference I think.

EViews Gareth
Fe ddaethom, fe welon, fe amcangyfrifon
Posts: 12295
Joined: Tue Sep 16, 2008 5:38 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

I think there are a few bugs in the implementation of this code when entering periods, rather than fractions of pi. The code:

Code: Select all

if (s > 1) then  'Allow the user to use periods when calling the routine.
s1 = 2/e
e   =  2/s
s   = 2/e
endif

is definitely wrong, since s1 has not been declared anywhere. It isn't clear to me what this code is meant to be doing, so I can't quite see how to fix it. I think perhaps the (s>1) condition should also be s>=1, since it seems to me that you should be able to enter s=1 as a valid starting period.

Entering fractions of pi rather than periods should mean that this part of the code is skipped. However the couple of tries I had at entering fractions caused error messages elsewhere to appear...

terrya
Posts: 107
Joined: Wed Aug 26, 2009 2:37 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

QMS Gareth wrote:I think there are a few bugs in the implementation of this code when entering periods, rather than fractions of pi. The code:

Code: Select all

if (s > 1) then  'Allow the user to use periods when calling the routine.
s1 = 2/e
e   =  2/s
s   = 2/e
endif

is definitely wrong, since s1 has not been declared anywhere. It isn't clear to me what this code is meant to be doing, so I can't quite see how to fix it. I think perhaps the (s>1) condition should also be s>=1, since it seems to me that you should be able to enter s=1 as a valid starting period.

Entering fractions of pi rather than periods should mean that this part of the code is skipped. However the couple of tries I had at entering fractions caused error messages elsewhere to appear...

Thanks

souliaris
Posts: 7
Joined: Thu Oct 29, 2009 5:17 am

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

There were a few bugs in the original code that have (hopefully) been fully addressed in the attached version (fdfilter.prg).

First, as Gareth noted, the initial code checking for "periods" rather than fractions was incorrect (s1 was not declared).
Second, "p!" (invalid syntax) was replaced with "!p" in subroutine dfte().

You can test the revised code on the "Demo.wf1" workfile that is sent out with EViews.

The correct output for the "gdpr" series in the "Demo" dataset is provided in gdp.xls file (attached). The test run (see the last few statements of fdfilter.prg) calculates the FD filter for (6,32) periods -- the usual setting for US real GDP.

Postscript: Gareth pointed out that the s>=1 condition didn't make sense given the subsequent logic for calculating the interval over (0,pi). The correct condition is s>=2. Code has been revised to show s>=2.
Attachments
fdfilter.prg
Eviews code for Ideal Band Pass Filter (Frequency Domain).
gdp.xls
Expected output from the FD filter for "GDPR", which appears in the "DEMO.WF1" demonstration file sent out with Eviews (Chapter 2).
Last edited by souliaris on Thu Jan 21, 2010 3:48 pm, edited 1 time in total.

EViews Gareth
Fe ddaethom, fe welon, fe amcangyfrifon
Posts: 12295
Joined: Tue Sep 16, 2008 5:38 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

This does seem to be an improvement. I do have one further note though.

I'll admit this is not my area of expertise, so I might be asking silly questions here, but can the start period be equal to 1? I note that you changed your program to allow s>=1, which might indicate that s=1 is a valid input. BUT, if s=1, then e=2/1=2, and you don't allow e>1, so the program still fails...

souliaris
Posts: 7
Joined: Thu Oct 29, 2009 5:17 am

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

Thanks Gareth. Code has been changed in the original post. The correct condition is (s>=2), allowing a maximum e (in terms of fraction of [0 to pi]) of 1. Apologies. Sam

EViews Gareth
Fe ddaethom, fe welon, fe amcangyfrifon
Posts: 12295
Joined: Tue Sep 16, 2008 5:38 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

Great. For EViews 7 users, I've added a graphical user interface to Sam's program. You'll have to download both my interface program, and Sam's program, since mine is an add-on to his, not a replacement. Note you'll have to remove the two "call" lines at the bottom of Sam's program when running my add-on program.
Attachments
fdfilterui.prg
User Interface for FDFILTER.prg

terrya
Posts: 107
Joined: Wed Aug 26, 2009 2:37 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

QMS Gareth wrote:Great. For EViews 7 users, I've added a graphical user interface to Sam's program. You'll have to download both my interface program, and Sam's program, since mine is an add-on to his, not a replacement. Note you'll have to remove the two "call" lines at the bottom of Sam's program when running my add-on program.

I've run the programs and they're great. What would be greater if you could add the GUI to allow a number of series to be filtered in a way that's similar to the Stock prog. I'm doing a study of the results of a number of different filters that involves at least 6 different data series.

EViews Gareth
Fe ddaethom, fe welon, fe amcangyfrifon
Posts: 12295
Joined: Tue Sep 16, 2008 5:38 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

Simple enough.
Attachments
fdfilterui.prg

terrya
Posts: 107
Joined: Wed Aug 26, 2009 2:37 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

QMS Gareth wrote:Simple enough.

Hi
Just run the program. Stops with error msg: "Syntax error in Series {%S_name} = resids. I was using 2 series.

terrya
Posts: 107
Joined: Wed Aug 26, 2009 2:37 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

terrya wrote:
QMS Gareth wrote:Simple enough.

Hi
Just run the program. Stops with error msg: "Syntax error in Series {%S_name} = resids. I was using 2 series.

Oops! My mistake. Didn't read the instructions properly. Works alright.

randomlearning
Posts: 8
Joined: Thu Jun 03, 2010 5:28 am

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

Is there a way of using this filter in a loop?

I am trying:

Code: Select all

for %i a b c
call IdealBandPass({%i},0.0625,0.33333,"ibf_{%i}","1950q1","2070q1")
next

but it gives me an error that output name that was specified is reserved.

EViews Gareth
Fe ddaethom, fe welon, fe amcangyfrifon
Posts: 12295
Joined: Tue Sep 16, 2008 5:38 pm

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

try:

Code: Select all

for %i a b c
%temp = "ibf_" + %i
call IdealBandPass({%i},0.0625,0.33333,%temp,"1950q1","2070q1")
next

randomlearning
Posts: 8
Joined: Thu Jun 03, 2010 5:28 am

Re: Ideal Band Pass Filter For Stationary/Non-Stationary Series

Thanks. Now I'm getting the error "Cannot assign string expression to numeric variable in "SERIES {%S_NAME}=resid".