时间序列预测实战：ARIMA、灰色预测与多元回归模型在MATLAB中的实现与比较

1. 项目概述：从一道赛题看预测模型的实战选择

2016年第五届数学建模国际赛（小美赛）的A题“臭氧消耗预测”，对于很多初次接触数学建模，尤其是时间序列预测的同学来说，是一道非常经典的入门与进阶练习题。这道题的核心并不在于提出一个惊世骇俗的全新理论，而在于如何将经典的预测模型，如ARMA、灰色预测、多元回归等，在一个具体的、有现实背景的问题中进行合理的应用、比较与融合。很多人在拿到题目后，第一反应是去搜索“标准答案”或“最优模型”，但建模的魅力恰恰在于没有唯一的“标准答案”，只有基于数据特征、问题背景和模型假设的“更优解”。今天，我就结合自己多年指导建模和数据分析的经验，来深度拆解这道题的解题全流程，重点不是给出一个死的程序，而是分享一套活的、可复用的建模思维与实操心法。

这道题通常会提供一段时期内（比如若干年）的月度或年度臭氧消耗相关数据，可能包括臭氧层空洞面积、特定区域臭氧浓度、以及一些潜在的影响因子，如氟氯烃（CFCs）排放量、太阳活动指数、大气环流指标等。我们的核心任务就是利用历史数据，构建数学模型，对未来一段时间（比如未来5年或10年）的臭氧消耗趋势进行预测。这本质上是一个时间序列预测问题，但可能融合了影响因素分析。关键词中的ARMA（自回归移动平均模型）、灰色预测、多元回归，正是解决此类问题的三把经典钥匙，而MATLAB则是挥舞这些钥匙最顺手的工具之一。接下来，我将从解题思路、模型详解、MATLAB实操到结果分析，一步步展开。

2. 解题核心思路与模型选型逻辑

面对一个预测问题，尤其是赛题，最忌讳的就是不加思考地套用模型。一个成熟的建模者，第一步永远是“读懂数据”和“定义问题”。

2.1 问题分析与数据预处理

首先，我们需要明确“臭氧消耗”这个指标是什么。是南极臭氧空洞的面积最大值？是全球平均臭氧柱总量？题目给出的具体指标决定了我们后续处理数据的粒度。假设我们拿到的是1990年至2015年南极臭氧空洞的年最大面积数据。

第一步：数据可视化与初步诊断。在MATLAB里，第一件事绝不是跑模型，而是画图。

% 假设数据已加载，year为年份向量，area为臭氧空洞面积向量 figure; plot(year, area, ‘o-‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 8); xlabel(‘年份‘); ylabel(‘臭氧空洞面积 (百万平方公里)‘); title(‘南极臭氧空洞年最大面积变化趋势‘); grid on;

这张图能告诉我们很多信息：趋势是上升还是下降？是否存在周期性（虽然年数据周期不明显）？有没有明显的异常点？数据是否平稳？通过看图，我们可能发现1990年代面积快速增长，21世纪初达到峰值后，由于《蒙特利尔议定书》生效，CFCs被限制，面积增长放缓甚至出现下降趋势。这提示数据可能具有趋势性和外部政策干预点。

第二步：数据平稳化处理。很多时间序列模型（如ARMA）要求数据是平稳的，即均值和方差不随时间变化。如果我们的数据有明显趋势，就需要进行差分。

% 计算一阶差分，消除线性趋势 diff_area = diff(area); % 画出差分后的序列图，观察是否平稳 figure; plot(year(2:end), diff_area, ‘s-‘); title(‘一阶差分后序列‘);

如果一阶差分后序列围绕0值波动，无明显趋势，则可认为基本平稳。此外，还需要进行**单位根检验（如ADF检验）**来统计上确认平稳性。在MATLAB中，可以使用Econometrics Toolbox中的adftest函数。

第三步：影响因素考量。如果题目除了时间序列数据，还提供了诸如CFCs浓度、太阳黑子数等协变量，那么问题就从单纯的时间序列预测，转变为时间序列回归或多元分析问题。我们需要分析这些因素与臭氧消耗的相关性。

