时间序列预测实战：ARMA、灰色预测与多元回归在臭氧消耗建模中的应用

1. 项目概述：从一道赛题到预测方法论的实战复盘

2016年第五届数学建模国际赛（小美赛）的A题“臭氧消耗预测”，对于当年参赛的我们来说，不仅仅是一道题目，更像是一次将多种经典预测模型置于真实环境数据下进行“同台竞技”的绝佳机会。这道题的核心，是要求参赛者基于给定的历史臭氧层消耗相关数据，构建数学模型，对未来趋势进行预测。这听起来像是典型的时间序列预测问题，但当你真正深入数据，会发现它融合了环境科学、统计学和计算编程的多重挑战。如今回头看，解题过程远不止是提交一篇论文和几行代码，它完整地呈现了从问题理解、数据勘探、模型选型、对比验证到结果分析的全链条建模思维。无论是当时主流的ARMA模型、灰色预测GM(1,1)，还是基础的多元线性回归，每一种方法的选择与调整背后，都有一套严密的逻辑和大量“踩坑”得来的经验。本文将基于当年的解题全流程文档，结合后续多年的建模与数据分析心得，为你拆解这道题的完整解决路径，并附上可复现的MATLAB程序核心代码与深度解析。无论你是正在备战数模的新手，还是希望巩固预测方法的数据分析从业者，这篇复盘都将提供从理论到实践的扎实参考。

2. 赛题核心与数据特征解析

2.1 问题重述与目标拆解

原题通常提供一段关于臭氧层消耗现象的背景描述，以及一份包含多个变量的时间序列数据集。变量可能包括年度或月度数据，例如：特定区域的臭氧柱总量、消耗性气体（如CFC-11, CFC-12）的浓度、太阳活动指数、大气温度等。题目的核心要求可以拆解为以下几个层次：

1. 趋势分析与预测：利用历史数据，建立臭氧消耗量（或相关指标）的预测模型，并给出未来若干时间点的预测值。
2. 关键因素识别：分析并量化不同因素（如各类消耗性气体浓度、环境因子）对臭氧消耗的影响程度。
3. 模型评估与比较：可能需要尝试多种预测模型，并比较其性能，说明各自的优缺点及适用条件。
4. 政策建议（通常为拓展部分）：基于预测结果，提出有针对性的环境保护建议。

我们的核心任务是完成前三点，形成一个闭环：理解数据 -> 选择并建立模型 -> 评估模型 -> 输出预测。

2.2 数据预处理与探索性分析（EDA）

拿到数据后的第一步绝不是直接套模型，而是花时间“认识”你的数据。这一步直接决定了后续模型的有效性。

2.2.1 数据清洗与缺失值处理原始数据可能存在记录错误、异常值或缺失值。对于时间序列数据，常见的处理方法包括：

• 线性插值：适用于缺失较少且数据趋势平缓的情况。
• 前向填充（Forward Fill）或后向填充（Backward Fill）：适用于连续性较强的监测数据。
• 剔除：如果缺失数据占比较小，且位于序列起始或末尾，可以考虑直接剔除该记录。

注意：处理缺失值的方法需要记录在论文中，并简要说明理由。粗暴地删除或随意插值都可能引入偏差。

2.2.2 平稳性检验与变换这是时间序列分析（如ARMA）的关键前提。我们使用MATLAB的adftest(Augmented Dickey-Fuller test) 来检验序列的平稳性。

% 假设 ozone_data 是臭氧浓度的时序向量 [h, pValue, stat, cValue] = adftest(ozone_data, 'model', 'TS'); if h == 0 disp('序列非平稳，需要进行差分处理。'); ozone_data_diff = diff(ozone_data); % 一阶差分 % 再次检验差分后序列的平稳性 [h_diff, ~] = adftest(ozone_data_diff); end

如果序列非平稳（h=0），通常需要进行差分运算，直到得到一个平稳序列。对于有明显趋势或季节性的数据，可能需要多次差分或进行对数变换等。

2.2.3 可视化分析利用MATLAB绘制时序图、自相关图（ACF）和偏自相关图（PACF），直观感受数据的趋势、周期性和模型初步定阶。

