博客專欄：
https://www.cnblogs.com/xuruilong100
5 用 Keras 構(gòu)建狀態(tài) LSTM 模型
首先，我們將在回測策略的某個樣本上用 Keras 開發(fā)一個狀態(tài) LSTM 模型。然后，我們將模型套用到所有樣本，以測試和驗證模型性能。
5.1 單個 LSTM 模型
對單個 LSTM 模型，我們選擇并可視化最近一期的分割樣本（Slice11），這一樣本包含了最新的數(shù)據(jù)。
split <- rolling_origin_resamples$splits[[11]]
split_id <- rolling_origin_resamples$id[[11]]
5.1.1 可視化該分割樣本
我么可以用 plot_split() 函數(shù)可視化該分割，設(shè)定 expand_y_axis = FALSE 以便將橫坐標縮放到樣本本身的范圍。
plot_split(
split,
expand_y_axis = FALSE,
size = 0.5)
theme(legend.position = 'bottom')
ggtitle(glue('Split: {split_id}'))
5.1.2 數(shù)據(jù)準備
首先，我們將訓(xùn)練和測試數(shù)據(jù)集合成一個數(shù)據(jù)集，并使用列 key 來標記它們來自哪個集合（training 或 testing）。請注意，tbl_time 對象需要在調(diào)用 bind_rows() 時重新指定索引，但是這個問題應(yīng)該很快在 dplyr 包中得到糾正。
df_trn <- training(split)
df_tst <- testing(split)
df <- bind_rows(
df_trn %>% add_column(key = 'training'),
df_tst %>% add_column(key = 'testing')) %>%
as_tbl_time(index = index)
df
## # A time tibble: 720 x 3
## # Index: index
## index value key
## <date> <dbl> <chr>
## 1 1949-11-01 144. training
## 2 1949-12-01 118. training
## 3 1950-01-01 102. training
## 4 1950-02-01 94.8 training
## 5 1950-03-01 110. training
## 6 1950-04-01 113. training
## 7 1950-05-01 106. training
## 8 1950-06-01 83.6 training
## 9 1950-07-01 91.0 training
## 10 1950-08-01 85.2 training
## # ... with 710 more rows
5.1.3 用 recipe 做數(shù)據(jù)預(yù)處理
LSTM 算法要求輸入數(shù)據(jù)經(jīng)過中心化并標度化。我們可以使用 recipe 包預(yù)處理數(shù)據(jù)。我們用 step_sqrt 來轉(zhuǎn)換數(shù)據(jù)以減少異常值的影響，再結(jié)合 step_center 和 step_scale 對數(shù)據(jù)進行中心化和標度化。最后，數(shù)據(jù)使用 bake() 函數(shù)實現(xiàn)處理轉(zhuǎn)換。
rec_obj <- recipe(value ~ ., df) %>%
step_sqrt(value) %>%
step_center(value) %>%
step_scale(value) %>%
prep()
df_processed_tbl <- bake(rec_obj, df)
df_processed_tbl
## # A tibble: 720 x 3
## index value key
## <date> <dbl> <fct>
## 1 1949-11-01 1.25 training
## 2 1949-12-01 0.929 training
## 3 1950-01-01 0.714 training
## 4 1950-02-01 0.617 training
## 5 1950-03-01 0.825 training
## 6 1950-04-01 0.874 training
## 7 1950-05-01 0.777 training
## 8 1950-06-01 0.450 training
## 9 1950-07-01 0.561 training
## 10 1950-08-01 0.474 training
## # ... with 710 more rows
接著，記錄中心化和標度化的信息，以便在建模完成之后可以將數(shù)據(jù)逆向轉(zhuǎn)換回去。平方根轉(zhuǎn)換可以通過乘方運算逆轉(zhuǎn)回去，但要在逆轉(zhuǎn)中心化和標度化之后。
center_history <- rec_obj$steps[[2]]$means['value']
scale_history <- rec_obj$steps[[3]]$sds['value']
c('center' = center_history, 'scale' = scale_history)
## center.value scale.value
## 7.549526 3.545561
5.1.4 規(guī)劃 LSTM 模型
我們需要規(guī)劃下如何構(gòu)建 LSTM 模型。首先，了解幾個 LSTM 模型的專業(yè)術(shù)語：
張量格式（Tensor Format）：
預(yù)測變量（X）必須是一個 3 維數(shù)組，維度分別是：samples、timesteps 和 features。第一維代表變量的長度；第二維是時間步（滯后階數(shù)）；第三維是預(yù)測變量的個數(shù)（1 表示單變量，n 表示多變量）
輸出或目標變量（y）必須是一個 2 維數(shù)組，維度分別是：samples 和 timesteps。第一維代表變量的長度；第二維是時間步（之后階數(shù)）
訓(xùn)練與測試：
批量大?。˙atch Size）：
時間步（Time Steps）:
周期（Epochs）：
考慮到這一點，我們可以提出一個計劃。我們將預(yù)測窗口或測試集的長度定在 120 個月（10年）。最優(yōu)相關(guān)性發(fā)生在 125 階，但這并不能被預(yù)測范圍整除。我們可以增加預(yù)測范圍，但是這僅提供了自相關(guān)性的最小幅度增加。我們選擇批量大小為 40，它可以整除測試集和訓(xùn)練集的觀察個數(shù)。我們選擇時間步等于 1，這是因為我們只使用 1 階滯后（只向前預(yù)測一步）。最后，我們設(shè)置 epochs = 300，但這需要調(diào)整以平衡偏差與方差。
# Model inputs
lag_setting <- 120 # = nrow(df_tst)
batch_size <- 40
train_length <- 440
tsteps <- 1
epochs <- 300
5.1.5 2 維與 3 維的訓(xùn)練、測試數(shù)組
下面將訓(xùn)練集和測試集數(shù)據(jù)轉(zhuǎn)換成合適的形式（數(shù)組）。記住，LSTM 模型要求預(yù)測變量（X）是 3 維的，輸出或目標變量（y）是 2 維的。
# Training Set
lag_train_tbl <- df_processed_tbl %>%
mutate(value_lag = lag(value, n = lag_setting)) %>%
filter(!is.na(value_lag)) %>%
filter(key =='training') %>%
tail(train_length)
x_train_vec <- lag_train_tbl$value_lag
x_train_arr <- array(
data = x_train_vec, dim = c(length(x_train_vec), 1, 1))
y_train_vec <- lag_train_tbl$value
y_train_arr <- array(
data = y_train_vec, dim = c(length(y_train_vec), 1))
# Testing Set
lag_test_tbl <- df_processed_tbl %>%
mutate(
value_lag = lag(
value, n = lag_setting)) %>%
filter(!is.na(value_lag)) %>%
filter(key =='testing')
x_test_vec <- lag_test_tbl$value_lag
x_test_arr <- array(
data = x_test_vec,
dim = c(length(x_test_vec), 1, 1))
y_test_vec <- lag_test_tbl$value
y_test_arr <- array(
data = y_test_vec,
dim = c(length(y_test_vec), 1))
5.1.6 構(gòu)建 LSTM 模型
我們可以使用 keras_model_sequential() 構(gòu)建 LSTM 模型，并像堆磚塊一樣堆疊神經(jīng)網(wǎng)絡(luò)層。我們將使用兩個 LSTM 層，每層都設(shè)定 units = 50。第一個 LSTM 層接收所需的輸入形狀，即[時間步，特征數(shù)量]。批量大小就是我們的批量大小。我們將第一層設(shè)置為 return_sequences = TRUE 和 stateful = TRUE。第二層和前面相同，除了 batch_size（batch_size只需要在第一層中指定），另外 return_sequences = FALSE 不返回時間戳維度（從第一個 LSTM 層返回 2 維數(shù)組，而不是 3 維）。我們使用 layer_dense(units = 1)，這是 Keras 序列模型的標準結(jié)尾。最后，我們在 compile() 中使用 loss ='mae' 以及流行的 optimizer = 'adam'。
model <- keras_model_sequential()
model %>%
layer_lstm(
units = 50,
input_shape = c(tsteps, 1),
batch_size = batch_size,
return_sequences = TRUE,
stateful = TRUE) %>%
layer_lstm(
units = 50,
return_sequences = FALSE,
stateful = TRUE) %>%
layer_dense(units = 1)
model %>%
compile(loss = 'mae', optimizer = 'adam')
model
## Model
## ______________________________________________________________________
## Layer (type) Output Shape Param #
## ======================================================================
## lstm_1 (LSTM) (40, 1, 50) 10400
## ______________________________________________________________________
## lstm_2 (LSTM) (40, 50) 20200
## ______________________________________________________________________
## dense_1 (Dense) (40, 1) 51
## ======================================================================
## Total params: 30,651
## Trainable params: 30,651
## Non-trainable params: 0
## ______________________________________________________________________
5.1.7 擬合 LSTM 模型
下一步，我們使用一個 for 循環(huán)擬合狀態(tài) LSTM 模型（需要手動重置狀態(tài)）。有 300 個周期要循環(huán)，運行需要一點時間。我們設(shè)置 shuffle = FALSE 來保存序列，并且我們使用 reset_states() 在每個循環(huán)后手動重置狀態(tài)。
for (i in 1:epochs) {
model %>%
fit(x = x_train_arr,
y = y_train_arr,
batch_size = batch_size,
epochs = 1,
verbose = 1,
shuffle = FALSE)
model %>% reset_states()
cat('Epoch: ', i)
}
5.1.8 使用 LSTM 模型預(yù)測
然后，我們可以使用 predict() 函數(shù)對測試集 x_test_arr 進行預(yù)測。我們可以使用之前保存的 scale_history 和 center_history 轉(zhuǎn)換得到的預(yù)測，然后對結(jié)果進行平方。最后，我們使用 reduce() 和自定義的 time_bind_rows() 函數(shù)將預(yù)測與一列原始數(shù)據(jù)結(jié)合起來。
# Make Predictions
pred_out <- model %>%
predict(x_test_arr, batch_size = batch_size) %>% .[,1]
# Retransform values
pred_tbl <- tibble(
index = lag_test_tbl$index,
value = (pred_out * scale_history center_history)^2)
# Combine actual data with predictions
tbl_1 <- df_trn %>%
add_column(key = 'actual')
tbl_2 <- df_tst %>%
add_column(key = 'actual')
tbl_3 <- pred_tbl %>%
add_column(key = 'predict')
# Create time_bind_rows() to solve dplyr issue
time_bind_rows <- function(data_1,
data_2, index) {
index_expr <- enquo(index)
bind_rows(data_1, data_2) %>%
as_tbl_time(index = !! index_expr)
}
ret <- list(tbl_1, tbl_2, tbl_3) %>%
reduce(time_bind_rows, index = index) %>%
arrange(key, index) %>%
mutate(key = as_factor(key))
ret
## # A time tibble: 840 x 3
## # Index: index
## index value key
## <date> <dbl> <fct>
## 1 1949-11-01 144. actual
## 2 1949-12-01 118. actual
## 3 1950-01-01 102. actual
## 4 1950-02-01 94.8 actual
## 5 1950-03-01 110. actual
## 6 1950-04-01 113. actual
## 7 1950-05-01 106. actual
## 8 1950-06-01 83.6 actual
## 9 1950-07-01 91.0 actual
## 10 1950-08-01 85.2 actual
## # ... with 830 more rows
5.1.9 評估單個分割樣本上 LSTM 模型的表現(xiàn)
我們使用 yardstick 包里的 rmse() 函數(shù)評估表現(xiàn)，rmse() 返回均方誤差平方根（RMSE）。我們的數(shù)據(jù)以“長”格式的形式存在（使用 ggplot2 可視化的最佳格式），所以需要創(chuàng)建一個包裝器函數(shù) calc_rmse() 對數(shù)據(jù)做預(yù)處理，以適應(yīng) yardstick::rmse() 的要求。
calc_rmse <- function(prediction_tbl) {
rmse_calculation <- function(data) {
data %>%
spread(key = key, value = value) %>%
select(-index) %>%
filter(!is.na(predict)) %>%
rename(
truth = actual,
estimate = predict) %>%
rmse(truth, estimate)
}
safe_rmse <- possibly(
rmse_calculation, otherwise = NA)
safe_rmse(prediction_tbl)
}
我們計算模型的 RMSE。
calc_rmse(ret)
## [1] 31.81798
RMSE 提供的信息有限，我們需要可視化。注意：當我們擴展到回測策略中的所有樣本時，RMSE 將在確定預(yù)期誤差時派上用場。
5.1.10 可視化一步預(yù)測
下一步，我們創(chuàng)建一個繪圖函數(shù)——plot_prediction()，借助 ggplot2 可視化單一樣本上的結(jié)果。
# Setup single plot function
plot_prediction <- function(data,
id,
alpha = 1,
size = 2,
base_size = 14) {
rmse_val <- calc_rmse(data)
g <- data %>%
ggplot(aes(index, value, color = key))
geom_point(alpha = alpha, size = size)
theme_tq(base_size = base_size)
scale_color_tq()
theme(legend.position = 'none')
labs(
title = glue(
'{id}, RMSE: {round(rmse_val, digits = 1)}'),
x = '', y = '')
return(g)
}
我們設(shè)置 id = split_id，在 Slice11 上測試函數(shù)。
ret %>%
plot_prediction(id = split_id, alpha = 0.65)
theme(legend.position = 'bottom')
LSTM 模型表現(xiàn)相對較好！ 我們選擇的設(shè)置似乎產(chǎn)生了一個不錯的模型，可以捕捉到數(shù)據(jù)中的趨勢。預(yù)測在下一個上升趨勢前搶跑了，但總體上好過了我的預(yù)期?，F(xiàn)在，我們需要通過回測來查看隨著時間推移的真實表現(xiàn)！
5.2 在 11 個樣本上回測 LSTM 模型
一旦我們有了能在一個樣本上工作的 LSTM 模型，擴展到全部 11 個樣本上就相對簡單。我們只需創(chuàng)建一個預(yù)測函數(shù)，再套用到 rolling_origin_resamples 中抽樣計劃包含的數(shù)據(jù)上。
5.2.1 構(gòu)建一個 LSTM 預(yù)測函數(shù)
這一步看起來很嚇人，但實際上很簡單。我們將 5.1 節(jié)的代碼復(fù)制到一個函數(shù)中。我們將它作為一個安全函數(shù)，對于任何長時間運行的函數(shù)來說，這是一個很好的做法，可以防止單個故障停止整個過程。
predict_keras_lstm <- function(split,
epochs = 300,
...) {
lstm_prediction <- function(split,
epochs,
...) {
# 5.1.2 Data Setup
df_trn <- training(split)
df_tst <- testing(split)
df <- bind_rows(
df_trn %>% add_column(key = 'training'),
df_tst %>% add_column(key = 'testing')) %>%
as_tbl_time(index = index)
# 5.1.3 Preprocessing
rec_obj <- recipe(value ~ ., df) %>%
step_sqrt(value) %>%
step_center(value) %>%
step_scale(value) %>%
prep()
df_processed_tbl <- bake(rec_obj, df)
center_history <- rec_obj$steps[[2]]$means['value']
scale_history <- rec_obj$steps[[3]]$sds['value']
# 5.1.4 LSTM Plan
lag_setting <- 120 # = nrow(df_tst)
batch_size <- 40
train_length <- 440
tsteps <- 1
epochs <- epochs
# 5.1.5 Train/Test Setup
lag_train_tbl <- df_processed_tbl %>%
mutate(
value_lag = lag(value, n = lag_setting)) %>%
filter(!is.na(value_lag)) %>%
filter(key =='training') %>%
tail(train_length)
x_train_vec <- lag_train_tbl$value_lag
x_train_arr <- array(
data = x_train_vec, dim = c(length(x_train_vec), 1, 1))
y_train_vec <- lag_train_tbl$value
y_train_arr <- array(
data = y_train_vec, dim = c(length(y_train_vec), 1))
lag_test_tbl <- df_processed_tbl %>%
mutate(
value_lag = lag(value, n = lag_setting)) %>%
filter(!is.na(value_lag)) %>%
filter(key =='testing')
x_test_vec <- lag_test_tbl$value_lag
x_test_arr <- array(
data = x_test_vec, dim = c(length(x_test_vec), 1, 1))
y_test_vec <- lag_test_tbl$value
y_test_arr <- array(
data = y_test_vec, dim = c(length(y_test_vec), 1))
# 5.1.6 LSTM Model
model <- keras_model_sequential()
model %>%
layer_lstm(
units = 50,
input_shape = c(tsteps, 1),
batch_size = batch_size,
return_sequences = TRUE,
stateful = TRUE) %>%
layer_lstm(
units = 50,
return_sequences = FALSE,
stateful = TRUE) %>%
layer_dense(units = 1)
model %>%
compile(loss = 'mae', optimizer = 'adam')
# 5.1.7 Fitting LSTM
for (i in 1:epochs) {
model %>%
fit(x = x_train_arr,
y = y_train_arr,
batch_size = batch_size,
epochs = 1,
verbose = 1,
shuffle = FALSE)
model %>% reset_states()
cat('Epoch: ', i)
}
# 5.1.8 Predict and Return Tidy Data
# Make Predictions
pred_out <- model %>%
predict(x_test_arr, batch_size = batch_size) %>%
.[,1]
# Retransform values
pred_tbl <- tibble(
index = lag_test_tbl$index,
value = (pred_out * scale_history center_history)^2)
# Combine actual data with predictions
tbl_1 <- df_trn %>%
add_column(key = 'actual')
tbl_2 <- df_tst %>%
add_column(key = 'actual')
tbl_3 <- pred_tbl %>%
add_column(key = 'predict')
# Create time_bind_rows() to solve dplyr issue
time_bind_rows <- function(data_1, data_2, index) {
index_expr <- enquo(index)
bind_rows(data_1, data_2) %>%
as_tbl_time(index = !! index_expr)
}
ret <- list(tbl_1, tbl_2, tbl_3) %>%
reduce(time_bind_rows, index = index) %>%
arrange(key, index) %>%
mutate(key = as_factor(key))
return(ret)
}
safe_lstm <- possibly(lstm_prediction, otherwise = NA)
safe_lstm(split, epochs, ...)
}
我們測試下 predict_keras_lstm() 函數(shù)，設(shè)置 epochs = 10。返回的數(shù)據(jù)為長格式，在 key 列中標記有 actual 和 predict。
predict_keras_lstm(split, epochs = 10)
## # A time tibble: 840 x 3
## # Index: index
## index value key
## <date> <dbl> <fct>
## 1 1949-11-01 144. actual
## 2 1949-12-01 118. actual
## 3 1950-01-01 102. actual
## 4 1950-02-01 94.8 actual
## 5 1950-03-01 110. actual
## 6 1950-04-01 113. actual
## 7 1950-05-01 106. actual
## 8 1950-06-01 83.6 actual
## 9 1950-07-01 91.0 actual
## 10 1950-08-01 85.2 actual
## # ... with 830 more rows
5.2.2 將 LSTM 預(yù)測函數(shù)應(yīng)用到 11 個樣本上
既然 predict_keras_lstm() 函數(shù)可以在一個樣本上運行，我們現(xiàn)在可以借助使用 mutate() 和 map() 將函數(shù)應(yīng)用到所有樣本上。預(yù)測將存儲在名為 predict 的列中。注意，這可能需要 5-10 分鐘左右才能完成。
sample_predictions_lstm_tbl <- rolling_origin_resamples %>%
mutate(predict = map(splits, predict_keras_lstm, epochs = 300))
現(xiàn)在，我們得到了 11 個樣本的預(yù)測，數(shù)據(jù)存儲在列 predict 中。
sample_predictions_lstm_tbl
## # Rolling origin forecast resampling
## # A tibble: 11 x 3
## splits id predict
## * <list> <chr> <list>
## 1 <S3: rsplit> Slice01 <tibble [840 x 3]>
## 2 <S3: rsplit> Slice02 <tibble [840 x 3]>
## 3 <S3: rsplit> Slice03 <tibble [840 x 3]>
## 4 <S3: rsplit> Slice04 <tibble [840 x 3]>
## 5 <S3: rsplit> Slice05 <tibble [840 x 3]>
## 6 <S3: rsplit> Slice06 <tibble [840 x 3]>
## 7 <S3: rsplit> Slice07 <tibble [840 x 3]>
## 8 <S3: rsplit> Slice08 <tibble [840 x 3]>
## 9 <S3: rsplit> Slice09 <tibble [840 x 3]>
## 10 <S3: rsplit> Slice10 <tibble [840 x 3]>
## 11 <S3: rsplit> Slice11 <tibble [840 x 3]>
5.2.3 評估回測表現(xiàn)
通過將 calc_rmse() 函數(shù)應(yīng)用到 predict 列上，我們可以得到所有樣本的 RMSE。
sample_rmse_tbl <- sample_predictions_lstm_tbl %>%
mutate(rmse = map_dbl(predict, calc_rmse)) %>%
select(id, rmse)
sample_rmse_tbl
## # Rolling origin forecast resampling
## # A tibble: 11 x 2
## id rmse
## * <chr> <dbl>
## 1 Slice01 48.2
## 2 Slice02 17.4
## 3 Slice03 41.0
## 4 Slice04 26.6
## 5 Slice05 22.2
## 6 Slice06 49.0
## 7 Slice07 18.1
## 8 Slice08 54.9
## 9 Slice09 28.0
## 10 Slice10 38.4
## 11 Slice11 34.2
sample_rmse_tbl %>%
ggplot(aes(rmse))
geom_histogram(
aes(y = ..density..),
fill = palette_light()[[1]], bins = 16)
geom_density(
fill = palette_light()[[1]], alpha = 0.5)
theme_tq()
ggtitle('Histogram of RMSE')
而且，我們可以總結(jié) 11 個樣本的 RMSE。專業(yè)提示：使用 RMSE（或其他類似指標）的平均值和標準差是比較各種模型表現(xiàn)的好方法。
sample_rmse_tbl %>%
summarize(
mean_rmse = mean(rmse),
sd_rmse = sd(rmse))
## # Rolling origin forecast resampling
## # A tibble: 1 x 2
## mean_rmse sd_rmse
## <dbl> <dbl>
## 1 34.4 13.0
5.2.4 可視化回測的結(jié)果
我們可以創(chuàng)建一個 plot_predictions() 函數(shù)，把 11 個回測樣本的預(yù)測結(jié)果繪制在一副圖上?。?！
plot_predictions <- function(sampling_tbl,
predictions_col,
ncol = 3,
alpha = 1,
size = 2,
base_size = 14,
title = 'Backtested Predictions') {
predictions_col_expr <- enquo(predictions_col)
# Map plot_split() to sampling_tbl
sampling_tbl_with_plots <- sampling_tbl %>%
mutate(
gg_plots = map2(
!! predictions_col_expr, id,
.f = plot_prediction,
alpha = alpha,
size = size,
base_size = base_size))
# Make plots with cowplot
plot_list <- sampling_tbl_with_plots$gg_plots
p_temp <- plot_list[[1]] theme(legend.position = 'bottom')
legend <- get_legend(p_temp)
p_body <- plot_grid(plotlist = plot_list, ncol = ncol)
p_title <- ggdraw()
draw_label(
title,
size = 18,
fontface = 'bold',
colour = palette_light()[[1]])
g <- plot_grid(
p_title,
p_body,
legend,
ncol = 1,
rel_heights = c(0.05, 1, 0.05))
return(g)
}
結(jié)果在這里。在一個不容易預(yù)測的數(shù)據(jù)集上，這是相當令人印象深刻的！
sample_predictions_lstm_tbl %>%
plot_predictions(
predictions_col = predict,
alpha = 0.5,
size = 1,
base_size = 10,
title = 'Keras Stateful LSTM: Backtested Predictions')
