clc;clear;
warning off all
data=xlsread('data.xlsx');
%预报步长
step=24;
TempData=data;
output=data;
%去趋势线
TempData=detrend(TempData);
%趋势函数
TrendData=data-TempData;
%差分，平稳化时间序列
[H,PValue,TestStat,CriticalValue] = adftest(TempData);
difftime=0;
SaveDiffData=[];
while ~H
%差分，平稳化时间序列
SaveDiffData=[SaveDiffData,TempData(1,1)];
TempData=diff(TempData);
%差分次数
difftime=difftime+1;
%adf检验，判断时间序列是否平稳化
[H,PValue,TestStat,CriticalValue] = adftest(TempData);
end
%模型定阶或识别
test = [];
%自回归对应PACF,给定滞后长度上限p和q
for p = 0:5
%移动平均对应ACF
for q = 0:5
if(p + q ~= 0)
m = armax(TempData,[p q]);
%armax(p,q),计算AIC
AIC = aic(m);
test = [test;p q AIC];
end
end
end
for k = 1:size(test,1)
%选择AIC值最小的模型
if test(k,3) == min(test(:,3))
p_test = test(k,1);
q_test = test(k,2);
break;
end
end
%armax(p,q),[p_test q_test]对应AIC值最小
m = armax(TempData,[p_test q_test]);
ARIMA_Predict = predict(m,TempData,step);
ARIMA_Forcast = forecast(m,TempData,step);
PreR=[ARIMA_Predict' ARIMA_Forcast'];
%还原差分
if size(SaveDiffData,2)~=0
for index=size(SaveDiffData,2):-1:1
%差分还原
PreR=cumsum([SaveDiffData(index),PreR]);
end
end
%预测趋势并返回结果
TrendData=TrendData';
mp1=polyfit(1:size(TrendData,2),TrendData,1);
xt=[];
for j=1:step
xt=[xt,size(TrendData,2)+j];
end
TrendResult=polyval(mp1,xt);
ARIMAoutput=[TrendData,TrendResult]+PreR;
ARIMAoutput_Predict=ARIMAoutput(1:end-step)';
ARIMAoutput_Forcast=ARIMAoutput(end-step+1:end)';
%% 数据输出
error=ARIMAoutput_Predict-output;
pererror=error./output;
avererror=sum(abs(error))/length(error);
averpererror=sum(abs(pererror))/length(pererror);
RMSE = sqrt(mean((error).^2));
disp('ARIMA模型预测平均绝对误差MAE');
disp(avererror)
disp('ARIMA模型预测平均绝对误差百分比MAPE');
disp(averpererror)
disp('ARIMA模型预测均方根误差RMSE')
disp(RMSE)
disp('未来年预报')
disp(ARIMAoutput_Forcast)
%% 可视化
figure
plot(ARIMAoutput,':.','linewidth',2);
hold on
plot(data,'k--','linewidth',2);
legend('预测值', '实际值','Location','NorthWest')
title('ARIMA模型预测')
xlabel('时间')
ylabel('结果')
axis tight
hold off
%-------------------------------------------------------------------------------------
figure
subplot(2,1,1)
autocorr(data)
subplot(2,1,2)
parcorr(data)
%-------------------------------------------------------------------------------------
figure()
plot(pererror,'-.','Color',[128 0 255]./255,'linewidth',2)
legend('ARIMA模型预测相对误差','Location','NorthEast','FontName','华文宋体')
title('ARIMA模型相对误差图','fontsize',12,'FontName','华文宋体')
ylabel('误差','fontsize',12,'FontName','华文宋体')
xlabel('样本','fontsize',12,'FontName','华文宋体')
axis tight
%-------------------------------------------------------------------------------------
figure();
histogram(error);
title(['Error Mean = ' num2str(mean(error)) ', Error StD = ' num2str(std(error))]);
xlabel('Error Histogram');
