简介：门限自回归（TAR）模型能为存在机制转换或临界效应的非线性时间序列提供灵活的分段建模方案，这份压缩包面向需要运用MATLAB完成TAR建模的研究者与数据分析学习者。包内共6个文件，以3个m脚本为核心，配合txt数据与Readme说明文件，整体仅70KB，轻量集中，便于快速运行与对照学习。目前已有357人学习下载。代码示例覆盖了TAR模型构建的主要环节，包括阈值识别、不同区间的自回归拟合、最大似然参数估计以及残差诊断；同时包含LR（似然比）图绘制逻辑，用以比较增加门限后的似然增益，辅助判断最佳门限数量、避免过拟合。配套的txt数据文件可直接用于实验演练，帮助读者从数据预处理、阈值检测到模型验证完整走通流程。对于研究经济波动、气象变化等存在状态切换的场景，这套MATLAB程序提供了清晰可复用的参考实现。

1. 门限自回归：当序列在不同状态下来回横跳，线性 AR 模型真的顶不住

jasa_03m 这类月度序列最常见的毛病，是均值、方差和自相关结构在某个临界点前后像换了个人。行情好的时候 GDP 增速的惯性很强，跌到某个阈值以下就变成另一种波动模式；用电量低峰和高峰的回归系数也不一样。普通 AR 模型只给一套全局系数，遇到这种 regime switching 只能把门槛两侧的差异平均掉，结果预测在拐点处永远慢半拍。门限自回归（Threshold Autoregressive，TAR）的思路简单粗暴：用门限变量把序列切成几段，每段各跑一个 AR，让系数自己随状态切换。这篇文章就把 TAR 从原理到落地讲透，适合手里有月度、季度观测序列、想改善拐点预测效果的数据分析或量化从业者。

2. 为什么序列一有门限效应线性回归就不够用：从模型结构说起

2.1 TAR 的数学表达：分段线性是它最核心的武器

门限自回归最早由 Tong 在 1970 年代末提出，后来和 Lim 一起整理成书，核心思想用一个分段函数就能说清。考虑一个单变量时间序列 $y_t$，如果门限变量取的是 $y_{t-d}$（也就是序列自己延迟 d 期的值），模型叫 SETAR（Self-Exciting Threshold Autoregressive），这就是最常用的一种形式。两段 SETAR 的写法是：

$$ y_t = \begin{cases} \phi_1^{(0)} + \sum_{i=1}^{p_1} \phi_1^{(i)} y_{t-i} + \varepsilon_t^{(1)}, & y_{t-d} \le c \ \phi_2^{(0)} + \sum_{i=1}^{p_2} \phi_2^{(i)} y_{t-i} + \varepsilon_t^{(2)}, & y_{t-d} > c \end{cases} $$

上面这个式子里，$c$ 是门限值，$d$ 是门限延迟阶数，$p_1$ 和 $p_2$ 是两段各自的 AR 阶数，两个残差序列 $\varepsilon_t^{(1)}$ 和 $\varepsilon_t^{(2)}$ 被假定为独立同分布的白噪声。也就是说，当前时刻落在哪个 regime，完全由 $y_{t-d}$ 和 $c$ 的大小关系决定。这个设计有个好处：门限变量是内生滞后值，预测的时候不需要额外估计未来的外生变量，递归预测非常方便。

要注意门限变量的选取不止滞后值一种。宏观实证里经常用"外生门限变量"，比如利率、油价、产出缺口，模型写成 TAR（广义叫 TVAR，Threshold Vector Autoregression，门限向量自回归）。这种写法的门限变量可以是另一个序列 $z_t$，门限条件写成 $z_{t-d} \le c$。两种做法的取舍会在 2.3 展开。至少从结构上看，TAR 不是黑匣子式的非线性模型，它本质是"多个线性模型 + 一个切换规则"，每一段都可以用 OLS 估计，解释起来跟 AR 差不多顺手。这在实际落地时是一个非常大的优势。

2.2 线性 AR 在门限场景下会翻车：一个最小例子

