门限自回归实战：月度序列区制切换建模与预测全解析

简介：当数据在不同区间表现出不同动态特征时，传统线性AR模型往往难以胜任，而门限自回归（TAR）模型提供了更灵活的处理方式。这份MATLAB代码包为TAR模型提供了可运行的实现示例，适合经济学、金融学及工程领域的研究者与学生按需参考。包内共6个文件，其中3个m脚本覆盖模型构建、参数估计、阈值检测与LR图绘制等核心流程，3个txt文件则提供数据与说明文档，整体压缩包仅70KB，非常轻量。读者通过运行示例，可以直观体会TAR模型将序列按阈值分段、在各区间分别拟合AR模型的思路，并借助自带数据验证最大似然估计和似然比检验的具体操作。目前已有356人学习/下载，无论是初步了解门限自回归原理，还是为实际序列建模寻找可复用的MATLAB参考，这份资源都能提供直接帮助。

1. 为什么月度序列的“拐点”让线性回归集体失灵：门限自回归的切入点

在做月度时间序列回归时，我最常遇到的不是趋势和季节性，而是“同一个模型忽然不灵了”：前两个月拟合得很好，第三个月开始残差连续同号，回头一看，序列整体换了一个运行状态。jasa_03m 就是这样一类月度序列的代表——它在线性自回归框架下表现平庸，不是因为数据乱，而是因为序列内部存在可辨识的状态门槛。门限自回归（Threshold Autoregressive Model，下文简称 TAR，也叫 Threshold model）的核心思路，就是把这个门槛找出来，将序列拆成两个（或更多）区制，在每个区制里分别做时间序列回归。这篇文章不绕理论，直接以 jasa_03m 为对象，把门限识别、参数估计、效果检验和落地避坑串成一条可复现的路径，适合正在处理销量、流量、设备状态等带阶段切换特征的月度数据的从业者。

2. 门限自回归的模型结构与选型：为什么“分两段回归”能赢

2.1 SETAR 与 Threshold model 的核心形式：门限变量从哪里来

门限自回归的基本形式可以写成：

y_t = (a1 + β11 y_{t-1} + ... + β1p y_{t-p}) · I(y_{t-d} ≤ C) + (a2 + β21 y_{t-1} + ... + β2p y_{t-p}) · I(y_{t-d} > C) + ε_t

其中 I(·) 是指示函数，条件成立取 1，否则取 0；d 是门限延迟，表示“当前时刻属于哪个区制由 d 期前的观测决定”；C 是待估的门限值。当门限变量使用序列自身的滞后项 y_{t-d} 时，这种特例叫 SETAR（Self-Exciting TAR），也是实际项目里最常用的起点。jasa_03m 这类月度序列，通常会先尝试 d=1，也就是“上一期的状态决定这一期的回归方程”，因为业务含义最直接：设备上个月过载，这个月的运行模式就变了。

门限变量也可以是时间趋势项 t、某个外生变量（比如宏观指标、气温），甚至另一个序列的滞后值。这种场景下模型严格叫 TAR 而不是 SETAR。我的经验是，除非你手里有一个强解释力的外生变量，否则先用内生滞后项做 SETAR，状态切分更好解释，也少一个变量选择的问题。等内生门限跑通后，再考虑把外生变量塞进门限函数里做对比，而不是一开始就上复杂度。

值得注意的一点是：门限自回归本质上仍然是一个线性分段模型，每个区制内回归系数是常数。它和“对全样本拟合一个多项式”的根本差异在于，TAR 的分段由可解释的门限值触发，而不是用高阶项去硬凑非线性曲线。所以它适合的序列特征很明确：存在若干个相对稳定的运行状态，状态之间切换较突然，且状态切换的方向和幅度有规律可循。

2.2 线性 AR 在什么条件下必须让位

线性 AR 模型假设整个样本期内自回归系数恒定。当序列存在区制切换时，拟合出来的系数是两套（甚至多套）真实系数的加权平均，看起来每个系数都“凑合能用”，但残差里会留下明显的结构信号——最常见的表现是残差自相关检验通不过，而增加 AR 阶数也压不下去。另一个信号是预测误差出现系统性偏置：低波动期预测尚可，高波动期预测一路滞后。

