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