*
* @RGSE( options ) series start end
* Computes Robinson's Gaussian Semiparametric Estimator (GSE) for the
* fractional differencing parameter.
*
* Parameters:
*   series     series to analyze
*   start end  range of <<series>> to use. By default, the defined range of <<series>>
*
* Options:
*   POWER=power of nobs to use as frequencies [.8]
*   [PRINT]/NOPRINT
*   TITLE=title of report ["Robinson GSE, Series xxx"]
*
* Variables Defined:
*   %%D    estimated value of D
*   %%DSE  estimated standard error of D
*
* Critical Info:
*   Most often you will need to difference the data before estimation.
*   d is constrained to (-1.5 < d < .5).
*
* Reference:
*   Robinson(1992), "Semiparametric analysis of long-memory time series",
*    Annals of Statistics, vol 22, 515-539.
*
* Revision Schedule:
*   12/1996 Written by Rob Schoen (translated from S-Plus).
*           email: rjschoen@earthlink.net or rjs2@stern.nyu.edu
*   12/2006 Updated to use SSTATS and FIND by Tom Doan, Estima
*   06/2012 Change to use %%D and %%DSE for compatiblity with other
*           procedures. Use REPORT. Add PRINT and TITLE options.
*           Use local PARMSET.
*   04/2018 Don't create REPORT if NOPRINT
*
procedure RGSE series start end
type series series
type integer start end
*
option real    power  0.8
option switch  print  1
option string  title
*
local integer  n m
local series   lam pgram
local real     m1 m2
local integer  startl endl
local real     rgseval
local parmset  dparm
local report   rreport
local string   ltitle
*
declare real 	%%d %%dse
*
inquire(series=series) startl<<start endl<<end
compute n = endl-startl+1
compute m = fix(n^power)
*
freq 1 n
rtoc startl endl 1
# series
# 1
fft 1
cmult(scale=1.0/(2*%pi*n)) 1 1
set pgram 1 m = %real(%z(t+1,1))
set lam   1 m = 2 * sin(%pi*t/n)
sstats(mean) 1 m log(lam)>>m2
nonlin(parmset=dparm) %%d
compute %%d=.5
find(noprint,parmset=dparm) min rgseval
   sstats(mean) 1 m pgram*lam^(2.0*%%d)>>m1
   compute rgseval=log(m1)-2.0*%%d*m2
end find
compute %%dse = sqrt((1/(4.0*m)))
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Robinson GSE, Series "+%label(series)
*
compute %nobs=n
if .not.print
   return

report(use=rreport,action=define,title=ltitle)
report(use=rreport,atrow=1,atcol=1,span) ltitle
report(use=greport,atrow=2,atcol=1,span) "Data From "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=rreport,atrow=4,atcol=1) "Observations" %nobs
report(use=rreport,atrow=5,atcol=1) "Frequencies" m
report(use=rreport,atrow=6,atcol=1) "d-hat" %%d
report(use=rreport,atrow=7,atcol=1) "std error" %%dse
report(use=rreport,action=format,atrow=6,picture="#.####")
report(use=rreport,action=show)
end
