简介这份源码包聚焦非平稳时间序列分析的R语言实现适合经济、金融、工程等领域需要处理股票价格、销售记录、气象数据等非平稳序列的数据分析初学者与进阶者。资源以R脚本为主共5个文件包含4个R代码文件和1个Rhistory命令历史文件压缩包仅5KB轻量实用方便直接对照运行。内容覆盖时间序列对象构建、可视化、描述性统计、差分平稳化、自相关与偏自相关分析、ADF单位根检验、ARIMA建模、季节性分解及预测等核心环节并配有上课笔记代码和老师课后更新内容便于按学习进度对照理解。目前已有492人学习下载通过实际运行这些脚本读者能快速上手R中非平稳时间序列的完整分析流程掌握从数据预处理、平稳性检验到模型构建与预测的完整链路同时加深对统计理论与实际数据操作之间关系的理解。1. 非平稳时间序列分析的第一个 R 代码步骤先承认数据可能不平稳非平稳时间序列平时最容易出现在两类场景里一类是带明显趋势的经济与金融数据比如股价对数收益外的原始收盘价、GDP 季度值另一类是带季节波动的工厂能耗、交通流量和气象观测。很多入门教程在讲到平稳性时往往只用一张“前后图像对比”带过但真正落到 R 代码上时先要回答的问题不是“能不能建模”而是“用什么检验去证明它不是平稳的”。R 里围绕这套判断的常用工具集中在tseries、forecast、urca三个包里对五年以上经验的工程师来说这些函数本身没有什么可记的但滞后阶数怎么选、两个方向相反的检验同时出现矛盾结果时该信谁、差分之后序列长度变化对后续预测代码的影响属于真正值得掰开的部分。这篇内容围绕“R代码_非平稳时间序列分析_源码”这个标题按一套可复现的源码路径把检验、变换、建模和滚动验证串起来讲。2. 用单位根检验先给非平稳性定性adf.test、kpss.test、ur.df 的滞后与输出解读2.1 ADF 检验在 R 中的代码形态与拒绝域边界先造一组明显带趋势的模拟数据目的是让后续每个函数的输出都能对照直觉。下面这段源码模拟的是一个从 2010 年开始、每月一采样的 200 期随机游走加漂移序列set.seed(2024) x_trend - ts(cumsum(rnorm(200, mean 0.15, sd 1)), frequency 12, start c(2010, 1)) plot(x_trend, main non-stationary series)rnorm的均值取 0.15 表示每一步都叠加正向漂移cumsum做累加后序列在视觉上一定是向上的直线形态任何看图的初判都会指向“非平稳”。接下来跑 ADF 检验library(tseries) adf_result - adf.test(x_trend) print(adf_result)adf.test的默认原假设是序列存在单位根也就是非平稳。如果输出里的p-value大于 0.05就没有理由拒绝原假设此时应把序列按非平稳处理。模拟代码跑出来的典型结果是Dickey-Fuller -2.31, Lag order 4, p-value 0.45这个数值足够说明问题ADF 的检测统计量必须比临界值更“负”才拒绝单位根而 -2.31 落在接受域里。这里要提醒一个 R 特有的细节adf.test的默认滞后阶数是floor(length(x) - 1)^(1/3)一类基于样本量的启发式结果但它不一定是数据生成过程真实的 AR 阶数。滞后太少残差的自相关会污染统计量滞后太多检验功效下降。更可控的做法是用urca::ur.df手动指定lagslibrary(urca) ur_df_trend - ur.df(x_trend, type trend, lags 5) summary(ur_df_trend)type trend表示回归式中同时包含常数项和趋势项这对应带漂移的非平稳序列。从summary输出中看tau3统计量即可它和adf.test里的Dickey-Fuller值含义一致。两层代码叠加还不是关键更关键的是 ADF 检验对“确定性趋势”和“随机趋势”并不敏感所以当type trend与默认检验结果矛盾时通常要去参考 KPSS 的结果。2.2 KPSS 检验与 ADF 方向相反的第二道确认tseries::kpss.test的原假设是序列平稳备择假设是存在单位根或趋势。把 2.1 的序列放进去kpss_result - kpss.test(x_trend) print(kpss_result)输出结果通常是KPSS Level 1.742, Truncation lag parameter 4, p-value 0.01。由于 p-value 很小拒绝“平稳”的原假设因此结论与 ADF 一致。两套检验真正有意思的地方在于它们方向和侧重点不同ADF 是“有单位根则拒绝平稳”KPSS 是“有单位根则拒绝平稳”但 KPSS 对缓慢变化的趋势更敏感。两者一起用可以参照下表组合判断ADF 结论KPSS 结论常见情况拒绝单位根不拒绝平稳序列本身平稳进入常规 ARMA 建模不拒绝单位根拒绝平稳典型非平稳先差分或做去趋势拒绝单位根拒绝平稳可能是结构性断点或长记忆过程不拒绝单位根不拒绝平稳样本量不足或两检验功效都弱需要更多数据提示同一份数据在多个检测函数之间出现不一致不是代码坏了而是备择假设不同。此时我一般以 KPSS 的结果为准因为 KPSS 对趋势型偏离的检测更直接而 ADF 的小样本功效确实偏弱。2.3 直接把三组检验封装成一个判断函数把常见的两步判定封装成函数方便后续在批量序列上复用。这个函数返回差分阶数建议同时也保留每一层检验的 p-value便于排查问题decide_diff_order - function(series, max_d 2) { require(tseries, quietly TRUE) for (d in 0:max_d) { adf_p - tryCatch(adf.test(series)$p.value, error function(e) NA) kpss_p - tryCatch(kpss.test(series)$p.value, error function(e) NA) if (is.na(adf_p) || is.na(kpss_p)) break if (adf_p 0.05 kpss_p 0.05) { return(list(diff_order d, adf_pvalue adf_p, kpss_pvalue kpss_p)) } series - diff(series) } return(list(diff_order max_d, note reached max diff)) } decide_diff_order(x_trend)max_d限制了最多差分几次防止对白噪声序列做无意义的过度差分每一次循环如果 ADF 的 p-value 小于 0.05 且 KPSS 的 p-value 大于 0.05就认为当前差分阶数已经让序列进入平稳状态。实际使用中如果遇到返回reached max diff要回过头检查原始数据是否有异常缺失、恒定值段或者突变跳跃而不是简单提高max_d。3. 差分、对数与 STL 分解非平稳序列在 R 代码里的三类转换路径3.1 用 ndiffs() 决定差分阶数再手动 diff() 确认长度变化检验只是给结论真正动手操作时第一选择是差分。差分在 R 里直接用diff()但差分的阶数不要靠肉眼猜。forecast::ndiffs()提供两种常见的自动化判断library(forecast) d_auto - ndiffs(x_trend, test adf) d_kpss - ndiffs(x_trend, test kpss) print(paste(d_auto, d_kpss))test adf时内部调用 ADF 检验test kpss时调用 KPSS 检验结合 2.3 的判断思路通常建议使用 KPSS 版本因为它对趋势的识别更主动。得到d 1后执行x_diff1 - diff(x_trend, lag 1, differences 1) plot(x_diff1)这里有一个容易忽略的 R 行为diff()不会保留时间序列的起始点x_trend有 200 个点x_diff1会变成 199 个点。如果后面要把差分序列和其他变量放进同一个data.frame长度不一致会导致模型报错。处理办法一般是先把原始时间标签对齐或者保留一个ts对象并显式重建x_diff_ts - ts(x_diff1, frequency frequency(x_trend), start time(x_trend)[2])frequency(x_trend)取回原始月份频率start从第二个观测点起算这样后续预测代码的h参数才能基于正确的时间索引工作。3.2 对数变换的适用边界方差稳定、负值与季节性幅度非平稳序列里如果振幅随水平值同步变大例如客流高峰月数值是平时的 3 到 5 倍且波动幅度也成倍增加这通常暗示方差与均值相关。此时先做对数变换再差分效果优于直接对原始值差分x_log - log(x_trend) x_log_diff - diff(x_log, lag 1, differences 1)对数变换的意义有两个层面。第一它把乘法关系变成加法关系对诸如“夏季销量是冬季的 2 倍”这类季节效应取对数后季节效应从倍数变成常量更适合用线性模型表达。第二对数差分近似等于收益率或变化率这能让序列在量纲上比原始值更容易比较。注意原始序列存在 0 或负数时log()会产生NaN或-InfR 不会自动报错而是静默传播后续acf()、arima()都会以莫名的方式失败。遇到这种情况需要先做any(x_trend 0)检查否则先考虑加常数平移或改用其他变换。3.3 STL 分解在 R 中的频率参数与稳健性设置差分适合趋势型非平稳但季节性强的序列建议先做 STL 分解把趋势、季节、剩余拆开。STL 的优势是不要求数据均匀方差且能够处理随时间变化的季节强度。R 中直接用stl()x_season_ts - ts(rnorm(120, 10, 2) seq(1, 120) * 0.1 rep(c(0, 2, -1, 3), each 30), frequency 12) decomposed - stl(x_season_ts, s.window periodic, robust TRUE) plot(decomposed)s.window periodic表示季节成分在整个序列上保持周期恒定适合年度规律稳定的情况。如果季节模式本身会逐年漂移需要给s.window一个奇数数值比如 13让局部回归窗口仅随邻近周期更新。robust TRUE会使用低权重抵御异常点这对传感器数据或活动数据里的尖峰特别有效。分解之后要建模的核心序列可以从分解对象中提取season_adj - decomposed$time.series[, trend] decomposed$time.series[, remainder]去掉 seasonal 列后剩余序列通常更容易通过 ADF 检验。不过要注意的是STL 分解得到的 “trend” 并不是平稳的所以这只能算预清洗不能替代差分。4. 非平稳序列进入 SARIMA 建模auto.arima 的 d、D 与残差诊断4.1 auto.arima 在源码层面的默认搜索逻辑与参数收紧经过检验和变换数据已经可以作为预测模型输入。R 里对非平稳序列最主流的建模入口是forecast::auto.arima()它的搜索逻辑大致是先用 KPSS 或 ADF 自动定d再通过 AICc 对 ARMA 部分做逐步搜索。直接跑一行代码虽然方便但生产环境里建议收紧参数library(forecast) fit_sarima - auto.arima(x_log_diff, max.p 5, max.q 5, stepwise FALSE, approximation FALSE, seasonal FALSE) summary(fit_sarima)stepwise FALSE让 R 不做逐步搜索而是把候选模型空间充分遍历代价是运行时间变长。对 200 到 500 个观测点的序列这个时间完全值得。approximation FALSE则关闭近似的快速估计保证信息准则值按最大似然计算。但这里有一个非常微妙的点当我们把x_log_diff传给auto.arima时实际上差分已经人工完成auto.arima不会再做一次差分。如果传原始序列并让auto.arima自己决定d它会内部调ndiffs()但这种隐式处理会在预测阶段自动还原差分。从预测准确性上看两条路径等价从调试角度看显式差分后建模更可控因为能看到残差里是否残留趋势。4.2 季节性非平稳序列的 D 参数同频差分与 SARIMA如果有明显的季节周期模型升级成 SARIMA对应参数是order和seasonal两个向量。对月度数据常见的设置如下fit_sarima - auto.arima(x_season_ts, order c(2, 1, 2), seasonal c(1, 1, 1, 12), stepwise TRUE, trace TRUE)seasonal c(1, 1, 1, 12)表示季节自回归阶数为 1、季节差分阶数为 1、季节移动平均阶数为 1周期为 12。这里的季节差分D 1作用是消除年度为单位重复出现的非平稳周期与普通差分d 1消除的短期趋势不同。trace TRUE会把每一步候选模型的 AICc 打印出来。这是 R 源码应用里常被忽略的调试利器当预测结果不理想时你能回看模型选择路径判断是 AR 阶数受限还是季节项搜索过早停止。如果auto.arima的结果残差仍然不白噪声就继续看下面的诊断代码。4.3 Ljung-Box 残差检验是模型边界判断的最后一环模型拟合完成后必须验证残差不再包含可建模的自相关结构。tsdiag提供三张图但更精准的是Box.testresid_sarima - residuals(fit_sarima) Box.test(resid_sarima, lag 12, type Ljung-Box, fitdf length(coef(fit_sarima)))fitdf指定消耗的自由度数量用来修正因估计参数带来的检验偏移这个参数在 R 里容易被漏掉。p-value 大于 0.05 意味着残差没有显著自相关模型信息被充分提取。如果 p-value 低于 0.05建议检查ACF哪个滞后阶还在置信区间之外然后手动增加对应p或q阶数。关于d和D的选择有一点要明确过度差分会降低序列方差使原本可识别的 AR 结构被抹平。auto.arima默认不会同时做超过一次的普通差分和季节差分但人工建模时有可能出现d 2, D 1的情况这时一定要对比差分前后的方差变化避免为了形式上的平稳牺牲预测方向。5. 在非平稳结构下做滚动原点预测验证单次划分训练集和测试集对非平稳序列的验证效果有限因为序列生成机制可能在中期就改变了。更可靠的方式是用forecast::tsCV()做滚动原点交叉验证每次只向前预测一步、两步或多步然后移动训练窗口重新拟合并记录误差。关键是理解tsCV的返回值是按“预测原点”对齐的误差矩阵不是直接给出汇总指标library(forecast) farima_model - function(y, h) { fit_temp - auto.arima(y, seasonal FALSE, stepwise FALSE) forecast(fit_temp, h h) } cv_errors - tsCV(x_log_diff, farima_model, h 1) rmse_rolling - sqrt(mean(cv_errors^2, na.rm TRUE)) print(rmse_rolling)cv_errors的长度与输入序列一致但最后一部分会因为窗口过短产生NA所以计算时使用na.rm TRUE。h 1时验证的是短期预测能力如果需要评估多步预测比如接下来 3 个月或 6 个月的趋势就把h调大并且在此基础上计算每个预测步长的 RMSE 平均值cv_errors_h3 - tsCV(x_log_diff, farima_model, h 3) rmse_h3 - apply(cv_errors_h3, 1, function(row) { sqrt(mean(row^2, na.rm TRUE)) })apply(cv_errors_h3, 1, ...)是按行计算每个预测原点的多步平均误差。这个结果能直接画成折线图观察误差是否随预测步长增长而快速膨胀——对非平稳序列来说误差增长的斜率比误差本身更有信息量斜率越陡说明模型结构稳定性越差。滚动验证之外还有一个常用技巧把差分阶数和季节参数固定在验证阶段不让每次循环重新搜索。因为auto.arima每做一次搜索都会消耗不少时间更关键的是它会针对不同训练窗口选择不同d这会造成误差指标失真。实务中先在整个样本上用一次ndiffs锁定d再把结论带入滚动函数判断才具备可比性。本文还有配套的精品资源点击获取