MATLAB数学建模——时间序列ARIMA及可视化

% 运行之前注意保存已有数据
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

运行结果如下图所示：