企业官网建设流程全解析

简介：门限自回归（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。两者的共性和差异可以用一个判断框架说清楚。

我个人的习惯是：先跑一遍 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_col=0, parse_dates=True).iloc[:, 0] # ADF 检验：p 值小于 0.05 视为平稳 adf_result = adfuller(series.dropna()) print(f"ADF 统计量: {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 波动，可以传入regression='n'去掉常数项，能提高检验功效。还有一点容易忽略——月度数据应该把autolag='AIC'设为默认的自动滞后选择，它会按信息准则挑最优滞后阶数，比固定滞后更稳。

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, lags=ar_order).fit() model_high = AutoReg(s_high, lags=ar_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, delay=1, threshold_candidates=threshold_candidates, ar_order=2) 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~5，d ∈ 1~3，组合数量不大，暴力枚举就行。

import warnings warnings.filterwarnings("ignore") def select_tar_order(series, p_range, d_range, threshold_candidates, min_samples=20): """网格搜索使 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, delay=d, threshold_candidates=threshold_candidates, ar_order=p) 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(key=lambda 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_range=range(1, 6), d_range=range(1, 4), threshold_candidates=series_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, lags=ar_order).fit() model_high = AutoReg(s_high, lags=ar_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_order+1): 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 选滞后阶数，可能在低体制选出 p=3，高体制选出 p=5，两套模型的复杂度加在一起远超数据能支撑的信息量，造成过拟合。更隐蔽的问题是：门限模型每个体制内的观测数少，相同阶数的 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_sim=1000): """残差自举法生成预测区间""" sim_paths = [] for _ in range(n_sim): # 有放回抽样残差 boot_resid = np.random.choice(residuals, size=horizon, replace=True) # 用不带噪声的模型预测初始路径 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, axis=0) upper = np.percentile(sim_paths, 97.5, axis=0) return np.array(base_forecast), lower, upper

代码里最关键的一行是np.random.choice(residuals, size=horizon, replace=True)，它决定了重抽样的随机性来源。如果残差序列存在自相关，直接重抽样会破坏时序结构，此时应该先对残差拟合一个 AR(1) 模型，再抽取白噪声部分做自举。我在月度数据上的经验是：残差如果有明显的序列相关（Ljung-Box 检验 p 值小于 0.05），先对残差做 ARIMA 过滤，再对过滤后的白噪声自举，区间会准很多。预测区间的宽度是模型不确定性的直观标尺：如果区间宽度超过预测值本身的 50%，说明模型在这个数据集上的信息量有限，后面要考虑引入外部变量或更高的采样频率。

个人习惯是：任何门限模型的预测报告里，一律把单点预测和区间一起交付。区间不是锦上添花，是判断模型可靠性的核心依据——如果一位工程师只给我单点预测却不给区间，我基本不会采用。希望这些内容能帮你在做门限自回归的路上少走些弯路，把这套非线性工具真正用起来。

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