[TOC] 翻譯自《Demo Week: Tidy Forecasting with sweep》 原文鏈接:www.business science.io/code tools/2017/10/25/demo_week_sweep.html 時間序列分析工具箱——sweep 的用途 正如 包之於 ...
目錄
翻譯自《Demo Week: Tidy Forecasting with sweep》
原文鏈接:www.business-science.io/code-tools/2017/10/25/demo_week_sweep.html
時間序列分析工具箱——sweep
sweep
的用途
正如 broom
包之於 stats
包,sweep
包用來簡化使用 forecast
包的工作流。本教程將逐一介紹常用函數 sw_tidy
、sw_glance
、sw_augment
和 sw_sweep
的用法。
sweep
和 timetk
帶來的額外好處是,如果 ts
對象是由 tbl
對象轉換來的,那麼在預測過程中日期和時間信息會以 timetk
索引的形式保留下來。一句話概括:這意味著我們最終可以在預測時使用日期,而不是 ts
類型數據使用的規則間隔數字日期。
載入包
本教程要使用到四個包:
sweep
:簡化forecast
包的使用forecast
:提供 ARIMA、ETS 和其他流行的預測演算法tidyquant
:獲取數據併在後臺載入tidyverse
系列工具timetk
:時間序列數據處理工具,用來將tbl
轉換成ts
# Load libraries
library(sweep) # Broom-style tidiers for the forecast package
library(forecast) # Forecasting models and predictions package
library(tidyquant) # Loads tidyverse, financial pkgs, used to get data
library(timetk) # Functions working with time series
數據
我們使用 timetk
教程中數據——啤酒、紅酒和蒸餾酒銷售數據,用 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")
現在你對我們要分析的時間序列有了直觀的感受,那麼讓我們繼續!
Demo:forecast
+ sweep
的簡化預測工作流
我們將聯合使用 forecast
和 sweep
來簡化預測分析。
關鍵想法:使用 forecast
包做預測涉及到 ts
對象,用起來並不簡潔。對於 stats
包來說有 broom
來簡化使用;forecast
包就用 sweep
。
目標:我們將用 ARIMA 模型預測未來 12 個月的數據。
STEP 1:創建 ts
對象
使用 timetk::tk_ts()
將 tbl
轉換成 ts
,從之前的教程可以瞭解到這個函數有兩點好處:
- 這是一個統一的方法,實現與
ts
對象的相互轉換。 - 得到的
ts
對象包含timetk_idx
屬性,是一個基於初始時間信息的索引。
下麵開始轉換,註意 ts
對象是規則時間序列,所以要設置 start
和 freq
。
# Convert from tbl to ts
beer_sales_ts <- tk_ts(
beer_sales_tbl,
start = 2010,
freq = 12)
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
檢查 ts
對象具有 timetk_idx
屬性。
# Check that ts-object has a timetk index
has_timetk_idx(beer_sales_ts)
## [1] TRUE
OK,這對後面要用的 sw_sweep()
很重要。下麵我們就要建立 ARIMA 模型了。
STEP 2A:ARIMA 模型
我們使用 forecast
包里的 auto.arima()
函數為時間序列建模。
# Model using auto.arima
fit_arima <- auto.arima(beer_sales_ts)
fit_arima
## Series: beer_sales_ts
## ARIMA(3,0,0)(1,1,0)[12] with drift
##
## Coefficients:
## ar1 ar2 ar3 sar1 drift
## -0.2498 0.1079 0.6210 -0.2817 32.1157
## s.e. 0.0933 0.0982 0.0925 0.1333 5.8882
##
## sigma^2 estimated as 175282: log likelihood=-535.49
## AIC=1082.97 AICc=1084.27 BIC=1096.63
STEP 2B:簡化模型
就像 broom
簡化 stats
包的使用一樣,我麽可以使用 sweep
的函數簡化 ARIMA 模型。下麵介紹三個函數:
sw_tidy()
:用於檢索模型參數sw_glance()
:用於檢索模型描述和訓練集的精確度度量sw_augment()
:用於獲得模型殘差
sw_tidy
sw_tidy()
函數以 tibble
對象的形式返回模型參數。
# sw_tidy - Get model coefficients
sw_tidy(fit_arima)
## # A tibble: 5 x 2
## term estimate
## <chr> <dbl>
## 1 ar1 -0.2497937
## 2 ar2 0.1079269
## 3 ar3 0.6210345
## 4 sar1 -0.2816877
## 5 drift 32.1157478
sw_glance
sw_glance()
函數以 tibble
對象的形式返回訓練集的精確度度量。可以使用 glimpse
函數美化顯示結果。
# sw_glance - Get model description and training set accuracy measures
sw_glance(fit_arima) %>%
glimpse()
## Observations: 1
## Variables: 12
## $ model.desc <chr> "ARIMA(3,0,0)(1,1,0)[12] with drift"
## $ sigma <dbl> 418.6665
## $ logLik <dbl> -535.4873
## $ AIC <dbl> 1082.975
## $ BIC <dbl> 1096.635
## $ ME <dbl> 1.189875
## $ RMSE <dbl> 373.9091
## $ MAE <dbl> 271.7068
## $ MPE <dbl> -0.06716239
## $ MAPE <dbl> 2.526077
## $ MASE <dbl> 0.4989005
## $ ACF1 <dbl> 0.02215405
sw_augument
sw_augument()
函數返回的 tibble
表中包含 .actual
、.fitted
和 .resid
列,有助於在訓練集上評估模型表現。註意,設置 timetk_idx = TRUE
返回初始的日期索引。
# sw_augment - get model residuals
sw_augment(fit_arima, timetk_idx = TRUE)
## # A tibble: 84 x 4
## index .actual .fitted .resid
## <date> <dbl> <dbl> <dbl>
## 1 2010-01-01 6558 6551.474 6.525878
## 2 2010-02-01 7481 7473.583 7.416765
## 3 2010-03-01 9475 9465.621 9.378648
## 4 2010-04-01 9424 9414.704 9.295526
## 5 2010-05-01 9351 9341.810 9.190414
## 6 2010-06-01 10552 10541.641 10.359293
## 7 2010-07-01 9077 9068.148 8.852178
## 8 2010-08-01 9273 9263.984 9.016063
## 9 2010-09-01 9420 9410.869 9.130943
## 10 2010-10-01 9413 9403.908 9.091831
## # ... with 74 more rows
我們可以可視化訓練數據上的殘差,看一下數據中有沒有遺漏的模式沒有被髮現。
# Plotting residuals
sw_augment(fit_arima, timetk_idx = TRUE) %>%
ggplot(aes(x = index, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, color = "red") +
labs(title = "Residual diagnostic") +
scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
theme_tq()
STEP 3:預測
使用 forecast()
函數做預測。
# Forecast next 12 months
fcast_arima <- forecast(fit_arima, h = 12)
一個問題是,預測結果並不“tidy”。我們需要數據框形式的預測結果,以便應用 tidyverse
的功能,然而預測結果是 forecast
類型的,一種基於 ts
的對象。
class(fcast_arima)
## [1] "forecast"
STEP 4:用 sweep
簡化預測
我們使用 sw_sweep()
簡化預測結果,一個額外的好處是,如果 forecast
對象有 timetk
索引,我們可以用它返回一個日期時間索引,不同於 ts
對象的規則索引。
首先要確認 forecast
對象有 timetk
索引,這需要在使用 sw_sweep()
時設置 timetk_idx
參數。
# Check if object has timetk index
has_timetk_idx(fcast_arima)
## [1] TRUE
現在,使用 sw_sweep()
來簡化預測結果,它會在內部根據 time_tk
構造一條未來時間序列索引(這一步總是會被執行,因為我們在第 1 步中用 tk_ts()
構造了 ts
對象)註意:這意味著我們最終可以在 forecast
包中使用日期(不同於 ts
對象中的規則索引)!。
# sw_sweep - tidies forecast output
fcast_tbl <- sw_sweep(fcast_arima, timetk_idx = TRUE)
fcast_tbl
## # A tibble: 96 x 7
## index key price lo.80 lo.95 hi.80 hi.95
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010-01-01 actual 6558 NA NA NA NA
## 2 2010-02-01 actual 7481 NA NA NA NA
## 3 2010-03-01 actual 9475 NA NA NA NA
## 4 2010-04-01 actual 9424 NA NA NA NA
## 5 2010-05-01 actual 9351 NA NA NA NA
## 6 2010-06-01 actual 10552 NA NA NA NA
## 7 2010-07-01 actual 9077 NA NA NA NA
## 8 2010-08-01 actual 9273 NA NA NA NA
## 9 2010-09-01 actual 9420 NA NA NA NA
## 10 2010-10-01 actual 9413 NA NA NA NA
## # ... with 86 more rows
STEP 5:比較真實值和預測值
我們可以使用 tq_get()
來檢索實際數據。註意,我們沒有用於比較的完整數據,但我們至少可以比較前幾個月的實際值。
actuals_tbl <- tq_get(
"S4248SM144NCEN",
get = "economic.data",
from = "2017-01-01",
to = "2017-12-31")
註意,預測結果放在 tibble
中,可以方便的實現可視化。
# Visualize the forecast with ggplot
fcast_tbl %>%
ggplot(
aes(x = index, y = price, color = key)) +
# 95% CI
geom_ribbon(
aes(ymin = lo.95, ymax = hi.95),
fill = "#D5DBFF", color = NA, size = 0) +
# 80% CI
geom_ribbon(
aes(ymin = lo.80, ymax = hi.80, fill = key),
fill = "#596DD5", color = NA,
size = 0, alpha = 0.8) +
# Prediction
geom_line() +
geom_point() +
# Actuals
geom_line(
aes(x = date, y = price), color = palette_light()[[1]],
data = actuals_tbl) +
geom_point(
aes(x = date, y = price), color = palette_light()[[1]],
data = actuals_tbl) +
# Aesthetics
labs(
title = "Beer Sales Forecast: ARIMA", x = "", y = "Thousands of Tons",
subtitle = "sw_sweep tidies the auto.arima() forecast output") +
scale_x_date(
date_breaks = "1 year",
date_labels = "%Y") +
scale_color_tq() +
scale_fill_tq() +
theme_tq()
我們可以研究測試集上的誤差(真實值 vs 預測值)。
# Investigate test error
error_tbl <- left_join(
actuals_tbl,
fcast_tbl,
by = c("date" = "index")) %>%
rename(
actual = price.x, pred = price.y) %>%
select(date, actual, pred) %>%
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 8601.815 62.18469 0.007177365
## 2 2017-02-01 10017 10855.429 -838.42908 -0.083700617
## 3 2017-03-01 11960 11502.214 457.78622 0.038276439
## 4 2017-04-01 11019 11582.600 -563.59962 -0.051147982
## 5 2017-05-01 12971 12566.765 404.23491 0.031164514
## 6 2017-06-01 14113 13263.918 849.08191 0.060163106
## 7 2017-07-01 10928 11507.277 -579.27693 -0.053008504
## 8 2017-08-01 12788 12527.278 260.72219 0.020388035
並且,我們可以做簡單的誤差度量。MAPE 接近 4.3%,比簡單的線性回歸模型略好一點,但是 RMSE 變差了。
# Calculate 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> 6.588034
## $ rmse <dbl> 561.4631
## $ mae <dbl> 501.9144
## $ mape <dbl> 4.312832
## $ mpe <dbl> -0.3835956