简介：面向需要处理季节性时间序列的 MATLAB 用户，资源以 SARIMA（季节性差分自回归滑动平均模型）为核心，覆盖从数据预处理、平稳性检验、ACF/PACF 与 AIC/BIC 参数寻优，到模型构建、拟合评估、预测及优化调整的完整流程，适合高校学生、数据分析师和科研人员快速上手。包体内共 13 个文件，以 .m 脚本为主，配以 .jpg 可视化图、.xlsx/.mat 数据文件及 .docx 说明文档，压缩包整体约 112KB，结构清晰便于对照学习。已有 1479 人学习下载。按步骤拆分的代码脚本和实验数据可直接运行，有助于理解季节性差分、参数定阶与预测效果评估等关键环节，也能为后续二次开发与业务预测节省建模时间。
1. SARIMA 时间序列预测分析：这份 MATLAB 资源包拆开看
拿到一组带明显季节波动的业务数据——航空客运量、月度零售额、电网负荷，第一反应就是用 SARIMA 试试。SARIMA（季节性差分自回归滑动平均模型）这类模型在月度、季度数据上可以说是性价比很高的选择，它比 ETS 更灵活，比 LSTM 可解释性强得多。你这份sarima.zip里的内容很标准：Demo_SARIMA.m是总入口，creatSARIMA.m负责建模，SARMA_Order_Select.m做定阶选参，Fun_SARIMA_Forecast.m管预测，还带了两套 Excel 数据和一份经典的Data_Airline.mat航空旅客数据，外加四张可视化图像和一份说明文档。也就是说，从数据检查到参数选择、模型拟合、预测评估，这条链路在这个包里是完整闭环的。适合谁？打算在 MATLAB 里落地时间序列预测、又不想完全靠 GUI 点选的人，以及刚接触 SARIMA、想找一套能直接跑通再改参数的参考代码的人。接下来我会按实际操作的顺序，把这个包里的逻辑和坑逐个拆开。
2. SARIMA 模型原理与参数定阶：先看懂季节周期，再谈调参
2.1 从 ARIMA 到 SARIMA：季节差分到底解决了什么
ARIMA(p,d,q) 处理的是非平稳序列，公式可以写成：
(1 - φ1B - ... - φpBp)(1-B)d Xt = (1 + θ1B + ... + θqBq) εt普通差分(1-B)d能消掉线性趋势，但对周期为 12 的月度数据，只做一阶差分后序列里依然残留明显的年度季节性——自相关图在滞后 12、24 的位置仍然显著。原因很简单：趋势是逐点变化的，而季节性是每隔 S 期重复的规律，需要在滞后 S 的位置再做一次差分。SARIMA 就是在此基础上加了季节差分项和季节 AR/MA 项：
(1-B)d (1-BS)D Xt = [φ(B) Φ(BS)]-1 [θ(B) Θ(BS)] εt其中(1-BS)D是季节差分，比如 D=1、S=12 时算的是Xt - Xt-12，意思是"今年 1 月比去年 1 月多了多少"。
在 MATLAB 里没有独立的sarima函数，正常的做法是用arima模型设置Seasonality属性来等价实现。示例代码如下：
% 构建 SARIMA(1,1,1)(1,1,1)_12 mdl = arima('Constant', 0, ... 'ARLags', 1, ... % 普通 AR 滞后阶数 p=1 'MALags', 1, ... % 普通 MA 滞后阶数 q=1 'SARLags', 12, ... % 季节 AR 滞后为 S=12 'SMALags', 12, ... % 季节 MA 滞后为 S=12 'D', 1, ... % 普通差分阶数 d=1 'Seasonality', 12); % 季节差分周期 S=12，等价 D=1逻辑说明：ARLags和MALags控制非季节部分的阶数，SARLags和SMALags必须填季节周期的整数倍（这里填 12），而不是填季节阶数 P 的值。比如你想做 SARIMA(1,1,1)(1,1,1)12，那SARLags就写 12 而不是 1。这一点和 R 语言里Arima(order=c(1,1,1), seasonal=c(1,1,1))的参数语义不太一样，刚转过来的人容易混。如果季节周期是季度数据 S=4，则SARLags写 4、SMALags写 4。Seasonality一旦设置成大于 0 的整数，模型会自动对数据做一阶季节差分，等价于公式里的 D=1。
另一个容易忽略的点：Constant要不要设成 0。如果做了差分（d≥1 或 D≥1），模型里通常不建议再加常数项，否则预测会带出人为的漂移趋势。creatSARIMA.m里如果看到'Constant', 0的写法，就是这个原因。
2.2 ACF/PACF 怎么读：滞后 12 处显著才是季节信号
定阶的核心逻辑分两层：先看低滞后（1～3）的 ACF/PACF 判断 p 和 q，再看滞后 S、2S 处的表现判断 P 和 Q。这属于标准流程，但很容易出现误读——特别是把季节信号当成普通自相关来处理，结果阶数越调越离谱。
普通部分的判断规则是：ACF 在 q 阶后截尾、PACF 拖尾，适合 MA(q)；ACF 拖尾、PACF 在 p 阶后截尾，适合 AR(p)；两个都拖尾，考虑 ARMA 或提高差分阶数。季节部分的规则类似，只是把观察点换到 S 的整数倍滞后上。举个例子：月度数据 S=12，如果 ACF 在滞后 12、24 处显著不为 0，之后快速衰减，而 PACF 在滞后 12 处有明显的单峰，这说明季节部分偏 MA(1)，也就是 Q=1、P=0。
SARMA_Order_Select.m这个脚本解决的是"AIC/BIC 网格搜索"那部分——你不必只靠眼睛看图猜，直接对一组候选 (p,q,P,Q) 组合逐个拟合并输出信息准则值。一个常用的网格搜索模板：
% 网格定阶：在 p,q,P,Q 的候选范围内遍历 pList = 0:2; qList = 0:2; PList = 0:1; QList = 0:1; results = []; for p = pList for q = qList for P = PList for Q = QList try mdl = arima('Constant', 0, ... 'ARLags', 1:p, ... 'MALags', 1:q, ... 'SARLags', 12*(1:P), ... 'SMALags', 12*(1:Q), ... 'D', 1, 'Seasonality', 12); [fitMdl, ~, logL] = estimate(mdl, Y, 'display', 'off'); [aic, bic] = aicbic(logL, fitMdl.NumEstimatedParameters, numel(Y)); results = [results; p, q, P, Q, aic, bic]; catch % 某些参数组合会收敛失败，跳过即可 end end end end end % 按 AIC 升序排列，看前几名 results = sortrows(results, 5); disp(results(1:5, :));逻辑说明：SARLags,1:p展开后是[1,2]，对应普通 AR 阶数 2；SARLags,12*(1:P)展开后如果 P=1 就是 12，P=2 则为[12,24]。estimate的第三输出logL是对数似然值，aicbic根据参数个数和样本量计算 AIC/BIC。try/catch在这里是必要的——有些组合比如 p=0 且 P=0 且 D=1 时，模型可能无法收敛，直接报错中断整个循环，跳过去是最省事的做法。
参数说明：AIC 适合预测导向的选模，BIC 对参数个数惩罚更重，倾向于给出更精简的模型。两者排名如果差别很大，说明候选模型里有明显的过拟合嫌疑，优先看 BIC 的结果。另外网格搜索耗时随候选数指数增长，这份数据样本量只有一百多个点，跑完没问题；如果是高频数据上千个点，建议把候选范围先缩到 p,q∈{0,1,2}, P,Q∈{0,1}，先粗筛一轮再细化。
2.3 数据中的季节周期 S 怎么确认
很多人一上来就默认月度数据 S=12、季度数据 S=4，这在多数场景下没问题，但拿到新数据还是先验证一下比较稳妥。一个简单的做法是直接对序列做频谱分析或画周期图：
% 用傅里叶变换找主周期 Y = Data_Airline; % 或者读入你自己的序列 Fs = 1; % 采样频率，月度数据为 1 个月 1 个点 nfft = 2^nextpow2(length(Y)); f = Fs * (0:(nfft/2))/nfft; Yfft = abs(fft(Y - mean(Y), nfft)); plot(f(1:nfft/2+1), Yfft(1:nfft/2+1)); xlabel('频率（周期/月）'); ylabel('幅值');逻辑说明：先减去均值再做 FFT，是为了避免直流分量把低频峰值盖掉。横轴频率如果看到 0.083 附近有明显峰值，即 1/12，那就确认周期为 12 个月。Data_Airline.mat这份航空旅客数据是月度流量，周期就是 12，直接按 S=12 处理没问题；但如果是周内季节性（比如电网日负荷数据周期为 7 天），S 就要改成 7。这个坑在 4.3 节还会展开说。
3. 在 MATLAB 里跑通 SARIMA：从数据预处理到模型拟合
3.1 数据预处理：先看图，再做平稳性检验
拿到data.xlsx或Data_Airline.mat之后，第一步不是建模，而是把序列画出来看形态。1.jpg到4.jpg大概率就是包作者跑出来的图：原始序列折线图、直方图、ACF/PACF 图，以及预测对比图。以航空数据为例，原始序列有明显的上升趋势和逐年振幅增大的季节性，直方图会呈现双峰或长尾分布，这说明序列既不平稳、方差也不稳定。
% 加载航空数据 load('Data_Airline.mat'); % 变量名通常是 Data 或 Y Y = Data(:); % 按列读取，转为列向量 % 快速三连画：原始序列、ACF、PACF tiledlayout(3,1); nexttile; plot(Y); title('原始序列'); nexttile; autocorr(Y); title('ACF'); nexttile; parcorr(Y); title('PACF');逻辑说明：Data(:)这种写法是为了避免 Excel 或 mat 文件里数据是行向量时后续代码出现维度不匹配——estimate和forecast都要求列向量输入。autocorr和parcorr是 Econometrics Toolbox 的函数，分别输出自相关图和偏自相关图。如果 ACF 图衰减非常缓慢（像一条从左上角滑下来的斜坡），说明序列有强趋势，必须先差分。
预处理阶段还有一步经常被跳过：方差稳定化。航空数据这种"振幅随水平升高而变大"的形态，比较规范的做法是先取对数再建模，把乘法性季节性变成加法性。判断方法很简单——看 3.1 第一张折线图，如果波动幅度和序列水平成正比，就做logY = log(Y)再继续。
3.2 平稳性检验：ADF 检验和差分阶数确定
差分阶数 d 的确定不是靠肉眼，而是靠单位根检验。MATLAB 的adftest和kpss可以一起用：ADF 的原假设是序列存在单位根（非平稳），KPSS 的原假设是序列平稳。两个检验可以互相印证——如果 ADF 的 p 值大于 0.05 且 KPSS 的 p 值小于 0.05，说明序列非平稳，需要差分。
% ADF 检验原始序列 [h1, pval1] = adftest(Y, 'model', 'TS'); % 含趋势项模型 fprintf('ADF 检验: h=%d, p=%.4f\n', h1, pval1); % 如果非平稳，做一阶差分后再检验 if h1 == 0 dY = diff(Y); [h2, pval2] = adftest(dY, 'model', 'AR'); fprintf('一阶差分后 ADF 检验: h=%d, p=%.4f\n', h2, pval2); end逻辑说明：adftest的model参数可选'AR'（只有截距）、'ARD'（截距加趋势）、'TS'（趋势平稳）。对有明显上升趋势的航空数据，原始序列检验用'TS'更合理；差分后用'AR'即可。这里h=0表示不能拒绝单位根原假设，即序列非平稳；h=1表示拒绝原假设，认为序列平稳。
注意一个常见问题：一阶差分后 ADF 的 p 值可能已经显著，但这只能说明普通差分够用，不代表季节性消掉了。Data_Airline.mat的序列一阶差分后，ACF 在滞后 12 处依然有明显峰值，这时就需要再加季节差分（也就是Seasonality, 12）。这个判断不能省——很多人只做了 d=1 就建模，结果残差在滞后 12 处依然显著自相关，整个模型等于白做。
3.3 模型拟合：estimate 的用法和结果解读
差分阶数确认后，就可以把模型放进estimate里拟合了。creatSARIMA.m里做的事情基本就是构建模型对象、传入数据、输出参数估计。一个完整的拟合流程结合参数解读如下：
load('Data_Airline.mat'); Y = log(Data(:)); % 先取对数再建模 % 1. 构建 SARIMA(1,1,1)(1,1,1)_12 mdl = arima('Constant', 0, ... 'ARLags', 1, 'MALags', 1, ... 'SARLags', 12, 'SMALags', 12, ... 'D', 1, 'Seasonality', 12); % 2. 拟合 estMdl = estimate(mdl, Y); % 3. 残差诊断：Ljung-Box Q 检验 [res, ~] = infer(estMdl, Y); [h, pValue] = lbqtest(res, 'Lags', [12, 24]); fprintf('Ljung-Box 检验: h=%d, p=%.4f\n', h, pValue);逻辑说明：infer用于从拟合好的模型中提取残差序列，这一步通常在预测前做，目的是确认残差是否为白噪声。lbqtest的Lags参数指定要检验的滞后阶数——对于月度季节数据，至少要看 12 和 24 这两个滞后，因为如果季节信息没提取干净，残差会在滞后 12 处表现出显著自相关，p 值小于 0.05，模型就需要修正。
参数说明：estimate估计出的结果里，AR{1}是 φ1、MA{1}是 θ1、SAR{12}是季节 AR 系数、SMA{12}是季节 MA 系数。如果某个系数估计值的 t 统计量不显著（|t|<1.96），说明这个阶数可能设多了，可以回 2.2 节把对应的 p 或 q 降一阶再拟合。这里有一个实战经验：季节 MA 系数SMA{12}如果估计值接近 1（比如 0.98 以上），模型可能已经接近不可逆边界，预测方差会被放大，遇到这种情况优先尝试减少 Q 而不是继续调其他参数。
3.4 预测与评估：Fun_SARIMA_Forecast.m 这条链路
Fun_SARIMA_Forecast.m是个函数文件，负责把训练好的模型向前预测 N 步。MATLAB 的forecast函数在没有并行计算时需要指定计算方式，通常建议直接用默认的'exact'方法：
% 假设 estMdl 已经在 3.3 节拟合好 numPeriods = 12; % 预测未来 12 个月 [Ypred, YMSE] = forecast(estMdl, numPeriods, 'Y0', Y); % 如果前面取了对数，预测值要转回原始尺度 Ypred_orig = exp(Ypred); lower = exp(Ypred - 1.96 * sqrt(YMSE)); upper = exp(Ypred + 1.96 * sqrt(YMSE)); % 可视化：历史 + 预测 + 置信区间 t = 1:length(Y); tPred = length(Y) + (1:numPeriods); figure; plot(t, exp(Y), 'k-'); hold on; plot(tPred, Ypred_orig, 'r-', 'LineWidth', 1.5); plot(tPred, lower, 'r--', tPred, upper, 'r--'); legend({'历史数据', '预测值', '95% 置信区间'});逻辑说明：Y0必须是拟合模型时用的同一条序列——注意这里有个细节，如果建模用的是log(Y)，那么forecast也要传log(Y)进去。YMSE是预测方差，转换成置信区间时要开根号。取对数后得到的区间是乘法区间，回到原始尺度后表现为上下置信线不对称，这是正常现象，不是 bug。
参数说明：forecast的第三个参数'Y0'传的是预测起点的历史信息。如果你在拟合后提取了残差，forecast会自动使用模型内部的残差状态；但如果你的模型包含移动平均项（q>0 或 Q>0），预测起点前的残差需要用到infer的结果，这时最稳妥的写法是先把残差存下来再调用forecast。实际操作中很多人直接在Y0里只传最后几个月的值，这会导致预测起点不正确——正确做法是把全部历史数据都传进去。
评估预测精度时，常用指标是 RMSE、MAE 和 MAPE，公式不多说，直接上代码：
% 假设 testY 是真实值，Ypred_orig 是对应预测值 rmse = sqrt(mean((testY - Ypred_orig).^2)); mae = mean(abs(testY - Ypred_orig)); mape = mean(abs((testY - Ypred_orig) ./ testY)) * 100; fprintf('RMSE=%.2f, MAE=%.2f, MAPE=%.2f%%\n', rmse, mae, mape);如果 MAPE 在 5% 以内，说明模型对这批数据拟合得相当好；超过 15% 基本要回炉重做，优先检查季节性周期和差分阶数。
4. SARIMA 常见问题与避坑：参数选错、过拟合、数据格式翻车
4.1 数据顺序错乱：Excel 乱序让模型学了个寂寞
现象：data.xlsx读进来之后画图发现序列杂乱无章，预测值来回跳，完全不像有规律的时间序列。原因：Excel 里的行序不是严格按时间递增排列的——有时是跨 sheet 拼接，有时是复制粘贴时错位。解决：读入后立刻按时间列排序，同时确保Y是列向量。不排序的话，SARIMA 会把乱序数据当成时间顺序来学，所有滞后结构全部失效。
4.2 差分阶数过大：拟合优度很漂亮，预测结果却发散
现象：模型的 AIC 很低，残差看起来也干净，但预测几步之后曲线迅速发散成喇叭形。原因：过度差分——数据本身只需要 d=1，你为了"让 ACF 看起来更干净"多差分了两次，把序列里的长期信息也差掉了。过度差分会导致预测方差被放大，表现就是置信区间越来越宽、预测均值偏离真实轨迹。解决：差分阶数只看 ADF 检验的结果——p 值小于 0.05 就停，不要靠肉眼判断"再差一次效果更好"。另外可以对比 d=1 和 d=2 的 RMSE，如果 d=2 的样本内拟合好了但样本外预测明显更差，妥妥的过拟合。
4.3 季节周期 S 设错：预测曲线怎么都对不上
现象：预测结果的波峰波谷位置和实际数据完全错位，比如实际 7 月是旺季，模型预测 9 月才起量。原因：S 设成了 12 但实际周期是 12 的另一倍或半倍——比如数据是"月度"但统计口径按 13 期记账；或者实际是周内周期（S=7），你按月份频率建了 S=12。解决：回到 2.3 节看周期图的主峰位置，确认周期后再设Seasonality和SARLags。另外注意，如果数据是季度频率但用了 S=12，那模型会把相邻年份同季度当成同一个月来建模，结果必然偏。
4.4 forecast 的 Y0 参数传错导致预测全是 NaN
现象：forecast执行成功后返回值全是 NaN 或者长度不对。原因：Y0传入的数据长度太短，或者传成了差分后的序列（diff(Y)），而模型内部期望的是原始序列的滞后值。解决：Y0一定传拟合时用的完整原始序列，不是差分序列。如果模型包含季节差分Seasonality，forecast会自动做内部差分，你不需要手动处理。再一个容易翻车的是：Y0和拟合数据长度不一致时 MATLAB 不会报错，但预测结果会从错误起点开始——务必检查forecast输出的前几个值和历史数据最后几期的衔接是否自然。
4.5 残差诊断被跳过：Ljung-Box 通不过你还在往下走
现象：拟合完成后直接去预测，结果看似合理，但模型在下一批数据上失效严重。原因：残差里还残留自相关结构，说明模型没把信息提取干净——最常见的残留位置就是滞后 12 处的季节性。解决：养成每次拟合后都跑一遍infer加lbqtest的习惯，滞后参数至少设[12, 24]。如果 p 值小于 0.05，说明残差不是白噪声，回到 2.2 节的网格搜索微调 P、Q 阶数，或检查数据是否漏了周期性。这一步是"后悔药"，越早做返工成本越低。
5. 把 SARIMA 用到生产环境：滚动预测、误差控制与模型边界
5.1 用滚动预测代替一次性切分
一次性把数据切成训练集和测试集，在时间序列预测里容易高估模型表现——因为测试集紧邻训练集，模型等于"猜"得很近的未来。更接近真实业务的验证方式是滚动预测（rolling forecast）：每次只用截止到当前月的历史数据训练，预测下一个月，然后真实值到达后再把窗口推进一步。代码如下：
load('Data_Airline.mat'); Y = log(Data(:)); numPeriods = 12; % 滚动预测 12 个月 trainSize = length(Y) - numPeriods; predictions = zeros(numPeriods, 1); for t = 1:numPeriods trainY = Y(1:trainSize + t - 1); % 窗口逐步扩大 mdl = arima('Constant', 0, ... 'ARLags', 1, 'MALags', 1, ... 'SARLags', 12, 'SMALags', 12, ... 'D', 1, 'Seasonality', 12); estMdl = estimate(mdl, trainY, 'display', 'off'); [pred, ~] = forecast(estMdl, 1, 'Y0', trainY); predictions(t) = pred(1); % 只取第一步预测 end % 对比真实值 testY = Y(trainSize + 1:end); plot(testY, 'k-', 'LineWidth', 1.2); hold on; plot(exp(predictions), 'r--', 'LineWidth', 1.2); legend('真实值', '滚动预测', 'Location', 'best');逻辑说明：每次循环都重新拟合一次模型，模拟的是生产环境的真实节奏——每月月底用最新数据更新模型，然后预测下个月。'display', 'off'是为了让循环不刷屏。注意forecast的第二个参数这里填的是 1，即每次只预测一步，预测值存入数组，下个循环再推进。
这套流程跑完，你得到的 RMSE/MAPE 才是模型真实水平的估计。一次性切分测试集的误差通常会比滚动预测好看 20%～30%，但那是虚的。
5.2 多步预测的误差累积控制
如果业务需要一次预测未来 12 个月，forecast(estMdl, 12, 'Y0', Y)直接给出 12 步预测，但误差会随步长增大而累积——SARIMA 的预测方差是单调递增的，前几个月还稳，后面几个月置信区间会明显变宽。想控制误差累积，常规做法是"滚动再估计"：每多获得一个真实值就重新拟合一次，或者用递归预测（每一步把上一步的预测值当真实值喂回去）。后者在 MATLAB 里没有内置函数，得手动循环：
Ypred = zeros(12, 1); tempY = Y; % 历史数据 for j = 1:12 [pred, ~] = forecast(estMdl, 1, 'Y0', tempY); Ypred(j) = pred(1); tempY = [tempY; pred(1)]; % 把预测值拼到序列尾部 end逻辑说明：递归预测的优势是代码简单、能跑出长期走势，代价是误差逐步累积——预测步数越往后越依赖之前的预测值。如果模型的 MA 系数较大，累积误差会膨胀得很快，这时我更建议走 5.1 的滚动再估计，哪怕重拟合的耗时多一些，精度提升是实打实的。
5.3 SARIMA 的边界：什么时候该换模型
SARIMA 不是万能药。如果你的数据存在节假日效应（比如春节、双十一）、多日历效应（比如当月工作日天数影响销量），或者突变的结构性断点，SARIMA 会把这些当成随机噪声处理。一个常见的判断依据：模型残差诊断做干净了，但滚动预测的表现还是跟不上业务精度要求——这时要反思的不是参数，而是模型结构。我一般会转去试 Prophet（对节假日和突变更友好）或者带外部回归项的 ARIMAX。还有一个小技巧：SARIMA 的预测区间假定误差正态同方差，如果你的数据预测误差本身就有季节性波动（比如旺季误差远大于淡季），可考虑对残差序列再做一次 GARCH 建模，把波动率的季节性也拟合进去。
这套流程走下来的心得是：SARMA_Order_Select.m跑完、出 AIC/BIC 表，只是万里长征第一步；真正决定模型生死的是残差白噪声检验和滚动预测。从那以后，我每次拿到带季节性的业务序列，都强制走一遍"画图 → 周期确认 → ADF 检验 → 网格定阶 → 残差诊断 → 滚动验证"，这套流程至少帮我把返工率降低了一半。希望帮到你。
本文还有配套的精品资源，点击获取
