function [d, nobs, tasy, sigasy, tols, sigols] = gph (series, incl, excl)
%GPH Geweke & Porter-Hudak (1983) estimation of the fractional differencing parameter d
%
% [D, NOBS, TASY, SIGASY, TOLS, SIGOLS] = GPH (SERIES, INCL, EXCL) returns
%    - GPH frequency domain estimator D of (differenced) time series in vector SERIES,
%    - the number of observations NOBS used in the frequency domain regression,
%    - asymptotic t-ratio TASY and associated level of significance SIGASY (when the known
%      theoretical variance of the residuals (pi¾/6) is imposed on the calculation), and
%    - OLS t-ratio TOLS and corresponding significance level SIGOLS (when the standard
%      errors are estimated from standard OLS regression output).
%
%    SIGASY and SIGOLS result from evaluating the null hypothesis H0: d = 0 against the
%    two-sided alternative H1: d ~= 0. If TCDF.M of the MATLAB Statistics Toolbox is not
%    on the path, SIGASY and SIGOLS assume NaN.
%
%    Parameters INCL and EXCL are used to select the relevant range of frequencies:
%    - the GPH regression excludes the first N^EXCL-1 frequencies
%      (default is EXCL = 0, so all lower frequencies are actually included), and
%    - it includes the remaining lower frequencies up to the N^INCL th frequency
%      (default is INCL = 0.5, so the first SQRT(N) frequencies are included).
%    Either INCL or EXCL, but not both, can be vector input, in which case the estimation
%    results returned are all vectors corresponding to the varying values of INCL/EXCL.
%
%    The objective of the GPH estimator is to capture any fractal structure in the LOWER
%    frequencies, i.e. a correlation structure that is neither I(0) nor I(1) but I(d)
%    where d is the fractional differencing parameter to be estimated here. These LOWER
%    frequencies are covered EXCLUSIVELY (i.e. no specification of a short-memory process
%    is required), because only a fraction of the first frequencies is used - hence the
%    INCL restriction on the number of observations in the regression. Including too many
%    frequencies would blur the analysis through higher-frequency cycles, while including
%    too few frequencies would make the OLS regression less reliable. 0.5 is the standard
%    value suggested by Kunsch (1986) and usually used in the literature. As argued by
%    Kunsch (1986), also frequencies around the origin should be excluded to get a
%    consistent estimator of D; yet, this is not always done in practice. Note also that
%    if a seasonally-adjusted time series is used, another adjustment would be required to
%    eliminate seasonal frequencies before running the GPH regression (see Ooms & Hassler,
%    1997). See also the bottom of the script for a literature list.
%
% The author assumes no responsibility for errors or damage resulting from usage. All
% rights reserved. Usage of the programme in applications and alterations of the code
% should be referenced. This script may be redistributed if nothing has been added or
% removed and nothing is charged. Positive or negative feedback would be appreciated.

%                     Copyright (c) 10 March 1998 by Ludwig Kanzler
%                     Department of Economics, University of Oxford
%                     Postal: Christ Church,  Oxford OX1 1DP,  U.K.
%                     E-mail: ludwig.kanzler@economics.oxford.ac.uk
%                     Homepage:      http://users.ox.ac.uk/~econlrk
%                     $ Revision: 1.31 $$ Date: 15 September 1998 $

% Check input arguments and assign default values, if necessary:
if nargin < 1
   error('The GPH regression requires a time series input.')
elseif prod(size(series)) == 1
   error('The GPH regression cannot be run on a scalar.')
end

if nargin < 2
   incl = 0.5;
elseif isempty(incl)
   incl = 0.5;
end

if nargin < 3
   excl = 0;
elseif isempty(excl)
   excl = 0;
end

if incl >= 1
 error('The regression can only be run on fewer frequencies than there are observations.')
end

if excl < 0
   error('Input argument EXCL cannot be negative.')
end

lincl = length(incl);
lexcl = length(excl);

if lincl > 1 & lexcl > 1
   error('Only one of input arguments INCL and EXCL can be a vector, not both.')
