企业官网建设流程全解析

简介：一套基于MATLAB的门限自回归（TAR）模型实现代码，面向时间序列分析与非线性建模的研究者、学生及工程人员，解决传统线性AR模型无法刻画不同阈值区间动态特性的问题。RAR压缩包共6个文件，约70KB，含3个.m脚本（如beta、beta2、absentee）与3个.txt文本（Readme说明、演示数据等），结构精简，便于直接运行与二次修改。该资源已有356人学习下载，适合希望快速上手TAR建模流程的初学者，也适合需要参考MATLAB实现的进阶用户。整套代码覆盖模型构建的核心环节，包括阈值识别、分段自回归拟合并结合最大似然法估计参数，同时绘制LR似然比图辅助确定门限数量，搭配配套数据文件可直观复现完整分析过程，有助于深入理解非线性时间序列回归及门限选择机制，具有很高的实用参考价值。

1. 门限自回归是时间序列回归里最容易解释的非线性方案：jasa_03m 先别急着上黑匣子

如果手头这份 jasa_03m 的走势像很多经济序列一样，低区制时围绕均值反复震荡，高区制时持续偏离，普通自回归模型会把两段动态平均成一个中间值，预测结果总是慢半拍。门限自回归（Threshold Autoregressive，TAR）要做的，就是在经典时间序列回归里显式加入一个门限变量和阈值：门限变量越过某个值，截距和自回归系数整体切换。它是分段线性的，保留回归的可解释性，又比单一线性 AR 多出状态切换能力，适合在神经网络的复杂度之前先试一轮。下面用 jasa_03m 这份序列把完整路径走一遍，从选型判断到参数调优，再到最容易翻车的六个细节。

2. Threshold Model 的机制判断：jasa_03m 该用 SETAR 还是外生门限 TAR

2.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} 自身，也就是用序列自己的历史值做开关，这个模型叫 SETAR（Self-Exciting Threshold Autoregressive），这是单变量时间序列里最常见的门限自回归形式。如果 r_{t-d} 来自外生变量 z_{t-d}，比如库存水位、利率、价差，它就是更一般的 TAR，也叫外生门限自回归。拿到 jasa_03m 后先别急着写代码，要问自己两个问题：这份序列是不是只有一列观测？业务上是否存在一个比 y 自身更能解释状态切换的变量？

我一般会先画一张散点图：把线性 AR 的残差按候选门限变量排序，横轴是门限变量的取值，纵轴是残差。如果残差在某个值前后明显分叉，一侧均值抬升、另一侧方差变大，才值得做两区制。如果残差随机得没有任何结构，强行拟合门限只会得到一个人为的切点，样本外预测大概率不如线性 AR。门限自回归的卖点是“分段线性”，不是“万能非线性”，这一点从一开始就要想清楚。

2.2 门限变量的三种选法：滞后值、外生变量、滑动均值

门限变量不是随便选一列数据放进去就完事。常见做法有以下三种，适用场景差别很大：

第一轮建模我强烈建议先做 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 在 p=1 到 p=5 之间选一个起点。这一步得到的 p 是给门限模型用的初值，不是最终值。门限模型允许两个区制的滞后阶不同，比如低区制是 AR(2)、高区制是 AR(1)，但初选阶段统一用同一个 p 会让后续格点搜索简单很多。

第三道是门限非线性检验。Tsay 检验是经典做法，思路是先对样本按门限变量排序，做递归回归得到预测残差，再检验这些残差是否与回归自变量相关；Hansen 检验则是在门限候选值上做格点搜索，计算 LR 统计量再通过自助法得到近似分布。两者的共同局限是：检验显著只说明“线性被拒绝了”，不保证“门限模型一定预测得更好”。BDS 检验可以放在最后做，它检测残差里是否还有一般性非线性依赖，但无法告诉你在哪个位置存在门限。

如果三道检查都指向非线性，再做估计。如果检验不显著但业务上强烈预期存在切换，也可以做一轮门限回归，但必须在后面用滚动验证说话，不能拿样本内拟合优度当证据。

3. 手写门限自回归估计器：在 jasa_03m 上跑通 SETAR 的最小 Python 实现

3.1 先定两个必选参数：门限延迟 d 和滞后阶 p

