y=xlsread('f:\高锰酸钾指数.xls','A1:A72');
%figure(2)
%subplot(2,1,1)
%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] = adftest(DX); %是否是稳定序列
for i = 1:10
if H== 1
break;
else
DX=diff(y,i); %进行差分
end
end
figure(3);
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))内
%对差分后的序列做拟合和预测,求出最好的阶数
DX=reshape(y,72,1);
z=[DX;zeros(12,1)];
z=iddata(DX);%将DX转化为matlab接受的格式
test = [];
% p = [0 1 2 3]; %自回归对应PACF,给定滞后长度上限p和q,一般取为T/10、ln(T)或T^(1/2),这里取
% T/10=12;
% q = [0 1 2 3]; %移动平均对应ACF
for p = 1:3 %自回归对应PACF,给定滞后长度上限p和q,一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q = 1:3 %移动平均对应ACF
m = armax(z(1:60),[p,q]);
AIC = aic(m); %armax(p,q),选择对应FPE最小,AIC值最小的模型
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:60),[p_test q_test]); %armax(p,q),[p_test q_test]对应AIC值最小
figure(4)
e = resid(m1,z); %拟合做残差分析
plot(e);
grid;
title('残差序列');
%检验残差的自相关和偏相关函数
figure(5)
subplot(2,1,1)
autocorr(e.OutputData) %一阶差分序列z自相关函数图MA(q),置信水平0.95
subplot(2,1,2)
parcorr(e.OutputData) %一阶差分序列z偏相关函数图AR(p),置信水平0.95
%预测过程
p=predict(m1,z,1);%lcl m是任何idmodel对象(idpoly,idss,idgrey或idarx),z是包含输入输出数据的iddata对象。1是预测时间长度
% p是模型的预测输出,为包含预测值的iddata对象
这是程序只可以看到72个预测到的数据 要想预测更多数据比如84个、96个怎么写程序
希望高手们给指导一下
谢谢!