% 运行之前注意保存已有数据
clc;clear;
close all;
% 是否强制确定ARIMA模型(不理想的情况下选择1)
isKill = 1;
% 预测时间长度
PreTime = 100;
% 是否进行新时间段的预测,建议先用0再用1
isPre = 1;
% 导入数据集
data = readtable('data.xlsx');
y = data.Value;
% 拆分数据集为训练集和测试集
train_size = floor(0.8 * length(y));
train_data = y(1:train_size);
if isPre == 1
train_data = y;
end
test_data = y(train_size+1:end);
% 平稳性检验 - ADF检验
[h_adf, p_adf] = adftest(y);
fprintf('ADF检验结果:\n');
fprintf('是否拒绝原假设(序列非平稳): %d\n', h_adf);
fprintf('p值: %.4f\n', p_adf);
% 平稳性检验 - KPSS检验
[h_kpss, p_kpss] = kpsstest(y);
fprintf('KPSS检验结果:\n');
fprintf('是否拒绝原假设(序列平稳): %d\n', ~h_kpss);
fprintf('p值: %.4f\n', p_kpss);
if isKill == 0
% 差分确定d参数
if h_adf == 1 && h_kpss == 0
d = 0;
disp('序列已经平稳')
elseif h_adf == 0 && h_kpss == 1
d = 1;
disp('序列需要一阶差分')
elseif h_adf == 1 && h_kpss == 1
d = 2;
disp('序列需要二阶差分')
else
error('无法确定d参数');
end
% 差分处理
if d > 0
y_diff = diff(train_data, d);
else
y_diff = train_data;
end
else
y_diff = train_data;
end
% 自相关法确定p参数
figure
autocorr(y_diff)
title('Autocorrelation')
% 偏自相关法确定q参数
figure
parcorr(y_diff)
title('Partial Autocorrelation')
% 训练ARIMA模型
if isKill == 0
model = arima('D', d, 'ARLags', [1 2], 'MALags', 1); % ARIMA(p, d, q) model
else
model = arima(2, 1, 1);
end
model = estimate(model, y_diff);
% 测试未来时间预测情况
if isPre == 0
horizon = length(y) - train_size;
else
horizon = PreTime;
end
[forecastTest, YMSE, forecastCI] = forecast(model, horizon, 'Y0', y_diff);
ResultScale(:, 1) = forecastTest - 1.64*sqrt(YMSE); %90置信区间下限 %95 改成1.96
ResultScale(:, 2) = forecastTest + 1.64*sqrt(YMSE); %90置信区间上限
% 自动确定的可能修正方法
if isKill == 0
forecastTest = forecastTest + train_data(end) - forecastTest(1);
ResultScale(:, 1) = forecastTest - 1.64*sqrt(YMSE); %90置信区间下限
ResultScale(:, 2) = forecastTest + 1.64*sqrt(YMSE); %90置信区间上限
end
if isPre == 0
% 可视化测试结果
figure
plot(y, 'ko-')
hold on
plot(train_size+1:length(y), forecastTest, 'b-', 'LineWidth', 1.5)
xlabel('Time')
ylabel('Value')
legend('Observations', 'Predictions')
title('ARIMA Forecast Results with Prediction Range')
% 预测结果取值范围可视化
plot(train_size+1:length(y), ResultScale(:, 1), 'r--', 'LineWidth', 1.5)
hold on
plot(train_size+1:length(y), ResultScale(:, 2), 'r--', 'LineWidth', 1.5)
fill([train_size+1:length(y), fliplr(train_size+1:length(y))], [ResultScale(:, 1)', fliplr(ResultScale(:, 2)')], 'r', 'FaceAlpha', 0.3)
plot(train_size+1:length(y), forecastTest, 'b-', 'LineWidth', 1.5)
xlabel('Time')
ylabel('Value')
legend('Lower Bound', 'Upper Bound', 'Predictions')
else
train_size = length(y);
figure
plot(y, 'ko-')
hold on
plot(train_size+1:length(y)+PreTime, forecastTest, 'b-', 'LineWidth', 1.5)
xlabel('Time')
ylabel('Value')
legend('Observations', 'Predictions')
title('ARIMA Forecast Results with Prediction Range')
% 预测结果取值范围可视化
plot(train_size+1:length(y)+PreTime, ResultScale(:, 1), 'r--', 'LineWidth', 1.5)
hold on
plot(train_size+1:length(y)+PreTime, ResultScale(:, 2), 'r--', 'LineWidth', 1.5)
fill([train_size+1:length(y)+PreTime, fliplr(train_size+1:length(y)+PreTime)], [ResultScale(:, 1)', fliplr(ResultScale(:, 2)')], 'r', 'FaceAlpha', 0.3)
plot(train_size+1:length(y)+PreTime, forecastTest, 'b-', 'LineWidth', 1.5)
xlabel('Time')
ylabel('Value')
legend('Lower Bound', 'Upper Bound', 'Predictions')
end
if isPre == 0
% 误差分析
error = y(train_size+1:end) - forecastTest;
mse = mean(error.^2);
rmse = sqrt(mse);
mae = mean(abs(error));
fprintf('Mean Squared Error (MSE): %.4f\n', mse);
fprintf('Root Mean Squared Error (RMSE): %.4f\n', rmse);
fprintf('Mean Absolute Error (MAE): %.4f\n', mae);
% 残差检验
residuals = infer(model, y_diff);
figure
plot(residuals)
xlabel('Time')
ylabel('Residuals')
title('Residuals')
% 残差自相关检验
figure
autocorr(residuals)
title('Autocorrelation of Residuals')
% 残差偏自相关检验
figure
parcorr(residuals)
title('Partial Autocorrelation of Residuals')
end
运行结果如下图所示: