@  Written/Modified by: Junsoo Lee

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

/*  Zivot and Andrews (1991, JBES) test

program: za.txt

Note: This program first determines the optimal lag for each of all possible
  cases of break points, and then search for optimal break points.
  Perron (1994) used the same method.  On the other hand, some applied papers
  determine the optimal break point first with lag = 0, and then determine the 
  optimal lag later by using the estimated break point with lag = 0. 
  Such method is not recommended.

*/

/*  Options to choose */

 output file = za.out reset;
  "ZA Min t-stat ";

 maxk = 8;
     @ # maximum augmentation lags allowed for @

 n = 100;
     @ # of observations @

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

 T1_end = .1;   T2_end = .9;
     @ End points (10% of both ends are not considered) @

 crit_t = 1.645;
     @ critical value of the t-stat for lag determination @


/* End of Options */

 @ imodel  1 = Crash model (dummy on the intercepts)
           2 = Breaking trend (dummy on the trend)   @
 @ option  1 = TB(t) included in the testing regression
           2 = excluded as in ZA  @

 imodel=1;
 do while imodel <=2;

 option=1;
 do while option <=2;

    T = n;
    T1 = round(T1_end * rows(y));
    T2 = round(T2_end * rows(y));

    lag = 0;

    if T1 < maxk+2; T1 = maxk + 3;  endif;

    westd1 = zeros(T,1);
    wtd1 = zeros(T,1);     wdf1 = zeros(T,1);
    wtk1 = zeros(T,1);     woptk1 = zeros(T,1);


  ii = T1;
  Do while ii <= T2;
     output off; "*";;

     Estconf = ii;

     wtk = zeros(maxk+1,1);
     wdf = zeros(maxk+1,1);
     westd = zeros(maxk+1,1);
     wtd = zeros(maxk+1,1);

     kk = 0;  @ lag @
     do while kk <= maxk;

     Estconf = ii;

     If imodel == 1;
      if option == 1;
       {wtstat, wed, wtd0, wres, wbeta, wtval, sig} = p_alt(y, 1, EstConf,kk);
          @ z(t) = [y(t-1), lags, 1, t, B(t), D(t)] @
      endif;
      if option == 2;
       {wtstat, wed, wtd0, wres, wbeta, wtval, sig} = p_alt1(y, 1, EstConf,kk);
      endif;
     Endif;

     If imodel == 2;
      if option == 1;
       {wtstat, wed, wtd0, wres, wbeta, wtval, sig} = p_altC(y, 1, EstConf,kk);
          @ z(t) = [y(t-1), lags, 1, t, B(t), D(t)] @
      endif;
      if option == 2;
      {wtstat, wed, wtd0, wres, wbeta, wtval, sig} = p_altC1(y, 1, EstConf,kk);
      endif;
     Endif;

     westd[kk+1] = wed;
     wtd[kk+1] = wtd0;
     wdf[kk+1] = wtstat;
     if kk > 0; wtk[kk] = wtval[kk+1];  endif;
     @ t-stat of the coeff. of the last augmented term @

    kk = kk + 1;
    endo;

      jj = maxk;  isw = 0;
      do while isw ne 1;
        if (abs(wtk[jj]) > crit_t) or (jj == 0);
           woptk1[ii] = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

      tt = woptk1[ii]+1;
      wdf1[ii] = wdf[tt];

   ii = ii + 1;
 Endo;

    output on;

    wo4 = minindc1(wdf1[T1:T2]) + T1 - 1;
    optb = wo4; optk = woptk1[optb];

    format /rd/m1 10,4;
    ""; "-----------------------------------------"; "";

    "Model  (1=A, 2=C)      = " imodel;
    "Option (1=with, 2=w/o) = " option;

    if option == 1;
       If imodel == 1;
         {wts,wd,wtd,wres,wb,wt,sig} = p_alt(y,1,optb,optk);
       endif;
       If imodel == 2;
        {wts,wd,wtd,wres,wb,wt,sig} = p_altC(y,1,optb,optk);
       endif;
     endif;

    if option == 2;
       If imodel == 1;
         {wts,wd,wtd,wres,wb,wt,sig} = p_alt1(y,1,optb,optk);
       endif;
       If imodel == 2;
        {wts,wd,wtd,wres,wb,wt,sig} = p_altC1(y,1,optb,optk);
       endif;
     endif;

     "";
     "Min. test statistic    = "  wts;
     "Estimated break point  = "  optb;
     "Selected lag           = "  optk;

     "";
     "Estimated coeff. of dummy var. = "  wd;
     "Its t-stat                     = "  wtd;
     "Standard error ..              = "  sig;
     "Standardized dummy coeff       = "  wd/sig;

     if imodel == 1;
      ""; "Coeff and t-stat";
      if option == 1;
        "X(t) = [y(t-1), (lags..omitted), 1, t, B(t), D(t)] ";
      endif;
      if option == 2;
        "X(t) = [y(t-1), (lags..omitted), 1, t, D(t)] ";
      endif;
     endif;

     if imodel == 2;
      ""; "Coeff and t-stat";
      if option == 1;
       "X(t) = [y(t-1), (lags..omitted), 1, t, B(t), D(t), DT(t)] ";
      endif;
      if option == 2;
       "X(t) = [y(t-1), (lags..omitted), 1, t, D(t), DT(t)] ";
      endif;
     endif;

     tt1 = wb[1] | wb[2+optk:rows(wb)];
     tt2 = wt[1] | wt[2+optk:rows(wt)];
     tt3 = tt1 ~ tt2;
     tt3; "";



