@  Written/Modified by: Junsoo Lee  (updated February 8, 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 one structural break.


Reference paper is: Lee, Junsoo, and Mark C. Strazicich, 2013, 
“Minimum LM Unit Root Test with One Structural Break.”
Economics Bulletin, 33(4): 2483-2492.

See also, 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.

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 possible
break point.  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 break point 
is determined, the program searches for the break point where the unit root test statistic is
minimized (i.e., the most negative).

Model 1 includes one structural change in the intercept or level of the time series.

Model 2 includes one structural change in the 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 being 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_one_break_high_school_grad_rate_1957_2015_w_max_8_lags.out reset;
 
 maxk = 8;
     @ maximum number augmentation lags allowed @

/* Not needed:
 n = 100;
     @ # 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];
 @ y = ln(y);    ... if log is to be taken @
 
 "";"";
 "jjj = "  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]';


 T1_end = .1;   T2_end = .9;
     @ end points (10% of both ends are not considered in break search) @

 crit_t = 1.645;
     @ critical value of the t-stat to determine number of augmented terms @

/* End of Options */

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

 imodel=1;
 do while imodel <=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;

    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;
     output off; "*";;

     Estconf = ii;

     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;

     Estconf = ii;

     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;

     lmestd[kk+1] = lmed;
     lmtd[kk+1] = lmtd0;
     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);
           lmoptk1[ii] = jj;  isw = 1;
        endif;
        jj = jj - 1;
      endo;

      tt = lmoptk1[ii]+1;
      lmlm1[ii] = lmlm[tt];

   ii = ii + 1;
 endo;

    output on;

    lmo4 = minindc1(lmlm1[T1:T2]) + T1 - 1;
    optb = lmo4; optk = lmoptk1[optb];

    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;
     "";
     "Estimated coeff. of dummy var. = "  lmed;
     "Its t-stat                     = "  lmtd0;

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

     if imodel == 1;
      ""; "Coeff and t-stat";
      "Z(t) = [S(t-1), (lags..omitted), 1, B(t)] ";
     endif;
     if imodel == 2;
      ""; "Coeff and t-stat";
      "Z(t) = [S(t-1), (lags..omitted), 1, B(t), D(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, stat;

     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));
     tt     = sqrt(sse*inv(ind'ind));
     tau    = alpha[1] ./ tt[1,1];

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

     stat = tau;
   
   /*  if q > 0;
     lr = lrvar(res,q);
     ratios = sse / lr;
     Za     = rho / ratios;
     Zt     = tau / sqrt(ratios);
     stat = zt;
    endif;
   */
   
     estd = alpha[rows(alpha)];
     td = estd ./ tt[rows(alpha),rows(alpha)];

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


/* Model C */

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, stat;

     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));
     tt     = sqrt(sse*inv(ind'ind));
     tau    = alpha[1] ./ tt[1,1];

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

     stat = tau;

   /*  if q > 0;
     lr = lrvar(res,q);
     ratios = sse / lr;
     Za     = rho / ratios;
     Zt     = tau / sqrt(ratios);
     stat = zt;
    endif;
   */
   
     estd = alpha[rows(alpha)];
     td = estd ./ tt[rows(alpha),rows(alpha)];

     retp(stat, 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;