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")