function [d, seasy, seprd] = EstimateGPH (mSeries,iTrunc)
% PURPOSE: The log-periodogram regression method for estimating d.
% -----------------------------------------------------
% USAGE:     [d, seasy, seprd] = EstimateGPH (mSeries,iTrunc)
%  where:    mSeries = a series matrix
%            iTrunc  = The number of remaining lower frequencies up to iTrunc, the default value is SQRT(N).
% -----------------------------------------------------
% RETURNS:   d= GPH frequency domain estimator D of (differenced) time series matrix
%            seasy,seprd =  result from evaluating the null hypothesis 
% -----------------------------------------------------             
% NOTE:      the function will be used for the initialization 
% -----------------------------------------------------


% Written and changed by
% Ying Ma
% yingma@email.com


if nargin < 1
   error('Error:The GPH regression requires a time mSeries input.')
elseif prod(size(mSeries)) == 1
   error('Error:The GPH regression cannot be run on a scalar.')
end

[ct n]    = size(mSeries);

if ct == 1 & n > 1
    mSeries = mSeries(:);
    [ct n]    = size(mSeries);
end

if nargin < 2
   iTrunc = sqrt(ct);
end

iTrunc = floor(iTrunc);
   
cp     = floor ( ct / 2);
  
if ct < 1024
    ymean   = mean(mSeries); 
    mSeries  = mSeries - ones(ct,1)*ymean;
    dft     = fft(mSeries);
    dft(1,:)  = [];
    dft     = abs(dft(1:cp,:)).^2;
    cper    = ct;
else
    dft   = fft(mSeries);
    dft(1,:)  = [];
    cper  =  size(dft,1);
    dft    = dft(1:cp,:);
    dft    = (abs(dft).^2) / (2*pi*ct) ;
end

if  iTrunc > 0 & iTrunc < cp
    dft = dft(1:iTrunc,:);
    cp  = iTrunc;
end

Freq    =  2*pi*[1:cp] / cper;

y       = log(dft);

x       = [ones(cp,1) -log(2*(1-cos(Freq')))];

beta    = olsv(y, x);

d       =  beta (2,:);

seasy   = pi*sqrt(diag(inv(x'*x)) /6 );
seasy   =   seasy(2:end)';

seprd   = 2*(1-cdf('Normal',d./seasy,0,1));