figure; subplot(2,2,1); plot(ozone_data, 'b-', 'LineWidth', 1.5); title('臭氧浓度原始时序图'); xlabel('时间'); ylabel('浓度'); grid on; subplot(2,2,2); autocorr(ozone_data); % 绘制自相关图 title('原始序列ACF'); subplot(2,2,3); plot(ozone_data_diff, 'r-'); title('一阶差分后时序图'); grid on; subplot(2,2,4); parcorr(ozone_data_diff); % 绘制偏自相关图 title('差分后序列PACF');

通过看图，我们可以初步判断：原始序列是否有上升/下降趋势（需差分）？ACF是否拖尾（可能为AR模型）？PACF是否截尾（可能为MA模型）？这为后续ARMA模型的定阶（p, q值）提供重要依据。

3. 三大预测模型的原理、实现与对比

针对本题，我们重点部署了三种经典模型：适用于单变量时间序列的ARMA和灰色预测，以及能处理多变量影响的多元回归。

3.1 ARMA模型：捕捉序列的内在记忆

ARMA（自回归移动平均）模型是处理平稳时间序列的利器。其思想是当前值由过去若干期的值（AR部分）和过去若干期的误差（MA部分）共同决定。

3.1.1 模型定阶与建立在确保序列平稳后，我们需要确定AR阶数p和MA阶数q。除了观察ACF/PACF图，更可靠的方法是结合信息准则（如AIC, BIC）进行网格搜索。

% 假设平稳序列为 stationary_data max_p = 5; % 预设最大AR阶数 max_q = 5; % 预设最大MA阶数 logL = zeros(max_p+1, max_q+1); % 对数似然值 numParams = zeros(max_p+1, max_q+1); % 参数个数 for p = 0:max_p for q = 0:max_q if p==0 && q==0 continue; % 跳过ARMA(0,0) end try mdl = arima(p, 0, q); % 创建ARMA(p,q)模型，d=0因已平稳 [estMdl, ~, logL(p+1, q+1)] = estimate(mdl, stationary_data, 'Display', 'off'); numParams(p+1, q+1) = p + q + 1; % +1为常数项参数 catch logL(p+1, q+1) = -Inf; % 模型拟合失败 end end end % 计算AIC和BIC，值越小越好 aic = 2*numParams - 2*logL; bic = log(length(stationary_data))*numParams - 2*logL; % 找到AIC/BIC最小的(p,q)组合 [minAIC, idxAIC] = min(aic(:)); [minBIC, idxBIC] = min(bic(:)); [p_aic, q_aic] = ind2sub(size(aic), idxAIC); [p_bic, q_bic] = ind2sub(size(bic), idxBIC); p_aic = p_aic - 1; q_aic = q_aic - 1; % 调整索引 fprintf('AIC推荐阶数: ARMA(%d, %d)\n', p_aic, q_aic); fprintf('BIC推荐阶数: ARMA(%d, %d)\n', p_bic, q_bic);

在实际操作中，AIC和BIC推荐的阶数可能不同。BIC对参数惩罚更重，倾向于选择更简单的模型。我们通常以BIC为准，兼顾模型的简洁性与预测能力。

3.1.2 模型诊断与预测模型建立后，必须进行残差诊断，检验残差是否为白噪声（均值为0、方差恒定、无自相关）。

% 使用选定的(p, q)拟合模型 bestMdl = arima(p_bic, 0, q_bic); estMdl = estimate(bestMdl, stationary_data); % 残差诊断 res = infer(estMdl, stationary_data); % 获取残差 figure; subplot(2,2,1); plot(res); title('残差序列图'); subplot(2,2,2); histogram(res, 20); title('残差直方图'); subplot(2,2,3); autocorr(res); title('残差ACF图'); subplot(2,2,4); parcorr(res); title('残差PACF图'); % 进行Ljung-Box检验，原假设为残差是白噪声 [h, pValue] = lbqtest(res, 'Lags', [10, 15]);

如果残差通过白噪声检验（h=0），说明模型已充分提取了序列信息。随后可以进行预测：

numSteps = 10; % 预测未来10期 [yF, yMSE] = forecast(estMdl, numSteps, 'Y0', stationary_data); % yF为预测值，yMSE为预测均方误差 lowerBound = yF - 1.96*sqrt(yMSE); % 95%置信区间下限 upperBound = yF + 1.96*sqrt(yMSE); % 95%置信区间上限

