/************************************************************************
       ADFtest.txt
       Performs Monte Carlo studies of unit root test
       Version: Nov 17, 2016
*************************************************************************/
output file=ADFmc.out reset ;

Niter = 1000; "Number of iterations: ";; Niter ;

Length = 100 ; "Length of series: ";; Length ;

ARorder = 3 ; yzero = zeros(Arorder,1) ;

Arcoeff = {.5,.3,.02} ; "AR coeff: " ;; Arcoeff' ;

/***** order selection for ADF test ************/

ordertp = 1 ;

/* ordertp = 0 for mean estimated; = 1 for trend and mean estimated */

ahat = zeros(Niter,1);
tstat = zeros(Niter,1);

i = 0 ;
do until i == Niter ;
  i = i+1 ;
  /**************** simulate AR-process ********************/
  xsr =recserar(rndn(Length+1,1),yzero,Arcoeff)  ;


  /************* ADF test **********************************/
 @ {ahat,tstat[i],c_t} =ADF(xsr,3) ;@
  {ahat[i],tstat[i]} =ADF(xsr,3) ;


  i ;
endo ;
frac1 = sortc(ahat,1) ;
"fractiles for estimated parameter" ;
"1%";;  (frac1[.01*Niter])' ;
"5%";;  (frac1[.05*Niter])' ;
"95%";; (frac1[.95*Niter])' ;
"99%";; (frac1[.99*Niter])' ;

"fractiles for ADF-test " ;
frac5 = sortc(tstat,1) ;
"1%";;  (frac5[.01*Niter])' ;
"5%";;  (frac5[.05*Niter])' ;
"95%";; (frac5[.95*Niter])' ;
"99%";; (frac5[.99*Niter])' ;

c_t=-3.45; @5 percent critical value, one sided@

 "ADF 5% critical value: " ;; c_t ;

 if sumc(Arcoeff)==1 ;
"size of ADF-test (per cent): ";; sumc(c_t .GE frac5)/Niter ; else
;

"power of ADF-test (per cent): ";; sumc(c_t .GE frac5)/Niter ;
endif ;

@========================================================@
/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@ @ !! ADF(Y,ADFLAG): 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 @ @ 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 ADF(Y,ADFLAG) @ @ Y: The data to be analyzed.
Null is that Y has a unit root @ @ ADFLAG: Lags to put in the ADF
regression @
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/

 proc(2)=ADF(y,adflag) ;

 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; */

 retp(phi[2,1]+1,tau) ;
@-------------------------------------------------------------@
endp;
@=============================================================@

output off ;