由于數(shù)據(jù)科學機器學習和深度學習的發(fā)展，時間序列預測在預測準確性方面取得了顯著進展。隨著這些 ML/DL 工具的發(fā)展，企業(yè)和金融機構現(xiàn)在可以通過應用這些新技術來解決舊問題，從而更好地進行預測。在本文中，我們展示了使用稱為 LSTM（長短期記憶）的特殊類型深度學習模型，該模型對涉及自相關性的序列預測問題很有用。我們分析了一個名為“太陽黑子”的著名歷史數(shù)據(jù)集（太陽黑子是指太陽表面形成黑點的太陽現(xiàn)象)。我們將展示如何使用 LSTM 模型預測未來 10 年的太陽黑子數(shù)量。

教程概覽

此代碼教程對應于 2018 年 4 月 19 日星期四向 SP Global 提供的 Time Series Deep Learning 演示文稿?？梢韵螺d補充本文的幻燈片。

這是一個涉及時間序列深度學習和其他復雜機器學習主題（如回測交叉驗證）的高級教程。如果想要了解 R 中的 Keras，請查看：Customer Analytics: Using Deep Learning With Keras To Predict Customer Churn。

本教程中，你將會學到：

• 用 keras 包開發(fā)一個狀態(tài) LSTM 模型，該 R 包將 R TensorFlow 作為后端。

• 將狀態(tài) LSTM 模型應用到著名的太陽黑子數(shù)據(jù)集上。

• 借助 rsample 包在初始抽樣上滾動預測，實現(xiàn)時間序列的交叉檢驗。

• 借助 ggplot2 和 cowplot 可視化回測和預測結果。

• 通過自相關函數(shù)（Autocorrelation Function，ACF）圖評估時間序列數(shù)據(jù)是否適合應用 LSTM 模型。

本文的最終結果是一個高性能的深度學習算法，在預測未來 10 年太陽黑子數(shù)量方面表現(xiàn)非常出色！這是回測后狀態(tài) LSTM 模型的結果。

商業(yè)應用

時間序列預測對營收和利潤有顯著影響。在商業(yè)方面，我們可能有興趣預測每月、每季度或每年的哪一天會發(fā)生大額支出，或者我們可能有興趣了解消費者物價指數(shù)（CPI）在未來六年個月如何變化。這些都是在微觀經(jīng)濟和宏觀經(jīng)濟層面影響商業(yè)組織的常見問題。雖然本教程中使用的數(shù)據(jù)集不是“商業(yè)”數(shù)據(jù)集，但它顯示了工具-問題匹配的強大力量，意味著使用正確的工具進行工作可以大大提高準確性。最終的結果是預測準確性的提高將對營收和利潤帶來可量化的提升。

長短期記憶（LSTM）模型

長短期記憶（LSTM）模型是一種強大的遞歸神經(jīng)網(wǎng)絡（RNN）。博文《Understanding LSTM Networks》（翻譯版)以簡單易懂的方式解釋了模型的復雜性機制。下面是描述 LSTM 內(nèi)部單元架構的示意圖，除短期狀態(tài)之外，該架構使其能夠保持長期狀態(tài)，而這是傳統(tǒng)的 RNN 處理起來有困難的：

來源：Understanding LSTM Networks

LSTM 模型在預測具有自相關性（時間序列和滯后項之間存在相關性）的時間序列時非常有用，因為模型能夠保持狀態(tài)并識別時間序列上的模式。在每次處理過程中，遞歸架構能使狀態(tài)在更新權重時保持或者傳遞下去。此外，LSTM 模型的單元架構在短期持久化的基礎上實現(xiàn)了長期持久化，進而強化了 RNN，這一點非常吸引人！

在 Keras 中，LSTM 模型可以有“狀態(tài)”模式，Keras 文檔中這樣解釋：

索引 i 處每個樣本的最后狀態(tài)將被用作下一次批處理中索引 i 處樣本的初始狀態(tài)

在正常（或“無狀態(tài)”）模式下，Keras 對樣本重新洗牌，時間序列與其滯后項之間的依賴關系丟失。但是，在“狀態(tài)”模式下運行時，我們通?？梢酝ㄟ^利用時間序列中存在的自相關性來獲得高質量的預測結果。

