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 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. if (s >= 2) then 'Allow the user to use periods when calling the routine. !s1 = e e = 2/s s = 2/!s1 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 call IdealBandPass(gdpr, 0.0625,0.3333,"f_gdp","1952q1","2003q4") 'This should give the same result as above --- here we use periods to define the range. call IdealBandPass(gdpr, 6,32,"p_gdp","1952q1","2003q4")