不过我也不建议一上来就把“非线性”当成救命稻草。时间序列回归领域有个常被忽略的原则：先用线性模型把趋势、季节性、异常点处理干净，再看残差。jasa_03m 这类月度序列尤其如此，因为月度数据天然带有年度周期性，如果季节性没有剔除，很可能被门限模型误判成“两个区制”，一个区制对应旺季，另一个区制对应淡季，这在业务上没意义。我会在第三章演示怎么在建模前先排除这个干扰。

判断是否该上 TAR，我的标准有三条：第一，序列或经过预处理后的残差在业务上能讲出“状态”概念，比如正常运行与故障运行、促销期与平销期；第二，样本量至少能支撑每个区制各自完成 AR 系数估计，经验值是每个区制至少要有 p+5 个观测；第三，线性 AR 残差的非线性检验给出明显信号。如果三条都不满足，门限自回归大概率只是给训练集做了一次华丽的过拟合。

2.3 模型设定的关键参数：门限个数、AR 阶数与门限延迟

落地 TAR 时实际要拍板的参数有三组：门限个数、每个区制的 AR 阶数 p、门限延迟 d。门限个数一般从 2 区制起步，也就是一个门限值。两个区制都解释不通时再试 3 区制，两个门限值分别对应“低、中、高”或“冷、正常、热”三种状态。三个区制意味着要在两维区间上做网格搜索，样本量需求陡增，月度数据如果少于 200 个点我基本不会尝试。

AR 阶数 p 可以先看偏自相关函数（PACF）截尾位置，再用 AIC/BIC 做小范围比较。月度序列的 p 通常不会太大，1～4 之间居多，因为门限本身已经消化了一部分自相关结构。门限延迟 d 的意义容易被低估，它描述的是“状态传导速度”：d=1 表示今天的区制由昨天决定，d=2 表示由前天决定。对 jasa_03m 这类数据，我不建议把 d 设成与 p 相同，二者物理意义完全不同，后面避坑章节会专门讲这个混淆点。

这三组参数的最优组合没有解析解，常见做法是把 p 和 d 放进小网格穷举，对每个组合做一次门限网格搜索，再按 AIC 或 BIC 选优。这个流程在第二章展开，但先记住一个原则：参数搜索范围宁可保守，不要贪大。p 到 4、d 到 2、门限个数固定为 1，这个范围已经能覆盖绝大多数月度业务序列的实际需求。

3. jasa_03m 的门限建模全流程：从平稳性到网格搜索代码

3.1 读取 jasa_03m 并完成平稳性预处理

拿到 jasa_03m 后，第一步不是建模，而是确认序列是否满足门限自回归的基本前提。TAR 和线性 AR 一样，要求序列是平稳的或经过变换后平稳。现实中月度销量、流量序列往往带趋势和季节性，直接套 TAR 会把确定性成分错误地归入门限切换。我一般先对原始序列做 ADF 检验，p 值大于 0.05 就先做一阶差分，差分后如果还有固定周期，再做一步季节差分。

import pandas as pd from statsmodels.tsa.stattools import adfuller # 读取 jasa_03m，假设文件里有 date 和 value 两列 df = pd.read_csv('jasa_03m.csv', parse_dates=['date'], index_col='date') y = df['value'].astype(float).dropna() # 对原始序列做 ADF 检验 res = adfuller(y, autolag='AIC') print(f'ADF = {res[0]:.3f}, p 值 = {res[1]:.3f}') # p 值大于 0.05 时做一阶差分后再检验 if res[1] > 0.05: dy = y.diff().dropna() res2 = adfuller(dy, autolag='AIC') print(f'一阶差分后 ADF = {res2[0]:.3f}, p 值 = {res2[1]:.3f}')