门限模型里有两个参数比门限值本身更难选：门限延迟 d 和自回归滞后阶 p。d 表示门限变量取第 t-d 期的观测值，它决定“用多久之前的状态来切换当前回归”；p 表示自回归阶数，它决定“用多少期历史来解释当前值”。d=1 是最常用的起点，因为上一期的状态往往最能解释这一期的行为。p 的初值按 2.3 节的 AIC 或 PACF 结果取，一般先从 p=1 或 p=2 开始。

一个常见的错误是上来就做三区制。三区制 SETAR 的待估参数几乎是两区制的两倍，jasa_03m 如果只有几百个观测，每个区制里很容易只剩几十个点，系数的置信区间会宽到没有业务参考价值。我的习惯是：第一轮只做两区制，跑通数据流，再根据残差诊断决定是否需要增加区制。门限回归的建模顺序应该是“先结构后参数”，结构错了，后面调什么都是白费。

3.2 格点搜索门限值：条件最小二乘的完整实现

门限模型的参数里，β_1 和 β_2 是线性的，但门限值 c 以分段方式出现，不能直接用梯度下降求解。标准做法是条件最小二乘：固定 c，模型退化成两个独立的线性回归，用最小二乘直接算出 β_1、β_2 和残差平方和；然后遍历 c 的候选值，取残差平方和最小的那组作为估计结果。候选门限不要均匀地取在 y 的取值范围内，而是取 y 的分位数，这样才能保证门限两侧都有足够的样本，避免门限跑到数据边界上去。

下面是我在实际项目中用的一个极简实现，适合拿 jasa_03m 直接跑通流程。

import numpy as np def setar_grid(y, p=2, d=1, qs=(0.15, 0.85), n_th=20): y = np.asarray(y, dtype=float) 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)) # 单侧最少样本数：至少要能估计出 p+1 个参数，再留一点余量 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], rcond=None)[0] beta2 = np.linalg.lstsq(X[high], yt[high], rcond=None)[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) 的参数个数是 (1+p1) + (1+p2) + 1，比同等滞后阶的线性 AR 明显多。如果只看 AIC 的最小值，高区制模型很容易获胜，因为它用更多的参数去拟合了噪声。

BIC 的惩罚项是 k * log(n)，比 AIC 重得多，在样本只有几百个点时，它会强烈地偏好低区制模型。实际项目中我的经验是：把线性 AR、两区制 SETAR、三区制 SETAR 放同一张表里比较 AIC 和 BIC，如果 BIC 明确选择线性 AR，而业务上又坚信存在切换，这时候不要急着放弃门限模型，先去做滚动样本外验证。样本内信息准则解决的是“模型复杂度合不合理”的问题，解决不了“门限是否在业务时间点上真实存在”的问题。

这里还有一个隐蔽的坑：格点搜索本身是一种模型选择。当你在 20 个候选门限上挑最小 SSE 时，不管门限真实存不存在，你都有机会挑到一个让样本内 SSE 很小的切点。所以用 AIC 比较模型之前，候选门限的搜索范围和搜索密度必须事先固定，不能为了让某个模型赢而事后调整 n_th。我一般会把所有候选组合记录在一张表里：p 取 1 到 3，d 取 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, p=p, d=d)[1] def bootstrap_threshold(y, p, d, beta1, beta2, best_th, n_boot=199, seed=7): 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% 分位数之间搜索，且单侧样本量不得少于 p+2 的上界，比如 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，就认为两个模型在统计上有差异，否则只把差异描述为“信号偏弱”。

最后一步是把门限翻译成业务语言。不要写“模型在 c=2.317 处切换”，要写“上一期指标低于 2.3 时，序列会以 0.2 的速度向均值回落；高于 2.3 时，序列会以 0.9 的惯性继续上行”。门限自回归在时间序列回归里最强的优势就是这个：它能让你指着一条预测曲线说清楚现在处于哪个状态、为什么切换、切换后会朝向哪里。这个解释能力，很多时候比那一点 RMSE 的下降更值钱。

我做这个模型的一个习惯是：跑完门限估计后，先打印门限置信区间，再看两区制系数差异，最后才看预测误差。门限位置不稳的模型，系数再好看也不可信。希望这些路径和坑能帮你把 jasa_03m 这类序列的门限自回归做稳，少走一轮弯路。

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