跪求Matlab高手，在菲利普斯-配荣法做单位根检验时遇到问题。-经管之家官网！

小弟在看陆懋祖的《高等时间序列计量经济学》时，在做单位根检验数据模拟时，想把书后面附录的p336页的临界值表，实验一下。现在的情况是DF的情况一、二、四都已实现，但在做pp检验法时，表二情况四的临界值也没问题是 -3.3559，样本量是200.（可以查看p337页），但是在备择假设下mu由0变为0.5做检验时，结构很奇怪利用数据模拟了5000次，但是落入拒绝域的只有255，而非5000.现在将程序拷在下文，希望哪位大哥帮忙解决，谢谢诶，纠结很久了。
%PP检验法的的二种情况。
w=zeros(1,5000);
ru=1;
arfa=1;
darta=0;
for t=1:5000
e=randn(1,200);
u(1)=e(1);
y(1)=u(1);
for i=2:200
u(i)=0.01*u(i-1)+e(i);
y(i)=arfa+ru*y(i-1)+darta*i+u(i);
end
%对于大yipusilong的表示
for i=1:200
E(i)=sum(e(1,1:i));
end
%参数估计向量
%A=[arfahat,ruhat-1=ruhat001,dartahat-arfa];
B=[200 sum(E(1,1:199)) sum(1:200);
sum(E(1,1:199)) sum(E(1,1:199).^2) sum((2:200).*E(1,1:199));
sum(1:200) sum((2:200).*E(1,1:199)) sum((1:200).^2)];
C=[sum(e(1:200)),sum(E(1:199).*e(2:200)),sum((1:200).*e(1:200))];
A=inv(B)*C';
R=[200.^(1/2) 0 0;
0 200 0;
0 0 200.^(3/2)];
D=R*A;
yitasquare=[0 1 0]*inv(B)*[0 1 0]';
ruhat001=D(2,1)/200;
w(t)=ruhat001/(yitasquare)^(1/2);
end
a=sort(w);
b=a(fix(5000*0.05))
%在备择假设下的检验
z=zeros(1,5000);
ru=0.2;
arfa=0;
darta=0;
for t=1:5000
e=randn(1,200);
u(1)=e(1);
y(1)=u(1);
for i=2:200
u(i)=0.01*u(i-1)+e(i);
y(i)=arfa+ru*y(i-1)+darta*i+u(i);
end
%对于大yipusilong的表示
for i=1:200
E(i)=sum(e(1,1:i));
end
%参数估计向量
%A=[arfahat,ruhat-1=ruhat001,dartahat-arfa];
B=[200 sum(E(1,1:199)) sum(1:200);
sum(E(1,1:199)) sum(E(1,1:199).^2) sum((2:200).*E(1,1:199));
sum(1:200) sum((2:200).*E(1,1:199)) sum((1:200).^2)];
C=[sum(e(1:200)),sum(E(1:199).*e(2:200)),sum((1:200).*e(1:200))];
A=inv(B)*C';
R=[200.^(1/2) 0 0;
0 200 0;
0 0 200.^(3/2)];
D=R*A;
yitasquare=[0 1 0]*inv(B)*[0 1 0]';
ruhat001=D(2,1)/200;
z(t)=ruhat001/(yitasquare)^(1/2);
end
s=z<b;%等于1全部拒绝
v=sum(s)
运行结果：
b =
-3.3559
v =
255