ARIMA数据分析.pdf

ARIMA数据分析

clear;

%原序列的图形
figure(1);
subplot(2,1,1)
plot(y1);
title('洛阳电厂08年到煤量');
grid;
subplot(2,1,2)
plot(y2);
title('洛阳电厂09年到煤量');
grid;
年的到煤量
%figure(11)
%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(y);%是否是稳定序
列
fori=1:10
ifH~=1
break;
else
DX=diff(y,i);%进行差分
[H,PValue,TestStat,CriticalValue]=ADFTEST(DX);
end
end

figure(2);
subplot(211)
plot(DX(1:365));
title('洛阳电厂08年到煤量差分之后的序列');
grid;
subplot(212)
plot(DX(366:700));
grid;
title('洛阳电厂09年到煤量差分之后的序列');
%plot(DX);%进行差分之后的数列

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))内

%对差分后的序列做拟合和预测，求出最好的阶数
z=iddata(DX);%将DX转化为matlab接受的格式
test=[];
%p=[0123456789101112];%自回归对应PACF,给定滞
后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取
%T/10=12;
%q=[0123456789101112];%移动平均对应ACF
forp=1:10%自回归对应PACF,给定滞后长度上限p和q，一般
取为T/10、ln(T)或T^(1/2),这里取T/10=12
forq=1:10%移动平均对应ACF
m=armax(z(1:500),[p,q]);
AIC=aic(m);%armax(p,q),选择对应FPE最小，AIC值最小的模
型
%[H,P,Qstat,CV]=lbqtest(z,[p;q],0.05)%Ljung-BoxQ-
statisticlack-of-fithypothesistest
test=[test;pqAIC];
end
end
fork=1:size(test,1)
iftest(k,3)==min(test(:,3))%选择AIC值最小的模型
p_test=test(k,1);
q_test=test(k,2);
break;
end
end

%拟合过程
m1=armax(z(1:500),[p_testq_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);%lclm是任何idmodel对象(idpoly，idss，
idgrey或idar
x)，z是包含输入输出数据的iddata对象。1是预测时间长度
%p是模型的预测输出，为包含预测值的iddata对象
figure(6)
subplot(211)
plot(z(501:531),'r*',z(501:531));
title('原始序列');
grid;
subplot(212)
plot(p(501:531),'b*',p(501:531));
title('预测序列');
grid;
po=p.OutputData;
zo=z.Outp