@  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 for Panel LM test with breaks

  program: Panel_LS.prg

  Panel LM test using
  Lee and Strazicich Min LM test with one break.

  See also, LS.txt, panel_sp.txt (no break), panel_tw.txt (two breaks)

    Determining k and T(b) under the alternative for Min LM tests
    and compute the panel LM test statistic.

    (i) Four different methods are used to determine Tb.
    (ii) t-test is used to determine k.

*/

@ ================================================================ @
    new;
    out1 = "panel.out";
    out2 = "panelsum.out";

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


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

/* Options  */

    crit_t = 1.645;
        @ for lag selection @

    maxk = 8;
        @ max # of lags @

    End1 = 0.1;  End2 = 0.9;
        @ excluding end points for grid search for breaks @

    npanel = 19;
        @ # of panels @

    timedum = 0;
        @ 1 to control cross country time effect by adding time dummies @

   @ optons printing @

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

    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;
     "End points               : " End1  End2;


   @ set up @

     panel_t = zeros(4,npanel);
     panel_k = zeros(4,npanel);
     panel_b = zeros(4,npanel);


Timedum = 0;
Do while Timedum <= 1;

     output file = ^out1 on;
     "";"";"";
     "=======================================";
     " Time Dummy (1=yes, 0=no) : " timedum;
     "=======================================";

     output file = ^out2 on;
     "";"";"";
     "=======================================";
     " Time Dummy (1=yes, 0=no) : " timedum;
     "=======================================";

     output file = ^out1 on;


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

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

