ARIMA时间序列预测MATLAB代码模板（无需调试） |
您所在的位置：网站首页 › 时间序列模型构建图 › ARIMA时间序列预测MATLAB代码模板（无需调试） |
ARIMA时间序列预测MATLAB代码模板（无需调试）
|
小白专用,直接改成自己的数据运行即可完成预测并画图 我的数据在评论区自取， 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;设定未来预测的步数。您可以根据需要更改这个数字。 |
【本文地址】