Page 1 of 2

Ideal Band Pass Filter For Stationary/Non-Stationary Series

Posted: Thu Oct 29, 2009 11:43 am
by souliaris
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

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

Posted: Sat Jan 16, 2010 5:05 pm
by terrya
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.

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

Posted: Mon Jan 18, 2010 2:19 pm
by EViews Gareth
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...

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

Posted: Mon Jan 18, 2010 2:57 pm
by terrya
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

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

Posted: Thu Jan 21, 2010 3:02 pm
by souliaris
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.

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

Posted: Thu Jan 21, 2010 3:10 pm
by EViews Gareth
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...

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

Posted: Thu Jan 21, 2010 3:54 pm
by souliaris
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

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

Posted: Thu Jan 21, 2010 4:02 pm
by EViews Gareth
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.

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

Posted: Fri Jan 22, 2010 12:29 am
by terrya
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.

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

Posted: Fri Jan 22, 2010 3:27 pm
by EViews Gareth
Simple enough.

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

Posted: Fri Jan 22, 2010 4:49 pm
by terrya
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.

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

Posted: Fri Jan 22, 2010 4:57 pm
by terrya
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.

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

Posted: Tue Jun 15, 2010 10:43 am
by randomlearning
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.

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

Posted: Tue Jun 15, 2010 10:50 am
by EViews Gareth
try:

Code: Select all

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

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

Posted: Tue Jun 15, 2010 11:12 am
by randomlearning
Thanks. Now I'm getting the error "Cannot assign string expression to numeric variable in "SERIES {%S_NAME}=resid".