门限自回归（TAR/SETAR）建模全流程与避坑指南

简介一套基于MATLAB的门限自回归TAR模型实现代码面向时间序列分析与非线性建模的研究者、学生及工程人员解决传统线性AR模型无法刻画不同阈值区间动态特性的问题。RAR压缩包共6个文件约70KB含3个.m脚本如beta、beta2、absentee与3个.txt文本Readme说明、演示数据等结构精简便于直接运行与二次修改。该资源已有356人学习下载适合希望快速上手TAR建模流程的初学者也适合需要参考MATLAB实现的进阶用户。整套代码覆盖模型构建的核心环节包括阈值识别、分段自回归拟合并结合最大似然法估计参数同时绘制LR似然比图辅助确定门限数量搭配配套数据文件可直观复现完整分析过程有助于深入理解非线性时间序列回归及门限选择机制具有很高的实用参考价值。1. 门限自回归是时间序列回归里最容易解释的非线性方案jasa_03m 先别急着上黑匣子如果手头这份 jasa_03m 的走势像很多经济序列一样低区制时围绕均值反复震荡高区制时持续偏离普通自回归模型会把两段动态平均成一个中间值预测结果总是慢半拍。门限自回归Threshold AutoregressiveTAR要做的就是在经典时间序列回归里显式加入一个门限变量和阈值门限变量越过某个值截距和自回归系数整体切换。它是分段线性的保留回归的可解释性又比单一线性 AR 多出状态切换能力适合在神经网络的复杂度之前先试一轮。下面用 jasa_03m 这份序列把完整路径走一遍从选型判断到参数调优再到最容易翻车的六个细节。2. Threshold Model 的机制判断jasa_03m 该用 SETAR 还是外生门限 TAR2.1 门限自回归与线性 AR 的本质区别系数不是连续变而是跳变线性 AR 模型假设 y_t 的动态结构在样本期内保持不变写成 y_t α φ_1 y_{t-1} … φ_p y_{t-p} ε_t。门限自回归把这个假设放宽成两段当门限变量 r_{t-d} ≤ c 时回归系数是 (α_1, φ_{1,1}, …, φ_{1,p})当 r_{t-d} c 时回归系数变为 (α_2, φ_{2,1}, …, φ_{2,p})。注意这不是两个独立的模型因为门限变量是可观测的每个样本点在建模时就能被明确归入某一个区制而不是像马尔科夫转换模型那样依赖不可观测的状态概率。门限变量的选择决定了模型的具体名字。如果 r_{t-d} 取的是 y_{t-d} 自身也就是用序列自己的历史值做开关这个模型叫 SETARSelf-Exciting Threshold Autoregressive这是单变量时间序列里最常见的门限自回归形式。如果 r_{t-d} 来自外生变量 z_{t-d}比如库存水位、利率、价差它就是更一般的 TAR也叫外生门限自回归。拿到 jasa_03m 后先别急着写代码要问自己两个问题这份序列是不是只有一列观测业务上是否存在一个比 y 自身更能解释状态切换的变量我一般会先画一张散点图把线性 AR 的残差按候选门限变量排序横轴是门限变量的取值纵轴是残差。如果残差在某个值前后明显分叉一侧均值抬升、另一侧方差变大才值得做两区制。如果残差随机得没有任何结构强行拟合门限只会得到一个人为的切点样本外预测大概率不如线性 AR。门限自回归的卖点是“分段线性”不是“万能非线性”这一点从一开始就要想清楚。2.2 门限变量的三种选法滞后值、外生变量、滑动均值门限变量不是随便选一列数据放进去就完事。常见做法有以下三种适用场景差别很大门限变量形式模型称呼适合场景自身滞后值 y_{t-d}SETAR只有单变量状态切换由上一期的走势决定外生可观测变量 z_{t-d}TAR有明确的业务开关变量如库存阈值、政策门槛自身滑动均值 (y_t…y_{t-L1})/L平滑门限 TAR单点噪声太大怕门限被个别异常值触发第一轮建模我强烈建议先做 SETAR也就是用 y_{t-d} 做门限变量。原因很实际它不引入额外数据也不增加新的待估窗口参数 L。滑动均值会多出一个平滑长度 L等于给模型再加了一层调参负担如果 jasa_03m 的样本量不大L 的变化很可能比门限 c 的变化更敏感最后你根本分不清是门限效应在起作用还是平滑参数在起作用。外生门限 TAR 的诱惑在于解释性更强你可以在报告里直接写“当库存低于 1200 万时进入低增长区制高于 1200 万时进入高增长区制”。但前提是 z 和 y 之间不能存在明显的同期反馈关系否则门限变量本身也是内生的估计出来的 c 会偏。先做 SETAR 还有一个好处门限变量的滞后阶 d 和自回归滞后阶 p 在同一套参数体系里后面切到外生 TAR 时思路完全一致只是把 r 的取值来源换掉。2.3 非线性检验与参数初选ADF、Tsay、BDS 怎么配合在正式估计门限模型之前至少要过三道检查。第一道是平稳性。门限自回归的区制内通常要求近似平稳如果 jasa_03m 的水平值明显带趋势先做差分或者取对数收益别把趋势项留给门限去切。对月度序列还要看季节项有固定季节峰值的序列最好先做季节调整否则门限很容易被季节效应劫持这一点在第 5 章会专门讲。第二道是线性 AR 定阶。用 PACF 看截尾阶数或者用 AIC 在 p1 到 p5 之间选一个起点。这一步得到的 p 是给门限模型用的初值不是最终值。门限模型允许两个区制的滞后阶不同比如低区制是 AR(2)、高区制是 AR(1)但初选阶段统一用同一个 p 会让后续格点搜索简单很多。第三道是门限非线性检验。Tsay 检验是经典做法思路是先对样本按门限变量排序做递归回归得到预测残差再检验这些残差是否与回归自变量相关Hansen 检验则是在门限候选值上做格点搜索计算 LR 统计量再通过自助法得到近似分布。两者的共同局限是检验显著只说明“线性被拒绝了”不保证“门限模型一定预测得更好”。BDS 检验可以放在最后做它检测残差里是否还有一般性非线性依赖但无法告诉你在哪个位置存在门限。参数初选方法实操建议pPACF 截尾 / AIC先统一 p建模后再让两区制分别定阶d对 1 到 p 分别跑 Tsay 检验选最小 p 值对应的 d同时看区制样本量c15% 到 85% 分位数格点搜索禁止在序列边界附近找门限区制数AIC / BIC 比较两区制优先三区制只在样本充足时试如果三道检查都指向非线性再做估计。如果检验不显著但业务上强烈预期存在切换也可以做一轮门限回归但必须在后面用滚动验证说话不能拿样本内拟合优度当证据。3. 手写门限自回归估计器在 jasa_03m 上跑通 SETAR 的最小 Python 实现3.1 先定两个必选参数门限延迟 d 和滞后阶 p门限模型里有两个参数比门限值本身更难选门限延迟 d 和自回归滞后阶 p。d 表示门限变量取第 t-d 期的观测值它决定“用多久之前的状态来切换当前回归”p 表示自回归阶数它决定“用多少期历史来解释当前值”。d1 是最常用的起点因为上一期的状态往往最能解释这一期的行为。p 的初值按 2.3 节的 AIC 或 PACF 结果取一般先从 p1 或 p2 开始。一个常见的错误是上来就做三区制。三区制 SETAR 的待估参数几乎是两区制的两倍jasa_03m 如果只有几百个观测每个区制里很容易只剩几十个点系数的置信区间会宽到没有业务参考价值。我的习惯是第一轮只做两区制跑通数据流再根据残差诊断决定是否需要增加区制。门限回归的建模顺序应该是“先结构后参数”结构错了后面调什么都是白费。3.2 格点搜索门限值条件最小二乘的完整实现门限模型的参数里β_1 和 β_2 是线性的但门限值 c 以分段方式出现不能直接用梯度下降求解。标准做法是条件最小二乘固定 c模型退化成两个独立的线性回归用最小二乘直接算出 β_1、β_2 和残差平方和然后遍历 c 的候选值取残差平方和最小的那组作为估计结果。候选门限不要均匀地取在 y 的取值范围内而是取 y 的分位数这样才能保证门限两侧都有足够的样本避免门限跑到数据边界上去。下面是我在实际项目中用的一个极简实现适合拿 jasa_03m 直接跑通流程。import numpy as np def setar_grid(y, p2, d1, qs(0.15, 0.85), n_th20): y np.asarray(y, dtypefloat) n len(y) # 建模样本从第 p 期开始因为要用前 p 期做滞后特征 idx np.arange(p, n) # 门限变量取滞后 d 期的 y 值 r y[idx - d] # 设计矩阵常数项 p 列滞后值 X np.column_stack( [np.ones_like(y[idx])] [y[idx - i] for i in range(1, p 1)] ) yt y[idx] # 候选门限取分位数避免在样本边界处搜索 th_candidates np.quantile(r, np.linspace(qs[0], qs[1], n_th)) # 单侧最少样本数至少要能估计出 p1 个参数再留一点余量 min_cnt max(p 2, 10) best None for th in th_candidates: low r th high ~low # 某一侧样本太少时直接跳过不参与比较 if low.sum() min_cnt or high.sum() min_cnt: continue # 固定门限后两个区制分别做 OLS 求解系数 beta1 np.linalg.lstsq(X[low], yt[low], rcondNone)[0] beta2 np.linalg.lstsq(X[high], yt[high], rcondNone)[0] # 总残差平方和两段各自的 SSE 相加 sse np.sum((yt[low] - X[low] beta1) ** 2) \ np.sum((yt[high] - X[high] beta2) ** 2) if best is None or sse best[0]: best (sse, th, beta1, beta2) return best代码的逻辑分成四步构造滞后特征生成门限变量固定候选门限做两段 OLS最后挑选最小 SSE 对应的门限。参数里的 p 是自回归阶数d 是门限延迟qs 规定了候选门限搜索的采样区间n_th 是候选门限个数。n_th 不宜太大20 到 50 足够太大了会让模型在样本噪声上找到一个偶然达到最小 SSE 的门限样本外表现反而更差。返回结果里 best[0] 是 SSE、best[1] 是门限值、best[2] 和 best[3] 是两段回归系数。3.3 用 R 的 tsDyn 封装做同款估计什么时候不需要自己写自己手写这段代码的好处是每个参数都暴露在眼前适合做第 4 章的 bootstrap 置信区间也适合排查为什么门限值会跑到边界上。坏处是它默认两区制共用同一个 p 阶滞后且不支持外生门限变量。如果 jasa_03m 的分析场景更复杂需要每侧独立定阶、需要外生门限变量、需要自动比较多区制模型直接换 R 的 tsDyn 包会更省事。tsDyn 里做 SETAR 的入口函数是 setar参数名里 m 对应滞后阶d 对应门限延迟thDelay 对应门限变量的额外延迟具体参数名建议以你自己环境里 ?setar 的帮助文档为准因为不同版本的默认值有过调整。我的习惯是先用手写版本跑出初始门限和系数再换封装包做正式估计。这样有两个好处一是能发现数据里明显的样本量问题候选门限搜索如果跳过了大量点说明某个区制根本撑不住二是换到 tsDyn 之后看到它给出的门限值和我手写结果一致我对结果的信任度才会从“模型跑通”提升到“模型可信”。4. 门限值、滞后阶与置信区间门限自回归调优的三件套4.1 AIC/BIC 选择区制数与滞后阶的注意事项门限模型的 AIC 比线性 AR 的 AIC 更容易被误解。先回忆一下 AIC 的表达AIC n * log(σ²) 2k其中 k 是参数个数。在门限模型里k 不是简单的一列系数而是两个区制各自截距和滞后系数之和再加上一个门限值 c。换句话说两区制 SETAR(2, p1, p2) 的参数个数是 (1p1) (1p2) 1比同等滞后阶的线性 AR 明显多。如果只看 AIC 的最小值高区制模型很容易获胜因为它用更多的参数去拟合了噪声。BIC 的惩罚项是 k * log(n)比 AIC 重得多在样本只有几百个点时它会强烈地偏好低区制模型。实际项目中我的经验是把线性 AR、两区制 SETAR、三区制 SETAR 放同一张表里比较 AIC 和 BIC如果 BIC 明确选择线性 AR而业务上又坚信存在切换这时候不要急着放弃门限模型先去做滚动样本外验证。样本内信息准则解决的是“模型复杂度合不合理”的问题解决不了“门限是否在业务时间点上真实存在”的问题。这里还有一个隐蔽的坑格点搜索本身是一种模型选择。当你在 20 个候选门限上挑最小 SSE 时不管门限真实存不存在你都有机会挑到一个让样本内 SSE 很小的切点。所以用 AIC 比较模型之前候选门限的搜索范围和搜索密度必须事先固定不能为了让某个模型赢而事后调整 n_th。我一般会把所有候选组合记录在一张表里p 取 1 到 3d 取 1 到 p区制数取 1 和 2然后统一用分位数 0.15 到 0.85 做候选门限。这样后续任何人复现结果都可比。4.2 用残差 bootstrap 给门限值画置信区间门限值 c 的分布不是标准正态分布直接用回归系数那套标准误去推根本不成立。常见做法是用残差 bootstrap 给 c 一个置信区间代码在上一节 setar_grid 的基础上可以这样写。def reestimate_threshold(yb, p, d): # 对 bootstrap 样本重新搜索门限返回门限值 return setar_grid(yb, pp, dd)[1] def bootstrap_threshold(y, p, d, beta1, beta2, best_th, n_boot199, seed7): rng np.random.default_rng(seed) n len(y) # 用原始拟合残差作为误差池 resid [] for t in range(p, n): Xt np.r_[1, [y[t - i] for i in range(1, p 1)]] mu Xt (beta1 if y[t - d] best_th else beta2) resid.append(y[t] - mu) resid np.asarray(resid) boot_ths [] for _ in range(n_boot): # 以原始前 p 期为初始值用残差重抽样递推生成新序列 yb np.empty(n) yb[:p] y[:p] for t in range(p, n): Xt np.r_[1, [yb[t - i] for i in range(1, p 1)]] e rng.choice(resid) if yb[t - d] best_th: yb[t] Xt beta1 e else: yb[t] Xt beta2 e # 每条 bootstrap 序列都重新搜索门限 th_b reestimate_threshold(yb, p, d) if th_b is not None: boot_ths.append(th_b) # 95% 分位数区间 return np.quantile(boot_ths, [0.025, 0.975])这段代码是固定滞后阶、固定门限参数下的一种近似残差 bootstrap它假设残差是近似不相关的先重抽样残差再按门限回归的递推关系生成新序列然后重新搜索门限。需要注意两点。第一如果残差存在明显的自相关或条件异方差这个近似会偏乐观更稳妥的是改用 block bootstrap 或 sieve bootstrap。第二如果门限置信区间跨度和序列本身的数量级差不多甚至把 10% 分位数和 90% 分位数都包进去了说明门限位置极不稳定这时候报告一个“平均门限”没有意义应该回到数据检查是否有异常值或缺失段。4.3 输出可解释结果两区制系数和门限的业务含义门限模型最终交付的不是一个预测值而是一个可解释的状态机制。拿到两区制系数后先分别计算两个区制的隐含长期均值公式是 α_i / (1 - Σφ_{i,j})它代表门限两侧序列各自趋向的水平。如果低区制长期均值是 1.2高区制长期均值是 4.8门限 c 是 2.3业务表达就是当上一期观测低于 2.3 时序列倾向于回落到 1.2 附近一旦越过 2.3序列会转向 4.8 附近的新水平。还要看两个区制的自回归系数和。低区制系数和接近 0.2说明进入低区制后衰减很快高区制系数和接近 0.9说明越过门限后序列有很强的惯性。这种差异才是门限模型最值钱的信息。两区制系数如果几乎没有差别门限只是切了一下截距那本质上更接近分段线性回归而不是门限自回归预测收益会非常有限。5. 排查与避坑门限自回归建模中六个容易翻车的细节下面六条是我在 jasa_03m 这类序列上反复踩过的坑按出现频率排序。5.1 门限值被拽到序列边界现象最优门限落在序列 5% 分位数以下或 95% 分位数以上另一区制只有零星几个样本。原因当数据里其实只有一个区制时SSE 最小化会把几个离群点切出去人为制造一个“高区制”来吸收异常值。解决候选门限必须在 15% 到 85% 分位数之间搜索且单侧样本量不得少于 p2 的上界比如 10 个。如果按这个约束搜索后最优门限还是贴着边界说明数据不支持两区制别硬凑。5.2 门限非线性检验显著但门限模型预测更差现象Tsay 检验 p 值小于 0.01滚动验证里 SETAR 的 RMSE 反而比线性 AR 高。原因检验拒绝线性只说明存在某种非线性结构不一定是硬门限也可能是平滑转换STAR或条件异方差导致的伪拒绝。解决把门限模型和 STAR 模型都跑一遍比较滚动预测误差如果两个模型都打不过线性 AR就老实承认这门限在这个数据集上没有实用价值。5.3 门限附近预测值来回跳变现象真实值在门限附近波动时预测结果一会用低区制系数、一会用高区制系数预测序列出现明显锯齿。原因门限变量是滞后的当它恰好卡在 c 附近时微小变化就会触发区制切换。解决报告门限的置信区间在区间内的预测不要只取一个区制的结果而是对两个区制预测做加权平均权重按门限变量与 c 的距离衰减。这个“犹豫区”在业务上比单一门限点更真实。5.4 季节成分劫持门限现象门限值总落在某个月份或某个固定时间点附近区制切换看起来像日历效应。原因月度序列没做季节调整季节波动比周期切换更强门限模型把“每年固定月份的变化”解释成了状态切换。解决先做季节调整或用季节哑变量参与回归门限变量一定用季节调整后的循环成分。jasa_03m 如果带月份标签这一步不能省。5.5 残差自相关导致门限检验虚高现象Tsay 检验显著但把 AR 阶数从 1 加到 3 后检验又变不显著了。原因门限检验的统计量通常假定残差是 iid 的线性 AR 定阶不足时残差里的自相关会被门限结构吸收造成伪非线性。解决在跑门限检验之前先对线性 AR 残差做 Ljung-Box 白噪声检验和 ARCH 检验。残差不白就先加滞后阶或换模型处理方差再回来看门限。5.6 用样本内 SSE 对比模型得出过拟合幻觉现象SETAR 的样本内 SSE 比 AR 低 20%滚动预测却落后。原因格点搜索在 20 到 50 个候选门限里挑最小值本质上做了一次隐式模型选择样本内 SSE 一定会偏乐观。解决门限候选搜索方案事先固定模型比较一律用滚动样本外误差。有一个简单办法把 jasa_03m 前一半做训练、后一半做测试训练段选出最优门限后测试段不再重估门限只按训练时的 c 做预测。这样能筛掉大量靠噪声拟合出来的门限。6. 验证门限自回归的进阶做法滚动胜率、DM 检验与门限的业务化门限模型好不好不能只看一两次预测的误差。我习惯的验证方式是把两个模型放在同一条滚动跑道上窗口从序列中间开始每一步用当前窗口重估门限和系数预测下一期然后窗口向后推一期。滚动结束后记录每一期的 AR 误差和 SETAR 误差先看总体 RMSE再看“SETAR 胜出的期数占比”。占比比平均误差更能说明问题如果两模型只在个别区间有差异平均 RMSE 会把这种差异稀释掉。两模型误差之间的差异是否显著可以用 Diebold-Mariano 检验来回答。它检验的原假设是两个模型的预测误差相等统计量由误差差值的均值除以长期方差得到对一步预测直接计算差值的均值与标准差就能做一个近似判断。实际操作中我会把滚动误差差值的均值除以标准差如果绝对值大于 1.96就认为两个模型在统计上有差异否则只把差异描述为“信号偏弱”。最后一步是把门限翻译成业务语言。不要写“模型在 c2.317 处切换”要写“上一期指标低于 2.3 时序列会以 0.2 的速度向均值回落高于 2.3 时序列会以 0.9 的惯性继续上行”。门限自回归在时间序列回归里最强的优势就是这个它能让你指着一条预测曲线说清楚现在处于哪个状态、为什么切换、切换后会朝向哪里。这个解释能力很多时候比那一点 RMSE 的下降更值钱。我做这个模型的一个习惯是跑完门限估计后先打印门限置信区间再看两区制系数差异最后才看预测误差。门限位置不稳的模型系数再好看也不可信。希望这些路径和坑能帮你把 jasa_03m 这类序列的门限自回归做稳少走一轮弯路。本文还有配套的精品资源点击获取