ARIMA预测模型

1. 自回归模型AR( p)

   • 描述当前值和历史值之间的关系，用变量自身的历史数据对自身进行预测，其必须要满足平稳性要求，只适用于预测与自身前期相关的现象(时间序列的自相关性)；
   • $p$ 阶自回归过程的公式定义：$$y_t=\mu +\sum _{i=1}^p{\gamma }_i{y}_{t-i}+{\epsilon }_t$$
   • $p$ 表示用几期的历史值来预测， $y_t$ 是当前值、 $\mu$ 是常数项、 $p$ 是阶数、 ${\gamma }_i$ 是自相关系数、${\epsilon }_t$ 为误差项；

2. 移动平均模型MA(q)

   • 移动平均模型关注的是自回归模型中误差项的累计；
   • $q$ 阶自回归过程的公式定义：$$y_t=\mu +\sum _{i=1}^q{\theta }_i{\epsilon }_{t-i}+{\epsilon }_t$$
   • 即时间序列当前值与历史值没有关系，而只依赖于历史白噪声的线性组合；
   • 移动平均法能有效地消除预测中的随机波动。

3. 自回归移动平均模型ARMA(p,q)

   • 自回归与移动平均的结合；
   • 公式定义：$$y_t=\mu +\sum _{i=1}^p{\gamma }_i{y}_{t-i}+{\epsilon }_t+\sum _{i=1}^q{\theta }_i{\epsilon }_{t-i}$$
   • 该式表明：
     • 一个随机时间序列可以通过一个自回归移动平均模型来表示，即该序列可以由其自身的过去或滞后值以及随机扰动项来解释；
     • 如果该序列是平稳的，即它的行为并不会随着时间的推移而变化，那么我们就可以通过该序列过去的行为来预测未来。

4. 差分整合移动平均自回归模型ARIMA(p,q,d)

   • 将 自回归模型AR( p) 、移动平均模型MA(q) 和 差分I(d) 结合，我们就得到了差分整合移动平均自回归模型ARIMA(p,q,d)；
   • $p$ 是自回归项，$q$ 为移动平均项数，$d$ 为时间序列成为平稳时所做的差分次数；
   • 原理：将非平稳时间序列转化为平稳时间序列然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。

1.算法步骤

1. 数据平稳性检验

   • 平稳性就是要求经由样本时间序列所得到的拟合曲线在未来的一段时间内仍然能够按照现有的形态延续下去，平稳性要求序列的均值和方差不发生明显变化：
     • 严平稳：序列所有的统计性质（期望，方差）都不会随着时间的推移而发生变化；
     • 宽平稳：期望与相关系数（依赖性）不变，就是说 t 时刻的值 X 依赖于过去的信息。
   • 对序列绘图，进行平稳性检验，观察序列是否宽平稳；
   • 方法：ADF 检验（Augmented Dickey-Fuller Test）；
     • ADF 大致的思想就是基于随即游走（不平稳的一个特殊序列）的，对其进行回归，如果发现 $p_{value}=1$，说明序列满足随机游走，就是非平稳的；
     • 判断标准：p-value < 0.05 则拒绝非平稳假设；
   • 若不平稳需差分处理（ d 次），d 一般不超过 3 。

2. 模型识别与定阶

   • ACF（自相关函数）确定 MA 阶数 $q$；
   • PACF（偏自相关函数）确定 AR 阶数 $p$；
   • 或者使用 AIC、BIC 遍历确定 $p、q$。

3. 参数估计

   • 极大似然估计或最小二乘法求解参数。

4. 模型诊断

   • 残差检验：Ljung-Box 检验残差是否为白噪声；
     • 

5. 预测应用

   • 使用最优模型进行未来值预测。

2.参数选择

(1)拖尾截尾判断法

• ACF（自相关函数）确定 MA 阶数 $q$：

  • 有序的随机变量序列与其自身相比较。自相关系数反映了统一序列在不同时序的取值之间的相关性，对于时间序列 $y_t$，$y_t$ 与 $y_{t-k}$ 的相关系数称为 $y_t$ 间隔 $k$ 的自相关系数，取值范围为 [-1 ,1]；
  • 公式：$$ACF(k)={\rho }_k=\frac{Cov(y_t,y_{t-k})}{Var(y_t)}$$