实操心得：ARMA模型对序列的平稳性要求严格。如果原始序列有很强的趋势或季节性，仅靠差分可能不够，可能需要考虑ARIMA（加入差分项）或SARIMA（加入季节性项）。对于臭氧数据，其长期趋势（如受政策影响下降）和可能的年度周期都需要仔细处理。

3.2 灰色预测GM(1,1)：小样本、贫信息的利器

当数据量较少（通常少于20个），且序列呈现近似指数增长或衰减趋势时，灰色预测模型往往能发挥奇效。它通过累加生成（AGO）弱化原始序列的随机性，挖掘其内在规律。

3.2.1 模型建立过程GM(1,1)是灰色预测中最基础的模型，其建模步骤如下：

1. 原始序列：X(0) = [x(0)(1), x(0)(2), ..., x(0)(n)]
2. 一次累加生成（1-AGO）：X(1)(k) = sum_{i=1}^{k} X(0)(i)，得到新序列X(1)。
3. 建立灰微分方程：dX(1)/dt + aX(1) = u，其中a为发展系数，u为灰色作用量。
4. 求解参数：利用最小二乘法估计参数a和u。
5. 得到时间响应式（预测公式）：X^(1)(k+1) = [X(0)(1) - u/a] * exp(-a*k) + u/a。
6. 累减还原：X^(0)(k+1) = X^(1)(k+1) - X^(1)(k)，得到原始序列的预测值。

3.2.2 MATLAB实现代码

