*! cfitzrw  1.0.0  CF Baum/M Lopez 30jun2006
* 1.0.0: CF Mata version from
*        http://ideas.repec.org/c/wpa/wuwppr/0507001.html
*        corr for tsrevar, byable(recall)

program define cfitzrw, rclass byable(recall,noheader)
        version 9.2

syntax varlist(ts) [if] [in], PLO(string) PHI(string) STUB(string)

        marksample touse
        _ts timevar panelvar if `touse', sort onepanel
        markout `touse' `timevar'
        tsreport if `touse', report
        if r(N_gaps) {
                di in red "sample may not contain gaps"
                exit
                }
        if `plo' >= `phi' {
                di in red "Lower bound (plo) must be less than upper bound (phi)"
                error 198
                }
        if `plo' < 2     {
                di in red "Lower bound (plo) must be at least 2"
                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'
                gen float `stub'_`v'_`suf' = .
                local varlist3 "`varlist3' `stub'_`v'_`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: cfrwmata("`varlist1'","`varlist3'","`touse'",`plo', `phi')

        return local rawvars "`varlist'"
        return local filtvars "`varlist3'"
        return local phi "`phi'"
        return local plo "`plo'"
        end


mata:
// from Pawel Kowal's CF function for MATLAB
       void cfrwmata(string scalar vname,
                   string scalar vname3,
                   string scalar touse,
                   real scalar plo,
                   real scalar phi)
       {
        real scalar a, b, n2, n1, B_end, B_start,T
        real matrix X, Y
        real vector B_main, B_tmp, i, B, BB, dX, J
        string rowvector vars, vars3
        string scalar v, v3

// access the Stata variables in varlist, varlist3, honoring touse
// use st_data, not st_view, as X will be modified
        vars = tokens(vname)
        v = vars[|1,.|]
		X = st_data(.,v,touse)
        vars3 = tokens(vname3)
        v3 = vars3[|1,.|]
        st_view(Y,.,v3,touse)
     	T = rows(X)

        plo = max((2, plo))

        a = 2*pi()/phi
        b = 2*pi()/plo

// remove OLS trends from X
        dX = (X[T,.] - X[1,.]) :/ (T-1)
        X = X - ( (0::(T-1)) * dX)

        J   = (1::T)
        B   = (sin(b:*J) - sin(a:*J)):/ (pi():*J)
        B  = ( (b-a)/pi() \ B )

        B_main = B[1::T]
        B_tmp  = B[1::T]

        n2 = T-1
        n1 = 0

        for(i=1; i<=T; i++){
        	B_end  = sum(B_tmp[(1+n2)::T])
       		B_start = sum(B_tmp[(1+n1)::T])

            BB = B_main
         	BB[1] = B_start
       		BB[T] = B_end

         	Y[i,.] =  colsum(X :* BB)

       		if(i<T) {
       			B_main = (B[i+1] \ B_main[1::(T-1)])
				n2 = n2-1
				n1 = n1+1
			}
   		}
}

    end