*
* @AGFRACTD(options) series  start  end
* Estimates fractional difference power for series using the bias-reduced
* technique from Andrews and Guggenberger.
*
* Options:
*   POWER=Lowest T^power frequencies are used in running the
*         regression [not used]
*   M    =Number of frequencies to use in running the regression [Andrews plug-in value]
*
*   These are alternative for choosing the number of frequencies to use

*   R  = Number of even powers of the frequencies to use [R=0 gives Geweke-Porter-Hudak]
*
*   [PRINT]/NOPRINT
*   TITLE=title of report ["Andrews-Guggenberger Regression, Series xxx"]
*
* Variables Defined:
*   %%D     estimated fractional difference
*   %%DSE   standard error of estimate for %%D
*
* 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
*
option   integer  m
option   real     power
option   integer 	r     	0
option   switch  	print 	1
option   string  	title
*
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
local    string   ltitle
local    report   areport
*
declare 	real     %%d %%dse
*
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.0*i+2*j+1)*(2.0*i+1)*(2.0*j+1))
ewise mu(i)=2.0*i/(2.0*i+1)^2
ewise xi(i,j)=2.0*j*(3.0+2.0*i)/((2.0*i+2.0*j+3.0)*(2.0*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.0*%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
if %defined(power)
   compute ords = fix(obs^power)
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)
compute %nobs=obs
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Andrews-Guggenberger Regression, Series "+%label(series)

report(use=areport,action=define,title=ltitle)
report(use=areport,atrow=1,atcol=1,span) ltitle
report(use=areport,atrow=2,atcol=1,span) "Data From "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=areport,atrow=4,atcol=1) "Even Powers"           r
report(use=areport,atrow=5,atcol=1) "Regression Ordinates" 	ords
report(use=areport,atrow=6,atcol=1) "Estimated d"           %%d
report(use=areport,atrow=7,atcol=1) "Standard Error"        %%dse
report(use=areport,atrow=6,action=format,picture="#.####")
if print
   report(use=areport,action=show)
end
