@  Written/Modified by: Junsoo Lee

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

/* program: Panel_tw.txt

  Panel LM test using
  Lee and Strazicich Min LM test with two breaks.

  See also, LStwo.txt, panel_sp.txt (no break), panel_LS.txt (one break)

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

*/

@ ================================================================ @
     out2 = "panel_tw.out";

 "SP t-stat with the data-dependent choice of # lags ";
 "Also, panel test ";

 load crit[27,19] = c:\work\gauss\panelcri.dat;
@ ================================================================ @

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

 npanel = 19;
      @ # of panels @

    output file = ^out2 reset;
    format /rd/m1 4,0;
     "Date                     : ";; tt=date; tt[1:3]'; "";
     "max p                    : " maxk;
     "# panels                 : " npanel;
     format /rd/m1 4,1;

   @ set up @

     panel_t = zeros(npanel,1);
     panel_k = zeros(npanel,1);
     panel_b = zeros(npanel,2);

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

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

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

/* End of Options */

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

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

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

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

@ Reading critical values @

    @ row location @
        if T < 10;  row_n = 1; endif;
        if T > 10 and T < 25;  row_n = T-10+1; endif;
        if T > 24 and T < 28;  row_n = 16; endif;
        if T > 27 and T < 33;  row_n = 17; endif;
        if T > 32 and T < 38;  row_n = 18; endif;
        if T > 37 and T < 43;  row_n = 19; endif;
        if T > 42 and T < 48;  row_n = 20; endif;
        if T > 47 and T < 53;  row_n = 21; endif;
        if T > 54 and T < 65;  row_n = 22; endif;
        if T > 64 and T < 75;  row_n = 23; endif;
        if T > 74 and T < 85;  row_n = 24; endif;
        if T > 84 and T < 95;  row_n = 25; endif;
        if T > 94 and T < 150; row_n = 26; endif;
        if T > 149;            row_n = 27; endif;

@ Calculating the panel LM test statistic @

""; "";
"--------------------------------";
" PANEL LM statistic and summary ";
"--------------------------------";

output file = ^out2 on;

  format /rd/m1 7,3;

  "";"";

  LM_t = panel_t;
  meanLM = meanc(LM_t);

    meanp = 0; varp = 0;

    mmm = 1;
    do while mmm <= npanel;

        @ column location @
        if panel_k[mmm] > 8; optk = 8;
        else; optk = panel_k[mmm];  endif;

        ct1 = crit[row_n,2+2*optk];
        ct2 = crit[row_n,3+2*optk];
        ct = ct1 ~ ct2 ; ct';

        meanp = meanp + crit[row_n,2+2*optk];
        varp  = varp + crit[row_n,3+2*optk];

        mmm = mmm + 1;
    endo;

    meanp = meanp / npanel;
    varp  = varp / npanel;

    panel_LM = sqrt(npanel) * (meanLM - meanp) / sqrt(varp);
    "Panel LM test statistic = "  panel_LM;  "";

    outt = seqa(1,1,npanel)~panel_t~panel_k~panel_b;
    "Individual (i) LM statistic (ii) opt p and (iii) opt break point"; outt;

    ""; t1 = seqa(1,1,19); t1'; $name[2:20]';


/*
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;