@  Written/Modified by: Junsoo Lee

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

/*  Lumsdain and Papell (1997, Restat) minimum LM unit root test with two breaks

program: LPtwo.txt

  Also see LStwo.txt;
  Lee and Strazicich Min LM test with two breaks.

  Note: This program first determines the optimal lag for each of all possible
  cases of break points, and then search for optimal break points.

*/

/*  Options to choose */

 output file = LPtwo.out reset;

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

 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 @

 n = 50;
     @ # of observations @

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

 /*
    load data[43,1] = c:\work\help\mark\egypt.txt;
    @load data[27,1] = c:\work\help\mark\israel.txt;@

    let name = egypt;
    @let name = israel;@

    data = ln(data);

    jjj = 1;
    do while jjj <= cols(data);

    y = data[.,jjj];

    format /m1/rd 20,0 ;
    ""; "";
    "|========>  ";;
    "Series (logged) : "  $name[jjj,1];   "";

     n = rows(y);
     format /m1/rd 5,0;
    "# of obs = " n;
*/


/* End of Options */

 @ imodel  1 = Crash model (dummy on the intercepts)
           2 = Breaking trend (dummy on the trend)   @

 imodel=1;
 do while imodel <=2;

    cls;
    "";
    "---------------------------------------------------";
    "Lumsdain and Papell Min t-stat with two breaks ";
    "Model = " imodel;

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

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


    @ any big value @
    minlp1 = 1000;
    minlp2 = 1000;

    wtk1 = zeros(maxk+1,1);
    wdf1 = zeros(maxk+1,1);
    wtk2 = zeros(maxk+1,1);
    wdf2 = zeros(maxk+1,1);

  ii = T1;
  Do while ii <= T2;

     if imodel == 1;  iij = ii + 2;  endif;
     if imodel == 2;  iij = ii + 3;  endif;

     Do while iij <= T2;

      output off;
      format /rd/m1 5,0;
      Estconf = ii | iij;
      locate 5,1; Estconf';

     lmtk = zeros(maxk+1,1);
     lmlm = zeros(maxk+1,1);
     lmestd = zeros(maxk+1,1);
     lmtd = zeros(maxk+1,1);

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


     If imodel == 1;
        {wtstat1, wed, wtd0, wres, wbeta, wtval1, sig}
          = p_alt(y, 1, EstConf,kk);
             @ z(t) = [y(t-1), lags, 1, t, B(t), D(t)] @
        {wtstat2, wed, wtd0, wres, wbeta, wtval2, sig}
          = p_alt1(y, 1, EstConf,kk);
             @ z(t) = [y(t-1), lags, 1, t, D(t)] @
     Endif;

     If imodel == 2;
        {wtstat1, wed, wtd0, wres, wbeta, wtval1,sig}
          = p_altC(y, 1, EstConf,kk);
             @ z(t) = [y(t-1), lags, 1, t, B(t), D(t), DT(t)] @
        {wtstat2, wed, wtd0, wres, wbeta, wtval2,sig}
          = p_altC1(y, 1, EstConf,kk);
             @ z(t) = [y(t-1), lags, 1, t, D(t), DT(t)] @
     Endif;

     wdf1[kk+1] = wtstat1;
     if kk > 0; wtk1[kk] = wtval1[kk+1];  endif;
          @ t-stat of the coeff. of the last augmented term @

     wdf2[kk+1] = wtstat2;
     if kk > 0; wtk2[kk] = wtval2[kk+1];  endif;
          @ t-stat of the coeff. of the last augmented term @

    kk = kk + 1;
    endo;

    @ LP test (I) @
      jj = maxk;  isw = 0;
      do while isw ne 1;
        if (abs(wtk1[jj]) > crit_t) or (jj == 0);
           woptk1 = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

      lp1_k = wdf1[woptk1+1];
      if lp1_k < minlp1;
         minlp1 = lp1_k;
         lp1tb = Estconf;
         lp1k = woptk1;
         locate 7,10;
         "LP with B(t) ";; Estconf';; lp1k;
      endif;

   @ LP test (II) @
      jj = maxk;  isw = 0;
      do while isw ne 1;
        if (abs(wtk2[jj]) > crit_t) or (jj == 0);
           woptk2 = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

      lp2_k = wdf2[woptk2+1];
      if lp2_k < minlp2;
         minlp2 = lp2_k;
         lp2tb = Estconf;
         lp2k = woptk2;
         locate 8,10;
         "LP w/o B(t)  ";; Estconf';; lp2k;
      endif;



    iij = iij + 1;
  Endo;

  ii = ii + 1;
 Endo;

    output on;

    optb = lmtb; optk = lmk;

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

    "Model  (1=A, 2=C)      = " imodel;  "";

    "1. LP with B(t) "; "";
     optb = lp1tb; optk = lp1k;
     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;

     "LP Min. test with B(t) = "  wts;
     "Estimated break point  = "  optb';
     "Selected lag           = "  optk;
     "";
     "Est. coeff. of dummy   = "  wd';
     "Its t-stat             = "  wtd';

     "Standard error ..      = "  sig;
     "Standardized dummy     = "  wtd'./sig;

     if imodel == 1;
      ""; "Coeff and t-stat";
      "Z(t) = [y(t-1),(lags..omit), 1,t,B1(t),B2(t),D1(t),D2(t)] ";
     endif;
     if imodel == 2;
      ""; "Coeff and t-stat";
      "Z(t) = [y(t-1),(lags..omit), 1,t,B1(t),B2(t),D1(t),D2(t),DT1,DT2] ";
     endif;

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


    "2. LP without B(t) "; "";
     optb = lp2tb; optk = lp2k;
     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;

     "LP Min. test without B(t) = "  wts;
     "Estimated break point  = "  optb';
     "Selected lag           = "  optk;
     "";
     "Est. coeff. of dummy   = "  wd';
     "Its t-stat             = "  wtd';

     "Standard error ..      = "  sig;
     "Standardized dummy     = "  wtd'./sig;

     if imodel == 1;
      ""; "Coeff and t-stat";
      "Z(t) = [y(t-1),(lags..omit), 1,t,D1(t),D2(t)] ";
     endif;
     if imodel == 2;
      ""; "Coeff and t-stat";
      "Z(t) = [y(t-1),(lags..omit), 1,t,D1(t),D2(t),DT1,DT2] ";
     endif;

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


