clear?all;
clc;

%從MA模型X（t）=a（t）-1.3a（t-1）+0.4a（t-2）中提取200個樣本點
K=200;
a1=randn（1K）;
x1=zeros（1K）;
x1（1）=1;
x1（2）=1;

for?i=3:K
????x1（i）=a1（i）-1.3*a1（i-1）+0.4*a1（i-2）;
end

X1=x1;

%從AR模型X（t）+0.8X（t-1）=a（t）中提取1000個樣本點
N=1000;
a2=randn（1N）;
x2=zeros（1N）;
x2（1）=0.1;
x2（2）=0.1;
x2（3）=1;

for?i=3:N
????x2（i）=1.5*x2（i-1）-0.5*x2（i-2）+a2（i）;
end

%取x2的前800個數據作為模型識別的樣本觀測值
L=800;
X2=x2（11:L）;

%p1=autocorr（X1L/10）;
%?autocorr（X1）;

%計算樣本自相關函數和樣本偏相關函數
r=zeros（1L）;???????????%樣本協方差矩陣
for?k=0:L-1?????????????%求樣本協方差
???rk=0;
???for?i=1:L-k
??????rk=rk+X2（i）*X2（i+k）;
???end
??r（k+1）=（1/L）*rk;
end

p=zeros（1L）;
for?i=1:L????????????%求樣本自相關函數
????p（i）=r（i）/r（1）;
end



%ppp=parcorr（X2L/10）;

?M=L/10;
pp=zeros（MM）;???????????%樣本偏相關系數矩陣
pp（11）=p（2）;
p（1）=[];??????

for?k=1:M-1??????????????%求樣本偏相關函數
????s1=0;
????s2=0;
????for?j=1:k
????????s1=s1+p（k+1-j）*pp（kj）;
????????s2=s2+p（j）*pp（kj）;
????end
????pp（k+1k+1）=（p（k+1）-s1）/（1-s2）;
????for?j=1:k
????pp（k+1j）=pp（kj）-pp（k+1k+1）*pp（kk+1-j）;
????end
end??

pp1=diag（pp）;????????%提取偏相關系數

subplot（211）
plot（p）;
title（‘自相關圖‘）;
subplot（212）
plot（pp1）;
title（‘偏相關圖‘）;


%識別模型，返回模型的階數
%取M=L/10個自相關和偏相關系數

count=0;
for?f1=1:M
????if?abs（p（f1））<=1/sqrt（L）
????????count=count+1;
????end
????if??count/M>=0.683
???????break;
????end
end


count_1=0;
for?f2=1:M
????if?abs（pp1（f2））<=1/sqrt（L）
????????count_1=count_1+1;
????end
????if?（count_1）/M>=0.683
????????break
????end
end


if?f1==M&&f2????disp（‘建立AR模型‘）;
????
????%確定模型階數
?????for?i=1:M???????????????????????
????????if??abs（pp1（i））<=1/sqrt（L）
????????????jieshu=i-1;
????????????break
????????end
????end
????fprintf（‘其階數為:%d\n‘jieshu）;
????
????%fai=ar（X2jieshu）;?????????????%AR模型求系數函數
????
????%最小二乘法求模型系數
????for?n=1:jieshu???????????????????
????????Y=X2‘;
????????Y（1:n）=[];
????????m=L-n;X=[];
????????for?i=1:m
????????????for?j=1:n
????????????????X（ij）=X2（n+i-j）;
????????????end
????????end
????????xishu=inv（X‘*X）*（X‘*Y）;
????end
????
????%用后200個數據檢驗所得模型是否正確
????X3=x2（1L+1:end）;???%提取后200個樣本點
????T=N-L;
????a3=zeros（1T）;????????%由觀測樣本點計算得到的白噪聲點
????ra=ones（1T+1）;??????%白噪聲的協方差
????
????a3（1）=X3（1）;
????a3（2）=X3（2）-xishu（1）*X3（1）;
????
???if?jieshu>1?????
????for?i=1:T