@  Written/Modified by: Junsoo Lee

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

/*  Amsler and Lee (1995, Econometric Theory) LM unit root test

program: AL.txt

*/

/*  Options to choose */

 output file = AL.out reset;
 "Amsler and Lee (Econometric Theory, 1995) LM t-stat ";

 maxk = 8;
     @ # maximum augmentation lags allowed for @

 n = 100;
     @ # of observations @

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

 pbreak = 50;
     @ break point between 0 1 and n-1 @

 crit_t = 1.645;

/* End of Options */

     lmtk = zeros(maxk+1,1);
     lmlm = zeros(maxk+1,1);

     kk = 0;  @ lag @
     do while kk <= maxk;

     Estconf = pbreak;

       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMk(y, EstConf,kk,0);

     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 = jj;  isw = 1;
        endif;
        jj = jj - 1;
    endo;

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

    output on;

    format /rd/m1 10,4;
    ""; "-----------------------------------------"; "";

       {lmtstat, lmed, lmtd0, lmres, lmbeta, lmtval, lmsig} =
           nLMk(y, EstConf,lmoptk1, 0);

     "LM test statistic      = "  lmtstat;
     "Selected lag           = "  lmoptk1;
     "Given break point      = "  pbreak;
     "";
     "";
     "Estimated coeff. of dummy var. = "  lmed;
	 "Its t-stat                     = "  lmtd0;

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


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

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;