这段代码里的 adfuller 是 statsmodels 的官方接口，autolag='AIC' 表示由函数自动选择用于检验的滞后阶数，比手动指定更稳。dropna() 必须加，因为 diff() 会引入第一个缺失值。需要说明的是，如果检验对象是差分后序列，后面网格搜索和预测都在差分序列上进行，最终预测结果要再反向累加回原始值，否则预测的绝对值会整体偏低或偏高，这是差分建模最容易翻车的地方。

对月度数据，我还会额外看一眼季节性。一个快速判断方法是画出逐月分组箱线图，如果各月均值差异明显，就在建模前做季节差分或把季节虚拟变量放进模型。jasa_03m 如果存在明显周期而未处理，门限搜索几乎一定会捕捉到周期的高低谷，门限值落在某个固定月份区间，这种模型换个年份就失效。

3.2 用线性 AR 残差确认非线性信号

不要跳过线性基线直接跑 TAR。先拟合一个合适的线性 AR 模型，把残差留下来做非线性诊断，才能证明门限结构是数据里真实存在的东西，而不是建模者的执念。最常用的诊断是检验残差是否仍有自相关，以及残差平方序列是否出现自相关——后者是 ARCH 效应和非线性依赖的常见信号。

from statsmodels.tsa.ar_model import AutoReg from statsmodels.stats.diagnostic import acorr_ljungbox # 在差分序列 dy 上拟合 AR(2)，含截距项 model_ar = AutoReg(dy, lags=2, trend='c').fit() resid = model_ar.resid.dropna() # 残差白噪声检验 lb_resid = acorr_ljungbox(resid, lags=[6, 12], return_df=True) print('残差 Ljung-Box:\n', lb_resid) # 残差平方的白噪声检验（McLeod-Li 思路） lb_sq = acorr_ljungbox(resid ** 2, lags=[6, 12], return_df=True) print('残差平方 Ljung-Box:\n', lb_sq)

lags=[6, 12] 对月度数据是常见配置，分别对应半年和一年的记忆长度。return_df=True 让结果以 DataFrame 形式返回，方便直接看 lb_stat 和 lb_pvalue 两列。残差检验里如果 p 值小于 0.05，说明线性 AR 没吃干净信息；残差平方检验 p 值很小，则提示存在非线性依赖，门限模型就值得尝试。两者的区别要分清：前者说明还有线性自相关没提取完，后者说明信号之间可能存在状态依赖的交互结构。

如果这两个检验都干净，我不会继续往下做 TAR，而是直接改用更高阶的线性 AR 或用 ARIMA。门限自回归不是万金油，强行建模只会得到一个靠 AIC 选中、但滚动预测一塌糊涂的模型。这个“先证明非线性存在”的习惯，能帮你挡掉很多自欺欺人的建模实验。

3.3 门限网格搜索与条件最小二乘估计

确认残差里有非线性信号后，进入 TAR 建模的核心步骤：门限估计。门限值 C 没有显式解，工程上最常用的是网格搜索加条件最小二乘：把门限变量排序，在候选区间内逐点尝试，每个 C 下把样本分成两个区制，分别做一次最小二乘回归，记录残差平方和，最终选残差平方和最小的那个 C。这个搜索过程在数学上等价于对有结构突变的分段回归做最大似然估计，只是门限参数通过搜索而非梯度求解。