• PACF（偏自相关函数）确定 AR 阶数 $p$：

  • 对于一个平稳 AR§ 模型，求出滞后 $k$ 自相关系数 ${\rho }_k$ 时，实际上得到并不是 $y_t$ 与 $y_{t-k}$ 之间单纯的相关关系；
  • 因为 $y_t$ 同时还会受到中间 $k-1$ 个随机变量 $y_{t-1}、y_{t-2}、\dots 、y_{t-k+1}$ 的影响，而这 $k-1$ 个随机变量又都和 $y_{t-k}$ 具有相关关系，所以自相关系数里面实际掺杂了其他变量对 $y_t$ 与 $y_{t-k}$ 的影响；
  • 为了能单纯测度 $y_{t-k}$ 对 $y_t$ 的影响，引进偏自相关系数的概念。对于平稳时间序列 { y t } \{y_t\} {yt​} ，所谓滞后 $k$ 偏自相关系数指在剔除了中间 $k-1$ 个随机变量 $y_{t-1}、y_{t-2}、\dots 、y_{t-k+1}$ 的干扰之后，$x(t-k)$对 { y t } \{y_t\} {yt​}影响的相关程度；
  • 公式：$$PACF(k)=\frac{Cov[(y_t-{\bar{y}}_t),(y_{t-k}-{\bar{y}}_{t-k})]}{\sqrt{Var(y_t-{\bar{y}}_t)}\sqrt{Var(y_{t-k}-{\bar{y}}_{t-k})}}$$

• 拖尾、截尾判断法：拖尾指序列以指数率单调递减或震荡衰减，而截尾指序列从某个时点变得非常小；

  • 拖尾（出现以下情况，通常视为（偏）自相关系数拖尾）：

    • 如果有超过 5% 的样本（偏）自相关系数都落入 2 倍标准差范围之外；
    • 或者是由显著非 0 的（偏）自相关系数衰减为小值波动的过程比较缓慢或非常连续；
    • 

  • 截尾（出现以下情况，通常视为（偏）自相关系数 d 阶截尾）：

    • 在最初的d阶明显大干 2 倍标准差范围；
    • 之后几乎 95% 的（偏）自相关系数都落在 2 倍标准差范围以内；
    • 且由非零自相关系数衰减为在零附近小值波动的过程非常突然。
    • 

  • 

(2) AIC 准则

通过拖尾和截尾对模型定阶，具有很强的主观性。对于模型参数估计的方法，是通过对损失和正则项的加权评估。在参数选择的时候，需要平衡预测误差与模型复杂度。可以根据信息准则函数法，来确定模型的阶数。
AIC 准则全称是最小化信息量准则（AIC，Akaike Information Criterion）： $$AIC=-2ln(L)+2K$$其中 $L$ 表示模型的极大似然函数，$K$ 表示模型参数个数。
AIC准则存在一定的不足。当样本容量很大时，在 AIC 准则中拟合误差提供的信息就要受到样本容量的放大，而参数个数的惩罚因子却和样本容量没关系（一直是2），因此当样本容量很大时使用 AIC 准则的模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多。

(3) BIC 准则

贝叶斯信息准则（BIC，Bayesian Information Criterion）贝叶斯信息准则弥补了 AIC 的不足：$$BIC=-2ln(L)+Kln(n)$$其中 $n$ 表示样本容量。
AIC、BIC 这两个评价指标越小越好。可以通过网格搜索，确定 AIC、BIC 最优的模型(p、q)。

3.MATLAB 实现