只看数学可能觉得 TAR 只是多了个 if-else，但它在实证上和线性 AR 的差别立竿见影。设想一个两段 SETAR(2, 1, 1) 过程，门限值 $c = 0$，两段系数分别是 $\phi_1 = 0.9$（低位 regime 强惯性）和 $\phi_2 = -0.7$（高位 regime 强反转），噪声两个 regime 里方差也不一样。把这段序列用普通 AR(1) 去拟合，OLS 会把 0.9 和 -0.7 平均成一个接近 0.1 的系数。结果是什么？模型预测的是一堆"温和的随机游走"，而真实数据是"低位持续走低、高位反复打脸"。用线性 AR 做出来的残差里会残留明显的自相关，Ljung-Box 检验大概率拒绝白噪声原假设——这就是典型的模型设定偏误。

实际做单变量序列预测时，我发现线性 AR 在拐点附近的误差方向常常是一致的：真实值向上突破时线性模型反应不足，真实值急跌时线性模型还停在原处。原因不是 AR 本身错，而是它的假设太强——自相关结构必须全周期稳定。TAR 允许每个 regime 有自己的均值和波动率，相当于把"惯性模式 A"和"波动模式 B"分开建模，拐点处切换。在 JASA_03M 这种月度数据上，如果存在明显的"扩张期—收缩期"交替，TAR 的样本外预测在转折点附近的优势会非常明显。判断方法也直观：对线性 AR 的残差做 BDS 检验或 RESET 检验，检验统计量显著就说明残差里还有非线性结构没提干净。

2.3 门限值、门限延迟和 regime 个数：三个要一起决定的参数

TAR 模型在参数上比 AR 多了一层选择困难。除了每段的 AR 阶数 $p_1$、$p_2$，还要定门限个数 k、门限值 $c$ 和延迟阶数 $d$。k 一般取 2 就够用，样本量不够大时取 3 段会导致每段样本太少，估计方差爆炸。延迟 d 表示"用多少期之前的值来决定现在的状态"，如果序列的真实前瞻时间是 3 个月，那用 $y_{t-3}$ 做门限效果最好，$d=1$ 反而会晚两期才切换。tsDyn 包里setar()函数的thDelay参数就是干这个的。

门限值 c 是 TAR 参数的"灵魂"。它跟 AR 系数不一样，不能直接用 OLS 解出来，原因在于 c 出现在分段函数的边界上——给定 c 可以用 OLS 估计两段系数，但 c 本身是一个取值连续的非凸参数。实际操作中不会直接求解析解，而是用 Chan (1993) 的网格搜索法：把候选门限值和序列的分位数一一对应，用两个 regime 的残差平方和最小作为准则去搜。搜索时一般只考虑 $y_{t-d}$ 的 15% 到 85% 分位数，留出足够的样本量给两段估计。注意这个网格搜索过程在 R 里几秒钟跑完，没必要自己优化，后面第 4 章会给出一个完整可复现的脚本。

这里得提一个常见误区：有人把门限值解释成"预测的目标水平"，比如 GDP 增速到 3% 就切换。这是不对的，门限值是门限变量分布的相对位置，它反映的是序列处于"低区间"还是"高区间"的自然分界，不是政策目标或业务目标。门限值要靠数据说话，不要先入为主设定。

3. 从数据到模型：用 tsDyn 跑通 jasa_03m 的最小完整流程

3.1 数据准备：平稳性、结构断点的前处理

在把月度序列喂给setar()之前，有几个前处理步骤必须做。首先是缺失值，TAR 模型估计时用的是普通过去值的滞后结构，中间有空洞会让延迟阶数错位，季度月度数据常用插补法填掉；如果缺失集中在尾部，宁可砍掉也不要硬填。其次是季节性：如果序列像用电量、零售额那样有明显的年度周期，不处理的话门限估计会被季节波动带走，常见的做法是做季节差分：y_t - y_{t-12}（月度数据），或者用tslm分离出季节成分后取残差。

平稳性是第三个要过的关卡。严格说 TAR 的每段 AR 都要求平稳，如果原始序列是带趋势的 I(1) 过程，两段 OLS 估计出来系数可能虚假。至少用 ADF 检验确认序列是否平稳：

