原文鏈接：http://tecdat.cn/?p=9024

用GAM進行建模時間序列

我已經(jīng)準備了一個文件，其中包含四個用電時間序列來進行分析。數(shù)據(jù)操作將由data.table程序包完成。

將提及的智能電表數(shù)據(jù)讀到data.table。

DT <- as.data.table(read_feather("DT_4_ind"))

使用GAM回歸模型。將工作日的字符轉(zhuǎn)換為整數(shù)，并使用recode包中的函數(shù)重新編碼工作日：1.星期一，…，7星期日。

1. DT[, week_num := as.integer(car::recode(week,

2. "'Monday'='1';'Tuesday'='2';'Wednesday'='3';'Thursday'='4';

3. 'Friday'='5';'Saturday'='6';'Sunday'='7'"))]

將信息存儲在日期變量中，以簡化工作。

1. n_type <- unique(DT[, type])

2. n_date <- unique(DT[, date])

3. n_weekdays <- unique(DT[, week])

4. period <- 48

讓我們看一下用電量的一些數(shù)據(jù)并對其進行分析。

1. data_r <- DT[(type == n_type[1] & date %in% n_date[57:70])]

2. 


3. ggplot(data_r, aes(date_time, value)) +

4. geom_line() +

5. theme(panel.border = element_blank(),

6. panel.background = element_blank(),

7. panel.grid.minor = element_line(colour = "grey90"),

8. panel.grid.major = element_line(colour = "grey90"),

9. panel.grid.major.x = element_line(colour = "grey90"),

10. axis.text = element_text(size = 10),

11. axis.title = element_text(size = 12, face = "bold")) +

12. labs(x = "Date", y = "Load (kW)")

在繪制的時間序列中可以看到兩個主要的季節(jié)性：每日和每周。我們在一天中有48個測量值，在一周中有7天，因此這將是我們用來對因變量–電力負荷進行建模的自變量。

訓(xùn)練我們的第一個GAM。通過平滑函數(shù)s對自變量建模，對于每日季節(jié)性，使用三次樣條回歸，對于每周季節(jié)性，使用P樣條。

1. gam_1 <- gam(Load ~ s(Daily, bs = "cr", k = period) +

2. s(Weekly, bs = "ps", k = 7),

3. data = matrix_gam,

4. family = gaussian)

首先是可視化。

1. layout(matrix(1:2, nrow = 1))

2. plot(gam_1, shade = TRUE)

我們在這里可以看到變量對電力負荷的影響。在左圖中，白天的負載峰值約為下午3點。在右邊的圖中，我們可以看到在周末負載量減少了。

讓我們使用summary函數(shù)對第一個模型進行診斷。

1. ##

2. ## Family: gaussian

3. ## Link function: identity

4. ##

5. ## Formula:

6. ## Load ~ s(Daily, bs = "cr", k = period) + s(Weekly, bs = "ps",

7. ## ? ? k = 7)

8. ##

9. ## Parametric coefficients:

10. ## ? ? ? ? ? ? Estimate Std. Error t value Pr(>|t|)

11. ## (Intercept) ?2731.67 ? ? ?18.88 ? 144.7 ? <2e-16 ***

12. ## ---

13. ## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

14. ##

15. ## Approximate significance of smooth terms:

16. ## ? ? ? ? ? ? ?edf Ref.df ? ? F p-value

17. ## s(Daily) ?10.159 12.688 119.8 ?<2e-16 ***

18. ## s(Weekly) ?5.311 ?5.758 130.3 ?<2e-16 ***

19. ## ---

20. ## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

21. ##

22. ## R-sq.(adj) = ?0.772 ? Deviance explained = 77.7%

23. ## GCV = 2.4554e+05 ?Scale est. = 2.3953e+05 ?n = 672

EDF：估計的自由度–可以像對給定變量進行平滑處理那樣來解釋（較高的EDF值表示更復(fù)雜的樣條曲線）。P值：給定變量對因變量的統(tǒng)計顯著性，通過F檢驗進行檢驗（越低越好）。調(diào)整后的R平方（越高越好）。我們可以看到R-sq.（adj）值有點低。

讓我們繪制擬合值：

我們需要將兩個自變量的交互作用包括到模型中。

第一種交互類型對兩個變量都使用了一個平滑函數(shù)。

1. gam_2 <- gam(Load ~ s(Daily, Weekly),

2. 


3. 


4. summary(gam_2)$r.sq

## [1] 0.9352108

R方值表明結(jié)果要好得多。

summary(gam_2)$s.table

1. ## ? ? ? ? ? ? ? ? ? ? edf ? Ref.df ? ? ? ?F p-value

2. ## s(Daily,Weekly) 28.7008 28.99423 334.4754 ? ? ? 0

似乎也很好，p值為0，這意味著自變量很重要。擬合值圖：

現(xiàn)在，讓我們嘗試上述張量積交互。這可以通過function完成te，也可以定義基本函數(shù)。

## [1] 0.9268452

與以前的模型相似gam_2。

summary(gam_3)$s.table

1. ## ? ? ? ? ? ? ? ? ? ? ? edf ? Ref.df ? ? ? ?F p-value

2. ## te(Daily,Weekly) 23.65709 23.98741 354.5856 ? ? ? 0

非常相似的結(jié)果。讓我們看一下擬合值：

與gam_2模型相比，只有一點點差異，看起來te擬合更好。

## [1] 0.9727604summary(gam_4)$sp.criterion

1. ## ? GCV.Cp

2. ## 34839.46

summary(gam_4)$s.table

1. ## ? ? ? ? ? ? ? ? ? ? ? edf ? Ref.df ? ? ? ?F p-value

2. ## te(Daily,Weekly) 119.4117 149.6528 160.2065 ? ? ? 0

我們可以在這里看到R方略有上升。
讓我們繪制擬合值：

這似乎比gam_3模型好得多。

## [1] 0.965618summary(gam_4_fx)$s.table

1. ## ? ? ? ? ? ? ? ? ?edf Ref.df ? ? ? ?F ? ? ? p-value

2. ## te(Daily,Weekly) 335 ? ?335 57.25389 5.289648e-199

我們可以看到R平方比模型gam_4低，這是因為我們過度擬合了模型。證明GCV程序（lambda和EDF的估計）工作正常。

因此，讓我們在案例（模型）中嘗試ti方法。

## [1] 0.9717469summary(gam_5)$sp.criterion

1. ## ? GCV.Cp

2. ## 35772.35

summary(gam_5)$s.table

1. ## ? ? ? ? ? ? ? ? ? ? ? ?edf ? ? Ref.df ? ? ? ? ?F p-value

2. ## s(Daily) ? ? ? ? 22.583649 ?27.964970 ?444.19962 ? ? ? 0

3. ## s(Weekly) ? ? ? ? 5.914531 ? 5.995934 1014.72482 ? ? ? 0

4. ## ti(Daily,Weekly) 85.310314 110.828814 ? 41.22288 ? ? ? 0

然后使用t2。

## [1] 0.9738273summary(gam_6)$sp.criterion

1. ## ? GCV.Cp

2. ## 32230.68

summary(gam_6)$s.table

1. ## ? ? ? ? ? ? ? ? ? ? ? edf ? Ref.df ? ? ? ?F p-value

2. ## t2(Daily,Weekly) 98.12005 120.2345 86.70754 ? ? ? 0

我還輸出了最后三個模型的GCV得分值，這也是在一組擬合模型中選擇最佳模型的良好標準。我們可以看到，對于t2相應(yīng)模型gam_6，GCV值最低。

在統(tǒng)計中廣泛使用的其他模型選擇標準是AIC（Akaike信息準則）。讓我們看看三個模型：

AIC(gam_4, gam_5, gam_6)

1. ## ? ? ? ? ? ? df ? ? ?AIC

2. ## gam_4 121.4117 8912.611

3. ## gam_5 115.8085 8932.746

4. ## gam_6 100.1200 8868.628

最低值在gam_6模型中。讓我們再次查看擬合值。

我們可以看到的模型的擬合值gam_4和gam_6非常相似?？梢允褂密浖母嗫梢暬湍Ｐ驮\斷功能來比較這兩個模型。

第一個是function?gam.check，它繪制了四個圖：殘差的QQ圖，線性預(yù)測變量與殘差，殘差的直方圖以及擬合值與因變量的關(guān)系圖。讓我們診斷模型gam_4和gam_6。

gam.check(gam_4)

1. ##

2. ## Method: GCV ? Optimizer: magic

3. ## Smoothing parameter selection converged after 7 iterations.

4. ## The RMS GCV score gradiant at convergence was 0.2833304 .

5. ## The Hessian was positive definite.

6. ## The estimated model rank was 336 (maximum possible: 336)

7. ## Model rank = ?336 / 336

8. ##

9. ## Basis dimension (k) checking results. Low p-value (k-index<1) may

10. ## indicate that k is too low, especially if edf is close to k'.

11. ##

12. ## ? ? ? ? ? ? ? ? ? ? ?k' ? ?edf k-index p-value

13. ## te(Daily,Weekly) 335.00 119.41 ? ?1.22 ? ? ? 1

gam.check(gam_6)

1. ##

2. ## Method: GCV ? Optimizer: magic

3. ## Smoothing parameter selection converged after 9 iterations.

4. ## The RMS GCV score gradiant at convergence was 0.05208856 .

5. ## The Hessian was positive definite.

6. ## The estimated model rank was 336 (maximum possible: 336)

7. ## Model rank = ?336 / 336

8. ##

9. ## Basis dimension (k) checking results. Low p-value (k-index<1) may

10. ## indicate that k is too low, especially if edf is close to k'.

11. ##

12. ## ? ? ? ? ? ? ? ? ? ? ?k' ? ?edf k-index p-value

13. ## t2(Daily,Weekly) 335.00 ?98.12 ? ?1.18 ? ? ? 1

我們可以再次看到模型非常相似，只是在直方圖中可以看到一些差異。

?

1. layout(matrix(1:2, nrow = 1))

2. plot(gam_4, rug = FALSE, se = FALSE, n2 = 80, main = "gam n.4 with te()")

3. plot(gam_6, rug = FALSE, se = FALSE, n2 = 80, main = "gam n.6 with t2()")

?該模型gam_6?有更多的“波浪形”的輪廓。因此，這意味著它對因變量的擬合度更高，而平滑因子更低。?

?

1. vis.gam(gam_6, n.grid = 50, theta = 35, phi = 32, zlab = "",

2. ticktype = "detailed", color = "topo", main = "t2(D, W)")

我們可以看到最高峰值是Daily變量的值接近30（下午3點），而Weekly變量的值是1（星期一）。

?

1. vis.gam(gam_6, main = "t2(D, W)", plot.type = "contour",

2. color = "terrain", contour.col = "black", lwd = 2)

再次可以看到，電力負荷的最高值是星期一的下午3:00，直到星期四都非常相似，然后負荷在周末減少。

?

最受歡迎的見解

1.在python中使用lstm和pytorch進行時間序列預(yù)測

2.python中利用長短期記憶模型lstm進行時間序列預(yù)測分析

3.使用r語言進行時間序列（arima，指數(shù)平滑）分析

4.r語言多元copula-garch-模型時間序列預(yù)測

5.r語言copulas和金融時間序列案例

6.使用r語言隨機波動模型sv處理時間序列中的隨機波動

7.r語言時間序列tar閾值自回歸模型

8.r語言k-shape時間序列聚類方法對股票價格時間序列聚類

9.python3用arima模型進行時間序列預(yù)測

?