matlab之 时间序列MA（q）模型

理学院真的是一个伟大的学院，什么时间序列啊，这么难的东西，那么厚的一本书，我们竟然要把它学会，还要考试（说好的人与人之间的信任呢/(ㄒoㄒ)/~~）

已知 零均值平稳时间序列MA（q）序列的自协方差函数，求模型的系数。

function [k,cur_pii,sigma2,b]=time_serise(r)
%MA(2)
q=2;    #MA（q）
guest=100; % should be changed,最小迭代次数
eps=1e-8;
%r0=12.4168;
%r=[12.4168,-4.7520,5.2];%r1是原数据r0，自协方差函数
%r=[7.4084,-2.664,3.4];
r0=r(1);
for i=4:guest*10
    r(i)=0;
end
 
pre_pii=[];cur_pii=[];omg=[];gama=[];
%自协方差阵
 
for k=1:guest
    for i=1:k
        for j=1:k
            gama(i,j)=r(abs(i-j)+1);%公式，+1是因为r（1）实际是r0
        end
    end
    for i=1:q
        for j=i:(k+i-1)
            omg(i,j-i+1)=r(j+1);
        end
    end
    if k==1
        pre_pii=omg*inv(gama)*omg';
    else
        cur_pii=omg*inv(gama)*omg';
        sum=0;
        for i=1:q
            for j=1:q
                sum=sum+(pre_pii(i,j)-cur_pii(i,j))^2;
            end
        end
        if sum < eps
            k;
            cur_pii;
            break
        end  
        pre_pii=cur_pii;
    end
end
for i=q:q
    tmp=eye(i+1);
    A=tmp(2:i+1,1:i);
    C=zeros(i,1);
    C(1)=1;
    Tmp=C'*cur_pii*C;
    sigma2=r0-C'*cur_pii*C;
    Gama=[];
    
    for j=1:i
        Gama(j)=r(1+j);
    end
    b=[];
    b=(Gama'-A*cur_pii*C)/sigma2;
end
end
r=[7.4084,-2.664,3.4];#自协方差函数
[k,cur_pii,sigma2,b]=time_serise(r)%b是系数，sigma2是sigma平方。

向量b就是模型的系数。