library(urca) adf_test <- ur.df(jasa_03m, type = "trend", selectlags = "AIC") summary(adf_test)

说明：type="trend"允许检验回归里带趋势项，selectlags="AIC"自动选择 ADF 回归的滞后阶数。如果检验统计量大于 5% 临界值（不能拒绝单位根），先做一阶差分再继续。注意，差分后的序列做出来的是"增长率/变化量"的门限模型，解释口径要跟着变；但差分也是消除趋势手段里最稳的一种，比做线性去趋势更不容易残留伪结构。

3.2 用 setar 拟合两段门限模型：核心参数与输出解读

R 里最常用的门限自回归实现是tsDyn包的setar()函数。它内部自动完成门限网格搜索、分段 OLS 估计，输出模型系数、门限值和评价指标。我建议的最小调用方式是：

library(tsDyn) # 假设 jasa_03m 是月度 ts 对象 / 数值向量，先转成 ts y <- ts(jasa_03m, start = c(1990, 1), frequency = 12) # 自动选择延迟阶数，再拟合 SETAR(2, p1, p2) tar_fit <- setar(y, m = 3, thDelay = 1, trim = 0.15, trace = TRUE)

说明：m = 3表示两段 AR 的候选最大滞后阶数都是 3，建模时每段会按 AIC 再筛选（如果m是单个数字，默认用select.order自动定阶）。thDelay = 1是门限变量的滞后阶数 d=1，即用y_{t-1}做状态判据；如果业务经验判断"拐点提前三个月就能看出来"，改成 3。trim = 0.15是关键参数，意思是门限值只在序列分布的 15%~85% 分位数之间搜索，防止门限落在样本边缘导致某段只有几个观测值。跑完后用print(tar_fit)看结果：

print(tar_fit, digits = 5)

输出里主要关注四块内容。第一块是门限值Threshold var和Threshold value——它告诉你最优切分点在哪。第二块是两段各自的系数表，细看低位段和高位段的自回归系数结构是不是真的不同；如果两段系数非常接近，说明序列可能没有实际的门限效应。第三块是 AIC/BIC 和残差方差，用于和线性 AR 直接对比。第四块是Non-linearity test：对"两段系数相同"这个原假设做检验，p 值小于 0.05 才能说明门限效应显著。我用这类输出判断 TAR 是否值得用，通常 AIC 比线性 AR 低 2 以上 + 非线检验显著，才算有真东西。

3.3 把 fit 变成预测：递归预测与样本外评估

setar()拟合好之后，预测和评估是落地重点。predict()可以按给定的超前步数做递归预测，tsDyn还提供了样本外评估的辅助函数。下面这段代码展示标准的滚动预测评估流程：

# 滚动窗口样本外预测：固定起点，递推预测 12 个月 h <- 12 n <- length(y) pred_tar <- numeric(h) for (i in 1:h) { refit <- setar(y[1:(n - h + i - 1)], m = 3, thDelay = 1) pred_tar[i] <- predict(refit, n.ahead = 1)[1] }

这段代码的要点是每预测一个月就重估一次模型，避免把未来数据卷进参数估计造成前视偏差。predict(refit, n.ahead = 1)只推进一步，下一步会把上一步真实观测加进估计窗口——这是严格的滚动原点评估，不是递归多步预测，两者口径要区分清楚。真正做长期预测时，用predict(tar_fit, n.ahead = h)一次性迭代，注意中间路径状态切换不可观测，后几步误差会累积变大；因此样本外比较建议多用滚动 h 步方式，而不是一次跑 36 个月。

做完预测，把 TAR 的结果和线性 AR 的预测放在一张图里：横轴是时间、纵轴是实际值和两条预测线。判断标准是 RMSE 和 MAE，特别注意拐点月份的重合质量——TAR 在转折方向上的领先是它真正的价值，如果优势只来自整体 RMSE 小幅改善，可能只是过拟合。

4. 自己写网格搜索实现 TAR：不依赖黑匣子，把门限估值过程拆开看

