@  Written/Modified by: Junsoo Lee

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

/* LM test using one break.

  See also, SP.g,  LS_two.g, panel_sp.g (no break), panel_one.g (one break), panel_tw.g (two breaks)

  Lee and Strazicich (2004) Minimum LM unit root test with one structural break.

Reference paper: Lee, Junsoo, and Mark C. Strazicich, 2004, Manuscript,
Department of Economics, Appalachian State University.  

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_ONE.g

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 t-test statistic is minimized.

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, you must name the output file, choose the maximum lag length k 
(8 in the example below), enter the number of observations in each time series (51 in the example below),
enter the number of time series tested (43 in the example below), enter the name of the
data source file, and give each time series a name.
 
*/

@ ================================================================ @
  output file = "LS_one_gdp.out" reset;

 "SP-LM t-stat with the data-dependent choice of # lags ";
@ ================================================================ @

/* Options  */

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

     ""; "";
     
    format /rd/m1 4,0;
     "Date                     : ";; tt=date; tt[1:3]'; "";
     "max p                    : " maxk;
     format /rd/m1 4,1;

@ ==== DATA ====================================================== @

/*
   load data[45,20] = c:\work\hyster\unemp2.txt;
     "Data: Unemp2.txt";

     let name =
     Year Australia Austria Belgium Canada Denmark Finland
     France Germany Ireland Italy Japan Netherlands Norway
     Newzealand Spain Sweden Switzerland UK USA;

     jjj = 2;
     do while jjj <= 20;

     y = data[.,jjj];

     @ Due to missing observations @
     if jjj == 16;  y = y[2:rows(data)];  endif;
     if jjj == 20;  y = y[2:rows(data)];  endif;
     if jjj == 11;  y = y[1:rows(data)-1];  endif;
     if jjj == 14;  y = y[1:rows(data)-1];  endif;
     if jjj == 18;  y = y[1:rows(data)-1];  endif;
     if jjj == 19;  y = y[2:rows(data)-1];  endif;

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

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

*/

@ ================================================================ @

   load data[38,21] = gdp.txt;
   npanel = 20;
      @ # of panels @

     let name =
        Year  Australia  Austria  Belgium  Canada  Denmark  Finland  France
        Germany  Greece  Iceland  Ireland  Italy  Japan  Netherlands  Norway
        Spain  Sweden  Switzerland  UK  USA;

    data = ln(data);

 Timedum = 0;
 Do while Timedum <= 1;

     "";"";"";
     "=======================================";
     " Time Dummy (1=yes, 0=no) : " timedum;
     "=======================================";

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

    y = data[.,jjj];

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

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

     if cols(data) == npanel;  mmm = jjj;           fstcol = 0;
     elseif cols(data) == npanel+1;  mmm = jjj-1;   fstcol = 1;
     else; "Error in # of panels !!!";  stop;  endif;

     format /m1/rd 6,4;
     "Panel # = ";; mmm;  "";
     
        @ To control cross-country time effect @
     
         if timedum == 1;
             if fstcol == 0;
                 newdata = data - meanc(data');
             else;
                 tt1 = meanc(data[.,2:cols(data)]');
                 newdata = data[.,2:cols(data)] - tt1;
             endif;
             y = newdata[.,mmm];
         endif;
    
@ ================================================================ @


/* End of Options */

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

 imodel=1;
 do while imodel <= 2;
 @ once for each model for the panel test @

    cls;
    "---------------------------------------------------";
    "Lee and Strazicich Min LM t-stat with ONE break ";
    "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;

      output off;
      format /rd/m1 5,0;
      Estconf = ii;
      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 break = ";; Estconf;; "   lag = ";; lmk;
       endif;

  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;

      panel_t[mmm] = lmtstat;
      panel_k[mmm] = optk;
      panel_b[mmm] = optb;

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

Timedum = Timedum + 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;