@  Written/Modified by: Junsoo Lee

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

/*  Perron (1989, Econometrica) test

program: perron.txt

*/

/*  Options to choose */

 output file = perron.out reset;
  "Perron's 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;  @

 pbreak = 50;
     @ break point; between 1 and n-1 @

/* 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
           (it will be a mistake to exclude the term, though) @

 imodel=1;
 do while imodel <=2;

 option=1;
 do while option <=2;

     Estconf = 0;

     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 = pbreak;

     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 = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

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

    optk = woptk1;
	optb = pbreak;

    output on;

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

    "Model  (1=A, 2=C)      = " imodel;
    "Option (1=w/o, 0=with) = " 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;
     "Given break point      = "  optb;
     "Selected lag           = "  optk;
     "";
     "Estimated coeff. of dummy var. = "  wd;
	 "Its t-stat                     = "  wtd;

     "Standard error ..      = "  sig;
     "Standardized dummy coeff "  wd/sig;

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;