[学科前沿] gauss 编程应用于学科前沿二：非线性单位根 [推广有奖]

11楼

xuelida

发表于 2013-7-12 11:49:28

下面看看几个其他的检验方法Sup test，Exp test，Ave test，这三大检验很多人知道，相关理论可以参考相关的论文，如，Andrews, D. W. K. (1993). "Tests for parameter instability and structural change with unknown change point." Econometrica, 61: 821-856.
该方法也可以适用于未知突变情形，程序见下面。

12楼

xuelida

发表于 2013-7-13 15:55:27

下面的程序需要大家把论文看一下。
new;
library pgraph;
graphset;
load simul[200,14] =  simul.txt;
timex = seqa(1926,0.25,198);

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

/* Andrews (1993), "Tests for Parameter Instability and Structural Change with Unknown Change Point",  Econometrica, 61(4), p. 821-856.                            */
/* ---------  LM Version  -------------------*/
proc (2) = AndrewsLM(y, piezero, bw);
local T, Pie1MPie, LRvar, S, supLM, andrwsV, end1, eMu, H;

   T = rows(y);
   supLM = 0;                /* 给 supLM赋值为0，为了存储后面的计算结果 */
   andrwsV = zeros(T,1);    /* 给 andrwsV赋值为zeros(T,1)，为了存储后面的计算结果，并寻找“最大值”计算supLM */
   end1 = floor(T*piezero); /* floor求整，读者可以参考help，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  -----------------*/
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;
retp(andrws, andrwsV);
endp;

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