2.2 三大预测模型的核心思想与适用场景

为什么是ARMA、灰色预测和多元回归？因为它们代表了三种不同的预测哲学，适用于不同的数据特征。

1. ARMA模型：挖掘数据内部规律

   • 核心思想：认为当前值可以由其过去若干期的历史值（自回归AR部分）以及过去若干期的随机扰动（移动平均MA部分）的线性组合来解释。它不关心外部因素，只专注于序列自身的内在结构。
   • 适用场景：适用于平稳的、具有短期相关性的时间序列。对于臭氧数据，如果我们通过差分消除了趋势，得到的平稳序列就可以尝试用ARMA或其扩展模型ARIMA（自回归积分移动平均，I指差分）来建模。
   • 关键步骤：模型识别（确定AR阶数p和MA阶数q）→ 参数估计 → 模型检验 → 预测。

2. 灰色预测GM(1,1)：小样本、趋势性数据的利器

   • 核心思想：针对“部分信息已知，部分信息未知”的“灰色系统”。它通过对原始数据进行一次累加生成（1-AGO），弱化随机性，凸显指数增长趋势，然后建立一阶微分方程进行预测，最后再累减还原得到预测值。
   • 适用场景：数据量少（通常只需4个以上数据点）、具有明显指数增长或衰减趋势的序列。臭氧消耗在特定历史阶段可能呈现近似指数变化，适合用灰色预测。但它对波动大的数据预测效果较差。
   • 关键步骤：原始数据累加 → 构建灰微分方程 → 求解发展系数和灰色作用量 → 时间响应函数预测 → 累减还原。

3. 多元线性回归：建立因果驱动模型

   • 核心思想：假设因变量（臭氧消耗）与多个自变量（CFCs排放、太阳活动等）之间存在线性关系。通过历史数据拟合出一个线性方程，然后代入自变量的未来值（或预测值）来预测因变量。
   • 适用场景：当影响因子明确、且能获取其历史及未来（或可假设）数据时。这需要较强的先验知识或题目提供协变量数据。它从“因果”角度进行预测，解释性强。
   • 关键步骤：变量选择与共线性诊断 → 建立回归方程 → 显著性检验（R², F检验, t检验） → 利用自变量进行预测。

模型选择心法：在实际解题中，不要只用一个模型。一个稳健的策略是：先用灰色预测做快速趋势捕捉（尤其数据少时），再用ARIMA对去趋势后的细节进行建模，如果有多元数据，则构建回归模型。最后，可以比较不同模型的预测结果，甚至进行组合预测（如加权平均），往往能提升预测的稳健性和精度。这正是数学建模竞赛中体现分析深度的地方。

3. 模型实现与MATLAB实操详解

理论说得再多，不如一行代码。下面我将分别展示三大模型在MATLAB中的核心实现步骤和代码片段，并穿插大量实操中才会遇到的“坑”和技巧。

3.1 ARIMA模型建模全流程

ARIMA模型是ARMA模型对非平稳序列的扩展，记作ARIMA(p,d,q)，其中d是差分阶数。我们以消除趋势后的序列为例。

步骤1：平稳性检验与差分阶数d确定

% 假设原始序列为 area，年份为 year % 1. 画图观察 figure; subplot(2,1,1); plot(year, area); title(‘原始序列‘); subplot(2,1,2); plot(year(2:end), diff(area)); title(‘一阶差分序列‘); % 2. 单位根检验 (需Econometrics Toolbox) [h1, pValue1] = adftest(area); % 检验原始序列 [h2, pValue2] = adftest(diff(area)); % 检验一阶差分序列 % h=0表示不拒绝原假设（存在单位根，非平稳），h=1表示拒绝（平稳） fprintf(‘原始序列ADF检验p值: %.4f，是否平稳: %d\n‘, pValue1, h1); fprintf(‘一阶差分序列ADF检验p值: %.4f，是否平稳: %d\n‘, pValue2, h2);

如果一阶差分后序列平稳，则设置d=1。

步骤2：模型识别与定阶(p, q)定阶是ARIMA建模的难点和核心。常用方法是观察**自相关函数(ACF)和偏自相关函数(PACF)**图。

