new;
library pgraph;
graphset;

load simul[200,14] = c:\docs\bilkent\gradmetr\511\strucchg\simul.txt;

timex = seqa(1926,0.25,198);


/* ---------  The Andrews Tests for a Structural Break  ----------------*/

/*      Andrews (1993), "Tests for Parameter Instability and Structural 
           Change with Unknown Change Point", 
           Econometrica, 61(4), p. 821-856.                             */


/* ---------  LM Version  -------------------
**
** { supLM, andrwsV } = AndrewsLM( y, piezero, bw );
*/

proc (2) = AndrewsLM( y, piezero, bw );
   local T, Pie1MPie, LRvar, S, supLM, andrwsV, end1, eMu, H;

   /* Initialize variables */
   T = rows(y);
   supLM = 0;  andrwsV = zeros(T,1);
   end1 = floor(T * piezero);
   Pie1MPie = ( seqa(1,1,T)/T ) .* ( 1 - seqa(1,1,T)/T );
   eMu = y - meanc(y);
   LRvar = s2l( eMu, bw );
   S = cumsumc(eMu);
   H = (1/T)*(S)^2 ./ LRvar;

   /* ----  Calculate the Tests  --------- */
   andrwsV[end1:T-end1] = (1/Pie1MPie[end1:T-end1]) .* H[end1:T-end1];
   supLM  = maxc(andrwsV);
   retp( supLM, andrwsV );
endp;



/* -----------  Wald Version  -----------------
**
** { andrws, andrwsV } = AndrewsW( y, piezero, bw );
*/

proc (2) = AndrewsW( y, piezero, bw );
  local T, i, end1, n2, a1, a2, andrws1, andrws, andrwsV;
  local V1, V2;
  T = rows(y);
  andrws = 0;  andrwsV = zeros(T,1);
  end1 = floor(T * piezero);
  i=end1-1; do while i<(T-end1); i=i+1;
    n2=i;
    a1=meanc(y[1:n2]); a2=meanc(y[n2+1:T]);
    if bw==0;
      V1=((n2-1)/n2)*stdc(y[1:n2])^2;  V2=((T-n2-1)/(T-n2))*stdc(y[n2+1:T])^2;
    else;
      V1 = Vhat( bw, y[1:n2]-a1);  V2 = Vhat( bw, y[n2+1:T]-a2);
    endif;
    andrws1 = (a1-a2)^2/( (V1/i) + (V2/(T-i)) );
    if andrws1 > andrws;
      andrws=andrws1;
    endif;
    andrwsV[n2]=andrws1;
  endo; @ i loop @
  retp( andrws, andrwsV );
endp;

/* ---------  Other procs  ------------- */



/*
**   s2l
*/
proc (1) = s2l( e, l );
   local w, s, i, li, s2li, T;
   T = rows(e);
   s2li = 0;
   i=0; do while i<rows(l); i=i+1;
      li = floor( l[i]*(T/100)^0.25 );
      w=0; s=0; do while (s<li); s=s+1;
        w = w + ( 1 - s/(li+1) ) * e[s+1:T]'e[1:T-s];
      endo;
      s2li = (s2li)|( (e'e)/T + (2/T)* w );
   endo;
   retp( trimr(s2li,1,0) );
endp;

proc (1) = Vhat( bw, e );
  local S, V;
  S = Shat( bw, e );
  V = S;
  retp(V);
endp;

proc (1) = Shat( bw, e );
  local n, q, v, S0, Sv, Gv;
  n = rows(e);
  q = floor( bw*(n/100)^0.25 );
  S0 = (1/n)*e[1:n]'e[1:n];
  Sv=0; v=0; do while (v<q); v=v+1;
    Gv = (1/n)*e[v+1:n]'e[1:n-v];
    Sv = Sv + (1 - v/(q+1))*2*Gv;
  endo;
  retp( S0 + Sv );
endp;

@diff = simul[3:201,5]-simul[2:200,5];@
diff = simul[3:200,6];
{supLM,AndrewsV} = AndrewsLM(diff,0.10,1);
xy(timex,AndrewsV);supLM;