简介：面向时间序列分析与计量经济学研究者，这份资源提供基于MATLAB的门限自回归（TAR）模型实现示例，旨在解决数据存在阈值效应时线性AR模型拟合不足的问题。压缩包共6个文件，3个m脚本分别承担阈值检测、分段建模与参数估计、似然比（LR）图绘制等任务，3个txt文件则提供详细说明和示例数据，整套代码约70KB，轻量精悍。目前已有356人学习下载，适合课程教学、论文复现或入门非线性时间序列建模。通过该示例可完整走通TAR建模流程：先进行数据预处理与异常值排查，再运用Quandt-Andrews或CUSUM检验初筛候选阈值，随后按不同区间拟合独立AR模型并利用最大似然法求参，最后通过LR图对比不同门限数量的拟合增益，辅助判断是否值得增加门限以平衡复杂度与拟合度；同时附带残差独立性、均值与方差稳定性等诊断思路，方便读者迁移到其他非线性时序数据。
1. 门限自回归到底在解决什么问题：从 jasa_03m 这类会切换机制的序列说起
拿到 jasa_03m 这份月度数据时，我的第一反应是直接用 ARIMA 硬拟合——大部分单变量时间序列回归，用 ARIMA 收尾都够用。但残差图在几个区间里出现成段的同号偏离，说明线性假设在局部失效。门限自回归（Threshold Autoregressive Model，TAR）就是为这种机制切换设计的：用一个滞后观测值做门限变量，把样本切成两个 regime，各自跑一套 AR。它解决的问题很具体——序列在某个触发点前后呈现截然不同的行为，普通时间序列回归会把这种切换平均掉，TAR 则把它当作模型主体。适合水文流量、用电负荷、宏观指标这类门槛效应明显的数据。这篇笔记从建模逻辑讲起，落到 R 实现、参数与踩坑。
2. 门限自回归的建模逻辑：TAR 是分段的时间序列回归，不是 AR 的补丁
2.1 一个门限变量把样本切成两个 regime
TAR 最朴素的形式是两 regime 的 AR(p)。模型写出来是这样：
y_t = β10 + β11 y_{t-1} + … + β1p y_{t-p} + ε_t，当 y_{t-d} ≤ r
y_t = β20 + β21 y_{t-1} + … + β2p y_{t-p} + ε_t，当 y_{t-d} > r
两个 regime 的滞后阶数可以不同，回归变量甚至也可以不同，但最常用的版本是「同阶数、同回归变量、只切截距与系数」。这里的关键是门限变量取的是序列自身的滞后值 y_{t-d}，而不是时间 t，也不是外部变量。这就是为什么它叫 threshold model 而不叫分段线性回归——切分依据是序列状态，不是时间墙。jasa_03m 这类数据里最常见的情形是：低值区平稳、高值区波动放大，或者反过来，而恰恰在中间某个触发点上行为突变。
两个 regime 各自都是线性回归，所以参数估计在给定 r 之后就是两个独立的 OLS 问题。麻烦的是 r 本身也是未知的，它不能用 OLS 直接解出来，得靠搜索。这个「先搜门限、再估系数」的过程，决定了 TAR 的参数估计和普通时间序列回归有本质区别——普通 AR 的估计是一个凸优化，TAR 的估计是一个非光滑的网格搜索。
2.2 三个要素：滞后阶数 p、延迟参数 d、门限值 r
建模 TAR 要固定三个东西。p 是每个 regime 内部的 AR 阶数，它决定模型的动态记忆，通常用 PACF 截尾或者 AIC 选。d 是延迟参数，决定用哪一期的滞后观测做触发变量。d=1 表示「上一期的状态决定本期落在哪个 regime」，这是最常见的设定；d 变大意味着切换存在滞后，比如金融波动率里常见 d=2、d=3 的延迟响应。r 是门限值，决定样本如何被划分。
Chan (1993) 给出的估计思路是：在门限变量的取值范围内逐个候选点搜索，每个候选点把样本切成两段，分别做 OLS，取两段残差平方和最小的那个点作为 r 的估计。这个做法看着简单，实际是 TAR 建模里最容易被忽略自由度的地方——因为 r 一旦参与搜索，模型的复杂性就不是「两个 regime 各一套系数」这么简单了。
三者的辨识度是不一样的。p 和 d 是离散的，可以用信息准则去选；r 是连续的，但搜索结果对样本量极其敏感。实操上我会把 d 固定在一个业务上说得通的值（比如 1 或与差分阶数对齐），再用搜索定 r，最后回头微调 p，而不是三个一起满世界搜。还有一个容易被忽略的细节：如果门限变量滞后太多，切换本身会呈现周期性——比如月度数据里 d=12，门限触发几乎等于「上一年的同月状态」，模型会退化成一个季节性模型，那就不该叫门限效应了。
2.3 和马尔可夫切换、线性 AR、分解类方法的选型边界
TAR 不是唯一能处理机制切换的模型。选型时常见的是这几类：
| 模型 | 切换机制 | 门限/状态可否解释 | 典型场景 |
|---|---|---|---|
| 线性 AR/ARIMA | 无切换 | 不涉及 | 平稳无突变序列 |
| TAR | 由 y_{t-d} 与 r 的比较触发 | 可以，r 有业务含义 | 水位、负荷、宏观指标分段 |
| 马尔可夫切换 MS-AR | 隐状态按转移概率切换 | 难解释，只给概率 | 金融波动、景气周期 |
| STL/小波分解 | 无切换，按频率/趋势分解 | 不涉及 | 强季节性与趋势叠加 |
判断依据很朴素：先画 y_t 对 y_{t-d} 的散点，看是否呈现明显的分段线性结构。如果有两组点云且分界清晰，TAR 值得做；如果只是波动率忽大忽小但均值连续，MS-AR 更合适；如果只是季节性和趋势的叠加，先 STL 分解再做线性 AR 就够，没必要引入门限。选 TAR 而不是其他模型，最终要回答的问题是：你的门限值能不能被业务解释。水文上的「警戒水位」、用电上的「峰谷切换点」、宏观上的「阈值政策生效点」，都是 r 可以直接对应的东西。如果搜出来的 r 在业务上说不通，那大概率是模型选错了，而不是参数没调好。模型定完后的第一件事都是看残差是否白噪声，TAR 的残差诊断和线性 AR 完全一样，Ljung-Box、QQ 图、异方差检验都可以直接用；残差还有明显结构，说明门限没把机制切干净，或者切点选得太少。
3. 在 R 里跑通 jasa_03m 的门限自回归：数据清洗到 setar 拟合
3.1 数据读入、序列化与平稳性检查
jasa_03m 这类命名常见于按月采样的单变量数据，03m 一般指 03 号月度序列或月度粒度。读进来之后第一件事不是建模，是确认序列的形态和统计性质。我一般用 read.csv 读入，再用 ts() 声明频率，方便后面画图和做季节诊断。
## 读入数据：文件里至少要有时间列和取值列 jasa <- read.csv("jasa_03m.csv", header = TRUE) y <- ts(jasa$value, start = c(2004, 1), frequency = 12) ## 看一眼形态和自相关 plot(y) acf(y, lag.max = 24) ## 平稳性：ADF 原假设是非平稳，p 值不显著就先差分 library(tseries) adf.test(y) dy <- diff(y) adf.test(dy)plot 和 acf 是帮你定「要不要预处理」的。acf 在 lag 1 上拖尾衰减，大概率有单位根；lag 12 和 24 上明显突出，说明有季节性。最常见的处理是取一阶差分 dy 把趋势拿掉；如果季节性占主导，可以先做 12 阶差分或 STL 分解再进入 TAR。门限变量也建议用差分后的 dy，而不是原始 y，否则门限经常切在趋势段上而不是状态段上，这一点在第 5 章会重点展开。
3.2 线性检验：先证明数据确实存在门限效应
跳过线性检验直接 setar 是新手最常见的翻车点，拟合完你会发现门限值落在样本边缘、regime 大小失衡，模型纯粹是硬搜出来的。正规做法是先做线性检验。常见可用的是 Tsay (1989) 的排序回归检验，思路是：把样本按门限变量 y_{t-d} 升序排序，按顺序做递推回归，得到一组递推残差；如果模型是线性 AR，递推残差与回归变量无关；如果存在门限切换，递推残差会被回归变量系统性解释，F 统计量显著。
## Tsay 排序回归检验的示意实现 tsay_f <- function(y, p = 2, d = 1) { mat <- embed(y, p + 1) # 列：y_t, y_{t-1}, ..., y_{t-p} yt <- mat[, 1] X <- cbind(1, mat[, -1]) # 截距 + p 个滞后 z <- mat[, d + 1] # 门限变量 y_{t-d} ord <- order(z) yt <- yt[ord]; X <- X[ord, , drop = FALSE] e <- rep(NA, nrow(X)) for (i in (p + 2):nrow(X)) { b <- solve(crossprod(X[1:i, ]), crossprod(X[1:i, ], yt[1:i])) e[i] <- yt[i] - sum(X[i, ] * b) } e_valid <- e[(p + 2):nrow(X)] X_valid <- X[(p + 2):nrow(X), , drop = FALSE] fit <- lm(e_valid ~ X_valid - 1) sse_full <- sum(residuals(fit)^2) sse0 <- sum(e_valid^2) k <- ncol(X) n <- length(e_valid) F <- ((sse0 - sse_full) / k) / (sse_full / (n - k)) c(f_stat = F, p_value = 1 - pf(F, k, n - k)) } tsay_f(dy, p = 2, d = 1)代码里 crossprod 是纯矩阵运算，比循环里调 lm 快很多，小样本无所谓，但习惯上我会一直用这种写法。输出的 p 值如果小于 0.05，拒绝线性 AR 原假设，说明存在门限效应，可以继续；如果 p 值大于 0.1，停止，回去做 ARIMA 或别的模型，不要硬上 TAR。这一步才是 TAR 建模真正的门槛。
3.3 setar 拟合与结果解读：门限、系数、AIC
确认存在门限效应后，用 tsDyn 包的 setar 拟合。这是 R 里门限自回归（Threshold model）最常用的实现，支持自动搜索门限。
library(tsDyn) ## m 是 regime 内 AR 阶数，d 是延迟参数 ## thDelay 是门限变量的滞后，通常与 d 对齐 ## th = "search" 表示让模型在候选门限中自动搜索 mod_tar <- setar(dy, m = 2, d = 1, thDelay = 1, th = "search", trace = FALSE) summary(mod_tar)summary 输出里重点看三块。第一是 Threshold 那一行，给出最终选中的门限值 r，以及落入两个 regime 的样本数——两个 regime 的样本量如果比例悬殊（比如 95:5），说明这个门限没有太多实际意义。第二是两套 regime 系数，注意截距和滞后系数在 regime 间的符号、量级差异，这直接对应业务上的切换；第三是 AIC，后面模型对比要靠它。强烈建议把 trace 改成 TRUE 重新跑一遍，打印模型尝试过的门限候选和对应 SSR，你能看到 RSS 在 r 附近是不是一个清晰的谷——如果是很平的沟，说明 r 辨识度差，要在第 5 章的处理里认真对待。thDelay 控制门限变量取哪一阶滞后，不同版本对基准期定义略有差异，跑通后把 summary 显示的门限变量和 lag(dy, 1) 画在一起核对，确保两个序列一致。
setar 拟合完之后，先做残差自相关检查再谈预测。残差诊断和普通时间序列回归完全一样，Ljung-Box 检验直接跑：
Box.test(residuals(mod_tar), lag = 12, type = "Ljung-Box")p 值大于 0.05 说明残差没有显著自相关，模型基本的动态结构是对的；显著则说明门限只切了均值或截距，滞后结构还没吸干净。先不要急着调 p，回到 2.1 的结构图再确认 regime 定义对不对。
4. 门槛值搜索、延迟参数 d 与滞后阶数：三个必调参数
4.1 d 的选择：从线性检验 p 值反推
d 是最容易被忽视的参数。很多人默认 d=1，但如果切换不是隔一期就触发，模型会把门限嵌进错误的位置。我一般这样排：先针对 d=1 到 d=6 分别跑一次 tsay_f，记录 p 值；p 值最小（最显著）的 d 表示在这个延迟下门限结构最容易辨识，就拿它做初值；然后用对应的 setar 拟合看 AIC，候选的几个 d 之间差距小于 2 个 AIC 单位就选业务上更说得通的那个。
这里有个细节需要提醒：d 和 p 的关系。d 越大，门限变量离当期越远，切换滞后越明显，但有效样本量也会减少（embed 会吃掉前 max(p, d) 个样本）。月度数据里 d=12 等于拿「去年同月」的状态做门限，这会让模型退化成季节性对比，绝大多数情况下不是你想要的门限效应。我在金融波动率数据上见过 d=2 或 d=3 的合理案例，但 jasa_03m 这类水文、负荷数据，d=1 通常是首选。
4.2 门限搜索的样本边界与最小 regime 约束
setar 的自动搜索虽然方便，但它默认的门限候选范围不见得适合你手里的样本。门限搜索的本质是：对门限变量的排序序列，逐个尝试候选切点。Chan 的经典建议是在门限变量的 10%-90% 分位数之间搜索，因为靠近两端时某个 regime 样本太少，参数估计方差会爆炸。
## 手动门限网格搜索，理解 setar 内部在干什么 grid_search_r <- function(y, p = 2, d = 1, probs = seq(0.15, 0.85, by = 0.05)) { mat <- embed(y, p + 1) yt <- mat[, 1] X <- mat[, -1, drop = FALSE] # p 个滞后列 z <- mat[, d + 1] # y_{t-d} 作为门限变量 cand <- as.numeric(quantile(z, probs)) aic_vec <- numeric(length(cand)) for (j in seq_along(cand)) { r1 <- z <= cand[j] r2 <- z > cand[j] n_min <- p + 2 # 每个 regime 至少 p+2 个样本 if (sum(r1) < n_min || sum(r2) < n_min) { aic_vec[j] <- Inf next } fit1 <- lm(yt[r1] ~ X[r1, , drop = FALSE]) fit2 <- lm(yt[r2] ~ X[r2, , drop = FALSE]) rss <- sum(residuals(fit1)^2) + sum(residuals(fit2)^2) k <- 2 * (p + 1) + 1 # 两套回归系数 + 门限本身 n <- nrow(mat) aic_vec[j] <- n * log(rss / n) + 2 * (k + 1) } data.frame(threshold = cand, aic = aic_vec) } grid_search_r(dy, p = 2, d = 1)跑出来的 threshold 列是候选门限值，aic 列是对应 AIC。你会看到 AIC 随门限变化的曲线：要么有一个清晰的谷（门限辨识度高），要么一片平坦（门限辨识度低）。这个函数虽然简陋，但它帮你把「搜索结果」变成「搜索风景」，很多判断比单看一个最优 r 值可靠。n_min 用 p+2 而不是 p+1，因为两个 regime 各估 p+1 个参数，要留出一点自由度。把输出画成折线图会更直观：
res_r <- grid_search_r(dy, p = 2, d = 1) plot(res_r$threshold, res_r$aic, type = "b")提示：如果曲线在谷底附近有多个候选点 AIC 几乎一样，说明 r 的辨识度低。此时不要执着于那个「最优」门限，先检查样本量，再用 bootstrap 给出 r 的区间。
4.3 模型对比：AR、TAR 与滚动回归，用 AIC 和预测 MSE 定胜负
模型定完形态之后，要拿它和基准模型比。我的对比表长这样：
| 比较维度 | 线性 AR(p) | TAR(p,d,r) | 滚动窗口回归 |
|---|---|---|---|
| 参数个数 | p+1 | 2(p+1)+1 | 随窗口变化 |
| 切换能力 | 无 | 两个 regime | 窗口外全丢弃 |
| AIC 对比 | 基准 | 更低才算赢 | 不可直接比 |
| 切换期预测表现 | 稳定但滞后 | 切换期明显占优 | 滞后于切换 |
对比方法上，我会把样本最后 6 到 12 个点留出来，不参与训练和门限搜索，然后做多步预测，比较均方误差。这一步能揭穿很多「AIC 漂亮但预测翻车」的模型——特别是一旦 d 和 r 是联合搜出来的，AIC 往往带乐观偏误。滚动回归做基准有一个好处：它对切换的响应是无模型的，如果 TAR 连滚动回归都比不过，说明门限结构没有转化为预测能力，应回头检查 4.2 的搜索风景。
这些参数设定最后要落成一句话：搜索越少越好。固定 d，限制 r 的搜索范围，p 用 AIC 在一个小集合里选，这是 TAR 建模里最不容易翻车的配置。如果你只有一百多个观测，却想同时搜出最优的 p、d、r，那模型自由度已经被数据撑爆了，任何一个信息准则都会骗你。
5. 门限自回归的避坑清单：从跑通到跑对的 5 个教训
5.1 伪门限：门限切在趋势段而不是状态段
现象：模型摘要很漂亮，门限值显著，但样本外预测比线性 AR 还差。折线图上看，门限把前段和后段的时间切开了。
原因：数据用原始水平值建模，趋势没处理干净。门限变量 y_{t-d} 在该数据里单调上升，门限实际在切时间区间，而不是切状态区间。
解决：用差分序列 dy 做门限变量和回归变量，或者在建模前先做 STL 分解并只对余项建 TAR。判断标准：把门限值 r 对应的位置标在 y_t 折线图上，如果 r 把时间轴切成几段，而不是把同一时段内的观测按状态分开，基本就是伪门限。
5.2 门限估计不稳定：r 的值随 p、d 变化剧烈
现象：把 p 从 2 改成 3，或者把 d 从 1 改成 2，门限估计值从 0.5 跳到 2.1，业务解释从「低水位」变成「中水位」。
原因：RSS 在门限 r 附近非常平坦，很多候选点给出的残差平方和相差无几；或者样本量小，个别点主导了最优选择。门限参数的得分面不是光滑的，极值附近没有解析标准误可用。
解决：画 grid_search_r 的 AIC 曲线看谷是否清晰；用 bootstrap 重抽样重新估计 r，得到 r 的分布区间。汇报结果时只报点估计是不够的，我会标出 r 的 90% bootstrap 区间，区间宽度超过序列标准差的 1/3 就该对 r 的精确值保持谨慎。
5.3 拟合完残差仍有自相关：门限没把结构切干净
现象：summary 的 AIC 很低，但 residuals(mod_tar) 的 Ljung-Box p 值小于 0.05，ACF 图在 lag 12 或 lag 24 突出。
原因：常见两种。一种是季节性没预处理，TAR 本身不处理周期成分，门限只是把季节性分成了两段；另一种是 regime 内滞后阶数 p 太小，动态记忆不够。
解决：先确认预处理（差分、STL、季节哑变量），再做 TAR；排除预处理后仍显著，把 regime 内 p 增大一阶重估，但注意 p 变大对两个 regime 样本量的挤压——这时要回头收紧 4.2 的最小 regime 约束。不要试图靠多加一个门限去吸收季节性，那会把模型参数翻倍，换来的却是季节性伪拟合。
5.4 AIC 联合搜索过拟合：训练集好看，测试集翻车
现象：p、d、r 全部交给自动搜索，训练 AIC 远优于线性 AR，但回测的滚动 MSE 更差。
原因：联合搜索的自由度太大，r 被搜到恰好拟合训练残差的位置；AIC 的惩罚项对非线性参数的复杂性补偿不够，门限又是一个让似然面很不光滑的参数。
解决：把 d 固定为 1 或用线性检验选定，只让 r 和 p 参与搜索；做至少 6 步的样本外预测对比；在两个 regime 样本量接近 50:50 时才承认门限的平均效应可靠，比例悬殊时优先怀疑过拟合。
5.5 一个门限不够：残差在 regime 内仍分层
现象：双 regime TAR 拟合后，把两个 regime 的残差分别画在 y_{t-d} 的排序图上，还能看出分段结构；把 regime 1 的样本单独做一次线性检验又显著。
原因：数据存在两个以上的切换点，比如「低位平稳—中位波动—高位回落」三段结构，双 regime 模型被迫把后两段合并。
解决：先做残差再检验：用 regime 1 样本重跑 tsay_f，若仍拒绝线性，考虑三 regime TAR，或换 LSTAR（平滑转移）避免硬门限。三 regime 模型的参数是双 regime 的 1.5 倍，样本没有 300 个观测以上我不建议硬上，否则会掉进 5.4 的过拟合坑。
6. 残差诊断与滚动预测对比：验证门限自回归是不是真的更优
6.1 残差诊断：不只是看 Ljung-Box
模型拟合完，我至少跑三组检查：Ljung-Box 看自相关，QQ 图看正态性，以及把残差按门限变量排序后画散点，看是否还有分段残留。最后这一步是 TAR 专属的验证方法，普通时间序列回归不需要，但门限模型必须做——如果原假设是「门限把机制切干净了」，那么残差对门限变量排序后应该没有趋势和分层。
## 残差白噪声检查 Box.test(residuals(mod_tar), lag = 12, type = "Ljung-Box") ## 按门限变量排序后的残差散点：不应再有分段结构 z <- lag(dy, 1)[-1] e <- residuals(mod_tar) plot(z, e) lines(lowess(z, e), col = "red")lowess 曲线如果还是明显的 V 形或阶梯形，说明当前门限没切到位。先别急着换模型，回到第 4 章重新审视 d 和 r 的搜索范围，往往比直接加 regime 更有效。
6.2 滚动预测对比的完整骨架
验证 TAR 是否真的更优，我常用的做法是把最后 6 个观测留出做测试集，AR 和 TAR 都用前段训练，比较多步预测的 MSE。
n <- length(dy) cut <- n - 6 train <- dy[1:cut] test <- dy[(cut + 1):n] ## 线性 AR 基准 mod_ar <- ar(train, order.max = 3) pred_ar <- as.numeric(predict(mod_ar, n.ahead = 6)$pred) ## TAR 模型：训练集重新拟合，不接触测试集 mod_tar <- setar(train, m = 2, d = 1, thDelay = 1, th = "search") pred_tar <- as.numeric(predict(mod_tar, n.ahead = 6)) mse_ar <- mean((test - pred_ar)^2) mse_tar <- mean((test - pred_tar)^2) cat("AR MSE:", mse_ar, "TAR MSE:", mse_tar, "\n")注意测试集不能参与门限搜索，否则 r 会被测试集的信息污染。TAR 的 MSE 如果只是略低，我更看重它在切换点附近的单步表现——把预测误差按时间画出来，看 TAR 是否在机制切换后的第一二步明显收窄。我这两年做门限建模养成的习惯是：先画 y_t 对 y_{t-1} 的散点，手动标出可能的切点，再让模型去搜。这个习惯帮我挡住了至少一半的伪门限。门限自回归好用，但它的每一步都得对数据有交代，参数少搜一个是一个，希望帮到你。
本文还有配套的精品资源，点击获取
