簡(jiǎn)介這份資源面向具備一定R語言基礎(chǔ)、希望深入掌握時(shí)間序列建模的數(shù)據(jù)分析與預(yù)測(cè)從業(yè)者聚焦加法與乘法過程兩類核心思路并延伸至廣義可加模型GAM的非線性建模場(chǎng)景。壓縮包內(nèi)共1個(gè)文件為R腳本整體約1KB可直接在R環(huán)境中加載運(yùn)行便于對(duì)照代碼理解建模流程。內(nèi)容圍繞時(shí)間序列的基本構(gòu)成展開涉及趨勢(shì)、季節(jié)性與隨機(jī)成分的線性相加與相互影響兩種假設(shè)并串聯(lián)ARIMA、SARIMA、STL分解以及基于mgcv包的gam()擬合等實(shí)現(xiàn)路徑同時(shí)涵蓋數(shù)據(jù)加載、模型選擇、參數(shù)調(diào)優(yōu)與殘差診斷等環(huán)節(jié)。目前已有618人學(xué)習(xí)下載適合作為經(jīng)濟(jì)、金融、氣象等領(lǐng)域預(yù)測(cè)任務(wù)的實(shí)操參考幫助讀者把抽象模型轉(zhuǎn)化為可復(fù)用的代碼框架并借助GAM樣條項(xiàng)捕捉復(fù)雜非線性趨勢(shì)。1. 時(shí)間序列的加法與乘法為什么你的預(yù)測(cè)總在趨勢(shì)拐點(diǎn)翻車很多人做時(shí)間序列預(yù)測(cè)時(shí)拿到數(shù)據(jù)第一件事就是往 ARIMA 或 Prophet 里一塞跑出來一條平滑曲線看著擬合得不錯(cuò)一到真實(shí)業(yè)務(wù)里就發(fā)現(xiàn)趨勢(shì)一變預(yù)測(cè)全廢。問題往往不在模型本身而在于你默認(rèn)了一個(gè)假設(shè)——這個(gè)序列是「加法結(jié)構(gòu)」還是「乘法結(jié)構(gòu)」。這兩個(gè)詞聽起來像統(tǒng)計(jì)學(xué)課本里的老古董但在 R 語言里用 GAM廣義可加模型做時(shí)間序列分解和預(yù)測(cè)時(shí)它直接決定了你的模型是「能跟著趨勢(shì)走」還是「一到旺季就低估」。加法模型假設(shè)各成分獨(dú)立疊加觀測(cè)值 趨勢(shì) 季節(jié) 殘差。乘法模型假設(shè)成分之間是比例關(guān)系觀測(cè)值 趨勢(shì) × 季節(jié) × 殘差?，F(xiàn)實(shí)業(yè)務(wù)里銷售額隨趨勢(shì)增長(zhǎng)時(shí)季節(jié)性波動(dòng)幅度也在放大這就是典型的乘法結(jié)構(gòu)。如果你硬用加法模型去擬合殘差會(huì)呈現(xiàn)明顯的異方差預(yù)測(cè)區(qū)間在高峰期窄得離譜低谷期又寬得沒用。R 語言里的 GAM 通過mgcv包提供了靈活的平滑項(xiàng)可以讓你在同一個(gè)框架下處理這兩種結(jié)構(gòu)甚至混合結(jié)構(gòu)。這篇文章面向的是已經(jīng)會(huì)用 R 做基本數(shù)據(jù)分析、但想在時(shí)間序列上把 GAM 用對(duì)的從業(yè)者。我會(huì)把加法與乘法的選擇邏輯、GAM 的平滑項(xiàng)設(shè)置、以及實(shí)際跑通一套分解加預(yù)測(cè)的流程講清楚順帶把那些讓我翻過車的參數(shù)坑標(biāo)出來。2. 加法與乘法結(jié)構(gòu)的判定先看殘差再談模型2.1 從業(yè)務(wù)場(chǎng)景反推結(jié)構(gòu)類型拿到一條時(shí)間序列不要急著畫圖。先問一個(gè)業(yè)務(wù)問題當(dāng)趨勢(shì)上升 10% 時(shí)季節(jié)波動(dòng)的絕對(duì)幅度是保持不變還是也跟著漲 10%如果是前者加法結(jié)構(gòu)更合理如果是后者乘法結(jié)構(gòu)更合理。這個(gè)判斷不需要任何統(tǒng)計(jì)檢驗(yàn)靠業(yè)務(wù)常識(shí)就能定調(diào)。舉個(gè)例子某電商平臺(tái)的日訂單量。過去兩年整體趨勢(shì)從每天 1000 單漲到 3000 單同時(shí)「雙十一」當(dāng)天的峰值從 5000 單漲到 15000 單。峰值與均值的比例大致穩(wěn)定在 5 倍左右這就是乘法結(jié)構(gòu)的典型信號(hào)。反過來如果是一個(gè)成熟產(chǎn)品的日活用戶趨勢(shì)基本平緩周末比工作日固定少 2000 人那加法結(jié)構(gòu)就夠了。在 R 里你可以用forecast包的decompose()函數(shù)分別跑加法和乘法分解然后對(duì)比殘差圖。乘法分解的殘差如果比加法分解更接近白噪聲且殘差絕對(duì)值不隨趨勢(shì)增大而膨脹那就選乘法。library(forecast) # 假設(shè) ts_data 是一個(gè) ts 對(duì)象頻率為 7周季節(jié) # 加法分解 decomp_add - decompose(ts_data, type additive) # 乘法分解 decomp_mult - decompose(ts_data, type multiplicative) # 對(duì)比殘差的標(biāo)準(zhǔn)差隨趨勢(shì)的變化 # 如果乘法分解的殘差更穩(wěn)定選乘法 sd_add - sd(decomp_add$random, na.rm TRUE) sd_mult - sd(decomp_mult$random, na.rm TRUE) # 更直觀的方法看殘差與趨勢(shì)的相關(guān)系數(shù) trend_add - decomp_add$trend resid_add - decomp_add$random cor_add - cor(trend_add, abs(resid_add), use complete.obs) trend_mult - decomp_mult$trend resid_mult - decomp_mult$random cor_mult - cor(trend_mult, abs(resid_mult), use complete.obs) cat(加法殘差與趨勢(shì)相關(guān)系數(shù):, cor_add, \n) cat(乘法殘差與趨勢(shì)相關(guān)系數(shù):, cor_mult, \n)這段代碼的邏輯是如果殘差的絕對(duì)值和趨勢(shì)有正相關(guān)說明波動(dòng)幅度隨趨勢(shì)增長(zhǎng)加法模型就不合適。參數(shù)上type指定分解類型frequency在ts()里設(shè)定。注意decompose()只適合季節(jié)性固定的序列如果季節(jié)模式在變得用stl()或 GAM。2.2 用 GAM 的平滑項(xiàng)同時(shí)捕捉趨勢(shì)與季節(jié)GAM 的優(yōu)勢(shì)在于它不預(yù)設(shè)趨勢(shì)是線性的也不預(yù)設(shè)季節(jié)是固定的。mgcv包里的gam()函數(shù)可以用s()指定平滑項(xiàng)用te()指定交互項(xiàng)。對(duì)于時(shí)間序列我一般這樣建library(mgcv) # 構(gòu)造時(shí)間索引和季節(jié)索引 n - length(ts_data) time_idx - 1:n season_idx - cycle(ts_data) # 1 到 frequency # 加法結(jié)構(gòu) GAM gam_add - gam(ts_data ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian()) # 乘法結(jié)構(gòu)先取對(duì)數(shù)再跑加法 GAM gam_mult - gam(log(ts_data) ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian())這里的關(guān)鍵參數(shù)是k它控制平滑項(xiàng)的最大自由度。k太小會(huì)欠擬合趨勢(shì)被抹平k太大會(huì)過擬合把噪聲當(dāng)信號(hào)。我一般從 10 到 30 之間試用gam.check()看殘差和k的顯著性。bs cc是循環(huán)平滑專門用于季節(jié)索引這種周期性變量不加這個(gè)12 月和 1 月之間會(huì)出現(xiàn)斷點(diǎn)。乘法結(jié)構(gòu)取對(duì)數(shù)后跑加法 GAM等價(jià)于在原始尺度上做乘法分解。預(yù)測(cè)時(shí)記得exp()回來但要注意偏差修正——直接exp()會(huì)低估均值因?yàn)閷?duì)數(shù)正態(tài)分布的均值是exp(mu sigma^2/2)。這個(gè)坑我后面會(huì)細(xì)說。2.3 模型選擇的量化依據(jù)AIC 與殘差診斷加法 GAM 和乘法 GAM 跑完后不能只看圖。AIC()可以比較但注意乘法模型是在對(duì)數(shù)尺度上算的AIC 不能直接和加法模型比。正確做法是把兩個(gè)模型的預(yù)測(cè)值都變換回原始尺度算 RMSE 或 MAE用交叉驗(yàn)證比較。# 樣本外交叉驗(yàn)證留出最后 30 個(gè)點(diǎn) train_n - n - 30 train_data - ts_data[1:train_n] test_data - ts_data[(train_n 1):n] # 重新擬合加法模型 gam_add_train - gam(train_data ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 重新擬合乘法模型 gam_mult_train - gam(log(train_data) ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 預(yù)測(cè) pred_add - predict(gam_add_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult_log - predict(gam_mult_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult - exp(pred_mult_log) # 計(jì)算 RMSE rmse_add - sqrt(mean((test_data - pred_add)^2)) rmse_mult - sqrt(mean((test_data - pred_mult)^2)) cat(加法 RMSE:, rmse_add, \n) cat(乘法 RMSE:, rmse_mult, \n)這段代碼的核心是「用樣本外誤差說話」。參數(shù)上k在訓(xùn)練集上重新選不要用全量數(shù)據(jù)調(diào)好的k直接套。cycle()提取季節(jié)位置predict()的newdata必須包含和訓(xùn)練時(shí)同名的變量。如果乘法模型的 RMSE 明顯更低且殘差沒有異方差那就選乘法。3. 在 R 里跑通 GAM 時(shí)間序列從數(shù)據(jù)到預(yù)測(cè)的完整鏈路3.1 數(shù)據(jù)準(zhǔn)備與 ts 對(duì)象構(gòu)造R 里做時(shí)間序列第一步是把數(shù)據(jù)轉(zhuǎn)成ts對(duì)象。很多人從 CSV 讀進(jìn)來直接跑結(jié)果cycle()返回 NULL季節(jié)平滑項(xiàng)直接報(bào)錯(cuò)。正確做法是明確頻率和起始時(shí)間。# 假設(shè) raw_data 是數(shù)據(jù)框date 列是日期value 列是觀測(cè)值 raw_data$date - as.Date(raw_data$date) # 按日期排序 raw_data - raw_data[order(raw_data$date), ] # 構(gòu)造 ts 對(duì)象頻率為 7周數(shù)據(jù) ts_data - ts(raw_data$value, frequency 7, start c(year(min(raw_data$date)), as.numeric(format(min(raw_data$date), %j)))) # 檢查 cat(頻率:, frequency(ts_data), \n) cat(周期數(shù):, length(ts_data) / frequency(ts_data), \n) head(cycle(ts_data))frequency 7表示每周 7 個(gè)觀測(cè)start參數(shù)指定起始年份和年內(nèi)第幾天。如果頻率設(shè)錯(cuò)比如日數(shù)據(jù)設(shè)成 365cycle()會(huì)返回 1 到 365季節(jié)平滑項(xiàng)k就得設(shè)得很大計(jì)算量爆炸且容易過擬合。常見做法是日數(shù)據(jù)用frequency 7捕捉周內(nèi)模式或者frequency 365.25捕捉年內(nèi)模式但后者需要k至少 50 以上我一般先用周頻率跑通再考慮年頻率。3.2 GAM 平滑項(xiàng)的參數(shù)設(shè)置與 gam.check 解讀gam()函數(shù)里最關(guān)鍵的三個(gè)參數(shù)k、bs、family。k是基函數(shù)維度決定平滑曲線的最大彎曲次數(shù)。bs是基函數(shù)類型tp是薄板回歸樣條默認(rèn)cc是循環(huán)三次樣條用于周期變量cr是三次回歸樣條計(jì)算更快。family指定分布族高斯用于連續(xù)值泊松用于計(jì)數(shù)負(fù)二項(xiàng)用于過離散計(jì)數(shù)。# 完整模型擬合 gam_fit - gam(ts_data ~ s(time_idx, k 25, bs tp) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 REML 估計(jì)平滑參數(shù) # 檢查 gam.check(gam_fit)gam.check()輸出四張圖殘差 vs 擬合值、QQ 圖、直方圖、響應(yīng) vs 擬合值。重點(diǎn)看第一張如果殘差呈現(xiàn)漏斗形隨擬合值增大而擴(kuò)散說明方差非恒定加法高斯模型不合適要么換乘法取對(duì)數(shù)要么換family Gamma。QQ 圖如果尾部偏離說明殘差非正態(tài)預(yù)測(cè)區(qū)間會(huì)不準(zhǔn)。method REML是我強(qiáng)烈建議的默認(rèn)的 GCV 在樣本量小的時(shí)候容易過擬合REML 更穩(wěn)健。這個(gè)參數(shù)不設(shè)gam()會(huì)用 GCV殘差診斷經(jīng)常顯示k不夠其實(shí)換了 REML 就好了。3.3 預(yù)測(cè)與置信區(qū)間對(duì)數(shù)變換后的偏差修正乘法模型預(yù)測(cè)時(shí)predict()返回的是對(duì)數(shù)尺度上的值exp()回去之后得到的是中位數(shù)不是均值。如果業(yè)務(wù)要的是均值預(yù)測(cè)必須加偏差修正。# 乘法模型預(yù)測(cè)對(duì)數(shù)尺度 pred_log - predict(gam_mult, newdata data.frame( time_idx (n 1):(n 30), season_idx rep(1:7, length.out 30) ), se.fit TRUE) # 偏差修正exp(mu sigma^2/2) sigma2 - summary(gam_mult)$scale pred_mean - exp(pred_log$fit sigma2 / 2) # 置信區(qū)間近似 lower - exp(pred_log$fit - 1.96 * pred_log$se.fit) upper - exp(pred_log$fit 1.96 * pred_log$se.fit) # 注意這個(gè)區(qū)間是條件均值區(qū)間不是預(yù)測(cè)區(qū)間 # 預(yù)測(cè)區(qū)間還要加殘差方差 pred_lower - exp(pred_log$fit - 1.96 * sqrt(pred_log$se.fit^2 sigma2)) pred_upper - exp(pred_log$fit 1.96 * sqrt(pred_log$se.fit^2 sigma2))summary(gam_mult)$scale提取的是殘差方差估計(jì)。偏差修正項(xiàng)sigma2 / 2看起來小但當(dāng)sigma2大的時(shí)候不修正會(huì)低估均值 5% 到 10%在庫(kù)存補(bǔ)貨場(chǎng)景里就是真金白銀的差距。置信區(qū)間和預(yù)測(cè)區(qū)間是兩回事前者是均值的區(qū)間后者是單次觀測(cè)的區(qū)間。業(yè)務(wù)上要「明天可能賣多少」用預(yù)測(cè)區(qū)間要「明天平均賣多少」用置信區(qū)間。4. 避坑與排查那些讓我重跑模型的參數(shù)陷阱4.1 現(xiàn)象gam.check 顯示 k 不夠調(diào)大后過擬合原因k的默認(rèn)值是 10對(duì)于趨勢(shì)變化劇烈的序列10 個(gè)基函數(shù)不夠彎曲。但直接調(diào)到 50模型會(huì)把噪聲也擬合進(jìn)去樣本外誤差反而上升。解決先用gam.check()看殘差是否還有模式。如果殘差已經(jīng)像白噪聲k就夠了不用管 p 值。如果殘差還有趨勢(shì)每次加 5直到殘差無模式。同時(shí)用method REML它比 GCV 更不容易過擬合。我一般從k 20開始趨勢(shì)項(xiàng)和季節(jié)項(xiàng)分開調(diào)。4.2 現(xiàn)象乘法模型預(yù)測(cè)值在低谷期為負(fù)原因?qū)?shù)變換后預(yù)測(cè)再exp()回去理論上不會(huì)負(fù)。但如果用了family gaussian()且沒有取對(duì)數(shù)直接跑乘法結(jié)構(gòu)比如把季節(jié)項(xiàng)寫成乘積形式線性預(yù)測(cè)子可能為負(fù)。解決乘法結(jié)構(gòu)必須在對(duì)數(shù)尺度上建模。log(ts_data)后跑加法 GAMexp()回來。如果原始數(shù)據(jù)有零或負(fù)值先加一個(gè)常數(shù)再取對(duì)數(shù)或者改用family Gamma(link log)后者直接在對(duì)數(shù)鏈接函數(shù)上建模不需要手動(dòng)變換。4.3 現(xiàn)象季節(jié)平滑項(xiàng) bs cc 報(bào)錯(cuò)「A term has fewer unique covariate combinations than specified maximum degrees of freedom」原因k設(shè)得比季節(jié)周期的唯一值數(shù)量還大。比如周頻率frequency 7season_idx只有 7 個(gè)唯一值k設(shè)成 10 就報(bào)錯(cuò)。解決k必須小于等于唯一值數(shù)量。周頻率用k 7月頻率用k 12日頻率用k 7周內(nèi)模式或k 365年內(nèi)模式但計(jì)算量大。如果確實(shí)需要更大的k改用bs tp并手動(dòng)構(gòu)造傅里葉項(xiàng)。4.4 現(xiàn)象預(yù)測(cè)區(qū)間在趨勢(shì)上升段明顯偏窄原因GAM 的predict(se.fit TRUE)只給了平滑項(xiàng)的不確定性沒有包含殘差方差。而且如果模型是加法結(jié)構(gòu)但數(shù)據(jù)實(shí)際是乘法結(jié)構(gòu)殘差方差會(huì)隨趨勢(shì)增大區(qū)間估計(jì)系統(tǒng)性偏窄。解決預(yù)測(cè)區(qū)間用sqrt(se.fit^2 sigma2)不要只用se.fit。同時(shí)檢查殘差是否異方差如果是換乘法結(jié)構(gòu)或family Gamma。我習(xí)慣在預(yù)測(cè)后畫殘差 vs 擬合值圖確認(rèn)沒有漏斗形才發(fā)結(jié)果。4.5 現(xiàn)象時(shí)間索引 time_idx 從 1 到 n預(yù)測(cè)時(shí)新數(shù)據(jù)的時(shí)間索引對(duì)不上原因訓(xùn)練時(shí)time_idx是 1 到 n預(yù)測(cè)時(shí)新數(shù)據(jù)的time_idx必須從 n1 開始。如果直接用1:30模型會(huì)以為你在預(yù)測(cè)歷史區(qū)間。解決預(yù)測(cè)時(shí)構(gòu)造newdata的time_idx用(n1):(nh)season_idx用cycle()或手動(dòng)指定。更穩(wěn)妥的做法是把時(shí)間索引存成數(shù)據(jù)框的一列訓(xùn)練和預(yù)測(cè)用同一套構(gòu)造邏輯。5. 進(jìn)階技巧用 te() 捕捉趨勢(shì)與季節(jié)的交互效應(yīng)5.1 什么時(shí)候需要交互項(xiàng)加法 GAM 假設(shè)趨勢(shì)和季節(jié)是獨(dú)立作用的趨勢(shì)上升不影響季節(jié)波動(dòng)的形狀。但現(xiàn)實(shí)中很多序列的季節(jié)模式會(huì)隨趨勢(shì)變化。比如一個(gè)產(chǎn)品剛上市時(shí)只有周末賣得好隨著市場(chǎng)成熟工作日銷量也上來了季節(jié)波動(dòng)的幅度和形狀都變了。這時(shí)候s(time_idx) s(season_idx)就不夠了需要用te(time_idx, season_idx)張量積平滑來捕捉交互。# 張量積交互模型 gam_te - gam(ts_data ~ te(time_idx, season_idx, k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) # 對(duì)比無交互模型 gam_no_te - gam(ts_data ~ s(time_idx, k 15) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 AIC 比較同一尺度下 cat(交互模型 AIC:, AIC(gam_te), \n) cat(無交互 AIC:, AIC(gam_no_te), \n)te()的參數(shù)k是一個(gè)向量分別對(duì)應(yīng)兩個(gè)變量的基函數(shù)維度。bs也是向量對(duì)應(yīng)各自的基函數(shù)類型。交互模型的自由度是k[1] * k[2]的量級(jí)計(jì)算量比加法模型大很多樣本量少于 200 時(shí)慎用容易過擬合。我一般先跑加法模型如果殘差在特定時(shí)間段比如旺季還有系統(tǒng)偏差再上交互項(xiàng)。5.2 用 vis.gam 可視化交互效應(yīng)mgcv自帶的vis.gam()可以畫三維透視圖或等高線圖直觀看到季節(jié)模式如何隨趨勢(shì)變化。# 三維透視圖 vis.gam(gam_te, view c(time_idx, season_idx), theta 30, phi 20, color heat, plot.type persp, ticktype detailed) # 等高線圖 vis.gam(gam_te, view c(time_idx, season_idx), plot.type contour, color heat)view指定兩個(gè)變量theta和phi控制視角plot.type選persp或contour。等高線圖里如果等高線是平行的說明沒有交互如果彎曲或交叉說明季節(jié)模式隨趨勢(shì)變了。這個(gè)圖我每次做完交互模型都會(huì)看比 AIC 更直觀。5.3 一個(gè)我常犯的錯(cuò)誤交互項(xiàng)加了但沒做預(yù)測(cè)對(duì)比剛用te()的時(shí)候我看到 AIC 降了就以為交互模型更好直接拿去預(yù)測(cè)結(jié)果樣本外 RMSE 反而更高。后來養(yǎng)成習(xí)慣不管 AIC 多低一定做樣本外交叉驗(yàn)證。交互模型的方差更大樣本量不夠時(shí)AIC 的懲罰項(xiàng)不足以抵消過擬合。# 樣本外對(duì)比 train_n - n - 30 gam_te_train - gam(ts_data[1:train_n] ~ te(1:train_n, cycle(ts_data[1:train_n]), k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) pred_te - predict(gam_te_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(ts_data[(train_n 1):n]) )) rmse_te - sqrt(mean((ts_data[(train_n 1):n] - pred_te)^2)) cat(交互模型樣本外 RMSE:, rmse_te, \n)如果rmse_te比加法模型的樣本外 RMSE 還大說明交互項(xiàng)在擬合噪聲果斷回退。這個(gè)習(xí)慣讓我少了很多「AIC 好看但業(yè)務(wù)翻車」的情況。5.4 最后的習(xí)慣先畫圖再跑模型我現(xiàn)在拿到任何時(shí)間序列第一件事是plot(ts_data)第二件事是plot(decompose(ts_data))或plot(stl(ts_data, s.window periodic))。圖上看一眼趨勢(shì)和季節(jié)的形態(tài)比任何統(tǒng)計(jì)檢驗(yàn)都快。加法還是乘法很多時(shí)候圖上一眼就能定如果季節(jié)波動(dòng)的幅度隨趨勢(shì)明顯放大乘法如果幅度穩(wěn)定加法。GAM 只是把這個(gè)判斷量化并給出預(yù)測(cè)區(qū)間。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取