*! fdfilter 1.00 03jan2010 * Define program program define fdfilter, rclass byable(recall, noheader) version 9.2 * Syntax requires: variables, start and end points, and stub to name new variables syntax varlist(ts) [if] [in], S(string) E(string) STUB(string) * Mark the sample to use marksample touse * Verify that data have been tsset _ts timevar panelvar if `touse', sort onepanel * Exclude observations if time variable is missing markout `touse' `timevar' * Check for gaps tsreport if `touse', report * Error message is sample contains gaps if r(N_gaps) { di in red "sample may not contain gaps" exit } * Allow the user to use periods when calling the subroutine if `s'>=2 { local s1=`e' local e=2/`s' local s=2/`s1' } * Check end point larger than starting point if `s'>`e' { di in red "Starting point must be smaller than ending point" error 198 } if `s'<=0 | `e'<=0 | `e'>1 { di in red "Invalid e or s" error 198 } * validate each new varname defined by stub() local kk: word count `varlist' local varlist2: subinstr local varlist "." "_", all local suf = _byindex() qui forval i = 1/`kk' { local v: word `i' of `varlist2' confirm new var `stub'_`v'_`suf' confirm new var `stub'_`v'_sm_`suf' gen double `stub'_`v'_`suf' = . gen double `stub'_`v'_sm_`suf' = . local varlist3 "`varlist3' `stub'_`v'_`suf'" local varlist4 "`varlist4' `stub'_`v'_sm_`suf'" } * create temp vars for any ts operators in the varlist * pass the resulting varlist1 to Mata fn tsrevar `varlist' local varlist1 `r(varlist)' mata: idbp("`varlist1'","`varlist3'","`varlist4'","`touse'",`s',`e') return local rawvars "`varlist'" return local filtvars "`varlist3'" return local s "`s'" return local e "`e'" end mata /* Function to compute Fourier transform */ complex matrix dft( transmorphic matrix x) { real matrix yr, yi, dftr, dfti, range real scalar n, c, i, ae complex matrix df yr=Re(C(x)) yi=Im(C(x)) n=rows(x) c=cols(x) dftr=J(n,c,.) dfti=J(n,c,.) range=range(1,n,1) ae=-2*pi()/n for (i=0; i<=n-1; i++) { dftr[i+1,.]=colsum(yr:*cos((range:-1):*ae:*i)-yi:*sin((range:-1):*ae:*i)) dfti[i+1,.]=colsum(yr:*sin((range:-1):*ae:*i)+yi:*cos((range:-1):*ae:*i)) } df=C(dftr,dfti) return(df/n) } /* Function to compute the inverse of Fourier transform */ complex matrix inv_dft(complex matrix x) { complex matrix df real matrix dftr, dfti, dftinvr, dftinvi, range real scalar n, c, ae, i n=rows(x) c=cols(x) dftr=Re(x) dfti=Im(x) dftinvr=J(n,c,.) dftinvi=J(n,c,.) range=range(1,n,1) ae=2*pi()/n for (i=0; i<=n-1; i++) { dftinvr[i+1,.]=colsum(dftr:*cos((range:-1):*ae:*i)-dfti:*sin((range:-1):*ae:*i)) dftinvi[i+1,.]=colsum(dftr:*sin((range:-1):*ae:*i)+dfti:*cos((range:-1):*ae:*i)) } df=C(dftinvr,dftinvi) return(df) } /* Function to extract the relevant frequencies */ real matrix dfte(transmorphic matrix x, real scalar s, real scalar e) { real matrix dftr, dfti real scalar n,m, j, c, k complex matrix df, dfte n=rows(x) c=cols(x) df=dft(C(x)) dftr=Re(df) dfti=Im(df) m=(n+mod(n,2))/2 for (j=1; j<=m; j++) { k=j/m if (ke) { if (j==1) { dftr[j,.]=J(1,c,0) dfti[j,.]=J(1,c,0) } else { dftr[j,.]=J(1,c,0) dftr[n-(j-2),.]=J(1,c,0) dfti[j,.]=J(1,c,0) dfti[n-(j-2),.]=J(1,c,0) } } } if (s>0) { dftr[1,.]=J(1,c,0) dfti[1,.]=J(1,c,0) } if (mod(n,2)==0 & e<1) { dftr[m+1,.]=J(1,c,0) dfti[m+1,.]=J(1,c,0) } dfte=Re(inv_dft(C(dftr,dfti))) return(dfte) } void idbp(string scalar vname, string scalar vname3, string scalar vname4, string scalar touse, real scalar s, real scalar e) { real matrix trend, t, dfftx, M, X, Y, Z string rowvector vars, vars3, vars4 string scalar v,v3,v4 real scalar i vars = tokens(vname) v = vars[|1,.|] X=st_data(., v,touse) vars3 = tokens(vname3) v3 = vars3[|1,.|] st_view(Y,.,v3,touse) vars4 = tokens(vname4) v4 = vars4[|1,.|] st_view(Z,.,v4,touse) trend=range(1,rows(X),1)/rows(X) t=dfte(trend,s,e) dfftx=dfte(X,s,e) M=I(rows(X),rows(X))-t*invsym(t'*t)*t' for (i=1; i<=cols(X); i++) { Y[.,i]=M*dfftx[.,i] } Z[.,.]=X-Y } end