function [out]=gph(data, m, plotit)
% calculate GPH estimator as in Perron and Qu (2008)
% Use FFT to calculate
% [d]=gph(data, m)
% m is the maximum number of frequencies for GPH
if nargin < 3, plotit=1; end
if nargin < 2, m = length(data)/2; end
% Check the data; data should be row vector
[d1 d2] =size(data);
if d1 > 1, data = data'; end;
n = length(data);
dft = fft(data); dft(1) = [];
omega = 2*pi*(1:floor(n/2))./n;
power = (abs(dft).^2)/(2*pi*n);
y = log(power);
k = min(m, length(omega));
if k > 3000;
lk = linspace(10, k, 3000);
lk = floor(lk);
else lk = 10:1:k;
end;
d=[];
%d1=[];
for j=1:1:length(lk)
X =[ones(lk(j),1) log(2*sin(omega(1:lk(j))./2))'];
beta = X\y(1:lk(j))';
d(j) = -beta(2)/2;
%X1 =[ones(j,1) log(omega(1:j))'];
%beta1 = X1\y(1:j)';
%d1(j-9) = -beta1(2)/2;
end;
ac = acf(data, length(data), 0);
ylabelv = [min(ac)-.1 max(ac)+.1];
out.d = d;
out.ac = ac;
out.k = lk;
mopt = floor(sqrt(n));
if mopt >= 10;
out.opt = d(lk == mopt);
out.optci = [out.opt-1.96*pi/sqrt(24*mopt) out.opt+1.96*pi/sqrt(24*mopt)];
end;
if mopt < 10;
out.opt = d(lk == floor(n^(.67)));
out.optci = [out.opt-1.96*pi/sqrt(24*floor(n^(.67))) out.opt+1.96*pi/sqrt(24*floor(n^(.67)))];
end;
if plotit==1
clf();
subplot(2,2,1),
plot(data);
xlim([1 length(data)]);
title('Data plot','fontsize', 12);
% set(gca,'PlotBoxAspectRatio',[2 ,1 ,1]);
subplot(2,2,2)
plot((1:1:length(ac))./length(ac), ac);
hold on;
plot([0 1], [0 0], 'b');
xlabel('Lag fraction \kappa = (h/T)');
ylim(ylabelv);
ylabel('$\widehat{\rho}$', 'interpreter', 'latex', 'fontsize', 12);
title(['Sample ACF T=', num2str(length(data))]);
hold on;
plot([.293 .293], ylabelv, 'r', 'Linewidth', 1);
hold on;
plot([.592 .592], ylabelv, 'r', 'Linewidth', 1);
subplot(2,2,3),
id = [n^(1/3) n^(1/3) n^(1/2) n^(1/2) n^(2/3) n^(2/3)];
id = floor(id);
id = omega(id);
periody = power(1:length(omega));
pym = min(periody); pyM = max(periody);
loglog(omega, periody);
hold on;
plot(id(1:2), [pym pyM] ,'r', 'Linewidth', 1);
hold on;
plot(id(3:4), [pym pyM] ,'r', 'Linewidth', 1);
hold on;
plot(id(5:6), [pym pyM] ,'r', 'Linewidth', 1);
title('Periodogram', 'fontsize', 12);
xlabel('$\log(\omega_j)$', 'interpreter', 'latex', 'fontsize', 12);
ylabel('$\log(I_{x,T})$', 'interpreter', 'latex', 'fontsize', 12);
axis tight; % axis square;
subplot(2,2,4),
plot(lk, d, 'Linewidth', 1.5);
% plot(10:1:k, d, 10:1:k, d1);
hold on;
plot([n^(1/3) n^(1/3)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
hold on;
plot([n^(1/2) n^(1/2)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
hold on;
plot([n^(2/3) n^(2/3)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
axis tight; % axis square;
ylim([min(d)-.05 max(d)+.05]);
xlim([0 max(lk)+1]);
xlabel('Number of frequencies');
ylabel('GPH $\widehat{d}$', 'interpreter', 'latex', 'fontsize', 12);
title('GPH estimates of d ', 'fontsize', 12);
end;
if plotit ==2
plot(lk, d, 'Linewidth', 1);
% plot(10:1:k, d, 10:1:k, d1);
hold on;
plot([n^(1/3) n^(1/3)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
hold on;
plot([n^(1/2) n^(1/2)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
hold on;
plot([n^(2/3) n^(2/3)], [min(d)-.05 max(d)+.05] ,'r', 'Linewidth', 1);
axis tight; % axis square;
% ylim([min(d)-.05 max(d)+.05]);
xlim([0 max(lk)+1]);
xlabel('Number of frequencies');
ylabel('GPH $\widehat{d}$', 'interpreter', 'latex', 'fontsize', 12);
title('GPH estimates of d ', 'fontsize', 12);
end;