data_stationary = diff(area); % 使用平稳化后的数据 figure; subplot(2,1,1); autocorr(data_stationary, ‘NumLags‘, 20); % ACF图 subplot(2,1,2); parcorr(data_stationary, ‘NumLags‘, 20); % PACF图

• ACF拖尾，PACFp阶截尾-> AR(p)模型。
• ACFq阶截尾，PACF拖尾-> MA(q)模型。
• 两者都拖尾-> ARMA(p,q)模型。

此外，可以使用AIC/BIC信息准则进行自动定阶，遍历一组(p,q)值，选择AIC/BIC最小的模型。

max_p = 3; % 假设最大尝试AR阶数 max_q = 3; % 假设最大尝试MA阶数 d = 1; % 差分阶数 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; % 跳过ARIMA(0,1,0)，即纯随机游走 end mdl = arima(p, d, q); [~,~,logL(p+1, q+1)] = estimate(mdl, area, ‘Display‘, ‘off‘); numParams(p+1, q+1) = p + q; % 不包括常数项和方差 end end % 计算AIC: AIC = -2*logL + 2*k, k为参数总数（包括常数项和方差，这里简化用numParams+2） AIC = -2*logL + 2*(numParams + 2); % 找到最小AIC对应的(p,q) [minAIC, idx] = min(AIC(:)); [p_opt, q_opt] = ind2sub(size(AIC), idx); p_opt = p_opt - 1; q_opt = q_opt - 1; fprintf(‘AIC准则推荐模型: ARIMA(%d, %d, %d)\n‘, p_opt, d, q_opt);

步骤3：模型估计与检验

% 使用推荐的阶数拟合模型 Mdl = arima(p_opt, d, q_opt); EstMdl = estimate(Mdl, area, ‘Display‘, ‘full‘); % 详细输出估计结果 % 残差检验：好的模型残差应近似为白噪声（无自相关） res = infer(EstMdl, area); % 获取残差 figure; subplot(2,2,1); plot(res); title(‘残差序列‘); subplot(2,2,2); histogram(res); title(‘残差直方图‘); subplot(2,2,3); autocorr(res); title(‘残差ACF‘); subplot(2,2,4); parcorr(res); title(‘残差PACF‘); % 进行Ljung-Box Q检验，检验残差是否自相关 [h, pValue] = lbqtest(res, ‘Lags‘, [5, 10, 15]); % h=0表示不拒绝原假设（残差是白噪声），模型通过检验。

步骤4：进行预测

numYears = 5; % 预测未来5年 [YF, YMSE] = forecast(EstMdl, numYears, ‘Y0‘, area); % YF: 预测值 % YMSE: 预测均方误差 lower = YF - 1.96*sqrt(YMSE); % 95%置信区间下限 upper = YF + 1.96*sqrt(YMSE); % 95%置信区间上限 % 绘制预测图 figure; plot(year, area, ‘b-o‘, ‘LineWidth‘, 1.5); hold on; futureYears = (year(end)+1):(year(end)+numYears); plot(futureYears, YF, ‘r–s‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 8); fill([futureYears, fliplr(futureYears)], [lower‘, fliplr(upper‘)], ‘r‘, ‘FaceAlpha‘, 0.2, ‘EdgeColor‘, ‘none‘); legend(‘历史数据‘, ‘预测值‘, ‘95%置信区间‘, ‘Location‘, ‘best‘); xlabel(‘年份‘); ylabel(‘臭氧空洞面积‘); title(‘ARIMA模型预测结果‘); grid on;

实操心得与避坑指南：

1. 差分过度：差分阶数d不是越大越好。过度差分会使序列方差变大，并可能引入不必要的结构。通常d=0,1,2，通过ADF检验和观察差分后序列图综合判断。
2. 常数项问题：arima模型默认包含常数项。如果原始序列差分后均值接近0，可以在estimate时指定‘Constant‘, 0来拟合无常数项模型，有时能简化模型。
3. 收敛问题：estimate可能不收敛，尤其是当初始参数设置不合理或数据量太少时。可以尝试使用‘Display‘, ‘full‘查看迭代过程，或使用arima的‘ARLags‘, ‘MALags‘等参数手动指定初始值。
4. 预测起点：forecast函数的‘Y0‘参数是用于条件预测的初始数据。通常传入整个历史序列area即可。预测的置信区间会随着预测步长增加而迅速变宽，这是时间序列预测的不确定性本质。

