@  Written/Modified by: Junsoo Lee  (updated March 4, 2005)

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

/*  Lee and Strazicich Minimum LM unit root test with two structural breaks.


Reference paper: Lee, Junsoo, and Mark C. Strazicich, 2003,
"Minimum Lagrange Multiplier Unit Root Test with Two Structural Breaks."
The Review of Economics and Statistics, 85(4): 1082-1089.

For information on the similar test with one break, see the reference paper:
Lee, Junsoo, and Mark C. Strazicich, 2013, 
“Minimum LM Unit Root Test with One Structural Break.”
Economics Bulletin, 33(4): 2483-2492.

program: LS.txt

Notes:

To correct for serial correlations by including k first differenced lagged (augmented) terms,
the program first determines the optimal lag length (k) at each combination of two
breaks.  The optimal lag length is determined by a general to specific procedure.
Starting with the maximum number of lags, maxk, the t-statistic on the maximum lagged term
is examined to see if significant at the asymptotic 10% level.  If not, the maximum lagged term 
is dropped and the test is repeated until the maximum lagged term is significant, or no lags 
are found (see, for example, Ng and Perron, 1995).  Once the optimal lag length at each combination
of two breaks is determined, the program searches for the two break points where
the unit root test statistic is minimized (i.e., the most negative).

Model 1 includes two changes in intercept or level of the time series.

Model 2 includes two changes in intercept and trend slope.

To run the program below, name the output file, choose the maximum lag length k 
(8 in the example below), enter the number of observations in each time series (59 in the example below),
enter the number of time series tested (3 in the example below), enter the name of the
data source file, and give each time series a name.
 
*/

/*  Options to choose */

 output file = c:\Minorities\LS_two_break_high_school_grad_rate_w_max_8_lags.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 @

 /* Not needed:
 n = 50;
     @ # of observations @

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

    load data[59,3] = c:\Minorities\High_School_Grad_Rate_1957_2015.txt;
    
    let name = LNBWHSGRAD LNBMWHSGRAD LNBFWHSGRAD;
    
    jjj = 1;
    do while jjj <= cols(data);

    y = data[.,jjj];

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

     n = rows(y);
     format /m1/rd 5,0;
    "# of obs = " n;
    
    format /m1/rd 10,4;
    " First 3 observations: "  y[1:3,1]';
    " Last 3 observations:  "  y[n-2:n,1]';

/* End of Options */

 @ imodel  1 = Crash model (dummy on the intercepts, i.e., changing intercepts)
           2 = Breaking trend (dummy on the trends, i.e., changing intercepts and trends) @

 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;