本期分享一个关于ARIMA电价预测的MATLAB代码实例。
ARIMA(Autoregressive Integrated Moving Average)模型是一种常用于时间序列分析和预测的统计模型。它结合了自回归(AR)模型、差分(Integrated)过程和移动平均(MA)模型的特性。ARIMA模型作为一种经典的时间序列预测方法,能够有效地捕捉时间序列的自相关性和趋势而广泛应用于实际问题中的预测和分析任务。
1. ARIMA的组成部分:
-
自回归(AR)部分:表示当前观测值与过去若干个观测值的线性关系。AR(p)模型中,当前值与p个滞后值相关。
-
差分(Integrated)部分:为了使时间序列平稳(stationary),需要进行差分操作,即对序列的当前值与之前的值相减,直到得到平稳序列。
-
移动平均(MA)部分:表示当前观测值与过去若干个随机误差的线性关系。MA(q)模型中,当前值与q个滞后的误差相关。
2. ARIMA的预测机理:
-
模型建立:通过观察时间序列的自相关图(ACF)和偏自相关图(PACF),确定合适的p(自回归阶数)和q(移动平均阶数),以及是否需要进行差分(d阶数)来使序列平稳。
-
参数估计:利用最大似然估计等方法,估计模型的参数。
-
预测:通过已有的历史数据,结合ARIMA模型中的参数,进行未来一段时间的预测。预测结果不仅依赖于历史数据的值,还依赖于残差的预测。
ARIMA模型的AIC(赤池信息准则)和BIC(贝叶斯信息准则)是模型选择中常用的两个准则,用于在多个候选模型中选择最优模型。在ARIMA模型拟合过程中,通常会尝试多个可能的ARIMA模型(不同的p、d、q值组合),然后根据AIC和BIC的数值选择最优模型。选取AIC或BIC较小的模型作为最终的ARIMA模型。通过使用AIC和BIC,可以帮助确定在平衡模型拟合优良程度与模型复杂度之间的最佳点,从而提高时间序列预测模型的准确性和可靠性。
本期代码以AIC和BIC为指标,计算ARIMA模型种不同的p, d, q组合,求出最佳的p,d,q值。
代码如下:
clc;clear;
% 1. 读取数据
data = readmatrix('electricity_dah_prices2.csv');
time_series_data = data(101:160,3); %取出100个数据点用于本次的训练预测测试
% 2. 划分训练集和测试集 - 这里使用80%的数据作为训练集
train_size = round(length(time_series_data) * 0.9);
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 = 1;
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
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 = 6; % 预测训练集之后的6个时刻的数据
[Y_pre,Y_pre_RMSE] = forecast(EstMdl, num_steps, 'Y0',train_data);
% 计算 95% 置信区间
z = norminv(0.95);
Y_pre_CI = [Y_pre - 1.96 .* sqrt(Y_pre_RMSE), Y_pre + 1.96.* sqrt(Y_pre_RMSE)];
% 8. 输出预测结果
disp(['预测结果(', num2str(num_steps), '个步长):']);
disp(Y_pre);
disp(['预测置信区间(', num2str(num_steps), '个步长):']);
disp(Y_pre_CI);
% 9. 可视化预测结果
figure;
fill([train_size+1:train_size+num_steps, train_size+num_steps : -1 : train_size+1], [Y_pre_CI(:,1)', flip(Y_pre_CI(:,2)')], ...
'r', 'FaceColor', [1, 0.8, 0.8], 'EdgeColor', 'none')
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, Y_pre, 'r', 'LineWidth', 1);hold on
xlim([1, length(time_series_data) + num_steps]);
title('ARIMA 时间序列预测');
xlabel('时间');
ylabel('电价');
legend('95%置信区间','历史数据', '真实数据', '预测数据', '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)]);
预测未来6个时刻的数据,并绘制置信区间。结果图如下:
代码获取
完整代码获取,后台回复关键词:
ARIMA预测