3.2 灰色预测GM(1,1)实现

灰色预测的MATLAB实现相对紧凑，我们可以自己编写一个函数。

步骤1：编写GM(1,1)核心函数

function [predict, a, b, C, P] = gm11(x0, num_predict) % GM(1,1)灰色预测模型 % 输入： % x0: 原始数据行向量 (1 x n) % num_predict: 需要预测的步数 % 输出： % predict: 预测值（包括历史拟合值和未来预测值）(1 x (n+num_predict)) % a: 发展系数 % b: 灰色作用量 % C: 后验差比 % P: 小误差概率 n = length(x0); % 1. 累加生成(1-AGO) x1 = cumsum(x0); % 2. 构造数据矩阵B和数据向量Y B = [-0.5*(x1(1:end-1) + x1(2:end))‘, ones(n-1, 1)]; Y = x0(2:end)‘; % 3. 最小二乘估计参数 a, b ab = (B‘ * B) \ (B‘ * Y); a = ab(1); b = ab(2); % 4. 建立时间响应函数（离散解） % 预测累加序列 x1(k+1) = (x0(1)-b/a)*exp(-a*k) + b/a k = 0:(n + num_predict - 1); x1_hat = (x0(1) - b/a) * exp(-a * k) + b/a; % 5. 累减还原，得到预测值 predict = [x1_hat(1), diff(x1_hat)]; % diff计算后项减前项 % 6. 模型检验（仅对历史数据部分） x0_hat = predict(1:n); % 历史拟合值 epsilon = x0 - x0_hat; % 残差 delta = abs(epsilon ./ x0); % 相对误差 % 计算后验差比C和小误差概率P S1 = std(x0); % 原始序列标准差 S2 = std(epsilon); % 残差标准差 C = S2 / S1; avg_epsilon = mean(epsilon); S0 = 0.6745 * S1; % 计算S0 e = abs(epsilon - avg_epsilon); P = sum(e < S0) / n; % 输出模型精度等级（参考） if (P > 0.95 && C < 0.35) level = ‘好‘; elseif (P > 0.80 && C < 0.50) level = ‘合格‘; elseif (P > 0.70 && C < 0.65) level = ‘勉强合格‘; else level = ‘不合格‘; end fprintf(‘发展系数 a = %.4f，灰色作用量 b = %.4f\n‘, a, b); fprintf(‘后验差比 C = %.4f，小误差概率 P = %.4f，模型精度等级: %s\n‘, C, P, level); end

步骤2：应用与预测

% 使用数据 x0 = area‘; % 确保是行向量 num_predict = 5; [predict_all, a, b, C, P] = gm11(x0, num_predict); % 分离历史拟合和未来预测 fit_values = predict_all(1:length(x0)); future_values = predict_all(length(x0)+1:end); % 绘图对比 figure; plot(year, x0, ‘bo-‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 8, ‘DisplayName‘, ‘原始数据‘); hold on; plot(year, fit_values, ‘r^--‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 6, ‘DisplayName‘, ‘拟合值‘); futureYears = (year(end)+1):(year(end)+num_predict); plot(futureYears, future_values, ‘gs--‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 10, ‘DisplayName‘, ‘预测值‘); xlabel(‘年份‘); ylabel(‘臭氧空洞面积‘); title(‘灰色预测GM(1,1)模型结果‘); legend(‘Location‘, ‘best‘); grid on;

实操心得与避坑指南：

