/*--------------------------------------------------------------------------
@                                                                           
@  ADFABIC.PRC                                                              
@  Purpose: Test the null of unit root process by ADF                       
@           Number of lags chosen by AIC or BIC      	                    
@  Written by Hyeongwoo Kim (Feb 29, 2008)                                  
@---------------------------------------------------------------------------
@                                                                           
@  Format: {lag,adf,b,std} = adfabic(y,pmax,c,crt);                          
@									    
@  Model1: dy(t) = b_0*y(t-1)                                 		       
@                   + b_1*dy(t-1) + ... + b_(p)*dy(t-p) + e(t)              
@  Model2: dy(t) = b_0*y(t-1) + b_1                                         
@                   + b_2*dy(t-1) + ... + b_(p+1)*dy(t-p) + e(t)            
@  Model3: dy(t) = b_0*y(t-1) + b_1 + b_2*t                                 
@                   + b_3*dy(t-1) + ... + b_(p+2)*dy(t-p) + e(t)            
@                                                                           
@  Input : y    (TX1) Vector of a Time Series                                
@          pmax (1x1) Maximum number of lags (differenced)  	    	    
@          c=0  (1x1) No constant                                            
@           =1        Constant                                              
@           =2        Constant and linear time trend                         
@          crt=1      Use the BIC rule
@             =2      Use the AIC rule
@                                                                           
@  Output: lag        Chosen number of lags                         	    
@          adf        ADF statistic                                         
@          b          Coefficient estimates                                 
@          std        Standard Errors                                       
@          It also prints ADF test results                                  
@                                                                           
--------------------------------------------------------------------------*/


proc(4) = adfabic(y,pmax,c,crt);
 local n1,y_1,dy,n2,p,j,x,non,tvb,cn,i,b,rsd,ssq,bic,sic,seb,vcb,tv,adf,rho,std,lag;

n1 = rows(y);

y_1 = y[1:n1-1];            @ first lagged y: y(t-1)        @
dy  = y[2:n1] - y_1;        @ differenced y: dy(t)          @

n2  = rows(dy);

bic = 100000; adf = 100; rho = 100; std = 100; lag = 100;		//initialization

j = pmax; do until j < 0;

if c eq 1;  x = y_1[j+1:n2]~ones(n2-j,1); 
 elseif c eq 2; x = y_1[j+1:n2]~ones(n2-j,1)~seqa(1,1,n2-j);
 else; x = y_1[j+1:n2];
endif;

  i = 1; do until i > j;
   x = x~dy[j+1-i:n2-i];
   i = i+1;
  endo;

b = inv(x'x)*x'dy[j+1:n2];
rsd = dy[j+1:n2] - x*b;
non = rows(x);

if crt eq 1; cn = ln(non); endif;
if crt eq 2; cn = 2; endif;

ssq = rsd'rsd;
sic = ln(ssq/non) + cn*(j+1)/non;
vcb = (ssq/non)*inv(x'x);
seb = sqrt(diag(vcb));                  @ standard errors of b's        @
tv  = b./seb;                           @ t-values                      @

if sic le bic; bic = sic; adf = tv[1]; rho = b; std = seb; lag = j; endif;

j = j - 1;

endo;

"";
"-------------------------------------------------------------------";"";
"One-sided test of H0: Unit root vs. H1: Stationary";
"Lag selection by general to specific";
"Approximate asymptotic critical values (t-ratio):";"";
"------------------------------------------------------------";
"  1%      5%      10%      Model";
"------------------------------------------------------------";
"-2.56   -1.94   -1.62";;
"     Simple ADF (no constant or trend)";
"-3.43   -2.86   -2.57";;
"     ADF with constant (no trend)";
"-3.96   -3.41   -3.13";;
"     ADF with constant & trend";
"------------------------------------------------------------";"";

format/M1/RDN 12,3;

"Chosen number of lags        ";;lag;
"ADF Statistic (t-ratio)      ";;adf;
"Rho and its standard error   ";;rho[1]~std[1];"";

"-------------------------------------------------------------------";"";

 retp(lag,adf,rho,std);

endp;