%  URREG.M



% Procedures to implement unit root testing with regression covariates.



% Bruce E. Hansen
% Department of Economics
% Social Science Building
% University of Wisconsin
% Madison, WI 53706-1393
% bhansen@ssc.wisc.edu
% http://www.ssc.wisc.edu/~bhansen/



% This file contains four procedures:
% UR_REG, LR_VAR, UR_CRITS, UR_ADF.




% (1) UR_REG produces the coefficient and studentized statistics
% for unit roots outlined in B. Hansen "Rethinking the univariate approach
% to unit root testing: Using covariates to increase power" 
% Econometric Theory (1995).



% The format is:



% {tstat,crits,rho2,sig2} = UR_REG(y,x,p,q,k1,k2);



% The inputs are:



% y   -  dependent variable, Tx1.
% x   -  right-hand-side variables, Txk.
% p   -  order of deterministic time trend.
%        if p=-1, then none included
%        if p=0,  then constant included
%        if p =1, then constant and time trend included.
% q   -  number of lagged dy terms
% k1  -  number of lagged x terms
% k2  -  number of lead x terms



% The outputs are:



% tstat   -  studentized coefficient test
% crits -  asymptotic critical values (1%, 5%, 10%) for tstat
% rho2    -  estimated "Rho-Squared"
% sig2    -  estimated error variance



% Notes:



% If k1=k2=0, then the variable x appears on the right-hand-side of the
% regression.  This is useful if the user wishes to directly control the
% context of the included covariates.



% The variable "x" needs to be stationary.  In most cases, it will be the
% first-difference of some levels variable.  (In the notation of the paper,
% it should be "dx").




% (2) LR_VAR  calculates a "Long-Run Variance" for a vector-valued
% series, using the suggestions of Andrews (1991).  The format is:



% {omega,bandw} = LR_VAR(e);



% where e is a Txk matrix of observations, omega is a kxk covariance
% matrix, and bandw is the bandwidth used in estimation of omega.




% Both UR_REG and LR_VAR depend on the following set of global variables.



% kernel_   -  indicator for kernel method
%              if kernel_ = 1, then Parzen kernel is used
%              if kernel_ = 2, then Quadratic Spectral kernel is used
%              if kernel_ = 3, then Bartlett (Newey-West) kernel is used
% band_     -  bandwidth used to calculate long-run variance matrix
%              if band_ = 0, then the bandwidth is calculated using
%              the suggestion of Andrews (1991)
% white_    -  pre-whiten error indicator
%              if white_ = 1, then the errors are pre-whitened with a VAR(1)
%              before kernel is applied.
% urprint_  -  print to screen supressor
%              if urprint_ = 0, then results (for UR_REG and UR_ADF) are printed
%              to screen
%              if urprint_ = 1, then no printing to screen
%             (useful for Monte Carlo runs)



% Global Defaults are set below
% Since they are globals, they may be changed in a program if used.




% (3)  UR_CRITS gives asymptotic critical values for the unit root tests.
% The format is



% crits = ur_crits(rho2,p);



% The inputs are:



% rho2  - R-squared measure, as given in UR_REG output
% p     -  order of deterministic time trend, as given in UR_REG input
%          if p=-1, then none included
%          if p=0,  then constant included
%          if p =1, then constant and time trend included.



% The outputs are:



% crits - asymptotic 1%, 5%, 10% critical values for tstat




% (4)  UR_ADF



% Format:



% {tstat,sig2} = UR_ADF(y,p,q);



% The inputs are:



% y   -  dependent variable, Tx1.
% p   -  order of deterministic time trend.
%        if p=-1, then none included
%        if p=0,  then constant included
%        if p =1, then constant and time trend included.
% q   -  number of lagged dy terms



% The outputs are:



% tstat  -  studentized coefficient test
% sig2   -  estimated error variance



% Print-to-screen is determined by the global variable _urprint,
% documented above.



% 



% PROCEDURE CODE: %
function [omega,bandw]=lr_var(u);
global kernel_;
global band_;
global urprint_;
global white_;
kernel_  = 1;
band_    = 0;
urprint_ = 0;
white_   = 0;



[tu,p]=size(u);
if white_==1
    te=tu-1;
    au=(u(2:tu,:)'/u(1:te,:)')';
    e=u(2:tu,:)-u(1:te,:)*au;
else
    e=u;
    te=tu;
end;



