@  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 = sp1.txt @
@  **  SP (Schmidt and Phillips, Oxford 1992)        @
@  **  SL (Schmidt and Lee, Economics Letters, 1991) @

@  Transformation tests using Fixed truncation lag   @
@  To use an automatic bandwidth selection, COINT 2.0 is needed @

new;

/*  Options to choose */

 output file = sp1.out reset;

 lag = 8;
     @ # truncation lag to estimate the longrun variance @

 n = 100;
     @ # of observations @

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

/* End of Options */


/* SP or SL Test */

lag = 0;
do while lag <= 8;

   "";"";
   "Truncation lag  = " lag; "";

   "  Schmidt-Phillips Tests"; "";
   select = 1;

   {alpha, Za, Zt} = SP(y, select, lag);

   format /rd/m1 12,4;
   "    T(a-1)          = " alpha;
   "    Za stat         = " Za;
   "    Zt stat         = " Zt;  "";


   "  Schmidt-Lee Tests"; "";
   select = 0;

   {alpha, Za, Zt} = SP(y, select, lag);

   format /rd/m1 12,4;
   "    T(a-1)          = " alpha;
   "    Za stat         = " Za;
   "    Zt stat         = " Zt;

   lag = lag + 1;
  endo;

end;

/*
**
**  SP (Schmidt and Phillips, Oxford 1992)
**  SL (Schmidt and Lee, Economics Letters, 1991)
**
** Purpose: Computes Z(rho) or Z(tau) statistics of SP (and SL) allowing for
**          deterministic trends of an arbitrary order.
** Transformation tests only (not augmented tests)
**
** Format:  {alpha, Za, Zt} = SP(x,p,l);
**
** Input:   x     -- time series
**
**          p     -- This is NOT "order of the time-polynomial to include
**                   in the fitted regression."  Lee re-defined p.
**                   Set p = 1 for Schmidt and Phillips test (w/ a constant).
**                       p = 0 for Schmidt and Lee (1991) test (w/o a const).
**
**          l     -- number of auto-correlations to use when computing the
**                   spectrum(@0).
**
**
** Output:  alpha       --  estimate of the autoregressive parameter;
**
**          Za          --  Z(rho) statistic
**
**          Zt          --  Z(tau) statistic
**
*/
proc (3) = SP(x,p,k);
     local lr, res, lhat, alpha, za, zt, nobs, fhat, xihat, st, dx, s1, rho;
     local sse, tau, ratios, sum2, is, it, wsl, lag1, e, n;

     if (p < 0);
        "Error: p cannot be set < 0.  Set p=1 for SP, and p=0 for SL.";
         retp(0,0,0,zeros(6,1),zeros(6,1));
     endif ;
     if (p > 1);
        "Error: p cannot be set > 1.  Set p=1 for SP, and p=0 for SL.";
         retp(0,0,0,zeros(6,1),zeros(6,1));
     endif ;

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

     nobs    = rows(x) ;
     fhat    = (x[nobs,1] - x[1,1]) / (nobs - 1);  /* mean of y */
     xihat   = x[1,1] - fhat;
     st      = x - xihat - fhat*seqa(1,1,nobs);

     if p == 0;   /* SL test */
         dx  = diff(st,1);
         s1  = trimr(lagn(st,1),1,0);
     else;        /* SP test */
         dx  = diff(x,1) - fhat;  /* demeaned */
         s1  = trimr(lagn(st,1),1,0);
         s1  = detrend(s1,p-1);
     endif;

     alpha   = dx/s1 ;
     rho     = nobs * alpha;
     res     = dx - (s1*alpha);
     sse     = res'res/(rows(res));
     tau     = alpha/sqrt(sse/(s1's1));


      e = res;
      lag1 = lag+1;
      n = rows(res);

      it=1; is=1;
      sum2=0;
      do while is < lag1;
        it=is+1;
        wsl=1-is/lag1;
        do while it < n+1;
          sum2=sum2+wsl.*e[it,1].*e[it-is,1];
          it=it+1;
        endo;
        is=is+1;
      endo;

      sum2=2*sum2/rows(res);
      lr = sse + sum2 ;

     /* lr      = lrvar(res,k);
     lhat    = covarf(res,k);  @ lr = sse + 2*lhat @ */

     ratios  = sse / lr;
     Za      = rho / ratios;
     Zt      = tau / sqrt(ratios);

     retp(alpha,za,zt);
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 ;

proc detrend(data,p) ;
local t, u, m, timep, it, iit, xmat, invx, beta, resid, nobs ;
     if (p == -1) ;
         retp(data);
     endif ;
     nobs     = rows(data) ;
     u        = ones(nobs,1) ;
     if p > 0 ;
       timep    = zeros(nobs,p) ;
       t        = seqa(1,1,nobs)/nobs ;
       m        = 1 ;
       do while m <= p ;
          timep[.,m] = t^m ;
          m = m + 1 ;
       endo ;
       xmat     = u~timep ;
     else ;
       xmat = u ;
     endif ;
     invx     = inv(xmat'xmat) ;
     beta     = invx*(xmat'data)  ;
     resid    = data - xmat*beta  ;
retp(resid) ;
endp ;