1. 数据量要求：GM(1,1)理论上最少需要4个数据点，但实际应用中，数据量太少（如少于7个）会导致模型对随机波动过于敏感，外推预测风险极高。建议至少使用7-10个以上数据点。
2. 发展系数a的符号：a的符号决定了趋势。a < 0时，模型呈指数增长趋势；a > 0时，呈指数衰减趋势。如果臭氧消耗数据在下降期，理论上a应为正。如果符号与物理意义相反，需要检查数据或模型适用性。
3. 精度检验：后验差比C越小越好（<0.35优秀），小误差概率P越大越好（>0.95优秀）。如果检验不合格，说明原始数据可能不适合直接用GM(1,1)。可以考虑对数据做平移变换（所有数据加上一个常数，使序列更平滑）或使用其他灰色模型（如DGM、灰色Verhulst模型）。
4. 长期预测风险：灰色预测基于指数规律，长期外推时，增长或衰减会非常迅速，可能与物理实际不符。因此它更适用于短期趋势预测。在赛题中，预测未来3-5年相对可靠，预测10年以上需非常谨慎，并必须在论文中说明此局限性。

3.3 多元线性回归模型实现

假设我们除了臭氧面积area，还有两个潜在影响因素：CFCs排放量cfc和太阳活动指数solar。目标是建立area = b0 + b1*cfc + b2*solar的模型。

步骤1：数据准备与可视化

% 假设数据为列向量: year, area, cfc, solar % 绘制散点图矩阵，观察两两关系 figure; plotmatrix([area, cfc, solar]); % 快速查看关系 title(‘变量间散点图矩阵‘); % 计算相关系数矩阵 corr_matrix = corrcoef([area, cfc, solar]); disp(‘相关系数矩阵:‘); disp(corr_matrix);

步骤2：建立回归模型与检验

% 使用fitlm函数建立线性模型 tbl = table(cfc, solar, area, ‘VariableNames‘, {‘CFC‘, ‘Solar‘, ‘Area‘}); mdl = fitlm(tbl, ‘Area ~ CFC + Solar‘); % 公式表示Area由CFC和Solar线性解释 % 显示详细的回归结果摘要 disp(mdl); % 关键输出解读： % R-squared: 决定系数，越接近1模型拟合越好。 % Adjusted R-squared: 调整后的R方，考虑了自变量个数，更可靠。 % F-statistic vs. constant model: F检验的p值，若小于0.05，说明模型整体显著。 % pValue for each coefficient: 每个系数（包括截距）的t检验p值，小于0.05说明该变量显著。

步骤3：模型诊断（非常重要！）一个合格的回归分析必须进行模型诊断，检查假设是否成立。

figure; subplot(2,2,1); plotResiduals(mdl, ‘fitted‘); % 残差 vs 拟合值图 title(‘残差-拟合值图‘); % 应随机分布，无规律，检验同方差性 xlabel(‘拟合值‘); ylabel(‘残差‘); subplot(2,2,2); plotResiduals(mdl, ‘lagged‘); % 残差 vs 滞后残差图（检验自相关） title(‘残差自相关图‘); xlabel(‘滞后残差‘); ylabel(‘残差‘); subplot(2,2,3); plotResiduals(mdl, ‘probability‘); % 正态概率图 title(‘正态概率图‘); % 点应近似在一条直线上，检验残差正态性 subplot(2,2,4); plotDiagnostics(mdl, ‘cookd‘); % Cook‘s距离，检测强影响点 title(‘Cook‘s距离‘); xlabel(‘观测序号‘); ylabel(‘Cook‘s距离‘); % 多重共线性诊断：计算方差膨胀因子(VIF) % VIF = 1 / (1 - R_i^2)， R_i^2是第i个自变量对其他自变量回归的R方。 % VIF > 5 或 10 表示存在严重共线性。 X = [cfc, solar]; [~, ~, ~, ~, stats] = regress(cfc, [ones(size(solar)), solar]); % 以cfc为因变量，solar为自变量 R2_cfc = stats(1); VIF_cfc = 1 / (1 - R2_cfc); % 同理计算VIF_solar [~, ~, ~, ~, stats] = regress(solar, [ones(size(cfc)), cfc]); R2_solar = stats(1); VIF_solar = 1 / (1 - R2_solar); fprintf(‘CFC的VIF: %.2f, Solar的VIF: %.2f\n‘, VIF_cfc, VIF_solar);

