@========================================================@
@::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @
@ !! UROOT(Y,ADFLAG,NWLAG): This proc computes                         @
@   1. Studentized coefficient for ADF test of unit root, no trend     @
@   2. Studentized coefficient for ADF test of unit root, with trend   @
@   3. t-test for Phillips-Perron test of unit root, no trend          @
@   4. t-test for Phillips-Perron test of unit root, with trend        @
@                                                                      @
@   ALL TESTS INCLUDE A CONSTANT                                       @
@     for no trend results consult Hamilton table B.6, case 2          @
@     for trend results, consult Hamilton table B.6, case 3            @
@                                                                      @
@ USAGE: CALL UROOT(Y,ADFLAG,NWLAG)                                    @
@   Y: The data to be analyzed. Null is that Y has a unit root         @
@   ADFLAG: Lags to put in the ADF regression                          @
@   NWLAG:  Lags to put in Newey-West Estimator for Phillips-Perron    @
@ The proc writes the results to the current output file               @
@   and does not return anything                                       @
@ :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @
 proc(0)= uroot(y,adflag,nwlag) ;

 local sig,sigsq,talpha,lam1,lam2,mbaryy,ztau,ft,nmom,n2,w,s,t,
       ser,a2,tau,sto,nobs,ser,varreg,resid,phi,i,j,k,t0,t1,dy,x,
       myy,mty,my,m,zalpha,c3;

 t1=rows(y);      dy=zeros(t1,1);

 dy[2:t1,1]=y[2:t1,1]-y[1:t1-1,1];     @ First difference @


 j=adflag+2;   @ CONSTRUCT X MATRIX FOR ADF REGRESSION--NO TREND @
  x=ones(rows(y[j:t1,1]),1)~y[j-1:t1-1,1]~dy[j-1:t1-1,1];
  k=2; do while k <= adflag;
     x=x~dy[j-k:t1-k,1];
     k=k+1;
  endo;


 nobs=rows(x);         @ RUNNING THE ADF REGRESSION     @
 phi=inv(x'x)*x'*dy[j:t1,1];
 resid=dy[j:t1,1]-x*phi;        @ ADF regression residuals @
 varreg=(resid'*resid)/(nobs-cols(x));
 ser=sqrt(varreg);              @ Std. Error of regression @
 sto=varreg*inv(x'x);
 tau=phi[2,1]/sqrt(sto[2,2]);

"ADF t-stat for Unit Root (No Trend)   ";;tau;; adflag;;" Lags";;
" rho";;phi[2,1]+1;



@====================================================@

 j=adflag+2;   @ CONSTRUCT X MATRIX FOR ADF REGRESSION--WITH TREND @
  x=ones(rows(y[j:t1,1]),1)~y[j-1:t1-1,1]~dy[j-1:t1-1,1]
         ~(seqa(1,1,nobs)-nobs/2);
  k=2; do while k <= adflag;
     x=x~dy[j-k:t1-k,1];
     k=k+1;
  endo;

 phi=inv(x'x)*x'dy[j:t1,1];    @ Run the ADF regression @
 resid=dy[j:t1,1]-x*phi;        @ ADF regression residuals @
 varreg=(resid'*resid)/(nobs-cols(x));
 ser=sqrt(varreg);              @ Std. Error of regression @
 sto=varreg*inv(x'x);
 tau=(phi[2,1])/sqrt(sto[2,2]);

"ADF t-stat for Unit Root (With Trend) ";;tau;; adflag;;" Lags";;
" rho";;phi[2,1]+1;

@---------------- COMPUTATIONS FOR PHILLIPS-PERRON--------------@
if nwlag>=0;
  j=2;  @ CONSTRUCT X MATRIX FOR PP REGRESSION--NO TREND @
  x=ones(rows(y[j:t1,1]),1)~y[j-1:t1-1,1];

 nobs=rows(x);         @ RUNNING THE PP REGRESSION     @
 phi=inv(x'x)*x'y[j:t1,1];
 resid=y[j:t1,1]-x*phi;        @ ADF regression residuals @
 varreg=resid'*resid/nobs;
 ser=sqrt(varreg);              @ Std. Error of regression @
 sto=meanc(y[j-1:t1-1,1]);
 a2=sumc((y[j-1:t1-1,1]-sto)^2);  a2=sqrt(a2);
 talpha=(phi[2,1]-1)*a2/ser;

ft=resid';                @ CONSTRUCT FT @
nmom=rows(ft);
n2=cols(ft);

   w=zeros(nmom,nmom) ;    @ DOING THE NEWEY AND WEST THING @
   t=1;                    @ COMPUTING SIGSQ @
   do while t <= n2 ;   w = ft[.,t]*ft[.,t]' + w ;  t=t+1 ; endo ;
   s=w;
  i=1 ;
  do while i <= nwlag ;  w=zeros(nmom,nmom) ;
   t=1; do while t <= n2-i;  w=ft[.,t]*ft[.,t+i]' + w ; t=t+1 ; endo ;
  s = s + (w + w')*(1-(i/(nwlag+1)));
  i=i+1 ;
  endo ; @-----------------------Spectral density at freq 0 computed---@


  sigsq=s/nobs;
  sig=sqrt(sigsq);
  sto=meanc(y);
  mbaryy=sumc((y-sto)^2); mbaryy=mbaryy/(nobs^2);
  lam1=(sigsq-varreg)/2;
  lam2=lam1/sigsq;
  ztau=(ser*talpha/sig) - lam2*sig/(sqrt(mbaryy));

"PP  t-stat for Unit Root (No Trend)   ";;ztau;; nwlag;;" Lags";;
" rho";;phi[2,1];

@-----------------------------------------------------------@

 j=2;  @ CONSTRUCT X MATRIX FOR PP REGRESSION--WITH TREND @
  x=ones(rows(y[j:t1,1]),1)~y[j-1:t1-1,1]~(seqa(1,1,nobs)-(nobs)/2);

      @ RUNNING THE PP REGRESSION     @
 phi=inv(x'x)*x'y[j:t1,1];
 resid=y[j:t1,1]-x*phi;        @ ADF regression residuals @
 varreg=resid'*resid/nobs;
 ser=sqrt(varreg);              @ Std. Error of regression @
 sto=inv(x'x);
 c3=sto[2,2];
 talpha=(phi[2,1]-1)/sqrt(ser*ser*c3);


ft=resid';                @ CONSTRUCT FT @
nmom=rows(ft);
n2=cols(ft);

   w=zeros(nmom,nmom) ;    @ DOING THE NEWEY AND WEST THING @
   t=1;                    @ COMPUTING SIGSQ @
   do while t <= n2 ;   w = ft[.,t]*ft[.,t]' + w ;  t=t+1 ; endo ;
   s=w;
  i=1 ;
  do while i <= nwlag ;  w=zeros(nmom,nmom) ;
   t=1; do while t <= n2-i;  w=ft[.,t]*ft[.,t+i]' + w ; t=t+1 ; endo ;
  s = s + (w + w')*(1-(i/(nwlag+1))) ;
  i=i+1 ;
  endo ; @-----------------------Spectral density at freq 0 computed---@

  sigsq=s/nobs;
  sig=sqrt(sigsq);
 myy=(nobs^(-2))*sumc(y[j:t1,1]^2);
 mty=(nobs^(-5/2))*sumc(seqa(1,1,nobs).*y[j:t1,1]);
 my=(nobs^(-3/2))*sumc(y[j:t1,1]);
 m=(1-nobs^(-2))*myy-12*mty^2+12*(1+(1/nobs))*mty*my
      -(4+(6/nobs)+2/(nobs^2))*my^2;

 lam1=(sigsq-varreg)/2;
 lam2=lam1/sigsq;
 ztau=(ser*talpha/sig) - lam2*sig/(sqrt(m));

"PP  t-stat for Unit Root (With Trend) ";;ztau;; nwlag;;" Lags";;
" rho";;phi[2,1];
endif;
@-------------------------------------------------------------@
endp;
@=============================================================@