*! version 1.2.4    22Dec2001       C F Baum/V Wiggins
*  Implement Robinson Log-Periodogram Regression, Annals of Statistics 1995
*  for univariate and multivariate cases
*  1.2.0: modified output, handles multiple powers for single var, constraints
*  1.2.1: published STB-57
*  1.2.2: trap j<1, add varname to output
*  1.2.3: observations miscounted when variables started at different periods
*  1.2.4: L option did not work with multiple variables

program define roblpr6, rclass
	version 6

	syntax varlist(ts) [if] [in] [ , Constraints(string) J(integer 1) /*
		*/ L(integer 0) Powers(string) REGRESSONLY ]

	_rmcoll `varlist' `if'
	local varlist `r(varlist)'
	local orgvlst `varlist'
	tsrevar `varlist'
	local varlist `r(varlist)'

	local nrvars : word count `varlist'
	if `nrvars' == 1 & "`powers'" != "" {
		numlist "`powers'", range(>0 < 1)
		local powers `r(numlist)'
	}
	else if `nrvars' > 1 & "`powers'" != "" {
		capture numlist "`powers'", range(>0 < 1) min(1) max(1)
		if _rc {
			di in red "with multiple variables option "	/*
				*/ "powers(#) accepts only a single # and it"
			di in red "must be between 0 and 1"
			exit 198
		}
		local powers `r(numlist)'
	}
	else	local powers .9

	local powct : word count `powers'
   	marksample touse

	if `j' < 1 {
		di in red "Error: j must be at least 1"
		exit 198
		}
			/* get time variables */
	_ts timevar, sort
	markout `touse' `timevar'
	tsreport if `touse', report
	if r(N_gaps) {
		di in red "sample may not contain gaps"
		exit
	}

	tempname xr xi n pg lpg ok vok vcv bs
	tempvar var vecy

	if `nrvars' != 1 | "`regressonly'" != "" {
		preserve
		gen long `n' = sum(`touse')-1
		tokenize `varlist'
		local i 1
		local last 0
	 	/* ensure even number of ordinates */
		qui count if `touse'
		local em=int(r(N)^`powers')+1
		if mod(`em',2)!=0 {local em=`em'+1}
		/* calculate number of obs needed, adjust if necessary */
		local newobs = `em'*`nrvars'
		if `nrvars' > 1  {
			qui count
			if `newobs' > r(N) { qui set obs `newobs' }
		}
		qui replace `touse'=0 if `touse'==.
		qui gen double `vecy' =.
		qui gen double `var'=.
		qui gen `ok'=_n
		local xlist "`vecy'"
* 1.2.4
		local sofar 1
		while `i'<=`nrvars' {
			local last=`last'+`em'
			local first=`last'-`em'+1
			local offset=(`i'-1)*`em'
			capt drop `xr' `xi' `lpg'
			qui replace `var'=``i''
			fft `var' if `touse', gen(`xr' `xi')
			if `j'==1 {
				qui gen double `lpg' = log((`xr'^2+`xi'^2)) /*
					*/ if `touse'
				}
			else 	{
				/* moving sum of periodogram over J elements */
				qui gen double `pg' = (`xr'^2+`xi'^2) if `touse'
				local ii=1
				qui gen double `lpg' = `pg'
				while `ii'<`j' {
					qui replace `lpg'=`lpg'+`pg'[_n-`ii']
					local ii=`ii'+1
					}
				qui replace `lpg'=log(`lpg')
				qui gen `vok'=(mod(_n,`j')==0)
				drop `pg'
				}
			qui replace `vecy'=`lpg'[_n-`offset'] 		/*
				*/ in `first'/`last'
			if `j'>1 {
				qui replace `ok'=`vok'[_n-`offset'] 	/*
					*/ in `first'/`last'
				drop `vok'
				}
			tempvar llam iota
			qui gen double `llam'=0 in 1/`newobs'
			qui gen double `iota'=0 in 1/`newobs'
			qui replace `llam' = 2.0 * 			   /*
				*/ log( 2.0 * _pi*(`n'[_n-`offset'])/`em') /*
				*/ in `first'/`last'
			qui replace `iota' = 1.0 in `first'/`last'
			/* implement L option via iota; adjust by one
			 * for zero freq, interact with J */
			if `l'>0 {
* 1.2.4
				local hi =`l'*`j'+`sofar'
				local lo = `sofar'+1
				qui replace `iota'=. in `lo'/`hi'
				local sofar = `sofar' + `em'
***				local hi =`l'*`j'+1
***				qui replace `iota'=. in 2/`hi'
				}
