%AR模型參數(shù)估計
clear
p=input('請輸入AR模型的階數(shù)p=')
for i=1:p
a(i)=input('請輸入AR模型的參數(shù)a=')
end
Q=input('請輸入方差q=')
L=input('請輸入樣本個數(shù)L=')

%產(chǎn)生L個w(n)，S(n)
w=randn(1,L);
b=sqrt(Q);
w=w*b;
power_w=var(w)
figure;
plot(w);
title('w(n)');

sum=zeros(1,L);
s=zeros(1,L);
s=w;
for n=1+p:L
    for i=1:p
        s(n)=w(n)-a(i)*s(n-i);
    end
end
figure;
plot(s)
title('s(n)');

%估計p+1個自相關函數(shù)，得a1....ap
fss=zeros(1,p+1);
sum=zeros(1,p+1);
for m=1:p+1
    for i=1:L-abs(m)+1
        sum(m)=sum(m)+s(i)*s(i+abs(m-1));
        d=1/(L-abs(m));
        fss(m)=d*sum(m);
    end
end
U=zeros(1,p+1);
U(1,1)=1;
for i=1:p+1
    for j=1:p+1
        Rss(i,j)=fss(abs(i-j)+1);
    end
end


B=inv(Rss)*U';
A=B/B(1);
a1=A(2)
Q1=B(1)

aex=mean(a1-a(1))^2
qex=mean(Q1-Q)^2