4.1 门限网格搜索的思路：残差平方和最小化

tsDyn封装得很好，但很多从业者（包括我）会在项目早期自己写一遍门限网格搜索，原因有两个：一是 setar 的参数组合有限，二是自己想调"门限搜索粒度"和"分段最小样本量"时，封装函数不一定给你接口。门限搜索的核心是一个两层循环：外层遍历候选门限值 $c$，内层把样本切成两段分别拟合 AR，记录每段残差平方和之和（相当于似然值），最后选残差平方和最小对应的 $c$。

这套思路对应 Chan 的一致性估计方法，计算量不大，月度序列几百个观测完全能循环完。自己实现最大的好处是可以在运行过程中打印每一组门限值对应的残差变化曲线，直观看到模型对门限位置的敏感程度——如果残差曲线在很大一段区间里都平平的，那门限值的"最优"就很脆弱，换一批样本可能就完全变了。

4.2 R 代码实现：带分位数裁剪与分段 OLS 的完整脚本

下面给出一段不依赖 tsDyn 的门限网格搜索 + 分段拟合的实现，适合 jasa_03m 这类单变量月度序列：

fit_tar_grid <- function(y, p = 2, d = 1, trim = 0.15) { n <- length(y) # 构造滞后矩阵：每行是 [y_t, y_{t-1}, ..., y_{t-p}] # 门限变量取 y_{t-d}，注意 d 和 p 的关系要保证样本不重叠错乱 lag_max <- max(p, d) y_lag <- sapply(1:lag_max, function(k) c(rep(NA, k), y[1:(n - k)])) colnames(y_lag) <- paste0("L", 1:lag_max) # 去掉有缺失的前 lag_max 行 df <- data.frame(y, y_lag) df <- df[complete.cases(df), ] # 候选门限：y_{t-d} 的分位数区间 [q15, q85]，步长为 1% 分位 th_var <- df[[paste0("L", d)]] qs <- quantile(th_var, probs = seq(trim, 1 - trim, by = 0.01)) best <- NULL for (c in qs) { low_idx <- th_var <= c high_idx <- th_var > c if (sum(low_idx) < p + 2 || sum(high_idx) < p + 2) next # 两段用同一个公式结构：y_t ~ L1 + L2 + ... + Lp fmla <- as.formula(paste("y ~", paste(paste0("L", 1:p), collapse = " + "))) fit_low <- lm(fmla, data = df[low_idx, ]) fit_high <- lm(fmla, data = df[high_idx, ]) rss <- sum(residuals(fit_low)^2) + sum(residuals(fit_high)^2) # AIC 近似：残差平方和 + 2 * 参数个数惩罚 k_total <- 2 * (p + 1) + 1 # 两段截距+系数 + 门限本身 aic_val <- n * log(rss / n) + 2 * k_total if (is.null(best) || aic_val < best$aic) { best <- list(c = c, rss = rss, aic = aic_val, fit_low = fit_low, fit_high = fit_high, low_n = sum(low_idx), high_n = sum(high_idx)) } } best } res <- fit_tar_grid(jasa_03m, p = 2, d = 1, trim = 0.15) res$c summary(res$fit_low) summary(res$fit_high)

这段代码的逻辑分三层。第一层是构造滞后矩阵，用sapply生成 $L_1$ 到 $L_{\max(p,d)}$，所有回归统一用这些滞后变量，避免每次拟合都重新拼数据。第二层是候选门限的生成，Trim 默认 15%，用分位数平方根等间隔取候选值；如果你想更细，把by=0.01改成by=0.002，搜索粒度变大，但这么做对单变量几百个样本意义不大，反而可能搜到局部过拟合点。第三层是分段 OLS 和评价，k_total = 2*(p+1)+1是我常用的惩罚项公式：两段各 p+1 个参数加上 1 个门限值参数，用 AIC 辅助比较、用 RSS 做网格选择主依据。

4.3 参数调优：搜索粒度、最小分段样本量和滞后阶数的取舍

自己实现的代码让你看清三组参数的相互作用。第一是搜索粒度与样本量的关系：样本量 n=300 时，1% 分位数步长大约只给 3 个观测的间距，门限估计可能过拟合到单个点；保守做法是seq(0.15, 0.85, by=0.05)只试 15 个候选，宁可粗糙不可过头。我通常在探索阶段用 by=0.01，确认门限效应存在后，用 by=0.05 重新估，看门限值稳不稳。

最小分段样本量的设置直接影响低 regime 的估计质量。p+2是最低限，实际建议至少留 10% 的样本，也就是 300 个观测里每段至少 30 个。tsDyn::setar的trim参数默认 0.15，就是干这个。如果最优门限本身就在 0.2 分位附近，那 trim=0.15 会导致搜索边界非常接近最优值，稳健性很差——调大 trim 到 0.25 再试一次，若门限偏移很大，说明数据里门限证据并不强。

滞后阶数的选择因果链比较隐蔽。p 太低会让残差自相关没清干净，把门限效应和线性动态混在一起；p 太高则每段参数浪费，小样本分段时方差爆炸。经验做法是先用线性 AR 的 AIC 定一个基础阶数 $p_{base}$，然后尝试 $p \in {p_{base}-1, p_{base}, p_{base}+1}$，分别跑网格搜索，选择 AIC 最低且两段系数差异显著的组合。延迟 d 用同样方式扫描1:6（月度数据最多试到 6 期）。如果多个 d 给出的门限值和残差差异不大，优先选较小的 d——简单意味着更稳定。

5. TAR 建模高频踩坑清单：门限漂移、过拟合和检验失效的实战记录

5.1 门限值搜出了"最优"，换样本却大幅漂移

现象：某个月度序列第一次跑网格搜索，门限值落在 $y_{t-d}$ 的 70% 分位数，AIC 很漂亮；把样本平移半年再跑一次，门限值跳到 30% 分位，模型完全变样。

原因：门限值的识别依赖足够多的观测跨越门限两侧。如果序列长时间停留在一个 regime，跨越门限的观测只有寥寥几个，残差平方和函数在对应区间非常平，门限不可识别。

解决：画出门限搜索过程的 RSS 曲线（横轴候选门限分位数，纵轴残差平方和），如果曲线没有清晰的最低点，不要相信最小 AIC。这时改用"门限值固定为序列中位数"的方式建模，或者增大 trim 压缩搜索区间，把门限估计问题转化成"门限范围敏感性分析"，而不是选点。

5.2 非平稳序列直接建模，两段系数全是伪回归

现象：对原始水平值（明显的上升趋势）直接跑 setar，输出结果里两段 AR(1) 系数都接近 1，门限值落在趋势中点，预测长期跟着漂移，看起来"合理"但误差灾难。

原因：带趋势或单位根的序列不满足分段平稳假设，OLS 估计量收敛速度不同，门限值会被趋势主导而不是被状态切换主导。

解决：先做 ADF 或 KPSS 检验。若存在单位根，取一阶差分后重新拟合；若只是确定性趋势，用tslm(y ~ trend)去趋势后对残差建模。门限自回归允许段内均值不同，但不允许段内非平稳——这是底线。

5.3 门限效应检验显著，但残差里还有强自相关

现象：非线性检验 p 值小于 0.01，模型看起来没问题，但残差 Ljung-Box 检验仍然拒绝白噪声，预测误差呈现明显的序列相关。

原因：门限设置吸收了状态切换，但每段 AR 阶数不足。setar 两段共享同一个 p，如果低位 regime 惯性长、高位 regime 惯性短，共用阶数必然有一方欠拟合。

解决：分段分别定阶。tsDyn里可以对两段分别给mL、mH参数，自己实现的网格搜索也可以给两段不同的滞后阶数 $p_1$、$p_2$。检验门槛：每段估计完之后分别做Box.test(residuals, type="Ljung-Box")，各自通过才算合格。如果加阶数后非线性检验变得不显著，说明原来的"门限"只是线性遗漏产生的伪影——这时老老实实线性 AR 加阶数就够了。

5.4 门限变量滞后阶数选错，切换总是慢半拍

现象：样本内拟合很好，但样本外预测在拐点上系统性偏晚，真实数据已经切换到低波动 regime，预测还停留在旧 regime 里。

原因：thDelay取了 1，而真实状态切换实际上由两周甚至更早的滞后期决定。预测越远期，切换滞后被放大，误差在拐点处积累。

解决：把延迟 d 当作超参，扫描 1:6 之后对比样本外预测的拐点误差而不是整体 RMSE。拐点检测误差可以用切换月份的预测误差绝对值来衡量，比如"真实切换月 ± 2 个月窗口内 MAE"。选择 d 的另一个参考是偏自相关函数的显著阶数——门限变量往往是序列自身信息最强的那个滞后值，PACF 截尾位置能给你重要线索。

5.5 两段系数看起来差异很大，但 t 检验不显著

现象：打印系数表看到低位段 AR(1)=0.6、高位段 AR(1)=-0.3，差别明显，心里觉得模型成立；再看检验结果 p 值 0.3，根本不显著。

原因：分段后每段样本量变小，特别是 trim=0.15 时低位段只有样本的 15%，系数标准误大幅膨胀。其实"系数差异的经济显著性"和"统计显著性"是两码事。

解决：不要只看系数表，跑一个正式的线性约束检验：构造交互项模型y ~ L1*I(y_{t-d} > c) + L2*I(y_{t-d} > c) + ...，对交互项整体做 F 检验。如果 F 不显著，两个 regime 的差异可能只是噪声驱动的。另一个实用做法是看两段的残差方差是否差异明显，TAR 的一个重要应用场景就是"波动率状态切换"——即使回归系数不显著，异方差也是可以用的信号。

6. 走出"拟合好看"：门限漂移稳健性检验与一张决策图

模型拟合完不是终点，TAR 最容易被质疑的就是"你找到的门限是不是过拟合出来的"。我养成了一个固定习惯：对门限值做 1000 次自举抽样，检验门限估计的分布。做法是把fit_tar_grid包在一个 bootstrap 循环里，对残差做有放回重抽，每次生成一条新的序列重新估计门限和两段系数，最后输出门限值的 90% 置信区间。如果区间宽度超过序列四分位距的一半，那这个门限就是"面条型的"，不具备实际业务含义。

set.seed(2024) boot_c <- numeric(1000) for (b in 1:1000) { resamp <- sample(residuals(res$fit_low), size = length(y), replace = TRUE) y_boot <- fitted(res$fit_low) + fitted(res$fit_high) + resamp # 简化重抽样 boot_c[b] <- fit_tar_grid(y_boot, p = 2, d = 1)$c } quantile(boot_c, probs = c(0.05, 0.5, 0.95))

上面的代码是 bootstrap 门限分布的一个简化版，实际用的时候要把整个模型拟合过程封装成函数，保证重抽样和重估不丢样本。如果门限分布的 IQR 很宽，说明门限不稳定——这时我会降级处理：不做门限判断，改用平滑转换回归（STAR），让制度切换是渐进的而不是跳跃的。STAR 和 TAR 的区别在于转换函数是连续的 logistic 函数而不是阶跃函数，门限效应变成"随门限变量缓慢变化"，对门限位置的敏感度低得多。tsDyn里的star()可以直接调用，参数调整和setar基本同构。这个取舍是 TAR 落地里最常见的一条分叉路：样本内门限清晰选 TAR；门限模糊但转换趋势明显，选 STAR 更稳。

最后分享一个我总结的决策流程：先跑线性 AR 做基线，记录残差；对残差做非线性检验（BDS 或 RESET），不显著就用线性 AR，显著再上 TAR；TAR 拟合后做门限 bootstrap 分布检验，稳定才进入正式预测。这个流程帮我避开了至少三次"看起来有门限、实际是过拟合"的翻车。门限自回归不是银弹，但对 JASA_03M 这类存在 regimes 的月度序列，它是性价比最高的非线性建模起点——参数少、可解释、分段线性的性质让你能跟业务方说清楚"哪个状态在用什么逻辑跑"。希望这篇落地笔记能帮你在实战里少走几步弯路。

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