简介本资源面向计算机、电子信息工程、数学等专业的大学生及算法入门者提供一套完整的ARIMA-SSA-LSTM时间序列预测实现方案可用于课程设计、期末大作业或毕业设计。代码基于Python与TensorFlow编写采用参数化编程思路参数可灵活调整且配有保姆级注释几乎一行一注释便于零基础读者理解建模流程。压缩包共3个文件包含2个CSV数据文件与1个Python源码文件整体约51KB数据与代码配套齐全可直接运行验证。资源将ARIMA、麻雀搜索算法SSA与LSTM相结合展示了从数据读取、参数优化到预测输出的完整链路读者可据此掌握组合模型的搭建方法与调参技巧。目前已有510人学习适合希望快速上手时间序列预测实战的读者参考。1. 从一条震荡的销量曲线说起ARIMA-SSA-LSTM 到底在解决什么去年帮一个做小家电的朋友看数据他手里有 36 个月的出货记录想预测未来 6 个月。我先用 ARIMA 跑了一版趋势抓得还行但每逢大促月份的尖峰全被抹平换成 LSTM 单模型尖峰能咬住一部分可样本量太小验证集上的曲线抖得像心电图。这不是模型不行是单一模型对「线性趋势 非线性突变」这种混合结构天然偏科。ARIMA-SSA-LSTM 这套组合就是冲这个痛点来的ARIMA 负责把线性成分剥出来SSA奇异谱分析把残差里的噪声和周期分量拆开、重构出更干净的序列LSTM 再去学剩下的非线性映射。它适合手里有几百到几千条单变量序列、又不想上重型深度模型的从业者比如销量、负荷、流量、GDP 这类按时间排队的数值预测。下面我按自己落地的顺序把数据怎么切、三个模块怎么串、参数怎么调、哪里最容易翻车讲清楚。2. 数据准备与平稳性处理ARIMA 的地基怎么打2.1 时间序列的四个前置检查拿到一列带时间戳的数值别急着往模型里灌。我一般先做四件事确认时间间隔是否等距、检查缺失值和异常点、看量纲是否需要归一化、画一遍原始曲线和差分曲线。等距是 ARIMA 的硬要求如果原始数据是「有单才记一笔」的事件流得先重采样成固定频率比如按天求和或按周取均值。缺失值不要直接dropna销量类数据缺失往往意味着当天真没卖出去填 0 比填均值更合理传感器类数据则适合线性插值。异常点用滚动中位数加 3 倍 MAD 判定超过阈值的先标记别急着删后面 SSA 重构时还能用上。平稳性决定了 ARIMA 的差分阶数 d。ADF 检验是常用手段p 值小于 0.05 认为平稳。但我要提醒一句ADF 对短序列和强周期序列经常误判最好配合 ACF/PACF 图一起看。如果 ACF 衰减极慢、PACF 在 1 阶后截尾基本就是需要差分的信号。差分不是越多越好过度差分会把序列里的有效信息也差掉一般 d 取 0、1、2 三档试优先 1 阶。import pandas as pd import numpy as np from statsmodels.tsa.stattools import adfuller # 读取数据假设两列ds(时间), y(数值) df pd.read_csv(series.csv, parse_dates[ds]).sort_values(ds) df df.set_index(ds).asfreq(D) # 统一成日频缺失处补 NaT df[y] df[y].interpolate(linear) # 线性插值补缺失 # 异常点标记滚动中位数 3倍MAD med df[y].rolling(7, centerTrue).median() mad (df[y] - med).abs().rolling(7, centerTrue).median() df[outlier] (df[y] - med).abs() 3 * 1.4826 * mad # ADF 检验 def adf_report(s, name): r adfuller(s.dropna(), autolagAIC) print(f{name}: stat{r[0]:.3f}, p{r[1]:.4f}, lags{r[2]}) adf_report(df[y], 原始序列) adf_report(df[y].diff(), 一阶差分)这段代码里asfreq(D)把不规则时间轴强制对齐到日频interpolate(linear)只补数值不补时间。MAD 判定用 1.4826 这个系数是为了让 MAD 在正态分布下与标准差可比这是稳健统计里的常规做法。ADF 的autolagAIC让函数自动选滞后阶数比手动指定稳。跑完看两个 p 值原始不平稳、一阶差分平稳d 就定 1。2.2 ARIMA 定阶别迷信 auto_arima定阶我踩过的坑最多。pmdarima的auto_arima很方便但它默认用 AIC 做准则短序列上容易选出过大的 p 和 q导致过拟合。我的习惯是先用auto_arima给一个候选范围再手动在 (p,d,q) 的邻域里用 BIC 和滚动验证比一遍。BIC 对参数个数惩罚更重小样本下更保守。季节性明显的序列比如月度数据还要加 (P,D,Q,s) 这组季节项s 取 12。from statsmodels.tsa.arima.model import ARIMA from sklearn.metrics import mean_absolute_error y df[y].values train, test y[:-12], y[-12:] best None for p in range(0, 4): for q in range(0, 4): try: m ARIMA(train, order(p, 1, q)).fit() pred m.forecast(len(test)) mae mean_absolute_error(test, pred) if best is None or mae best[0]: best (mae, p, q, m) except Exception: continue print(最优阶数:, best[1], best[2], MAE:, round(best[0], 3))这里用滚动外推的方式评估而不是只看训练集拟合优度。order(p,1,q)里的 1 就是前面定下的 d。forecast(len(test))做的是样本外预测比fittedvalues更能反映真实泛化。如果所有组合的 MAE 都很大说明线性部分解释力有限这恰恰是后面要上 SSA 和 LSTM 的信号不是失败。2.3 残差提取把线性部分「减」出去ARIMA 拟合完真正要交给下游的不是预测值而是残差。残差 真实值 − ARIMA 拟合值它承载了线性模型没抓住的非线性结构和周期波动。提取时注意对齐索引训练段和测试段要分开存别把未来信息漏进训练。残差序列通常均值接近 0如果还有明显趋势说明 d 选小了回去重定阶。resid_train train - best[3].fittedvalues resid_test test - best[3].forecast(len(test)) np.save(resid_train.npy, resid_train) np.save(resid_test.npy, resid_test)fittedvalues是样本内一步预测forecast是样本外。两者相减得到的两段残差拼起来才是完整残差序列。存成 npy 是为了下一步 SSA 直接读避免反复跑 ARIMA。到这里ARIMA 这层地基就算打完了接下来处理残差。3. SSA 分解与重构残差里的噪声和周期怎么拆3.1 奇异谱分析的轨迹矩阵与分组逻辑SSA 的核心是把一维序列嵌成轨迹矩阵再做奇异值分解按奇异值大小把分量分成「趋势」「周期」「噪声」几组最后对角平均重构回一维。它不像小波那样需要选基函数也不像 EMD 那样容易模态混叠对残差这种信噪比不高的序列比较友好。关键参数有两个窗口长度 L 和分组方式。L 一般取序列长度的 1/3 到 1/2太短拆不出周期太长计算量陡增且分量之间开始混。分组是 SSA 里最像「玄学」的一步。奇异值从大到小排前几个通常对应趋势和主周期中间一段对应次周期尾部小奇异值基本是噪声。我一般画一张奇异值谱找那个明显拐点拐点之后的重构时直接丢掉。如果拐点不明显就用分组重构后的序列做 ADF 和方差比检验看保留多少组能让残差既平滑又不丢周期。import numpy as np def ssa_decompose(series, L): N len(series) K N - L 1 # 轨迹矩阵 X np.column_stack([series[i:iL] for i in range(K)]) # 奇异值分解 U, s, Vt np.linalg.svd(X, full_matricesFalse) return U, s, Vt, L, K def ssa_reconstruct(U, s, Vt, L, K, groups): groups: 列表的列表每个子列表是一组要保留的分量索引 recon np.zeros(L K - 1) for g in groups: Xg np.zeros((L, K)) for i in g: Xg s[i] * np.outer(U[:, i], Vt[i, :]) # 对角平均 for k in range(L K - 1): idx [j for j in range(K) if 0 k - j L] recon[k] np.mean([Xg[k - j, j] for j in idx]) return recon resid np.load(resid_train.npy) L len(resid) // 3 U, s, Vt, L, K ssa_decompose(resid, L) print(奇异值谱:, np.round(s[:10], 3))ssa_decompose里np.column_stack构造的是 Hankel 型轨迹矩阵每一列是原序列的一个滑动窗口。np.linalg.svd返回的 s 是降序奇异值。ssa_reconstruct里的对角平均是 SSA 反变换的标准操作把矩阵元素沿反对角线取平均还原成一维。groups参数让你自由指定保留哪几组比如[[0,1],[2,3,4]]表示趋势组和周期组分开重构。3.2 分组策略用方差贡献率代替肉眼判断肉眼找拐点不稳定换个人看可能结论就变了。我后来改成算方差贡献率每个奇异值的平方除以总和累计到 85% 到 95% 之间的那一段就是主分量。剩下的尾部直接归入噪声组丢弃。这样每次跑结果可复现也方便写进脚本自动执行。var_ratio s**2 / np.sum(s**2) cum_ratio np.cumsum(var_ratio) keep np.where(cum_ratio 0.90)[0] # 保留累计贡献90%的分量 print(保留分量数:, len(keep), 累计贡献:, round(cum_ratio[keep[-1]], 3)) recon ssa_reconstruct(U, s, Vt, L, K, [list(keep)]) resid_clean resid - recon # 去掉主分量后的噪声残差这里keep取累计贡献不超过 90% 的索引意味着把 90% 的能量当作「有结构的部分」重构出来剩下的 10% 当噪声。resid_clean是去噪后的残差它比原始残差平滑更适合喂给 LSTM。注意recon长度是 LK-1正好等于原序列长度可以直接相减。3.3 重构序列的对齐与长度校验SSA 重构最容易翻车的地方是长度对不上。轨迹矩阵的维度是 L × K对角平均后长度是 L K − 1而 K N − L 1所以 L K − 1 N理论上严丝合缝。但如果你在分组时对某些分量做了截断或补零长度就会漂。我的做法是重构完立刻断言len(recon) len(resid)不等就报错别让它悄悄错位到 LSTM 阶段才暴露。assert len(recon) len(resid), f长度不匹配: {len(recon)} vs {len(resid)} resid_clean resid - recon np.save(resid_clean_train.npy, resid_clean)这个断言看着简单但帮我省过至少两次深夜排查。SSA 输出的resid_clean就是 LSTM 的输入。到这一步线性部分归 ARIMA周期和噪声归 SSA剩下的非线性映射交给 LSTM。4. LSTM 建模把 SSA 重构后的残差喂进网络4.1 滑窗构造监督样本LSTM 要的是监督学习格式得把一维序列切成 (样本, 时间步, 特征) 的三维张量。滑窗长度时间步是个关键参数太短学不到长依赖太长梯度容易消失且样本数骤减。经验上取序列主周期的 1 到 2 倍比如月度数据主周期 12窗口取 12 到 24。步长一般取 1做数据增强时可以取 2 或 3 增加样本多样性。def make_windows(series, window, horizon1): X, y [], [] for i in range(len(series) - window - horizon 1): X.append(series[i:iwindow]) y.append(series[iwindow:iwindowhorizon]) X np.array(X).reshape(-1, window, 1) y np.array(y).reshape(-1, horizon) return X, y resid_clean np.load(resid_clean_train.npy) WINDOW 12 X, y make_windows(resid_clean, WINDOW, horizon1) print(样本形状:, X.shape, y.shape)make_windows里reshape(-1, window, 1)的最后一个 1 是特征维度单变量预测就是 1。horizon1表示单步预测要多步就把 horizon 调大y 会变成二维。注意样本数 序列长度 − window − horizon 1窗口越大样本越少小数据集上要权衡。4.2 网络结构与训练参数结构上我用两层 LSTM 加一层全连接隐藏单元第一层 64、第二层 32dropout 0.2。优化器 Adam学习率 1e-3损失 MSE早停 patience 设 10。批大小 32epoch 上限 200。这些不是金科玉律但作为起点很稳。归一化用 MinMax 把残差缩到 [0,1]因为残差有正有负标准化到均值 0 方差 1 也可以我倾向 MinMax反归一化时不容易出错。import torch import torch.nn as nn class LSTMForecaster(nn.Module): def __init__(self, input_size1, hidden164, hidden232, dropout0.2): super().__init__() self.lstm1 nn.LSTM(input_size, hidden1, batch_firstTrue) self.lstm2 nn.LSTM(hidden1, hidden2, batch_firstTrue) self.drop nn.Dropout(dropout) self.fc nn.Linear(hidden2, 1) def forward(self, x): out, _ self.lstm1(x) out, _ self.lstm2(out) out self.drop(out[:, -1, :]) # 取最后一个时间步 return self.fc(out) model LSTMForecaster() opt torch.optim.Adam(model.parameters(), lr1e-3) loss_fn nn.MSELoss()batch_firstTrue让输入维度是 (batch, seq, feature)符合习惯。out[:, -1, :]只取最后一个时间步的隐藏状态因为预测目标在窗口末尾之后。两层 LSTM 堆叠能提取更抽象的时序模式但层数再多小数据上就开始过拟合了。dropout 加在 LSTM 输出和全连接之间别加在 LSTM 内部PyTorch 的 LSTM 自带 dropout 参数只在多层时生效容易混淆。4.3 训练循环与早停训练循环里我固定做三件事每个 epoch 记录训练和验证损失、验证损失连续 10 轮不降就停、保存验证损失最低的权重。这样即使训练震荡最后拿到的也是验证集上最好的那个模型而不是最后一个 epoch 的。from torch.utils.data import TensorDataset, DataLoader X_t torch.tensor(X, dtypetorch.float32) y_t torch.tensor(y, dtypetorch.float32) split int(len(X_t) * 0.8) train_ds TensorDataset(X_t[:split], y_t[:split]) val_ds TensorDataset(X_t[split:], y_t[split:]) train_dl DataLoader(train_ds, batch_size32, shuffleTrue) val_dl DataLoader(val_ds, batch_size32) best_val, patience, wait float(inf), 10, 0 for epoch in range(200): model.train() for xb, yb in train_dl: opt.zero_grad() loss loss_fn(model(xb), yb) loss.backward() opt.step() model.eval() with torch.no_grad(): val_loss np.mean([loss_fn(model(xb), yb).item() for xb, yb in val_dl]) if val_loss best_val: best_val, wait val_loss, 0 torch.save(model.state_dict(), best_lstm.pt) else: wait 1 if wait patience: print(f早停于 epoch {epoch}) breakshuffleTrue只在训练集开验证集保持顺序。torch.save存的是 state_dict加载时先实例化同结构模型再load_state_dict。早停的 patience 设 10 是经验值数据噪声大可以放宽到 15数据干净可以收紧到 5。训练完把验证集预测和真实值画在一起看是欠拟合还是过拟合再决定要不要调窗口或隐藏单元。4.4 三模型叠加预测与反归一化最终预测 ARIMA 预测 SSA 重构的 LSTM 预测。注意 LSTM 学的是resid_clean而resid_clean resid − recon所以 LSTM 输出的是「去噪后残差」的预测要加回recon再和 ARIMA 相加。顺序别搞反这是最容易出错的一步。model.load_state_dict(torch.load(best_lstm.pt)) model.eval() with torch.no_grad(): lstm_pred model(torch.tensor(X_test, dtypetorch.float32)).numpy().flatten() final_pred arima_pred recon_test lstm_predarima_pred是 ARIMA 的样本外预测recon_test是测试段 SSA 重构的主分量lstm_pred是网络对去噪残差的预测。三者相加才是最终值。如果前面归一化过这里要先反归一化再相加或者全程在归一化空间里加完再反归一化但三种量的归一化参数必须一致否则量纲对不上。5. 避坑与排查这套组合拳最容易翻车的五个地方5.1 现象验证集 MAE 比测试集还低以为模型很强原因滑窗构造时用了未来数据或者归一化参数在划分训练测试之前就用全量数据算了。这是时间序列里最经典的数据泄漏。解决所有统计量均值、方差、最大最小值只在训练段计算然后应用到验证和测试段。滑窗的切分要在归一化之后、但归一化参数必须来自训练段。我习惯写一个fit_scaler(train)和transform(scaler, data)分开的函数从代码结构上杜绝泄漏。5.2 现象SSA 重构后序列和原序列几乎一样LSTM 学不到东西原因分组时保留了太多分量把噪声也重构进去了等于没去噪。或者窗口 L 取太小轨迹矩阵秩不够分解不出有效结构。解决检查累计方差贡献率保留到 85% 到 90% 就停别贪多。L 至少取序列长度的 1/4短序列可以到 1/2。重构后算一下和原残差的相关系数如果高于 0.99说明基本没去掉东西回去调分组。5.3 现象LSTM 训练损失一直降验证损失从第 3 轮就开始涨原因模型容量相对数据量太大或者窗口太长导致样本数太少。小数据集上两层 LSTM 加 64 单元很容易过拟合。解决先减容量隐藏单元降到 32 和 16层数降到一层。再加大 dropout 到 0.3 到 0.4。还不行就缩短窗口增加样本数。数据量少于 500 条时我一般直接用单层 LSTM 加 32 单元起步别一上来就堆结构。5.4 现象最终预测曲线整体偏移所有点都偏高或偏低原因反归一化时用的参数不对或者 ARIMA 预测和 LSTM 预测的量纲没对齐。常见于 ARIMA 用了差分、LSTM 用了归一化两者还原顺序不一致。解决统一在原始量纲上做加法。具体做法是 ARIMA 预测先还原到原始尺度SSA 重构也在原始尺度LSTM 输出反归一化到原始尺度三者都还原完再相加。别在归一化空间里混着加。5.5 现象每次跑结果都不一样MAE 波动超过 20%原因PyTorch 默认随机初始化DataLoader 的 shuffle 也带随机性。深度学习部分本身就有方差。解决固定随机种子torch.manual_seed(42)、np.random.seed(42)DataLoader 的 generator 也设种子。然后跑 5 次取平均作为最终报告值。如果 5 次方差仍然很大说明模型不稳定回去检查数据质量和窗口设置而不是继续调种子。6. 进阶技巧用滚动重训和残差诊断把误差再压一档模型搭完能跑通只是及格线真正拉开差距的是滚动重训和残差诊断。我现在的习惯是每预测一步就把真实值并入训练集重新拟合 ARIMA 和微调 LSTM这叫 walk-forward validation。它比一次性切分训练测试更贴近生产环境因为生产上你每天都会拿到新数据。代价是计算量成倍增加但小数据集上完全扛得住。def walk_forward(y, window, horizon, n_splits): preds, truths [], [] for i in range(n_splits): end len(y) - n_splits i train, test y[:end], y[end:endhorizon] if len(test) horizon: break # 这里复用前面的 ARIMA SSA LSTM 流程 p full_pipeline(train, horizon) preds.append(p) truths.append(test) return np.array(preds), np.array(truths)full_pipeline把前面所有步骤封成一个函数输入训练序列和预测步长输出预测值。n_splits控制回测轮数一般取 5 到 10。跑完算 MAPE 和 RMSE和单次切分的结果对比如果滚动重训明显更好说明数据分布随时间漂移生产上必须用滚动策略。残差诊断是另一个被低估的环节。把最终预测误差画出来看有没有自相关。如果误差的 ACF 在滞后几阶还有显著尖峰说明模型没抓干净还有结构留在误差里。这时候可以回到 SSA 阶段多保留一组分量或者给 LSTM 加一层注意力。我一般用 Ljung-Box 检验p 值大于 0.05 就认为误差接近白噪声可以收手了。from statsmodels.stats.diagnostic import acorr_ljungbox err truths.flatten() - preds.flatten() lb acorr_ljungbox(err, lags[5, 10, 20], return_dfTrue) print(lb)acorr_ljungbox返回每个滞后阶的统计量和 p 值。如果 10 阶 p 值小于 0.05误差还有自相关模型没榨干。这时候别急着换模型先回去看 SSA 分组和 ARIMA 阶数往往问题出在前两层。参数搜索上我最后会做一轮粗网格ARIMA 的 (p,q) 在 0 到 3 之间、SSA 的 L 在 N/4 到 N/2 之间、LSTM 窗口在 6 到 24 之间各取三档组合用滚动验证的 MAPE 排序。组合数控制在 50 以内不然跑一天都跑不完。排序后取前 5 组人工看预测曲线和真实曲线的贴合度选那个不仅误差小、而且拐点也咬得住的。误差小但曲线平滑得像直线的生产上反而危险因为它把突变全抹了。最后说个我自己的习惯每次上线新数据前先拿历史数据做一次「假装不知道未来」的回测把预测值和真实值并排打出来逐月看偏差方向。如果连续三个月同向偏差说明模型有系统性偏移得检查是不是数据采集口径变了而不是模型参数问题。这套流程跑顺之后ARIMA-SSA-LSTM 在小样本单变量预测上确实比单模型稳但它不是万能药数据里的突变如果来自外部事件任何时序模型都抓不住那得靠外部特征。希望帮到你。本文还有配套的精品资源点击获取
