写这篇博文之前先说个实际感受很多做时间序列预测的朋友拿到非平稳数据直接上ARMA拟合出来R²倒是挺好看但一到预测段就漂得厉害。我自己踩过这个坑之后才认真把“小波分解ARMA”这套组合捡起来。这篇文章就围绕“基于小波分解和ARMA预测附Matlab代码”这个项目把原理、选型、完整代码和调试经验一次性讲透。内容偏实操代码可以直接抄但每个关键步骤我都会补一段“为什么这么做”的解释。1. 项目背景与整体思路拆解1.1 为什么单一ARMA模型搞不定非平稳数据ARMA模型的全称是自回归移动平均模型核心假设是序列平稳。所谓平稳粗略理解就是均值、方差在时间轴上不剧烈变化自相关结构只与时间间隔有关。现实中你拿到的数据无论是电力负荷、风速、交通流量还是设备振动信号几乎没有一路是天然平稳的。趋势项、周期项、随机冲击叠加在一起ARMA直接硬套残差里依然藏着强自相关结构参数估计的方差被严重低估预测区间完全失真。我当初用纯ARMA预测一组带明显季节性的负荷数据训练段拟合得停不下来滚动预测却只能在均值附近荡秋千。后来才意识到问题不在于模型阶数定得不够好而是数据本身含有多个不同频段的成分单一模型压根不具备同时刻画趋势和局部波动结构的能力。这个场景下用小波分解先把序列按频率拆开再分别用ARMA对各分量建模是一种比“整体硬建模”要稳妥得多的思路。1.2 小波分解在这套方案里扮演什么角色小波分解的核心作用可以理解成“分拣机”。原始信号是混在一起的混合体小波变换通过母小波的平移和伸缩把信号分解成若干不同频带的细节分量高频成分和一个近似分量低频趋势。对预测任务来说这一步的真正价值是原本非平稳的整体序列被拆成若干平稳性明显改善的子序列。近似分量虽然还带有缓慢趋势但它比原始序列的光滑度高很多各层细节分量则在零均值附近波动符合平稳假设的程度大幅提升。实际工程里常听到的EMD经验模态分解也能做类似的事情但EMD有模态混叠问题和端点效应分解结果依赖停止条件复现性差一些。小波分解有严格的数学框架分解层数可控而且通过逆变换可以无损重构原始信号。对需要稳定产出、重复实验的项目来说小波是更省心的选择。1.3 组合方案的架构设计整套方案的设计思路其实就三步分解、建模、重构。第一步用wavedec对训练数据做N层小波分解得到近似分量A_N和细节分量D_1到D_N。第二步对逐个分量建立ARMA模型这里要独立定阶、独立拟合并滚动预测。第三步把各分量预测结果直接相加还原成原始尺度上的预测值。有人会问为什么不直接对整体数据做ARMA模型反而绕这一圈原因是组合模型把“趋势跟随”和“波动捕获”两项任务分给不同类型的子模型去承担。近似分量负责趋势外推细节分量负责局部周期波动。分工明确参数解释也更清晰。实测中如果数据本身存在明显的多尺度特征这套组合通常比单一ARMA预测误差低20%到40%尤其在多步预测场景下提升更明显。2. 核心理论知识精讲小波基选择与ARMA定阶2.1 小波基怎么选db4为什么是默认选项Matlab的wavedec函数需要一个指定小波基的参数一般写成db4、sym4、coif3这样。理论上有无限多种小波基可选但工程实践中大多数情况下db4足够胜任。db系列是Daubechies小波紧凑支撑且正交能量集中度比较好。db4表示消失矩为4意思是可以精确逼近到三次多项式的信号趋势。对绝大多数负荷、风速、振动类信号来说db4的时频局部化能力与计算复杂度之间平衡得很好。如果信号形态有较强的对称性可以考虑sym4相位失真更小。如果信号本身包含较多尖锐突变或者你更关注奇异点位置可以试试haar小波它是最简单的分解但平滑性太差一般不做预测首选。选小波基还有一条经验不一定追求越高阶越好。db8的频带分割更锐利但边界效应也更明显尤其在数据长度有限的情况下高阶小波会产生明显的端点振荡反而污染预测精度。2.2 分解层数定多少三层四层还是更多分解层数直接影响分量数量和每个分量的频带宽度。层数太少分离不彻底层数太多高频细节分量会变成近乎纯噪声ARMA建不出有效模型。常用的经验规则是数据长度在几百到一两千之间分解层数取3到5层。具体可以观察近似分量曲线是否已经足够平滑如果A_5依然还有明显的周期波动说明该继续分解。也可以借助小波包能量占比来做层数选择但预测场景下不必搞这么复杂。我做负荷预测时通常先用三层分解近似分量A_3细节分量D_1最高频、D_2中高频、D_3较低频。A_3代表长期趋势和主要周期D_3体现短周期波动D_1和D_2则对应随机噪声和高频抖动。这样一个结构对大多数业务场景来说足够解释了。2.3 ARMA模型定阶别只盯AICARMA(p,q)的定阶是另一个容易翻车的点。常用方法有两种看ACF自相关函数和PACF偏自相关函数的截尾拖尾特性或者用AIC、BIC信息准则自动搜索。实际经验是信息准则给出的是统计意义上的“优”不是预测意义上的“好”。AIC偏低容易选过大的阶数处理不好就是过拟合。分量序列的样本量不大时AIC选出来的高阶模型往往在线外预测中表现惨淡。我建议的做法是用aicbic函数算一组候选阶数然后结合分量序列的实际特性取较小阶数。近似分量可以用稍高一点的阶数比如p在3到5之间细节分量尤其是高频分量阶数控制在1到2就好甚至可以退化为MA模型因为高频分量本身记忆性弱。还有个细节Matlab的armax函数输出的是离散时间多项式模型不是你手动定义的ARMA系数。写代码时最好统一用idpoly或者直接配合estimate函数避免在模型格式转换上浪费时间。3. 完整Matlab实现代码与逐步解析3.1 代码整体框架一览先把主要流程列出来后面逐段解释。%% 清空环境 clear; clc; close all; %% 生成模拟数据趋势 周期 噪声 rng(42); t (1:1000); data 0.05 * t 8 * sin(2 * pi * t / 120) 3 * sin(2 * pi * t / 35) randn(1000, 1) * 2; train_len 800; test_len 200; train_data data(1:train_len); test_data data(train_len1:end); %% 参数设置 wname db4; level 3; horizon test_len; % 预测步长模拟数据就按“线性趋势 两个不同周期的正弦 高斯噪声”叠加这样可以验证分解能否有效拆出不同频段成分。训练集八百个点测试集两百个点预测步长设成与测试集等长算是比较苛刻的多步预测场景了。3.2 小波分解wavedec与wrcoef的配合%% 小波分解 [C, L] wavedec(train_data, level, wname); A3 wrcoef(a, C, L, wname, level); D3 wrcoef(d, C, L, wname, level); D2 wrcoef(d, C, L, wname, level - 1); D1 wrcoef(d, C, L, wname, level - 2);wavedec返回的是小波分解结构C和长度记录L。C里装的是所有层的小波系数和最后一层尺度系数L记录了每一段系数的长度。直接拿C去建模没有物理意义必须用wrcoef按层重构回时域信号。wrcoef的a表示重构近似分量d表示重构细节分量。一个常见的错误是直接把C当成分量序列丢给ARMA模型这绝对不行。C是系数域不是时域。wrcoef重构回来才是与原始信号等长的分量序列。3.3 各分量分别建模循环里的自主定阶每个分量都要单独定阶、建模。为了演示我写了一个简单的定阶函数基于AIC同时限制阶数上限避免过拟合。%% 各分量ARMA建模函数 function model fitarma(series, maxp, maxq) best_aic Inf; best_model []; best_p 0; best_q 0; data_id iddata(series, [], 1); for p 0:maxp for q 0:maxq if p 0 q 0 continue; end try M armax(data_id, [p q]); [~, aic_val] aicbic(M); if aic_val best_aic best_aic aic_val; best_model M; best_p p; best_q q; end catch continue; end end end model best_model; end调用时对不同分量给不同的阶数上限%% 对近似分量和细节分量分别建模 [A_model] fitarma(A3, 5, 3); [D3_model] fitarma(D3, 3, 2); [D2_model] fitarma(D2, 2, 2); [D1_model] fitarma(D1, 1, 1);近似分量给5阶上限细节分量从3降到1高频分量限得最紧。我建议你实际跑数据时也这样处理高频分量模型越简单越稳硬要让它记住每一个毛刺预测效果反而像抽签。3.4 预测与重构forecast的维度陷阱%% 预测未来horizon个时间点 f_A forecast(A_model, A3, horizon); f_D3 forecast(D3_model, D3, horizon); f_D2 forecast(D2_model, D2, horizon); f_D1 forecast(D1_model, D1, horizon); f_A f_A(:); f_D3 f_D3(:); f_D2 f_D2(:); f_D1 f_D1(:); forecast_total f_A f_D3 f_D2 f_D1;forecast函数的输入第一个参数是模型对象第二个是历史数据第三个是预测步长。它输出的是一个列向量如果直接和行向量做加法会出现维度不匹配或隐式扩展代码里统一转置成行向量再累加。这里有一个容易被忽视的维度陷阱小波分解后重构出的分量序列长度等于原始序列长度。forecast对每个分量预测horizon个点得到的是各个分量未来段的估计。整体预测就是这几个分量预测值逐点相加。这个加法操作的意义是重构和wavedec无关。3.5 误差评估与可视化%% 误差评估 actual test_data; mae mean(abs(forecast_total - actual)); rmse sqrt(mean((forecast_total - actual).^2)); mape mean(abs((forecast_total - actual) ./ actual)) * 100; fprintf(MAE: %.3f\n, mae); fprintf(RMSE: %.3f\n, rmse); fprintf(MAPE: %.2f%%\n, mape); %% 绘图对比 figure; plot(1:test_len, actual, b-, LineWidth, 1.5); hold on; plot(1:test_len, forecast_total, r--, LineWidth, 1.5); legend(真实值, 小波分解ARMA预测); xlabel(时间步长); ylabel(数值); title(预测结果对比); grid on; %% 分解效果可视化 figure; subplot(5,1,1); plot(train_data); title(原始训练序列); subplot(5,1,2); plot(A3); title(近似分量 A3); subplot(5,1,3); plot(D3); title(细节分量 D3); subplot(5,1,4); plot(D2); title(细节分量 D2); subplot(5,1,5); plot(D1); title(细节分量 D1);误差评价用MAE、RMSE、MAPE三个指标一起看。MAPE对接近零的真实值特别敏感如果数据里存在接近零的采样点MAPE会显得很大此时建议以RMSE为主参考。做分解效果可视化有一个实际价值你能直观看到哪些分量包含了主要能量哪些分量几乎是噪声。如果D1基本像白噪声对预测的贡献很小建模时甚至可以省掉它只对A3、D2、D3建模。4. 常见问题与调试经验实录4.1 边界效应导致预测前段偏差过大小波分解本质上是对有限长序列做卷积变换在序列两端不可避免地出现边界效应。ARMA预测时边界效应主要影响的是重构出来的最后一个点附近的分量值而这个值恰好是滚动预测的起点。解决思路有三条第一分解前对数据两端做对称延拓延拓长度建议不低于小波支撑长度。预测完成后再截掉延拓部分。这个方法有效但会增加代码复杂度。第二换成sym小波基来缓解相位失真对边界问题有一定改善但不彻底。第三训练段和测试段之间留一段缓冲带比如训练集只用到800个点预测起点从801开始但实际业务上有个固定的预测启动点。简单说预测起点附近如果存在剧烈的边界振荡可以考虑给预测结果前几个点一个小幅修正或者忽略前几个点的误差统计。我在实际项目中测试过边界效应最严重的是高频分量D1和D2。如果这两个分量的预测值异常大检查一下重构后的端点值是不是偏离了正常范围可以用前后几个点的均值来平滑掉这个异常。4.2 高频分量预测误差拖累整体结果这是最常遇到的坑也是小波加ARMA方案最致命的弱点。D1代表最高频成分本质上就是随机噪声的提取物。你拿ARMA去拟合噪声训练集内可以拟合出一些参数但这些参数描述的是噪声的随机实现不是规律。预测时它不仅没有预测能力还会引入额外误差。我踩坑后总结的做法有三种第一种当D1的方差远小于其他分量时直接放弃预测它把它的预测值设为零。D1的均值本来就接近零预测为零是合理的期望。第二种给D1、D2建模时强制限制低阶比如D1用ARMA(1,0)或者干脆用AR(1)D2用ARMA(1,1)。模型只捕捉最基础的短期相关性不追求训练集内拟合优度。第三种在重构前将高频分量乘以一个衰减系数比如0.5。这个系数相当于软阈值降噪可以减少高频误差对整体的干扰。但系数需要根据验证集调不是固定的。4.3 forecast函数报错或被拒绝Invalid forecast horizon. Must be a positive scalar.这个问题多发生在horizon被错误定义成向量或者0时。检查一下horizon是否为正整数另外确认历史数据是列向量还是行向量。Matlab的forecast函数对维度要求比较严格历史数据是行向量时某些版本会报错。统一转成列向量输出再转置回来。另一个常见问题是模型对象格式不正确。armax返回的对象可以直接用forecast但如果你自己手动构建了idpoly模型要先检查属性是否完整最省事的方式还是直接对iddata对象用armax一步到位。4.4 分解层数与模型稳定性之间的平衡分解层数不是越大越好。层数增加后A_N会越来越平滑低频段的预测越来越保守同时D_N会不断分裂出新的频段导致每个分量包含的信息量变少。当某个分量的方差小到接近数值精度时ARMA建模会退化成一个低方差常数预测意义不大。还有一种判断层数是否合适的思路分解到D_N后如果D_N的自相关函数呈现衰减较慢的拖尾形态说明它里面依然含有值得建模的周期性成分可再分解一层。如果D_N的ACF在零附近快速徘徊说明已经是噪声为主。我一般会写个循环逐层检查把这部分逻辑自动掉化不过在小项目中手动看一下也就够了。4.5 表格问题现象、原因与解决办法速查问题现象可能原因解决方法预测前段偏离大小波边界效应污染起点附近分量值对数据进行延拓或改用sym小波基整体预测偏保守几乎平稳近似分量阶数太低趋势外推能力不够适当提高A_N的AR阶数或考虑对A_N做二次差分建模高频分量预测值异常波动D1/D2拟合噪声成分限制低阶或直接预测为零训练集拟合很好测试集很差过拟合模型阶数太高降低p、q上限用AIC加约束用滚动验证方式检查稳定性forecast函数报错horizon类型不对或输入是行向量确保horizon是正标量历史数据统一为列向量分解重构后曲线两端发散边界效应叠加了高频分量发散用wdenoise对高低频分离或裁剪边界区域5. 实操过程中的进阶心得分解层数定完模型建完预测画出来之后真正的工程问题才刚开始。真实数据不像模拟数据这么友好尤其是低频趋势项的处理需要灵活调整。有些业务场景比如电力负荷预测近似分量A_3依然带有明显的非平稳特征。这时可以不用ARMA直接预测而是对A_3先做一阶差分差分平稳后用ARMA建模预测结果再逆差分回累计量。这一招特别管用时序里的趋势项在未来段的走向往往依赖最近几期的变化率而不是长期均值。另外ARMA模型对数据的量纲和规范化程度不敏感但小波分解不是。如果你把风速数据按m/s输入wavedec的尺度系数和高频系数的能量比例会显得失调。建议在分解前做一次标准化减去均值除标准差分解完成后预测结果再反标准化。这能显著提高高频分量建模的稳定性。关于工具箱版本我用过Matlab R2021b到R2023bwavedec和armax的接口没有本质变化。但如果你用的是较新版本要留意forecast函数对horizon参数类型的要求更严格了以前可以传double没问题新版有些情况要求传int32。具体问题具体看报错提示代码里可以写成int32(horizon)保险。最后再说一个很多人忽视的小细节小波分解的结果与边界延拓模式有关默认是周期延拓。如果你的数据端点值差异很大周期延拓会在边界制造伪影。可以试一下dwtmode(sym)切换到对称延拓模式实测对预测起点影响很大尤其是高频率分量。这是我看过Matlab官方文档里一句不起眼的话实际项目里真帮我解决过预测起点跳变的问题。用完之后记得dwtmode(per)切回来避免影响后续别的代码。这套方案我已经在不同的数据集上跑过很多遍总体感受是遇到非平稳、多尺度、有明显趋势和周期混叠的数据时它的优势非常明显值得作为首选方案之一但如果你的数据本身已经是平稳序列小波分解就是在画蛇添足。所以动手之前先用adftest做一次单位根检验再决定要不要套用小波分解这套流程完全来得及。
