简介这份资源面向能源经济分析、金融量化及时间序列建模的学习者提供一套用R语言实现煤炭价格多元时序预测的完整案例。内容围绕VAR向量自回归模型展开涵盖平稳性检验ADF、PP、协整检验Johansen、滞后阶数确定AIC、BIC、HQIC、模型拟合、脉冲响应分析以及未来价格预测与置信区间评估等关键环节帮助读者理解多变量间的动态交互效应。资源包共3个文件包含1个R脚本、1个docx文档和1个csv数据文件压缩包约1.38MB分别对应建模代码、方法说明与历史价格数据便于直接复现与对照学习。目前已有361人学习下载。通过复现该案例读者可掌握R语言在时间序列分析中的实操流程理解煤炭价格的动态变化规律为政策研究、投资决策与企业运营提供量化依据。1. 煤炭价格多元时序预测从单变量到多因子的实战跨越很多做能源数据分析的朋友一开始都会陷入一个误区拿历史煤价跑个 ARIMA发现预测效果时好时坏一遇到政策调控或港口库存突变就彻底翻车。煤炭价格从来不是孤立波动的数字它背后是供需、运输、库存、下游电厂的日耗甚至包括进口煤到岸价的博弈。所谓多元时序预测就是把煤价这条主线和这些外生变量绑在一起建模让模型知道价格为什么涨、为什么跌。R language 在这个领域有天然优势forecast、tsibble、xgboost这些包能把特征工程和时序交叉验证串成一条流水线。这篇文章面向的是手里有煤价数据、想从单变量预测升级到多因子预测的从业者不管你是刚接触 R 的新手还是已经用过auto.arima的老手都能从下面的步骤里找到可复现的路径。2. 多元时序预测的底层逻辑与 R 生态选型2.1 为什么单变量模型在煤价预测上容易失效煤炭价格的形成机制决定了它不是一个自回归过程就能打发的。秦皇岛港的动力煤平仓价短期受日耗、库存、锚地船舶数影响中期受长协履约率和进口煤价差影响长期还受产能周期和宏观用电量牵引。单变量模型只从价格的历史序列里挖信息本质上是在假设“过去的波动模式会原样重复”。但煤价在迎峰度夏和供暖季的驱动因子完全不同一个在夏天有效的 ARIMA 参数到了冬天可能连残差白噪声检验都过不了。我自己的血泪经验是2021 年那波煤价从 500 多涨到 2500 多任何基于历史价格训练的纯时序模型都预测不出来因为驱动因素已经切换了。多元时序预测的核心价值就是把那些“切换”信号提前喂给模型。常见做法是构建一个包含目标变量煤价和多个外生变量的数据框外生变量包括港口库存、电厂日耗、进口煤到岸价、运费指数、期货升贴水等。这些变量不需要全部用上但至少要覆盖供需两端。R language 里处理这类问题的生态非常完整。tsibble负责把多列时间序列对齐成规整的 tbl 格式fable提供 ARIMA、ETS、VAR 等模型的统一接口xgboost和lightgbm则负责处理非线性关系和特征交互。选型的关键在于如果外生变量和煤价的关系近似线性且平稳用fable里的 ARIMA xreg 就够了如果关系复杂、存在阈值效应直接上梯度提升树做回归预测更稳。2.2 用 tsibble 把多源数据对齐成建模输入煤炭数据最头疼的是频率不一致。煤价是日频库存是日频但偶尔缺报日耗是日频但周末不更新进口煤价是周频。直接合并会出现大量 NA而大多数时序模型不接受缺失值。我的处理习惯是先统一到日频用tsibble的fill_gaps()补全时间索引再用interpolate()做线性插值最后用future_frame()生成预测期的外生变量占位。library(tsibble) library(dplyr) library(tidyr) # 假设已有三个数据框price_df煤价、inventory_df库存、consumption_df日耗 # 每个数据框包含 date 和 value 两列 # 统一为 tsibble 并补全日期索引 price_ts - price_df %% as_tsibble(index date) %% fill_gaps() %% mutate(price interpolate(price, date)) inventory_ts - inventory_df %% as_tsibble(index date) %% fill_gaps() %% mutate(inventory interpolate(inventory, date)) consumption_ts - consumption_df %% as_tsibble(index date) %% fill_gaps() %% mutate(consumption interpolate(consumption, date)) # 按 date 合并为多元 tsibble coal_ts - price_ts %% left_join(inventory_ts, by date) %% left_join(consumption_ts, by date) %% fill_gaps() %% mutate( inventory interpolate(inventory, date), consumption interpolate(consumption, date) ) # 查看缺失情况和时间跨度 coal_ts %% count_gaps() coal_ts %% head(10)这段代码的逻辑是先把每个数据源单独转成tsibble并补全日期避免合并时因为日期不连续产生错位。fill_gaps()会自动插入缺失的日期行interpolate()则对这些行做线性插值。合并之后再次fill_gaps()是为了处理不同数据源日期范围不一致的情况。参数上需要注意interpolate()默认对数值型列做线性插值如果你的数据有季节性可以改用model(ARIMA(value ~ season()))做插值预测但计算量会大很多。count_gaps()用来检查还有没有未补全的缺口如果有说明某个数据源的起始日期晚于其他源需要决定是截断还是向前填充。提示煤价数据在节假日经常停报插值会引入平滑效应建议在建模时加入一个“是否节假日”的哑变量来抵消这种影响。2.3 特征工程从原始外生变量到可入模的因子把数据对齐只是第一步真正决定预测上限的是特征质量。煤价的驱动因子往往不是原始值本身而是它们的变换形式。比如港口库存的绝对量意义不大库存同比变化率和库存可用天数才是交易员盯的指标。电厂日耗要看环比变化和与往年同期的偏离度。进口煤价要算内外贸价差因为价差决定进口煤是否冲击国内市场。我一般会构造以下几类特征差分特征一阶差分、七日差分、滚动统计特征7 日滚动均值、14 日滚动标准差、比率特征库存/日耗、内贸价/进口价、滞后特征外生变量的 1 到 7 阶滞后。在 R 里可以用dplyr的mutate配合lag()和rollmean()批量生成。library(zoo) coal_features - coal_ts %% mutate( # 差分特征 price_diff1 price - lag(price, 1), price_diff7 price - lag(price, 7), # 滚动统计 price_ma7 rollmean(price, 7, fill NA, align right), price_sd14 rollapply(price, 14, sd, fill NA, align right), # 比率特征 inventory_consumption_ratio inventory / consumption, # 外生变量滞后 inventory_lag1 lag(inventory, 1), consumption_lag1 lag(consumption, 1), # 同比偏离假设有年度数据这里用 365 天滞后近似 consumption_yoy consumption / lag(consumption, 365) - 1 ) %% drop_na() # 检查特征与煤价的相关系数 cor(coal_features %% select(-date), use complete.obs)[, price]这段代码的关键参数是滚动窗口长度。7 日对应一周的交易节奏14 日对应半个月的库存周期365 日对应同比。窗口不是越长越好煤价在旺季的波动周期可能只有 3 到 5 天窗口太长会把信号平滑掉。drop_na()会删掉因为滞后和滚动产生的缺失行如果数据量本身不大可以用fill()向前填充代替删除但要小心引入未来信息。相关系数矩阵用来快速筛查哪些特征和煤价线性关系强但不要只依赖这个非线性关系在 Pearson 系数上会显示为 0。3. 用 fable 和 xgboost 搭建两套可复现的预测流程3.1 ARIMA xreg 的多元回归预测流程fable包里的 ARIMA 模型支持通过xreg()传入外生变量这是最接近传统计量经济学的做法。它的假设是外生变量和煤价之间存在稳定的线性关系残差服从 ARIMA 过程。适合场景是数据量中等、外生变量维度不高3 到 5 个、关系相对平稳的情况。library(fable) library(tsibble) # 划分训练集和测试集最后 30 天作为测试 train - coal_features %% filter(date max(date) - 30) test - coal_features %% filter(date max(date) - 30) # 拟合 ARIMA xreg 模型 fit_arima - train %% model( arima_xreg ARIMA(price ~ inventory consumption inventory_consumption_ratio price_ma7 price_sd14) ) # 查看模型系数和残差诊断 report(fit_arima) fit_arima %% gg_tsresiduals() # 在测试集上预测 fc_arima - fit_arima %% forecast(new_data test) # 计算预测精度 accuracy(fc_arima, test)ARIMA()函数会自动选择差分阶数、AR 和 MA 的滞后阶数通过最小化 AICc 准则。xreg里的变量会作为回归项进入模型系数反映的是在其他变量不变时该变量每变化一个单位对煤价的影响。report()会输出系数估计和显著性检验如果某个外生变量的 p 值很大说明它和煤价的关系不显著可以考虑剔除。gg_tsresiduals()用来检查残差是否还有自相关如果残差 ACF 图有显著尖峰说明 ARIMA 部分没抓干净需要手动调整pdq参数。accuracy()输出的 RMSE 和 MAPE 是评估预测效果的核心指标煤价预测里 MAPE 能控制在 3% 以内就算不错了。注意xreg里的外生变量在预测期必须已知。如果预测未来 30 天你需要先预测这 30 天的库存和日耗或者用情景假设值。这是多元预测最容易被忽略的环节很多人模型拟合得很好一到实际预测就卡在外生变量的未来值上。3.2 xgboost 回归预测模型的调参与特征重要性当外生变量和煤价的关系呈现非线性、阈值效应或复杂交互时梯度提升树是更好的选择。xgboost在 R 里通过xgboost()函数调用核心是把时序问题转化为监督学习问题用过去 N 天的特征预测未来 1 天的煤价。library(xgboost) library(Matrix) # 构造监督学习矩阵用 t-1 到 t-7 的特征预测 t 时刻的 price build_lagged_matrix - function(data, target_col, feature_cols, n_lag 7) { result - data for (col in feature_cols) { for (i in 1:n_lag) { result[[paste0(col, _lag, i)]] - lag(result[[col]], i) } } result - result %% drop_na() return(result) } feature_cols - c(inventory, consumption, inventory_consumption_ratio, price_ma7, price_sd14, price_diff1) model_data - build_lagged_matrix(coal_features, price, feature_cols, n_lag 7) # 划分训练和测试 train_idx - 1:(nrow(model_data) - 30) train_x - model_data[train_idx, ] %% select(-date, -price) %% as.matrix() train_y - model_data$price[train_idx] test_x - model_data[-train_idx, ] %% select(-date, -price) %% as.matrix() test_y - model_data$price[-train_idx] # 训练 xgboost 模型 dtrain - xgb.DMatrix(data train_x, label train_y) dtest - xgb.DMatrix(data test_x, label test_y) params - list( objective reg:squarederror, eta 0.05, max_depth 4, subsample 0.8, colsample_bytree 0.8, min_child_weight 3 ) xgb_model - xgb.train( params params, data dtrain, nrounds 500, watchlist list(train dtrain, test dtest), early_stopping_rounds 30, verbose 1 ) # 预测与评估 pred_xgb - predict(xgb_model, dtest) rmse - sqrt(mean((pred_xgb - test_y)^2)) mape - mean(abs((pred_xgb - test_y) / test_y)) * 100 cat(RMSE:, rmse, MAPE:, mape, %\n) # 特征重要性 importance_matrix - xgb.importance(model xgb_model) print(importance_matrix[1:15, ]) xgb.plot.importance(importance_matrix[1:15, ])这段代码的核心参数需要重点解释。eta是学习率煤价数据噪声大建议设在 0.03 到 0.08 之间太大容易过拟合太小需要更多轮次。max_depth控制树的深度时序预测里一般 3 到 5 就够了太深会记住噪声。subsample和colsample_bytree是行采样和列采样比例0.7 到 0.9 之间能提升泛化能力。early_stopping_rounds配合watchlist用测试集监控如果 30 轮内验证误差不下降就停止这是防止过拟合的关键。min_child_weight控制叶子节点的最小样本权重煤价数据量通常不大设 3 到 5 比较稳妥。特征重要性输出里如果price_ma7_lag1排在第一说明煤价自身的短期动量最重要如果inventory_consumption_ratio_lag3排名靠前说明库存消费比的滞后效应显著。这些信息可以反过来指导 ARIMA 模型的外生变量选择。3.3 时序交叉验证与模型融合的实操单次划分训练测试集的结果波动很大煤价预测尤其如此。我一般用tsibble的stretch_tsibble()做滚动原点交叉验证模拟真实预测场景每次用截至 t 时刻的数据训练预测 t1 到 t7然后向前滚动。library(fabletools) # 滚动交叉验证初始训练 200 天每次滚动 7 天 cv_fit - coal_features %% stretch_tsibble(.init 200, .step 7) %% model( arima_xreg ARIMA(price ~ inventory consumption inventory_consumption_ratio price_ma7), snaive SNAIVE(price) ) %% forecast(h 7) # 计算交叉验证精度 cv_accuracy - cv_fit %% accuracy(coal_features) print(cv_accuracy) # 模型融合取 ARIMA 和 xgboost 预测的加权平均 # 假设已有 fc_arima 和 pred_xgb 对齐后的预测值 ensemble_pred - 0.4 * fc_arima$.mean 0.6 * pred_xgbstretch_tsibble()的.init参数控制初始训练窗口长度煤价数据建议至少 200 个交易日覆盖一个完整的季节性周期。.step是滚动步长7 天对应一周如果数据量充足可以设 1 天做更细粒度的验证。accuracy()会输出每个滚动窗口的 RMSE、MAE、MAPE看均值的同时也要看标准差如果某些窗口误差特别大说明模型在特定市场状态下失效。融合权重 0.4 和 0.6 不是固定的可以用验证集上的误差倒数做加权ARIMA 在平稳期表现好xgboost 在波动期更稳动态权重能兼顾两者。4. 煤炭价格预测的避坑与排查清单4.1 外生变量未来值未知导致的预测断裂现象模型在训练集和测试集上精度都很好但一到实际预测未来一个月就报错或输出离谱结果。原因xreg或特征矩阵里的外生变量在预测期没有真实值代码里用了NA或者错误地重复了历史值。解决提前建立外生变量的预测流程或者用情景分析法设定乐观、中性、悲观三套外生变量路径分别输出煤价预测区间。我一般会先用fable对库存和日耗单独做预测再把预测值作为new_data传入煤价模型。4.2 数据频率不一致导致的隐性错位现象合并后的数据框行数对不上或者某些日期的外生变量全是 NA。原因不同数据源的日期格式不统一有的是Date类型有的是字符串或者某个数据源用的是交易日历另一个用的是自然日。解决合并前统一用as.Date()转换日期列用tsibble的fill_gaps()检查缺口。如果某个数据源只有交易日数据需要决定是向前填充到自然日还是把煤价也截断到交易日。我的习惯是统一到自然日交易日缺失用前值填充但加一个is_trading_day哑变量标记。4.3 滞后特征引入的未来信息泄露现象交叉验证精度高得离谱MAPE 低于 1%但实盘预测一塌糊涂。原因构造滞后特征时用了lag()但方向搞反了或者滚动统计用了align center导致当前值参与了计算。解决所有滞后特征必须确保 t 时刻只能看到 t-1 及之前的数据。rollmean()的align参数必须设为rightlag()的n必须为正数。构造完特征后随机抽几行手动核对t 行的price_ma7应该等于 t-6 到 t 的均值而不是 t-3 到 t3。4.4 煤价政策干预导致的模型失效现象模型在正常市场环境下误差可控但遇到限价政策、进口配额调整等事件时预测完全偏离。原因政策干预是外生冲击不在历史数据模式里任何基于历史统计的模型都无法提前预知。解决在模型里加入政策哑变量比如“是否处于限价区间”“进口煤通关是否受限”用 0/1 标记。同时建立人工干预机制当监测到政策发布时手动调整预测值或切换到情景分析模式。不要指望模型能自动学会政策影响样本量根本不够。4.5 xgboost 参数过拟合的隐蔽表现现象训练集 RMSE 极低测试集 RMSE 是训练集的 3 倍以上。原因max_depth设得太大、nrounds没有用早停、特征维度太高但样本量不足。解决先用xgb.cv()做交叉验证确定最佳nrounds再把max_depth限制在 5 以内min_child_weight提高到 5 以上。如果特征数量超过样本量的十分之一先做特征筛选用xgb.importance()去掉排名后 50% 的特征再重新训练。5. 用滚动预测和情景分析把模型推到生产可用模型在测试集上跑出好看的 MAPE 只是第一步真正到生产环境里煤价预测的交付物不是一个点估计值而是一个带置信区间的预测路径。我现在的习惯是每次预测输出三列——乐观值、中性值、悲观值分别对应外生变量的三种情景。中性情景用模型预测的外生变量均值乐观和悲观用 ±1 个标准差。这样业务方拿到的不只是一个数字而是一个决策区间。具体操作上我会把fable的 ARIMA 预测和xgboost的预测做分位数融合。ARIMA 的forecast()输出自带 80% 和 95% 置信区间xgboost 则用分位数回归或者对残差做 bootstrap 得到区间。两者加权后用geom_ribbon()画出一条带阴影的预测带。这个图比任何精度指标都更能说服业务方因为它直观展示了“最坏情况可能到哪里”。还有一个容易被忽略的技巧把预测结果和实际值做滚动回测每周更新一次模型。煤价的市场结构在变三年前有效的特征今天可能已经失效。我一般每季度重新跑一次特征重要性排序如果发现某个外生变量的重要性持续下降就把它从模型里剔除。模型不是一次训练就一劳永逸的它需要跟着市场一起迭代。最后说一个我踩过的坑曾经花了很大精力把 MAPE 从 2.8% 优化到 2.3%结果业务方根本感知不到差别他们更关心的是“预测方向对不对”和“极端行情能不能提前预警”。所以后来我把评估指标从 MAPE 改成了方向准确率加上极端值召回率模型调优的目标也跟着变了。技术指标要服务于业务目标这个道理在煤价预测里尤其明显。希望帮到你。本文还有配套的精品资源点击获取