*		noi di in r "sofar, lo, hi `sofar' `lo' `hi'"
			drop ``i''
			rename `llam' ``i''
			local llam "``i''"
			local xlist `xlist' ``i'' `iota'
			local blist `blist' ``i''
			local i=`i'+1
			}

			/* log periodogram regression: use noc to get
			 * indiv constant terms */
		if "`constraints'" != "" {
			cons def 988 _cons = 0
			qui cnsreg `xlist' if `ok'>0 , 	/*
				*/ constraint(`constraints' 988)
		}
		else	qui regress `xlist' if `ok'>0, noc


		if "`regressonly'" != "" {
			return local depvar = "`orgvlst'"
			return scalar nord = e(N)
			return scalar rob = -_b[`llam']
			mat `vcv'=e(V)
			return scalar se = sqrt(`vcv'[1,1])
			return scalar t = return(rob)/return(se)
			return scalar p = tprob(`em',return(t))
			return scalar power = `powers'
			exit
		}
	}


	if `nrvars'==1 {
		di in gr _n "Robinson estimates of fractional differencing parameter for " /*
		*/ in ye "`orgvlst'"
		if `l' != 0 {
			di in gr _col(35) "Skipped ords" _col(48) "="	/*
				*/ in ye %7.0f `l'
		}
		if `j' != 1 {
		   di in gr _col(35) "Averaging" _col(48) "=" in ye %7.0f `j'
		}

		di in gr _dup(55) "-"
		di in gr "Power   Ords      Est d    Std Err  "		/*
			*/ "t(H0: d=0)    P>|t|"
		di in gr _dup(55) "-"

		if "`l'" != "" { local optbase l(`l') }
		if "`j'" != "" { local optbase `optbase' j(`j') }
		tokenize `powers'
		local i 1
		while "``i''" != "" {
			capture roblpr6 `varlist' if `touse',	/*
				*/ regressonly `optbase' power(``i'')
			if _rc {
				di in blue "  roblpr could not be " 	/*
					*/ " calculated for power = ``i''"
			}
			else {
				di in gr   " "    %4.2g ``i'' in ye  	/*
					*/ " "    %6.0f r(nord)		/*
					*/ "  "   %9.0g r(rob)		/*
					*/ "  "   %9.0g r(se) 		/*
					*/ "  "   %9.4f r(t)		/*
					*/ "  "   %8.3f r(p)
			}
			local i = `i' + 1
		}

		di in gr _dup(55) "-"
		return local depvar = "`orgvlst'"
		return scalar p = r(p)
		return scalar se = r(se)
		return scalar t = r(t)
		return scalar rob = r(rob)
		return scalar nord = r(nord)
		return scalar power = r(power)
	}
	else {
		return scalar nord = e(N)/`nrvars'
		di in gr _n "Robinson estimates of fractional differencing parameters"
		di in gr "Power = " %4.3g `powers'			/*
		   */ _col(41) "Ords" _col(54) "=" in ye %7.0f return(nord)
		if `l' != 0 {
			di in gr _col(41) "Skipped ords" _col(54) "="	/*
				*/ in ye %7.0f `l'
		}
		if `j' != 1 {
		   di in gr _col(41) "Averaging" _col(54) "=" in ye %7.0f `j'
		}
		di in gr _dup(61) "-"
		di in gr "Variable         |     Est d       Std Err      " /*
			*/ "t       P>|t|"
		di in gr "-----------------+-------------------------"	/*
			*/ "------------------"
		local i=1
		local j=1
		mat `bs'=e(b)
		mat `vcv'=e(V)
		tokenize `orgvlst'
		return local depvar = "`orgvlst'"
		while `i'<=`nrvars' {
			return scalar rob`i' = -`bs'[1,`j']
			return scalar se`i' = sqrt(`vcv'[`j',`j'])
			return scalar t`i' = return(rob`i')/return(se`i')
			return scalar p`i' = tprob(`em',return(t`i'))
			di in gr "``i''" _column(18) "|" in ye	 	/*
				*/ _column(20)  %9.0g -`bs'[1,`j']	/*
				*/ _column(34) %9.0g return(se`i') 	/*
				*/ _column(42) %9.4f return(t`i') 	/*
				*/ _column(56) %6.3f return(p`i')
			local i=`i'+1
			local j=`j'+2
			}
* logic from testparm v3.1
		tokenize `blist'
		local lhs "`1'"
        	quietly {
                	mac shift
                	test `1'=`lhs', notest
                	mac shift
                	while ("`1'"!="") {
                        	test `1'=`lhs', accum notest
                        	macro shift
                		}
        		test
        		}
		return scalar F = r(F)
		return scalar df = r(df)
		return scalar df_r = r(df_r)
		return scalar pF = r(p)
		di in gr _dup(61) "-"
		if r(df) > 0 {
		   di in gr "Test for equality of d coefficients:" _continue
		   di in gr "   F(" in ye r(df) in gr "," in ye r(df_r)	/*
			*/ in gr  ") = " in ye %7.0g r(F)		/*
			*/ in gr "   Prob > F = " in ye %6.4f r(p)
		}

	}
end

exit