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.