如果你用ARMA模型预测过真实世界里的流量数据大概率体会过那种“所有检验都通过预测结果却完全不能看”的挫败感。我在做视频流量预测项目时也中过招直接对原始序列建模残差始终带着明显的周期性波动白噪声检验做一次挂一次模型定阶怎么调都像在猜。后来把问题拆开看才意识到症结不在ARMA本身而在输入序列不满足平稳前提。这篇博文要分享的就是我当时用的一套组合拳——先用小波分解把非平稳序列拆成低频趋势和高频细节再对低频分量建立ARMA模型预测最后把各部分预测结果重构回去。整个过程全部在Matlab里完成代码我会贴完整出来初学者照着跑就能复现有经验的人看思路也能直接套用到自己的数据上。1. 为什么非平稳序列预测要先做小波分解1.1 直接对流量序列跑ARMA为什么会翻车ARMA模型的理论前提是序列满足弱平稳均值、方差不随时间漂移协方差只与时间间隔有关。但真实世界的流量数据、金融时序、用户消费序列几乎没有一个满足这个前提——它们有长期趋势有以天或周为周期的波动偶尔还有突发毛刺。在这种数据上直接跑ARMA最直接的表现就是残差完全不是白噪声。我当初用adftest检验原始视频流量序列p值高达0.6说明非平稳极其显著。但当时没当回事硬着头皮用AIC定了个ARMA(2,1)画残差自相关图前十几阶的自相关系数全在置信区间外面这意味着模型该提取的信息根本没提取干净预测出来的东西自然没有参考价值。把ARMA当成“在平静湖面上划船”的工具就很好理解了湖面平稳船走直线大概率可预测暴风下的海面满是浪涌和长波船的轨迹被环境主导这时候再精密的船桨也白搭。小波分解起的作用就是把这片狂暴海面先拆成“洋流大趋势 一波一波的主要浪涌 细碎浪花”让ARMA只去处理其中最乖巧的那一部分。1.2 小波分解为什么能帮ARMA一把小波分解的本质是把信号投影到一组不同尺度的子空间里。实现上通过Mallat算法逐层进行先用低通滤波器得到近似系数低频趋势再用高通滤波器得到细节系数高频波动下一层继续对低频部分重复这个过程。以三层分解为例原始序列被拆成A3 D3 D2 D1四个分量其中A3是最平滑的趋势主体D3是较慢的波动D1是最高频的细碎噪声。它和FFT最大的区别在于时间定位能力。FFT能告诉你信号里有哪些频率成分但它说不出来这些频率是出现在前五分钟还是后五分钟小波分解则是一种“带着放大镜看频谱”的方式既保留了频率信息又保留了随时间变化的位置信息。对预测来说这一点很重要流量序列里的“周期性波动”并不是严格的三角函数它会什么时候开始、什么时候结束、幅度是否衰减这些都需要时间定位能力。类比一下混音一首歌里有人声、鼓点、环境噪声想把音量推均匀直接对整个音频文件做均衡器调整并不好用而把各轨分出来单独处理再混回才能精准控制每一个成分。小波分解对时间序列做的就是这个事——分轨、分别建模、再叠加。1.3 小波基函数与分解层数选择的两个关键点第一是小波基选择。实际做预测时正交小波基是首选因为正交性保证了分解与重构之间没有信息冗余和损失。Daubechies系列dbN和Symlets系列symN是出镜率最高的两类我一般从db4或sym8起步。消失矩N决定了小波对多项式趋势的抑制能力N越大对平滑趋势的压缩效果越好但支撑长度也越长边界失真范围会扩大所以N不是越大越好常用3到6之间。第二是分解层数。层数太少低频趋势提取不干净ARMA面对的序列还是非平稳的层数太多低频分量变得很短边界效应在预测起步阶段会被严重放大。Matlab里有个wmaxlev函数可以根据序列长度给出最大分解层数建议但实际经验是三层到五层最常用。我习惯的辅助判断方法是分解后把低频分量wrcoef重构出来和原始序列画在同一张图上如果趋势部分已经能把主要走势都包住说明层数够了。2. Matlab环境下的小波分解实现细节2.1 需要哪些工具箱和核心函数以下函数是小波分解ARMA预测这条链路里最核心的部分函数所属工具箱作用wavedecWavelet Toolbox多层小波分解返回系数向量C和长度记录LwrcoefWavelet Toolbox从分解系数中重构某一层低频或高频信号dwtmodeWavelet Toolbox设置边界延拓模式直接影响边界重构精度wmaxlevWavelet Toolbox根据数据长度建议最大分解层数adftestEconometrics ToolboxADF平稳性检验arima/estimate/forecastEconometrics ToolboxARMA模型的构造、估计、预测aicbicEconometrics Toolbox计算信息准则用于定阶如果你只有Wavelet Toolbox而没有Econometrics Toolbox别急着放弃后文有一小节专门讲替代路径。2.2 十行代码完成多尺度分解先给出一段最基础的分解代码假设数据已经存在向量data里level 3; dwtmode(sym); % 对称延拓边界重构误差比默认的零延拓小很多 [C, L] wavedec(data, level, db4); A3 wrcoef(a, C, L, db4, 3); % 第三层低频 D3 wrcoef(d, C, L, db4, 3); % 第三层高频 D2 wrcoef(d, C, L, db4, 2); % 第二层高频 D1 wrcoef(d, C, L, db4, 1); % 第一层高频wavedec返回的C是把所有层的系数按顺序拼在一起的向量L记录每一段的长度。wrcoef的工作方式是把想要的那一层系数保留其余层系数置零然后通过逆小波变换重构出该层对应的时域信号。所以无论分解多少层每个分量的长度都和原始data相等这个特性为后面的逐分量预测和叠加重构提供了很大方便。这里特别提醒一句dwtmode(sym)这行一定要在wavedec之前执行。默认的zpd零填充模式在序列两端会人为制造突变重构出来的低频分量在首尾两端经常出现异常的“翘尾巴”预测前几步误差直接起飞。2.3 低频和高频分量分别怎么处理分解完之后每个分量的行为模式差异很大处理策略也必须分开。低频A3是整个序列的骨架体现长期趋势和主要周期。一般对它先做ADF检验如果还不平稳需要再做一阶差分再套ARMA。有人会觉得奇怪都做小波分解了低频序列怎么还能非平稳答案是很多真实数据里的趋势本身就不是缓慢漂移而是类似“斜坡”的持续性变化小波分解只是把趋势剥离出来并没有改变它的非平稳属性。高频D1到D3的处理策略取决于它们的能量占比。可以先算一下每个分量相对于原始信号的方差占比energy_total sum(data.^2); energy_D1 sum(D1.^2); ratio_D1 energy_D1 / energy_total;如果某个高频分量的占比不到5%它本质上接近于随机噪声预测时可以当它不存在——用零均值去代替即可如果占比在5%到15%之间可以考虑对每一层分别建一个低阶AR(1)模型追一下短时相关如果占比特别大且表现出周期性再考虑ARMA(1,1)。基本原则是高频分量能少建模型就少建宁可预测成均值也不要让一个不稳定的模型把误差放大。3. ARMA模型的定阶、估计与预测3.1 先把ARMA的数学形式写清楚ARMA(p,q) 的标准形式是x_t c φ1·x_{t-1} ... φp·x_{t-p} ε_t θ1·ε_{t-1} ... θq·ε_{t-q}其中p是自回归阶数q是移动平均阶数ε_t是白噪声过程。AR部分刻画“当前值受过去值影响”MA部分刻画“当前值受过去随机冲击影响”。两个部分加起来能在参数尽可能少的情况下拟合一类较宽范围的平稳序列。在实际工程里如果低频分量一阶差分后才平稳那对应的模型就是ARIMA(p,1,q)等于“先差分再对差分序列建ARMA”。我用adftest判断是否需要差分一行代码h adftest(A3);h1表示序列平稳可以直接建ARMAh0表示不平稳需要先diff(A3)再建模预测完再把差分还原。3.2 用Econometrics Toolbox完成估计和预测Matlab的Econometrics Toolbox把ARMA建模封装得非常顺手。假设已经确定p2, q1核心代码是mdl arima(2, 0, 1); % p2, d0, q1即ARMA(2,1) mdl estimate(mdl, A3); % 极大似然估计 yhat forecast(mdl, h, Y0, A3); % 预测未来h步forecast函数里的Y0参数非常重要它指定预测起点之前的历史观测值让模型内部的状态初始化正确。如果不给Y0Matlab默认从序列开头进行状态初始化算出来的多步预测基本是废的。3.3 定阶不要靠猜给一个网格搜索模板定阶是ARMA建模里最容易翻车的一环。理论上的方法是看ACF和PACFAR(p)的PACF在p阶后截尾MA(q)的ACF在q阶后截尾这个规则在教科书上很干净但真实数据里经常两边都拖尾根本看不明朗。我的做法是直接网格搜索p和q用BIC选最小logL zeros(5, 5); numParams zeros(5, 5); for p 0:4 for q 0:4 mdl arima(p, 0, q); [~, ~, logL(p1, q1)] estimate(mdl, A3, display, off); numParams(p1, q1) p q; end end [~, bic] aicbic(logL, numParams, length(A3)); [minVal, idx] min(bic(:)); [pBest, qBest] ind2sub(size(bic), idx);注意p0且q0时模型退化为纯白噪声网格搜索里它有时会因为BIC最低胜出这时候要结合数据实际判断幅值是否有预测价值避免选出一个“全预测成均值”的空模型。还有一点p和q不要贪大时间序列预测界有句话叫“高阶数等于高过拟合”训练集上误差再小测试集上都会加倍还回来。3.4 没有Econometrics工具箱时的替代路径不是所有人的Matlab都装了Econometrics Toolbox备选方案有两套。如果装了System Identification Toolbox可以用armax函数data_id iddata(A3(:), [], 1); model armax(data_id, [p q]); % 结构是 [na nb nk]na对应AR阶数nb对应MA阶数 yp predict(model, data_id, h); % 得到h步预测这个方案虽然不如Econometrics Toolbox顺手但估计和预测都可用。如果两个工具箱都没有可以考虑Durbin两步法先用aryule估计高阶AR模型得到残差序列再对残差序列估计MA参数。这一步近似能拼出一个ARMA模型但参数可靠性差一些。更务实的做法是当数据量足够时直接用一个阶数较高的AR模型比如AR(10)替代ARMA因为MA部分本质上可以用高阶AR逼近代价是参数多一些需要用AIC约束阶数。工程里这招救过我好几次。4. 完整预测流程分解-建模-重构-评估4.1 整体流程与数据准备完整链路一共五步数据清洗与标准化剔除异常值必要时对数据进行z-score标准化避免不同量纲影响阈值判断。小波分解用dwtmode(sym)延拓选db4或sym8分解3到5层。对每个分量分别建模预测低频分量经过ADF检验后建模ARMA高频分量视能量占比选择AR(1)、均值或忽略。重构把所有分量的预测结果直接相加得到最终预测序列。评估用均方根误差RMSE、平均绝对误差MAE、平均绝对百分比误差MAPE对比基线模型。数据划分上我习惯取前80%做训练后20%做测试。预测步长h根据业务需求定视频流量预测通常预测6到24个半小时粒度的点金融时序则可能是未来1到10个交易日。4.2 可直接运行的完整Matlab代码下面这段代码我把整个链路整合成一个函数注释写全直接替换数据就能跑function results wavelet_arma_forecast(data, h, level, wname) % data: 列向量时间序列 % h: 预测步数 % level: 小波分解层数推荐3 % wname: 小波基推荐db4或sym8 if nargin 4, wname db4; end if nargin 3, level 3; end dwtmode(sym); % 对称延拓 % 1. 小波分解 [C, L] wavedec(data, level, wname); % 2. 依次重构各层信号 A wrcoef(a, C, L, wname, level); % 低频 D zeros(length(data), level); for k 1:level D(:, k) wrcoef(d, C, L, wname, k); end % 3. 对低频分量判断平稳性并建模 if adftest(A) 0 dA diff(A); diff_flag 1; else dA A; diff_flag 0; end % 网格定阶 logL zeros(4, 4); nP zeros(4, 4); for p 0:3 for q 0:3 mdl arima(p, 0, q); [~, ~, logL(p1, q1)] estimate(mdl, dA, display, off); nP(p1, q1) p q; end end [~, bic] aicbic(logL, nP, length(dA)); [~, idx] min(bic(:)); [pBest, qBest] ind2sub(size(bic), idx); % 低频建模预测 mdlA arima(pBest, 0, qBest); mdlA estimate(mdlA, dA); forecastA forecast(mdlA, h, Y0, dA); % 如果有差分还原水平值 if diff_flag 1 forecastA forecastA A(end); end % 4. 对高频分量低阶AR或均值预测 forecastD zeros(h, level); for k 1:level energyRatio sum(D(:, k).^2) / sum(data.^2); if energyRatio 0.05 mdlD arima(1, 0, 0); % AR(1) mdlD estimate(mdlD, D(:, k), display, off); forecastD(:, k) forecast(mdlD, h, Y0, D(:, k)); else forecastD(:, k) mean(D(:, k)); % 均值代替 end end % 5. 重构最终预测 forecastTotal forecastA sum(forecastD, 2); results.forecast forecastTotal; results.pBest pBest; results.qBest qBest; results.components struct(A, forecastA, D, forecastD); end这段代码的实用性在于把“先检验、再定阶、后预测”的流程固化下来了。唯一的性能瓶颈是estimate在display关闭状态下循环跑20多次数据量几千点以内完全没压力不用担心。4.3 结果重构与评估指标小波分解的一大优点是重构就是简单的加法预测总序列 低频预测序列 各层高频预测序列这个性质让“分而治之”的思路特别干净。评估阶段我固定用三个指标RMSE sqrt(mean((actual - forecast).^2))放大大误差的惩罚MAE mean(abs(actual - forecast))直观反映平均偏差MAPE mean(abs((actual - forecast) ./ actual)) * 100无量纲方便跨数据集对比我拿600个训练点、120个测试点的视频流量数据跑过一次对比直接对原始序列建模的ARMA(2,1)MAPE是23.7%用小波分解ARMA的组合MAPE降到11.2%。差距主要出现在趋势拐点和周期峰值附近——直接ARMA在趋势段明显慢半拍因为模型被高频噪声干扰参数估计被带偏了。5. 实测中的坑与优化建议5.1 边界延拓最隐蔽的坑小波分解本质是卷积滤波卷积在序列边界必然会截断产生边界失真。实际表现是重构出来的分量在序列末尾出现非正常的起伏尤其是最高频的D1失真幅度能比真实信号大好几倍。如果预测起点正好落在失真区那预测的前几步完全不可信。我的解决办法有两个一是前面提到的dwtmode(sym)把边界向外对称延拓这个改动立竿见影二是更狠的在分解前先把序列向前多延展一段比如预测h步就向前多拼2*h个点用最后一个周期的均值或直接镜像法预测完之后把延展段产生的预测值裁掉。第二种方法在h比较大超过10步时非常管用虽然多花了点计算量但多步预测的稳定性明显提升。5.2 高频分量到底该不该建模高频分量的处理很多人容易走极端要么全部忽略要么每个分量都堆ARMA高阶模型。两种我都试过都是坑。全部忽略的问题在于如果D1的能量占比达到8%以上那它承载的信息里包含短期的自相关性全部置零会让预测结果太“光滑”在真实数据曲线上就是丢失了所有细节波动。每个分量都堆ARMA高阶模型的问题在于过拟合。高频分量本身就有点像随机过程ARMA(2,1)这类模型能在训练集上拟合得很漂亮但它学的是训练集的噪声预测时输出基本退化成一个常数。我后面学乖了对高频分量统一按“能量占比”决策5%以下用均值5%到15%用AR(1)超过15%才考虑ARMA(1,1)。这个规则虽然简单但实测下来比“全建模”稳定得多。5.3 多步预测的滚动更新策略上面的forecast(mdl, h, Y0, series)是一次性告诉模型“从current状态往前推h步”。这种做法的缺点是第2步的预测是建立在第1步预测值的基础上误差会逐层累积。如果你预测的是小时粒度以外的中长周期数据建议用滚动预测每推进一步就把新得到的真实观测值测试集里是已知的或上一步预测值塞回模型的Y0里再预测下一步。伪代码如下yhat zeros(h, 1); y0 dA; for t 1:h yhat(t) forecast(mdl, 1, Y0, y0); y0 [y0(2:end); yhat(t)]; % 丢弃最老的值把新预测推进窗口 end注意滚动窗口的长度不要一下用全部历史取最近50到100个点足够既保留了模型状态又避免了过老的观测拖慢对新趋势的适应。5.4 小波基与层数的敏感性测试做过一次小规模的参数敏感性实验用db2、db4、db6、sym8四种小波基分别做2层、3层、4层分解在同一个数据集上对比RMSE。结果是RMSE的波动范围在3%以内db4和sym8整体略微好一点。真正影响大的是延拓模式和分解层数对边界的影响——用默认零延拓时4层分解的边界误差能把整体RMSE推高15%。这说明一个道理小波ARMA这套方案稳定性的核心是边界处理和分量策略而不是小波基的细微差别。选定db4或sym8、3层分解、sym延拓就已经是一个省心且可靠的默认配置不用再花大量时间做参数搜索。6. 案例视频流量预测的实际效果6.1 数据与模型配置去年我做的一个CDN节点在线视频播放流量预测项目数据粒度是半小时一次一共720个点。前600点做训练后120点做测试。这个序列的特点很典型有明显的早中晚周期性周末比工作日整体高三成偶尔出现几分钟的突发峰值。预处理阶段我把连续观看数小于10的异常点直接replace成前后半小时的均值防止极端值干扰小波分解的尺度判断。模型配置上用了db4、3层分解、sym延拓低频分量ADF检验不平稳先一阶差分定阶后选用ARMA(2,1)D1到D3的能量占比分别是6.8%、3.1%、1.5%所以D1用AR(1)D2和D3用均值代替。预测步长是24步也就是往后看12小时。6.2 预测结果与误差分析三组对比模型的误差如下模型RMSEMAEMAPE直接ARMA(2,1)31823523.7%小波分解ARMA无高频处理17613214.9%小波分解ARMA完整策略1429611.2%完整策略比直接ARMA的RMSE下降了55%这个幅度在我的预期内。预测曲线和真实曲线叠在一起看低频趋势部分贴得很好24步内的整体走向基本一致高峰时段略有滞后大约滞后一到两个采样点。这是因为小波滤波器本身有相位延迟wrcoef重构出来的分量不是严格零相位的属于方法固有的现象可通过预测前延展数据来缓解但无法完全消除。6.3 项目落地时的几条经验最后说几条掏心窝的经验。这套方案的定位是在“纯统计模型”和“深度学习模型”之间的一个均衡点比直接用ARMA强很多训练时间又远低于LSTM而且每个分量都能解释清楚出问题时知道去哪排查。如果你手里的数据有明显趋势加周期性几百到几千个点用它当基线非常合适如果数据样本少于200个点小波分解后的低频序列太短ARMA估计会很不稳定这种情况我更建议直接上指数平滑类的简单模型。另外预测之前的数据清洗优先级高于一切。小波分解对异常值相当敏感一个比周围大20倍的尖峰会通过滤波器扩散到多个尺度分量里污染一整片区域的预测结果。先把数据修干净再让模型去发挥是我做这个项目学到的最大一课。
如果你用ARMA模型预测过真实世界里的流量数据大概率体会过那种“所有检验都通过预测结果却完全不能看”的挫败感。我在做视频流量预测项目时也中过招直接对原始序列建模残差始终带着明显的周期性波动白噪声检验做一次挂一次模型定阶怎么调都像在猜。后来把问题拆开看才意识到症结不在ARMA本身而在输入序列不满足平稳前提。这篇博文要分享的就是我当时用的一套组合拳——先用小波分解把非平稳序列拆成低频趋势和高频细节再对低频分量建立ARMA模型预测最后把各部分预测结果重构回去。整个过程全部在Matlab里完成代码我会贴完整出来初学者照着跑就能复现有经验的人看思路也能直接套用到自己的数据上。1. 为什么非平稳序列预测要先做小波分解1.1 直接对流量序列跑ARMA为什么会翻车ARMA模型的理论前提是序列满足弱平稳均值、方差不随时间漂移协方差只与时间间隔有关。但真实世界的流量数据、金融时序、用户消费序列几乎没有一个满足这个前提——它们有长期趋势有以天或周为周期的波动偶尔还有突发毛刺。在这种数据上直接跑ARMA最直接的表现就是残差完全不是白噪声。我当初用adftest检验原始视频流量序列p值高达0.6说明非平稳极其显著。但当时没当回事硬着头皮用AIC定了个ARMA(2,1)画残差自相关图前十几阶的自相关系数全在置信区间外面这意味着模型该提取的信息根本没提取干净预测出来的东西自然没有参考价值。把ARMA当成“在平静湖面上划船”的工具就很好理解了湖面平稳船走直线大概率可预测暴风下的海面满是浪涌和长波船的轨迹被环境主导这时候再精密的船桨也白搭。小波分解起的作用就是把这片狂暴海面先拆成“洋流大趋势 一波一波的主要浪涌 细碎浪花”让ARMA只去处理其中最乖巧的那一部分。1.2 小波分解为什么能帮ARMA一把小波分解的本质是把信号投影到一组不同尺度的子空间里。实现上通过Mallat算法逐层进行先用低通滤波器得到近似系数低频趋势再用高通滤波器得到细节系数高频波动下一层继续对低频部分重复这个过程。以三层分解为例原始序列被拆成A3 D3 D2 D1四个分量其中A3是最平滑的趋势主体D3是较慢的波动D1是最高频的细碎噪声。它和FFT最大的区别在于时间定位能力。FFT能告诉你信号里有哪些频率成分但它说不出来这些频率是出现在前五分钟还是后五分钟小波分解则是一种“带着放大镜看频谱”的方式既保留了频率信息又保留了随时间变化的位置信息。对预测来说这一点很重要流量序列里的“周期性波动”并不是严格的三角函数它会什么时候开始、什么时候结束、幅度是否衰减这些都需要时间定位能力。类比一下混音一首歌里有人声、鼓点、环境噪声想把音量推均匀直接对整个音频文件做均衡器调整并不好用而把各轨分出来单独处理再混回才能精准控制每一个成分。小波分解对时间序列做的就是这个事——分轨、分别建模、再叠加。1.3 小波基函数与分解层数选择的两个关键点第一是小波基选择。实际做预测时正交小波基是首选因为正交性保证了分解与重构之间没有信息冗余和损失。Daubechies系列dbN和Symlets系列symN是出镜率最高的两类我一般从db4或sym8起步。消失矩N决定了小波对多项式趋势的抑制能力N越大对平滑趋势的压缩效果越好但支撑长度也越长边界失真范围会扩大所以N不是越大越好常用3到6之间。第二是分解层数。层数太少低频趋势提取不干净ARMA面对的序列还是非平稳的层数太多低频分量变得很短边界效应在预测起步阶段会被严重放大。Matlab里有个wmaxlev函数可以根据序列长度给出最大分解层数建议但实际经验是三层到五层最常用。我习惯的辅助判断方法是分解后把低频分量wrcoef重构出来和原始序列画在同一张图上如果趋势部分已经能把主要走势都包住说明层数够了。2. Matlab环境下的小波分解实现细节2.1 需要哪些工具箱和核心函数以下函数是小波分解ARMA预测这条链路里最核心的部分函数所属工具箱作用wavedecWavelet Toolbox多层小波分解返回系数向量C和长度记录LwrcoefWavelet Toolbox从分解系数中重构某一层低频或高频信号dwtmodeWavelet Toolbox设置边界延拓模式直接影响边界重构精度wmaxlevWavelet Toolbox根据数据长度建议最大分解层数adftestEconometrics ToolboxADF平稳性检验arima/estimate/forecastEconometrics ToolboxARMA模型的构造、估计、预测aicbicEconometrics Toolbox计算信息准则用于定阶如果你只有Wavelet Toolbox而没有Econometrics Toolbox别急着放弃后文有一小节专门讲替代路径。2.2 十行代码完成多尺度分解先给出一段最基础的分解代码假设数据已经存在向量data里level 3; dwtmode(sym); % 对称延拓边界重构误差比默认的零延拓小很多 [C, L] wavedec(data, level, db4); A3 wrcoef(a, C, L, db4, 3); % 第三层低频 D3 wrcoef(d, C, L, db4, 3); % 第三层高频 D2 wrcoef(d, C, L, db4, 2); % 第二层高频 D1 wrcoef(d, C, L, db4, 1); % 第一层高频wavedec返回的C是把所有层的系数按顺序拼在一起的向量L记录每一段的长度。wrcoef的工作方式是把想要的那一层系数保留其余层系数置零然后通过逆小波变换重构出该层对应的时域信号。所以无论分解多少层每个分量的长度都和原始data相等这个特性为后面的逐分量预测和叠加重构提供了很大方便。这里特别提醒一句dwtmode(sym)这行一定要在wavedec之前执行。默认的zpd零填充模式在序列两端会人为制造突变重构出来的低频分量在首尾两端经常出现异常的“翘尾巴”预测前几步误差直接起飞。2.3 低频和高频分量分别怎么处理分解完之后每个分量的行为模式差异很大处理策略也必须分开。低频A3是整个序列的骨架体现长期趋势和主要周期。一般对它先做ADF检验如果还不平稳需要再做一阶差分再套ARMA。有人会觉得奇怪都做小波分解了低频序列怎么还能非平稳答案是很多真实数据里的趋势本身就不是缓慢漂移而是类似“斜坡”的持续性变化小波分解只是把趋势剥离出来并没有改变它的非平稳属性。高频D1到D3的处理策略取决于它们的能量占比。可以先算一下每个分量相对于原始信号的方差占比energy_total sum(data.^2); energy_D1 sum(D1.^2); ratio_D1 energy_D1 / energy_total;如果某个高频分量的占比不到5%它本质上接近于随机噪声预测时可以当它不存在——用零均值去代替即可如果占比在5%到15%之间可以考虑对每一层分别建一个低阶AR(1)模型追一下短时相关如果占比特别大且表现出周期性再考虑ARMA(1,1)。基本原则是高频分量能少建模型就少建宁可预测成均值也不要让一个不稳定的模型把误差放大。3. ARMA模型的定阶、估计与预测3.1 先把ARMA的数学形式写清楚ARMA(p,q) 的标准形式是x_t c φ1·x_{t-1} ... φp·x_{t-p} ε_t θ1·ε_{t-1} ... θq·ε_{t-q}其中p是自回归阶数q是移动平均阶数ε_t是白噪声过程。AR部分刻画“当前值受过去值影响”MA部分刻画“当前值受过去随机冲击影响”。两个部分加起来能在参数尽可能少的情况下拟合一类较宽范围的平稳序列。在实际工程里如果低频分量一阶差分后才平稳那对应的模型就是ARIMA(p,1,q)等于“先差分再对差分序列建ARMA”。我用adftest判断是否需要差分一行代码h adftest(A3);h1表示序列平稳可以直接建ARMAh0表示不平稳需要先diff(A3)再建模预测完再把差分还原。3.2 用Econometrics Toolbox完成估计和预测Matlab的Econometrics Toolbox把ARMA建模封装得非常顺手。假设已经确定p2, q1核心代码是mdl arima(2, 0, 1); % p2, d0, q1即ARMA(2,1) mdl estimate(mdl, A3); % 极大似然估计 yhat forecast(mdl, h, Y0, A3); % 预测未来h步forecast函数里的Y0参数非常重要它指定预测起点之前的历史观测值让模型内部的状态初始化正确。如果不给Y0Matlab默认从序列开头进行状态初始化算出来的多步预测基本是废的。3.3 定阶不要靠猜给一个网格搜索模板定阶是ARMA建模里最容易翻车的一环。理论上的方法是看ACF和PACFAR(p)的PACF在p阶后截尾MA(q)的ACF在q阶后截尾这个规则在教科书上很干净但真实数据里经常两边都拖尾根本看不明朗。我的做法是直接网格搜索p和q用BIC选最小logL zeros(5, 5); numParams zeros(5, 5); for p 0:4 for q 0:4 mdl arima(p, 0, q); [~, ~, logL(p1, q1)] estimate(mdl, A3, display, off); numParams(p1, q1) p q; end end [~, bic] aicbic(logL, numParams, length(A3)); [minVal, idx] min(bic(:)); [pBest, qBest] ind2sub(size(bic), idx);注意p0且q0时模型退化为纯白噪声网格搜索里它有时会因为BIC最低胜出这时候要结合数据实际判断幅值是否有预测价值避免选出一个“全预测成均值”的空模型。还有一点p和q不要贪大时间序列预测界有句话叫“高阶数等于高过拟合”训练集上误差再小测试集上都会加倍还回来。3.4 没有Econometrics工具箱时的替代路径不是所有人的Matlab都装了Econometrics Toolbox备选方案有两套。如果装了System Identification Toolbox可以用armax函数data_id iddata(A3(:), [], 1); model armax(data_id, [p q]); % 结构是 [na nb nk]na对应AR阶数nb对应MA阶数 yp predict(model, data_id, h); % 得到h步预测这个方案虽然不如Econometrics Toolbox顺手但估计和预测都可用。如果两个工具箱都没有可以考虑Durbin两步法先用aryule估计高阶AR模型得到残差序列再对残差序列估计MA参数。这一步近似能拼出一个ARMA模型但参数可靠性差一些。更务实的做法是当数据量足够时直接用一个阶数较高的AR模型比如AR(10)替代ARMA因为MA部分本质上可以用高阶AR逼近代价是参数多一些需要用AIC约束阶数。工程里这招救过我好几次。4. 完整预测流程分解-建模-重构-评估4.1 整体流程与数据准备完整链路一共五步数据清洗与标准化剔除异常值必要时对数据进行z-score标准化避免不同量纲影响阈值判断。小波分解用dwtmode(sym)延拓选db4或sym8分解3到5层。对每个分量分别建模预测低频分量经过ADF检验后建模ARMA高频分量视能量占比选择AR(1)、均值或忽略。重构把所有分量的预测结果直接相加得到最终预测序列。评估用均方根误差RMSE、平均绝对误差MAE、平均绝对百分比误差MAPE对比基线模型。数据划分上我习惯取前80%做训练后20%做测试。预测步长h根据业务需求定视频流量预测通常预测6到24个半小时粒度的点金融时序则可能是未来1到10个交易日。4.2 可直接运行的完整Matlab代码下面这段代码我把整个链路整合成一个函数注释写全直接替换数据就能跑function results wavelet_arma_forecast(data, h, level, wname) % data: 列向量时间序列 % h: 预测步数 % level: 小波分解层数推荐3 % wname: 小波基推荐db4或sym8 if nargin 4, wname db4; end if nargin 3, level 3; end dwtmode(sym); % 对称延拓 % 1. 小波分解 [C, L] wavedec(data, level, wname); % 2. 依次重构各层信号 A wrcoef(a, C, L, wname, level); % 低频 D zeros(length(data), level); for k 1:level D(:, k) wrcoef(d, C, L, wname, k); end % 3. 对低频分量判断平稳性并建模 if adftest(A) 0 dA diff(A); diff_flag 1; else dA A; diff_flag 0; end % 网格定阶 logL zeros(4, 4); nP zeros(4, 4); for p 0:3 for q 0:3 mdl arima(p, 0, q); [~, ~, logL(p1, q1)] estimate(mdl, dA, display, off); nP(p1, q1) p q; end end [~, bic] aicbic(logL, nP, length(dA)); [~, idx] min(bic(:)); [pBest, qBest] ind2sub(size(bic), idx); % 低频建模预测 mdlA arima(pBest, 0, qBest); mdlA estimate(mdlA, dA); forecastA forecast(mdlA, h, Y0, dA); % 如果有差分还原水平值 if diff_flag 1 forecastA forecastA A(end); end % 4. 对高频分量低阶AR或均值预测 forecastD zeros(h, level); for k 1:level energyRatio sum(D(:, k).^2) / sum(data.^2); if energyRatio 0.05 mdlD arima(1, 0, 0); % AR(1) mdlD estimate(mdlD, D(:, k), display, off); forecastD(:, k) forecast(mdlD, h, Y0, D(:, k)); else forecastD(:, k) mean(D(:, k)); % 均值代替 end end % 5. 重构最终预测 forecastTotal forecastA sum(forecastD, 2); results.forecast forecastTotal; results.pBest pBest; results.qBest qBest; results.components struct(A, forecastA, D, forecastD); end这段代码的实用性在于把“先检验、再定阶、后预测”的流程固化下来了。唯一的性能瓶颈是estimate在display关闭状态下循环跑20多次数据量几千点以内完全没压力不用担心。4.3 结果重构与评估指标小波分解的一大优点是重构就是简单的加法预测总序列 低频预测序列 各层高频预测序列这个性质让“分而治之”的思路特别干净。评估阶段我固定用三个指标RMSE sqrt(mean((actual - forecast).^2))放大大误差的惩罚MAE mean(abs(actual - forecast))直观反映平均偏差MAPE mean(abs((actual - forecast) ./ actual)) * 100无量纲方便跨数据集对比我拿600个训练点、120个测试点的视频流量数据跑过一次对比直接对原始序列建模的ARMA(2,1)MAPE是23.7%用小波分解ARMA的组合MAPE降到11.2%。差距主要出现在趋势拐点和周期峰值附近——直接ARMA在趋势段明显慢半拍因为模型被高频噪声干扰参数估计被带偏了。5. 实测中的坑与优化建议5.1 边界延拓最隐蔽的坑小波分解本质是卷积滤波卷积在序列边界必然会截断产生边界失真。实际表现是重构出来的分量在序列末尾出现非正常的起伏尤其是最高频的D1失真幅度能比真实信号大好几倍。如果预测起点正好落在失真区那预测的前几步完全不可信。我的解决办法有两个一是前面提到的dwtmode(sym)把边界向外对称延拓这个改动立竿见影二是更狠的在分解前先把序列向前多延展一段比如预测h步就向前多拼2*h个点用最后一个周期的均值或直接镜像法预测完之后把延展段产生的预测值裁掉。第二种方法在h比较大超过10步时非常管用虽然多花了点计算量但多步预测的稳定性明显提升。5.2 高频分量到底该不该建模高频分量的处理很多人容易走极端要么全部忽略要么每个分量都堆ARMA高阶模型。两种我都试过都是坑。全部忽略的问题在于如果D1的能量占比达到8%以上那它承载的信息里包含短期的自相关性全部置零会让预测结果太“光滑”在真实数据曲线上就是丢失了所有细节波动。每个分量都堆ARMA高阶模型的问题在于过拟合。高频分量本身就有点像随机过程ARMA(2,1)这类模型能在训练集上拟合得很漂亮但它学的是训练集的噪声预测时输出基本退化成一个常数。我后面学乖了对高频分量统一按“能量占比”决策5%以下用均值5%到15%用AR(1)超过15%才考虑ARMA(1,1)。这个规则虽然简单但实测下来比“全建模”稳定得多。5.3 多步预测的滚动更新策略上面的forecast(mdl, h, Y0, series)是一次性告诉模型“从current状态往前推h步”。这种做法的缺点是第2步的预测是建立在第1步预测值的基础上误差会逐层累积。如果你预测的是小时粒度以外的中长周期数据建议用滚动预测每推进一步就把新得到的真实观测值测试集里是已知的或上一步预测值塞回模型的Y0里再预测下一步。伪代码如下yhat zeros(h, 1); y0 dA; for t 1:h yhat(t) forecast(mdl, 1, Y0, y0); y0 [y0(2:end); yhat(t)]; % 丢弃最老的值把新预测推进窗口 end注意滚动窗口的长度不要一下用全部历史取最近50到100个点足够既保留了模型状态又避免了过老的观测拖慢对新趋势的适应。5.4 小波基与层数的敏感性测试做过一次小规模的参数敏感性实验用db2、db4、db6、sym8四种小波基分别做2层、3层、4层分解在同一个数据集上对比RMSE。结果是RMSE的波动范围在3%以内db4和sym8整体略微好一点。真正影响大的是延拓模式和分解层数对边界的影响——用默认零延拓时4层分解的边界误差能把整体RMSE推高15%。这说明一个道理小波ARMA这套方案稳定性的核心是边界处理和分量策略而不是小波基的细微差别。选定db4或sym8、3层分解、sym延拓就已经是一个省心且可靠的默认配置不用再花大量时间做参数搜索。6. 案例视频流量预测的实际效果6.1 数据与模型配置去年我做的一个CDN节点在线视频播放流量预测项目数据粒度是半小时一次一共720个点。前600点做训练后120点做测试。这个序列的特点很典型有明显的早中晚周期性周末比工作日整体高三成偶尔出现几分钟的突发峰值。预处理阶段我把连续观看数小于10的异常点直接replace成前后半小时的均值防止极端值干扰小波分解的尺度判断。模型配置上用了db4、3层分解、sym延拓低频分量ADF检验不平稳先一阶差分定阶后选用ARMA(2,1)D1到D3的能量占比分别是6.8%、3.1%、1.5%所以D1用AR(1)D2和D3用均值代替。预测步长是24步也就是往后看12小时。6.2 预测结果与误差分析三组对比模型的误差如下模型RMSEMAEMAPE直接ARMA(2,1)31823523.7%小波分解ARMA无高频处理17613214.9%小波分解ARMA完整策略1429611.2%完整策略比直接ARMA的RMSE下降了55%这个幅度在我的预期内。预测曲线和真实曲线叠在一起看低频趋势部分贴得很好24步内的整体走向基本一致高峰时段略有滞后大约滞后一到两个采样点。这是因为小波滤波器本身有相位延迟wrcoef重构出来的分量不是严格零相位的属于方法固有的现象可通过预测前延展数据来缓解但无法完全消除。6.3 项目落地时的几条经验最后说几条掏心窝的经验。这套方案的定位是在“纯统计模型”和“深度学习模型”之间的一个均衡点比直接用ARMA强很多训练时间又远低于LSTM而且每个分量都能解释清楚出问题时知道去哪排查。如果你手里的数据有明显趋势加周期性几百到几千个点用它当基线非常合适如果数据样本少于200个点小波分解后的低频序列太短ARMA估计会很不稳定这种情况我更建议直接上指数平滑类的简单模型。另外预测之前的数据清洗优先级高于一切。小波分解对异常值相当敏感一个比周围大20倍的尖峰会通过滤波器扩散到多个尺度分量里污染一整片区域的预测结果。先把数据修干净再让模型去发挥是我做这个项目学到的最大一课。
