简介门限自回归TAR模型是在传统线性AR模型基础上引入阈值机制的非线性时间序列方法适用于刻画不同状态下回归关系发生结构性变化的数据。这份MATLAB代码资源面向时间序列分析方向的学习者与研究人员围绕TAR模型构建与模型选择提供了一套可运行的完整示例。压缩包共6个文件含3个m脚本和3个txt文件体积仅70KB脚本覆盖数据预处理、阈值识别、分段AR拟合、最大似然参数估计及LR似然比图绘制等核心环节txt文件提供示例数据与使用说明便于对照运行。通过调试和运行这套代码读者可掌握TAR模型从建模到诊断的完整流程并学会利用LR图判断最佳门限数量提升对非线性时间序列的建模能力。目前已有356人学习下载适合具备基础回归知识、希望实践门限模型的MATLAB用户。1. 门限自回归时间序列回归里的“变脸”问题终于有解了做时间序列预测的人迟早会撞上一堵墙模型在某个时间点之前拟合得漂亮之后却集体失效。不是数据变脏了而是数据背后的“机制”变了——宏观数据遇到政策拐点风速序列遇到季节转换销售数据遇到大促前后。线性 AR 模型假设过去的关系永远不变但现实里很多序列存在“门限效应”当某个驱动变量跨过阈值序列的行为就会切换成另一套逻辑。门限自回归Threshold Autoregressive Model简称 TAR就是专门为这种“变关系”设计的模型。它不强行用一条曲线去拟合所有时期而是把序列按门限变量的取值切成若干区间每个区间内分别拟合自回归方程。jasa_03m 这类月度数据正是 TAR 模型最典型的用武之地月度观测里天然藏着季节、周期和突变点用门限模型能比常规 ARIMA 多抓住一层结构。本文不绕弯子直接讲清楚门限自回归的建模流程、Python 实现和参数设置给你一份能照着复现的实战路径。2. 门限自回归的核心逻辑为什么要“分段”做回归2.1 线性 AR 的假设盲区关系不是恒定的先看一个最基础的问题门限自回归到底在解决什么。经典的 AR(p) 模型写成 y_t c φ₁y_{t-1} … φ_p y_{t-p} ε_t它隐含的假设是从第一期到最后一期φ 这些系数是不变的。换句话说无论序列处于高位、低位还是中间状态过去值对当前值的影响力度都一样。这个假设在平稳序列上还能应付一旦序列存在结构性变化——比如某个月出台新政策、市场需求突然换挡——单一组的 φ 就完全不够用了。门限自回归的破局思路很直接把序列按某个变量的数值划分成不同的“体制”regime每个体制内各跑各的自回归。最经典的二体制 TAR 模型写成y_t (φ₁₀ φ₁₁y_{t-1} … φ₁ₚy_{t-p}) · I(z_t ≤ c) (φ₂₀ φ₂₁y_{t-1} … φ₂ₚy_{t-p}) · I(z_t c) ε_t其中 z_t 是门限变量c 是门限值I(·) 是指示函数。公式看着复杂意思是当 z_t 小于等于 c 时序列服从第一套自回归方程当 z_t 大于 c 时切换到第二套方程。z_t 可以是滞后值 y_{t-d}这叫 SETAR自激励门限自回归也可以是外生变量这叫 TAR。这个分段机制让模型天然具备“识别变结构”的能力。2.2 TAR 与 SETAR 怎么选门限变量是核心分水岭建模前必须先想清楚你的门限变量是什么。这个问题决定了模型的类别。如果门限变量是序列自身的滞后值比如 z_t y_{t-1}模型就叫 SETAR如果门限变量是另一个外生变量比如利率、汇率、气温模型就是 TAR。两者的共性和差异可以用一个判断框架说清楚。维度SETARTAR门限变量序列自身滞后值 y_{t-d}外生变量 z_t适用场景序列有自我强化的周期或波动聚集序列受外部变量驱动切换机制建模复杂度较低数据本身够用需要配套的门限变量数据典型应用太阳黑子、气温、股价波动销售额受政策利率影响、风速受气压影响我个人的习惯是先跑一遍 SETAR 做摸底因为不需要额外找数据如果检验出明显的门限效应但不稳定再考虑换外生门限变量。见过不少项目是数据明明有突变点却硬套 ARIMA 硬扛改造成门限模型后误差直接降了一截。选择 TAR 还是 SETAR本质是在问你更相信“过去的状态”还是“外部的力量”能解释当前的行为切换。2.3 门限效应的检验不要上来就分段这里有个新手常犯的错误看到序列图上有几个折点就急着定门限值。门限模型的第一个正式步骤应该是假说检验——先验证“是否存在门限效应”再估计门限值。统计上一般用 Chan (1993) 的超检验过程先按门限变量的取值排序对每一个可能作为门限的候选点把样本分成两段分别估计两段的自回归模型算残差平方和RSS取 RSS 最小的那个候选点作为门限估计值。这一步在 R 的 tsDyn 包和 Python 里都有成熟实现。但这还不够因为哪怕没有真正的门限效应你也能找到“使 RSS 最小”的点只是它没有统计显著意义。标准的做法是用自举bootstrap生成零分布在原假设“无门限效应”下模拟一大批序列计算每个模拟序列的超检验统计量然后看真实统计量落在分布的多极端位置算出 p 值。如果 p 值小于 0.05才算有统计依据做分段回归。3. 用 Python 实现门限自回归从数据检验到参数估计3.1 数据准备与平稳性预检验没有一个可靠的过程就没有一个可靠的模型任何门限模型的起点都是数据清洗和平稳性检验。不像线性回归可以容忍一些粗糙的预处理门限模型对数据的“分层”非常敏感如果序列里混着异常值或趋势项门限点会被强行拉到一个错误的位置。先做 ADF 检验确认平稳性不平稳就差分这步没有绕过空间。import numpy as np import pandas as pd from statsmodels.tsa.stattools import adfuller # 加载 jasa_03m 月度数据假设为单列时间序列 series pd.read_csv(jasa_03m.csv, index_col0, parse_datesTrue).iloc[:, 0] # ADF 检验p 值小于 0.05 视为平稳 adf_result adfuller(series.dropna()) print(fADF 统计量: {adf_result[0]:.4f}, p 值: {adf_result[1]:.4f}) # 若 p 值偏大做一阶差分后重新检验 if adf_result[1] 0.05: series_diff series.diff().dropna() print(原序列非平稳使用一阶差分序列) else: series_diff series代码的逻辑很直观先用adfuller做 ADF 单位根检验p 值小说明没有单位根、序列平稳p 值大就做一阶差分。这里有个参数注意点adfuller默认的回归项含常数项c如果序列均值明显非零这个设置是对的如果序列围绕 0 波动可以传入regressionn去掉常数项能提高检验功效。还有一点容易忽略——月度数据应该把autolagAIC设为默认的自动滞后选择它会按信息准则挑最优滞后阶数比固定滞后更稳。3.2 门限候选点的网格搜索把连续阈值变成可计算的离散候选门限值本身是连续的但实际计算时只能从观测值里选。常见做法是把门限变量 z_t 从小到大排序然后按一定的分位数区间比如 15% 到 85%作为候选范围去掉两端太少的样本保证每个体制内的观测数足够做回归。接下来对每个候选点把全样本切两半分别估计 AR 模型累加两个子模型的残差平方和。import itertools from statsmodels.tsa.ar_model import AutoReg def fit_tar_rss(series, delay, threshold_candidates, ar_order): 遍历候选门限点返回 RSS 最小的门限值及对应模型信息 best {rss: np.inf, threshold: None, models: None} for th in threshold_candidates: # 按门限变量这里是 y_{t-delay}分割样本 mask_low series.shift(delay) th mask_high series.shift(delay) th s_low, s_high series.loc[mask_low].dropna(), series.loc[mask_high].dropna() # 每个体制内要求最小样本量否则跳过 if len(s_low) 20 or len(s_high) 20: continue model_low AutoReg(s_low, lagsar_order).fit() model_high AutoReg(s_high, lagsar_order).fit() rss_total model_low.ssr model_high.ssr if rss_total best[rss]: best[rss] rss_total best[threshold] th best[models] (model_low, model_high, mask_low, mask_high) return best # 以滞后 1 期为门限变量候选门限设为序列 15%-85% 分位区间内的观测值 series_clean series_diff.dropna() threshold_candidates series_clean.shift(1).dropna().quantile([0.15, 0.25, 0.5, 0.75, 0.85]).values result fit_tar_rss(series_clean, delay1, threshold_candidatesthreshold_candidates, ar_order2) print(f最优门限值: {result[threshold]:.4f}, 最小 RSS: {result[rss]:.4f})这段代码的核心是fit_tar_rss函数里的双重循环逻辑外层遍历门限候选点内层对分段后的两个子样本分别拟合AutoReg模型。model_low.ssr和model_high.ssr分别是低体制和高体制模型的残差平方和加起来就是总 RSS。延迟参数delay决定门限变量用的是哪一期滞后——这个参数直接影响模型行为后面会专门讲怎么选。我在项目里的经验是候选门限别用手工指定用我这段的分位数切片方式更客观5 个分位点做粗筛找到最优区间后在区间内细化网格重跑一遍能有效避免漏掉真正的门限值。3.3 固定滞后阶数与门限延迟用 AIC/BIC 在模型空间里选择门限值定下来之后还需要确定两件事每个体制内的 AR 滞后阶数 p以及门限延迟 d如果用序列自身做门限变量就是 z_t y_{t-d} 里的 d。这两组参数的最佳组合通常靠信息准则搜索确定。常见的搜索范围是 p ∈ 1~5d ∈ 1~3组合数量不大暴力枚举就行。import warnings warnings.filterwarnings(ignore) def select_tar_order(series, p_range, d_range, threshold_candidates, min_samples20): 网格搜索使 AIC 最小的 (p, d) 组合 results [] for p, d in itertools.product(p_range, d_range): # 门限变量使用 y_{t-d}序列做回归时使用 1..p 期滞后 try: res fit_tar_rss(series, delayd, threshold_candidatesthreshold_candidates, ar_orderp) if res[threshold] is None: continue model_low, model_high, _, _ res[models] aic_total model_low.aic model_high.aic results.append((aic_total, p, d, res[threshold], res[rss])) except Exception: continue results.sort(keylambda x: x[0]) best_aic, best_p, best_d, best_th, best_rss results[0] print(f最优 AIC{best_aic:.2f}, p{best_p}, d{best_d}, 门限{best_th:.4f}) return results result_sorted select_tar_order( series_clean, p_rangerange(1, 6), d_rangerange(1, 4), threshold_candidatesseries_clean.shift(1).dropna().quantile(np.arange(0.2, 0.8, 0.1)).values )这里我用AIC而不是BIC—— TAR 模型本身分段估计参数已经比线性模型多了一倍再用 BIC 这种偏保守的准则容易选出过简模型AIC 在样本量中等的情况下更兼顾拟合与复杂度。不过如果你的样本量很大超过 500 个观测换成 BIC 问题也不大最终效果差异很小。还有个细节AutoReg的ssr属性和aic属性直接获取残差平方和与 Akaike 信息准则值省去手写公式的麻烦。跑完网格后注意检查最优组合和第二优组合的 AIC 差距如果差距很小小于 2说明模型选择不够稳健需要检查数据里有没有异常值扰动。4. jasa_03m 月度数据的实战建模从检验到预测的完整流程4.1 为什么月度数据特别适合门限模型我拿 jasa_03m 这批月度数据举例因为它非常有代表性。月度观测天然具有三种结构第一是季节性周期比如 12 个月的景气循环第二是趋势变化比如经济扩张期和收缩期第三是变点事件比如疫情冲击、政策调整、供需关系的结构性转变。这三者叠加在一个线性模型里往往互相掩盖——季节项把突变点抹平趋势项把短期的门限效应吞掉。门限模型处理月度数据的优势在于它能自动根据门限变量可以是滞后值也可以是外部月度指标把数据分成“高体制”和“低体制”比如物价高企的月份和物价低迷的月份两个体制内部分别拟合不同的自回归结构捕捉到的依赖关系就完全不一样。你可能已经注意到这有点像把“时间分段”换成“状态分段”——不是按时间轴切一刀而是按变量取值切一刀这样即便相同月份处于不同年份只要状态相似就会被分到同一个体制里估计参数信息的利用率高得多。用 jasa_03m 这类数据做实证时我通常先把序列拆成训练集和测试集比如前 80% 训练后 20% 验证在训练集上完成全部参数选择测试集只用于最终的滚动预测评估。4.2 带外生门限变量的 TAR把外部驱动因素引入模型做月度数据时门限变量只选序列自身的滞后值有时不够。比如你预测一个城市的月度用电量门限变量用“上月气温是否超过某个阈值”比用“上个月用电量是否超过某个值”更贴近物理实际。这就回到 2.2 节讲的 TAR 和 SETAR 的分野。带外生门限变量的实现方式与纯 SETAR 几乎一样只是门限变量的数据来源变了。def fit_tar_exog(series, exog_threshold, threshold_candidates, ar_order): 外生门限变量的 TAR 拟合 exog_threshold: 与 series 等长的外生门限变量如气温、利率 best {rss: np.inf, threshold: None, models: None} for th in threshold_candidates: mask_low exog_threshold th mask_high exog_threshold th s_low, s_high series.loc[mask_low].dropna(), series.loc[mask_high].dropna() if len(s_low) 20 or len(s_high) 20: continue model_low AutoReg(s_low, lagsar_order).fit() model_high AutoReg(s_high, lagsar_order).fit() rss_total model_low.ssr model_high.ssr if rss_total best[rss]: best[rss] rss_total best[threshold] th best[models] (model_low, model_high, mask_low, mask_high) return best和 3.2 的fit_tar_rss一比唯一的区别是第四行的exog_threshold替代了原有的series.shift(delay)其他一切照旧。看起来改动不大但这个选择的背后是建模思路的变化当门限变量来自外部时你得先确认它的“时间对齐”——是同步值当月气温还是滞后值上月气温这取决于业务上哪个变量“驱动”了序列切换机制。做外生门限变量时最容易翻车的就是对齐问题门限变量和序列的索引对应错位一位整个模型就串味了。我在项目里会先用print(pd.DataFrame({series: series, gate: exog_threshold}).dropna().head())检查对齐再做拟合这一步五分钟能省三小时的排查功夫。4.3 模型预测与双体制的样本外对比模型估计完最终要回答的问题是它比单一线性 AR 模型强在哪儿验证方法不复杂——用训练好的门限模型对测试集逐期做预测并和线性 AR 模型做误差对比。门限模型的预测逻辑是每一步预测前先看门限变量落在哪个体制再用对应体制的 AR 模型来预测。预测值可以直接用AutoReg的predict方法但门限模型的预测要小心滞后阶数的对齐推荐把预测函数封装好。def tar_forecast(series, exog_threshold, th, model_low, model_high, ar_order, horizon): 门限模型滚动预测, h 步预测 forecasts [] history list(series.values) gate_history list(exog_threshold.values) for _ in range(horizon): gate_val gate_history[-1] # 当前状态决定用哪个体制的模型预测 active_model model_low if gate_val th else model_high # 取最近 ar_order 个观测作为起点 recent history[-ar_order:] # 手动构造滞后特征 X np.array([recent[::-1]]).reshape(1, -1) yhat active_model.predict(X)[0] if hasattr(active_model, predict) else np.nan # 手动用参数算滞后系数 常数项 if np.isnan(yhat): yhat active_model.params.iloc[0] # 常数 for lag_i in range(1, ar_order1): yhat active_model.params.iloc[lag_i] * recent[-lag_i] forecasts.append(yhat) # 更新历史序列 history.append(yhat) gate_history.append(gate_val) return forecasts这段代码用的active_model.predict(X)是理论化写法实际用statsmodels的AutoReg对象做单步预测时需要先拟合一个AutoReg的预测包装器或者直接读取参数做滚动推算。我更推荐后一种方式把active_model.params按常数和滞后系数拆开手动加权计算下一期预测值逻辑透明、不容易因为时序索引错位而出错。horizon是预测步长月度数据一般做 3~12 步预测。预测结束后对比mean_absolute_error(forecasts, actual)和线性 AR 的同一指标如果门限模型没有显著优势就回头检查门限值是不是不够显著或者门限体制的样本量差异太大——这些在下一章详细说。5. 门限自回归避坑指南五个最容易翻车的地雷5.1 现象门限效应检验的 p 值不稳定换个样本区间结果就变原因门限效应检验依赖自举抽样的随机性不同的随机种子会产生不同的模拟分布p 值落在 0.04 到 0.06 之间浮动时结论就会摇摆。另一个原因是样本量不够自举分布的形状不够稳定。解决跑检验时固定随机种子np.random.seed(42)并且至少跑三次不同的种子观察 p 值是否稳定地在 0.05 以下。如果 p 值在临界点附近晃动别急着建模先用数据可视化确认门限变量的分布是否存在双峰或明显的跳跃结构。我在月度数据上踩过这个坑第一次跑 p 值是 0.03换了个子样本变成 0.08最后发现是样本里有一个极端的异常值把门限点拉偏了剔除后 p 值稳定在 0.01。5.2 现象模型在训练集上拟合得很好测试集上一塌糊涂原因这是门限模型最典型的过拟合陷阱——搜索门限值的时候用的是“残差平方和最小”这个目标但是候选门限点本身是从观测值里挑的你其实是在用测试信息反选门限。本质上你是在拿数据来拟合门限的位置而不是先验地指定它。这会导致门限点选择过度优化训练集内部再分段做回归自由度大幅上升泛化能力自然大幅下降。解决这一点我算交过学费的。后来固定做法是——门限值只从训练集的前 80% 里选后 20% 完全不参与门限值搜索选完门限后再用全训练集估参数。这样至少保证门限点的选择没有“偷看”最后一段信息。另外给每个体制设定最小观测数我通常设为总样本的 15% 以上防止门限切在边缘导致单体制样本太少。5.3 现象两个体制的滞后阶数都用 AIC 定但模型整体比线性 AR 还差原因每个体制单独用 AIC 选滞后阶数可能在低体制选出 p3高体制选出 p5两套模型的复杂度加在一起远超数据能支撑的信息量造成过拟合。更隐蔽的问题是门限模型每个体制内的观测数少相同阶数的 AR 模型在小样本下估计方差更大。解决我一般会先固定两个体制使用相同的滞后阶数再用联合 AIC两个体制 AIC 之和去选。如果不同体制的滞后阶数差异在业务上站不住脚就保持同构。还有一种更实用的策略残差诊断用 Ljung-Box 检验如果残差已经近似白噪声就不要再加滞后阶数了——AIC 在这个场景下容易钩住多余的主项。5.4 现象门限值估计出来正好落在数据的中位数附近但切换逻辑说不通原因这不一定是你选错了但也可能是门限变量选错了。比如用 y_{t-1} 做门限门限值切出来的两段刚好把数据按高低排序切开但这两段没有业务上的“机制差异”只是数学上的偶然。这种模型的预测能力往往很脆弱换一个数据集门限值就漂移。解决回到业务层面问一个问题这个门限变量的哪个数值在你的领域里意味着“机制切换”比如对销量数据门限可能是“上期销量高于历史均值 1.2 倍”——这是产能瓶颈点对风速数据门限可能是“切变层高度超过某个海拔”——这是风切变机制切换点。如果门限值给不出业务解释就换门限变量别硬撑着用数学结果。5.5 现象带外生门限变量时两种体制的参数看起来没有显著差异原因这说明门限变量选的没有区分度——它在两个区间内对应的序列行为本质上一样分段只是给模型额外塞了一堆参数。另一个常见原因是外生门限变量和序列同期相关导致门限效应被序列自身的惯性吸收门限切换本身没有增量信息。解决在拟合门限变量之前先分别统计两个体制内序列的均值、方差和自相关系数。如果三个统计量没有显著差异就不要做门限模型。可以用简单的 t 检验比较两个体制的均值差异是否显著。我自己有个铁律门限模型如果只是让 AIC 降低了不到 5%这个模型就没有实用价值不如老老实实用线性 ARIMA。6. 残差自举法计算预测区间门限模型的不确定性量化最后一节是门限模型里最实用但很多人没做到位的一环预测区间。区间的价值在于帮你区分“正常的波动”和“异常的偏移”。在月度数据上一个未来 3 个月的预测值如果没有区间业务方根本不知道数字是大概率事件还是赌运气。残差自举法的操作分四步第一步从门限模型的残差序列里随机抽样第二步把抽样残差叠加到预测值上生成一条模拟的未来路径第三步重复 1000 次收集所有模拟路径在每个未来时间点的值第四步取 2.5% 和 97.5% 分位数作为 95% 置信区间。这个方法的优势是不需要假设残差服从正态分布——月度数据经常有偏态用正态假设会算出过窄的区间。def tar_forecast_interval(residuals, series, exog_threshold, th, model_low, model_high, ar_order, horizon, n_sim1000): 残差自举法生成预测区间 sim_paths [] for _ in range(n_sim): # 有放回抽样残差 boot_resid np.random.choice(residuals, sizehorizon, replaceTrue) # 用不带噪声的模型预测初始路径 base_forecast tar_forecast(series, exog_threshold, th, model_low, model_high, ar_order, horizon) # 叠加残差生成模拟路径 sim_path np.array(base_forecast) boot_resid sim_paths.append(sim_path) sim_paths np.array(sim_paths) lower np.percentile(sim_paths, 2.5, axis0) upper np.percentile(sim_paths, 97.5, axis0) return np.array(base_forecast), lower, upper代码里最关键的一行是np.random.choice(residuals, sizehorizon, replaceTrue)它决定了重抽样的随机性来源。如果残差序列存在自相关直接重抽样会破坏时序结构此时应该先对残差拟合一个 AR(1) 模型再抽取白噪声部分做自举。我在月度数据上的经验是残差如果有明显的序列相关Ljung-Box 检验 p 值小于 0.05先对残差做 ARIMA 过滤再对过滤后的白噪声自举区间会准很多。预测区间的宽度是模型不确定性的直观标尺如果区间宽度超过预测值本身的 50%说明模型在这个数据集上的信息量有限后面要考虑引入外部变量或更高的采样频率。个人习惯是任何门限模型的预测报告里一律把单点预测和区间一起交付。区间不是锦上添花是判断模型可靠性的核心依据——如果一位工程师只给我单点预测却不给区间我基本不会采用。希望这些内容能帮你在做门限自回归的路上少走些弯路把这套非线性工具真正用起来。本文还有配套的精品资源点击获取