企业官网建设流程全解析

简介：这份源码包聚焦非平稳时间序列分析的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手动指定lags：

library(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 对缓慢变化的趋势更敏感。两者一起用，可以参照下表组合判断：

提示：同一份数据在多个检测函数之间出现不一致，不是代码坏了，而是备择假设不同。此时我一般以 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或-Inf，R 不会自动报错而是静默传播，后续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.test：

resid_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，再把结论带入滚动函数，判断才具备可比性。

本文还有配套的精品资源，点击获取