在完成本教程時，我們會進一步解釋。就目前而言，可以認為 LSTM 模型對涉及自相關性的時間序列問題可能非常有用，而且 Keras 有能力創(chuàng)建完美的時間序列建模工具——狀態(tài) LSTM 模型。

太陽黑子數(shù)據(jù)集

太陽黑子是隨 R 發(fā)布的著名數(shù)據(jù)集（參見 datasets 包）。數(shù)據(jù)集跟蹤記錄太陽黑子，即太陽表面出現(xiàn)黑點的事件。這是來自 NASA 的一張照片，顯示了太陽黑子現(xiàn)象。相當酷！

來源：NASA

本教程所用的數(shù)據(jù)集稱為 sunspots.month，包含了 265（1749 ～ 2013）年間每月太陽黑子數(shù)量的月度數(shù)據(jù)。

構建 LSTM 模型預測太陽黑子

讓我們開動起來，預測太陽黑子。這是我們的目標：

目標：使用 LSTM 模型預測未來 10 年的太陽黑子數(shù)量。

1 若干相關包

以下是本教程所需的包，所有這些包都可以在 CRAN 上找到。如果你尚未安裝這些包，可以使用 install.packages() 進行安裝。注意：在繼續(xù)使用此代碼教程之前，請確保更新所有包，因為這些包的先前版本可能與所用代碼不兼容。

# Core Tidyverse
library(tidyverse)
library(glue)
library(forcats)

# Time Series
library(timetk)
library(tidyquant)
library(tibbletime)

# Visualization
library(cowplot)

# Preprocessing
library(recipes)

# Sampling / Accuracy
library(rsample)
library(yardstick)

# Modeling
library(keras)

如果你之前沒有在 R 中運行過 Keras，你需要用 install_keras() 函數(shù)安裝 Keras。

# Install Keras if you have not installed before
install_keras()

2 數(shù)據(jù)

數(shù)據(jù)集 sunspot.month 隨 R 一起發(fā)布，可以輕易獲得。它是一個 ts 類對象（非 tidy 類），所以我們將使用 timetk 中的 tk_tbl() 函數(shù)轉換為 tidy 數(shù)據(jù)集。我們使用這個函數(shù)而不是來自 tibble 的 as.tibble()，用來自動將時間序列索引保存為zoo yearmon 索引。最后，我們將使用 lubridate::as_date()（使用 tidyquant 時加載）將 zoo 索引轉換為日期，然后轉換為 tbl_time 對象以使時間序列操作起來更容易。

sun_spots <- datasets::sunspot.month %>%
   tk_tbl() %>%    
   mutate(index = as_date(index)) %>%    
   as_tbl_time(index = index) sun_spots

## # A time tibble: 3,177 x 2##
# Index: index
##    index      value
##    <date>     <dbl>
##  1 1749-01-01  58.0
##  2 1749-02-01  62.6
##  3 1749-03-01  70.0
##  4 1749-04-01  55.7
##  5 1749-05-01  85.0
##  6 1749-06-01  83.5
##  7 1749-07-01  94.8
##  8 1749-08-01  66.3
##  9 1749-09-01  75.9
## 10 1749-10-01  75.5
## # ... with 3,167 more rows

3 探索性數(shù)據(jù)分析

時間序列很長（有 265 年！）。我們可以將時間序列的全部（265 年）以及前 50 年的數(shù)據(jù)可視化，以獲得該時間系列的直觀感受。

3.1 使用 COWPLOT 可視化太陽黑子數(shù)據(jù)

我們將創(chuàng)建若干 ggplot 對象并借助 cowplot::plot_grid() 把這些對象組合起來。對于需要縮放的部分，我們使用 tibbletime::time_filter()，可以方便的實現(xiàn)基于時間的過濾。

