时间序列matlab代码

function [yucezhi]=shijianxulie(data,cishu);
[hang lie]=size(data);
yucezhi=zeros(hang,lie);
z=zeros(hang,cishu);
for nt=1:hang
x=data(nt,:);
xx=x;
x=[0 diff(x)];
[NnN Nn]=size(x);
%%%%%%用AR模型
for n=1:sqrt(Nn)
x2=x((n+1):Nn);
for j=1:(Nn-n)
for ji=1:n
x1(j,ji)=x(n-ji+j);
end
end
dd=pinv((x1'*x1))*x1'*x2';

p=dd;%%%%%%%%p为行向量
for j=n+1:Nn
for ii=1:n
ax(ii)=p(ii)*x(j-ii);%%%%拟合出的xt
end
ct(j)=x(j)-sum(ax);%%%%%%%%白噪声序列
end
t=[];
t=ct(n+1:Nn).^2;
b(n)=sum(t)/(Nn-n);
AIC(n)=log(b(n))+2*n/Nn;%%%%%%%%AIC准则
x1=[];x2=[];
end
%%%%%%%%%%%%%%%%%求AIC的最小值
pp=AIC(1);
nn=1;
for i=1:sqrt(Nn)
if pp>AIC(i)
pp=AIC(i);
nn=i;
end
end
%%%%%%%%%%%%%%%AIC中最小值为pp，回归节数为nn
%%%%%%%%%求白噪声序列%%%%%%
for j=nn+1:Nn
for ii=1:nn
ax(ii)=p(ii)*x(j-ii);
end
ct(j)=x(j)-sum(ax);%%%%%%%%%白噪声序列%%%%%%
end
at0(nn+1:Nn)=ct(nn+1:Nn);
%%%%%%%%%%%%%估计参数aa bb
jies=12;
if jies>=Nn
jies=Nn-1;
end
for p=1:jies
for q=1:jies
hh=[p,q,nn];
L=max(hh);

X=[];e1=[];YY=[];Y=[];m=[];n=[];aabb=[];pa=[];pb=[];
for j=1:Nn-L
for km=1:p
X(j,km)=x(L+j-km);
end
end


for j=1:Nn-L
for jm=1:q
e1(j,jm)=at0(L-jm+j);
end
end

for i=1:Nn-L
YY(i)=x(L+i);
end
Y=YY';
XX=[X e1];


aabb=pinv(XX)*Y;



pa=aabb(1:p,1);%%%%%%%%%%%pa为自回归系数
pb=aabb(p+1:p+q,1);%%%%%%%%%%%pb为移动平均系数
%%%%%%%%%%%%计算AICB（p,q)
for j=L+1:Nn
for ii=1:p
ax(ii)=pa(ii)*x(j-ii);
end
for kk=1:q
bx(kk)=pb(kk)*x(j-kk);
end
abtt(j)=x(j)-sum(ax)-sum(bx);%%%%%%白噪声序列%%%%%%
end
t=[];
t=abtt(L+1:Nn).^2;
ab(p,q)=sum(t)/(Nn-L);
AICB(p,q)=log(ab(p,q))+2*(p+q)/Nn;
end
end
%%%%%%%%%%%%%%%%求AICB(p,q)的最小值
pppp=AICB(1,1);
nna=1;
nnb=1;
for ik=1:jies
for jk=1:jies
if pppp>AICB(ik,jk)
pppp=AICB(ik,jk);
nna=ik;
nnb=jk;
end
end
end
%%%%%%%%%%%%%%%%确定模型为ARMA（p,q)%%%%%%%%%%求自回归系数和移动平均系数

X=[];e1=[];
p=nna;
q=nnb;
hh=[p,q];
L=max(hh);


for j=1:Nn-L
for km=1:p
X(j,km)=x(L+j-km);
end
end

for j=1:Nn-L
for jm=1:q
e1(j,jm)=at0(L-jm+j);
end
end

for i=1:Nn-L
YY(i)=x(L+i);
end
Y=YY';

XX=[X e1];

aabb=pinv(XX)*Y;


abc=aabb;
pma=[];pmb=[];
pma=abc(1:p,1)';%%%%%%%%%自回归系数（行向量
pmb=abc(p+1:p+q,1)';%%%%%%%%%%移动平均系数（行向量
%%%%%%%%%%%模型ARMA（p,q)%%%%%%自回归系数pma%%%%%%%移动平均系数pmb
%%%

%%signa=sum( Y-X*pma'-e1*pmb').^2/(Nn-L)
%%%%%%%%%%%%%%%预测%%%%%%%%%%%%%%%

yyuce=[x(1:L)';X*pma'+e1*pmb'];


ee=(Y-X*pma'-e1*pmb');
ee=[zeros(L,1);ee; zeros(cishu,1)]';%%%%%%%白噪声（L+1:Nn)（行向量

c=zeros(1,p);
cc=zeros(1,q);
for i=1:p
c(i)=x(Nn-i+1);
end
for i=1:q
cc(i)=ee(Nn-i+1);
end
z1=c*pma'+cc*pmb';
z(nt,1)=z1;
for k=2:cishu
if p==1
z2=z1;
z1=z2*pma'+cc*pmb';
end
if p>1
for i=1:p-1
c(p-i+1)=c(p-i);
end
c(1)=z1;
for i=1:q
cc(i)=ee(Nn-i+k);
end
z1=c*pma'+cc*pmb';
end
z(nt,k)=z1;
end

z(nt,1)=z(nt,1)+xx(Nn);
for k=2:cishu
z(nt,k)=z(nt,k)+z(nt,k-1);
end
yyuce=(yyuce)';
yuce(1)=xx(1);
for i=1:Nn-1
yuce(i+1)=yyuce(i)+xx(i+1);
end

yucezhi(nt,:)=yuce;

end
yucezhi=[yucezhi z];

end