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就是模型的系数。