@ ===== End of reading data ====================================== @

     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 @

 @ Note: the following correction is made only for the unemp data due to
   missing observations @

   if timedum == 1;
        if fstcol == 0;
            newdata = data - meanc(data');
        else;
            tt1 = meanc(data[.,2:cols(data)]');
            tt1[1] = (sumc(data[2,.]')-data[2,16]-data[2,19]-data[2,20])/17;
            tt1[19] = (sumc(data[20,.]')-data[20,11]-data[20,14]
                      -data[20,18]-data[20,19])/16;
            newdata = data[.,2:cols(data)] - tt1;
        endif;
        y = newdata[.,mmm];
    endif;

    @ Set up @

    T = n;

    T1 = round(end1 * rows(y));
    T2 = round(end2 * rows(y));

    if T1 < maxk+2; T1 = maxk + 3;  endif;

    lmsbc1 = zeros(T,1);    lmestd1 = zeros(T,1);
    lmtd1 = zeros(T,1);     lmlm1 = zeros(T,1);
    lmtk1 = zeros(T,1);     lmoptk1 = zeros(T,1);


  @ ==== Loop ====@

  ii = T1;
  Do while ii <= T2;
     output off; "*";;

     Estconf = ii;

      lmtk = zeros(maxk+1,1);
      lmsbc = 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;

        {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
             nLMk(y, EstConf,kk,0);
         @ z(t) = [1, t, D(t)], Reg = [S(t-1), 1, B(t), lags] @

       Vars = (lmres'lmres)/T;
       lmsbc[kk+1] = T*LN(Vars) + (5 + kk).* LN(T);
       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;

    @ --- Choosing k for each of given break points --- @

      @  Select k when the last augmented coeff is significant
         starting from maxk to 0.  If all coeff are insig., k = 0.
         That is, general to specific approach @

      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;
      lmsbc1[ii] = lmsbc[tt];
      lmestd1[ii] = lmestd[tt];
      lmtd1[ii] = lmtd[tt];
      lmlm1[ii] = lmlm[tt];

   ii = ii + 1;
 Endo;
 @ === End of loop === @

   output on;


 @ =======  Now, for the testing model @

    lmo1 = minindc1(lmsbc1[T1:T2]) + T1 - 1;
    lmo2 = maxindc(abs(lmestd1[T1:T2])) + T1 - 1;
    lmo3 = maxindc(abs(lmtd1[T1:T2])) + T1 - 1;
    lmo4 = minindc1(lmlm1[T1:T2]) + T1 - 1;

    format /rd/m1 8,3;  "";

     "(1)    SBC      "; "";
      optb = lmo1; optk = lmoptk1[optb];
      {lmts,lmd,lmtd,lmres,lmb,lmt,lmsig} = nLMk(y,optb,optk,0);

      "  LM stat = " lmts;
      "  opt_k   = " optk;
      "  opt_b   = " optb; "";
      "  Est_b   = " lmb[1];
      "  Est_d   = " lmd;
      "  Est_td  = " lmtd;
      "  sig     = " lmsig;
      "  std sig = " lmd/lmsig;    "";

      panel_t[1,mmm] = lmts;
      panel_k[1,mmm] = optk;
      panel_b[1,mmm] = optb;


     "(2)    Max|d|   ";  "";
      optb = lmo2; optk = lmoptk1[optb];
      {lmts,lmd,lmtd,lmres,lmb,lmt,lmsig} = nLMk(y,optb,optk,0);

      "  LM stat = " lmts;
      "  opt_k   = " optk;
      "  opt_b   = " optb; "";
      "  Est_b   = " lmb[1];
      "  Est_d   = " lmd;
      "  Est_td  = " lmtd;
      "  sig     = " lmsig;
      "  std sig = " lmd/lmsig;   "";

      panel_t[2,mmm] = lmts;
      panel_k[2,mmm] = optk;
      panel_b[2,mmm] = optb;


     "(3)    Max|td|  ";  "";
      optb = lmo3; optk = lmoptk1[optb];
      {lmts,lmd,lmtd,lmres,lmb,lmt,lmsig} = nLMk(y,optb,optk,0);

      "  LM stat = " lmts;
      "  opt_k   = " optk;
      "  opt_b   = " optb; "";
      "  Est_b   = " lmb[1];
      "  Est_d   = " lmd;
      "  Est_td  = " lmtd;
      "  sig     = " lmsig;
      "  std sig = " lmd/lmsig; "";

      panel_t[3,mmm] = lmts;
      panel_k[3,mmm] = optk;
      panel_b[3,mmm] = optb;


     "(4)   Min t-stat";  "";
      optb = lmo4; optk = lmoptk1[optb];
      {lmts,lmd,lmtd,lmres,lmb,lmt,lmsig} = nLMk(y,optb,optk,0);

      "  LM stat = " lmts;
      "  opt_k   = " optk;
      "  opt_b   = " optb; "";
      "  Est_b   = " lmb[1];
      "  Est_d   = " lmd;
      "  Est_td  = " lmtd;
      "  sig     = " lmsig;
      "  std sig = " lmd/lmsig;

      panel_t[4,mmm] = lmts;
      panel_k[4,mmm] = optk;
      panel_b[4,mmm] = optb;

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;

ooo = 1;
do while ooo <= 4;

  format /rd/m1 7,3;

  "";"";
  if ooo == 1; "(1) Using SBC ";  ""; endif;
  if ooo == 2; "(2) Using max|d| ";  ""; endif;
  if ooo == 3; "(3) Using max|td| ";  ""; endif;
  if ooo == 4; "(4) Using Min LM test";  ""; endif;

  LM_t = panel_t[ooo,.]';
  meanLM = meanc(LM_t);

    meanp = 0; varp = 0;

    mmm = 1;
    do while mmm <= npanel;

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

        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[ooo,.]'~panel_k[ooo,.]'~panel_b[ooo,.]';
    "Individual (i) LM statistic (ii) opt p and (iii) opt break point"; outt;

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

  ooo = ooo + 1;
endo;

Timedum = Timedum + 1;
endo;


@ ====== Procedures for min Za and LM tests @

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 ;

/* *** Dec. 16, 1998 *** */

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 ;

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 (3) = LM(y, config, 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, 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);
       du = zeros(nobs,RR);

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

     retp(stat, za, res);
endp ;


proc (3) = LMC(y, config, 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, 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 ~ 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;

     retp(stat, za, res);
endp ;


proc (5) = nLM(y, config, 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, 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 ~ 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, za, res, estd, td);
endp ;


proc (5) = nLMC(y, config, 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, 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 ~ 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, za, res, estd, td);
endp ;


proc (6) = p_null(x,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;

     @dep     = trimr(x,1,0);@
     dep     = diff(x,1);
     ch      = diff(x,1);
     k       = 0     ;
     z       = ones(rows(dep),1);

     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));@

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2;

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

       z1 =  du;
       dz = diff(z1,1);

     Endif;

     @z1 = trimr(z1,1,0);@
     @z1 = trimr(z1,k,0);@
     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp(b, tval, estd, td, res, sqrt(so));
endp ;


proc (7) = p_alt(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;

     @ dep    = trimr(x,1,0);@
     dep     = diff(x,1);

     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;

     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z)) ;

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

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

       z1 =  du;
       dz = diff(z1,1);

     Endif;

    @ z1 = trimr(z1,1,0);
      z1 = z1[1:rows(z1)-k-1,.];@

     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     @z[1:10,.]';@

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;


proc (6) = p_nullC(x,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

    if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;

     @dep     = trimr(x,1,0);@
     dep     = diff(x,1);

     ch      = diff(x,1);
     k       = 0     ;
     z       = ones(rows(dep),1);

     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));@

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

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

       z1 =  du;
       dz = diff(z1,1);

     Endif;

     @z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);@
     z1 = z1[k+2:nobs,.];


     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp(b, tval, estd, td, res, sqrt(so));
endp ;



proc (7) = p_altC(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td, tval;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;

     @ dep    = trimr(x,1,0);@
     dep     = diff(x,1);

     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;@
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

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

       z1 =  du ~ dt;
       dz = diff(z1,1);

     Endif;

     @z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);@
     z1 = z1[k+2:nobs,.];

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     tval    = b ./ sqrt(diag(var_cov));
     estd    = b[rows(b)];
     td      = b[rows(b)] ./ diag(sqrt(var_cov[rows(b),rows(b)]));

retp((b[1,1])/sqrt(var_cov[1,1]), estd, td, res, b, tval, sqrt(so)) ;
endp ;


/*  ********** */

proc (6) = nperon2(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td;
     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;


     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ; @
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

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

       z1 =  du;
       dz = diff(z1,1);

     Endif;

     z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     estd    = b[cols(b)];
     td      = b[cols(b)] ./ diag(sqrt(var_cov[cols(b),cols(b)]));

retp(b[1,1],(b[1,1]-1)/sqrt(var_cov[1,1]), tvalk, estd, td, res) ;
endp ;


proc (6) = nperon2C(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f, tvalk, estd, td;

     local R, dz, du, pos, temp, kk, temp2, t2, dt, i;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;


     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     @if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ; @
     z = z ~ ones(rows(z),1) ~ seqa(1,1,rows(z));

  @ Here @

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
       @du = du ~ temp2;@
       du = temp2 ~ du;

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

       z1 =  du ~ dt;
       dz = diff(z1,1);

     Endif;

     z1 = trimr(z1,1,0);
     z1 = trimr(z1,k,0);

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z);
     tvalk   = b[2:2+k] ./ diag(sqrt(var_cov[2:2+k,2:2+k]));
     estd    = b[cols(b)];
     td      = b[cols(b)] ./ diag(sqrt(var_cov[cols(b),cols(b)]));

retp(b[1,1],(b[1,1]-1)/sqrt(var_cov[1,1]), tvalk, estd, td, res) ;
endp ;



/*
** Perron 2P
**
** Purpose: Compute the Perron's Augmented Dickey Fuller statistic, allowing
**          for deterministic polynomial time trends of an arbitrary
**          order.

   With B(t)
   Two step procedure is considered.

**
** Format:  {alpha, tstat} = perron2(x,p,config,l);
**
** Input:   x     -- time series variable
**
**          p     -- order of the time-polynomial to include in the
**                   ADF regression.  Set p = -1 for no deterministic
**                   part.
**
**          l     -- number of lagged changes of x to include in the
**                   fitted regression.
**
** Output:  alpha       --  estimate of the autoregressive parmaeter;
**
**          tstat       --  ADF t-statistic
**
*/

proc (2) = perron2P(x,p,config,l) ;
     local b,k,z,res,so,var_cov,xx,z1;
     local timep,t,m,xmat,nobs,dep,ch,f,i;

     local R, dz, du, pos, temp, kk, temp2, z0;

     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;

     nobs    = rows(x);

     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;


     dep = x;
     z = ones(rows(x),1) ~ seqa(1,1,rows(x));

     If (config ne 0);

       R = rows(config);

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

       temp2 = zeros(nobs,1);
       temp2[pos+1,1] = 1;
      @du = du ~ temp2;@

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

     Endif;

     z = z ~ du;

     b = dep / z;
     res = dep - z*b;

     x = res;


     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z      = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
           z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);

     z1 = trimr(temp2,1,0);
     z1 = trimr(z1,k,0);

     if (config ne 0); z = z ~ z1; endif;

     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z) ;

retp(b[1,1],(b[1,1]-1)/sqrt(var_cov[1,1]) ) ;
endp ;


/*
** ADF
**
** Purpose: Compute the Augmented Dickey Fuller statistic, allowing
**          for deterministic polynomial time trends of an arbitrary
**          order.
**
** Format:  {alpha, tstat} = adf(x,p,l);
**
** Input:   x     -- time series variable
**
**          p     -- order of the time-polynomial to include in the
**                   ADF regression.  Set p = -1 for no deterministic
**                   part.
**
**          l     -- number of lagged changes of x to include in the
**                   fitted regression.
**
** Output:  alpha       --  estimate of the autoregressive parmaeter;
**
**          tstat       --  ADF t-statistic
**
*/
proc (2) = adf(x,p,l) ;
     local b,k,z,res,so,var_cov,xx;
     local timep,t,m,xmat,nobs,dep,ch,f;
     if (p < -1);
        "Error: p cannot be set < -1";
         retp(0,0,zeros(6,1));
     endif ;
     if (cols(x) > 1);
        "Error: ADF cannot handle a data matrix; cols(x) > 1 (=" cols(x) ")";
         retp(0,0,zeros(6,1));
     endif ;
     nobs    = rows(x);
     if (nobs - (2*l) + 1 < 1) ;
        "Error: l is too large; negative degrees of freedom.";
         retp(0,0,zeros(6,1));
     endif ;
     dep     = trimr(x,1,0);
     ch      = diff(x,1);
     k       = 0     ;
     z       = trimr(lagn(x,1),1,0) ;
     If (l > 0) ;
        do until k >= l ;
           k = k + 1 ;
       z = z~lagn(ch,k) ;
        endo ;
     Endif ;
     z       = trimr(z,k,0);
     dep     = trimr(dep,k,0);
     if ( p > -1) ;
        z = z~ptrend(p,rows(z));
     endif ;
     b       = dep/z ;
     res     = dep - z*b ;
     so      = (res'res)/(rows(dep)-cols(z));
     var_cov = so*inv(z'z) ;

retp(b[1,1],(b[1,1]-1)/sqrt(var_cov[1,1]) ) ;
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;


proc (1) = ident(stat,tmax,tmin);
        @ program to count frequencies of each cell
          Use like this:
            ipos = ident(tau,tmax,tmin);
            itcnt(ipos) = itcnt(ipos) + 1;
          Refer to cdfgr.prg;
        @

        local temp, range, range1, ipos;
        ipos=0;
        temp = (stat - tmin)*100;
        if (temp < 0); temp=0;  endif;
        ipos =temp;
        range = (tmax-tmin) * 100;
        if (ipos > range);  ipos = range;  endif;
        ipos = ipos+1;
      retp(ipos);
endp;

     /*       subroutine ident(stat,tmax,tmin,itemp)
              implicit double precision (a-h,o-z)
              itemp=0
              temp = (stat - tmin)*100.0d0
              if(temp.lt.0) temp=0
              itemp=temp
              if(itemp.lt.0) itemp=0
              range=tmax-tmin
              irange=range
              itemax=irange*100
              if(itemp.gt.itemax) itemp=itemax
              itemp=itemp+1
            return
            end
      */