步骤4：进行预测要进行预测，我们需要未来年份的CFC和Solar数据。这在现实中很难获取，但在赛题中，有时会给出假设情景（如CFC按一定比例递减），或需要我们根据历史趋势外推（这又引入了新的预测误差）。

% 假设我们通过其他方法（如简单趋势外推）得到了未来5年的cfc_future和solar_future % cfc_future, solar_future 是长度为5的向量 future_data = table(cfc_future‘, solar_future‘, ‘VariableNames‘, {‘CFC‘, ‘Solar‘}); % 使用predict函数进行预测，并计算置信区间 [area_pred, area_ci] = predict(mdl, future_data); % area_pred: 点预测值 % area_ci: 置信区间（默认95%） % 绘图 figure; plot(year, area, ‘b-o‘, ‘LineWidth‘, 1.5); hold on; plot(futureYears, area_pred, ‘r–s‘, ‘LineWidth‘, 1.5, ‘MarkerSize‘, 8); plot(futureYears, area_ci(:,1), ‘g--‘, ‘LineWidth‘, 1); % 置信区间下限 plot(futureYears, area_ci(:,2), ‘g--‘, ‘LineWidth‘, 1); % 置信区间上限 legend(‘历史数据‘, ‘回归预测值‘, ‘95%置信区间‘, ‘Location‘, ‘best‘); xlabel(‘年份‘); ylabel(‘臭氧空洞面积‘); title(‘多元线性回归预测结果‘); grid on;

实操心得与避坑指南：

1. 共线性陷阱：如果CFC和Solar高度相关（VIF很大），回归系数的估计会变得不稳定，解释性变差。解决方法包括：剔除其中一个变量、使用主成分回归(PCR)或岭回归(Ridge Regression)等有偏估计方法。在MATLAB中，可以使用ridge函数进行岭回归。
2. 自相关问题：时间序列数据做回归，残差常常存在自相关（见‘lagged‘图），这违背了回归模型的独立同分布假设，会导致标准误低估，t检验失效。解决方法：使用时间序列回归模型，如加入滞后项，或使用arima模型拟合带有外生变量的ARIMAX模型（MATLAB中arima模型支持外生回归因子‘X‘参数）。
3. 预测的局限性：回归预测的准确性严重依赖于自变量的未来值是否准确。如果自变量的预测本身误差很大，那么因变量的预测误差会被放大。在论文中必须明确指出这一不确定性来源。
4. 交互项与非线性：如果怀疑CFC和Solar对臭氧的影响不是独立的（例如，高CFC时，Solar的影响更大），可以加入交互项CFC*Solar。如果关系是非线性的，可以考虑多项式回归或变量变换（如取对数）。

4. 模型比较、组合与结果分析

单一模型总有局限。在数学建模竞赛中，展示对不同模型的深刻理解和综合运用能力，是获得高分的关键。

4.1 模型评价指标

为了公平比较不同模型的预测性能，我们需要在历史数据上划分出一部分作为“测试集”（例如，用前20年数据建模，预测后5年，并与真实后5年数据对比）。常用指标有：

• 均方根误差 (RMSE):sqrt(mean((y_true - y_pred).^2))，衡量预测值与真实值的平均偏差，对较大误差更敏感。
• 平均绝对误差 (MAE):mean(abs(y_true - y_pred))，衡量平均绝对偏差，更稳健。
• 平均绝对百分比误差 (MAPE):mean(abs((y_true - y_pred) ./ y_true)) * 100，反映相对误差，便于不同量级序列比较。
• 决定系数 (R²): 在测试集上计算，越接近1越好。

在MATLAB中实现：