p1 <- sun_spots %>%
   ggplot(aes(index, value))    
   geom_point(
           color = palette_light()[[1]], alpha = 0.5)    
    theme_tq()    
    labs(title = 'From 1749 to 2013 (Full Data Set)') p2 <- sun_spots %>%
        filter_time('start' ~'1800') %>%    
        ggplot(aes(index, value))    
        geom_line(color = palette_light()[[1]], alpha = 0.5)    
        geom_point(color = palette_light()[[1]])    
        geom_smooth(method = 'loess', span = 0.2, se = FALSE)    
        theme_tq()    
        labs(
                title = '1749 to 1800 (Zoomed In To Show Cycle)',        
                caption = 'datasets::sunspot.month') p_title <- ggdraw()
                     draw_label(
                             'Sunspots',        
                             size = 18,        
                             fontface = 'bold',        
                             colour = palette_light()[[1]])
                             
plot_grid(    p_title, p1, p2,
   ncol = 1,
   rel_heights = c(0.1, 1, 1))

乍一看，這個時間序列應該很容易預測。但是，我們可以看到，周期（10 年）和振幅（太陽黑子的數(shù)量）似乎在 1780 年至 1800 年之間發(fā)生變化。這產(chǎn)生了一些挑戰(zhàn)。

接下來我們要做的是確定 LSTM 模型是否是一個適用的好方法。LSTM 模型利用自相關性產(chǎn)生序列預測。我們的目標是使用批量預測（一種在整個預測區(qū)域內(nèi)創(chuàng)建單一預測批次的技術，不同于在未來一個或多個步驟中迭代執(zhí)行的單一預測）產(chǎn)生未來 10 年的預測。批量預測只有在自相關性持續(xù) 10 年以上時才有效。下面，我們來檢查一下。

首先，我們需要回顧自相關函數(shù)（Autocorrelation Function，ACF），它表示時間序列與自身滯后項之間的相關性。stats 包庫中的 acf() 函數(shù)以曲線的形式返回每個滯后階數(shù)的 ACF 值。但是，我們希望將 ACF 值提取出來以便研究。為此，我們將創(chuàng)建一個自定義函數(shù) tidy_acf()，以 tidy tibble 的形式返回 ACF 值。

tidy_acf <- function(data,

                              value,

                              lags = 0:20) {

        value_expr <- enquo(value)

        acf_values <- data %>%

               pull(value) %>%

               acf(lag.max = tail(lags, 1), plot = FALSE) %>%

               .$acf %>%

               .[,,1]        

        ret <- tibble(acf = acf_values) %>%

               rowid_to_column(var = 'lag') %>%

               mutate(lag = lag - 1) %>%

               filter(lag %in% lags)        

        return(ret) 

}

接下來，讓我們測試一下這個函數(shù)以確保它按預期工作。該函數(shù)使用我們的 tidy 時間序列，提取數(shù)值列，并以 tibble 的形式返回 ACF 值以及對應的滯后階數(shù)。我們有 601 個自相關系數(shù)（一個對應時間序列自身，剩下的對應 600 個滯后階數(shù)）。一切看起來不錯。

max_lag <- 12 * 50

sun_spots %>%
   tidy_acf(value, lags = 0:max_lag)

## # A tibble: 601 x 2
##      lag   acf
##    <dbl> <dbl>
##  1    0. 1.00
##  2    1. 0.923
##  3    2. 0.893
##  4    3. 0.878
##  5    4. 0.867
##  6    5. 0.853
##  7    6. 0.840
##  8    7. 0.822
##  9    8. 0.809
## 10    9. 0.799
## # ... with 591 more rows

下面借助 ggplot2 包把 ACF 數(shù)據(jù)可視化，以便確定 10 年后是否存在高自相關滯后項。

sun_spots %>%
   tidy_acf(value, lags = 0:max_lag) %>%
   ggplot(aes(lag, acf))
   geom_segment(
       aes(xend = lag, yend = 0),        
       color = palette_light()[[1]])
   geom_vline(
           xintercept = 120, size = 3,        
           color = palette_light()[[2]])    
   annotate(
           'text', label = '10 Year Mark',        
           x = 130, y = 0.8,        
           color = palette_light()[[2]],        
           size = 6, hjust = 0)
   theme_tq()
       labs(title = 'ACF: Sunspots')

