matlab的arma程序-经管之家官网！

请各位高手看一下，这个程序有问题吗？邮箱：longzhu147@126.com
clear all;clc;
y=xlsread('E:\guo\A3.xls','A1:A96');%填写待处理数据的详细
ls=length(y);
figure(1);%第1张图
plot(y);
set(gca,'Xlim',[0 ls]);
figure(2);
subplot(2,1,1);
set(gca,'Xlim',[0 ls]);
autocorr(y);
%原序列的自相关函数图MA(q),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
subplot(2,1,2)
parcorr(y); %原序列的偏相关函数图AR(p)，观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
%如果该序列不是平稳的做差分图，否则跳过该步
DX=y;
[H,PValue,TestStat,CriticalValue]=dfARDTest(y,[],0.05,'T');%是否是稳定序列
for i = 1:10
if H == 1
break;
else
DX=diff(y,i);%进行差分
[H,PValue,TestStat,CriticalValue] = dfARDTest(DX,[],0.05,'T');
end
end
figure(3);
plot(DX);%进行差分之后的数列
set(gca,'Xlim',[0 ls]);
figure(4);
subplot(2,1,1)
autocorr(DX); %差分序列DX自相关函数图MA(q),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
subplot(2,1,2)
parcorr(DX); %差分序列DX偏相关函数图AR(p)，观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
set(gca,'Xlim',[0 ls]);
%对差分后的序列做拟合和预测，求出最好的阶数
z=iddata(DX);%将DX转化为matlab接受的格式
test = [];
% p = [0 1 2 3 4 5 6 7 8 9 10 11 12]; %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取
%T/10=12;
% q = [0 1 2 3 4 5 6 7 8 9 10 11 12]; %移动平均对应ACF
for p = 1:10 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q = 1:10 %移动平均对应ACF
m = armax(z(1:96),[p q]);
AIC = aic(m); %armax(p,q),选择对应FPE最小，AIC值最小的模型
%[H, P, Qstat, CV] = lbqtest(z, [p;q], 0.05) %Ljung-Box Q-statistic lack-of-fit hypothesis test
test = [test;p q AIC];
end
end
for k = 1:size(test,1)
if test(k,3) == min(test(:,3)) %选择AIC值最小的模型
p_test = test(k,1);
q_test = test(k,2);
break;
end
end
%拟合过程
m1 = armax(z(1:96),[p_test q_test]); %armax(p,q),[p_test q_test]对应AIC值最小
figure(5);
e = resid(m1,z); %拟合做残差分析
plot(e);
set(gca,'Xlim',[0 ls]);
%检验残差的自相关和偏相关函数
figure(6);
subplot(2,1,1)
autocorr(e.OutputData) %一阶差分序列z自相关函数图MA(q)，置信水平0.95
subplot(2,1,2)
parcorr(e.OutputData)
set(gca,'Xlim',[0 ls]);%一阶差分序列z偏相关函数图AR(p)，置信水平0.95
%预测过程
p=predict(m1,z,1);
x=[1:1:95];
po = p.OutputData
figure(7)
plot(x,po,'r');
hold on
x=[1:1:96];
plot(x,y,'b');
set(gca,'Xlim',[0 ls]);
hold off
[Pr,DWr] = dwtest(e.OutputData,z.OutputData);
if Pr<0.05
disp('can not use this model');
else disp('can use this model');
end