@  Written/Modified by: Junsoo Lee

   Disclaimer: This code is provided gratis without any guarantees 
               or warrantees.  Proprietary modifications of this code 
			   are not permitted.
@

@  file = adf.txt @
@  Augmented Dickey Fuller Tests @
@  Using fixed lags or data-driven lag selection procedures @

new;

/*  Options to choose */

 output file = adf.out reset;

 maxk = 8;
     @ max # lags allowed in the data-driven procedure @

 n = 100;
     @ # of observations @

 crit_t = 1.645;
     @ critical value used for significance of lagged terms @

 p = 1;
     @ 1 = with trend, 0 = with constant only, -1 = no constant or trend @

 select = 0;
     @ 1 = fixed lag, 0 = data dependent lag selection @

 y = cumsumc(rndn(n,1));
     @ Or, Read the dat file here.
        ->  load data[n,1] = data.txt;  @

/* End of Options */

/* Using a Fixed lag */
if select == 1;

   lag = 4;
   {beta,tval, alpha, tstat} = adf(y,p,lag);

   format /rd/m1 12,4;
   "Using a Fixed lag"
   "Fixed lag    = " lag;
   "T(a-1)       = " alpha;
   "ADF stat     = " tstat;
   "Estimated coefficients :
     (x(t-1), dx(t-1), dx(t-2), ... dx(t-k), const, trend) ";
   beta';
   "t-stat";
   tval';

   goto next;

endif;

/* Using data-driven lag selection procedure
   'general to specific' search
    see Ng and Perron (1995, JASA) or Lee (1997, WP) */

 @ maxk # regressions @

     statkeep = zeros(maxk+1,1);
     tvalkeep = zeros(maxk,1);

     {beta,tval, alpha, tstat} = adf(y,p,0);
     statkeep[1] = tval[1];

     j=1;
     do while j <= maxk;

        {beta,tval, alpha, tstat} = adf(y,p,j);

        statkeep[j+1] = tval[1];
        tvalkeep[j] = tval[j+1];

        j = j + 1;
     endo;

     /* statkeep  for k = 0, 1, 2, ..., maxk */
     /* tvalkeep  for k = 1, 2, ..., maxk */

@ selection using t-test @

     /* select the lag where the coeff is sig. from maxk to 0.
        if all coeff are insig., kk = 0  */

     j = maxk;
     do while j >= 0;
        if (abs(tvalkeep[j]) > crit_t) or (j == 0);
            kk = j;
            j = 0;
        endif;
        j = j - 1;
     endo;

   @ kk = selected lag @

   {beta, tval, alpha, tstat} = adf(y,p,kk);

   format /rd/m1 12,4;
   "Using Data dependent selection of the lag";
   "Selected lag = " kk;
   "T(a-1)       = " alpha;
   "ADF stat     = " tstat;
   "Estimated coefficients :
     (x(t-1), dx(t-1), dx(t-2), ... dx(t-k), const, trend) ";
   beta';
   "t-stat";
   tval';

next:
end;

/*
** ADF
**
** Purpose: Compute the Augmented Dickey Fuller statistic, allowing
**          for deterministic polynomial time trends of an arbitrary
**          order.
**
** Format:  {beta, tval, alpha, tstat} = adf(x,p,l);
**
** Input:   x     -- time series variable
**
**          p     -- order of the time-polynomial to include in the
**                   ADF regression.  Set p = -1 for no deterministic
**                   part.
**
**          l     -- number of lagged changes of x to include in the
**                   fitted regression.
**
** Output:  alpha       --  estimate of the autoregressive parameter;
**
**          tstat       --  ADF t-statistic
**
**          beta        --  coeff. of (x(t-1), dx(t-1), dx(t-2), ... dx(t-k),
**                                     constant, time trend)
**          tval        --  t-values
*/

proc (4) = adf(x,p,l) ;
     local b,k,z,res,so,var_cov,xx, tval;
     local timep,t,m,xmat,nobs,dep,ch,f;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;

     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);
     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;

     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);
     if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;
     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     b[1,1] = b[1,1] - 1;
     tval    = b ./ sqrt(diag(var_cov));
retp(b,tval,b[1,1],(b[1,1])/sqrt(var_cov[1,1])) ;
endp ;

proc lagn(x,n);
    local y;
    y = shiftr(x', n, (miss(0, 0))');
    retp(y');
endp;

proc diff(x,k) ;
     if ( k == 0) ;
        retp(x) ;
     endif ;
     retp(trimr(x,k,0)-trimr(lagn(x,k),k,0)) ;
endp ;