clear?all;
close?all;
N=1000;
sigma=0.0332;
a1=1.6;
a2=-1.46;
a3=0.616;
a4=-0.1525;
M=100;
wsum1=zeros（4N-4）;
w2sum1=zeros（4N-4）;
asum1=zeros（1N）;
a2sum1=zeros（1N）;
for?loop=1:M
v=randn（1N）*sqrt（sigma）;
u=zeros（1N）;
for?i=1:（N-4）
????u（i+4）=a1*u（i+3）+a2*u（i+2）+a3*u（i+1）+a4*u（i）+v（i+4）;
end???
u=u（5:N）;
x=0.98;
P2=eye（4）/0.004;
k2=zeros（41）;
w2=zeros（4N-4）;
aa2=zeros（1N）;
w=zeros（4N-4）;
K=zeros（4N）;
p=zeros（44N）;
aa=zeros（1N）;
p（::1）=eye（4）;
u1=zeros（4N）;
a=zeros（1N）;
u1（:1）=[u（5）?0?0?0]‘;
u1（:2）=[u（6）?u（5）?0?0]‘;
u1（:3）=[u（7）?u（6）?u（5）?0]‘;
for?i=4:（N-8）
????u1（:i）=[u（i+4）?u（i+3）?u（i+2）?u（i+1）]‘;
end
for?j=2:（N-6）
????w（:j）=w（:j-1）;
????aa（j）=u（j）-（u1（:j）.‘）*w（:j）;
????p（::j）=p（::j-1）;
????a（j）=（u1（:j）.‘）*?p（::j）*u1（:j）+0.005;
????K（:j）=?