目錄
翻譯自《Demo Week: Time Series Machine Learning with timetk》
原文鏈接：www.business-science.io/code-tools/2017/10/24/demo_week_timetk.html
時間序列分析工具箱——timetk
timetk 的主要用途
三個主要用途：
- 時間序列機器學習：使用回歸算法進行預測；
- 構造時間序列索引：基於時間模式提取、探索和擴展時間序列索引；
- 轉換不同類型的時間序列數據（例如
tbl、xts、zoo、ts之間）：輕松實現不同類型的時間序列數據之間的相互轉換。
我們今天將討論時間序列機器學習和數據類型轉換。第二個議題（提取和構造未來時間序列）將在時間序列機器學習中涉及，因為這對預測准確性非常關鍵。
加載包
需要加在兩個包：
tidyquant：用於獲取數據並在后台加載 tidyversetimetk：R 中用於處理時間序列的工具包
如果還沒有安裝過，請用下面的命令安裝：
# Install packages
install.packages("timetk")
install.packages("tidyquant")
加載包。
# Load libraries
library(timetk) # Toolkit for working with time series in R
library(tidyquant) # Loads tidyverse, financial pkgs, used to get data
數據
我們將使用 tidyquant 中的 tq_get() 函數從 FRED 獲取數據——啤酒、葡萄酒和蒸餾酒銷售數據。
# Beer, Wine, Distilled Alcoholic Beverages, in Millions USD
beer_sales_tbl <- tq_get(
"S4248SM144NCEN",
get = "economic.data",
from = "2010-01-01",
to = "2016-12-31")
beer_sales_tbl
## # A tibble: 84 x 2
## date price
## <date> <int>
## 1 2010-01-01 6558
## 2 2010-02-01 7481
## 3 2010-03-01 9475
## 4 2010-04-01 9424
## 5 2010-05-01 9351
## 6 2010-06-01 10552
## 7 2010-07-01 9077
## 8 2010-08-01 9273
## 9 2010-09-01 9420
## 10 2010-10-01 9413
## # ... with 74 more rows
可視化數據是一個好東西，這有助於幫助我們了解正在使用的是什么數據。可視化對於時間序列分析和預測尤為重要。我們將使用 tidyquant 畫圖工具：主要是用 geom_ma(ma_fun = SMA，n = 12) 來添加一個周期為 12 的簡單移動平均線來了解趨勢。我們還可以看到似乎同時存在着趨勢性（移動平均線以近似線性的模式增長）和季節性（波峰和波谷傾向於在特定月份發生）。
# Plot Beer Sales
beer_sales_tbl %>%
ggplot(aes(date, price)) +
geom_line(col = palette_light()[1]) +
geom_point(col = palette_light()[1]) +
geom_ma(ma_fun = SMA, n = 12, size = 1) +
theme_tq() +
scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
labs(title = "Beer Sales: 2007 through 2016")
現在你對我們要分析的時間序列有了直觀的感受，那么讓我們繼續！
timetk 教程：
教程分為兩部分。首先，我們將遵循時間序列機器學習的工作流程。其次，我們將介紹數據轉換工具。
PART 1：時間序列機器學習
時間序列機器學習是預測時間序列數據的一種很好的方法，但在我們開始之前，這里有幾個注意點：
- 關鍵洞察力：將時間序列簽名（時間戳信息按列擴展到特征集）用於執行機器學習。
- 目標：我們將使用時間序列簽名預測未來 12 個月的時間序列數據。
我們將遵循可用於執行時間序列機器學習的工作流程。你將看到幾個 timetk 函數如何幫助完成此過程。我們將使用簡單的 lm() 線性回歸進行機器學習，你將看到使用時間序列簽名會使機器學習更強大和准確。此外，你還應該考慮使用其他更強大的機器學習算法，例如 xgboost、glmnet（LASSO）等。
STEP 0：回顧數據
看看我們的起點，先打印出 beer_sales_tbl。
# Starting point
beer_sales_tbl
## # A tibble: 84 x 2
## date price
## <date> <int>
## 1 2010-01-01 6558
## 2 2010-02-01 7481
## 3 2010-03-01 9475
## 4 2010-04-01 9424
## 5 2010-05-01 9351
## 6 2010-06-01 10552
## 7 2010-07-01 9077
## 8 2010-08-01 9273
## 9 2010-09-01 9420
## 10 2010-10-01 9413
## # ... with 74 more rows
我們可以使用 tk_index() 來提取索引；使用 tk_get_timeseries_summary() 來檢索索引的摘要信息，進而快速了解時間序列。我們使用 glimpse() 輸出一個更漂亮的格式。
beer_sales_tbl %>%
tk_index() %>%
tk_get_timeseries_summary() %>%
glimpse()
## Observations: 1
## Variables: 12
## $ n.obs <int> 84
## $ start <date> 2010-01-01
## $ end <date> 2016-12-01
## $ units <chr> "days"
## $ scale <chr> "month"
## $ tzone <chr> "UTC"
## $ diff.minimum <dbl> 2419200
## $ diff.q1 <dbl> 2592000
## $ diff.median <dbl> 2678400
## $ diff.mean <dbl> 2629475
## $ diff.q3 <dbl> 2678400
## $ diff.maximum <dbl> 2678400
你可以看到一些重要的特征，例如起始、結束、單位等等。還有時間差的分位數（相鄰兩個觀察之間差距的秒數），這對評估規律性的程度很有用。由於時間尺度是月度的，因此每個月之間差距的秒數並不規則。
STEP 1：擴充時間序列簽名
tk_augment_timeseries_signature() 函數將時間戳信息逐列擴展到機器學習特征集中，並將時間序列信息列添加到初始數據表。
# Augment (adds data frame columns)
beer_sales_tbl_aug <- beer_sales_tbl %>%
tk_augment_timeseries_signature()
beer_sales_tbl_aug
## # A tibble: 84 x 30
## date price index.num diff year year.iso half quarter
## <date> <int> <int> <int> <int> <int> <int> <int>
## 1 2010-01-01 6558 1262304000 NA 2010 2009 1 1
## 2 2010-02-01 7481 1264982400 2678400 2010 2010 1 1
## 3 2010-03-01 9475 1267401600 2419200 2010 2010 1 1
## 4 2010-04-01 9424 1270080000 2678400 2010 2010 1 2
## 5 2010-05-01 9351 1272672000 2592000 2010 2010 1 2
## 6 2010-06-01 10552 1275350400 2678400 2010 2010 1 2
## 7 2010-07-01 9077 1277942400 2592000 2010 2010 2 3
## 8 2010-08-01 9273 1280620800 2678400 2010 2010 2 3
## 9 2010-09-01 9420 1283299200 2678400 2010 2010 2 3
## 10 2010-10-01 9413 1285891200 2592000 2010 2010 2 4
## # ... with 74 more rows, and 22 more variables: month <int>,
## # month.xts <int>, month.lbl <ord>, day <int>, hour <int>,
## # minute <int>, second <int>, hour12 <int>, am.pm <int>,
## # wday <int>, wday.xts <int>, wday.lbl <ord>, mday <int>,
## # qday <int>, yday <int>, mweek <int>, week <int>, week.iso <int>,
## # week2 <int>, week3 <int>, week4 <int>, mday7 <int>
STEP 2：模型
任何回歸模型都可以應用於數據，我們在這里使用 lm()。 請注意，我們刪除了 date 和 diff 列。大多數算法無法使用日期數據，而 diff 列對機器學習沒有什么用處（它對於查找數據中的時間間隔更有用）。
# linear regression model used, but can use any model
fit_lm <- lm(
price ~ .,
data = select(
beer_sales_tbl_aug,
-c(date, diff)))
summary(fit_lm)
##
## Call:
## lm(formula = price ~ ., data = select(beer_sales_tbl_aug, -c(date,
## diff)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -447.3 -145.4 -18.2 169.8 421.4
##
## Coefficients: (16 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.660e+08 1.245e+08 2.940 0.004738 **
## index.num 5.900e-03 2.003e-03 2.946 0.004661 **
## year -1.974e+05 6.221e+04 -3.173 0.002434 **
## year.iso 1.159e+04 6.546e+03 1.770 0.082006 .
## half -2.132e+03 6.107e+02 -3.491 0.000935 ***
## quarter -1.239e+04 2.190e+04 -0.566 0.573919
## month -3.910e+03 7.355e+03 -0.532 0.597058
## month.xts NA NA NA NA
## month.lbl.L NA NA NA NA
## month.lbl.Q -1.643e+03 2.069e+02 -7.942 8.59e-11 ***
## month.lbl.C 8.368e+02 5.139e+02 1.628 0.108949
## month.lbl^4 6.452e+02 1.344e+02 4.801 1.18e-05 ***
## month.lbl^5 7.563e+02 4.241e+02 1.783 0.079852 .
## month.lbl^6 3.206e+02 1.609e+02 1.992 0.051135 .
## month.lbl^7 -3.537e+02 1.867e+02 -1.894 0.063263 .
## month.lbl^8 3.687e+02 3.217e+02 1.146 0.256510
## month.lbl^9 NA NA NA NA
## month.lbl^10 6.339e+02 2.240e+02 2.830 0.006414 **
## month.lbl^11 NA NA NA NA
## day NA NA NA NA
## hour NA NA NA NA
## minute NA NA NA NA
## second NA NA NA NA
## hour12 NA NA NA NA
## am.pm NA NA NA NA
## wday -8.264e+01 1.898e+01 -4.353 5.63e-05 ***
## wday.xts NA NA NA NA
## wday.lbl.L NA NA NA NA
## wday.lbl.Q -7.109e+02 1.093e+02 -6.503 2.13e-08 ***
## wday.lbl.C 2.355e+02 1.336e+02 1.763 0.083273 .
## wday.lbl^4 8.033e+01 1.133e+02 0.709 0.481281
## wday.lbl^5 6.480e+01 8.029e+01 0.807 0.422951
## wday.lbl^6 2.276e+01 8.200e+01 0.278 0.782319
## mday NA NA NA NA
## qday -1.362e+02 2.418e+02 -0.563 0.575326
## yday -2.356e+02 1.416e+02 -1.664 0.101627
## mweek -1.670e+02 1.477e+02 -1.131 0.262923
## week -1.764e+02 1.890e+02 -0.933 0.354618
## week.iso 2.315e+02 1.256e+02 1.842 0.070613 .
## week2 -1.235e+02 1.547e+02 -0.798 0.428283
## week3 NA NA NA NA
## week4 NA NA NA NA
## mday7 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 260.4 on 57 degrees of freedom
## Multiple R-squared: 0.9798, Adjusted R-squared: 0.9706
## F-statistic: 106.4 on 26 and 57 DF, p-value: < 2.2e-16
STEP 3：構建未來（新）數據
使用 tk_index() 來擴展索引。
# Retrieves the timestamp information
beer_sales_idx <- beer_sales_tbl %>%
tk_index()
tail(beer_sales_idx)
## [1] "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01" "2016-11-01"
## [6] "2016-12-01"
通過 tk_make_future_timeseries 函數從現有索引構造未來索引。函數會在內部檢查索引的周期性，並返回正確的序列。我們在 whole vignette on how to make future time series 介紹了該主題更詳盡的內容。
# Make future index
future_idx <- beer_sales_idx %>%
tk_make_future_timeseries(
n_future = 12)
future_idx
## [1] "2017-01-01" "2017-02-01" "2017-03-01" "2017-04-01" "2017-05-01"
## [6] "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01"
## [11] "2017-11-01" "2017-12-01"
用 tk_get_timeseries_signature() 將未來索引轉換成時間序列簽名數據框。
new_data_tbl <- future_idx %>%
tk_get_timeseries_signature()
new_data_tbl
## # A tibble: 12 x 29
## index index.num diff year year.iso half quarter month
## <date> <int> <int> <int> <int> <int> <int> <int>
## 1 2017-01-01 1483228800 NA 2017 2016 1 1 1
## 2 2017-02-01 1485907200 2678400 2017 2017 1 1 2
## 3 2017-03-01 1488326400 2419200 2017 2017 1 1 3
## 4 2017-04-01 1491004800 2678400 2017 2017 1 2 4
## 5 2017-05-01 1493596800 2592000 2017 2017 1 2 5
## 6 2017-06-01 1496275200 2678400 2017 2017 1 2 6
## 7 2017-07-01 1498867200 2592000 2017 2017 2 3 7
## 8 2017-08-01 1501545600 2678400 2017 2017 2 3 8
## 9 2017-09-01 1504224000 2678400 2017 2017 2 3 9
## 10 2017-10-01 1506816000 2592000 2017 2017 2 4 10
## 11 2017-11-01 1509494400 2678400 2017 2017 2 4 11
## 12 2017-12-01 1512086400 2592000 2017 2017 2 4 12
## # ... with 21 more variables: month.xts <int>, month.lbl <ord>,
## # day <int>, hour <int>, minute <int>, second <int>, hour12 <int>,
## # am.pm <int>, wday <int>, wday.xts <int>, wday.lbl <ord>,
## # mday <int>, qday <int>, yday <int>, mweek <int>, week <int>,
## # week.iso <int>, week2 <int>, week3 <int>, week4 <int>,
## # mday7 <int>
STEP 4：預測新數據
將 predict() 應用於回歸模型。注意，和之前使用 lm() 函數時一樣，去掉 index 和 diff 列。
# Make predictions
pred <- predict(
fit_lm,
newdata = select(
new_data_tbl, -c(index, diff)))
predictions_tbl <- tibble(
date = future_idx,
value = pred)
predictions_tbl
## # A tibble: 12 x 2
## date value
## <date> <dbl>
## 1 2017-01-01 9509.122
## 2 2017-02-01 10909.189
## 3 2017-03-01 12281.835
## 4 2017-04-01 11378.678
## 5 2017-05-01 13080.710
## 6 2017-06-01 13583.471
## 7 2017-07-01 11529.085
## 8 2017-08-01 12984.939
## 9 2017-09-01 11859.998
## 10 2017-10-01 12331.419
## 11 2017-11-01 13047.033
## 12 2017-12-01 13940.003
STEP 5：比較實際值與預測值
我們可以使用 tq_get() 來檢索實際數據。注意，我們沒有用於比較的完整數據，但我們至少可以比較前幾個月的實際值。
actuals_tbl <- tq_get(
"S4248SM144NCEN",
get = "economic.data",
from = "2017-01-01",
to = "2017-12-31")
可視化我們的預測。
# Plot Beer Sales Forecast
beer_sales_tbl %>%
ggplot(aes(x = date, y = price)) +
# Training data
geom_line(color = palette_light()[[1]]) +
geom_point(color = palette_light()[[1]]) +
# Predictions
geom_line(
aes(y = value),
color = palette_light()[[2]],
data = predictions_tbl) +
geom_point(
aes(y = value),
color = palette_light()[[2]],
data = predictions_tbl) +
# Actuals
geom_line(
color = palette_light()[[1]],
data = actuals_tbl) +
geom_point(
color = palette_light()[[1]],
data = actuals_tbl) +
# Aesthetics
theme_tq() +
labs(
title = "Beer Sales Forecast: Time Series Machine Learning",
subtitle = "Using basic multivariate linear regression can yield accurate results")
我們可以檢查測試集上的錯誤（實際值 vs 預測值）。
# Investigate test error
error_tbl <- left_join(
actuals_tbl, predictions_tbl) %>%
rename(
actual = price, pred = value) %>%
mutate(
error = actual - pred,
error_pct = error / actual)
error_tbl
## # A tibble: 8 x 5
## date actual pred error error_pct
## <date> <int> <dbl> <dbl> <dbl>
## 1 2017-01-01 8664 9509.122 -845.1223 -0.097544127
## 2 2017-02-01 10017 10909.189 -892.1891 -0.089067495
## 3 2017-03-01 11960 12281.835 -321.8352 -0.026909301
## 4 2017-04-01 11019 11378.678 -359.6777 -0.032641592
## 5 2017-05-01 12971 13080.710 -109.7099 -0.008458092
## 6 2017-06-01 14113 13583.471 529.5292 0.037520667
## 7 2017-07-01 10928 11529.085 -601.0854 -0.055004156
## 8 2017-08-01 12788 12984.939 -196.9386 -0.015400265
接着，我們可以做一些誤差度量。預測值和實際值的 MAPE（平均絕對百分誤差）接近 4.5%，對於簡單的多元線性回歸模型來說是相當好的結果。更復雜的算法可以產生更精確的結果。
# Calculating test error metrics
test_residuals <- error_tbl$error
test_error_pct <- error_tbl$error_pct * 100 # Percentage error
me <- mean(test_residuals, na.rm=TRUE)
rmse <- mean(test_residuals^2, na.rm=TRUE)^0.5
mae <- mean(abs(test_residuals), na.rm=TRUE)
mape <- mean(abs(test_error_pct), na.rm=TRUE)
mpe <- mean(test_error_pct, na.rm=TRUE)
tibble(me, rmse, mae, mape, mpe) %>% glimpse()
## Observations: 1
## Variables: 5
## $ me <dbl> -349.6286
## $ rmse <dbl> 551.7818
## $ mae <dbl> 482.0109
## $ mape <dbl> 4.531821
## $ mpe <dbl> -3.593805
時間序列機器學習可能產生異常預測。對這個議題感興趣的人可以閱讀我們的短文 time series forecasting using timetk。
PART 2：轉換
- 問題：R 中不同類型的時間序列數據難以方便一致地實現相互轉換。
- 解決方案：
tk_tbl、tk_xts、tk_zoo、tk_ts
tk_xts
我們開始的時候用 tbl 對象，一個劣勢是有時候必須轉換成 xts 對象，因為要使用其他包（例如 xts、zoo、quantmod 等等）里基於 xts 對象的函數。
我們可以使用 tk_xts() 函數輕松地將數據轉換成 xts 對象。注意，tk_xts() 函數會自動檢測包含時間的列，並把該列當做 xts 對象的索引。
# Coerce to xts
beer_sales_xts <- tk_xts(beer_sales_tbl)
# Show the first six rows of the xts object
beer_sales_xts %>%
head()
## price
## 2010-01-01 6558
## 2010-02-01 7481
## 2010-03-01 9475
## 2010-04-01 9424
## 2010-05-01 9351
## 2010-06-01 10552
我們也可以從 xts 轉成 tbl 。我們設定 rename_index = "date" 讓索引的名字和開始的時候保持一致。這種操作在以前不太容易。
tk_tbl(beer_sales_xts, rename_index = "date")
## # A tibble: 84 x 2
## date price
## <date> <int>
## 1 2010-01-01 6558
## 2 2010-02-01 7481
## 3 2010-03-01 9475
## 4 2010-04-01 9424
## 5 2010-05-01 9351
## 6 2010-06-01 10552
## 7 2010-07-01 9077
## 8 2010-08-01 9273
## 9 2010-09-01 9420
## 10 2010-10-01 9413
## # ... with 74 more rows
tk_ts
有許多包用了另一種類型的時間序列數據——ts，其中最常見的可能就是 forecast 包。使用 tk_ts() 函數的優點有兩個：
- 與其他
tk_函數兼容，可以方便直接的實現數據的轉換和逆轉換。 - 更重要的是：當使用
tk_ts函數時，ts對象將初始的不規則時間索引（通常是具體的日期）轉換成一個索引屬性。這使得保留日期和時間信息成為可能。
一個例子，使用 tk_tbl(timetk_index = TRUE)函數轉換成 ts 對象。因為 ts 對象只能處理規則時間索引，我們需要添加參數 start = 2010 和 freq = 12。
# Coerce to ts
beer_sales_ts <- tk_ts(
beer_sales_tbl,
start = 2010,
freq = 12)
# Show the calendar-printout
beer_sales_ts
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2010 6558 7481 9475 9424 9351 10552 9077 9273 9420 9413
## 2011 6901 8014 9833 9281 9967 11344 9106 10468 10085 9612
## 2012 7486 8641 9709 9423 11342 11274 9845 11163 9532 10754
## 2013 8395 8888 10109 10493 12217 11385 11186 11462 10494 11541
## 2014 8559 9061 10058 10979 11794 11906 10966 10981 10827 11815
## 2015 8398 9061 10720 11105 11505 12903 11866 11223 12023 11986
## 2016 8540 10158 11879 11155 11916 13291 10540 12212 11786 11424
## Nov Dec
## 2010 9866 11455
## 2011 10328 11483
## 2012 10953 11922
## 2013 11139 12709
## 2014 10466 13303
## 2015 11510 14190
## 2016 12482 13832
有兩種方法轉換回 tbl：
- 直接使用
tk_tbl()，我們將得到YEARMON類型的“規則”的時間索引（來自zoo包）。 - 如果原始對象用
tk_ts()創建，並且有屬性timetk_index，我們可以通過命令tk_tbl(timetk_index = TRUE)轉換回去，並得到Date格式 “非規則”時間索引。
方法 1：注意，日期列是 YEARMON 類型的。
tk_tbl(beer_sales_ts, rename_index = "date")
## # A tibble: 84 x 2
## date price
## <S3: yearmon> <int>
## 1 Jan 2010 6558
## 2 Feb 2010 7481
## 3 Mar 2010 9475
## 4 Apr 2010 9424
## 5 May 2010 9351
## 6 Jun 2010 10552
## 7 Jul 2010 9077
## 8 Aug 2010 9273
## 9 Sep 2010 9420
## 10 Oct 2010 9413
## # ... with 74 more rows
方法 2：設置參數 timetk_idx = TRUE 來找回初始的日期或時間信息。
首先，用 has_timetk_idx() 檢查 ts 對象是否存在 timetk 索引。
# Check for timetk index.
has_timetk_idx(beer_sales_ts)
## [1] TRUE
如果返回值是 TRUE，在調用 tk_tbl() 時設定 timetk_idx = TRUE。現在可以看到日期列是 date 類型的，這在以往不容易做到。
# If timetk_idx is present, can get original dates back
tk_tbl(beer_sales_ts, timetk_idx = TRUE, rename_index = "date")
## # A tibble: 84 x 2
## date price
## <date> <int>
## 1 2010-01-01 6558
## 2 2010-02-01 7481
## 3 2010-03-01 9475
## 4 2010-04-01 9424
## 5 2010-05-01 9351
## 6 2010-06-01 10552
## 7 2010-07-01 9077
## 8 2010-08-01 9273
## 9 2010-09-01 9420
## 10 2010-10-01 9413
## # ... with 74 more rows