import numpy as np def ols_rss(X, y): """给设计矩阵 X 和目标 y，返回截距项下的残差平方和与系数""" X = np.column_stack([np.ones(len(X)), X]) beta, _, _, _ = np.linalg.lstsq(X, y, rcond=None) resid = y - X @ beta return resid @ resid, beta def tar_grid_search(y, p=1, d=1, trim=0.15): """对 y 做两区制 TAR 网格搜索，返回最优门限、最小 RSS 和系数""" y = np.asarray(y, dtype=float) n = len(y) m = max(p, d) if m >= n: raise ValueError('样本量不足以构造滞后项') # 构造回归矩阵：第 j 列是 y_{t-j} Y = y[m:] X = np.column_stack([y[m-j:n-j] for j in range(1, p+1)]) # 门限变量：y_{t-d} q = y[m-d:n-d] if d > 0 else y[m:n] # 按门限变量升序排序，便于逐点搜索 order = np.argsort(q) q_sorted, X_sorted, Y_sorted = q[order], X[order], Y[order] N = len(Y) lo = int(np.ceil(trim * N)) hi = int(np.floor((1 - trim) * N)) + 1 best_rss, best_c, best_info = np.inf, None, None for i in range(lo, hi): c = q_sorted[i] mask = q_sorted <= c # 每个区制至少要有 p 个样本，否则回归矩阵不满秩 if mask.sum() <= p or (N - mask.sum()) <= p: continue rss1, beta1 = ols_rss(X_sorted[mask], Y_sorted[mask]) rss2, beta2 = ols_rss(X_sorted[~mask], Y_sorted[~mask]) rss = rss1 + rss2 if rss < best_rss: best_rss = rss best_c = c best_info = (beta1, beta2, mask.sum(), N - mask.sum()) return best_c, best_rss, best_info

这段代码的参数要仔细解释。p 是每个区制内自回归阶数，d 是门限延迟步数，trim 是门限搜索范围的比例，默认 0.15 表示排序后最前和最后各 15% 的观测不参与门限候选。trim 的作用是防止门限被放到样本极值附近，导致某个区制只有几个观测，回归系数变成“三个点拟合一条线”的玄学。代码里 mask.sum() <= p 是第二道保险，确保两个区制的样本量都大于待估系数个数。

ols_rss 函数在 X 前拼了一列 1，表示每个区制的回归都带独立截距项。这是 TAR 的标准做法——如果强制两个区制共享截距，门限估计会被截距偏差带偏。网格搜索的复杂度是 O(N·p)，p 很小时几万条月度数据也秒出结果。如果你用的是 R，tsDyn 包也实现了类似的搜索，可以把两边的门限估计结果做交叉验证，防止语言实现细节导致结论差异。

3.4 用 AIC/BIC 在候选参数组合中定案

单个 (p, d) 组合下的门限搜索做完后，还需要在多个组合之间做选择。我的做法是写一个双层循环，外层遍历 p 和 d，内层调用 tar_grid_search，然后用 AIC 统一打分。打分公式为：

AIC = N * ln(RSS / N) + 2 * (2 * (p + 1) + 1)

其中 N 是参与回归的样本量，RSS 是网格搜索返回的最小残差平方和，2*(p+1) 是高低两个区制各 p+1 个系数，额外加 1 是给门限参数 C 的复杂度惩罚。BIC 的惩罚项则变成 (2*(p+1)+1) * ln(N)。月度序列样本量通常在 100～500 之间，AIC 和 BIC 差异不大，但 BIC 在样本量偏大时更保守，能帮你压住参数数量。

这个表是我每次建模前心里过一遍的参数清单。实际运行时，我会把每组 (p, d) 的门限值、低区制样本数、高区制样本数和 AIC 打在一张表里，直接选 AIC 最小的组合。选完还要做一件事：看两个区制的样本量比例是否悬殊。比如低区制 95 个点、高区制 8 个点，即使 AIC 最低也不建议采用，因为高区制的系数基本被几个极端点绑架，门限的置信区间会大得失去业务解释力。

3.5 门限模型预测：一步和多步

估计完参数后，预测逻辑也要按区制分流。预测时先用当前已知的门限变量（通常是 y_{t-d}）判断此刻落在哪个区制，再调用该区制的自回归方程生成下一步预测值。多步预测则把预测值不断追加到序列尾部，注意每追加一步都要重新判断区制。

def tar_predict_one(y, p, d, c, beta1, beta2): """单步预测：根据门限变量当前状态选择区制方程""" lag = y[-p:] # 最近 p 个观测，顺序为 y_{t-1}, ..., y_{t-p} q = y[-d] # 当前状态由 y_{t-d} 决定 if q <= c: return beta1[0] + np.dot(beta1[1:], lag) else: return beta2[0] + np.dot(beta2[1:], lag) def tar_predict_multi(y, p, d, c, beta1, beta2, h): """多步预测：预测值回填后继续递推""" y = list(y) preds = [] for _ in range(h): val = tar_predict_one(np.array(y), p, d, c, beta1, beta2) preds.append(val) y.append(val) return np.array(preds)