% 假设我们将最后5年数据作为测试集 train_idx = 1:(length(area)-5); test_idx = (length(area)-4):length(area); area_train = area(train_idx); area_test = area(test_idx); year_train = year(train_idx); year_test = year(test_idx); % 分别用ARIMA、灰色、回归模型在训练集上建模，并预测测试集年份 % ... (此处省略各模型在训练集上的拟合代码，与前面类似，只是数据换成area_train) % 假设得到了三个模型的测试集预测值：pred_arima, pred_gm11, pred_reg % 计算评价指标 metrics = table(); model_names = {‘ARIMA‘, ‘GM(1,1)‘, ‘Regression‘}; preds = {pred_arima, pred_gm11, pred_reg}; for i = 1:3 pred = preds{i}; rmse = sqrt(mean((area_test - pred‘).^2)); mae = mean(abs(area_test - pred‘)); mape = mean(abs((area_test - pred‘) ./ area_test)) * 100; % 计算测试集R² SS_res = sum((area_test - pred‘).^2); SS_tot = sum((area_test - mean(area_test)).^2); r2_test = 1 - SS_res / SS_tot; metrics = [metrics; table(model_names{i}, rmse, mae, mape, r2_test, ... ‘VariableNames‘, {‘Model‘, ‘RMSE‘, ‘MAE‘, ‘MAPE(%)‘, ‘R2_test‘})]; end disp(‘模型在测试集上的性能比较:‘); disp(metrics);

4.2 组合预测策略

如果多个单一模型各有优劣，可以采用组合预测来集成它们的智慧，降低风险。

• 等权平均：最简单，pred_comb = (pred_arima + pred_gm11 + pred_reg) / 3。
• 加权平均：根据各模型在测试集上的表现（如RMSE的倒数）分配权重。

  rmse_vec = [metrics.RMSE(1), metrics.RMSE(2), metrics.RMSE(3)]; weights = (1 ./ rmse_vec) / sum(1 ./ rmse_vec); % RMSE越小，权重越大 pred_comb_weighted = weights(1)*pred_arima + weights(2)*pred_gm11 + weights(3)*pred_reg;

• 最优线性组合：使用回归方法，以各单一模型的预测值为自变量，真实值为因变量，拟合一个线性组合模型。这相当于让数据自己决定如何组合。

4.3 结果分析与论文表述要点

在论文写作中，不能只罗列代码和图表，必须有深入的分析。

1. 模型优缺点分析：

   • ARIMA：优点在于严格的数据驱动，能捕捉序列的自相关结构，提供置信区间。缺点是对数据平稳性要求高，且纯时间序列模型无法利用外部影响因素信息。
   • 灰色预测：优点是小样本需求，对趋势数据拟合好，计算简单。缺点是假设数据服从指数规律，对波动数据适应性差，长期预测可靠性低。
   • 多元回归：优点是有明确的因果解释，如果自变量预测准，则预测说服力强。缺点是严重依赖自变量的准确未来值，且容易受共线性、自相关等问题干扰。

2. 预测结果解读：

   • 将最终预测结果（可能是组合模型的结果）以图表形式清晰展示。
   • 结合《蒙特利尔议定书》等现实背景，解释预测趋势的合理性。例如，预测显示臭氧空洞面积在未来五年呈缓慢下降趋势，这与国际社会持续减排CFCs的努力相符。
   • 必须讨论预测的不确定性：明确指出置信区间的含义（例如，我们有95%的把握认为真实值落在这个区间内），并说明随着预测时间变长，置信区间会变宽，预测不确定性增加。

3. 灵敏度分析：

   • 这是一个加分项。可以探讨如果某个关键参数（如CFC减排速率）发生变化，预测结果会如何改变。这展示了模型的鲁棒性和你对问题理解的深度。

   % 例如，在回归模型中，假设CFC未来以-2%/年或-5%/年的速度变化，重新预测臭氧面积。 cfc_scenario1 = cfc_future * 0.98.^(1:5)‘; % 每年减少2% cfc_scenario2 = cfc_future * 0.95.^(1:5)‘; % 每年减少5% % 然后分别代入模型预测，比较结果差异。

5. 常见问题、排查技巧与进阶思考

在实际操作和竞赛中，你会遇到各种各样的问题。这里记录一些典型问题和我的解决思路。

5.1 MATLAB环境与代码问题

• 问题：arima/estimate函数未定义？

  • 排查：确认是否安装了Econometrics Toolbox。在命令行输入ver查看已安装的工具箱列表。
  • 解决：如果没有，需要安装该工具箱。对于学生，可以使用学校提供的正版授权；对于个人，可以考虑MathWorks的Home版或试用版。