elseif sum(excl >= incl)
   error('There are more frequencies excluded than there are included.')
end

incl(1:lincl*lexcl) = incl;   % If one of EXCL and INCL is a vector, this transforms the
excl(1:lincl*lexcl) = excl;   % other into a vector of equal values for each element.


% Take the Discrete Fourier Transform, thus obtaining the complex vector DFT (see also
% bottom of script):
dft    = fft(series(:));
n      = length(dft);
dft(1) = [];

% The following FOR loop runs the GPH regression for each parameter in EXCL or INCL:
for i = 1 : lincl*lexcl
   first(i)= ceil (n^excl(i));       % first frequency (lowest) of the regression series
   last(i) = floor(n^incl(i));       % last frequency (highest) of the regression series
   nobs(i) = last(i) - first(i) + 1; % number  of  observations of the regression series

   % The periodogram (a plot of power versus frequency) is defined by...

   % ... the Fourier frequencies (harmonic ordinates) of the sample, defined as (see GPH):
   freq    = 2*pi*(first(i):last(i)) / n;

   % ... and the corresponding "power" observations of the periodogram, i.e. the squares
   % of the magnitude of DFT:
   power   = abs(dft(first(i):last(i))).^2;

   % The spectral regression on a transformation of the periodogram (see GPH) produces the
   % the vector of slope coefficients BETA:
   y       = log(power);
   X       = [ones(nobs(i),1), log(4*sin(freq'/2).^2)];
   beta    = X \ y;

   % A consistent estimate of fractional differencing parameter D is provided by the
   % negative slope coefficient on the frequency variable:
   d(i)    = - beta(2);

   % The known theoretical variance of the residuals pi¾/6 is used to compute the ASYMPTO-
   % TIC standard errors; the OLS standard errors are computed from the actual residuals:
   seasy   = sqrt(diag(inv(X'*X)) * pi^2/6 );
   tasy(i) = d(i)/seasy(2);

   if nargout > 4
      resid   = y - X*beta;
      seols   = sqrt(diag(inv(X'*X)) * resid'*resid/(nobs(i)-2));
      tols(i) = d(i)/seols(2);
   end
end

% If TCDF.M exists, evaluate H0 against the 2-sided alternative for all t-ratios at once:
if exist('tcdf.m','file')
   sigasy = min(1-tcdf(tasy, nobs-1), tcdf(tasy, nobs-1))*2;
   if nargout > 4
      sigols = min(1-tcdf(tols, nobs-1), tcdf(tols, nobs-1))*2;
   end
else
   sigasy(1:lincl*lexcl) = NaN;
   if nargout > 4
      sigols(1:lincl*lexcl) = NaN;
   end
end

% End of function.


% The following function is not actually used by the main function and only included here
% for the benefit of someone interested in translating this script into a programming
% language which does not have a built-in fast Fourier transform function. It may also be
% helpful to understand what the explicit Fourier transform equation computes and to check
% whether the built-in function FFT really produces the same results (it does! - to try
% simply replace FFT in the above code by SFT).

function dft = sft (series, n)
% SFT "Slow" Fourier Transform
%
%   SFT (SERIES) is the discrete Fourier transform (DFT) of vector SERIES, computed using
%   the explicit Fourier equation rather than using MATLAB's built-in function FFT.
%
%   SFT (SERIES, N) returns the n-point DFT. If the length of SERIES is less than N,
%   SERIES is padded with trailing zeros to length N. If the length of SERIES is greater
%   than N, the sequence SERIES is truncated.
%
%   Please note that unlike FFT, SFT only accepts vector input.
%
%                    Copyright (c)  1 April 1998  by Ludwig Kanzler
%                    Department of Economics,  University of Oxford
%                    Postal: Christ Church, Oxford OX1 1DP, England
%                    E-mail:  ludwig.kanzler@economics.oxford.ac.uk
%                    $ Revision: 1.11  $    $ Date: 29 April 1998 $

N      = length(series);
series = series(:);

if length(series) ~= N
   error('This Fourier transformation can only be performed on vector input.')
end

if nargin == 1
   n = N;
end

series(N+1:n) = 0;
dft   (1:n,:) = 0;

% Computation for each element K in DFT (i.e. k from 1 to n):
for k = 1 : n
   dft(k) = exp(-i*2*pi*(k-1)*(0:n-1)/n) * series(1:n);
end

% NOTE: The matrix multiplication used above avoids another loop for summation; it is
% equivalent to the following routine, which is easier to understand, but much slower:
%
% for k = 1 : n
%    dftk = 0;
%    for l = 1 : n
%       dftk = dftk + series(l)*exp(-i*2*pi*(k-1)*(l-1)/n);
%    end
%    dft(k) = dftk;
% end

% End of sub-function.


% RELATED LITERATURE:
%
% Baillie, Richard (1996), "Long Memory Processes and Fractional Integration in
%    Econometrics", Journal of Econometrics, vol. 73, no. 1 (July), pp. 5-59
%
% Diebold, Francis & Marc Nerlove (1990), "Unit Roots in Economic Time Series: A
%    Selective Survey", in George Rhodes Jr. & Thomas Fomby, eds., "Advances in
%    Econometrics: A Research Annual", vol. 8 ("Co-integration, Spurious Regressions, and
%    Unit Roots"), JAI Press, Greenwich, Connecticut, pp. 3-69
%
% Geweke, John & Susan Porter-Hudak (1983), "The Estimation and Application of Long Memory
%    Time Series Models", Journal of Time Series Analysis, vol. 4, no. 4, pp. 221-238
%
% Granger, Clive & Roselyne Joyeux (1980), "An Introduction to Long-Memory Time Series
%    Models and Fractional Differencing", Journal of Time Series Analysis, vol. 1, no. 1,
%    pp. 15-29
%
% Harvey, Andrew (1993), "Time Series Models", 2nd edition, Harvester Wheatsheaf, Hemel
%    Hempstead, England, Chapter 6 ("The Frequency Domain"), pp. 166-232
%
% Hosking, J.R.M. (1981), "Fractional Differencing", Biometrika, vol. 68, no. 1, pp.
%    165-176
%
% Kunsch, H. (1986), "Discrimination between Monotonic Trends and Long-Range Dependence",
%    Journal of Applied Probability, vol. 23, pp. 1025-1030
%
% Mandelbrot, BenoÓt (1977), "Fractals: Form, Chance, and Dimension", Free Press,
%    New York
%
% Ooms, Marius & Uwe Hassler (1997), "On the Effects of Seasonal Adjustment on the
%    Periodogram Regression", Economics Letters, vol. 56, no. 2, pp. 135-141


% This MATLAB script draws on a procedure used (and described) in the following papers:
%
% Barkoulas, John, Christopher Baum & Mustafa Caglayan (1998), "Persistent Dependence in
%    Foreign Exchange Rates?", Journal of International Money and Finance, forthcoming,
%    revised version of Boston College Working Papers in Economics, no. 377
%
% Barkoulas, John, Christopher Baum & Mustafa Caglayan (1998), "Long Memory or Level
%    Shifts: Can Either Explain Nonstationary Real Exchange Rates Under the Current
%    Float?", Boston College Working Papers in Economics, no. 380
%
% Barkoulas, John, Christopher Baum & Nickolaos Travlos (1998), "Long Memory in the Greek
%    Stock Market", Journal of Applied Financial Economics, forthcoming, revised version
%    of Boston College Working Papers in Economics, no. 356 (August 1997)
%
% Baum, Christopher & John Barkoulas (1996), "Testing for Fractal Structure in Stock
%    Returns", Boston College Working Papers in Economics, no. 314 (April)
%
% Berg, Lennart & Johan Lyhagen (1996), "Short and Long Run Dependence in Swedish Stock
%    Returns", Uppsala University, Department of Economics, working paper, no. 1996:19
%    (September)


% ACKNOWLEDGEMENT:
%
% The author is grateful to Christopher Baum for a discussion of the GPH procedure.


% End of file.