*
* @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
*
* 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.
*
* Reference:
*   Geweke and Porter-Hudak(1983), "The Estimation and Application of Long Memory
*   Time Series Models", J. of Time Series Analysis, pp. 221-238.
*
procedure gph  series  start  end
type series  series
type integer start  end
*
options real power 0.50
*
local integer startl endl ords obs
local vector[series] freqs
*
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"
display "Geweke-Porter-Hudak Regression, Series" %label(series)
linreg(noprint) freqs(1)
# constant freqs(2)
display "Power = " POWER  @+5 "Regression Ordinates = " %nobs
display "Estimated d = "     @30   #.##### -%beta(2)
display "Asymp Standard Error = "  @30   #.##### %pi*sqrt(%xx(2,2)/6.0)
display "OLS Standard Error = "  @30   #.##### sqrt(%seesq*%xx(2,2))
end