好消息。自相關系數(shù)在 120 階（10年標志）之后依然超過 0.5。理論上，我們可以使用高自相關滯后項來開發(fā) LSTM 模型。

sun_spots %>%
   tidy_acf(value, lags = 115:135) %>%    
   ggplot(aes(lag, acf))    
   geom_vline(
           xintercept = 120, size = 3,        
           color = palette_light()[[2]])    
   geom_segment(
           aes(xend = lag, yend = 0),        
           color = palette_light()[[1]])
   geom_point(
           color = palette_light()[[1]],        
           size = 2)    
   geom_label(
           aes(label = acf %>% round(2)),
                   vjust = -1,        
                   color = palette_light()[[1]])
    annotate(
            'text', label = '10 Year Mark',        
            x = 121, y = 0.8,        
            color = palette_light()[[2]],        
            size = 5, hjust = 0)    
     theme_tq()    
     labs(
             title = 'ACF: Sunspots',        
             subtitle = 'Zoomed in on Lags 115 to 135')

經(jīng)過檢查，最優(yōu)滯后階數(shù)位于在 125 階。這不一定是我們將使用的，因為我們要更多地考慮使用 Keras 實現(xiàn)的 LSTM 模型進行批量預測。有了這個觀點，以下是如何 filter() 獲得最優(yōu)滯后階數(shù)。

optimal_lag_setting <- sun_spots %>%
   tidy_acf(value, lags = 115:135) %>%    
   filter(acf == max(acf)) %>%    
   pull(lag) optimal_lag_setting

## [1] 125

4 回測：時間序列交叉驗證

交叉驗證是在子樣本數(shù)據(jù)上針對驗證集數(shù)據(jù)開發(fā)模型的過程，其目標是確定預期的準確度級別和誤差范圍。在交叉驗證方面，時間序列與非序列數(shù)據(jù)有點不同。具體而言，在制定抽樣計劃時，必須保留對以前時間樣本的時間依賴性。我們可以通過平移窗口的方式選擇連續(xù)子樣本，進而創(chuàng)建交叉驗證抽樣計劃。在金融領域，這種類型的分析通常被稱為“回測”，它需要在一個時間序列上平移若干窗口來分割成多個不間斷的序列，以在當前和過去的觀測上測試策略。

最近的一個發(fā)展是 rsample 包，它使交叉驗證抽樣計劃非常易于實施。此外，rsample 包還包含回測功能?！癟ime Series Analysis Example”描述了一個使用 rolling_origin() 函數(shù)為時間序列交叉驗證創(chuàng)建樣本的過程。我們將使用這種方法。

4.1 開發(fā)一個回測策略

我們創(chuàng)建的抽樣計劃使用 50 年（initial = 12 x 50）的數(shù)據(jù)作為訓練集，10 年（assess = 12 x 10）的數(shù)據(jù)用于測試（驗證）集。我們選擇 20 年的跳躍跨度（skip = 12 x 20），將樣本均勻分布到 11 組中，跨越整個 265 年的太陽黑子歷史。最后，我們選擇 cumulative = FALSE 來允許平移起始點，這確保了較近期數(shù)據(jù)上的模型相較那些不太新近的數(shù)據(jù)沒有不公平的優(yōu)勢（使用更多的觀測數(shù)據(jù)）。rolling_origin_resamples 是一個 tibble 型的返回值。

periods_train <- 12 * 50
periods_test  <- 12 * 10
skip_span     <- 12 * 20
rolling_origin_resamples <- rolling_origin(    sun_spots,
   initial    = periods_train,    
   assess     = periods_test,    
   cumulative = FALSE,    
   skip       = skip_span) rolling_origin_resamples

## # Rolling origin forecast resampling
## # A tibble: 11 x 2
##    splits       id    
##    <list>       <chr>  
##  1 <S3: rsplit> Slice01
##  2 <S3: rsplit> Slice02
##  3 <S3: rsplit> Slice03
##  4 <S3: rsplit> Slice04
##  5 <S3: rsplit> Slice05
##  6 <S3: rsplit> Slice06
##  7 <S3: rsplit> Slice07
##  8 <S3: rsplit> Slice08
##  9 <S3: rsplit> Slice09
## 10 <S3: rsplit> Slice10
## 11 <S3: rsplit> Slice11

