*
* @GPH( options)  series  start  end
* Estimates fractional difference power for series using the frequency
* domain regression techniques of Geweke and Porter-Hudak.
*
* Parameters:
*   series      series to analyze. This is usually the first difference
*               of the series actually being studied.
*   start end   range to analyze (by default, full range of series)
*
* Options:
*   POWER=Lowest T^power frequencies are used in running the
*         regression. Default is .5 (i.e. use square root of the number
*         of observations
*   [PRINT]/NOPRINT
*   TITLE=title for output ["Geweke-Porter-Hudak Regression: Series ...]
*
* Revision Schedule (by Estima):
*   12/1992 SIN function should be evaluated at (2(pi)j/T)/2,j=0,...,T-1.
*           Added calculation of both asymptotic and OLS standard errors
*
*   11/1997 Added CLEAR instruction to clear out FREQS(1) and (2). This avoids
*           problems that could occur if you re-executed GPH in a given session
*           with a shorter number of ordinates (i.e. by using a smaller POWER),
*           which could leave old data in the FREQS series.
*
*   07/2010 Switch to use REPORT. Add PRINT and TITLE options.
*   06/2012 Define %%D and %%DSE
*
* Variables Defined:
*   %%D    estimated value of D
*   %%DSE  estimated standard error of D
*   %NOBS  number of ordinates in the regression
*
* Reference:
*   Geweke and Porter-Hudak(1983), "The Estimation and Application of Long
*    Memory Time Series Models", J. of Time Series Analysis, vol 4, no 4,
*     221-238.
*
procedure gph  series  start  end
type series  series
type integer start  end
*
options 	real 				power 	0.50
option  	switch			print		1
option 	string			title
*
local 	integer 			startl endl ords obs
local 	vector[series] freqs
local 	string			ltitle
local 	report			greport
*
declare 	real         	%%d %%dse
*
inquire(series=series) startl<<start endl<<end
*
compute obs = endl-startl+1
*
freq 2 obs
rtoc startl endl 1
# series
#   1
*
fft 1
cmult 1 1
cln 1
*
cset 2 = log(4 * sin(%pi*(t-1)/obs)^2)
*
compute ords = fix(obs^power)
*
dim freqs(2)
clear freqs(1) freqs(2)
ctor 2 ords+1 1
#     1      2
#  freqs(1)  freqs(2)
labels freqs(1) freqs(2)
# "Log_Pdgm" "LogSin2"
linreg(noprint) freqs(1)
# constant freqs(2)
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Geweke-Porter-Hudak Regression, Series "+%label(series)

compute %%d  =-%beta(2)
compute %%dse=%pi*sqrt(%xx(2,2)/6.0)
report(use=greport,action=define,title=ltitle)
report(use=greport,atrow=1,atcol=1,span) ltitle
report(use=greport,atrow=2,atcol=1,span) "Data From "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=greport,atrow=4,atcol=1) "Power" 						power
report(use=greport,atrow=5,atcol=1) "Regression Ordinates" 	%nobs
report(use=greport,atrow=6,atcol=1) "Estimated d" 				%%d
report(use=greport,atrow=7,atcol=1) "Asymp Standard Error"  %%dse
report(use=greport,atrow=8,atcol=1) "OLS Standard Error"		sqrt(%seesq*%xx(2,2))
report(use=greport,atcol=2,action=format,picture="#.#####")
if print
   report(use=greport,action=show)
end