%% ARIMA时间序列预测（BIC自动定阶 + 训练集验证）
clc; clear; close all;%% 加载示例数据（航空乘客数据）
% 数据说明：1949-1960年每月航空乘客数（单位：千人）
data = [112 118 132 129 121 135 148 148 136 119 104 118 ...115 126 141 135 125 149 170 170 158 133 114 140 ...145 150 178 163 172 178 199 199 184 162 146 166 ...171 180 193 181 183 218 230 242 209 191 172 194 ...196 196 236 235 229 243 264 272 237 211 180 201 ...204 188 235 227 234 264 302 293 259 229 203 229 ...242 233 267 269 270 315 364 347 312 274 237 278 ...284 277 317 313 318 374 413 405 355 306 271 306 ...315 301 356 348 355 422 465 467 404 347 305 336 ...340 318 362 348 363 435 491 505 404 359 310 337 ...360 342 406 396 420 472 548 559 463 407 362 405 ...417 391 419 461 472 535 622 606 508 461 390 432]';
dates = datetime(1949,1:length(data),1, 'Format', 'MMM yyyy');%% 划分训练集和测试集（前80%训练，后20%测试）
trainRatio = 0.8;
n = length(data);
trainEnd = floor(n * trainRatio);
trainData = data(1:trainEnd);
testData = data(trainEnd+1:end);
trainDates = dates(1:trainEnd);
testDates = dates(trainEnd+1:end);% 可视化数据分割
figure('Color','white', 'Position',[100 100 1000 400]);
plot(trainDates, trainData, 'b', testDates, testData, 'm', 'LineWidth',1.5);
title('数据划分：训练集(蓝色) vs 测试集(品红)');
xlabel('日期'); ylabel('乘客数（千）'); 
legend('训练集', '测试集', 'Location','northwest'); grid on;%% 步骤1：数据平稳化处理
% ADF检验确定差分阶数d
d = 0;
[h, pValue] = adftest(trainData);
while pValue > 0.05 && d < 3d = d + 1;trainDataDiff = diff(trainData, d);[h, pValue] = adftest(trainDataDiff);
end
fprintf('自动确定差分阶数 d = %d\n', d);%% 步骤2：BIC准则自动选择最优p,q
maxP = 13;  % AR最大尝试阶数
maxQ = 13;  % MA最大尝试阶数
numModels = (maxP+1)*(maxQ+1); % 总模型数
bicMatrix = inf(maxP+1, maxQ+1); % 初始化BIC矩阵% 遍历所有p,q组合
for p = 0:maxPfor q = 0:maxQtry% 创建ARIMA模型model = arima(p, d, q);% 拟合模型（抑制命令行输出）[fitModel, ~, logL] = estimate(model, trainData, 'Display','off');% 计算参数总数（AR + MA + 常数项 + 方差）numParams = p + q + 1 + 1; % 计算BIC值：BIC = -2*logL + numParams*log(n)n = length(trainData) - d; % 有效样本量bic = -2*logL + numParams*log(n);bicMatrix(p+1, q+1) = bic;catch ME% 跳过无法收敛的模型fprintf('[警告] p=%d, q=%d 模型拟合失败: %s\n', p, q, ME.message);endend
end% 找到最小BIC对应的p,q
[minBIC, idx] = min(bicMatrix(:));
[bestP, bestQ] = ind2sub(size(bicMatrix), idx);
bestP = bestP - 1; % 调整索引偏移
bestQ = bestQ - 1;
fprintf('\n===== BIC自动定阶结果 =====\n');
fprintf('最优 p = %d, q = %d (BIC=%.2f)\n', bestP, bestQ, minBIC);% 可视化BIC热力图
figure('Color','white');
heatmap(0:maxQ, 0:maxP, bicMatrix, 'Colormap', parula);
xlabel('MA阶数 q'); ylabel('AR阶数 p'); 
title('BIC值热力图（颜色越深BIC越小）');%% 步骤3：使用最优参数建立ARIMA模型
model = arima(bestP, d, bestQ);
fitModel = estimate(model, trainData); % 训练最终模型
disp('最优ARIMA模型参数：');
disp(fitModel);%% 步骤4：模型诊断（残差检验）
residuals = infer(fitModel, trainData);figure('Color','white', 'Position',[100 100 1000 600]);
subplot(2,2,1); plot(residuals); title('残差序列');
subplot(2,2,2); autocorr(residuals); title('残差ACF');
subplot(2,2,3); parcorr(residuals); title('残差PACF');
subplot(2,2,4); qqplot(residuals); title('Q-Q图');% Ljung-Box检验
[h, p] = lbqtest(residuals, 'Lags', 20);
fprintf('残差白噪声检验p值 = %.4f\n', p);%% 步骤5：测试集预测与评估
numTest = length(testData);
[forecastData, forecastMSE] = forecast(fitModel, numTest, trainData);% 计算误差指标
actual = testData;
predicted = forecastData;
mae = mean(abs(actual - predicted));
rmse = sqrt(mean((actual - predicted).^2));
mape = mean(abs((actual - predicted)./actual)) * 100;
% 可视化预测对比
figure('Color','white', 'Position',[100 100 1200 500]);
hold on;
plot(trainDates, trainData, 'b', 'LineWidth',1.5);
plot(testDates, testData, 'm', testDates, forecastData, 'r', 'LineWidth',1.5);
plot(testDates, forecastData-1.96*sqrt(forecastMSE), 'k--');
plot(testDates, forecastData+1.96*sqrt(forecastMSE), 'k--');
text(testDates(end), forecastData(end)+30, ...sprintf('MAE=%.1f\nRMSE=%.1f', mae, rmse), ...'VerticalAlignment','bottom', 'HorizontalAlignment','right');
title(sprintf('ARIMA(%d,%d,%d) 测试集预测', bestP, d, bestQ));
xlabel('日期'); ylabel('乘客数（千）'); 
legend('训练集', '测试集实际值', '预测值', '95%置信区间', 'Location','northwest');
grid on; xlim([dates(1) dates(end)]);
fprintf('\n===== 测试集预测效果 =====\n');
fprintf('MAE: %.1f\nRMSE: %.1f\nMAPE: %.1f%%\n', mae, rmse, mape);

参考资料

[1] 终于弄明白ARIMA模型啦！包括确定pq值详细解释！！哔哩哔哩bilibili
[2] 时间序列模型（ARIMA）讲解(附matlab和python代码) 【数学建模快速入门】数模加油站 江北_哔哩哔哩_bilibili
[3] 「SPSSAU｜数据分析」：时间序列模型分析分析
[4] 5-参数选择_哔哩哔哩_bilibili