简介门限自回归TAR模型能为存在机制转换或临界效应的非线性时间序列提供灵活的分段建模方案这份压缩包面向需要运用MATLAB完成TAR建模的研究者与数据分析学习者。包内共6个文件以3个m脚本为核心配合txt数据与Readme说明文件整体仅70KB轻量集中便于快速运行与对照学习。目前已有357人学习下载。代码示例覆盖了TAR模型构建的主要环节包括阈值识别、不同区间的自回归拟合、最大似然参数估计以及残差诊断同时包含LR似然比图绘制逻辑用以比较增加门限后的似然增益辅助判断最佳门限数量、避免过拟合。配套的txt数据文件可直接用于实验演练帮助读者从数据预处理、阈值检测到模型验证完整走通流程。对于研究经济波动、气象变化等存在状态切换的场景这套MATLAB程序提供了清晰可复用的参考实现。1. 门限自回归当序列在不同状态下来回横跳线性 AR 模型真的顶不住jasa_03m 这类月度序列最常见的毛病是均值、方差和自相关结构在某个临界点前后像换了个人。行情好的时候 GDP 增速的惯性很强跌到某个阈值以下就变成另一种波动模式用电量低峰和高峰的回归系数也不一样。普通 AR 模型只给一套全局系数遇到这种 regime switching 只能把门槛两侧的差异平均掉结果预测在拐点处永远慢半拍。门限自回归Threshold AutoregressiveTAR的思路简单粗暴用门限变量把序列切成几段每段各跑一个 AR让系数自己随状态切换。这篇文章就把 TAR 从原理到落地讲透适合手里有月度、季度观测序列、想改善拐点预测效果的数据分析或量化从业者。2. 为什么序列一有门限效应线性回归就不够用从模型结构说起2.1 TAR 的数学表达分段线性是它最核心的武器门限自回归最早由 Tong 在 1970 年代末提出后来和 Lim 一起整理成书核心思想用一个分段函数就能说清。考虑一个单变量时间序列 $y_t$如果门限变量取的是 $y_{t-d}$也就是序列自己延迟 d 期的值模型叫 SETARSelf-Exciting Threshold Autoregressive这就是最常用的一种形式。两段 SETAR 的写法是$$ y_t \begin{cases} \phi_1^{(0)} \sum_{i1}^{p_1} \phi_1^{(i)} y_{t-i} \varepsilon_t^{(1)}, y_{t-d} \le c \ \phi_2^{(0)} \sum_{i1}^{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广义叫 TVARThreshold 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}$ 做门限效果最好$d1$ 反而会晚两期才切换。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)说明typetrend允许检验回归里带趋势项selectlagsAIC自动选择 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是门限变量的滞后阶数 d1即用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%用分位数平方根等间隔取候选值如果你想更细把by0.01改成by0.002搜索粒度变大但这么做对单变量几百个样本意义不大反而可能搜到局部过拟合点。第三层是分段 OLS 和评价k_total 2*(p1)1是我常用的惩罚项公式两段各 p1 个参数加上 1 个门限值参数用 AIC 辅助比较、用 RSS 做网格选择主依据。4.3 参数调优搜索粒度、最小分段样本量和滞后阶数的取舍自己实现的代码让你看清三组参数的相互作用。第一是搜索粒度与样本量的关系样本量 n300 时1% 分位数步长大约只给 3 个观测的间距门限估计可能过拟合到单个点保守做法是seq(0.15, 0.85, by0.05)只试 15 个候选宁可粗糙不可过头。我通常在探索阶段用 by0.01确认门限效应存在后用 by0.05 重新估看门限值稳不稳。最小分段样本量的设置直接影响低 regime 的估计质量。p2是最低限实际建议至少留 10% 的样本也就是 300 个观测里每段至少 30 个。tsDyn::setar的trim参数默认 0.15就是干这个。如果最优门限本身就在 0.2 分位附近那 trim0.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, typeLjung-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根本不显著。原因分段后每段样本量变小特别是 trim0.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显著再上 TARTAR 拟合后做门限 bootstrap 分布检验稳定才进入正式预测。这个流程帮我避开了至少三次看起来有门限、实际是过拟合的翻车。门限自回归不是银弹但对 JASA_03M 这类存在 regimes 的月度序列它是性价比最高的非线性建模起点——参数少、可解释、分段线性的性质让你能跟业务方说清楚哪个状态在用什么逻辑跑。希望这篇落地笔记能帮你在实战里少走几步弯路。本文还有配套的精品资源点击获取 门限自回归：时间序列状态切换的非线性预测方法
简介门限自回归TAR模型能为存在机制转换或临界效应的非线性时间序列提供灵活的分段建模方案这份压缩包面向需要运用MATLAB完成TAR建模的研究者与数据分析学习者。包内共6个文件以3个m脚本为核心配合txt数据与Readme说明文件整体仅70KB轻量集中便于快速运行与对照学习。目前已有357人学习下载。代码示例覆盖了TAR模型构建的主要环节包括阈值识别、不同区间的自回归拟合、最大似然参数估计以及残差诊断同时包含LR似然比图绘制逻辑用以比较增加门限后的似然增益辅助判断最佳门限数量、避免过拟合。配套的txt数据文件可直接用于实验演练帮助读者从数据预处理、阈值检测到模型验证完整走通流程。对于研究经济波动、气象变化等存在状态切换的场景这套MATLAB程序提供了清晰可复用的参考实现。1. 门限自回归当序列在不同状态下来回横跳线性 AR 模型真的顶不住jasa_03m 这类月度序列最常见的毛病是均值、方差和自相关结构在某个临界点前后像换了个人。行情好的时候 GDP 增速的惯性很强跌到某个阈值以下就变成另一种波动模式用电量低峰和高峰的回归系数也不一样。普通 AR 模型只给一套全局系数遇到这种 regime switching 只能把门槛两侧的差异平均掉结果预测在拐点处永远慢半拍。门限自回归Threshold AutoregressiveTAR的思路简单粗暴用门限变量把序列切成几段每段各跑一个 AR让系数自己随状态切换。这篇文章就把 TAR 从原理到落地讲透适合手里有月度、季度观测序列、想改善拐点预测效果的数据分析或量化从业者。2. 为什么序列一有门限效应线性回归就不够用从模型结构说起2.1 TAR 的数学表达分段线性是它最核心的武器门限自回归最早由 Tong 在 1970 年代末提出后来和 Lim 一起整理成书核心思想用一个分段函数就能说清。考虑一个单变量时间序列 $y_t$如果门限变量取的是 $y_{t-d}$也就是序列自己延迟 d 期的值模型叫 SETARSelf-Exciting Threshold Autoregressive这就是最常用的一种形式。两段 SETAR 的写法是$$ y_t \begin{cases} \phi_1^{(0)} \sum_{i1}^{p_1} \phi_1^{(i)} y_{t-i} \varepsilon_t^{(1)}, y_{t-d} \le c \ \phi_2^{(0)} \sum_{i1}^{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广义叫 TVARThreshold 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}$ 做门限效果最好$d1$ 反而会晚两期才切换。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)说明typetrend允许检验回归里带趋势项selectlagsAIC自动选择 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是门限变量的滞后阶数 d1即用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%用分位数平方根等间隔取候选值如果你想更细把by0.01改成by0.002搜索粒度变大但这么做对单变量几百个样本意义不大反而可能搜到局部过拟合点。第三层是分段 OLS 和评价k_total 2*(p1)1是我常用的惩罚项公式两段各 p1 个参数加上 1 个门限值参数用 AIC 辅助比较、用 RSS 做网格选择主依据。4.3 参数调优搜索粒度、最小分段样本量和滞后阶数的取舍自己实现的代码让你看清三组参数的相互作用。第一是搜索粒度与样本量的关系样本量 n300 时1% 分位数步长大约只给 3 个观测的间距门限估计可能过拟合到单个点保守做法是seq(0.15, 0.85, by0.05)只试 15 个候选宁可粗糙不可过头。我通常在探索阶段用 by0.01确认门限效应存在后用 by0.05 重新估看门限值稳不稳。最小分段样本量的设置直接影响低 regime 的估计质量。p2是最低限实际建议至少留 10% 的样本也就是 300 个观测里每段至少 30 个。tsDyn::setar的trim参数默认 0.15就是干这个。如果最优门限本身就在 0.2 分位附近那 trim0.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, typeLjung-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根本不显著。原因分段后每段样本量变小特别是 trim0.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显著再上 TARTAR 拟合后做门限 bootstrap 分布检验稳定才进入正式预测。这个流程帮我避开了至少三次看起来有门限、实际是过拟合的翻车。门限自回归不是银弹但对 JASA_03M 这类存在 regimes 的月度序列它是性价比最高的非线性建模起点——参数少、可解释、分段线性的性质让你能跟业务方说清楚哪个状态在用什么逻辑跑。希望这篇落地笔记能帮你在实战里少走几步弯路。本文还有配套的精品资源点击获取