% GE march 01 proc to compute the Q(g,k) statistic from the
%   initial condition paper.
%  this method does the usual OLS detrending as in the minitial paper
% inputs 
%  
% z    : Tx1 vector of data to be tested for a unit root
% g    : local point for POI test. Usual choice is g=10.
% k    : prior on initial condition under Ha in terms of unconditional 
%        variance of z, larger means larger variance 
% fdet : flag for deterministics, =1 for constant, =2 for time trend
% s0   : scaling parameter - spectral density at freq zero 
%        (this needs to be supplied to the procedure, any consistent
%        estimator is fine for the asymptotic theory).
%
% output
%
% q    : q statistic test for unit root
% format [q] = qstat0(z,g,k,fdet,s0);

function [q] = qstat0(z,g,k,fdet,s0);


% setup weights for M statistics 
if fdet==1;
 q0=-g;
 q1=-g*(1+g)*(k-2)/(2+g*k);
 q2=(2*g-g*k+g*g*k)/(2+g*k);
 q3=2*g*(k-2)/(2+g*k);
 q4=g*g;
 else;
 q0=-g;
 q1=(8*g*g+8*g*g*g-3*g*g*g*k-g*g*g*g*(k-2))/(24+24*g+8*g*g+k*g*g*g);
 q2=(8*g*g+g*g*g*(8-3*k)+g*g+g*g)/(24+24*g+8*g*g+k*g*g*g);
 q3=(8*g*g+2*g*g*g*(4-3*k))/(24+24*g+8*g*g+k*g*g*g);
 q4=g*g;
end;


% construct M statistics 
n=length(z);

if fdet==1;
    y=z-mean(z);
elseif fdet==2;
  x=[ones(n,1) (1:n)'];
  y=z-x*inv(x'*x)*x'*z;
end;

m1=(y(n,1)/sqrt(s0*n));
m0=(y(1,1)/sqrt(s0*n));
im2=(y'*y)/(n*n*s0);

% construct test statistic 

q=q0+q1*m0*m0+q2*m1*m1+q3*m1*m0+q4*im2;