戰(zhàn)：Python分段建模時(shí)間序列突變)
簡(jiǎn)介門限自回歸TAR模型是在傳統(tǒng)線性AR模型基礎(chǔ)上引入閾值機(jī)制的非線性時(shí)間序列方法適用于刻畫不同狀態(tài)下回歸關(guān)系發(fā)生結(jié)構(gòu)性變化的數(shù)據(jù)。這份MATLAB代碼資源面向時(shí)間序列分析方向的學(xué)習(xí)者與研究人員圍繞TAR模型構(gòu)建與模型選擇提供了一套可運(yùn)行的完整示例。壓縮包共6個(gè)文件含3個(gè)m腳本和3個(gè)txt文件體積僅70KB腳本覆蓋數(shù)據(jù)預(yù)處理、閾值識(shí)別、分段AR擬合、最大似然參數(shù)估計(jì)及LR似然比圖繪制等核心環(huán)節(jié)txt文件提供示例數(shù)據(jù)與使用說明便于對(duì)照運(yùn)行。通過調(diào)試和運(yùn)行這套代碼讀者可掌握TAR模型從建模到診斷的完整流程并學(xué)會(huì)利用LR圖判斷最佳門限數(shù)量提升對(duì)非線性時(shí)間序列的建模能力。目前已有356人學(xué)習(xí)下載適合具備基礎(chǔ)回歸知識(shí)、希望實(shí)踐門限模型的MATLAB用戶。1. 門限自回歸時(shí)間序列回歸里的“變臉”問題終于有解了做時(shí)間序列預(yù)測(cè)的人遲早會(huì)撞上一堵墻模型在某個(gè)時(shí)間點(diǎn)之前擬合得漂亮之后卻集體失效。不是數(shù)據(jù)變臟了而是數(shù)據(jù)背后的“機(jī)制”變了——宏觀數(shù)據(jù)遇到政策拐點(diǎn)風(fēng)速序列遇到季節(jié)轉(zhuǎn)換銷售數(shù)據(jù)遇到大促前后。線性 AR 模型假設(shè)過去的關(guān)系永遠(yuǎn)不變但現(xiàn)實(shí)里很多序列存在“門限效應(yīng)”當(dāng)某個(gè)驅(qū)動(dòng)變量跨過閾值序列的行為就會(huì)切換成另一套邏輯。門限自回歸Threshold Autoregressive Model簡(jiǎn)稱 TAR就是專門為這種“變關(guān)系”設(shè)計(jì)的模型。它不強(qiáng)行用一條曲線去擬合所有時(shí)期而是把序列按門限變量的取值切成若干區(qū)間每個(gè)區(qū)間內(nèi)分別擬合自回歸方程。jasa_03m 這類月度數(shù)據(jù)正是 TAR 模型最典型的用武之地月度觀測(cè)里天然藏著季節(jié)、周期和突變點(diǎn)用門限模型能比常規(guī) ARIMA 多抓住一層結(jié)構(gòu)。本文不繞彎子直接講清楚門限自回歸的建模流程、Python 實(shí)現(xiàn)和參數(shù)設(shè)置給你一份能照著復(fù)現(xiàn)的實(shí)戰(zhàn)路徑。2. 門限自回歸的核心邏輯為什么要“分段”做回歸2.1 線性 AR 的假設(shè)盲區(qū)關(guān)系不是恒定的先看一個(gè)最基礎(chǔ)的問題門限自回歸到底在解決什么。經(jīng)典的 AR(p) 模型寫成 y_t c φ?y_{t-1} … φ_p y_{t-p} ε_(tái)t它隱含的假設(shè)是從第一期到最后一期φ 這些系數(shù)是不變的。換句話說無論序列處于高位、低位還是中間狀態(tài)過去值對(duì)當(dāng)前值的影響力度都一樣。這個(gè)假設(shè)在平穩(wěn)序列上還能應(yīng)付一旦序列存在結(jié)構(gòu)性變化——比如某個(gè)月出臺(tái)新政策、市場(chǎng)需求突然換擋——單一組的 φ 就完全不夠用了。門限自回歸的破局思路很直接把序列按某個(gè)變量的數(shù)值劃分成不同的“體制”regime每個(gè)體制內(nèi)各跑各的自回歸。最經(jīng)典的二體制 TAR 模型寫成y_t (φ?? φ??y_{t-1} … φ??y_{t-p}) · I(z_t ≤ c) (φ?? φ??y_{t-1} … φ??y_{t-p}) · I(z_t c) ε_(tái)t其中 z_t 是門限變量c 是門限值I(·) 是指示函數(shù)。公式看著復(fù)雜意思是當(dāng) z_t 小于等于 c 時(shí)序列服從第一套自回歸方程當(dāng) z_t 大于 c 時(shí)切換到第二套方程。z_t 可以是滯后值 y_{t-d}這叫 SETAR自激勵(lì)門限自回歸也可以是外生變量這叫 TAR。這個(gè)分段機(jī)制讓模型天然具備“識(shí)別變結(jié)構(gòu)”的能力。2.2 TAR 與 SETAR 怎么選門限變量是核心分水嶺建模前必須先想清楚你的門限變量是什么。這個(gè)問題決定了模型的類別。如果門限變量是序列自身的滯后值比如 z_t y_{t-1}模型就叫 SETAR如果門限變量是另一個(gè)外生變量比如利率、匯率、氣溫模型就是 TAR。兩者的共性和差異可以用一個(gè)判斷框架說清楚。維度SETARTAR門限變量序列自身滯后值 y_{t-d}外生變量 z_t適用場(chǎng)景序列有自我強(qiáng)化的周期或波動(dòng)聚集序列受外部變量驅(qū)動(dòng)切換機(jī)制建模復(fù)雜度較低數(shù)據(jù)本身夠用需要配套的門限變量數(shù)據(jù)典型應(yīng)用太陽(yáng)黑子、氣溫、股價(jià)波動(dòng)銷售額受政策利率影響、風(fēng)速受氣壓影響我個(gè)人的習(xí)慣是先跑一遍 SETAR 做摸底因?yàn)椴恍枰~外找數(shù)據(jù)如果檢驗(yàn)出明顯的門限效應(yīng)但不穩(wěn)定再考慮換外生門限變量。見過不少項(xiàng)目是數(shù)據(jù)明明有突變點(diǎn)卻硬套 ARIMA 硬扛改造成門限模型后誤差直接降了一截。選擇 TAR 還是 SETAR本質(zhì)是在問你更相信“過去的狀態(tài)”還是“外部的力量”能解釋當(dāng)前的行為切換。2.3 門限效應(yīng)的檢驗(yàn)不要上來就分段這里有個(gè)新手常犯的錯(cuò)誤看到序列圖上有幾個(gè)折點(diǎn)就急著定門限值。門限模型的第一個(gè)正式步驟應(yīng)該是假說檢驗(yàn)——先驗(yàn)證“是否存在門限效應(yīng)”再估計(jì)門限值。統(tǒng)計(jì)上一般用 Chan (1993) 的超檢驗(yàn)過程先按門限變量的取值排序?qū)γ恳粋€(gè)可能作為門限的候選點(diǎn)把樣本分成兩段分別估計(jì)兩段的自回歸模型算殘差平方和RSS取 RSS 最小的那個(gè)候選點(diǎn)作為門限估計(jì)值。這一步在 R 的 tsDyn 包和 Python 里都有成熟實(shí)現(xiàn)。但這還不夠因?yàn)槟呐聸]有真正的門限效應(yīng)你也能找到“使 RSS 最小”的點(diǎn)只是它沒有統(tǒng)計(jì)顯著意義。標(biāo)準(zhǔn)的做法是用自舉bootstrap生成零分布在原假設(shè)“無門限效應(yīng)”下模擬一大批序列計(jì)算每個(gè)模擬序列的超檢驗(yàn)統(tǒng)計(jì)量然后看真實(shí)統(tǒng)計(jì)量落在分布的多極端位置算出 p 值。如果 p 值小于 0.05才算有統(tǒng)計(jì)依據(jù)做分段回歸。3. 用 Python 實(shí)現(xiàn)門限自回歸從數(shù)據(jù)檢驗(yàn)到參數(shù)估計(jì)3.1 數(shù)據(jù)準(zhǔn)備與平穩(wěn)性預(yù)檢驗(yàn)沒有一個(gè)可靠的過程就沒有一個(gè)可靠的模型任何門限模型的起點(diǎn)都是數(shù)據(jù)清洗和平穩(wěn)性檢驗(yàn)。不像線性回歸可以容忍一些粗糙的預(yù)處理門限模型對(duì)數(shù)據(jù)的“分層”非常敏感如果序列里混著異常值或趨勢(shì)項(xiàng)門限點(diǎn)會(huì)被強(qiáng)行拉到一個(gè)錯(cuò)誤的位置。先做 ADF 檢驗(yàn)確認(rèn)平穩(wěn)性不平穩(wěn)就差分這步?jīng)]有繞過空間。import numpy as np import pandas as pd from statsmodels.tsa.stattools import adfuller # 加載 jasa_03m 月度數(shù)據(jù)假設(shè)為單列時(shí)間序列 series pd.read_csv(jasa_03m.csv, index_col0, parse_datesTrue).iloc[:, 0] # ADF 檢驗(yàn)p 值小于 0.05 視為平穩(wěn) adf_result adfuller(series.dropna()) print(fADF 統(tǒng)計(jì)量: {adf_result[0]:.4f}, p 值: {adf_result[1]:.4f}) # 若 p 值偏大做一階差分后重新檢驗(yàn) if adf_result[1] 0.05: series_diff series.diff().dropna() print(原序列非平穩(wěn)使用一階差分序列) else: series_diff series代碼的邏輯很直觀先用adfuller做 ADF 單位根檢驗(yàn)p 值小說明沒有單位根、序列平穩(wěn)p 值大就做一階差分。這里有個(gè)參數(shù)注意點(diǎn)adfuller默認(rèn)的回歸項(xiàng)含常數(shù)項(xiàng)c如果序列均值明顯非零這個(gè)設(shè)置是對(duì)的如果序列圍繞 0 波動(dòng)可以傳入regressionn去掉常數(shù)項(xiàng)能提高檢驗(yàn)功效。還有一點(diǎn)容易忽略——月度數(shù)據(jù)應(yīng)該把a(bǔ)utolagAIC設(shè)為默認(rèn)的自動(dòng)滯后選擇它會(huì)按信息準(zhǔn)則挑最優(yōu)滯后階數(shù)比固定滯后更穩(wěn)。3.2 門限候選點(diǎn)的網(wǎng)格搜索把連續(xù)閾值變成可計(jì)算的離散候選門限值本身是連續(xù)的但實(shí)際計(jì)算時(shí)只能從觀測(cè)值里選。常見做法是把門限變量 z_t 從小到大排序然后按一定的分位數(shù)區(qū)間比如 15% 到 85%作為候選范圍去掉兩端太少的樣本保證每個(gè)體制內(nèi)的觀測(cè)數(shù)足夠做回歸。接下來對(duì)每個(gè)候選點(diǎn)把全樣本切兩半分別估計(jì) AR 模型累加兩個(gè)子模型的殘差平方和。import itertools from statsmodels.tsa.ar_model import AutoReg def fit_tar_rss(series, delay, threshold_candidates, ar_order): 遍歷候選門限點(diǎn)返回 RSS 最小的門限值及對(duì)應(yīng)模型信息 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() # 每個(gè)體制內(nèi)要求最小樣本量否則跳過 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 期為門限變量候選門限設(shè)為序列 15%-85% 分位區(qū)間內(nèi)的觀測(cè)值 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最優(yōu)門限值: {result[threshold]:.4f}, 最小 RSS: {result[rss]:.4f})這段代碼的核心是fit_tar_rss函數(shù)里的雙重循環(huán)邏輯外層遍歷門限候選點(diǎn)內(nèi)層對(duì)分段后的兩個(gè)子樣本分別擬合AutoReg模型。model_low.ssr和model_high.ssr分別是低體制和高體制模型的殘差平方和加起來就是總 RSS。延遲參數(shù)delay決定門限變量用的是哪一期滯后——這個(gè)參數(shù)直接影響模型行為后面會(huì)專門講怎么選。我在項(xiàng)目里的經(jīng)驗(yàn)是候選門限別用手工指定用我這段的分位數(shù)切片方式更客觀5 個(gè)分位點(diǎn)做粗篩找到最優(yōu)區(qū)間后在區(qū)間內(nèi)細(xì)化網(wǎng)格重跑一遍能有效避免漏掉真正的門限值。3.3 固定滯后階數(shù)與門限延遲用 AIC/BIC 在模型空間里選擇門限值定下來之后還需要確定兩件事每個(gè)體制內(nèi)的 AR 滯后階數(shù) p以及門限延遲 d如果用序列自身做門限變量就是 z_t y_{t-d} 里的 d。這兩組參數(shù)的最佳組合通?？啃畔?zhǔn)則搜索確定。常見的搜索范圍是 p ∈ 1~5d ∈ 1~3組合數(shù)量不大暴力枚舉就行。import warnings warnings.filterwarnings(ignore) def select_tar_order(series, p_range, d_range, threshold_candidates, min_samples20): 網(wǎng)格搜索使 AIC 最小的 (p, d) 組合 results [] for p, d in itertools.product(p_range, d_range): # 門限變量使用 y_{t-d}序列做回歸時(shí)使用 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最優(yōu) 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 模型本身分段估計(jì)參數(shù)已經(jīng)比線性模型多了一倍再用 BIC 這種偏保守的準(zhǔn)則容易選出過簡(jiǎn)模型AIC 在樣本量中等的情況下更兼顧擬合與復(fù)雜度。不過如果你的樣本量很大超過 500 個(gè)觀測(cè)換成 BIC 問題也不大最終效果差異很小。還有個(gè)細(xì)節(jié)AutoReg的ssr屬性和aic屬性直接獲取殘差平方和與 Akaike 信息準(zhǔn)則值省去手寫公式的麻煩。跑完網(wǎng)格后注意檢查最優(yōu)組合和第二優(yōu)組合的 AIC 差距如果差距很小小于 2說明模型選擇不夠穩(wěn)健需要檢查數(shù)據(jù)里有沒有異常值擾動(dòng)。4. jasa_03m 月度數(shù)據(jù)的實(shí)戰(zhàn)建模從檢驗(yàn)到預(yù)測(cè)的完整流程4.1 為什么月度數(shù)據(jù)特別適合門限模型我拿 jasa_03m 這批月度數(shù)據(jù)舉例因?yàn)樗浅Ｓ写硇?。月度觀測(cè)天然具有三種結(jié)構(gòu)第一是季節(jié)性周期比如 12 個(gè)月的景氣循環(huán)第二是趨勢(shì)變化比如經(jīng)濟(jì)擴(kuò)張期和收縮期第三是變點(diǎn)事件比如疫情沖擊、政策調(diào)整、供需關(guān)系的結(jié)構(gòu)性轉(zhuǎn)變。這三者疊加在一個(gè)線性模型里往往互相掩蓋——季節(jié)項(xiàng)把突變點(diǎn)抹平趨勢(shì)項(xiàng)把短期的門限效應(yīng)吞掉。門限模型處理月度數(shù)據(jù)的優(yōu)勢(shì)在于它能自動(dòng)根據(jù)門限變量可以是滯后值也可以是外部月度指標(biāo)把數(shù)據(jù)分成“高體制”和“低體制”比如物價(jià)高企的月份和物價(jià)低迷的月份兩個(gè)體制內(nèi)部分別擬合不同的自回歸結(jié)構(gòu)捕捉到的依賴關(guān)系就完全不一樣。你可能已經(jīng)注意到這有點(diǎn)像把“時(shí)間分段”換成“狀態(tài)分段”——不是按時(shí)間軸切一刀而是按變量取值切一刀這樣即便相同月份處于不同年份只要狀態(tài)相似就會(huì)被分到同一個(gè)體制里估計(jì)參數(shù)信息的利用率高得多。用 jasa_03m 這類數(shù)據(jù)做實(shí)證時(shí)我通常先把序列拆成訓(xùn)練集和測(cè)試集比如前 80% 訓(xùn)練后 20% 驗(yàn)證在訓(xùn)練集上完成全部參數(shù)選擇測(cè)試集只用于最終的滾動(dòng)預(yù)測(cè)評(píng)估。4.2 帶外生門限變量的 TAR把外部驅(qū)動(dòng)因素引入模型做月度數(shù)據(jù)時(shí)門限變量只選序列自身的滯后值有時(shí)不夠。比如你預(yù)測(cè)一個(gè)城市的月度用電量門限變量用“上月氣溫是否超過某個(gè)閾值”比用“上個(gè)月用電量是否超過某個(gè)值”更貼近物理實(shí)際。這就回到 2.2 節(jié)講的 TAR 和 SETAR 的分野。帶外生門限變量的實(shí)現(xiàn)方式與純 SETAR 幾乎一樣只是門限變量的數(shù)據(jù)來源變了。def fit_tar_exog(series, exog_threshold, threshold_candidates, ar_order): 外生門限變量的 TAR 擬合 exog_threshold: 與 series 等長(zhǎng)的外生門限變量如氣溫、利率 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一比唯一的區(qū)別是第四行的exog_threshold替代了原有的series.shift(delay)其他一切照舊?？雌饋砀膭?dòng)不大但這個(gè)選擇的背后是建模思路的變化當(dāng)門限變量來自外部時(shí)你得先確認(rèn)它的“時(shí)間對(duì)齊”——是同步值當(dāng)月氣溫還是滯后值上月氣溫這取決于業(yè)務(wù)上哪個(gè)變量“驅(qū)動(dòng)”了序列切換機(jī)制。做外生門限變量時(shí)最容易翻車的就是對(duì)齊問題門限變量和序列的索引對(duì)應(yīng)錯(cuò)位一位整個(gè)模型就串味了。我在項(xiàng)目里會(huì)先用print(pd.DataFrame({series: series, gate: exog_threshold}).dropna().head())檢查對(duì)齊再做擬合這一步五分鐘能省三小時(shí)的排查功夫。4.3 模型預(yù)測(cè)與雙體制的樣本外對(duì)比模型估計(jì)完最終要回答的問題是它比單一線性 AR 模型強(qiáng)在哪兒驗(yàn)證方法不復(fù)雜——用訓(xùn)練好的門限模型對(duì)測(cè)試集逐期做預(yù)測(cè)并和線性 AR 模型做誤差對(duì)比。門限模型的預(yù)測(cè)邏輯是每一步預(yù)測(cè)前先看門限變量落在哪個(gè)體制再用對(duì)應(yīng)體制的 AR 模型來預(yù)測(cè)。預(yù)測(cè)值可以直接用AutoReg的predict方法但門限模型的預(yù)測(cè)要小心滯后階數(shù)的對(duì)齊推薦把預(yù)測(cè)函數(shù)封裝好。def tar_forecast(series, exog_threshold, th, model_low, model_high, ar_order, horizon): 門限模型滾動(dòng)預(yù)測(cè), h 步預(yù)測(cè) forecasts [] history list(series.values) gate_history list(exog_threshold.values) for _ in range(horizon): gate_val gate_history[-1] # 當(dāng)前狀態(tài)決定用哪個(gè)體制的模型預(yù)測(cè) active_model model_low if gate_val th else model_high # 取最近 ar_order 個(gè)觀測(cè)作為起點(diǎn) recent history[-ar_order:] # 手動(dòng)構(gòu)造滯后特征 X np.array([recent[::-1]]).reshape(1, -1) yhat active_model.predict(X)[0] if hasattr(active_model, predict) else np.nan # 手動(dòng)用參數(shù)算滯后系數(shù) 常數(shù)項(xiàng) if np.isnan(yhat): yhat active_model.params.iloc[0] # 常數(shù) 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)是理論化寫法實(shí)際用statsmodels的AutoReg對(duì)象做單步預(yù)測(cè)時(shí)需要先擬合一個(gè)AutoReg的預(yù)測(cè)包裝器或者直接讀取參數(shù)做滾動(dòng)推算。我更推薦后一種方式把a(bǔ)ctive_model.params按常數(shù)和滯后系數(shù)拆開手動(dòng)加權(quán)計(jì)算下一期預(yù)測(cè)值邏輯透明、不容易因?yàn)闀r(shí)序索引錯(cuò)位而出錯(cuò)。horizon是預(yù)測(cè)步長(zhǎng)月度數(shù)據(jù)一般做 3~12 步預(yù)測(cè)。預(yù)測(cè)結(jié)束后對(duì)比mean_absolute_error(forecasts, actual)和線性 AR 的同一指標(biāo)如果門限模型沒有顯著優(yōu)勢(shì)就回頭檢查門限值是不是不夠顯著或者門限體制的樣本量差異太大——這些在下一章詳細(xì)說。5. 門限自回歸避坑指南五個(gè)最容易翻車的地雷5.1 現(xiàn)象門限效應(yīng)檢驗(yàn)的 p 值不穩(wěn)定換個(gè)樣本區(qū)間結(jié)果就變?cè)蜷T限效應(yīng)檢驗(yàn)依賴自舉抽樣的隨機(jī)性不同的隨機(jī)種子會(huì)產(chǎn)生不同的模擬分布p 值落在 0.04 到 0.06 之間浮動(dòng)時(shí)結(jié)論就會(huì)搖擺。另一個(gè)原因是樣本量不夠自舉分布的形狀不夠穩(wěn)定。解決跑檢驗(yàn)時(shí)固定隨機(jī)種子np.random.seed(42)并且至少跑三次不同的種子觀察 p 值是否穩(wěn)定地在 0.05 以下。如果 p 值在臨界點(diǎn)附近晃動(dòng)別急著建模先用數(shù)據(jù)可視化確認(rèn)門限變量的分布是否存在雙峰或明顯的跳躍結(jié)構(gòu)。我在月度數(shù)據(jù)上踩過這個(gè)坑第一次跑 p 值是 0.03換了個(gè)子樣本變成 0.08最后發(fā)現(xiàn)是樣本里有一個(gè)極端的異常值把門限點(diǎn)拉偏了剔除后 p 值穩(wěn)定在 0.01。5.2 現(xiàn)象模型在訓(xùn)練集上擬合得很好測(cè)試集上一塌糊涂原因這是門限模型最典型的過擬合陷阱——搜索門限值的時(shí)候用的是“殘差平方和最小”這個(gè)目標(biāo)但是候選門限點(diǎn)本身是從觀測(cè)值里挑的你其實(shí)是在用測(cè)試信息反選門限。本質(zhì)上你是在拿數(shù)據(jù)來擬合門限的位置而不是先驗(yàn)地指定它。這會(huì)導(dǎo)致門限點(diǎn)選擇過度優(yōu)化訓(xùn)練集內(nèi)部再分段做回歸自由度大幅上升泛化能力自然大幅下降。解決這一點(diǎn)我算交過學(xué)費(fèi)的。后來固定做法是——門限值只從訓(xùn)練集的前 80% 里選后 20% 完全不參與門限值搜索選完門限后再用全訓(xùn)練集估參數(shù)。這樣至少保證門限點(diǎn)的選擇沒有“偷看”最后一段信息。另外給每個(gè)體制設(shè)定最小觀測(cè)數(shù)我通常設(shè)為總樣本的 15% 以上防止門限切在邊緣導(dǎo)致單體制樣本太少。5.3 現(xiàn)象兩個(gè)體制的滯后階數(shù)都用 AIC 定但模型整體比線性 AR 還差原因每個(gè)體制單獨(dú)用 AIC 選滯后階數(shù)可能在低體制選出 p3高體制選出 p5兩套模型的復(fù)雜度加在一起遠(yuǎn)超數(shù)據(jù)能支撐的信息量造成過擬合。更隱蔽的問題是門限模型每個(gè)體制內(nèi)的觀測(cè)數(shù)少相同階數(shù)的 AR 模型在小樣本下估計(jì)方差更大。解決我一般會(huì)先固定兩個(gè)體制使用相同的滯后階數(shù)再用聯(lián)合 AIC兩個(gè)體制 AIC 之和去選。如果不同體制的滯后階數(shù)差異在業(yè)務(wù)上站不住腳就保持同構(gòu)。還有一種更實(shí)用的策略殘差診斷用 Ljung-Box 檢驗(yàn)如果殘差已經(jīng)近似白噪聲就不要再加滯后階數(shù)了——AIC 在這個(gè)場(chǎng)景下容易鉤住多余的主項(xiàng)。5.4 現(xiàn)象門限值估計(jì)出來正好落在數(shù)據(jù)的中位數(shù)附近但切換邏輯說不通原因這不一定是你選錯(cuò)了但也可能是門限變量選錯(cuò)了。比如用 y_{t-1} 做門限門限值切出來的兩段剛好把數(shù)據(jù)按高低排序切開但這兩段沒有業(yè)務(wù)上的“機(jī)制差異”只是數(shù)學(xué)上的偶然。這種模型的預(yù)測(cè)能力往往很脆弱換一個(gè)數(shù)據(jù)集門限值就漂移。解決回到業(yè)務(wù)層面問一個(gè)問題這個(gè)門限變量的哪個(gè)數(shù)值在你的領(lǐng)域里意味著“機(jī)制切換”比如對(duì)銷量數(shù)據(jù)門限可能是“上期銷量高于歷史均值 1.2 倍”——這是產(chǎn)能瓶頸點(diǎn)對(duì)風(fēng)速數(shù)據(jù)門限可能是“切變層高度超過某個(gè)海拔”——這是風(fēng)切變機(jī)制切換點(diǎn)。如果門限值給不出業(yè)務(wù)解釋就換門限變量別硬撐著用數(shù)學(xué)結(jié)果。5.5 現(xiàn)象帶外生門限變量時(shí)兩種體制的參數(shù)看起來沒有顯著差異原因這說明門限變量選的沒有區(qū)分度——它在兩個(gè)區(qū)間內(nèi)對(duì)應(yīng)的序列行為本質(zhì)上一樣分段只是給模型額外塞了一堆參數(shù)。另一個(gè)常見原因是外生門限變量和序列同期相關(guān)導(dǎo)致門限效應(yīng)被序列自身的慣性吸收門限切換本身沒有增量信息。解決在擬合門限變量之前先分別統(tǒng)計(jì)兩個(gè)體制內(nèi)序列的均值、方差和自相關(guān)系數(shù)。如果三個(gè)統(tǒng)計(jì)量沒有顯著差異就不要做門限模型。可以用簡(jiǎn)單的 t 檢驗(yàn)比較兩個(gè)體制的均值差異是否顯著。我自己有個(gè)鐵律門限模型如果只是讓 AIC 降低了不到 5%這個(gè)模型就沒有實(shí)用價(jià)值不如老老實(shí)實(shí)用線性 ARIMA。6. 殘差自舉法計(jì)算預(yù)測(cè)區(qū)間門限模型的不確定性量化最后一節(jié)是門限模型里最實(shí)用但很多人沒做到位的一環(huán)預(yù)測(cè)區(qū)間。區(qū)間的價(jià)值在于幫你區(qū)分“正常的波動(dòng)”和“異常的偏移”。在月度數(shù)據(jù)上一個(gè)未來 3 個(gè)月的預(yù)測(cè)值如果沒有區(qū)間業(yè)務(wù)方根本不知道數(shù)字是大概率事件還是賭運(yùn)氣。殘差自舉法的操作分四步第一步從門限模型的殘差序列里隨機(jī)抽樣第二步把抽樣殘差疊加到預(yù)測(cè)值上生成一條模擬的未來路徑第三步重復(fù) 1000 次收集所有模擬路徑在每個(gè)未來時(shí)間點(diǎn)的值第四步取 2.5% 和 97.5% 分位數(shù)作為 95% 置信區(qū)間。這個(gè)方法的優(yōu)勢(shì)是不需要假設(shè)殘差服從正態(tài)分布——月度數(shù)據(jù)經(jīng)常有偏態(tài)用正態(tài)假設(shè)會(huì)算出過窄的區(qū)間。def tar_forecast_interval(residuals, series, exog_threshold, th, model_low, model_high, ar_order, horizon, n_sim1000): 殘差自舉法生成預(yù)測(cè)區(qū)間 sim_paths [] for _ in range(n_sim): # 有放回抽樣殘差 boot_resid np.random.choice(residuals, sizehorizon, replaceTrue) # 用不帶噪聲的模型預(yù)測(cè)初始路徑 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代碼里最關(guān)鍵的一行是np.random.choice(residuals, sizehorizon, replaceTrue)它決定了重抽樣的隨機(jī)性來源。如果殘差序列存在自相關(guān)直接重抽樣會(huì)破壞時(shí)序結(jié)構(gòu)此時(shí)應(yīng)該先對(duì)殘差擬合一個(gè) AR(1) 模型再抽取白噪聲部分做自舉。我在月度數(shù)據(jù)上的經(jīng)驗(yàn)是殘差如果有明顯的序列相關(guān)Ljung-Box 檢驗(yàn) p 值小于 0.05先對(duì)殘差做 ARIMA 過濾再對(duì)過濾后的白噪聲自舉區(qū)間會(huì)準(zhǔn)很多。預(yù)測(cè)區(qū)間的寬度是模型不確定性的直觀標(biāo)尺如果區(qū)間寬度超過預(yù)測(cè)值本身的 50%說明模型在這個(gè)數(shù)據(jù)集上的信息量有限后面要考慮引入外部變量或更高的采樣頻率。個(gè)人習(xí)慣是任何門限模型的預(yù)測(cè)報(bào)告里一律把單點(diǎn)預(yù)測(cè)和區(qū)間一起交付。區(qū)間不是錦上添花是判斷模型可靠性的核心依據(jù)——如果一位工程師只給我單點(diǎn)預(yù)測(cè)卻不給區(qū)間我基本不會(huì)采用。希望這些內(nèi)容能幫你在做門限自回歸的路上少走些彎路把這套非線性工具真正用起來。本文還有配套的精品資源點(diǎn)擊獲取 門限自回歸實(shí)戰(zhàn)：Python分段建模時(shí)間序列突變
簡(jiǎn)介門限自回歸TAR模型是在傳統(tǒng)線性AR模型基礎(chǔ)上引入閾值機(jī)制的非線性時(shí)間序列方法適用于刻畫不同狀態(tài)下回歸關(guān)系發(fā)生結(jié)構(gòu)性變化的數(shù)據(jù)。這份MATLAB代碼資源面向時(shí)間序列分析方向的學(xué)習(xí)者與研究人員圍繞TAR模型構(gòu)建與模型選擇提供了一套可運(yùn)行的完整示例。壓縮包共6個(gè)文件含3個(gè)m腳本和3個(gè)txt文件體積僅70KB腳本覆蓋數(shù)據(jù)預(yù)處理、閾值識(shí)別、分段AR擬合、最大似然參數(shù)估計(jì)及LR似然比圖繪制等核心環(huán)節(jié)txt文件提供示例數(shù)據(jù)與使用說明便于對(duì)照運(yùn)行。通過調(diào)試和運(yùn)行這套代碼讀者可掌握TAR模型從建模到診斷的完整流程并學(xué)會(huì)利用LR圖判斷最佳門限數(shù)量提升對(duì)非線性時(shí)間序列的建模能力。目前已有356人學(xué)習(xí)下載適合具備基礎(chǔ)回歸知識(shí)、希望實(shí)踐門限模型的MATLAB用戶。1. 門限自回歸時(shí)間序列回歸里的“變臉”問題終于有解了做時(shí)間序列預(yù)測(cè)的人遲早會(huì)撞上一堵墻模型在某個(gè)時(shí)間點(diǎn)之前擬合得漂亮之后卻集體失效。不是數(shù)據(jù)變臟了而是數(shù)據(jù)背后的“機(jī)制”變了——宏觀數(shù)據(jù)遇到政策拐點(diǎn)風(fēng)速序列遇到季節(jié)轉(zhuǎn)換銷售數(shù)據(jù)遇到大促前后。線性 AR 模型假設(shè)過去的關(guān)系永遠(yuǎn)不變但現(xiàn)實(shí)里很多序列存在“門限效應(yīng)”當(dāng)某個(gè)驅(qū)動(dòng)變量跨過閾值序列的行為就會(huì)切換成另一套邏輯。門限自回歸Threshold Autoregressive Model簡(jiǎn)稱 TAR就是專門為這種“變關(guān)系”設(shè)計(jì)的模型。它不強(qiáng)行用一條曲線去擬合所有時(shí)期而是把序列按門限變量的取值切成若干區(qū)間每個(gè)區(qū)間內(nèi)分別擬合自回歸方程。jasa_03m 這類月度數(shù)據(jù)正是 TAR 模型最典型的用武之地月度觀測(cè)里天然藏著季節(jié)、周期和突變點(diǎn)用門限模型能比常規(guī) ARIMA 多抓住一層結(jié)構(gòu)。本文不繞彎子直接講清楚門限自回歸的建模流程、Python 實(shí)現(xiàn)和參數(shù)設(shè)置給你一份能照著復(fù)現(xiàn)的實(shí)戰(zhàn)路徑。2. 門限自回歸的核心邏輯為什么要“分段”做回歸2.1 線性 AR 的假設(shè)盲區(qū)關(guān)系不是恒定的先看一個(gè)最基礎(chǔ)的問題門限自回歸到底在解決什么。經(jīng)典的 AR(p) 模型寫成 y_t c φ?y_{t-1} … φ_p y_{t-p} ε_(tái)t它隱含的假設(shè)是從第一期到最后一期φ 這些系數(shù)是不變的。換句話說無論序列處于高位、低位還是中間狀態(tài)過去值對(duì)當(dāng)前值的影響力度都一樣。這個(gè)假設(shè)在平穩(wěn)序列上還能應(yīng)付一旦序列存在結(jié)構(gòu)性變化——比如某個(gè)月出臺(tái)新政策、市場(chǎng)需求突然換擋——單一組的 φ 就完全不夠用了。門限自回歸的破局思路很直接把序列按某個(gè)變量的數(shù)值劃分成不同的“體制”regime每個(gè)體制內(nèi)各跑各的自回歸。最經(jīng)典的二體制 TAR 模型寫成y_t (φ?? φ??y_{t-1} … φ??y_{t-p}) · I(z_t ≤ c) (φ?? φ??y_{t-1} … φ??y_{t-p}) · I(z_t c) ε_(tái)t其中 z_t 是門限變量c 是門限值I(·) 是指示函數(shù)。公式看著復(fù)雜意思是當(dāng) z_t 小于等于 c 時(shí)序列服從第一套自回歸方程當(dāng) z_t 大于 c 時(shí)切換到第二套方程。z_t 可以是滯后值 y_{t-d}這叫 SETAR自激勵(lì)門限自回歸也可以是外生變量這叫 TAR。這個(gè)分段機(jī)制讓模型天然具備“識(shí)別變結(jié)構(gòu)”的能力。2.2 TAR 與 SETAR 怎么選門限變量是核心分水嶺建模前必須先想清楚你的門限變量是什么。這個(gè)問題決定了模型的類別。如果門限變量是序列自身的滯后值比如 z_t y_{t-1}模型就叫 SETAR如果門限變量是另一個(gè)外生變量比如利率、匯率、氣溫模型就是 TAR。兩者的共性和差異可以用一個(gè)判斷框架說清楚。維度SETARTAR門限變量序列自身滯后值 y_{t-d}外生變量 z_t適用場(chǎng)景序列有自我強(qiáng)化的周期或波動(dòng)聚集序列受外部變量驅(qū)動(dòng)切換機(jī)制建模復(fù)雜度較低數(shù)據(jù)本身夠用需要配套的門限變量數(shù)據(jù)典型應(yīng)用太陽(yáng)黑子、氣溫、股價(jià)波動(dòng)銷售額受政策利率影響、風(fēng)速受氣壓影響我個(gè)人的習(xí)慣是先跑一遍 SETAR 做摸底因?yàn)椴恍枰~外找數(shù)據(jù)如果檢驗(yàn)出明顯的門限效應(yīng)但不穩(wěn)定再考慮換外生門限變量。見過不少項(xiàng)目是數(shù)據(jù)明明有突變點(diǎn)卻硬套 ARIMA 硬扛改造成門限模型后誤差直接降了一截。選擇 TAR 還是 SETAR本質(zhì)是在問你更相信“過去的狀態(tài)”還是“外部的力量”能解釋當(dāng)前的行為切換。2.3 門限效應(yīng)的檢驗(yàn)不要上來就分段這里有個(gè)新手常犯的錯(cuò)誤看到序列圖上有幾個(gè)折點(diǎn)就急著定門限值。門限模型的第一個(gè)正式步驟應(yīng)該是假說檢驗(yàn)——先驗(yàn)證“是否存在門限效應(yīng)”再估計(jì)門限值。統(tǒng)計(jì)上一般用 Chan (1993) 的超檢驗(yàn)過程先按門限變量的取值排序?qū)γ恳粋€(gè)可能作為門限的候選點(diǎn)把樣本分成兩段分別估計(jì)兩段的自回歸模型算殘差平方和RSS取 RSS 最小的那個(gè)候選點(diǎn)作為門限估計(jì)值。這一步在 R 的 tsDyn 包和 Python 里都有成熟實(shí)現(xiàn)。但這還不夠因?yàn)槟呐聸]有真正的門限效應(yīng)你也能找到“使 RSS 最小”的點(diǎn)只是它沒有統(tǒng)計(jì)顯著意義。標(biāo)準(zhǔn)的做法是用自舉bootstrap生成零分布在原假設(shè)“無門限效應(yīng)”下模擬一大批序列計(jì)算每個(gè)模擬序列的超檢驗(yàn)統(tǒng)計(jì)量然后看真實(shí)統(tǒng)計(jì)量落在分布的多極端位置算出 p 值。如果 p 值小于 0.05才算有統(tǒng)計(jì)依據(jù)做分段回歸。3. 用 Python 實(shí)現(xiàn)門限自回歸從數(shù)據(jù)檢驗(yàn)到參數(shù)估計(jì)3.1 數(shù)據(jù)準(zhǔn)備與平穩(wěn)性預(yù)檢驗(yàn)沒有一個(gè)可靠的過程就沒有一個(gè)可靠的模型任何門限模型的起點(diǎn)都是數(shù)據(jù)清洗和平穩(wěn)性檢驗(yàn)。不像線性回歸可以容忍一些粗糙的預(yù)處理門限模型對(duì)數(shù)據(jù)的“分層”非常敏感如果序列里混著異常值或趨勢(shì)項(xiàng)門限點(diǎn)會(huì)被強(qiáng)行拉到一個(gè)錯(cuò)誤的位置。先做 ADF 檢驗(yàn)確認(rèn)平穩(wěn)性不平穩(wěn)就差分這步?jīng)]有繞過空間。import numpy as np import pandas as pd from statsmodels.tsa.stattools import adfuller # 加載 jasa_03m 月度數(shù)據(jù)假設(shè)為單列時(shí)間序列 series pd.read_csv(jasa_03m.csv, index_col0, parse_datesTrue).iloc[:, 0] # ADF 檢驗(yàn)p 值小于 0.05 視為平穩(wěn) adf_result adfuller(series.dropna()) print(fADF 統(tǒng)計(jì)量: {adf_result[0]:.4f}, p 值: {adf_result[1]:.4f}) # 若 p 值偏大做一階差分后重新檢驗(yàn) if adf_result[1] 0.05: series_diff series.diff().dropna() print(原序列非平穩(wěn)使用一階差分序列) else: series_diff series代碼的邏輯很直觀先用adfuller做 ADF 單位根檢驗(yàn)p 值小說明沒有單位根、序列平穩(wěn)p 值大就做一階差分。這里有個(gè)參數(shù)注意點(diǎn)adfuller默認(rèn)的回歸項(xiàng)含常數(shù)項(xiàng)c如果序列均值明顯非零這個(gè)設(shè)置是對(duì)的如果序列圍繞 0 波動(dòng)可以傳入regressionn去掉常數(shù)項(xiàng)能提高檢驗(yàn)功效。還有一點(diǎn)容易忽略——月度數(shù)據(jù)應(yīng)該把a(bǔ)utolagAIC設(shè)為默認(rèn)的自動(dòng)滯后選擇它會(huì)按信息準(zhǔn)則挑最優(yōu)滯后階數(shù)比固定滯后更穩(wěn)。3.2 門限候選點(diǎn)的網(wǎng)格搜索把連續(xù)閾值變成可計(jì)算的離散候選門限值本身是連續(xù)的但實(shí)際計(jì)算時(shí)只能從觀測(cè)值里選。常見做法是把門限變量 z_t 從小到大排序然后按一定的分位數(shù)區(qū)間比如 15% 到 85%作為候選范圍去掉兩端太少的樣本保證每個(gè)體制內(nèi)的觀測(cè)數(shù)足夠做回歸。接下來對(duì)每個(gè)候選點(diǎn)把全樣本切兩半分別估計(jì) AR 模型累加兩個(gè)子模型的殘差平方和。import itertools from statsmodels.tsa.ar_model import AutoReg def fit_tar_rss(series, delay, threshold_candidates, ar_order): 遍歷候選門限點(diǎn)返回 RSS 最小的門限值及對(duì)應(yīng)模型信息 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() # 每個(gè)體制內(nèi)要求最小樣本量否則跳過 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 期為門限變量候選門限設(shè)為序列 15%-85% 分位區(qū)間內(nèi)的觀測(cè)值 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最優(yōu)門限值: {result[threshold]:.4f}, 最小 RSS: {result[rss]:.4f})這段代碼的核心是fit_tar_rss函數(shù)里的雙重循環(huán)邏輯外層遍歷門限候選點(diǎn)內(nèi)層對(duì)分段后的兩個(gè)子樣本分別擬合AutoReg模型。model_low.ssr和model_high.ssr分別是低體制和高體制模型的殘差平方和加起來就是總 RSS。延遲參數(shù)delay決定門限變量用的是哪一期滯后——這個(gè)參數(shù)直接影響模型行為后面會(huì)專門講怎么選。我在項(xiàng)目里的經(jīng)驗(yàn)是候選門限別用手工指定用我這段的分位數(shù)切片方式更客觀5 個(gè)分位點(diǎn)做粗篩找到最優(yōu)區(qū)間后在區(qū)間內(nèi)細(xì)化網(wǎng)格重跑一遍能有效避免漏掉真正的門限值。3.3 固定滯后階數(shù)與門限延遲用 AIC/BIC 在模型空間里選擇門限值定下來之后還需要確定兩件事每個(gè)體制內(nèi)的 AR 滯后階數(shù) p以及門限延遲 d如果用序列自身做門限變量就是 z_t y_{t-d} 里的 d。這兩組參數(shù)的最佳組合通?？啃畔?zhǔn)則搜索確定。常見的搜索范圍是 p ∈ 1~5d ∈ 1~3組合數(shù)量不大暴力枚舉就行。import warnings warnings.filterwarnings(ignore) def select_tar_order(series, p_range, d_range, threshold_candidates, min_samples20): 網(wǎng)格搜索使 AIC 最小的 (p, d) 組合 results [] for p, d in itertools.product(p_range, d_range): # 門限變量使用 y_{t-d}序列做回歸時(shí)使用 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最優(yōu) 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 模型本身分段估計(jì)參數(shù)已經(jīng)比線性模型多了一倍再用 BIC 這種偏保守的準(zhǔn)則容易選出過簡(jiǎn)模型AIC 在樣本量中等的情況下更兼顧擬合與復(fù)雜度。不過如果你的樣本量很大超過 500 個(gè)觀測(cè)換成 BIC 問題也不大最終效果差異很小。還有個(gè)細(xì)節(jié)AutoReg的ssr屬性和aic屬性直接獲取殘差平方和與 Akaike 信息準(zhǔn)則值省去手寫公式的麻煩。跑完網(wǎng)格后注意檢查最優(yōu)組合和第二優(yōu)組合的 AIC 差距如果差距很小小于 2說明模型選擇不夠穩(wěn)健需要檢查數(shù)據(jù)里有沒有異常值擾動(dòng)。4. jasa_03m 月度數(shù)據(jù)的實(shí)戰(zhàn)建模從檢驗(yàn)到預(yù)測(cè)的完整流程4.1 為什么月度數(shù)據(jù)特別適合門限模型我拿 jasa_03m 這批月度數(shù)據(jù)舉例因?yàn)樗浅Ｓ写硇?。月度觀測(cè)天然具有三種結(jié)構(gòu)第一是季節(jié)性周期比如 12 個(gè)月的景氣循環(huán)第二是趨勢(shì)變化比如經(jīng)濟(jì)擴(kuò)張期和收縮期第三是變點(diǎn)事件比如疫情沖擊、政策調(diào)整、供需關(guān)系的結(jié)構(gòu)性轉(zhuǎn)變。這三者疊加在一個(gè)線性模型里往往互相掩蓋——季節(jié)項(xiàng)把突變點(diǎn)抹平趨勢(shì)項(xiàng)把短期的門限效應(yīng)吞掉。門限模型處理月度數(shù)據(jù)的優(yōu)勢(shì)在于它能自動(dòng)根據(jù)門限變量可以是滯后值也可以是外部月度指標(biāo)把數(shù)據(jù)分成“高體制”和“低體制”比如物價(jià)高企的月份和物價(jià)低迷的月份兩個(gè)體制內(nèi)部分別擬合不同的自回歸結(jié)構(gòu)捕捉到的依賴關(guān)系就完全不一樣。你可能已經(jīng)注意到這有點(diǎn)像把“時(shí)間分段”換成“狀態(tài)分段”——不是按時(shí)間軸切一刀而是按變量取值切一刀這樣即便相同月份處于不同年份只要狀態(tài)相似就會(huì)被分到同一個(gè)體制里估計(jì)參數(shù)信息的利用率高得多。用 jasa_03m 這類數(shù)據(jù)做實(shí)證時(shí)我通常先把序列拆成訓(xùn)練集和測(cè)試集比如前 80% 訓(xùn)練后 20% 驗(yàn)證在訓(xùn)練集上完成全部參數(shù)選擇測(cè)試集只用于最終的滾動(dòng)預(yù)測(cè)評(píng)估。4.2 帶外生門限變量的 TAR把外部驅(qū)動(dòng)因素引入模型做月度數(shù)據(jù)時(shí)門限變量只選序列自身的滯后值有時(shí)不夠。比如你預(yù)測(cè)一個(gè)城市的月度用電量門限變量用“上月氣溫是否超過某個(gè)閾值”比用“上個(gè)月用電量是否超過某個(gè)值”更貼近物理實(shí)際。這就回到 2.2 節(jié)講的 TAR 和 SETAR 的分野。帶外生門限變量的實(shí)現(xiàn)方式與純 SETAR 幾乎一樣只是門限變量的數(shù)據(jù)來源變了。def fit_tar_exog(series, exog_threshold, threshold_candidates, ar_order): 外生門限變量的 TAR 擬合 exog_threshold: 與 series 等長(zhǎng)的外生門限變量如氣溫、利率 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一比唯一的區(qū)別是第四行的exog_threshold替代了原有的series.shift(delay)其他一切照舊?？雌饋砀膭?dòng)不大但這個(gè)選擇的背后是建模思路的變化當(dāng)門限變量來自外部時(shí)你得先確認(rèn)它的“時(shí)間對(duì)齊”——是同步值當(dāng)月氣溫還是滯后值上月氣溫這取決于業(yè)務(wù)上哪個(gè)變量“驅(qū)動(dòng)”了序列切換機(jī)制。做外生門限變量時(shí)最容易翻車的就是對(duì)齊問題門限變量和序列的索引對(duì)應(yīng)錯(cuò)位一位整個(gè)模型就串味了。我在項(xiàng)目里會(huì)先用print(pd.DataFrame({series: series, gate: exog_threshold}).dropna().head())檢查對(duì)齊再做擬合這一步五分鐘能省三小時(shí)的排查功夫。4.3 模型預(yù)測(cè)與雙體制的樣本外對(duì)比模型估計(jì)完最終要回答的問題是它比單一線性 AR 模型強(qiáng)在哪兒驗(yàn)證方法不復(fù)雜——用訓(xùn)練好的門限模型對(duì)測(cè)試集逐期做預(yù)測(cè)并和線性 AR 模型做誤差對(duì)比。門限模型的預(yù)測(cè)邏輯是每一步預(yù)測(cè)前先看門限變量落在哪個(gè)體制再用對(duì)應(yīng)體制的 AR 模型來預(yù)測(cè)。預(yù)測(cè)值可以直接用AutoReg的predict方法但門限模型的預(yù)測(cè)要小心滯后階數(shù)的對(duì)齊推薦把預(yù)測(cè)函數(shù)封裝好。def tar_forecast(series, exog_threshold, th, model_low, model_high, ar_order, horizon): 門限模型滾動(dòng)預(yù)測(cè), h 步預(yù)測(cè) forecasts [] history list(series.values) gate_history list(exog_threshold.values) for _ in range(horizon): gate_val gate_history[-1] # 當(dāng)前狀態(tài)決定用哪個(gè)體制的模型預(yù)測(cè) active_model model_low if gate_val th else model_high # 取最近 ar_order 個(gè)觀測(cè)作為起點(diǎn) recent history[-ar_order:] # 手動(dòng)構(gòu)造滯后特征 X np.array([recent[::-1]]).reshape(1, -1) yhat active_model.predict(X)[0] if hasattr(active_model, predict) else np.nan # 手動(dòng)用參數(shù)算滯后系數(shù) 常數(shù)項(xiàng) if np.isnan(yhat): yhat active_model.params.iloc[0] # 常數(shù) 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)是理論化寫法實(shí)際用statsmodels的AutoReg對(duì)象做單步預(yù)測(cè)時(shí)需要先擬合一個(gè)AutoReg的預(yù)測(cè)包裝器或者直接讀取參數(shù)做滾動(dòng)推算。我更推薦后一種方式把a(bǔ)ctive_model.params按常數(shù)和滯后系數(shù)拆開手動(dòng)加權(quán)計(jì)算下一期預(yù)測(cè)值邏輯透明、不容易因?yàn)闀r(shí)序索引錯(cuò)位而出錯(cuò)。horizon是預(yù)測(cè)步長(zhǎng)月度數(shù)據(jù)一般做 3~12 步預(yù)測(cè)。預(yù)測(cè)結(jié)束后對(duì)比mean_absolute_error(forecasts, actual)和線性 AR 的同一指標(biāo)如果門限模型沒有顯著優(yōu)勢(shì)就回頭檢查門限值是不是不夠顯著或者門限體制的樣本量差異太大——這些在下一章詳細(xì)說。5. 門限自回歸避坑指南五個(gè)最容易翻車的地雷5.1 現(xiàn)象門限效應(yīng)檢驗(yàn)的 p 值不穩(wěn)定換個(gè)樣本區(qū)間結(jié)果就變?cè)蜷T限效應(yīng)檢驗(yàn)依賴自舉抽樣的隨機(jī)性不同的隨機(jī)種子會(huì)產(chǎn)生不同的模擬分布p 值落在 0.04 到 0.06 之間浮動(dòng)時(shí)結(jié)論就會(huì)搖擺。另一個(gè)原因是樣本量不夠自舉分布的形狀不夠穩(wěn)定。解決跑檢驗(yàn)時(shí)固定隨機(jī)種子np.random.seed(42)并且至少跑三次不同的種子觀察 p 值是否穩(wěn)定地在 0.05 以下。如果 p 值在臨界點(diǎn)附近晃動(dòng)別急著建模先用數(shù)據(jù)可視化確認(rèn)門限變量的分布是否存在雙峰或明顯的跳躍結(jié)構(gòu)。我在月度數(shù)據(jù)上踩過這個(gè)坑第一次跑 p 值是 0.03換了個(gè)子樣本變成 0.08最后發(fā)現(xiàn)是樣本里有一個(gè)極端的異常值把門限點(diǎn)拉偏了剔除后 p 值穩(wěn)定在 0.01。5.2 現(xiàn)象模型在訓(xùn)練集上擬合得很好測(cè)試集上一塌糊涂原因這是門限模型最典型的過擬合陷阱——搜索門限值的時(shí)候用的是“殘差平方和最小”這個(gè)目標(biāo)但是候選門限點(diǎn)本身是從觀測(cè)值里挑的你其實(shí)是在用測(cè)試信息反選門限。本質(zhì)上你是在拿數(shù)據(jù)來擬合門限的位置而不是先驗(yàn)地指定它。這會(huì)導(dǎo)致門限點(diǎn)選擇過度優(yōu)化訓(xùn)練集內(nèi)部再分段做回歸自由度大幅上升泛化能力自然大幅下降。解決這一點(diǎn)我算交過學(xué)費(fèi)的。后來固定做法是——門限值只從訓(xùn)練集的前 80% 里選后 20% 完全不參與門限值搜索選完門限后再用全訓(xùn)練集估參數(shù)。這樣至少保證門限點(diǎn)的選擇沒有“偷看”最后一段信息。另外給每個(gè)體制設(shè)定最小觀測(cè)數(shù)我通常設(shè)為總樣本的 15% 以上防止門限切在邊緣導(dǎo)致單體制樣本太少。5.3 現(xiàn)象兩個(gè)體制的滯后階數(shù)都用 AIC 定但模型整體比線性 AR 還差原因每個(gè)體制單獨(dú)用 AIC 選滯后階數(shù)可能在低體制選出 p3高體制選出 p5兩套模型的復(fù)雜度加在一起遠(yuǎn)超數(shù)據(jù)能支撐的信息量造成過擬合。更隱蔽的問題是門限模型每個(gè)體制內(nèi)的觀測(cè)數(shù)少相同階數(shù)的 AR 模型在小樣本下估計(jì)方差更大。解決我一般會(huì)先固定兩個(gè)體制使用相同的滯后階數(shù)再用聯(lián)合 AIC兩個(gè)體制 AIC 之和去選。如果不同體制的滯后階數(shù)差異在業(yè)務(wù)上站不住腳就保持同構(gòu)。還有一種更實(shí)用的策略殘差診斷用 Ljung-Box 檢驗(yàn)如果殘差已經(jīng)近似白噪聲就不要再加滯后階數(shù)了——AIC 在這個(gè)場(chǎng)景下容易鉤住多余的主項(xiàng)。5.4 現(xiàn)象門限值估計(jì)出來正好落在數(shù)據(jù)的中位數(shù)附近但切換邏輯說不通原因這不一定是你選錯(cuò)了但也可能是門限變量選錯(cuò)了。比如用 y_{t-1} 做門限門限值切出來的兩段剛好把數(shù)據(jù)按高低排序切開但這兩段沒有業(yè)務(wù)上的“機(jī)制差異”只是數(shù)學(xué)上的偶然。這種模型的預(yù)測(cè)能力往往很脆弱換一個(gè)數(shù)據(jù)集門限值就漂移。解決回到業(yè)務(wù)層面問一個(gè)問題這個(gè)門限變量的哪個(gè)數(shù)值在你的領(lǐng)域里意味著“機(jī)制切換”比如對(duì)銷量數(shù)據(jù)門限可能是“上期銷量高于歷史均值 1.2 倍”——這是產(chǎn)能瓶頸點(diǎn)對(duì)風(fēng)速數(shù)據(jù)門限可能是“切變層高度超過某個(gè)海拔”——這是風(fēng)切變機(jī)制切換點(diǎn)。如果門限值給不出業(yè)務(wù)解釋就換門限變量別硬撐著用數(shù)學(xué)結(jié)果。5.5 現(xiàn)象帶外生門限變量時(shí)兩種體制的參數(shù)看起來沒有顯著差異原因這說明門限變量選的沒有區(qū)分度——它在兩個(gè)區(qū)間內(nèi)對(duì)應(yīng)的序列行為本質(zhì)上一樣分段只是給模型額外塞了一堆參數(shù)。另一個(gè)常見原因是外生門限變量和序列同期相關(guān)導(dǎo)致門限效應(yīng)被序列自身的慣性吸收門限切換本身沒有增量信息。解決在擬合門限變量之前先分別統(tǒng)計(jì)兩個(gè)體制內(nèi)序列的均值、方差和自相關(guān)系數(shù)。如果三個(gè)統(tǒng)計(jì)量沒有顯著差異就不要做門限模型。可以用簡(jiǎn)單的 t 檢驗(yàn)比較兩個(gè)體制的均值差異是否顯著。我自己有個(gè)鐵律門限模型如果只是讓 AIC 降低了不到 5%這個(gè)模型就沒有實(shí)用價(jià)值不如老老實(shí)實(shí)用線性 ARIMA。6. 殘差自舉法計(jì)算預(yù)測(cè)區(qū)間門限模型的不確定性量化最后一節(jié)是門限模型里最實(shí)用但很多人沒做到位的一環(huán)預(yù)測(cè)區(qū)間。區(qū)間的價(jià)值在于幫你區(qū)分“正常的波動(dòng)”和“異常的偏移”。在月度數(shù)據(jù)上一個(gè)未來 3 個(gè)月的預(yù)測(cè)值如果沒有區(qū)間業(yè)務(wù)方根本不知道數(shù)字是大概率事件還是賭運(yùn)氣。殘差自舉法的操作分四步第一步從門限模型的殘差序列里隨機(jī)抽樣第二步把抽樣殘差疊加到預(yù)測(cè)值上生成一條模擬的未來路徑第三步重復(fù) 1000 次收集所有模擬路徑在每個(gè)未來時(shí)間點(diǎn)的值第四步取 2.5% 和 97.5% 分位數(shù)作為 95% 置信區(qū)間。這個(gè)方法的優(yōu)勢(shì)是不需要假設(shè)殘差服從正態(tài)分布——月度數(shù)據(jù)經(jīng)常有偏態(tài)用正態(tài)假設(shè)會(huì)算出過窄的區(qū)間。def tar_forecast_interval(residuals, series, exog_threshold, th, model_low, model_high, ar_order, horizon, n_sim1000): 殘差自舉法生成預(yù)測(cè)區(qū)間 sim_paths [] for _ in range(n_sim): # 有放回抽樣殘差 boot_resid np.random.choice(residuals, sizehorizon, replaceTrue) # 用不帶噪聲的模型預(yù)測(cè)初始路徑 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代碼里最關(guān)鍵的一行是np.random.choice(residuals, sizehorizon, replaceTrue)它決定了重抽樣的隨機(jī)性來源。如果殘差序列存在自相關(guān)直接重抽樣會(huì)破壞時(shí)序結(jié)構(gòu)此時(shí)應(yīng)該先對(duì)殘差擬合一個(gè) AR(1) 模型再抽取白噪聲部分做自舉。我在月度數(shù)據(jù)上的經(jīng)驗(yàn)是殘差如果有明顯的序列相關(guān)Ljung-Box 檢驗(yàn) p 值小于 0.05先對(duì)殘差做 ARIMA 過濾再對(duì)過濾后的白噪聲自舉區(qū)間會(huì)準(zhǔn)很多。預(yù)測(cè)區(qū)間的寬度是模型不確定性的直觀標(biāo)尺如果區(qū)間寬度超過預(yù)測(cè)值本身的 50%說明模型在這個(gè)數(shù)據(jù)集上的信息量有限后面要考慮引入外部變量或更高的采樣頻率。個(gè)人習(xí)慣是任何門限模型的預(yù)測(cè)報(bào)告里一律把單點(diǎn)預(yù)測(cè)和區(qū)間一起交付。區(qū)間不是錦上添花是判斷模型可靠性的核心依據(jù)——如果一位工程師只給我單點(diǎn)預(yù)測(cè)卻不給區(qū)間我基本不會(huì)采用。希望這些內(nèi)容能幫你在做門限自回歸的路上少走些彎路把這套非線性工具真正用起來。本文還有配套的精品資源點(diǎn)擊獲取