*
* @AGFRACTD(options) series  start  end
* Estimates fractional difference power for series using the bias-reduced
* technique from Andrews and Guggenberger.
*
* Options:
*     M  = Number of frequencies to use in running the regression [sqrt(number of obs)]
*     R  = Number of even powers of the frequencies to use [R=0 gives Geweke-Porter-Hudak]
*     [PRINT]/NOPRINT
*
*     At this point, this isn't set up to choose M using Andrews' formulas.
*
* Reference:
*    Andrews and Guggenberger, "A Biased-Reduced Log-periodogram Regression
*    Estimator for the Long-Memory Parameter", Econometrica 2003, vol 71, no. 2,
*    pp 675-712.
*
* Revision Schedule:
*    08/2003 Written by Tom Doan, Estima
*
procedure agfractd series  start  end
type series  series
type integer start  end
*
options integer m
options integer r  0
options switch  print 1
*
local integer startl endl ords i j
local integer obs
local series diff logpdgm
local vector[series] freqs
local symm gamma
local rect xi
local vect mu c tau
*
inquire(series=series) startl<<start endl<<end
*
compute startl = startl + 1 , obs = endl-startl+1
diff series startl endl diff
*
freq 2 obs
rtoc startl endl 1
# series
#   1
*
fft 1
cmult 1 1
cln 1
*
* Define the quantities needed for computing plug-in m
*
dim gamma(r+1,r+1) xi(r+1,r+1) mu(r+1) c(r+1) tau(r+1)
ewise gamma(i,j)=4.0*i*j/((2*i+2*j+1)*(2*i+1)*(2*j+1))
ewise mu(i)=2.0*i/(2*i+1)**2
ewise xi(i,j)=2.0*j*(3+2*i)/((2*i+2*j+3)*(2*j+1))
ewise c(i)=1.0/(1-%qform(inv(%xsubmat(gamma,1,i,1,i)),%xsubvec(mu,1,i)))
ewise tau(i)=-(2*%pi)**(2+2*i)*(2+2*i)*c(i)/(2*%factorial(3+2*i)*(3+2*i)) * $
  (1-%dot(%xsubvec(mu,1,i),inv(%xsubmat(gamma,1,i,1,i))*tr(%xsubmat(xi,i,i,1,i))))
*
if %defined(m)
   compute ords = m
else
   compute ords = fix(sqrt(obs))
*
dim freqs(1+r)
ctor 2 ords+1 1
#     1
#  logpdgm
set freqs(1) 1 ords = -2*log(2*%pi*t/obs)
do i=1,r
   set freqs(1+i) 1 ords = (2.*%pi*t/obs)**(2*i)
   labels freqs(1+i)
   # "l**"+(2*i)
end do i
labels freqs(1)
# "LOG_LAMBDA**-2"
*
* This uses the finite sample estimate for the standard error, described on the
* top of page 686.
*
linreg(noprint) logpdgm 1 ords
# constant freqs
compute %%d=%beta(2),%%dse=%pi*sqrt(%xx(2,2)/6.0)
if .not.print
   return
display "Andrews-Guggenberger Regression, Series" %label(series)
display "m  = " ords @+5 "r = " r
display "estimated d = "     @30   #.##### %%d
display "standard error = "  @30   #.##### %%dse
end