4.2 可視化回測策略

我們可以用兩個自定義函數(shù)來可視化再抽樣。首先是 plot_split()，使用 ggplot2 繪制一個再抽樣分割圖。請注意，expand_y_axis 參數(shù)默認將日期范圍擴展成整個 sun_spots 數(shù)據(jù)集的日期范圍。當我們將所有的圖形同時可視化時，這將變得有用。

# Plotting function for a single split

plot_split <- function(split,
                      expand_y_axis = TRUE,                      
                      alpha = 1,                      
                      size = 1,                      
                      base_size = 14) {        # Manipulate data    train_tbl <- training(split) %>%        
                      add_column(key = 'training')        test_tbl  <- testing(split) %>%        
                      add_column(key = 'testing')        data_manipulated <- bind_rows(        train_tbl, test_tbl) %>%
       as_tbl_time(index = index) %>%        
       mutate(            
               key = fct_relevel(                    key, 'training', 'testing'))        
   
   # Collect attributes    train_time_summary <- train_tbl %>%
           tk_index() %>%        
           tk_get_timeseries_summary()        test_time_summary <- test_tbl %>%        
           tk_index() %>%        
           tk_get_timeseries_summary()
                   # Visualize    g <- data_manipulated %>%            ggplot(                aes(x = index,

                       y = value,                

                       color = key))        

            geom_line(size = size, alpha = alpha)

            theme_tq(base_size = base_size)

            scale_color_tq()

            labs(

                title    = glue('Split: {split$id}'),

                subtitle = glue(

                        '{train_time_summary$start} to {test_time_summary$end}'),

                                 y = '', x = '')

          theme(legend.position = 'none')        if (expand_y_axis) {                sun_spots_time_summary <- sun_spots %>%

            tk_index() %>%

            tk_get_timeseries_summary()                g <- g

           scale_x_date(

               limits = c(                    sun_spots_time_summary$start,                    sun_spots_time_summary$end))    }        return(g) }

plot_split() 函數(shù)接受一個分割（在本例中為 Slice01），并可視化抽樣策略。我們使用 expand_y_axis = TRUE 將橫坐標范圍擴展到整個數(shù)據(jù)集的日期范圍。

rolling_origin_resamples$splits[[1]] %>%
   plot_split(expand_y_axis = TRUE)
   theme(legend.position = 'bottom')

第二個函數(shù)是 plot_sampling_plan()，使用 purrr 和 cowplot 將 plot_split() 函數(shù)應用到所有樣本上。

# Plotting function that scales to all splits

plot_sampling_plan <- function(sampling_tbl,                               expand_y_axis = TRUE,                               ncol = 3,

                              alpha = 1,

                              size = 1,

                              base_size = 14,                               title = 'Sampling Plan') {        # Map plot_split() to sampling_tbl    sampling_tbl_with_plots <- sampling_tbl %>%

       mutate(

           gg_plots = map(                splits, plot_split,                expand_y_axis = expand_y_axis,

               alpha = alpha,

               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) }

現(xiàn)在我們可以使用 plot_sampling_plan() 可視化整個回測策略！我們可以看到抽樣計劃如何平移抽樣窗口逐漸切分出訓練和測試子樣本。

rolling_origin_resamples %>%
   plot_sampling_plan(
           expand_y_axis = T,        
           ncol = 3, alpha = 1,        
           size = 1, base_size = 10,            title = 'Backtesting Strategy: Rolling Origin Sampling Plan')

此外，我們可以讓 expand_y_axis = FALSE，對每個樣本進行縮放。

rolling_origin_resamples %>%
   plot_sampling_plan(
           expand_y_axis = F,        
           ncol = 3, alpha = 1,       
           size = 1, base_size = 10,         
           title = 'Backtesting Strategy: Zoomed In')

當在太陽黑子數(shù)據(jù)集上測試 LSTM 模型準確性時，我們將使用這種回測策略（來自一個時間序列的 11 個樣本，每個時間序列分為 50/10 兩部分，并且樣本之間有 20 年的偏移）。