/*--------------------------------------------------------------------------
                                                                           
  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];
    dy  = y[2:n1] - y_1;
    n2  = rows(dy);
    bic = 100000; 
    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));  
        tv  = b./seb;           
        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;