imodel = imodel + 1;
endo;

/*jjj = jjj + 1;
endo;*/


/*
Format: {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMCk(y, EstConf,kk,0);

Input:  y  = series
        Estconf = position of the structural break
        kk = # of augmentation lags
        q = # of truncation lags (ignore here)

Output: lmtstat  = ADF t-statistic
        lmd   = estimated coefficients of break dummies
        lmtd0  = t-stat for coefficients of break dummies
        lmres = residuals
        lmbeta   = all estimated coefficients
        lmtval   = t-stat of estimated coefficients
        lmsig  = 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);
     tval    = b ./ sqrt(diag(var_cov));

     if R == 0; estd = 0; td = 0; endif;
     if R > 0;
        estd = b[rows(b)-R+1:rows(b)];
         td  = b[rows(b)-R+1:rows(b)] ./
               diag(sqrt(var_cov[rows(b)-R+1:rows(b),rows(b)-R+1:rows(b)]));
     endif;

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) ;
     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;

     @ without B(t) @

     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 @

     R = 0;
     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);
     tval    = b ./ sqrt(diag(var_cov));

     if R == 0; estd = 0; td = 0; endif;
     if R > 0;
        estd = b[rows(b)-R+1:rows(b)];
         td  = b[rows(b)-R+1:rows(b)] ./
               diag(sqrt(var_cov[rows(b)-R+1:rows(b),rows(b)-R+1:rows(b)]));
     endif;

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



proc (6) = p_nullC(x,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 (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       = ones(rows(dep),1);

     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 = 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);
     tval    = b ./ sqrt(diag(var_cov));

     if R == 0; estd = 0; td = 0; endif;
     if R > 0;
        estd = b[rows(b)-R+1:rows(b)];
         td  = b[rows(b)-R+1:rows(b)] ./
               diag(sqrt(var_cov[rows(b)-R+1:rows(b),rows(b)-R+1:rows(b)]));
     endif;

retp(b, tval, estd, td, res, 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);
     tval    = b ./ sqrt(diag(var_cov));

     if R == 0; estd = 0; td = 0; endif;
     if R > 0;
        estd = b[rows(b)-R+1:rows(b)];
         td  = b[rows(b)-R+1:rows(b)] ./
               diag(sqrt(var_cov[rows(b)-R+1:rows(b),rows(b)-R+1:rows(b)]));
     endif;

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) ;
     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;

     @ without B(t) @

     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 @~ 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);
     tval    = b ./ sqrt(diag(var_cov));

     if R == 0; estd = 0; td = 0; endif;
     if R > 0;
        estd = b[rows(b)-R+1:rows(b)];
         td  = b[rows(b)-R+1:rows(b)] ./
               diag(sqrt(var_cov[rows(b)-R+1:rows(b),rows(b)-R+1:rows(b)]));
     endif;

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
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 ptrend(p,nobs) ;
local data, t, u, m, timep, it, iit, xmat, invx, beta, resid ;
     u        = ones(nobs,1) ;
     if p > 0 ;
       timep    = zeros(nobs,p) ;
       t        = seqa(1,1,nobs);
       m        = 1 ;
       do while m <= p ;
          timep[.,m] = t^m ;
          m = m + 1 ;
       endo ;
       xmat     = u~timep ;
     else ;
       xmat = u ;
     endif ;
retp(xmat) ;
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;