if band_==0
    eb=e(1:te-1,:);
    ef=e(2:te,:);
    ae=sum(eb.*ef)'./sum(eb.^2)';
    ee=ef-eb.*(ones(length(eb(:,1)),1)*ae');
    se=mean(ee.^2)';
    temp=(1-ae).^2;
    ad=sum((se./((1-ae).^2)).^2)';
    aese=ae.*se;
    a1=4*sum((aese./((((1-ae).^3).*(1+ae))*ones(1,length(aese(1,:))))).^2)'/ad;
    a2=4*sum((aese./(((1-ae).^4)*ones(1,length(aese(1,:))))).^2)'/ad;
    (aese./(((1-ae).^4)*ones(1,length(aese(1,:))))).^2
    
    if kernel_==2
        bandw = 1.3221*((a2*te)^.2); % Quadratic Spectral %
    elseif kernel_==1
        bandw = 2.6614*((a2*te)^.2); % Parzen %
        if bandw>(te-2)
            bandw=te-2;
        end;
    elseif kernel_==3;               % Bartlett %
        bandw=1.1447*((a1*te)^.333);
        if bandw>(te-2)
            bandw=te-2;
        end;
    end;
else
    bandw=band_;
end;



% Estimate Covariances %        
if kernel_==2  % Quadratic Spectral Kernel %
    tm=te-1;
    jb=((1:tm)/bandw')';
    jband = jb*1.2*pi;
    kern = ((sin(jband)./jband - cos(jband))./(jband.^2)).*3;
elseif kernel_==1 % Parzen kernel %
    tm=floor(bandw);
    if tm>0
        jb=((1:tm)/bandw')';
        kern=(1-(jb.^2)*6+(jb.^3)*6).*(jb <=.5);
        kern=kern + ((1-jb).^3).*(jb>.5)*2;
    end;
elseif kernel_==3 % Bartlett kernel %
    tm=floor(bandw);
    if tm>0
        kern=1-((1:tm)/bandw')';
    end;
end;



lam=zeros(p,p);
j=1;
while j<=tm
    kj=kern(j);
    lam=lam+(e(1:te-j,:)'*e(1+j:te,:))*kj;
    j=j+1;
end;
omega=(e'*e+lam+lam')/te;



if white_==1
    eau=inv(eye(p)-au);
    omega=eau'*omega*eau;
end;






function ct=ur_crits(r2,p);
global kernel_;
global band_;
global urprint_;
global white_;
if p==-1;
    crt=[-2.4611512       -1.7832090       -1.4189957
         -2.4943410       -1.8184897       -1.4589747
         -2.5152783       -1.8516957       -1.5071775
         -2.5509773       -1.8957720       -1.5323511
         -2.5520784       -1.8949965       -1.5418830
         -2.5490848       -1.8981677       -1.5625462
         -2.5547456       -1.9343180       -1.5889045
         -2.5761273       -1.9387996       -1.6020210
         -2.5511921       -1.9328373       -1.6128210
         -2.5658          -1.9393          -1.6156];;
elseif p==0
    crt=[-2.7844267       -2.1158290       -1.7525193
         -2.9138762       -2.2790427       -1.9172046
         -3.0628184       -2.3994711       -2.0573070
         -3.1376157       -2.5070473       -2.1680520
         -3.1914660       -2.5841611       -2.2520173
         -3.2437157       -2.6399560       -2.3163270
         -3.2951006       -2.7180169       -2.4085640
         -3.3627161       -2.7536756       -2.4577709
         -3.3896556       -2.8074982       -2.5037759
         -3.4336          -2.8621          -2.5671];
elseif p==1
    crt=[-2.9657928       -2.3081543       -1.9519926
         -3.1929596       -2.5482619       -2.1991651
         -3.3727717       -2.7283918       -2.3806008
         -3.4904849       -2.8669056       -2.5315918
         -3.6003166       -2.9853079       -2.6672416
         -3.6819803       -3.0954760       -2.7815263
         -3.7551759       -3.1783550       -2.8728146
         -3.8348596       -3.2674954       -2.9735550
         -3.8800989       -3.3316415       -3.0364171
         -3.9638          -3.4126          -3.1279];
end;



if r2<.1
    ct=crt(1,:);
else
    r210=r2*10;
    if r210>=10
        ct=crt(10,:);
    else
        r2a=floor(r210);
        r2b=ceil(r210);
        wa=r2b-r210;
        ct=wa*crt(r2a,:)+(1-wa)*crt(r2b,:);
    end;
end;



function [tstat,crits,rho2,sig2]=ur_reg(y,x,p,q,k1,k2);
global kernel_;
global band_;
global urprint_;
global white_;



t=length(y(:,1));
m=length(x(1,:));
ts=1+max([q;(k1-1)])';
tt=t-ts-k2;
dy=y(2:t,:)-y(1:t-1,:);
dv=dy(ts:length(dy(:,1))-k2,:);
z=y(ts:t-k2-1,:);
if p>=0
    z=[z,ones(tt,1)];
    if p==1
        z=[z,(1:tt)'];
    end;
end;



qq=1;
while qq<=q
    z=[z,dy(ts-qq:length(dy(:,1))-qq-k2,:)];
    qq=qq+1;
end;



z=[z,x(ts+1:t-k2,:)];



kk=1;
while kk<=k1
    z=[z,x(ts-kk+1:m-kk-k2,:)];
    kk=kk+1;
end;



kk=1;
while kk<=k2
    z=[z,x(ts+kk+1:t-k2+kk,:)];
    kk=kk+1;
end;



zzi=inv(z'*z);
beta=zzi*(z'*dv);
e=dv-z*beta;
sig2=(e'*e)/(tt - 2 - (1+k1+k2)*m - q - p);
cov=zzi*sig2;
se=sqrt(diag(cov));
tstat=beta(1)/se(1);
v=e+z(:,3+p+q:length(z(1,:)))*beta(3+p+q:length(z(1,:)));
[omega,bandw]=lr_var([v,e]);
r2=omega(2,2)/omega(1,1);
rho2=(omega(2,1)^2)/(omega(1,1)*omega(2,2));



crits=ur_crits(rho2,p);



if urprint_==0
    bs=[beta,se];
    fprintf('Regression Results\n\n');
    fprintf('Time Trend Order  %u\n',p);
    fprintf('Number AR Lags  %u\n',q);
    fprintf('Number of Lag DXs  %u\n',k1);
    fprintf('Number of Lead DXs %u\n',k2);
    if kernel_==1
        fprintf('Kernel = Parzen\n\n');
    elseif kernel_==2
        fprintf('Kernel = Quadratic Spectral\n\n');
    elseif kernel_==3
        fprintf('Kernel = Bartlett\n\n');
    end;
    fprintf('Y(t-1):\n');
    for j=1:length(bs(1,:))
            fprintf('%f   ',bs(1,j));
        end
        fprintf('\n');
    if p>=0
        fprintf('Constant, Trend:\n');
        for i=2:p+2
            for j=1:length(bs(1,:))
                 fprintf('%f   ',bs(i,j));
            end
        fprintf('\n');
        end;
        fprintf('\n');
    end;
    fprintf('DY(t-1) ... DY(t-q):\n');
    for i=p+3:p+2+q
        for j=1:length(bs(1,:))
            fprintf('%f   ',bs(i,j));
        end
        fprintf('\n');
    end;
    fprintf('\n');
    
    fprintf('X(t):\n');
    for i=p+3+q:p+2+q+m
        for j=1:length(bs(1,:))
            fprintf('%f   ',bs(i,j));
        end
        fprintf('\n');
    end;
    fprintf('\n');
    
    if k1>0
        fprintf('X(t-1) ... X(t-k1):\n');
        for i=3+p+q+m:2+p+q+m+k1*m
            for j=1:length(bs(1,:))
            fprintf('%f   ',bs(i,j));
        end
        fprintf('\n');
        end;
        fprintf('\n');
    end;
    if k2>0
        fprintf('X(t+1) ... X(t+k2):\n');
        for i=3+p+q+m+k1*m:length(bs(:,1))
            for j=1:length(bs(1,:))
            fprintf('%f   ',bs(i,j));
        end
        fprintf('\n');
        end;
        fprintf('\n');
    end;
    fprintf('\n');
    fprintf('Error Variance %f\n',sig2);
    fprintf('Studentized Unit Root Test Statistic %f\n\n',tstat);
    fprintf('1%, 5%, 10% Asymptotic Critical Values %f   %f   %f\n',crits(1),crits(2),crits(3));
    fprintf('Squared Correlation Measure (Rho^2) %f\n',rho2);
    fprintf('Variance Ratio (R^2)  %f\n',r2);
    fprintf('Bandwidth  %f\n\n',bandw)
end;





function [tstat,sig2]=ur_adf(y,p,q);
global kernel_;
global band_;
global urprint_;
global white_;
t=length(y(:,1));
tt=t-q-1;
dy=y(2:t,:)-y(1:t-1,:);
dv=dy(q+1:length(dy(:,1)),:);
z=y(q+1:t-1,:);
if p>=0
    z=[z,ones(tt,1)];
    if p==1
        z=[z,(1:tt)'];
    end;
end;
qq=1;
while qq<=q
    z=[z,dy(q-qq+1:length(dy(:,1))-qq,:)];
    qq=qq+1;
end;



zzi=inv(z'*z);
beta=zzi*(z'*dv);
e=dv-z*beta;
sig2=(e'*e)/(tt-2-q-p);
cov=zzi*sig2;
se=sqrt(diag(cov));
tstat=beta(1)/se(1);



if urprint_==0;
    bs=[beta,se];
    fprintf('Regression Results\n\n');
    fprintf('Y(t-1):\n');
    fprintf('%f\n\n',bs(1,:));
    if p>=0
        fprintf('Constant, Trend:\n');
        for i=2:p+2
            for j=1:length(bs(1,:))
                 fprintf('%f   ',bs(i,j));
            end
            fprintf('\n');
        end;
        fprintf('\n');
    end;
    if q>0
        fprintf('DY(t-1) ... DY(t-q):\n');
        for i=p+3:p+2+q
            for j=1:length(bs(1,:))
                 fprintf('%f   ',bs(i,j));
            end
            fprintf('\n');
        end;
    fprintf('\n');
    end;
    fprintf('Error Variance  %f\n\n',sig2);
    fprintf('t Statistic %f\n\n',tstat);
end;