时间序列建模完整闭环：从平稳性检验到SARIMA预测

简介人大王燕版《时间序列分析》课后习题参考答案覆盖第二章与第三章核心内容并包含上机操作题目解答适合统计学、数据科学及相关专业学生备考与自学。文档针对自相关系数计算、偏自相关图识别、AR模型、MA模型、ARIMA模型建模等高频考点逐题推演既有理论推导也有R语言与Excel实现细节可帮助读者补齐课堂练习中的易错环节理解平稳性检验与模型定阶思路。资源包仅含1个doc文档文件约1.2MB便于直接查看、打印或导入笔记目前已有350人学习下载是时间序列课程复习与考研准备中较为实用的一份补充材料。配套内容从第二章序列平稳性判断、纯随机序列检验到第三章模型平稳性与可逆性证明、Green函数计算再到指数平滑与趋势拟合的上机示例基本串起主干知识点能有效提升课后习题的完成效率。1. 从一套习题答案看时间序列建模的完整闭环拿到这份《人大(王燕)时间序列课后习题答案(2_5,含上机的)》时我正在帮人调一个 R 脚本原始序列明明在 adf.test 里 P 值小于 0.05可自相关图上的系数却拖了十几阶都不衰减。后来发现问题不是检验方法选错而是没分清“原序列”和“差分后序列”到底该用哪一个。这份答案的价值恰好在这里它把第二章的平稳性判断、第三章的 ARMA 模型识别、第四章的差分与 SARIMA 上机代码串成了一条完整链路。对准备复试上机、考研复习或者用 R 做时间序列预测但总在定阶上翻车的人来说值得逐题过一遍。2. 平稳性判断自相关系数、LB 统计量与单位根检验的配合2.1 自相关系数衰减慢说明什么文档第二章第 3 题的原话是“因为序列具有明显的趋势所以序列非平稳”。这句批注看着简单却是整个建模的起点。弱平稳要求均值、方差不随时间改变趋势的存在直接破坏了“均值恒定”这条假设。更直接的证据是自相关系数0.850、0.702、0.556、0.415、0.280、0.153衰减速度很慢到第 8 阶才勉强回到正负两倍标准差范围内。这种慢衰减意味着当前观测与很久以前的观测仍然相关典型的一阶单整序列特征。在 R 里复现这个过程时不要直接对原始序列调用 acf应该先画时序图再决定是否差分。我一般会这样写z - scan(习题4.6数据.txt) ts.plot(z) acf(z, lag.max 20, main 原始序列自相关图)acf() 返回各阶自相关系数和置信区间置信区间默认按 1.96/sqrt(n) 计算。如果前几阶相关系数远在区间之外而且衰减缓慢基本可以断定序列非平稳。注意文档里的括号数值来自近似公式比如 0.85 旁边的 0.85 用的是递推近似实际用 R 算出来会因为样本量和中心化处理不同而有小数点后的差异没必要强行对齐。这里最容易犯的错误是看到自相关图“有相关性”就急着建 ARMA 模型却忘了先处理趋势。趋势没去掉之前自相关图的显著相关只是趋势的投影不是真正可建模的短期相关结构。2.2 LB 统计量判断纯随机序列第 4 题直接给出 LB(6)1.6747、LB(12)4.9895对应 0.05 显著性水平下的临界值 12.59 和 21.0。比较后结论是“该序列为纯随机序列”。这里有个隐含前提这个 LB 检验是针对已经处理过的序列做的不是原始趋势序列。如果你拿带趋势的原始序列去跑 Box.test结果几乎必然显著但那个“非白噪声”只是趋势造成的假象。R 命令是Box.test(z, lag 6, type Ljung-Box) Box.test(z, lag 12, type Ljung-Box)type 参数选 Ljung-Box 而不是 Box-Pierce因为小样本下 Ljung-Box 的卡方逼近更准。lag 对应要检验的最大延迟阶数习题里取 6 和 12是为了覆盖短期相关性和一个完整季节周期的长度。P 值大于 0.05 时没有充分证据拒绝“纯随机”的原假设说明序列已经是白噪声不需要再建模。2.3 adf.test 与肉眼判读的边界文档 4.8 题开始引入单位根检验library(tseries) adf.test(z)adf.test 的原假设是序列有单位根即非平稳。P 值小于 0.05 时拒绝原假设可以认为序列平稳。这里要留意两个问题。第一adf.test 默认的回归形式包含常数项和趋势项如果序列本身没有确定趋势检验功效会下降结论更保守这种情况我一般会再用 urca 包的 ur.df(z, typedrift) 确认一次。第二adf.test 的滞后阶数 k 默认是 trunc((length(z)-1)^(1/3))样本量变化时 k 会跟着变不同 k 下结论可能不一致。所以我把 ACF 图、时序图和 adf.test 放在一起看三个证据互相印证而不是只看一个 P 值。工具原假设判断标准适用场景时序图无肉眼判断均值是否恒定建模前第一步ACF 图无自相关系数是否快速衰减到置信区间内判断相关性结构Ljung-Box序列是白噪声P 值 0.05 接受白噪声残差诊断adf.test存在单位根P 值 0.05 拒绝原假设定量判断平稳性在文档 4.8 题里adf.test 的结果与自相关图一致所以后文直接开始拟合。实际操作中如果两者矛盾优先检查是否漏了差分、数据里是否有缺失值或者 adf.test 的滞后阶数是否选得太大。这些是上机时最常见的翻车点。3. ARMA 模型识别ACF/PACF 截尾拖尾与可逆性验证3.1 AR(2) 平稳域与特征根第三章第 3 题给出了一个 AR(2) 模型第 4 题问参数 c 满足什么条件时模型平稳答案是“当 -1c0 时该 AR(2) 模型平稳”。这个结论不是靠猜的而是来自特征方程。AR(2) 模型 x_t φ1 x_{t-1} φ2 x_{t-2} ε_t 的特征方程为 λ² - φ1 λ - φ2 0只有当两个特征根的模都小于 1 时模型才平稳。R 里可以直接用 polyroot 求根phi1 - 0.5 phi2 - 0.3 roots - polyroot(c(1, -phi1, -phi2)) Mod(roots)polyroot(c(1, -phi1, -phi2)) 返回特征方程 λ² - φ1 λ - φ2 0 的两个根。c() 里第一个元素是常数项 1第二个是 -φ1第三个是 -φ2顺序不能写反。Mod() 取复数的模只要两个模都小于 1模型就平稳。习题答案里出现“0.450.2693i”这样的输出就是在说复根复根的模同样按实部平方加虚部平方开根号判断。很多人在这一步把平稳性条件记成“φ1φ21, φ2-φ11, |φ2|1”但用的时候不注意方向。其实这三个条件等价于特征根都在单位圆内直接求根更不容易错。文档第 5 题还证明“不论 c 取何值都会有一个特征根等于 1”这种情况下模型必然非平稳这就是单位根问题。3.2 MA 模型可逆性判别第 8 题是一道经典的 MA(1) 参数题模型可写成某种算子形式当 1-0.5B 能够整除相关多项式时模型等价于 MA(2)解得 C0.275。这里涉及的核心概念是“可逆性”。MA 模型本身总是平稳的但只有满足可逆性条件才能写成等价的自回归形式预测和参数估计才稳定。MA(1) 模型 x_t ε_t θ ε_{t-1} 的可逆性要求 |θ|1。判断高维 MA 是否可逆同样可以用 polyroot只是根要落在单位圆外theta - c(1, -0.275) # 对应算子多项式 1 - 0.275B roots - polyroot(theta) Mod(roots)这里的 c(1, -0.275) 对应算子多项式 1 - 0.275B 的系数。可逆性要求所有根的模都大于 1。文档第 11 题里“模型可逆”的结论都是先算特征根再比较模得到的。实际建模时如果 arima() 拟合出的 MA 系数绝对值接近 1就要警惕模型可能不可逆需要换阶数或者对序列做进一步变换。3.3 ACF/PACF 定阶与 Green 函数预测区间习题 4.7 中有一段非常典型的定阶过程观察 12 步差分后的 ACF看到 12 阶处显著、24 阶处回到置信区间判定季节部分用 MA(1)同时 ACF 在 1 阶处显著后立刻截尾PACF 在第 5、6 阶出界综合得到乘积模型。这段逻辑很多人看不明白核心是两条ACF 截尾对应 MA 项PACF 截尾对应 AR 项季节阶数要看 lag12、24、36 处的表现是否构成规律观察位置ACF 表现PACF 表现建议模型非季节部分1 阶截尾缓慢拖尾MA(1)非季节部分拖尾1 阶截尾AR(1)季节部分12 阶截尾12、24、36 递减季节 MA(1)季节部分12、24 递减12 阶截尾季节 AR(1)定阶之后文档还用 Green 函数计算了第三章第 16 题的预测区间。AR(1) 模型的 Green 函数是 G_j φ^j预测误差方差是 σ² 乘以 G_j 的平方和。所以预测区间就是点预测加减 1.96 倍标准差。R 里不用手动算predict() 会直接给fit - arima(z, order c(1, 0, 0)) p - predict(fit, n.ahead 5) p$pred 1.96 * p$se p$pred - 1.96 * p$sepredict() 返回的 se 是预测标准差按正态分布近似计算区间。文档里出现 9.9892±1.96×3.1419 这类结果就是套这个公式。n.ahead 表示向前预测几步超过 1 步后 AR 系数会通过递推放大不确定性区间会变宽这是正常现象。4. 非平稳序列建模差分、季节分解与 SARIMA 乘积模型4.1 趋势拟合与移动平滑的选择习题 4.6 给的是速度数据时序图呈现二次曲线形状。文档用了两种做法一种是直接做二次曲线拟合 lm(income~tt2)另一种是 N5 的简单移动平均 filter(income, rep(1/5,5), sides1)。两种方法都能提取趋势但逻辑不同。回归拟合适用于趋势形态明确的序列移动平滑更通用不需要预设函数形式但会损失首尾数据。R 代码income - scan(习题4.6数据.txt) t - 1:length(income) t2 - t^2 fit_lm - lm(income ~ t t2) summary(fit_lm) income_fil - filter(income, rep(1/5, 5), sides 1) ts.plot(income) lines(fit_lm$fitted.values, col 2) lines(income_fil, col 3)lm() 里的 t 和 t2 是时间趋势项和二次项summary() 看 t2 的 P 值是否显著若显著说明二次项有必要。filter() 的 rep(1/5,5) 构造等权重窗口sides1 表示只用过去值计算避免未来数据泄漏。这里有个细节sides1 时前 4 个点无法计算R 默认给 NA后续画线和预测都要处理这些 NA否则下标会错位。回到建模流程如果趋势明确且有理论背景优先用回归如果只是想快速看趋势走势移动平滑更稳。习题 4.7 里既有季节又有趋势单靠这两种方法都不够所以才引入了季节分解。4.2 decompose 加法模型与季节系数提取习题 4.7 的序列是月度数据时序图上能同时看到长期递增和以年为周期的季节波动。文档选择加法模型 z trend seasonal random并用 decompose(z, typeadditive) 分解。加法模型适合季节波动幅度不随趋势水平变化的数据如果波幅随着水平增大而增大应该改用 typemultiplicative。分解后的关键一步是“减掉季节系数再预测”z - ts(z, start c(1962,1), frequency 12) z_dec - decompose(z, type additive) z_adj - z - z_dec$seas t - 1:length(z_adj) r3 - lm(z_adj ~ t I(t^2) I(t^3)) pt - (length(z_adj)1) : (length(z_adj)12) pred - r3$coef[2:4] %*% rbind(pt, pt^2, pt^3) r3$coef[1] pred_season - pred z_dec$sea[1:12]decompose() 返回的 seas 是季节系数长度为一个周期加法模型里它围绕 0 波动。z_adj 去掉季节成分后只剩趋势和随机部分再用三次多项式拟合最后把季节系数加回去。这里的矩阵乘法 %*% 左边是三次项系数右边是未来 12 个时刻的 t、t²、t³拼成 3×12 矩阵。coef 向量里第 1 个是截距后面三个是 t、t²、t³ 的系数。这个做法不是万能的。decompose() 假定季节成分固定不变如果季节模式随年份缓慢变化这种“一次性提取”会引入滞后。文档里最后又用 SARIMA 重新建模本质上就是为了缓解这个问题。4.3 乘积 SARIMA 定阶与参数乘积模型写成 arima(z, orderc(0,1,1), seasonallist(orderc(0,1,1), period12))含义是先做一阶差分消除趋势再做 12 步差分消除季节然后对差分后的序列拟合非季节 MA(1) 和季节 MA(1)。order 的三个数分别对应 (p,d,q)seasonal 里的 order 是季节部分的 (P,D,Q)period12 表示月度周期。为什么不直接对原序列拟合 ARMA因为既有趋势又有季节效应的序列直接拟合残差很难通过白噪声检验。文档里明确写了“尝试拟合 ARMA 模型但拟合效果均不理想拟合残差均通不过白噪声检验”这是非常典型的教训。先差分再建模比强行加大 ARMA 阶数更干净。定阶依据来自差分后序列的 ACF/PACFsq - diff(diff(z), lag 12) acf(sq, lag.max 50) pacf(sq, lag.max 50)diff(z) 做一阶差分diff(..., lag12) 再做 12 步差分两步把趋势和季节同时去掉。文档里 ACF 在 1 阶处截尾说明非季节部分用 MA(1)季节部分在 12 阶截尾说明用季节 MA(1)最终阶数是 (0,1,1)×(0,1,1)_12。实际定阶时“截尾”和“拖尾”的判断很主观我的建议是列出一组候选阶数逐个拟合并比较 AIC不要只看图。4.4 arima() 的 fixed 参数处理不显著项习题 4.8 的序列经过差分后仍有些不明显的周期性文档尝试了几个 SARIMA 组合最终选中 orderc(4,1,0)、seasonallist(orderc(1,0,1), period12)并且用 fixed 把第三个 AR 系数固定为 0result - arima(z, order c(4,1,0), seasonal list(order c(1,0,1), period 12), fixed c(NA, NA, 0, NA, NA, NA))fixed 向量的长度必须和模型参数个数一致。这个模型的参数顺序是 ar1、ar2、ar3、ar4、sar1、sma1共 6 个所以 c(NA,NA,0,NA,NA,NA) 表示只把 ar3 强制为 0其余 5 个参数自由估计。NA 代表“不固定”数字代表“固定为该值”。使用 fixed 的场景是某个中间滞后阶数不显著但删掉整个 AR(4) 会丢失其他显著项。这样操作后自由度减少AIC 会变化必须重新比较。文档里说“第三个的 aic 值最小即模型拟合效果最好”就是这个流程。注意 fixed 是 R 的 arima() 特有的参数在 Python 的 statsmodels 里没有直接对应项需要手动构造约束所以跨语言复现时容易搞混。5. 上机题的验证技巧残差白噪声检验与滚动预测5.1 残差白噪声检验的坑拟合完模型只看得分不验残差等于白干。tsdiag(result) 会输出三张图标准化残差、残差 ACF、Ljung-Box 检验的 P 值。文档里强调“延迟 6 阶和 12 阶的 P 值均大于 0.05可以认为该残差序列即为白噪声”。但很多上机题会在这里扣分原因是直接调 Box.test 时没有设置 fitdf。模型消耗了几个参数自由度就要减几个否则 P 值会偏大本来不白噪声的残差也可能被误判成白噪声。我一般会写成Box.test(result$residuals, lag 12, type Ljung-Box, fitdf length(result$coef))result$coef 是模型估计出的所有系数fitdf 等于系数个数。注意 lag 必须大于 fitdf否则卡方自由度变成负数R 会直接报错。如果模型参数比较多就把 lag 提到 24。5.2 用自定义函数对比指数平滑与移动平均习题 4.8 里给了一个手写的指数平滑函数核心递推是 y[i1] a*z[i] (1-a)*y[i]。a 是平滑系数越接近 1平滑值越贴近最新观测序列越灵敏但也越容易跟着噪声走。文档用 a0.6配合简单移动平均一起看走势。实际上 R 有内置的 HoltWinters 和 ets考试上机允许用但手写一遍能帮助理解递推公式。移动平均部分要注意sma - filter(z, rep(1/12, 12), sides 1) sma - c(NA, sma)filter() 的窗口长度必须与数据周期一致月度数据用 12季度数据用 4。sides1 导致开头 11 个点全是 NA所以补一个 NA 让序列长度对齐。这里常见的坑是直接 lines(sma) 会因为 NA 位置不对而错位画出来的曲线整体右移一个点。5.3 从拟合到预测的完整校验流程最后一个上机习惯每拟合一个模型都要把“差分还原”这一步写清楚。arima() 在内部会自动完成差分还原所以 predict() 给出的是原始尺度下的预测值但如果手动做了 diff() 再用 lm() 预测就必须手工累积差分否则预测结果停留在“增量”而不是“水平”。一个可复用的校验流程是fit - arima(z, order c(2,1,0), seasonal list(order c(1,0,1), period 12)) tsdiag(fit) p - predict(fit, n.ahead 12) plot(z, xlim c(1, length(z) 12), ylim range(c(z, p$pred))) lines(p$pred, col 2) lines(p$pred 1.96 * p$se, col 2, lty 2) lines(p$pred - 1.96 * p$se, col 2, lty 2)先诊断再预测预测图必须带置信区间。上机题最后通常会问“下一年的预测值是多少”只报点预测不画区间会丢分。习题 4.7 的预测结果就是这么来的先用三次拟合外推 12 个月加回季节系数再通过 SARIMA 做二次校验两份结果对照才填表。写预测函数时把每一步的输入输出打印出来至少在 filter() 错位时能看出是哪一步开始偏差。本文还有配套的精品资源点击获取