这里的 beta1[0] 是低区制截距，beta1[1:] 依次对应滞后 1～p 阶的系数。由于 tar_grid_search 里回归矩阵的第一列是 y_{t-1}，与 lag 数组的首元素对齐，所以不需要翻转 lag 顺序。真正常踩的坑是预测时用全样本的门限 C 直接判断新数据的区制——如果新数据所在的业务环境已经变化，门限值本身也应该滚动重估，这一点在后面的滚动验证章节详细说。

4. 门限效应的显著性检验：三个动作确认模型不是自嗨

4.1 线性 AR 原假设下的 SupF 统计量

网格搜索找到了残差平方和最小的门限，但残差平方和变小是必然的——多了一整套区制参数，拟合效果自然会变好。要让门限模型站得住脚，必须检验“门限效应是否显著”。经典做法是对每个候选门限计算一个 Chow 型 F 统计量，然后取所有 F 值中的最大值，这个最大值叫 SupF 统计量。它的抽样分布不是标准 F 分布，因为门限本身是搜索出来的，不是事先给定的。实务中常见做法是用 bootstrap 计算 p 值，代码里则可以先算一个粗略的 F 值做方向判断。

def rough_sup_f(y, p=1, d=1, trim=0.15): """基于最优门限计算粗略 F 统计量，用于快速判断门限效应强弱""" best_c, rss_tar, info = tar_grid_search(y, p, d, trim) beta1, beta2, n1, n2 = info m = max(p, d) Y = y[m:] X = np.column_stack([y[m-j:len(y)-j] for j in range(1, p+1)]) rss_lin, _ = ols_rss(X, Y) N = len(Y) k_lin = p + 1 k_tar = 2 * (p + 1) F = ((rss_lin - rss_tar) / (k_tar - k_lin)) / (rss_tar / (N - k_tar)) return best_c, F

这个 F 统计量只能当参考值，不能直接查 F 分布表。更严谨的做法是把 tar_grid_search 的循环改成在每个候选门限下计算 F，记录最大值，然后对原序列做 block bootstrap 重采样，每次重采样都重新搜索门限并计算 SupF，最后看原始 SupF 落在 bootstrap 分布的哪个分位上。月度数据有自相关结构，bootstrap 必须用 block 方式而不能用逐点重采样，否则会破坏序列的依赖结构。

4.2 区制系数差异的直观检验

SupF 是一场“大考”，但对业务解释来说，还有一个更直观的检验：比较两个区制的回归系数是否显著不同。可以构造一个交互模型，把门限指示变量与所有滞后项相乘，放进一个回归里，然后对交互项做联合 F 检验。如果交互项整体不显著，说明两个区制的动态结构没有本质差异，门限切分只是在拟合噪声。

这个检验的 p 值同样要谨慎解读，因为门限 C 的估计不确定性没有纳入其中，p 值会偏乐观。我的习惯是把它当“筛选器”而不是“判决书”：交互项 F 检验 p 值大于 0.05 的模型直接放弃；小于 0.01 的进入下一轮滚动验证。真正决定模型去留的，永远是样本外预测表现，而不是任何单一检验的 p 值。

4.3 RSS 曲线与门限置信区间

网格搜索过程中，每个候选门限都对应一个残差平方和，把这些点连成曲线，能直观看到门限估计的稳定性。如果曲线在最低点附近形成一个尖锐的 V 形，说明门限辨识度高；如果谷底平坦得像碗底，说明一个区间范围内的门限值效果都差不多，此时报告单个门限值会误导业务方，不如给出置信区间。

门限置信区间可以用 bootstrap 近似：对残差做 block bootstrap，多次重估门限，取 2.5% 和 97.5% 分位数。实际操作里我会在 tar_grid_search 的循环里记录所有 (c, rss) 点，然后用 matplotlib 画出曲线。看到平坦谷底时，我会主动把业务方拉来讨论：与其纠结门限是 0.31 还是 0.38，不如定义“以 0.35 为中心的门限带”，在带内给两个区制的预测结果做加权平均。这个思路和预测控制里的软切换很像，能显著提高预测稳定性。

