ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

SSA-ARIMA-LSTM组合模型实战:Python实现时间序列预测

SSA-ARIMA-LSTM组合模型实战:Python实现时间序列预测 简介Python实现的ARIMA-SSA-LSTM组合时间序列预测完整源码与数据包面向计算机、电子信息、数学等专业学生的课程设计、期末大作业与毕业设计场景也适合刚接触深度学习与时间序列预测的开发者。资源压缩包共3个文件包含2个CSV数据集和1个Python主程序整体大小仅51KB轻量易部署。代码采用参数化编程参数集中配置、调整方便并几乎一行一注释保姆级注释可帮助初学者逐行理解数据预处理、模型构建、训练与预测的完整流程。数据文件针对焦作相关指标整理可直接运行验证并复现预测结果同时便于替换成自己的数据开展对比实验。目前已有509人学习下载这套源码兼顾完整性与易读性既能满足课程设计、毕业设计的完整方案需求也是学习ARIMA、SSA、LSTM组合预测模型的实用参考。1. 时间序列预测的配合战ARIMA、SSA 与 LSTM 各自解决什么单一模型做时间序列预测总是顾此失彼。ARIMA 线性趋势抓得稳碰到非线性波动就露怯LSTM 什么模式都能拟合但在序列短、信噪比低时会把噪声也记住SSA 能把序列拆成趋势、周期和噪声但拆完并不直接产出预测。水文径流、电力负荷、销量这类数据恰好是三种成分混在一起所以工程上越来越多地把它仨串成流水线SSA 负责信号分离ARIMA 负责线性主趋势LSTM 拾取非线性残差。下面用 Python 从零实现这条链路覆盖 SSA 分解、ARIMA 定阶、LSTM 训练与预测融合代码均可直接改用真实数据运行。2. 先做信号分离SSA 把趋势、周期与噪声拆成可预测成分时间序列预测的第一步往往不是选模型而是先看清数据里到底混了几种成分。SSA奇异谱分析的定位就是信号预处理器它不假设序列服从某种分布也不要求周期固定只凭矩阵分解就能把序列拆成若干个可解释的分量。对于长度足够、周期相对稳定的序列拆开后通常能直接读出三类东西缓慢变化的主趋势、几个能量递减的周期分量以及散乱的噪声。后面 ARIMA 和 LSTM 各自只处理其中一部分模型压力会小很多。同类工具里经验模态分解EMD和小波分解也常被拿来拆序列但 EMD 的模态混叠问题要靠集合平均去缓解小波则需要预先选好基函数和分解层数。SSA 的优势在于只有两个超参数窗口长度 L 和保留分量数 r且分解结果对参数扰动相对平缓工程上更容易调稳。这也是在 ARIMA-LSTM 组合之前先接 SSA 的最直接理由。2.1 轨迹矩阵、SVD 与对角平均SSA 的四步骨架SSA 的流程固定为四步嵌入、奇异值分解、分组和对角平均。第一步把一维序列“铺”成二维轨迹矩阵固定窗口长度 L从序列开头依次滑动取值得到形状为 (L, K) 的矩阵其中 K n − L 1。第二步对轨迹矩阵做 SVD得到按大小排列的奇异值大奇异值对应高能量成分。第三步分组把前 r 个奇异值分量划入信号其余归为噪声。第四步沿矩阵反对角线取平均把每个秩一矩阵还原成与原序列等长的一维分量。第四步最容易实现出错因为轨迹矩阵对角线上的元素个数不一致矩阵中间的对角线有 L 个元素两端的元素依次递减。用双层循环也能算但边界条件容易写错。用 bincount 一趟完成累加和计数边界问题交给函数自己处理。import numpy as np def ssa_decompose(series: np.ndarray, L: int): SSA 分解返回分量序列与奇异值。 series : 一维原始序列 L : 窗口长度嵌入维数需满足 1 L n/2 n len(series) K n - L 1 # 1) 嵌入轨迹矩阵 X形状 (L, K) X np.stack([series[i:i L] for i in range(K)], axis1) # 2) SVD U, s, Vh np.linalg.svd(X) # 3) 4) 每个奇异值对应一个秩一矩阵对角平均还原为长度 n 的序列 idx np.arange(L)[:, None] np.arange(K)[None, :] counts np.bincount(idx.ravel(), minlengthn) components np.zeros((L, n)) for i in range(L): Ei s[i] * np.outer(U[:, i], Vh[i, :]) components[i] np.bincount(idx.ravel(), weightsEi.ravel(), minlengthn) / counts return components, s这个函数把四步压缩成三个关键操作np.stack 完成嵌入np.linalg.svd 完成分解bincount 完成对角平均。L 是唯一的必填参数返回值 components 里每一行是一条长度为 n 的分量序列sv 是按大小排列的奇异值序列下一步选 r 时要用到它。2.2 窗口长度 L 和保留分量 r 怎么定L 和 r 是 SSA 仅有的两个超参数它们共同决定了“信号”和“噪声”的边界。先说 L。它表征嵌入维数直接控制了轨迹矩阵能携带多少时间延迟信息L 太小时信息量不足周期成分分不开L 太大时矩阵接近满秩SVD 得到的奇异值分布变得均匀信号和噪声又糊在一起。实践上 L 落在 n/4 到 n/2 之间比较稳若已知序列有明确周期优先取一个完整周期或周期倍数。场景L 的常用取值说明月度数据带年周期12 或 24对齐周期后周期成分集中在前几个奇异值日度数据带周周期7 或 14周内波动是主要周期成分无明显周期n/4 ~ n/2按数据长度试探看奇异值谱找拐点序列过短n 100不超过 n/2窗口过大导致轨迹矩阵近似满秩分解失效r 的选取看奇异值谱。把奇异值平方后归一化看前 r 个分量的累计能量占比超过 95% 即可截断。这个阈值允许浮动数据噪声大就放宽到 90%周期结构强就收紧到 98%。经验上径流、负荷这类序列前 35 个分量就能装下主趋势和一个主周期奇异值谱会出现明显拐点拐点之后的分量能量占比断崖式下跌基本是噪声。注意r 宁可取大不要取小。多保留一个低能量分量最多是让 ARIMA 多拟合一点周期尾巴少保留一个高能量分量趋势被截断后整条流水线的预测都会带上系统性偏差。2.3 用代码把信号和噪声拆出来数据文件建议整理成两列 CSVdate 是日期或序号value 是观测值。缺失值先插值异常值先剔除不要直接填 0时间序列模型对断点非常敏感。读取后先调用 ssa_decompose再按累计能量占比自动选 r。import pandas as pd # 数据格式两列 CSVdate 为日期或序号value 为观测值 df pd.read_csv(series.csv, parse_dates[date]) y df[value].values.astype(float) # SSA 分解L 取一个周期长度 12 comps, sv ssa_decompose(y, L12) # 按累计能量占比自动选择信号分量数 r energy sv ** 2 cum np.cumsum(energy) / energy.sum() r int(np.argmax(cum 0.95)) 1 # 信号 前 r 个分量之和噪声 原始序列 - 信号 signal comps[:r].sum(axis0) noise y - signal print(fsignal energy ratio: {cum[r-1]:.3f}, r {r}) print(fnoise mean: {noise.mean():.4f}, noise std: {noise.std():.4f})分解效果要从两个角度核对signal 曲线要平滑且形状和原始序列大体一致noise 的均值要接近 0如果噪声里还看得出周期性说明 r 偏小回去把阈值从 0.95 提到 0.98 再试一轮。这一步调稳了后面 ARIMA 和 LSTM 的参数量可以维持得很小。3. 线性主趋势交给 ARIMA平稳性检验、定阶与前向预测SSA 拆出的 signal 是平滑的主序列噪声已被剥离。这个序列的特点是自相关结构简单、周期性明确正是 ARIMA(p, d, q) 的舒适区。三个参数各有分工p 是自回归阶数描述“当前值由过去 p 个值解释多少”d 是差分阶数用来消除趋势带来的非平稳q 是滑动平均阶数描述“当前值由过去 q 个预测误差修正多少”。d 阶差分之后拿到的是平稳残差AR 和 MA 项只有在平稳序列上才有统计意义。3.1 d 阶差分与 ADF 检验先确认序列是否平稳d 的确定属于统计检验问题。ADF 检验的原假设是序列存在单位根p 值小于 0.05 时拒绝原假设认为序列平稳。实际代码里从 d0 开始逐次差分重测直到平稳或达到最大阶数。d 上限控制在 2超过 2 的差分不仅损失样本量还会把周期性信息一起差掉。from statsmodels.tsa.stattools import adfuller def find_d(signal, max_d2): 逐阶差分做 ADF 检验返回使序列平稳的最小 d。 s signal.copy() d 0 p_val adfuller(s)[1] while p_val 0.05 and d max_d: s np.diff(s) d 1 p_val adfuller(s)[1] return d注意差分只用于“判定阶数”真正建模时不手动差分化数据而是把 d 原样传给 ARIMA 的 order 参数。statsmodels 内部会对原始序列完成差分和差分回退预测值因此保持在原始量纲上。手动差分后再喂数据预测结果会整体偏移一个常数多步预测时误差被逐步放大。3.2 p、q 怎么定ACF/PACF 先看、AIC 网格再算定 p、q 有两条路。第一条是看图计算差分后序列的 ACF 和 PACF根据截尾、拖尾特征做初判第二条是算在 p、q 的小范围内网格搜索以 AIC 最小为准则。真实数据的相关图通常都不干净我更信任网格搜索图只用来限定搜索范围。观察到的自相关形态模型倾向原因PACF 在 p 阶截尾ACF 拖尾AR(p)过去 p 个值足以解释当前值ACF 在 q 阶截尾PACF 拖尾MA(q)当前值是前 q 个误差的线性组合两者都拖尾衰减ARMA(p,q)需要同时引入 AR 和 MA 项两者都快速截尾近似白噪声signal 分解过度趋势被拆没了网格搜索有两个细节p0 且 q0 的组合对应纯白噪声直接跳过某些 (p, q) 组合在短序列上会收敛失败用 try/except 捕获不要让一个坏组合拖垮整个搜索。signal 比原始序列平滑p、q 范围给到 04 通常足够。import itertools from statsmodels.tsa.arima.model import ARIMA def select_order(signal, d, max_p4, max_q4): AIC 网格搜索返回最优 (p, d, q)。 best_aic, best_order float(inf), None for p, q in itertools.product(range(max_p 1), range(max_q 1)): if p 0 and q 0: continue try: fit ARIMA(signal, order(p, d, q)).fit() if fit.aic best_aic: best_aic, best_order fit.aic, (p, d, q) except Exception: continue return best_orderAIC 是拟合优度和参数数量的权衡它不会直接给出“最好”的物理解释但在候选模型里选一个稳健够用的网格搜索比人眼盯图可靠得多。3.3 用最优阶数拟合 signal 并做多步预测定阶完成后用 statsmodels.tsa.arima.model 里的 ARIMA 类做最终拟合。statsmodels 0.12 之后新增的 ARIMA 类在接口和内部实现上都与旧版 ARIMA 模块不同跑代码前先确认 import 路径和版本。拟合后调用 forecast(steps) 得到预测序列返回值是原始量纲不需要手动逆差分。拟合之后要检查残差是否还有线性信息残留。用 Ljung-Box 检验原始残差p 值小于 0.05 说明残差仍存在显著自相关。此时先回到 SSA 那一步把 r 调大而不是盲目加大 p、q——信号分量没吃干净时ARIMA 再怎么调阶都只是在一个残缺信号上做文章。def arima_forecast(signal, d, steps): 按 AIC 网格选阶并预测 steps 步。返回预测序列与最优阶数。 order select_order(signal, d) fit ARIMA(signal, orderorder).fit() return fit.forecast(stepssteps), order, fit整个 ARIMA 环节有一条原则拟合对象永远是 SSA 拆出的 signal不是原始序列。后面与 LSTM 融合时ARIMA 预测信号部分LSTM 预测噪声部分两边不同量纲时相加没有意义。保持各自独立训练、最后同尺度相加是这套流水线的核心约定。4. 非线性残差交给 LSTM构窗口、定结构、控训练SSA 拆出的 noise 序列仍有内容里面是缓慢变化的非线性波动、突发脉冲和测量误差。ARIMA 的线性假设在这里基本失效硬套只会得到一条接近零的预测。LSTM 这类循环网络正好接盘它用门结构控制信息在时间步之间的流动遗忘门决定上一时刻记忆保留多少输入门写入新信息输出门决定放出什么天然适合“用过去一段窗口预测下一个值”的任务。做 LSTM 时间序列预测 Python 实现通常也是从两层 LSTM 加 Dropout 的标配结构起步。4.1 残差序列切片把一维数据变成监督学习样本LSTM 不能直接吃一维序列输入必须切成 (样本数, lookback, 特征数) 的三维张量。lookback 是“看多远”物理意义是记忆窗口通常取一个周期长度周期不明显时从 512 里试。切片前要做归一化LSTM 内部激活函数有界输入尺度跨度过大时梯度不稳定训练很难收敛。MinMaxScaler 把数据压到 [0, 1]并且只用训练段拟合 scaler验证段和测试段复用同一套 min/max这是避免信息泄漏的基本纪律。from sklearn.preprocessing import MinMaxScaler def make_windows(data, lookback): 滑动窗口切分用前 lookback 个点预测下一个点。 X, y [], [] for i in range(len(data) - lookback): X.append(data[i:i lookback]) y.append(data[i lookback]) return np.array(X), np.array(y) # 先按时间切分再做归一化防止测试集信息进入训练 split_idx int(len(noise) * 0.8) train_noise, test_noise noise[:split_idx], noise[split_idx:] scaler MinMaxScaler() scaler.fit(train_noise.reshape(-1, 1)) train_scaled scaler.transform(train_noise.reshape(-1, 1)).ravel() test_scaled scaler.transform(test_noise.reshape(-1, 1)).ravel() lookback 12 X_train, y_train make_windows(train_scaled, lookback) X_test, y_test make_windows(test_scaled, lookback) X_train X_train.reshape(X_train.shape[0], X_train.shape[1], 1) X_test X_test.reshape(X_test.shape[0], X_test.shape[1], 1)切片后样本数等于 len(noise) − lookback每多 1 个 lookback 就少 1 个样本。序列短几百个点以内时优先减小 lookback而不是压缩测试集比例测试集至少要留 20%否则指标置信度太低。4.2 网络结构两层 LSTM 加 Dropout 的默认配置第一层 LSTM 设置 return_sequencesTrue把每个时间步的隐状态都传给下一层第二层 return_sequencesFalse只输出最后一个时间步的隐状态再接 Dense(1) 产出预测值。units 的取值范围在 32128太小欠拟合太大会在噪声成分上过度记忆。Dropout 放在每层 LSTM 之后随机丢弃一部分神经元连接是控制过拟合最直接的手段。Keras 里第一层显式传入 input_shape后面的层会自动推断维度。from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout from tensorflow.keras.callbacks import EarlyStopping model Sequential([ LSTM(64, return_sequencesTrue, input_shape(lookback, 1)), Dropout(0.2), LSTM(32, return_sequencesFalse), Dropout(0.2), Dense(1) ]) model.compile(optimizeradam, lossmse, metrics[mae])损失函数用均方误差优化器选 adam两者都是时间序列回归的默认选项。监控指标加 mae 是为了对照参考训练过程主要盯 val_loss不是训练集上的 loss。4.3 训练超参早停、批大小和平移预测的识别训练阶段最该盯的是验证集 loss 曲线。数据量小时训练 loss 持续下降而验证 loss 回升是典型的过拟合信号。EarlyStopping 监听 val_losspatience 给 1015 个 epoch 的缓冲验证 loss 不再下降时自动回滚到最优权重。batch_size 决定梯度估计的稳定性噪声数据上用 1632 的小批量更容易跳出局部最优代价是训练时间变长。判断 LSTM 有没有真正学到模式最实用的办法是把预测值和真实值画在一起看如果预测曲线像把真实曲线整体右移了一个 lookback 窗口说明模型学会了“复制最近的输入”而不是提取模式。出现这种平移预测时先确认输入特征切的是 noise 而不是原始序列再尝试减小 units、加大 Dropout。early_stop EarlyStopping( monitorval_loss, patience15, restore_best_weightsTrue ) history model.fit( X_train, y_train, validation_split0.1, # 从训练集再划 10% 做验证 epochs200, batch_size32, callbacks[early_stop], verbose1 )超参数常用范围调整方向lookback13 个周期长度周期不明显时取小值序列短时优先减小units32128验证 loss 不平滑就加大出现过拟合就减小dropout0.10.3噪声成分大时加大batch_size1664数据量小时用 16训练太慢再加大patience1020序列波动大时适当加大避免太早停5. 预测融合落地评价指标与 SSA-LSTM 组合的三个常见坑5.1 多步预测拼接与指标对照融合前先统一步数。ARIMA 用 forecast(stepshorizon) 一次给出 horizon 步LSTM 做多步递归把最后一个 lookback 窗口作为种子每预测一步就把新值回填进窗口逐点滚动。def recursive_lstm_predict(model, scaler, seed_window, steps): 用最后一个 lookback 窗口做种子递归多步预测。 window seed_window.reshape(1, lookback, 1).astype(float32) results [] for _ in range(steps): nxt model.predict(window, verbose0)[0, 0] results.append(nxt) window np.concatenate( [window[:, 1:, :], np.array([[[float(nxt)]]], dtypefloat32)], axis1 ) return scaler.inverse_transform(np.array(results).reshape(-1, 1)).ravel() horizon 12 arima_pred, _, _ arima_forecast(signal, dfind_d(signal), stepshorizon) lstm_pred recursive_lstm_predict(model, scaler, test_scaled[-lookback:], horizon) final_pred arima_pred lstm_pred # 两部分都是原始量纲直接相加指标上同时看 RMSE 和 MAPERMSE 对大的偏差敏感能暴露模型在某一段上的失控MAPE 给出相对误差水平便于跨数据集对比。评估要做三路对照单独 ARIMA、单独 LSTM、组合模型各跑同一测试段。模型组合平稳段signal 主导波动段noise 主导实现成本单独 ARIMA好差低单独 LSTM中好高SSA ARIMA LSTM好好中5.2 三个容易出现偏差的地方第一个坑是 SSA 参数泄漏。正式回测时L 和 r 必须由训练段数据决定再用同一组参数处理测试段先对全序列分解再切分等于让模型提前看到了未来信息回测指标被系统性美化。严格的做法是滚动窗口每一轮只取当前窗口内的数据重新做 SSA 定参、ARIMA 定阶、LSTM 训练再预测下一步。第二个坑是 LSTM 陷入“复制平移”的局部最优特征切错序列是高频原因验证集曲线比训练集曲线更值得盯。第三个坑是把高频分量一律当噪声丢掉径流日数据里的汛期脉冲、销量数据里的突发促销往往藏在较高阶分量里r 宁可多留一个也不要一刀切。每次融合预测后把 signal 预测、noise 预测和最终预测三条曲线画在同一张图上若 final 与 signal 几乎重合说明 LSTM 本次没有贡献优先复查 lookback 和 r若 final 在波峰处系统性偏低多半是 LSTM 递归预测累积了误差缩短预测步数或用直接多步输出都可以缓解。本文还有配套的精品资源点击获取
返回列表