option = option + 1;
endo;

imodel = imodel + 1;
endo;

/*
Format: {wts,wd,wtd,wres,wb,wt,sig} = p_alt(y,p,optb,optk);

Input:  x  = series
        p  = 1 (with trend), 0 (with drift), -1 (no drift, no trend)
        optb = position of the structural break
        optk = # of augmentation lags

Output: wts  = ADF t-statistic
        wd   = estimated coefficients of break dummies
        wtd  = t-stat for coefficients of break dummies
        wres = residuals
        wb   = all estimated coefficients
        wt   = t-stat of estimated coefficients
        sig  = standard error of regression
*/

proc (7) = p_alt(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     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);@
     dep     = diff(x,1);

     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 ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z)) ;

  @ Here @

     If (config ne 0);

       R = rows(config);

       pos = config[1];
       du = zeros(pos,1) | ones(nobs-pos,1);

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

       i=2;
       do while i <= R;
          pos = config[i];
          temp = zeros(pos,1) | ones(nobs-pos,1);
          temp2 = zeros(nobs,1);
          temp2[pos+1,1] = 1;
          @du = du ~ temp ~ temp2;@
          du = du ~ temp2 ~ temp;
          i = i + 1;
       endo;

       z1 =  du;
       dz = diff(z1,1);

     Endif;

    @ z1 = trimr(z1,1,0);
      z1 = z1[1:rows(z1)-k-1,.];@

     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     @z[1:10,.]';@

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;

proc (7) = p_alt1(x,p,config,l) ;

     @ Without B(t) @
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     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);@
     dep     = diff(x,1);

     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 ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z)) ;

  @ Here @

     If (config ne 0);

       R = rows(config);

       pos = config[1];
       du = zeros(pos,1) | ones(nobs-pos,1);
       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       @du = temp2 ~ du;@

       i=2;
       do while i <= R;
          pos = config[i];
          temp = zeros(pos,1) | ones(nobs-pos,1);
          temp2 = zeros(nobs,1);
          temp2[pos+1,1] = 1;
          @du = du ~ temp ~ temp2;@
          @du = du ~ temp2 ~ temp;@
          du = du ~ temp;
          i = i + 1;
       endo;

       z1 =  du;
       dz = diff(z1,1);

     Endif;

    @ z1 = trimr(z1,1,0);
      z1 = z1[1:rows(z1)-k-1,.];@

     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     @z[1:10,.]';@

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;

proc (7) = p_altC(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     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);@
     dep     = diff(x,1);

     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 ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));

  @ Here @

     If (config ne 0);

       R = rows(config);

       pos = config[1];
       du = zeros(pos,1) | ones(nobs-pos,1);
       dt = zeros(pos,1) | seqa(1,1,nobs-pos);

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

       i=2;
       do while i <= R;
          pos = config[i];
          temp = zeros(pos,1) | ones(nobs-pos,1);
          temp2 = zeros(nobs,1);
          temp2[pos+1,1] = 1;
          t2 = zeros(pos,1) | seqa(1,1,nobs-pos);
          @du = du ~ temp ~ temp2;@
          du = du ~ temp2 ~ temp;
          dt = dt ~ t2;
          i = i + 1;
       endo;

       z1 =  du ~ dt;
       dz = diff(z1,1);

     Endif;

     @z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);@
     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;

proc (7) = p_altC1(x,p,config,l) ;

     @ Without B(t) @
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     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);@
     dep     = diff(x,1);

     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 ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));

  @ Here @

     If (config ne 0);

       R = rows(config);

       pos = config[1];
       du = zeros(pos,1) | ones(nobs-pos,1);
       dt = zeros(pos,1) | seqa(1,1,nobs-pos);

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       @du = temp2 ~ du;@

       i=2;
       do while i <= R;
          pos = config[i];
          temp = zeros(pos,1) | ones(nobs-pos,1);
          temp2 = zeros(nobs,1);
          temp2[pos+1,1] = 1;
          t2 = zeros(pos,1) | seqa(1,1,nobs-pos);
          @du = du ~ temp ~ temp2;@
          @du = du ~ temp2 ~ temp;@
          du = du ~ temp;
          dt = dt ~ t2;
          i = i + 1;
       endo;

       z1 =  du ~ dt;
       dz = diff(z1,1);

     Endif;

     @z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);@
     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;


proc (1) = minindc1(x);
  local d,m,i,minx;

  m = 1;
  d = x[1];

  i=1;
  do while i <=rows(x);
    if x[i] < d;  d = x[i];  endif;
    i=i+1;
  endo;

 i=2;
  do while i <= rows(x);
    if d == x[i]; m = i;  goto stops; endif;
    i = i + 1;
  endo;
  stops:
retp(m);
endp;