ARTICLE DETAIL

资讯详情

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

小波分解+ARMA组合预测非平稳时间序列:Matlab完整实现

小波分解+ARMA组合预测非平稳时间序列:Matlab完整实现 做时间序列预测这一行说实话大部分精力都耗在“预处理”和“模型选择”的拉扯上。单拿一个ARMA模型来说它对平稳性要求极高而现实中的风速、电价、交通流量、设备振动数据哪一个不是带着趋势和噪声的“刺头”硬拿原始数据去喂ARMA结果往往是残差一堆、预测曲线滞后好几拍。后来我把小波分解加进去做前置处理用Matlab把整个流程跑通之后感觉一下子把这个问题理顺了——先让不同频率的成分各回各家再对分解后的平稳分量分别建模精度和稳定性都有肉眼可见的提升。这篇博文就把我完整的实现思路、参数选择逻辑、Matlab代码细节、以及我踩过的坑写出来。如果你正被非平稳时间序列折磨或者想找一个比纯ARIMA更细腻、比深度学习更轻量的折中方案这篇内容应该能给你省下不少试错时间。1. 项目概述与模型选型思路1.1 核心需求解析这个项目的目标非常明确用“小波分解ARMA”的组合拳完成对非平稳时间序列的预测。标题里“附Matlab代码”说明这不是一个纯理论探讨而是要能落地跑出结果的。先说为什么要拆成两步走。ARMA模型也就是自回归滑动平均模型它的数学假设是序列平稳——均值、方差、自协方差在时间上保持稳定。但真实业务数据几乎不可能满足这个条件。拿电网负荷来说早上和晚上有高峰中午有低谷周末和工作日又是完全不同的用电节奏。这种数据直接建模ACF和PACF图就乱成一锅粥模型参数估计出来也不可靠。小波分解在这里扮演的角色有点像“拆解师”。它能把原始信号拆成不同频率的子序列——低频部分代表趋势和整体形态高频部分代表噪声和局部波动。这样一来原本混在一起的复杂信号变成了几条“性格更单纯”的子序列。对趋势项可以做差分或直接建模型对高频噪声项可以用更灵活的模型或者直接剥离。经过这一层处理ARMA面对的输入就变得乖巧多了。这个方案的适用人群很广。做金融数据分析的、搞设备故障预测的、研究气象风速的、做交通流预测的只要你手里的数据是非平稳的、带噪声的这套流程都值得试一试。特别是那些不想上LSTM、不想碰GPU、希望在一个脚本里快速跑出结果的场景小波加ARMA的组合优势非常明显。1.2 为什么不用纯ARIMA或深度学习每个方案都有它的脾气。纯ARIMA虽然通过差分也能处理非平稳性但它有个致命弱点一次差分处理不了“多尺度”问题。比如一个振动信号既有缓慢漂移的趋势又有叠加在里面的周期性冲击还有高频的测量噪声。差分成I(1)序列只能把趋势抹掉高频噪声还在建模时会被噪声带偏。而且ARIMA对异常值非常敏感一个突发的尖峰能直接影响差分后的序列形态。深度学习方案比如LSTM、GRU确实能端到端解决非平稳问题我在实际项目里也用过。但它的代价是需要大量样本、需要调参、需要GPU加速。样本量少几千条数据训练出来的LSTM可能还没有一个简单的线性模型靠谱。而且LSTM的可解释性很差业务方问你“为什么这个点预测高了”你很难给出清晰的解释。小波分解加ARMA的做法正好卡在中间。小波分解输出的每个分量都有明确的物理含义趋势是趋势波动是波动ARMA又提供了自回归系数和滑动平均系数可以直观看出近期值受历史哪些点影响。这种“可解释性”在工业场景里太重要了。而且整个计算过程在普通CPU上几秒钟就能完成部署成本极低。我后来对比过一组测试数据纯ARIMA预测的MAPE大约在8%左右LSTM在样本充足时能压到6%但小波加ARMA能做到5%上下而且波动更小、跑得更快。当然不同数据特性下结果会有差异但这个组合确实是一个值得优先尝试的强基线方案。2. 小波分解实操要点与参数选择2.1 小波基函数怎么选小波分解的第一步就是要选一个母小波。这个问题我一开始也纠结过后来通过实际测试发现不同的母小波对结果影响很大选错会让分解出的细节分量失真。常见的选择有db2、db4、sym4、coif3等。db系列是Daubechies小波特点是紧支撑、正交matlab里最常用。sym系列是symlet小波比同阶的db小波对称性更好相位失真小一点。coif系列更平滑但计算量稍大。我的选择经验是信号平滑、趋势明显的用db4或db6信号有剧烈突变的用sym4或coif3。Db2虽然计算量小但它的消失矩阶数低对平滑信号的分解效果偏粗糙。既然目标是预测分量质量直接影响后续模型精度没必要省这点计算时间。还要考虑分解层数。层数太少趋势和噪声没有完全分开层数太多低层分量会被过度平滑丢掉有用的高频信息。matlab里wmaxlev函数可以基于信号长度给出最大分解层数建议但我的经验是不要直接拉满设3到5层最实用。选择依据很简单每增加一层低频分量越来越平滑高频分量能量越来越小。当你发现某一层的细节分量能量占比低于1%时那就是极限了再加层只会产生一堆无效的边界效应。2.2 分解层数与重构策略在Matlab里小波分解用的函数是wavedec重构用的是waverec但做预测时我们一般不用直接重构原始信号而是把每个分量单独当成一条时间序列来处理。以三层分解为例% 原始信号 x长度N采样率fs x load(signal.mat); % 或者任何你的数据源 N length(x); % 三层小波分解使用db4小波 [C, L] wavedec(x, 3, db4); % 提取各层系数重构成与原始信号等长的分量 cA3 wrcoef(a, C, L, db4, 3); % 第三层低频近似 cD3 wrcoef(d, C, L, db4, 3); % 第三层高频细节 cD2 wrcoef(d, C, L, db4, 2); % 第二层高频细节 cD1 wrcoef(d, C, L, db4, 1); % 第一层高频细节这里的wrcoef是“waverec reconstruction coefficients”的缩写它能把指定层的系数重构为与原始信号等长的序列。每一层分量相加正好等于原始信号这是小波分解的完美重构性质。我在实操中的做法是低频分量cA3直接进ARMA模型高频分量先看能量占比再决定处理方式。如果cD1那条序列基本接近零轴线上的毛刺那它就是白噪声性质的成分可以单独用一个白噪声模型模拟或者干脆在重构时不加进去。如果是风速数据或者电力负荷数据cD2和cD3往往含有有意义的周期成分需要各建一个模型。分层有个很实用的判断标准观察各分量的自相关图。如果某个高频分量的ACF几乎全部落入置信带内说明它已经接近纯随机可以直接当白噪声处理。反之如果ACF显示出明显的周期性波动这个分量必须建模否则会丢失信息。2.3 边界效应处理心得所有基于卷积的信号处理方法都有边界问题小波分解也不例外。信号两侧因为要和滤波器卷积但滤波器超出信号范围的部分没有数据只能通过延拓方式补齐。matlab里wavedec默认使用周期性延拓但如果你用默认参数往往会在信号两端看到明显的“蝴蝶效应”——分解重构后首尾两端出现异常的波形。处理这个问题有几个办法。最省事的是在Matlab里指定sym延拓方式对称延拓。对称延拓对大多数数据形态都友好首尾不会出现突变。具体设置是[C, L] wavedec(x, 3, db4, mode, sym);但更好的做法是提前对信号做“边缘削波”。我会在分解前把数据头尾各切除一小段等预测结果出来后再映射回真实时间点。比如一段1000点的数据我分解时用第51到第950点预测出的结果前50点和后50点因为是边界区可靠性稍弱我会打上一个置信度降低的标记不让下游调用的人误读。这个细节在实时预测系统里极其关键。电力和气象预测里边界点往往对应刚发生的时刻决策者最关心的就是当前附近时段的预测值。如果你没处理边界效应模型给出的“最新值”和真实值差一大截系统上线第一天就会被投诉。3. ARMA建模与预测核心步骤3.1 平稳性检验与差分处理每个分量分解出来后接下来就是走ARMA建模的标准流程。第一步永远是平稳性检验这一步跳过去后面的定阶和参数估计全都悬在半空中。平稳性检验最常用的是ADF检验Augmented Dickey-Fuller TestMatlab的Econometrics Toolbox里直接调用adftest即可h adftest(cA3); % h 1 表示拒绝单位根假设序列平稳如果h返回0说明序列非平稳需要差分。但注意小波分解出的低频分量往往不是一阶差分就能搞定的。比如趋势项可能是近似线性的一阶差分会引入强烈的负自相关也有可能是缓慢曲线的二阶差分才是合适的。我的经验是每阶差分后都做一次ADF检验直到平稳为止但不要过度差分。过一次差分后要是还不平稳先看一下是不是原序列带明显趋势比如单调增长。如果是想想是不是数据本身的业务逻辑就不需要“去趋势”。比如预测一个缓慢爬坡的电池SOC荷电状态曲线趋势本身就是预测的核心你把它差分掉再预测最后还要积分回来误差反而累积了。这时候可以试试在ARMA模型里加入趋势项或者对低频分量采用“ARMA外生变量”的形式。高阶差分之后的序列信息量会大幅衰减。差分次数越多信噪比越低。数据量小、只有几百个点的时候差分超过两阶预测效果基本就崩了。小波分解的意义在这里再次凸显趋势项被单独分离出来后差分处理变得非常简单直接甚至很多数据一阶差分就完全平稳了。3.2 定阶方法ACF/PACF与AIC准则结合ARMA模型有两个核心参数p自回归阶数和q滑动平均阶数。定阶方法有两条路看自相关图ACF和偏自相关图PACF或者靠信息准则自动搜索。看图定阶这件事新手容易盲目自信。ACF拖尾、PACF截尾定AR模型ACF截尾、PACF拖尾定MA模型都拖尾就ARMA模型。但真实数据里“拖尾”和“截尾”的分界线不明显尤其在样本量少的时候。我的经验是用看图法确定一个大致的范围然后在这个范围内用AIC或BIC自动寻优。Matlab里可以用aicbic函数或者自己写个循环穷举。穷举的过程不复杂maxP 5; maxQ 5; AIC zeros(maxP, maxQ); for p 1:maxP for q 1:maxQ mdl arima(p, 0, q); % 默认不带差分 mdls(p, q) estimate(mdl, cA3, Display, off); [aic, bic] aicbic(logL, numParams); AIC(p, q) aic; end end [minVal, idx] min(AIC(:)); [pBest, qBest] ind2sub(size(AIC), idx);这个循环看起来简单但要注意几个细节。第一estimate函数在序列非平稳时会报错所以必须在循环之前确保cA3已经通过ADF检验。第二AIC的值本身是负的对数似然加上惩罚项所以越小越好但相邻几个p、q组合的AIC差距很小时建议选p、q之和更小的那个更保守的模型避免过拟合。另一个问题是定阶用的数据量不能太少。ARMA模型的参数估计通常用极大似然法样本量低于100时估计出来的参数方差会非常大AIC选出来的阶数也极不稳定。所以如果你手里只有500个数据点分解层数就不要超过4层否则每个分量的样本量被切得只剩一百多个建模基础就薄弱了。3.3 参数估计与残差白噪声检验定阶完成之后模型参数由Matlab直接给出。你会看到AR系数、MA系数、方差等估计值。这些系数能不能用不是看t统计量是不是全显著而是看残差是不是白的。残差白噪声检验是ARMA建模里最容易糊弄过去的一步。很多人估计完参数看一眼拟合图就觉得完事了这是不对的。我常用的检验方式是res infer(mdl, cA3); [h, pValue] lbqtest(res, Lags, 20);lbqtest是Ljung-Box Q检验。如果返回的p值小于0.05说明残差里还有未被模型捕捉的自相关需要调整阶数。这时候别急着改p、q先画一下残差的自相关图看看是哪一阶还在置信带之外。如果第12阶突出说明可能有月度周期的成分没被提取可以考虑加季节性ARMA或者回到小波分解那边增加分解层数。我做故障预测的时候遇到过一个有意思的情况残差在某个特定滞后阶上总是显著反复调p和q都没有改善。后来我画了残差的功率谱发现里面有一个固定频率的振动分量。回头检查小波分解参数发现我用了db4三层分解而数据本身在更低频段还有能量泄漏。换成db6五层分解后这个频率成分被成功分离到一个独立的细节层ARMA残差立刻变白。这就是“分解不够彻底”导致ARMA模型背黑锅的典型案例。4. 完整Matlab流程与关键代码实现4.1 主流程串联整个流程串起来的思路是数据加载 → 小波分解 → 逐个分量建ARMA模型 → 各自预测 → 重构叠加 → 输出结果。下面我给出一份完整的核心代码骨架你可以直接修改数据源后跑通。%% 基于小波分解与ARMA的时间序列预测主流程 % 数据准备x为原始序列step为预测步数 load(yourdata.mat); % x是N×1的double列向量 trainLen floor(0.8 * length(x)); % 80%训练20%测试 xTrain x(1:trainLen); xTest x(trainLen1:end); nPred length(xTest); %% 1. 小波分解 wName db4; nLevel 3; [C, L] wavedec(xTrain, nLevel, wName, mode, sym); % 重构每个分量 cA wrcoef(a, C, L, wName, nLevel); cD zeros(trainLen, nLevel); for k 1:nLevel cD(:, k) wrcoef(d, C, L, wName, k); end %% 2. 逐个分量平稳化并建模 components [cA, cD]; modelPreds zeros(nPred, nLevel1); for i 1:nLevel1 comp components(:, i); % 二阶平稳化检查 if adftest(comp) 0 comp diff(comp); if adftest(comp) 0 comp diff(comp); end end % 用AIC选p、q这里简化为固定值展示流程 mdl arima(2, 0, 2); mdl estimate(mdl, comp, Display, off); % 预测并还原差分 [pred, ~] forecast(mdl, nPred, Y0, comp(end-20:end)); % 若做过差分需要积分还原 modelPreds(:, i) restore_diff(pred, components(:, i), nPred); end %% 3. 叠加各分量预测 yPred zeros(nPred, 1); for i 1:nLevel1 yPred yPred modelPreds(:, i); end %% 4. 评估 mseVal mean((xTest - yPred).^2); mapeVal mean(abs((xTest - yPred)./xTest)) * 100; fprintf(MSE: %.4f, MAPE: %.2f%%\n, mseVal, mapeVal);注意代码里我留了一个restore_diff自定义函数的位置这是处理差分还原的地方。比如一阶差分预测出来的值代表的是变化量要还原成原尺度就把预测值从最后一个真实值开始累加function yRestored restore_diff(pred, original, nPred) lastVal original(end); yRestored zeros(nPred, 1); yRestored(1) lastVal pred(1); for k 2:nPred yRestored(k) yRestored(k-1) pred(k); end end4.2 高频噪声分量的特殊处理前面一直强调不是每个分量都需要建ARMA模型。对于纯噪声级别的高频分量叫cD1或者cD2强行建模只会引入过拟合。我用的判断规则是如果某个高频分量的标准差小于原始信号标准差的2%——这分量直接当零均值白噪声处理预测值为0不参与重构叠加。但你会问预测值为0不会损失信息吗这里要注意小波分解里高频细节分量的能量本来就占比很小对整体预测的贡献微乎其微。如果它里面真的有周期性信息它的标准差不会小成这样。如果标准差小到可以忽略那预测阶段把它当成常量忽略是合理的。还有一种情况高频分量有周期性但周期不稳。比如振动信号受转速波动影响频率有轻微漂移。这种情况ARMA建出来效果也不好我常用的是用Hilbert变换先提取瞬时频率再对瞬时频率建模。但这个已经超出“小波ARMA”的标准范围属于进阶玩法。基础流程里我建议遇到周期性高频分量时先把标准差和周期算出来如果周期稳定就用简单季节项预测周期不稳定才考虑用Hilbert那套方案。4.3 训练测试集划分与滚动预测时间序列预测和普通回归不一样数据不能随机打乱必须保持时间顺序。我通常用两种划分方式固定划分和滚动划分。固定划分简单前80%训练、后20%测试适合离线研究。但真实业务里数据是不断流入的模型需要定期更新。滚动预测的做法是每预测一步就把真实值追加到训练集里重新估计参数再预测下一步。horizon 1; % 单步滚动 for t 1:nPred % 每次用现有全部数据建模型 % 预测下一步 % 把真实值更新到训练集 end滚动预测的计算量更大但效果通常更好。我实测过一个项目固定划分的MAPE是7.2%改成单步滚动后降到5.8%。原因很简单模型能不断吸收最新的状态变化追赶数据的动态漂移。代价是每步都要重新分解、重新定阶、重新估计参数如果是三层分解加ARMA单步耗时大约几十毫秒完全可接受。如果预测步长较长比如24步有两种策略多步递归预测或直接多步预测。ARMA本质上是递归预测型的——预测第2步时把第1步的预测值当输入。这意味着误差会随着步长增加而累积。针对这个问题我会在分解阶段就减少高频分量对多步预测的影响。因为多层递归预测时高频噪声的误差会在每一层放大。低频趋势模型多做几步没问题高频细节模型只预测一到两步后面用零填充。这个“混步长”策略在电力负荷预测里帮我把24步的MAPE压低了将近3个百分点。5. 常见问题与排查技巧实录5.1 分解后分量相加对不上原始数据很多初学的人会遇到三个分量加在一起和原始信号差了一段数值。这通常不是小波函数的问题而是补边方式不一致带来的后果。wavedec默认补边是周期性延拓而你可能又用了sym或者在调用wrcoef时没保持和小波分解一致的延拓模式。注意分解和重构只要用了同一个[C, L]结构wrcoef内部会自动处理对应关系一般不会错。真正会错的是——你把不同模式下分解出来的系数混用或者把一个分量单独用了另一套分解参数。排查方法很简单直接验证完美重构。xRec waverec(C, L, wName); % waverec默认重构全部 max(abs(xRec - xTrain))这个误差应该是1e-10级别的。如果误差很大说明分解参数有问题。另外注意wrcoef和waverec的重构结果首尾会有轻微差异但内部点的误差在浮点精度内。做预测时不要因为这一丁点差异去调参浪费时间。5.2 ARMA建模时报错“Data must be stationary”我用的是Econometrics Toolbox里的arima和estimate组合它对平稳性检查很严格。报这个错通常是差分没做够。但还有一种隐蔽情况差分本身没问题但小波分量包含非零均值而估计器会因为均值的存在而报错。这时候可以先用detrend或者直接减去均值再建模。comp cA3; comp detrend(comp); % 去除均值/线性趋势 if adftest(comp) 0 comp diff(comp); end mdl arima(p, 0, q); mdl estimate(mdl, comp, Display, off);预测完记得把均值加回去。这一步很容易被忽略结果预测曲线整体偏移一个常数MAPE爆炸。如果数据里含有明显异常值传感器掉线导致的零值、脉冲尖刺ARMA的参数估计会被这些点严重拉偏。我的做法是预测前先用中值滤波或用一个简单阈值算法把尖刺标出来替换成插值值。做预测的不是做实时监控预处理阶段稍微花点时间后面模型稳很多。5.3 加了小波分解后效果反而变差这种情况真的存在而且不是个例。你要知道小波分解会延展边界、会产生滤波延迟而且分解后的分量在各层的“相位”可能不完全对齐。如果你数据本身的信噪比很低分解后噪声进入多个分量ARMA在每个分量上都会学到错误的信息叠加起来反而比直接建模更差。我的建议是在任何数据上都要先跑一个“无分解”基线。用ARMA直接对原始数据建模记录指标。然后跑小波分解版本对比再下结论。如果分解版本更差先检查分解层数是不是过多。一个快速有效的方法是计算每一层分量的能量占比画出能量分布柱状图。如果能量分布均匀地散在所有层没有哪一层明显突出说明这个信号不适合用小波分解做层级分离你可以考虑换EEMD集合经验模态分解试试。另外一个容易踩的坑是分解层数定了3层但信号长度只有200个点。第三层低频分量的长度变成了几十个点ARMA在这个长度上建模完全不够看。这种情况我果断把层数降到2层或者干脆只分解一层把高频噪声剥掉就算了。样本量少于等于50个点时不建议再用ARMA直接上指数平滑或均值回归更稳。5.4 p、q阶数自动寻优太慢怎么办AIC穷举法在小数据量下没有问题但数据量稍大、每步滚动预测时反复计算几十个模型的代价也不小。有一个优化思路不要每次都跑全阶数搜索而是只在p、q的一个小邻域内搜索。比如上一轮滚动训练得到的最优参数是p2、q2下一轮训练就在p∈[1,3]、q∈[1,3]这个范围内搜索。因为模型参数在相邻时间窗口内大概率变化不大这种做法不会牺牲多少精度但计算量直接减到原来的1/4。另外Matlab里arima模型的估计有内置优化器迭代次数上限默认比较高。对计算速度要求高的场景可以设置较小的Optionsopts optimoptions(fmincon, Display, off, MaxIterations, 50); mdl estimate(mdl, comp, Display, off, Options, opts);5.5 预测结果的置信区间为什么宽得离谱ARMA模型自带置信区间如果你发现预测区间很宽几乎覆盖了整个测试区间先检查模型是不是过度差分或过度拟合了。过度差分会导致预测区间呈喇叭形扩散而且中心预测值波动巨大。另一个原因是高频分量在预测时被当成随机噪声处理了但它的方差被纳入置信区间计算。如果高频分量本身噪声确实大宽区间是合理的——模型诚实地告诉你这件事它把握不大。但在业务场景里太宽的置信区间没有指导意义。我的处理办法在重构预测值时用低频趋势分量的置信区间作为最终置信区间不把高频噪声的方差叠进去。因为高频噪声本质上不可预测告诉决策者“这部分不可预测”比给一个宽到没有信息量的区间更有价值。这个方法让我在很多项目汇报里保住了模型的“可靠性”人设——毕竟置信区间窄业务方就会更信任你。写在最后的一点实操体会小波分解加ARMA这个组合我用了三年多从电网负荷到电池SOC预测都跑过。最初我以为它就是“预处理加了个滤镜”后来才发现它的价值不只是滤噪——它让你能把不同时间尺度的规律分开研究。趋势是什么节奏、波动是什么节奏、噪声是什么量级一目了然。这种分层视角比纯粹调box-jenkins参数更能帮你理解数据本身。如果后续想继续扩展可以在ARMA这一层换成其他模型。比如对低频趋势分量用多输出回归对高频周期分量用含外生变量的ARIMAX甚至对某一层用bilstm解析更复杂的非线性结构。小波分解给出的分量天然适合做这种“一器一模型”的分配这也是它比“一个模型吃所有”更先进的原因。你在实际项目里遇到特别难啃的数据不妨先分解看层再决定每一步怎么处理。先让小波把问题分层分清楚后面的事情能顺很多。
返回列表