*
* GPH_SEAS.SRC
*     Estimates fractional difference power for series using the
*     frequency domain regression techniques of Geweke and Porter-Hudak,
*     with adjustments proposed by Ooms and Hassler (1996) to prevent
*     inclusion of ordinates corresponding to seasonal frequencies (which,
*     for seasonally adjusted data,  contain log-periodogram values of
*     minus infinity; large negative values will be outliers in the
*     log-periodogram regression).
*
*     References: 
*     Geweke and Porter-Hudak, "The Estimation and Application
*     of Long Memory Time Series Models", J. of Time Series Analysis,
*     1983, (221-238).
*
*     Marius Ooms and Uwe Hassler, "On the Effects of Seasonal Adjustment
*     on the Periodogram Regression," Economics Letters, 1997, 56:2, 135-141.
*
*  Syntax:
*     @GPH_SEAS  series  start  end
*     series usually is the first difference of the series actually
*     being studied.
*
*  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
*
*     SADJ = [TRUE] to adjust sample size to even number of years by
*            zero padding and omit seasonal frequencies from log-periodogram
*            regression
*
*  Revision Schedule (originally Estima's GPH.SRC):
*    12/92 SIN function should be evaluated at (2(pi)j/T)/2,j=0,...,T-1.
*          Added calculation of both asymptotic and OLS standard errors
*    01/98 added explicit limits to LINREG (baum@bc.edu)
*    01/98 added print of t-stats (baum@bc.edu)
*    02/98 modified GPH.SRC to incorporate handling of seasonal frequencies
*          (baum@bc.edu)
*
PROCEDURE GPH_SEAS  SERIES  START  END
TYPE SERIES  SERIES
TYPE INTEGER START  END
*
OPTIONS REAL POWER 0.50
options switch  sadj 1
*
LOCAL INTEGER STARTL ENDL BASE ORDS seasfq nunit
LOCAL VECTOR[SERIES] FREQS
*
INQUIRE(SERIES=SERIES) STARTL>>START ENDL>>END
inquire(seasonal) seasfq
if sadj .and. seasfq ==1 {
dis 'Warning: seasonal = 1'
}
COMPUTE OBS = ENDL-STARTL+1
*
if sadj .and. seasfq>1 {
  comp nunit=fix(obs/seasfq)
  comp obs=(nunit+1)*seasfq
}
*
FREQ 3 OBS
RTOC STARTL ENDL 1
# SERIES
#   1
*
FFT 1
CMULT 1 1
CLN 1
*
cset 3 = 1.0
cset 3 1 1 = 0.0
CSET 2 = LOG(4 * SIN(%PI*(T-1)/OBS)**2)
if sadj {
cmask 3 / obs/seasfq 1 
} 
* cprint(lc=periods) / 1 3
*
COMPUTE ORDS = FIX(OBS**POWER)
*
DIM FREQS(3)
CTOR 1 ORDS+1 1
#     1      2         3
#  FREQS(1)  FREQS(2) freqs(3)
LABELS FREQS(1) FREQS(2) freqs(3)
# 'LOGPER' 'LOGSIN2' 'ssmpl'
DISPLAY 'Geweke-Porter-Hudak Regression, Series' %LABEL(SERIES)
* print / freqs(1) freqs(3)
LINREG(NOPRINT,smpl=freqs(3)) FREQS(1) 1 ords+1
# CONSTANT FREQS(2)
clear freqs(1) freqs(2) freqs(3)
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) $
 @40 '(' ##.### -%beta(2)/(%pi*sqrt(%xx(2,2)/6.0)) ')'
DISPLAY 'OLS Standard Error = '  @30   #.##### sqrt(%seesq*%xx(2,2)) $
 @40 '(' ##.### -%beta(2)/sqrt(%seesq*%xx(2,2)) ')'
END
