@  Written/Modified by: Junsoo Lee

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

/*  Lee and Strazicich (1999, WP) minimum LM unit root test with two breaks

program: LStwo.txt

  Lee and Strazicich Min LM test with two breaks.

  The program for the case with one break is given as LS.txt, and this program
  program provides the procedure for the case of 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 = LStwo.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;
    "---------------------------------------------------";
    "Lee and Strazicich Min LM 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 @
        minlm = 1000;

    lmestd1 = zeros(T,1);
    lmtd1 = zeros(T,1);     lmlm1 = zeros(T,1);
    lmtk1 = zeros(T,1);     lmoptk1 = zeros(T,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;
       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMk(y, EstConf,kk,0);
     Endif;

     If imodel == 2;
       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMCk(y, EstConf,kk,0);
     Endif;

     lmlm[kk+1] = lmtstat;
     if kk > 0; lmtk[kk] = lmtval[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(lmtk[jj]) > crit_t) or (jj == 0);
           lmoptk = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

      tt = lmoptk+1;
      lm_k = lmlm[tt];
      if lm_k < minlm;
         minlm = lm_k;
         lmtb = Estconf;
         lmk = lmoptk;
         format /rd/m1 5,0;
         locate 9,10;
         "Updated two breaks = ";; Estconf';; "   lag = ";; lmk;
       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;

    If imodel == 1;
       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMk(y, optb, optk, 0);
    Endif;

    If imodel == 2;
       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMCk(y, optb, optk, 0);
    Endif;

     "";
     "Min. test statistic    = "  lmtstat;
     "Estimated break point  = "  optb';
     "Selected lag           = "  optk;
     "";
     "Est. coeff. of dummy   = "  lmed';
     "Its t-stat             = "  lmtd0';

     "Standard error ..      = "  lmsig;
     "Standardized dummy     = "  lmed'./lmsig;

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

     tt1 = lmbeta[1] | lmbeta[2+optk:rows(lmbeta)];
     tt2 = lmtval[1] | lmtval[2+optk:rows(lmtval)];
     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) = nLMk(y, config, lag, q);
     local lr, res, lhat, alpha, za, zt, nobs, dy, rho;
     local sse, tau, ratios, i, n, z, dz, pos, st, du, s1, ind, tt, temp, sai;
     local RR, tdel, whvar, nwvar, temp2, dt, estd, td, ch, k;

     nobs = rows(y);

     if (config == 0);
         RR = 0;
         z = seqa(1,1,nobs);
         dz = ones(nobs-1,1);
     endif;

     if (config ne 0);
       RR = rows(config);
       pos = config[1];
       du = zeros(pos,1) | ones(nobs-pos,1);
       i=2;
       do while i <= RR;
          pos = config[i];
          temp = zeros(pos,1) | ones(nobs-pos,1);
           du = du ~ temp;
           i = i + 1;
       endo;
       z = seqa(1,1,nobs) ~ du;
       dz = diff(z,1);
     endif;

     dy = diff(y,1);
     tdel = dy / dz;
     sai = y[1] - z[1,.]*tdel;
     st = y - sai - z*tdel;
     s1 = st[1:rows(st)-1,1];

     ind = s1;
     ch  = diff(st,1);

     k       = 0     ;
     If (lag > 0) ;
        do until k >= lag ;
           k = k + 1 ;
           ind = ind ~ lagn(ch,k) ;
        endo ;
     Endif ;

     dy  = trimr(dy,k,0);
     ind = trimr(ind,k,0);
     dz  = trimr(dz,k,0);

     ind = ind ~ dz;

     alpha = dy / ind;
     rho = nobs * alpha[1];
     res    = dy - ind*alpha;

     sse    = res'res/(rows(res)-cols(ind));
     tt     = sqrt(sse*inv(ind'ind));
     tau    = alpha[1] ./ tt[1,1];

     @WhVar = WhiteEst(res,q);
     NWVar = NewWEST(res,q);@
     @LEE@

     @lr = lrvar(res,q);
     ratios = sse / lr;
     Za     = rho / ratios;
     Zt     = tau / sqrt(ratios);@

     if RR == 0; estd = 0; td = 0; endif;
     if RR > 0;
       estd = alpha[rows(alpha)-RR+1:rows(alpha)];
         td = estd ./
           diag(tt[rows(alpha)-RR+1:rows(alpha),rows(alpha)-RR+1:rows(alpha)]);
     endif;

     retp(tau, estd, td, res, alpha, alpha./diag(tt), sqrt(sse));
endp ;



proc (7) = nLMCk(y, config, lag, q);
     local lr, res, lhat, alpha, za, zt, nobs, dy, rho;
     local sse, tau, ratios, i, n, z, dz, pos, st, du, s1, ind, tt, temp, sai;
     local RR, tdel, whvar, nwvar, temp2, dt, estd, td, ch, k;

     nobs = rows(y);

     if (config == 0);
         RR = 0;
         z = seqa(1,1,nobs);
         dz = ones(nobs-1,1);
     endif;

     if (config ne 0);

       RR = rows(config);

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

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

       z = seqa(1,1,nobs) ~ du ~ dt;
       dz = diff(z,1);

     endif;

     dy = diff(y,1);
     tdel = dy / dz;
     sai = y[1] - z[1,.]*tdel;
     st = y - sai - z*tdel;
     s1 = st[1:rows(st)-1,1];

     ind = s1;
     ch      = diff(st,1);

     k       = 0     ;
     If (lag > 0) ;
        do until k >= lag ;
           k = k + 1 ;
           ind = ind ~ lagn(ch,k) ;
        endo ;
     Endif ;

     dy  = trimr(dy,k,0);
     ind = trimr(ind,k,0);
     dz  = trimr(dz,k,0);

     ind = ind ~ dz;

     alpha = dy / ind;
     rho = nobs * alpha[1];
     res    = dy - ind*alpha;

     sse    = res'res/(rows(res)-cols(ind));
     tt     = sqrt(sse*inv(ind'ind));
     tau    = alpha[1] ./ tt[1,1];

     @WhVar = WhiteEst(res,q);
     NWVar = NewWEST(res,q);@
     @LEE@

     @lr = lrvar(res,q);
     ratios = sse / lr;
     Za     = rho / ratios;
     Zt     = tau / sqrt(ratios);@

     if RR == 0; estd = 0; td = 0; endif;
     if RR > 0;
       estd = alpha[rows(alpha)-RR+1:rows(alpha)];
         td = estd ./
           diag(tt[rows(alpha)-RR+1:rows(alpha),rows(alpha)-RR+1:rows(alpha)]);
     endif;

     retp(tau, estd, td, res, alpha, alpha./diag(tt), sqrt(sse));
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;