%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This code implements the Perron-Yabu (2009) procedure  "Estimating Deterministic Trends with an Integrated or
%Stationary Noise Component,"  Journal of Econometrics  151 (2009), 56-69.
function [tstat1,stdb1,b1,ams,tstat2,stdb2,b2,aus,ahat,am,au,p]=tcode(dep,kmax)

[p] = mic1(dep,1,0,kmax,0);
[ahat,am,au] = RF(dep,p);
[tstat1,stdb1,b1,ams] = QFGLS(dep,p,am);
[tstat2,stdb2,b2,aus] = QFGLS(dep,p,au);
%********************************************************************%
%lagn(x,n) => x{t-n}*/
function y = lagn(x,n)
	if n>0;
		y = [zeros(n,size(x,2)) ; trimr(x,0,n)];
    else
		y = [trimr(x,abs(n),0) ; zeros(abs(n),size(x,2))];
    end


%********************************************************************%
%diff(x,n) with n>0 => x{t}-x{t-n}*/  
function y = diff(x,n)
	if n==0
		y = x;
    else
		y = [zeros(n,size(x,2)) ; (trimr(x,n,0)-trimr(lagn(x,n),n,0))];
    end

%********************************************************************%
%median-unbiased and upper-biased estimator
function [rhows,rhomd1,rhomd2]= RF(x,p)
    T=size(x,1);
    dep=x; reg=[ones(T,1),seqa(1,1,T)];
    %weighted symmetric least squares for AR(p)%
    wt=zeros(T,1);
    i=p+1;
    while i<=T;
        if i<=T-p+1;
            wt(i)=(i-p)/(T-2*p+2);
        elseif i>T-p+1;
	        wt(i)=1;
        end
        i=i+1;
    end
    w=diagrv(zeros(2*(T-p),2*(T-p)),[wt(p+1:T);(1-rev(wt(2:T-p+1)))]);         
    u=(eye(T)-reg*inv(reg'*reg)*reg')*dep;				
    depa=[u(p+1:T);rev(u(1:T-p))];
    rega=[u(p:T-1);rev(u(2:T-p+1))];
    i=2;
    while i<=p;
        rega=[rega,[u(p-i+1:T-i);rev(u(1+i:T-p+i))]]; 
        i=i+1;
   end
    bhat=inv(rega'*w*rega)*rega'*w*depa; 
    ss=(depa'*w*depa-bhat'*rega'*w*depa)/(T-p-2);
    vbhat=ss*inv(rega'*w*rega);
    vec1=ones(p,1);
    rhows = vec1'*bhat;
    vrhows = vec1'*vbhat*vec1;
    tau1=(rhows-1)/sqrt(vrhows);

    %Median unbiased estimator%
    t05=-1.96;
    k1=5; k2=sqrt(3*T);Ip=fix((p+1)/2);
    cont=(3-t05*t05*(1+Ip/T))/(t05*(t05+k1)*(1+Ip/T));
    if (tau1>t05); 
	    ctau1=-tau1; 
    elseif ((tau1<= t05)*(tau1>-k1))==1; 
	    ctau1=tau1/T - 3/(tau1+cont*(tau1+k1));
    elseif ((tau1<= -k1)*(tau1>-k2))==1; 
	    ctau1=tau1/T-3/tau1;
    elseif (tau1<= -k2); 
	    ctau1=0;
    end
    rhomd1 = rhows + ctau1*sqrt(vrhows(1,1));
    rhomd1 = min(rhomd1|1);

   %Upper biased estimator%
    t05=-2.85;
    cont=(3-t05*t05*(1+Ip/T))/(t05*(t05+k1)*(1+Ip/T));
    if (tau1>t05); 
	    ctau1=-tau1; 
    elseif ((tau1<= t05)*(tau1>-k1))==1; 
	    ctau1=tau1/T - 3/(tau1+cont*(tau1+k1));
    elseif ((tau1<= -k1)*(tau1>-k2))==1; 
	    ctau1=tau1/T-3/tau1;
    elseif (tau1<= -k2); 
	    ctau1=0;
    end
    rhomd2 = rhows + ctau1*sqrt(vrhows(1,1));
    rhomd2 = min(rhomd2|1);

%********************************************************************%
% Quasi Feasible GLS*/
function [t_QFGLS,stdb,b,amus]= QFGLS(x,p,ahat)
	T=size(x,1); 
    dep=x; dx=diff(dep,1);
    reg=[ones(T,1),seqa(1,1,T)];
	CR=(T^0.5)*abs(ahat-1);	
	amus=ahat*(CR>1)+1*(CR<=1);%supper efficient estimator%
	depg=([dep(1),(dep(2:T)-amus*dep(1:T-1))'])';
	regg=([ones(1,2)',(reg(2:T,:)-amus*reg(1:T-1,:))'])';	
	[b,v,~]=olsqr2(depg,regg);		

	
	
	if (p==1);%Estimation of long-run variance%
		h0=v'*v/(size(v,1));
	elseif (p>1)
		j=1;
		if amus==1 %judged as a unit root%
			while j<=p-1
				if (j==1)
					regv=lagn(v,j);
                else
					regv=[regv,lagn(v,j)];
				end
				j=j+1;
            end		
			depv=trimr(v,p-1,0);
			regv=trimr(regv,p-1,0);		
			[phi,e,~]=olsqr2(depv,regv);				
			h0=(e'*e/(T-(p-1)))/((1-sum(phi))^2);			
		elseif amus<1 %judged as a stationary process
	            while j<=p-1;
				if j==1
					regx=[reg,lagn(dx,j)];
                else
					regx=[regx,lagn(dx,j)];
				end
				j=j+1;
                end			
			depx=trimr(depg,p-1,0);	
			regx=trimr(regx,p-1,0);		
			[~,e,~]=olsqr2(depx,regx);					
			h0=e'*e/(T-(p-1));			
		end
	end

	VCV=(h0)*inv(regg'*regg);
    stdb=sqrt(diag(VCV));
	t_QFGLS=(b(2)-0)/stdb(2);%t-statistic%


%********************************************************************%
%This is the code written by NG and perron. See the following papers:
%Ng, S. and P. Perron (1997): "Lag Length Selection and Constructing unit
%root tests with good size and power," manuscript, Boston College and
%Boston University. BC WP 369.
%See also Perron and Ng (1996), Review of Economic Studies, Useful
%Modifications to Unit root tests with dependent errors and their local
%asymptotic properties, vol 63, p.435-464.*/
%********************************************************************%
function krule=mic1(y,p,penalty,kmax,kmin)
%------------------------------------------------------------
%y: a univariate time series of dimension nt.
%p: the order of the polynomial trend; either
%   a) p = 0 with only a constant in z, or
%   b) p = 1, with a constant and a linear trend in z.
%penalty: 0 for mic and 1 for bic.
%---------------------------------------------------------------------*/

nt=size(y,1);
if p == 0;
z=ones(nt,1);
cbar=-7.0;
end
if p == 1;
z=[ones(nt,1),seqa(1,1,nt)];
cbar=-13.5;
end
cbar=cbar;

[yt,~]=glsd(y,z,cbar); % transforming the data*/
[~,r,~]=olsqr2(yt(2:nt,1),yt(1:nt-1,1));  % estimate the alpha-hat*/
s2u=r'*r/(nt-1);
sumyt2=sum(yt(1:nt-1,1).^2)/(nt-1).^2;
krule=s2ar(yt,penalty,kmax,kmin); % estimate s2ar */


function kopt=s2ar(yts,penalty,kmax,kmin)
nt=size(yts,1);
min=9999999999;
tau=zeros(kmax+1,1);
s2e=999*ones(kmax+1,1);
dyts=diff(yts,1);
reg=lagn(yts,1);
i=1;
while i <= kmax
reg=[reg,lagn(dyts,i)];
i=i+1;
end
dyts0=dyts;
reg0=reg;
%loop over k*/
dyts0=trimr(dyts,kmax+1,0);
reg0=trimr(reg,kmax+1,0);
sumy=sum(reg0(:,1).*reg0(:,1));
nef=nt-kmax-1;
k=kmin;%kmin=0 means that the process is simply AR(1)%
while k <= kmax
%[b,e,fit]=olsqr2(dyts0,reg0(:,1:k+1));
b=reg0(:,1:k+1)\dyts0;
e=dyts0-reg0(:,1:k+1)*b;
s2e(k+1)=e'*e/nef;
tau(k+1)=(b(1)*b(1))*sumy/s2e(k+1);
k=k+1;
end
kk=seqa(0,1,kmax+1);
if penalty == 0
mic=log(s2e)+2.0*(kk+tau)./nef;
else
mic=log(s2e)+log(nef)*(kk)./nef;
end
kopt=minindc(mic);


function [yt,ssr]=glsd(y,z,cbar)
nt=size(y,1);
abar=1+cbar/nt;
ya=zeros(nt,1);
za=zeros(nt,size(z,2));
ya(1:1,1)=y(1:1,1);
za(1:1,:)=z(1:1,:);
ya(2:nt,1)=y(2:nt,1)-abar*y(1:nt-1,1);
za(2:nt,:)=z(2:nt,:)-abar*z(1:nt-1,:);
%constructing the gls detrended series*/
bhat=inv(za'*za)*za'*ya;
yt=y-z*bhat;
ssr=(ya-za*bhat)'*(ya-za*bhat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ahat,r,fit]=olsqr2(y,x)
ahat=(x'*x)\(x'*y);
r=y-x*ahat;
fit=x*ahat;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         AUXILIARY GAUSS TO MATLAB FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z=trimr(x,n1,n2)
% PURPOSE: return a matrix (or vector) x stripped of the specified rows.
% -----------------------------------------------------
% USAGE: z = trimr(x,n1,n2)
% where: x = input matrix (or vector) (n x k)
%       n1 = first n1 rows to strip
%       n2 = last  n2 rows to strip
% NOTE: modeled after Gauss trimr function
% -----------------------------------------------------
% RETURNS: z = x(n1+1:n-n2,:)
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl%jpl.econ.utoledo.edu

  [n, ~] = size(x);
  if (n1+n2) >= n; 
     error('Attempting to trim too much in trimr');
  end;
  h1 = n1+1;   
  h2 = n-n2;
  z = x(h1:h2,:);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function seq=seqa(a,b,c)
% PURPOSE: produce a sequence of values
% -----------------------------------------------------
% USAGE: y = seqa(a,b,c)
%  where    a = initial value in sequence 
%           b = increment
%           c = number of values in the sequence  
% -----------------------------------------------------
% RETURNS: a sequence, (a:b:(a+b*(c-1)))' in MATLAB notation
% ----------------------------------------------------- 
% NOTE: a Gauss compatability function
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu
       
% seqa Gauss eqivalent of seqa(a,b,c)
seq=(a:b:(a+b*(c-1)))';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = diagrv(x,v)
% PURPOSE: replaces main diagonal of a square matrix
% -----------------------------------------
% USAGE: y - diagrv(x,v)
% where: x = input matrix
%        v = vector to replace main diagonal
% -----------------------------------------
% RETURNS: y = matrix x with v placed on main diagonal
% -----------------------------------------
% NOTE: a Gauss compatability function
% -----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu
  
[r,c] = size(x);
if r ~= c;
  error('x matrix not square')
end;
if length(v) ~= r;
  error('v is not conformable with x')
end;
y = x - diag(diag(x)) + diag(v);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c=minindc(A)
    [~,n]=size(A);
    c=zeros(n,1);
    for i=1:n
        [~,index]=min(A(:,i));
        c(i)=index;
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function re=rev(x)
% REV 
% REV(x) reverses the elements of x
% It is equivalent to the Matlab functin flipud
d=size(x); li=d(1,1);
ind=li:-1:1;
re=x(ind,:);