• 问题：运行灰色预测函数，预测值出现NaN或Inf？

  • 排查：检查发展系数a是否非常接近0。在时间响应函数x1_hat = (x0(1)-b/a)*exp(-a*k) + b/a中，如果a接近0，b/a会趋于无穷大。
  • 解决：这通常意味着原始数据几乎没有增长或衰减趋势，不满足灰色模型的基本假设。应考虑换用其他模型，或对数据进行预处理（如分析其是否平稳，尝试ARMA）。

• 问题：回归诊断图中残差呈现明显的“漏斗形”或“弯曲形”？

  • 排查：这违反了同方差性假设，意味着误差方差随预测值增大而改变（漏斗形），或存在非线性关系（弯曲形）。
  • 解决：
    1. 对于异方差（漏斗形），可以对因变量进行变换，如取对数log(area)，或使用加权最小二乘法。
    2. 对于非线性（弯曲形），可以在模型中添加自变量的平方项（多项式回归），或使用更复杂的非线性回归方法。

5.2 模型与方法层面的问题

• 问题：数据既有趋势又有季节性（如月度数据），怎么办？

  • 解决：这是更经典的时间序列问题。可以使用季节性ARIMA模型 (SARIMA)。MATLAB中arima模型可以通过‘Seasonality‘参数指定季节性周期。例如，对于月度数据，周期s=12。模型记为SARIMA(p,d,q)×(P,D,Q)_s。这需要更复杂的模型识别过程。

• 问题：影响因素很多，如何选择进入回归模型的自变量？

  • 解决：不要一股脑全放进去。可以使用逐步回归（stepwiselm函数）来自动筛选变量。或者基于领域知识先筛选，然后通过方差膨胀因子(VIF)和显著性检验(p值)来剔除不显著或共线性严重的变量。原则是：在保证解释力的前提下，模型越简洁越好（奥卡姆剃刀原理）。

• 问题：感觉ARIMA和回归模型都挺好，到底该选哪个作为最终模型？

  • 心法：没有绝对的最好，只有最合适。如果外部驱动因素明确且可预测，回归模型的理论解释性更强。如果系统复杂，影响因素难以量化或获取，那么从数据自身出发的ARIMA可能更稳健。在竞赛中，展示你尝试了多种模型，并基于数据特征和检验结果进行了有理有据的选择和组合，这个过程比单纯给出一个“正确”答案更重要。

5.3 竞赛论文写作要点

1. 问题重述不是照抄：用自己的话简洁概括问题背景、目标和已知条件。
2. 模型假设要合理且明确：例如，“假设未来五年没有新的强效臭氧消耗物质被大量排放”、“假设太阳活动周期遵循历史规律”等。这界定了你模型的适用范围。
3. 符号说明要清晰：在模型建立前，用表格列出所有用到的主要符号及其含义。
4. 流程图是利器：用清晰的框图展示你的整体建模步骤（数据预处理→模型选择与建立→模型检验→预测→分析），能让评委快速抓住你的思路。
5. 图表规范美观：MATLAB出图后，注意调整线条粗细、标记大小、图例位置、坐标轴标签字体，保证在论文中清晰可读。一张专业的图表能极大提升印象分。
6. 模型检验部分不可或缺：无论是ARIMA的残差白噪声检验、灰色预测的后验差检验，还是回归的多种诊断图，都必须展示并简要说明结果，证明你的模型是可靠的。
7. 优缺点与推广实事求是：在结论部分，客观总结所用模型的优点和局限性，并提出可能的改进方向（例如，“未来可收集更多影响因素，尝试机器学习模型如随机森林进行预测”）。

这道“臭氧消耗预测”题，就像一把钥匙，打开了时间序列预测和回归分析的大门。它教会我们的，不是某个特定模型的死记硬背，而是一套面对数据、提出问题、选择工具、验证结果、解释现实的完整科学思维流程。真正有价值的不只是那几行MATLAB代码，而是在反复调试模型、对比结果、解读图表过程中，培养出的那种对数据的感觉和对模型局限性的清醒认知。当你下次遇到销售预测、能源需求预测、传染病传播预测等问题时，你会发现，思路都是相通的。