4.4 残差白噪声与区制业务含义

最后一关是检查每个区制的残差是否还残留自相关。把两个区制的样本和残差分别取出，各自运行 3.2 里的 Ljung-Box 检验。如果某个区制残差仍有明显自相关，说明这个区制内的 AR 阶数不够，或者区制内还存在子状态。此时优先增加该区制的 p，而不是增加门限个数。门限增多会让每个区制的样本量骤减，副作用比调高 p 更严重。

业务含义的交叉验证同样不能省。模型给出低区制系数为正、高区制系数为负时，要回到业务里找证据：是否低状态有惯性、高状态有反转？如果业务上完全讲不通，哪怕统计检验全部通过，我也会怀疑是伪门限。TAR 的优点恰恰在于每个区制都对应一段可观察的历史区间，把门限切出的两段时间序列分别打上标签，直接交付给业务方核对，比给出一堆检验统计量更有说服力。

5. 门限自回归的五个高频踩坑：现象、原因和处理办法

5.1 门限值“换 trim 就变”：RSS 曲线太平坦

现象：trim 从 0.10 调到 0.15，门限从 -0.32 跳到 -0.61，两个模型 AIC 却几乎一样。原因：RSS 曲面在谷底附近存在一个平坦区间，多个门限值的拟合优度差异很小，网格搜索只是挑了其中一个点。这在月度样本量不足 150 个时尤其常见。

解决：先把 trim 固定，在 tar_grid_search 里记录所有 (c, rss) 点并画曲线。如果曲线谷底平坦，改用门限带策略——以最优门限为中心，上下各取一个标准差范围，预测时对两个区制的结果做加权平均。更规范的做法是用 block bootstrap 给出门限置信区间，报告中写“门限估计为 -0.45，95% 置信区间 [-0.61, -0.32]”，而不是写一个孤零零的点值。

5.2 把自回归阶数 p 和门限延迟 d 当成一回事

现象：代码里把 p 和 d 设成同一个值，比如 p=2、d=2，结果门限估计异常，预测效果也明显偏差。原因：p 控制的是每个区制内“用哪几期滞后做回归”，d 控制的是“用哪一期的状态判断当前区制”，两者物理意义完全不同。p=2 意味着回归项包含 y_{t-1} 和 y_{t-2}；d=2 意味着当前时刻属于哪个区制由 y_{t-2} 决定，而不是 y_{t-1}。

解决：参数搜索时把 p 和 d 放在两层独立循环里，允许 p=2、d=1 这类组合存在。选模型时不要只看 AIC，还要看门限变量的业务解释：如果 d=2 显著优于 d=1，说明状态传导有滞后，需要确认业务机制是否支持“前两期的状态决定本期行为”。

5.3 月度序列没预处理，季节性被识别成伪门限

现象：门限估计恰好落在每年某个月份对应的数值附近，两个区制分别对应“旺季”和“淡季”，业务解释听起来很合理，但换一年数据就失效。原因：月度数据自带的年度周期性没有剔除，门限模型把确定性周期误判为状态切换。线性 AR 有同样问题，但门限模型的区制结构会让这种误判更隐蔽。

解决：建模前先做季节差分或季节调整。对 jasa_03m，我一般先跑 STL 分解，把趋势和季节分量剥掉，只用余项做门限建模。如果业务上必须保留季节性，就把季节虚拟变量直接放进每个区制的回归方程里，而不是让门限去承担周期切分的任务。

5.4 区制样本量悬殊，少数几个点绑架了高区制估计

现象：低区制样本 180 个，高区制只有 12 个，但 AIC 显示这个模型最优。原因：高区制的 12 个点拟合残差极小，显著拉低了整体 RSS，但这 12 个点的系数估计方差极大，换个样本区间就完全变样。trim 参数设得过小是直接原因。

