简介一篇面向时间序列预测、设备状态监测及组合模型研究人员的学术论文主题是基于加权马尔可夫链修正的ARIMA预测模型。该研究针对ARIMA模型在设备状态参数预测中的偏差与不稳定性引入加权马尔可夫链对残差序列建模结合状态特征值与线性插值把预测残差状态转化为具体值再以残差修正值作为最终预测结果并以船舶海水出口温度预测为例验证了修正模型较单一ARIMA模型预测精度显著提升。资源仅1个PDF文件约767KB内含论文全文、摘要、图表、公式与参考文献能直接用于算法、机器学习、人工智能方向的技术研读。已有514人学习下载既适合组合预测与设备状态分析的研究者借鉴思路也可为视情维修、状态监控等工程应用提供参考。1. 用ARIMA做GDP季度预测最怕的不是不收敛而是对“拐弯”没反应我用ARIMA(1,1,1)拟合一组80期的季度数据训练集MAPE能压到1.5%可一到预测段增速切换出现连续两期偏移误差直接跳到6%以上。原因很直接ARIMA用差分和自回归把趋势线性化状态跳跃这类信息它天生不敏感。后来拆到一份PDF《基于加权马尔可夫链修正的ARIMA预测模型的研究》思路说穿了就一句话ARIMA先兜底预测再把残差划分成状态用加权马尔可夫链修正系统性偏差。这套组合适合GDP、销量、客流这类有趋势又有波动的中短期预测也适合ARIMA精度卡住又不想换XGBoost整套框架的人。这篇笔记把状态划分、转移矩阵、权重设计和落地代码过一次最后记录我复现时踩过的坑。2. ARIMA与加权马尔可夫链为什么偏要两个模型叠一起2.1 ARIMA的强项与边界趋势之外的信息去哪了ARIMA(p,d,q)的含义不复杂p是自回归阶数d是差分数q是移动平均阶数。它拟合的是线性相关的结构输出是过去观测值和过去误差的线性组合。对GDP这类有明确趋势的序列d取1或2之后差分序列大体平稳ARIMA表现不错。但线性组合决定了它对状态切换不敏感政策冲击、节日效应、经营节奏切换这类信号会成块地留在残差里。换句话说ARIMA的残差不等于白噪声它往往带着状态信息只是模型没能力把它吃进去。我习惯在拟合完ARIMA之后立刻跑一次Ljung-Box检验看残差是否还有自相关。如果p值小于0.05说明残差里还有东西可挖这时候直接上马尔可夫链修正是有意义的如果残差已经是白噪声强行加修正反而会把噪声当成信号后面第6章会专门讲这个验证步骤。2.2 加权马尔可夫链修正状态划分、转移概率和权重普通马尔可夫链的基本假设是下一个状态只依赖当前状态用一步转移概率矩阵描述。对ARIMA残差来说这个假设有点浪费残差状态不仅和上一期有关还可能和更早的状态存在周期性关联。加权马尔可夫链做的事情是把滞后1阶、2阶、3阶的状态转移概率都拿出来分别查对应的转移矩阵再按权重线性相加得到当前残差的状态分布。状态划分是准备工作。把ARIMA的训练残差按分位数切成3到5个状态常见命名是负大偏差、负小偏差、正常、正小偏差、正大偏差。每个状态对应一个状态中心通常是该状态下训练残差的均值。修正的时候先算当前时刻处于各状态的概率再对状态中心做加权平均得到期望残差加到ARIMA的基准预测上。权重在这里是核心。滞后1阶的转移样本最多、短期依赖最强滞后3阶的样本减少但能反映季节周期或政策滞后效应。固定权重0.5、0.3、0.2是稳妥的起步值但要拿更好的精度就得根据数据特性重新标定第4章会给具体的定权方法。2.3 完整的五步建模路线从ARIMA到组合修正整套流程我拆成五个阶段每阶段都有明确的输入输出。第一个阶段是数据切分时序数据必须按时间顺序切不能随机打散否则状态转移矩阵会学到未来信息。第二个阶段是ARIMA建模在训练集上用AIC或BIC扫描定阶同时检查残差平稳性。第三个阶段是提取残差训练残差用来划分状态和统计状态中心预测残差是修正的对象。第四个阶段是状态划分与转移矩阵构建用分位数定阈值再用平滑拉普拉斯避免零概率。第五个阶段是加权修正把多阶转移概率加权求和得到修正期望并叠加到基准预测上。阶段输入输出关键参数数据切分时间序列train / val / test窗口长度、验证段比例ARIMA建模train序列拟合模型、残差p、d、qAIC准则状态划分ARIMA训练残差状态序列、分位数阈值状态数、分位数位置转移矩阵状态序列1到K阶转移矩阵最大滞后阶K、平滑系数加权修正基预测、转移矩阵最终预测值各阶权重、修正阈值这个五步流程在论文里是理论基础落到代码上就是第3章的内容。3. 把论文模型落地成可复现代码状态划分、转移矩阵与权重计算3.1 数据准备与ARIMA基线用小样本季度数据演示我用一份模拟的季度数据来做演示80期带线性趋势和周期性波动量级接近GDP数据。选线性趋势是为了让一阶差分后序列平稳避免模拟数据本身不满足ARIMA假设那会干扰后面马尔可夫修正的效果判断。import numpy as np import pandas as pd from statsmodels.tsa.arima.model import ARIMA from statsmodels.tsa.stattools import adfuller np.random.seed(42) t np.arange(80) trend 100 0.15 * t season 4 * np.sin(2 * np.pi * t / 4) noise np.random.normal(0, 1.2, sizet.shape) y trend season noise df pd.DataFrame({t: t, y: y}) train df.iloc[:70] test df.iloc[70:] adf_res adfuller(np.diff(train[y].values), autolagAIC) print(一阶差分后ADF p值:, adf_res[1]) grid {} for p in range(0, 4): for q in range(0, 4): try: m ARIMA(train[y], order(p, 1, q)).fit() grid[(p, q)] m.aic except Exception: continue best_order min(grid, keygrid.get) model ARIMA(train[y], orderbest_order).fit() print(最优阶数:, best_order, AIC:, grid[best_order])这段代码把ADF检验放在训练集内做只判断一阶差分后是否平稳正常p值应小于0.05如果大于这个值说明d1不够需要把差分阶数调到2。AIC扫描同样是只在训练集上跑避免测试集信息泄漏到模型选择里。这里扫描了p和q在0到3之间的全部组合样本量80期完全够用。3.2 残差状态划分与K阶转移矩阵ARIMA拟合之后取训练集内从第1期开始的拟合值计算残差序列。状态划分我优先用分位数而不是等宽阈值因为残差分布不一定对称分位数能让每个状态都有足够的样本后续转移矩阵不容易出现整行零。fit_line model.predict(start1, endlen(train) - 1, typlevels) resid train[y].iloc[1:].values - fit_line.values quants np.quantile(resid, [0.15, 0.35, 0.65, 0.85]) states np.digitize(resid, quants) def build_transition_matrix(states_seq, max_lag3, n_states5, alpha0.1): mats [] for lag in range(1, max_lag 1): cnt np.full((n_states, n_states), alpha) for i in range(len(states_seq) - lag): cnt[states_seq[i], states_seq[i lag]] 1 mats.append(cnt / cnt.sum(axis1, keepdimsTrue)) return mats trans_mats build_transition_matrix(states) print(1阶转移矩阵:) print(np.round(trans_mats[0], 3))状态划分的分位数取了15%、35%、65%、85%四个位置正好把残差切成5个区间。np.digitize返回0到4的编号对应负大偏差到正大偏差。build_transition_matrix里alpha是拉普拉斯平滑系数给每个转移对加了0.1的伪计数这样即使某些状态对在历史里一次都没出现过矩阵也不会出现零概率行后面的加权求和不会产生除零错误。3.3 加权马氏修正逐期输出最终预测修正的核心是每个测试点都用历史状态序列预测当前残差状态分布。注意测试期真实值不会立刻拿到滚动修正时需要等真实值到达再回填状态如果下一期真实值还没公布一般先用ARIMA基准预测值顶替残差计算状态等真实值出来再修正历史状态序列。state_centers np.array([ resid[states s].mean() for s in range(5) ]) weights np.array([0.5, 0.3, 0.2]) base_forecast model.predict(startlen(train), endlen(df) - 1, typlevels).values true_y test[y].values history list(states) corrected [] for i in range(len(true_y)): p np.zeros(5) for k, w in enumerate(weights): lag k 1 if len(history) lag: p w * trans_mats[k][history[-lag], :] p p / p.sum() exp_resid np.dot(p, state_centers) corrected.append(base_forecast[i] exp_resid) real_resid true_y[i] - base_forecast[i] real_state int(np.digitize(real_resid, quants)) history.append(real_state) def mape(y_true, y_pred): return np.mean(np.abs((y_true - y_pred) / y_true)) * 100 print(ARIMA MAPE:, round(mape(true_y, base_forecast), 3)) print(ARIMAMarkov MAPE:, round(mape(true_y, np.array(corrected)), 3))这里weights是滞后1阶、2阶、3阶的权重0.5、0.3、0.2是先验值含义是短期状态依赖比长期更强。循环里每次用history最后1个、2个、3个真实状态分别查对应阶数的转移矩阵得到三个概率向量再加权求和。归一化后与状态中心做点积得到期望残差修正。需要特别注意state_centers完全来自训练残差测试期真实状态只用来更新history不会重新计算阈值和状态中心这是防止数据泄漏的关键边界。3.4 关键参数速查表参数常见取值说明p、d、qAIC扫描确定d1对季度趋势多数够用残差非平稳就升d状态数3到5个太少没区分度太多转移矩阵稀疏分位数位置15%, 35%, 65%, 85%让各状态样本量大致均衡最大滞后阶K2到4太大样本减少且权重难定平滑系数alpha0.1到0.5防止零转移概率对短序列尤其重要权重组合0.5, 0.3, 0.2起步再用第4章方法按数据重估4. 权重到底怎么定从频率到贴近度的多种定权方式4.1 为什么权重不能随便拍滞后阶的信息量并不相等滞后1阶的转移直接反映相邻两期残差状态的关系样本量最大统计上最可靠。滞后2阶、3阶能捕捉到周期性和政策滞后效应但构造转移矩阵时每增加一阶可用样本就减少一个矩阵估计的方差变大。如果拍脑袋给它们相同权重等于把噪声和信号同等对待。这就是为什么论文里的加权设计值得单独研究。固定权重0.5、0.3、0.2适合快速验证流程通不通但要拿更好的预测精度得在验证段上比较各个滞后阶的独立预测能力再确定权重。4.2 三种可落地的定权方法频数、区分度、贴近度频数权重最简单直接认为样本量越大的滞后阶越可靠。区分度权重衡量转移矩阵各行之间的差异程度如果某阶矩阵从不同状态出发的概率分布差别很大说明这阶矩阵能有效区分状态路径权重应该更高。贴近度权重则是在验证段上实际跑一遍每种滞后阶的单独修正效果MAPE小的给高权重。from scipy.spatial.distance import pdist def frequency_weight(max_lag, states_len): counts np.array([states_len - lag for lag in range(1, max_lag 1)]) return counts / counts.sum() def discrimination_weight(trans_mats): scores [] for m in trans_mats: scores.append(pdist(m, metriceuclidean).mean()) scores np.array(scores) 1e-8 return scores / scores.sum() def closeness_weight(trans_mats, centers, val_y, val_base, train_states, quants): hist np.array( list(train_states) [np.digitize(v, quants) for v in (np.array(val_y) - np.array(val_base))] ) errors [] for k in range(1, len(trans_mats) 1): m trans_mats[k - 1] pred [] for t in range(len(val_y)): idx len(train_states) t - k if idx 0: pred.append(0.0) continue pred.append(np.dot(m[int(hist[idx]), :], centers)) adj np.array(val_base) np.array(pred) errors.append(mape(val_y, adj)) inv 1 / (np.array(errors) 1e-6) return inv / inv.sum()frequency_weight计算的是每个滞后阶实际可用样本对总数只看数据和滞后阶长度。discrimination_weight用pdist计算转移矩阵各行间的欧氏距离均值距离越大说明该阶矩阵的状态区分能力越强。closeness_weight把每个滞后阶单独当成完整修正器跑验证段统计MAPE后取倒数归一化本质上是在衡量“这个阶数是不是真的能预测”。closeness_weight是三种里面最贴近实际预测场景的但要注意它用的是真实残差状态序列来计算误差相当于近似评估完整做法应该滚动更新状态。对小样本序列如果这三种方法给出的权重差异很大优先相信频数权重因为高阶矩阵在小样本下方差太大区分度和贴近度都会被噪声主导。4.3 定权后的交叉验证权重不是越复杂越好我拿训练集最后15期当验证段分别用三种权重做滚动修正对比权重差异和MAPE变化。经验是如果三种权重下的MAPE差异小于5%直接选频数权重省事且稳定如果某个滞后阶在区分度或贴近度下权重超过0.6要回头检查状态划分是否过细而不是急着调权重。定权方式适用场景注意事项频数权重小样本、快速起步完全忽略预测能力区分度权重状态依赖明显、矩阵稳定稀疏矩阵下会高估噪声阶贴近度权重验证段充足、追求精度计算量最大注意状态泄漏5. 避坑与常见问题ARIMA马尔可夫链的四个翻车现场5.1 状态个数拍脑袋设成8个转移矩阵稀疏成什么样现象80期样本状态设8个每个状态平均不到10个样本二阶以上的转移矩阵大量位置是零加权修正结果基本等于ARIMA原预测加了等于没加。原因状态划分过细之后每个状态出现的频次不够转移矩阵变成了稀疏矩阵马尔可夫链根本学不到稳定的转移规律。解决状态数控制在3到5个用分位数划分保证样本均衡并在计数矩阵上做拉普拉斯平滑。def laplace_smooth(count_mat, alpha0.1): smoothed count_mat alpha return smoothed / smoothed.sum(axis1, keepdimsTrue)5.2 残差非平稳直接套马尔可夫修正项变成趋势跟随现象ARIMA用d0直接拟合带趋势的GDP残差序列ADF检验p值到0.3状态序列几乎只在负大偏差和正大偏差两个状态之间来回跳修正结果变成一条跟随原序列趋势的斜线。原因残差本身非平稳状态序列不具备马尔可夫性转移矩阵反映的是趋势漂移而不是真实状态切换。解决拟合后立刻对残差跑adfuller非平稳就升d或做对数变换另一个实用技巧是改用相对误差作为状态输入因为GDP量级逐年增大早期绝对残差小、后期绝对残差大用固定阈值会失效。adf_resid adfuller(resid, autolagAIC) print(残差ADF p值:, adf_resid[1]) if adf_resid[1] 0.05: print(残差非平稳先调ARIMA的d或换log变换)5.3 只看MAPE调参数高波动样本把模型带偏现象测试段某季度增速异常这个点的MAPE贡献了总误差一大半按MAPE最小化去调状态阈值所有阈值都往负偏差方向移动后面预测全部被往负数修正。原因MAPE对极端值和接近零的真实值过于敏感一个离群点就能主导调参方向。解决同时看RMSE和方向命中率定权重时用中位数误差而不是均值误差或者对极端修正做截断。这是典型的血泪经验单指标调参在时序预测里特别容易翻车。5.4 修正过头残差被修成恒等于零不是模型变强而是数据泄漏现象修正后MAPE下降90%画预测图和真实值发现曲线几乎完全贴合但换个时间窗口重新跑精度立刻崩回原样。原因状态中心是从训练残差均值算的如果验证段真实残差恰好落在某个固定状态附近期望残差回填会不断逼近真实误差相当于把测试期的残差分布信息偷偷透传进了修正项。解决严格保证状态中心和分位数阈值只来自训练段滚动回测里每个窗口只允许看到窗口之前的信息修正后检查残差序列的ACF如果出现明显负自相关说明修正项在吃真实残差模型已经退化成回填器。6. 验证与进阶白噪声检验、滚动窗口与修正阈值止损6.1 先跑白噪声检验什么时候马尔可夫修正没有意义马尔可夫修正的前提是ARIMA残差里还有信号。拟合后第一件事就是Ljung-Box检验p值大于0.05表示残差已经是白噪声没必要加修正p值小于0.05才有继续的价值。修正完成后同样要跑一次确认修正后的残差不再有自相关。from statsmodels.stats.diagnostic import acorr_ljungbox lb_arima acorr_ljungbox(model.resid, lags[6, 12], return_dfTrue) print(ARIMA残差Ljung-Box:, lb_arima[lb_pvalue].values) corrected_resid true_y - np.array(corrected) lb_markov acorr_ljungbox(corrected_resid, lags[6, 12], return_dfTrue) print(修正后残差Ljung-Box:, lb_markov[lb_pvalue].values)如果修正后p值比修正前更小说明马尔可夫链引入了新的自相关这不是修正过头就是状态划分有问题。我一般要求修正后残差至少不劣于修正前。6.2 滚动窗口复现避免一次性切分带来的运气成分一次性切分只能验证一个窗口容易把窗口选择的好坏当成模型好坏。滚动窗口的做法是每次用最近36期训练预测下一期然后窗口向后滑动一期重复执行最后把全部单步预测拼成完整预测序列。def rolling_arima_markov(series, order(1, 1, 1), window36, max_lag3, weightsNone): if weights is None: weights np.array([0.5, 0.3, 0.2]) preds [] for end in range(window, len(series)): train series.iloc[:end] m ARIMA(train, orderorder).fit() fit m.predict(start1, endend - 1, typlevels) resid train.iloc[1:].values - fit.values quants np.quantile(resid, [0.15, 0.35, 0.65, 0.85]) states np.digitize(resid, quants) mats build_transition_matrix(states, max_lagmax_lag) centers np.array([resid[states s].mean() for s in range(5)]) base m.predict(startend, endend, typlevels).values[0] prob np.zeros(5) for k in range(max_lag): if len(states) k 1: prob weights[k] * mats[k][states[-1 - k], :] prob / prob.sum() preds.append(base np.dot(prob, centers)) return np.array(preds)这段代码把第3章的完整流程压缩进了滚动迭代。注意states[-1 - k]取的是训练段末尾往前数第k1个状态mats[k]对应第k1阶转移矩阵两者匹配保证预测当前期时只用已知信息。窗口太短时状态划分不稳定我建议窗口至少36期并且每个滚动窗口都重新计算分位数阈值。6.3 阈值止损给修正加一道保险丝马尔可夫修正的问题是偶尔一个状态的期望残差特别大修正值直接把预测拉偏。我的做法是给修正值加一个上限上限取训练残差标准差的1到2倍超过就截断。这样既保留趋势修正能力又避免单期异常状态主导结果。def clip_correction(corr, resid_std, threshold1.5): limit threshold * resid_std if abs(corr) limit: return np.sign(corr) * limit return corr # 使用时 exp_resid np.dot(prob, state_centers) final_pred base_forecast[i] clip_correction(exp_resid, resid.std(), threshold1.5)threshold在1.0到2.0之间调。阈值太小修正能力被压制阈值太大又回到修正过头的风险。我一般先用1.5跑一轮滚动回测看修正项被截断的比例超过20%说明状态划分可能有问题。从那以后我每次跑ARIMA加马尔可夫链的组合都强制在日志里打印三行ARIMA残差Ljung-Box p值、修正后残差p值、修正项标准差与训练残差标准差之比。三者同时正常我才敢把预测结果交给业务方。希望帮到你。本文还有配套的精品资源点击获取