小白专用,直接改成自己的数据运行即可完成预测并画图
我的数据在评论区自取，
clear; clc
%小白专用，"*********《需要自己输入》**********"仅在有这种注释的地方改成自己的数据即可，一共有4个地方
DD=readmatrix("B.xlsx");%这里输入自己的单序列数据，要求行向量*********《需要自己输入》**********
P=DD(1:500,2)';
N=length(P);
n=490;%自己选取训练集个数*********《需要自己输入》**********
F = P(1:n+2);
%----------------------由于时间序列有不平稳趋势，进行两次差分运算，消除趋势性----------------------%
Yt=[0,diff(P,1)];
L=diff(P,2);%全体，比原始数据少2个，因为做了差分
Y=L(1:n); %输入
a=length(L)-length(Y);%单步预测步数
aa=a;%多步预测步数
% %画图
% figure;
% plot(P);
% title('原数据序列图');
% hold on;
% figure;
% plot(Y,'*');
% title('两次差分后的序列图和原数对比图');
%%
%--------------------------------------对数据标准化处理----------------------------------------------%
%处理的算法 : (data - 期望)/方差
Ux=sum(Y)/n ; % 求序列均值
yt=Y-Ux;
b=0;
for i=1:n
b=yt(i)^2/n+b;
end
v=sqrt(b); % 求序列方差
% Y=yt/v; % 标准化处理公式
Y=zscore(Y);
f=F(1:n);
t=1:n;
%画图
% figure;
% plot(t,f,t,Y,'r')
% title('原始数据和标准化处理后对比图');
% xlabel('时间t'),ylabel('油价y');
% legend('原始数据 F ','标准化后数据Y ',"Location","best");
%%
%--------------------------------------对数据标准化处理----------------------------------------------%
%------------------------检验预处理后的数据是否符合AR建模要求，计算自相关和偏相关系数---------------%
%---------------------------------------计算自相关系数-----------------------------------%
%%
R0=0;
for i=1:n
R0=Y(i)^2/n+R0; %标准化处理后的数据的方差
end
for k=1:20
%R 协方差
R(k)=0;
for i=k+1:n
R(k)=Y(i)*Y(i-k)/n+R(k);
end
end
x=R/R0 ; %自相关系数x = 协方差/方差
%画图
% figure;
% plot(x)
% title('自相关系数分析图');
%%
%-----------------------------------计算自相关系数-------------------------------------%
%-----------------------解Y-W方程，其系数矩阵是Toeplitz矩阵(多普里兹矩阵)。求得偏相关函数X-------------------
X1=x(1);xx(1,1)=1;X(1,1)=x(1);B(1,1)=x(1);
K=0;T=X1;
for t=2:n
at=Y(t)-T(1)*Y(t-1);
K=(at)^2+K;
end
U(1)=K/(n-1) ; % 1阶模型残差方差
for i =1:19
B(i+1,1)=x(i+1);
xx(1,i+1)=x(i);
A=toeplitz(xx);
XX=A\B; %x=a\b是方程a*x =b的解
XXX=XX(i+1);
X(1,i+1)=XXX;
K=0;T=XX;
for t=i+2:n
r=0;
for j=1:i+1
r=T(j)*Y(t-j)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(i+1)=K/(n-i+1); %计算20阶以内的模型残差方差
end
%-----------------------------------解Y-W方程，得偏相关函数X-------------------------------------%
% figure;
% plot(X);
% title('偏相关函数图');%自己要根据图先判断阶次
%%
q=20;%猜测阶数,通过看上面偏相关图,*********《需要自己输入》**********
%-----根据偏相关函数截尾性，初判模型阶次为5。用最小二乘法估计参数,计算20阶以内的模型残差方差和AIC值，应用AIC准则为模型定阶--%
S(1,1)=R0;
for i = 1:q-1
S(1,i+1)=R(i);
end
G=toeplitz(S);
%inv(G)返回G的反函数
W=inv(G)*[R(1:q)]' ; % 参数W(i) 与X5相同 G*W = [R(1:5)]'
U=20*U ;
for i=1:20
AIC2(i)=n*log(U(i))+2*(i) ;
end
% 比如AIC2值为：172.6632 165.4660 153.2087 145.1442 140.7898 141.6824 142.9944 144.5601 146.3067 148.7036
%-----------------取使AIC值为最小值的阶次，判断模型阶次为5。用最小二乘法估计参数--------------------%
%%
q=20;%确定阶数 ，通过看AIC2值最小的位置,*********《需要自己输入》**********
%------------------检验{at}是否为白噪声。求{at}的自相关系数，看其是否趋近于零-----------------------%
C=0;K=0;
for t=q+2:n
at=Y(t)+Y(q+1);
for i=1:q
at=-W(i)*Y(t-i)-W(i)*Y(q-i+1)+at;
end
at1=Y(t-1);
for i=1:q
at1=-W(i)*Y(t-i-1)+at1;
end
%at1=Y(t-1)-W(1)*Y(t-2)-W(2)*Y(t-3)-W(3)*Y(t-4)-W(4)*Y(t-5)-W(5)*Y(t-6);
C=at*at1+C;
K=(at)^2+K;
end
p=C/K ; %若p接近于零，则{at}可看作是白噪声
%--------------------------------{at}的自相关系数,趋近于零，模型适用---------%
%%
%------------AR(5)模型方程为-----------------------------------------------%
% X(t)=W(1)*X(t-1)-W(2)*X(t-2)-W(3)*X(t-3)-W(4)*X(t-4)-W(5)*X(t-5)+at
%------------------------------------------后六年的数据 进行预测和效果检验------------------------------------%
%注意注意注意a为测试集的元素个数
%-----------------------------单步预测 预测当前时刻后的a个数据---------------%
XT=[L(n-q+1:n+a)];
for t=q+1:q+a
m(t)=0;
for i=1:q
m(t)=W(i)*XT(t-i)+m(t);
end
end
m=m(q+1:q+a);
%-------------预测值进行反处理---------------%
for i =1:a
m(i)=Yt(n+i+1)+m(i); %一次反差分
z1(i)=P(n+i+1)+m(i);%二次反差分
end
% z1 ; % 单步预测的向后6个预测值
%---------------------------绘制数据模型逼近曲线----------------------------%
for t=q+1:n
r=0;
for i=1:q
r=W(i)*Y(t-i)+r;
end
at= Y(t)-r;
end
figure;
for t=q+1:n
y(t)=0;
for i=1:q
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(n+2:n+a+1,z1,'r-*'); %向后a布预测
hold on;
plot(P,"--"); % 原样本曲线
title('AR(q)模型样本逼近预测曲线');
legend("训练样本预测值","测试集预测值","真实值","Location","best");
%-------------------------检测单步预测误差
D_a=P(n+2:end-1);
for i=1:a
e6_a(i)=D_a(i)-z1(i);
PE6_a(i)= (e6_a(i)/D_a(i))*100;
end
e6_a; % 多步预测的绝对误差
PE6_a; % 多步预测的相对误差
1-abs(PE6_a); % 准确率
%------多步预测平均绝对误差
mae6_a=sum(abs(e6_a)) /6 ;
%------多步预测平均绝对百分比误差
MAPE6_a=sum(abs(PE6_a))/6
%------绘制预测结果和实际值的比较图
figure;
plot(1:a,D_a,'-+')
hold on;
plot(z1,'r-*');
title('单步，向后a步预测值和实际值对比图');
legend("真实值","预测值","Location","best");
hold off;
%%
%-----------------------------绘制数据模型逼近曲线--------------------------%
%-------------------------预测误差分析（多步）------------------------%
%----------------------------------多步预测 目的是向后aa步预测--------------%
Z(1)=0;Xt=0;
for i =1:q
Xt(1,i)=Y(n-q+i);
end
%Xt=[ Y(n-4) Y(n-3) Y(n-2) Y(n-1) Y(n)]; %取当前时刻之前的q个数据
for i =1:q
Z(1)=W(i)*Xt(q-i+1)+Z(1);
end
%Z(1)=W(1)*Xt(5)+W(2)*Xt(4)+W(3)*Xt(3)-W(4)*Xt(2)-W(5)*Xt(1) ;
%------求向前l步的预测值
%预测步数小于q时
for l=2:q
K(l)=0;
for i=1:l-1
K(l)=W(i)*Z(l-i)+K(l);
end
G(l)=0;
for j=l:q
G(l)=W(j)*Xt(q+l-j)+G(l);
end
Z(l)=K(l)+G(l);
end
%预测步数大于q时（向前aa步预测）
for l=q+1:aa
K(l)=0;
for i=1:q
K(l)=W(i)*Z(l-i)+K(l);
end
Z(l)=K(l);
end
%----预测值进行反标准化处理
r=Z*v+Ux ;
r(1)=Yt(n+2)+r(1); %一次反差分
z(1)=P(n+2)+r(1) ; %二次反差分
for i=2:aa
r(i)=r(i-1)+r(i);
z(i)=z(i-1)+r(i) ;
end
%%
%---------------------------- 预测误差分析 ------------------------------%
%-------计算绝对误差和相对误差
D=P(n+2:end-1);
for i=1:aa
e6(i)=D(i)-z(i);
PE6(i)= (e6(i)/D(i))*100;
end
e6 ; % 多步预测的绝对误差
PE6 ; % 多步预测的相对误差
1-abs(PE6) ; % 准确率
%------多步预测平均绝对误差
mae6=sum(abs(e6)) /6 ;
%------多步预测平均绝对百分比误差
MAPE6=sum(abs(PE6))/6
%------绘制预测结果和实际值的比较图
figure;
plot(1:aa,D,'-+')
hold on;
plot(z,'r-*');
title('多步，向后aa步预测值和实际值对比图');
legend("真实值","预测值","Location","best");
hold off;原程序运行图窗结果：
命令行输出结果：
MAPE6_a =
3.6147
MAPE6 =
12.5103由上可知：单步预测准确率约为96.4%
多步（这里是8步）预测准确率约为87.5%
方法2：用arima函数实现时间序列预测
clc;clear;
% 1. 读取数据 - 请将'B.xlsx'替换为您的数据文件名，并将'data(:,2)'根据要预测的列确定
data = readmatrix('B.xlsx');
time_series_data = data(:,2);
% 2. 划分训练集和测试集 - 这里使用80%的数据作为训练集，您可以根据需要调整比例
train_size = round(length(time_series_data) * 0.8);
train_data = time_series_data(1:train_size);
test_data = time_series_data(train_size+1:end);
% 3. 初始化最小AIC和BIC以及最优参数 - 选择模型参数的范围（p、d、q的最大值）
max_p = 5;
max_d = 2;
max_q = 5;
min_aic = Inf;
min_bic = Inf;
best_p = 0;
best_d = 0;
best_q = 0;
% 4. 循环遍历不同的p, d, q值，尝试拟合ARIMA模型，并计算AIC和BIC
for p = 0:max_p
for d = 0:max_d
for q = 0:max_q
% 创建ARIMA模型
Mdl = arima(p, d, q);
% 拟合模型，并计算AIC和BIC
try
[EstMdl,~,logL] = estimate(Mdl, train_data, 'Display', 'off');
[aic, bic] = aicbic(logL, p + q + 1, length(train_data));
catch
continue;
end
% 更新最优参数
if bic < min_bic
min_aic = aic;
min_bic = bic;
best_p = p;
best_d = d;
best_q = q;
end
end
end
end
% 5. 使用最优参数创建ARIMA模型
best_mdl = arima(best_p, best_d, best_q);
% 6. 拟合模型
EstMdl = estimate(best_mdl, train_data);
% 7. 对测试集数据后的值进行预测 - 设定预测步长
num_steps = 20; % 预测测试集之后的20天数据
[forecast,forecast_RMSE] = forecast(EstMdl, num_steps, 'Y0', train_data);
% 计算 95% 置信区间
z = norminv(0.975);
forecast_CI = [forecast - z * forecast_RMSE, forecast + z * forecast_RMSE];
% 8. 输出预测结果
disp(['预测结果（', num2str(num_steps), '个步长）:']);
disp(forecast);
disp(['预测置信区间（', num2str(num_steps), '个步长）:']);
disp(forecast_CI);
% 9. 可视化预测结果
figure;
hold on;
plot(time_series_data, 'k', 'LineWidth', 1);hold on
plot(train_size+1:train_size+length(test_data), test_data, 'b', 'LineWidth', 1); hold on% 绘制测试集数据
plot(train_size+1:train_size+num_steps, forecast, 'r', 'LineWidth', 1);hold on
xlim([1, length(time_series_data) + num_steps]);
title('ARIMA 时间序列预测');
xlabel('时间');
ylabel('值');
legend('实际数据', '测试集数据', '预测', 'Location', 'best');
% 10. 输出模型参数
disp(['最优模型参数: p = ', num2str(best_p), ', d = ', num2str(best_d), ', q = ', num2str(best_q)]);
disp(['最小 AIC: ', num2str(min_aic)]);
disp(['最小 BIC: ', num2str(min_bic)]);这个代码将首先使用AIC和BIC选择最优ARIMA模型参数，然后使用这些参数拟合训练集数据。接下来，它将预测测试集数据后的10个步长。预测结果和预测置信区间将分别输出。您可以根据需要调整预测步长。
以下参数可以根据您的实际情况进行设定：
- 数据文件名和数据列名： 将
your-data.csv替换为您的实际数据文件名，并将your_column替换为您的数据列名。 - 训练集和测试集的划分比例： 代码中使用
train_size = round(length(time_series_data) * 0.8);将数据划分为80%的训练集和20%的测试集。您可以根据需要更改0.8这个比例。 - ARIMA模型参数： 代码中使用
Mdl = arima(p,d,q);创建一个ARIMA(p,d,q)模型。这里的数字分别代表AR（自回归）、I（差分）和MA（移动平均）的阶数。您可以根据您的时间序列数据特点选择不同的参数。 - 预测步数： 代码中使用
forecast_steps = 10;设定未来预测的步数。您可以根据需要更改这个数字。
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
