%-----------------------------单步预测 预测当前时刻后的六个数据----------------------------------%
XT=[L(84:94)];
for t=6:11
m(t)=0;
for i=1:5
m(t)=W(i)*XT(t-i)+m(t);
end
end
%---------------------------绘制数据模型逼近曲线-----------------------------------%
for t=6:88
r=0;
for i=1:5
r=W(i)*Y(t-i)+r;
end
at= Y(t)-r;
end
figure;
for t=6:88
y(t)=0;
for i=1:5
y(t)=W(i)*Y(t-i)+y(t);
end
y(t)=y(t)+at;
y(t)=Yt(t+1)-y(t);
y(t)=P(t+1)-y(t);
end
plot(y,'r-*'); % 样本数据模型逼近曲线
hold on;
plot(91:96,z1,'r-*');
hold on;
plot(P,'--'); % 原样本曲线
title('AR(5)模型样本逼近预测曲线');
pause
%-----------------------------绘制数据模型逼近曲线-----------------------------------%
Z(1)=W(1)*Xt(5)+W(2)*Xt(4)+W(3)*Xt(3)-W(4)*Xt(2)-W(5)*Xt(1)
%------求向前l步的预测值
%预测步数小于5时
for l=2:5
K(l)=0;
for i=1:l-1
K(l)=W(i)*Z(l-i)+K(l);
end
G(l)=0;
for j=l:5
G(l)=W(j)*Xt(5+l-j)+G(l);
end
Z(l)=K(l)+G(l);
end
%预测步数大于5时(向前6步预测)
for l=6:6
K(l)=0;
for i=1:5
K(l)=W(i)*Z(l-i)+K(l);
end
Z(l)=K(l);
end
%----预测值进行反标准化处理
r=Z*v+Ux % 0.0581 0.0844 0.0156 0.0319 0.0632 0.0652
r(1)=Yt(90)+r(1); %一次反差分
z(1)=P(90)+r(1) %二次反差分
for i=2:6
r(i)=r(i-1)+r(i);
z(i)=z(i-1)+r(i)
end
%----------------------由于时间序列有不平稳趋势,进行两次差分运算,消除趋势性----------------------%
for i=2:96
Yt(i)=P(i)-P(i-1);
end
for i=3:96
L(i)=Yt(i)-Yt(i-1);
end
figure;
L=L(3:96);
Y=L(1:88);
plot(P);
title('原数据序列图');
hold on;
pause
plot(Y,'r');
title('两次差分后的序列图和原数对比图');
pause
%--------------------------------------对数据标准化处理----------------------------------------------%
Ux=sum(Y)/88 % 求序列均值
yt=Y-Ux;
b=0;
for i=1:88
b=yt(i)^2/88+b;
end
v=sqrt(b) % 求序列方差
Y=(Y-Ux)/v; % 标准化处理公式
f=F(1:88);
t=1:88;
figure;
plot(t,f,t,Y,'r')
title('原始数据和标准化处理后对比图');
xlabel('时间t'),ylabel('油价y');
legend('原始数据 F ','标准化后数据Y ');
pause
%--------------------------------------对数据标准化处理----------------------------------------------%
%------------------------检验预处理后的数据是否符合AR建模要求,计算自相关和偏相关系数---------------%
%---------------------------------------计算自相关系数-----------------------------------%
R0=0;
for i=1:88
R0=Y(i)^2/88+R0;
end
R0
for k=1:20
R(k)=0;
for i=k+1:88
R(k)=Y(i)*Y(i-k)/88+R(k);
end
R %自协方差函数R
end
x=R/R0 %自相关系数x
figure;
plot(x)
title('自相关系数分析图');
pause
%-----------------------------------计算自相关系数-------------------------------------%