解决：强制检查 tar_grid_search 返回的 n1 和 n2，任一区制样本数少于 p+5 就丢弃该组合。同时把 trim 下限设到 0.15。还有一种思路是让两个区制使用不同的 AR 阶数：样本多的区制用 p=2，样本少的区制用 p=1，减少待估参数量。这个功能需要修改网格搜索的代码结构，但收益非常明显。

5.5 预测时区制来回横跳，门限附近预测值抖动

现象：预测点正好落在门限附近，前一步按低区制预测，后一步按高区制预测，预测序列出现锯齿状跳变。原因：门限 C 本质上是一个点估计，真实状态切换在 C 附近存在模糊带，但模型做了硬切分。单步预测对门限位置高度敏感。

解决：引入滞回带（hysteresis）机制——进入高区制需要门限变量连续 k 期大于 C+δ，退出高区制需要连续 k 期小于 C-δ。δ 可以取门限置信区间宽度的一半，k 取 2 或 3。这样预测序列在门限附近不会频繁切换状态，代价是对真实突变点的响应会滞后一到两期，但整体预测稳定性大幅提升。这个技巧在工业设备预警场景里非常实用。

6. 用滚动时间窗验证门限模型的真实增益：一个预测步长的对比实验

门限模型在训练集上几乎总是比线性 AR 好看，所以我的最后一个动作永远是滚动时间窗对比验证。做法很朴素：固定一个窗口长度（比如 120 个月），从第 120 期开始逐期滚动，每一期只用窗口内的数据重新估计两个模型，然后各自预测下一期，对比预测误差。这样既检验了模型的样本外能力，也模拟了真实业务中“模型定期重训”的落地姿势。

from statsmodels.tsa.ar_model import AutoReg def rolling_compare(y, p, d, window=120, horizon=1): """滚动对比线性 AR 和 TAR 的一步预测 MSE""" y = np.asarray(y, dtype=float) n = len(y) ar_errors, tar_errors = [], [] for end in range(window, n - horizon + 1): train = y[:end] test = y[end:end + horizon] # 线性 AR 基线 ar = AutoReg(train, lags=p, trend='c').fit() ar_pred = ar.forecast(horizon)[-1] # TAR：滚动窗口内重新估计门限 try: c, _, info = tar_grid_search(train, p, d, trim=0.15) beta1, beta2, _, _ = info tar_pred = tar_predict_one(train, p, d, c, beta1, beta2) except Exception: continue ar_errors.append((test[0] - ar_pred) ** 2) tar_errors.append((test[0] - tar_pred) ** 2) print('AR MSE:', np.mean(ar_errors)) print('TAR MSE:', np.mean(tar_errors)) return np.mean(ar_errors), np.mean(tar_errors)

这段代码里 TAR 的门限 C 在每个窗口都重新估计，这是刻意为之。有些团队为了省事，用全样本训一次门限，然后对后续所有窗口固定使用，这在序列结构稳定时问题不大，但一旦业务状态漂移，固定门限会让预测迅速失真。滚动窗口内的重估成本不高，因为 tar_grid_search 的时间复杂度很低，逐期滚动对月度数据完全吃得消。

对比结果要留意一个常见现象：整体 MSE 差异不大，但分时段看差异明显。我会把 AR 和 TAR 的逐期绝对误差画在同一张图里，标出误差差值最大的时间段，再回到业务日历里找对应事件。如果 TAR 的优势集中在某个特定营业状态切换期，说明门限模型的价值是“在关键时刻不出大错”，这种价值用平均 MSE 看不出来，但业务部门往往最在意。如果滚动验证里 TAR 的 MSE 只是偶尔低于 AR，我会选择更保守的线性 AR，因为简单的模型在维护成本和异常态表现上都更可控。

我现在的习惯是，把滚动验证纳入每个时间序列项目的固定交付物。门限模型最大的价值不是证明“非线性存在”，而是让预测在状态切换时少交学费。真正经过滚动验证的东西，哪怕只是比 AR 好一点点，也值得留在线上。希望帮到你。

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