作者简介：
徐瑞龙，量化分析师
博客专栏：https://www.cnblogs.com/xuruilong100
本文翻译自《Demo Week: Time Series Machine Learning with timetk》
原文链接：
www.business-science.io/code-tools/2017/10/24/demo_week_timetk.html
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-16STEP 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.003STEP 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 rowstk_ts
有许多包用了另一种类型的时间序列数据——ts，其中最常见的可能就是 forecast 包。使用 tk_ts() 函数的优点有两个：
与其他
tk_函数兼容，可以方便直接的实现数据的转换和逆转换。更重要的是：当使用
tk_ts函数时，ts对象将初始的不规则时间索引（通常是具体的日期）转换成一个索引属性。这使得保留日期和时间信息成为可能。
一个例子，使用 tk_tbl(timetk_index = TRUE)函数转换成 ts 对象。因为 ts 对象只能处理规则时间索引，我们需要添加参数 start = 2010 和 freq = 12。
# Coerce to tsbeer_sales_ts <- tk_ts(
beer_sales_tbl,
start = 2010,
freq = 12)# Show the calendar-printoutbeer_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大家都在看
公众号后台回复关键字即可学习
回复 爬虫 爬虫三大案例实战
回复 Python 1小时破冰入门回复 数据挖掘 R语言入门及数据挖掘
回复 人工智能 三个月入门人工智能
回复 数据分析师 数据分析师成长之路
回复 机器学习 机器学习的商业应用
回复 数据科学 数据科学实战
回复 常用算法 常用数据挖掘算法
扫一扫
1377
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