function [predict, a, u] = gm11(data, predict_step) % data: 原始行向量，如 [data1, data2, ..., datan] % predict_step: 预测步长 n = length(data); % 1. 累加生成 X1 = cumsum(data); % 2. 构造数据矩阵B和常数向量Y B = [-0.5*(X1(1:end-1)+X1(2:end))', ones(n-1,1)]; Y = data(2:end)'; % 3. 最小二乘求解参数 a, u parameters = (B'*B) \ (B'*Y); a = parameters(1); u = parameters(2); % 4. 计算预测值（累加序列） X1_predict = zeros(1, n+predict_step); X1_predict(1) = data(1); for k = 1:(n+predict_step-1) X1_predict(k+1) = (data(1) - u/a) * exp(-a*k) + u/a; end % 5. 累减还原，得到原始序列预测值 predict = [data(1), X1_predict(2:end) - X1_predict(1:end-1)]; predict = predict(n+1:end); % 只返回未来的预测值 end % 调用示例 ozone_historical = [数据序列]; % 替换为实际数据 future_steps = 5; [pred_vals, a_coef, u_coef] = gm11(ozone_historical, future_steps); fprintf('发展系数 a = %.4f，灰色作用量 u = %.4f\n', a_coef, u_coef); disp('未来预测值：'); disp(pred_vals);

3.2.3 模型检验灰色预测模型必须进行精度检验，常用方法有：

• 后验差检验：计算后验差比值C和小误差概率P。
  • C = 原始序列标准差 S1 / 残差标准差 S2。C越小越好（<0.35优秀，<0.5合格，<0.65勉强合格）。
  • P = 概率 {|残差 - 残差均值| < 0.6745*S1}。P越大越好（>0.95优秀，>0.8合格）。
• 相对误差检验：计算模型拟合值与历史真实值的相对误差。

注意事项：GM(1,1)默认适用于具有指数趋势的序列。如果原始序列非常平缓或波动剧烈，预测效果可能不佳。此外，灰色预测是“滚动的”，用最新数据重新建模预测下一步通常比一次性预测多步更准。对于臭氧数据，如果其下降趋势符合近似指数衰减，GM(1,1)会是一个简洁有效的选择。

3.3 多元线性回归：量化多因素影响

如果题目提供了除时间外的其他变量（如各类气体浓度、温度），那么多元线性回归可以帮助我们量化这些因素对臭氧消耗的具体影响，并基于影响因素的变化进行预测。

3.3.1 模型构建与变量筛选模型形式为：Ozone = β0 + β1*X1 + β2*X2 + ... + βk*Xk + ε。 在MATLAB中，可以使用fitlm函数。

% 假设数据表 T 包含变量：Ozone, CFC11, CFC12, SolarIndex, Year % Year 可能作为控制变量或用于捕捉线性趋势 mdl_lm = fitlm(T, 'Ozone ~ CFC11 + CFC12 + SolarIndex + Year'); disp(mdl_lm); % 显示详细的回归结果摘要

回归结果摘要会显示每个系数的估计值、标准误差、t统计量和p值。p值用于判断该变量是否对因变量有显著影响（通常以p<0.05为显著）。

3.3.2 模型诊断与共线性处理建立回归模型后，必须进行诊断：

1. 残差分析：检查残差是否满足独立性、正态性、同方差性。绘制残差图。

   figure; subplot(2,2,1); plotResiduals(mdl_lm, 'fitted'); % 残差与拟合值图 subplot(2,2,2); plotResiduals(mdl_lm, 'probability'); % 正态概率图 subplot(2,2,3); plotResiduals(mdl_lm, 'lagged'); % 残差与滞后残差图（查自相关）

2. 多重共线性诊断：如果自变量之间高度相关，会导致系数估计不稳定。计算方差膨胀因子（VIF）。

   vif = diag(inv(corrcoef(table2array(T(:, {'CFC11', 'CFC12', 'SolarIndex'}))))); % 计算除截距项外自变量的VIF disp('VIF值:'); disp(vif);

   VIF > 10 通常认为存在严重共线性。解决方法包括剔除相关性高的变量、使用主成分回归（PCR）或岭回归（Ridge Regression）。

3.3.3 预测与置信区间使用训练好的模型对新自变量数据进行预测。

% 创建新观测值表格 newData = table([CFC11_new], [CFC12_new], [SolarIndex_new], [Year_new], ... 'VariableNames', {'CFC11', 'CFC12', 'SolarIndex', 'Year'}); [pred_y, pred_ci] = predict(mdl_lm, newData); % pred_ci为预测区间

实操心得：多元回归的强大之处在于可解释性。你可以明确说出“CFC-11浓度每增加1单位，臭氧浓度平均减少β1单位”。但它的预测精度严重依赖于自变量的未来值是否已知或可准确预测。在本题中，如果需要预测未来臭氧，你必须先有或先预测出未来CFC浓度等变量的值，这构成了一个“嵌套预测”问题，增加了不确定性。

4. 模型比较、评估与综合策略

单一模型往往有局限性，在实际解题中，我们通常会运行多个模型，并进行系统比较。

4.1 评估指标选择

我们使用以下指标在历史数据上（如留出最后几年的数据作为验证集）评估模型：

• 均方根误差（RMSE）：衡量预测值与真实值之间的偏差，对较大误差更敏感。

  rmse = sqrt(mean((y_true - y_pred).^2));

• 平均绝对百分比误差（MAPE）：相对误差，易于理解。

  mape = mean(abs((y_true - y_pred) ./ y_true)) * 100;

• 决定系数（R²）：反映模型对数据波动的解释能力，越接近1越好。

  ss_res = sum((y_true - y_pred).^2); ss_tot = sum((y_true - mean(y_true)).^2); r2 = 1 - (ss_res / ss_tot);

4.2 对比结果分析与模型融合

将ARMA、GM(1,1)和多元回归在验证集上的表现填入下表进行对比：

基于对比，我们可能采取以下策略：

1. 择优选择：选择一个在验证集上RMSE和MAPE最小、R²最高的模型作为最终预测模型。
2. 组合预测：如果各模型各有优劣，可以采用加权平均的方式组合预测结果。权重可以根据各模型在验证集上的误差倒数或其他方法确定。

   % 假设有三个模型的预测结果：pred1, pred2, pred3 % 计算在验证集上的权重（例如，根据RMSE的倒数） w1 = 1/rmse1; w2 = 1/rmse2; w3 = 1/rmse3; total_w = w1 + w2 + w3; final_pred = (w1/total_w)*pred1 + (w2/total_w)*pred2 + (w3/total_w)*pred3;

3. 情景分析：提交多个模型的预测结果，并说明各自适用的条件和假设。例如，“在假设未来消耗性气体浓度保持当前下降趋势的前提下，采用ARIMA模型预测...；若考虑极端太阳活动事件，则基于多元回归模型的情景分析如下...”。

5. 完整解题流程复盘与避坑指南

回顾整个解题过程，以下几个关键环节最容易出问题，也是新手需要特别注意的地方：

5.1 数据预处理阶段的坑

• 忽视平稳性检验直接上ARMA：这是最常见的错误。非平稳序列拟合ARMA会导致伪回归，预测毫无意义。务必先画图、做ADF检验，必要时进行差分或变换。
• 对缺失值处理不当：时间序列的缺失值处理需要谨慎。简单用全局均值填充会破坏时间依赖性。优先使用时序方法（如插值、前向填充）或基于模型的插补。
• 未考虑季节性：臭氧数据可能具有年度周期性。如果ACF图在滞后12（月度数据）或1（年度数据）的倍数处出现峰值，提示存在季节性。这时应考虑SARIMA模型或引入季节性虚拟变量。

5.2 模型建立与诊断阶段的坑

• ARMA模型定阶过度依赖自动函数：MATLAB的auto.arima（需Econometrics Toolbox）或相关自动定阶函数很方便，但不能完全替代人工判断。一定要结合ACF/PACF图和信息准则，并检查残差。有时，一个更简洁的模型（如AR(1)）可能比复杂的ARMA(p,q)预测效果更稳健。
• 灰色预测不进行精度检验：建完GM(1,1)模型后，务必计算后验差比C和小误差概率P。如果检验不合格（如C>0.65, P<0.7），说明该模型不适用于当前数据，预测结果不可信。
• 多元回归忽视假设检验：回归不是把变量扔进去得到系数就完了。必须检查残差是否独立、正态、同方差，必须检验多重共线性。如果残差存在自相关（时间序列数据常见），则需要改用时间序列回归模型（如加入滞后项）或广义最小二乘法。

5.3 预测与报告阶段的坑

• 混淆预测区间与置信区间：forecast或predict函数给出的区间通常是预测区间，它比置信区间（对条件均值的估计区间）更宽，因为它包含了模型误差和观测误差。在报告中要表述清楚。
• 外推风险：所有模型都是基于历史数据建立的。用它们预测未来，本质上是假设“历史规律在未来持续”。对于臭氧消耗这种受国际公约（如《蒙特利尔议定书》）强烈影响的过程，政策突变点前后的数据规律可能完全不同。在报告中必须强调这一外推风险，并进行讨论。
• 只谈结果，不谈不确定性：优秀的数模论文不仅要给出预测值，还要给出预测的不确定性范围（如95%预测区间），并讨论哪些因素（如数据误差、模型假设、未来情景）会影响预测的准确性。

5.4 MATLAB编程实操技巧

• 代码模块化与注释：将数据读取、预处理、模型拟合、评估分别写成函数或独立脚本节，并添加清晰注释。这便于调试和报告复现。
• 保存中间结果与图形：使用save命令保存工作区变量，使用saveas或exportgraphics高分辨率保存生成的图表。避免每次修改代码都要重新运行耗时长的部分。
• 利用循环进行批量操作：比如需要尝试多组模型参数时，使用for或parfor（并行循环）可以极大提高效率。

  % 示例：批量测试不同训练集长度对预测误差的影响 train_ratios = 0.6:0.05:0.9; % 训练集比例 rmse_results = zeros(size(train_ratios)); for i = 1:length(train_ratios) ratio = train_ratios(i); split_idx = floor(ratio * length(data)); train_data = data(1:split_idx); test_data = data(split_idx+1:end); % ... 在此训练模型并预测测试集 ... rmse_results(i) = calculateRMSE(test_pred, test_data); end figure; plot(train_ratios, rmse_results, 'o-'); xlabel('训练集比例'); ylabel('RMSE');

这道“臭氧消耗预测”赛题，本质上是一个经典的时间序列预测案例。它教会我们的，不是某个特定模型的用法，而是一套面对预测问题时的标准化建模流程：从数据洞察出发，严谨地选择、建立、诊断、比较模型，最后审慎地解读和报告预测结果。这个过程里用到的工具（平稳性检验、ACF/PACF、AIC/BIC、残差诊断、VIF、RMSE/MAPE）和思维框架，适用于绝大多数预测场景。最后分享一个个人体会：在数模竞赛和实际工作中，没有“最好”的模型，只有“最合适”的模型。模型的复杂性应与数据量和问题背景相匹配。有时，一个精心构建的简单模型，其可靠性和可解释性远胜于一个黑箱般的复杂模型。