Dự báo nhiều chuỗi thời gian (phần 2)
Nguyễn Ngọc Bình
2022-03-03
Phân tích chuỗi thời gian là một phương pháp được sử dụng rộng rãi trong kinh doanh để có được những thông tin hữu ích như dự báo nhu cầu, dự báo giá cổ phiếu, giá tiền mã hóa và các phân tích dự báo khác.
Trong phần này, tôi sẽ giải thích cách dự báo hàng nghìn chuỗi (1000+) sử dụng mô hình Thống kê / Máy học / Học sâu để dự đoán giá trị trong tương lai) & phân loại mô hình nhu cầu (phân loại sản phẩm dựa trên số lượng và thời gian).
Chuỗi phân tích này bao gồm 5 phần:
Phần 1: Làm sạch dữ liệu và tạo features đơn giản.
Phần 2: Mô hình Chuỗi thời gian với phương pháp thống kê (ARIMA, ETS, CROSTON, v.v.) sử dụng gói fable R.
Phần 3: Time Series Feature Engineering sử dụng gói timetk R.
Phần 4: Các mô hình Machine Learning (XGBoost, Random Forest, v.v.) & điều chỉnh Hyperparameter bằng cách sử dụng các gói modeltime & tidymodels R.
Phần 5: Mô hình Deeplearning (NBeats & DeepAR) & điều chỉnh Hyperparameter bằng cách sử dụng các gói modeltime, modeltime.gluonts R.
# install.packages("devtools")
# devtools::install_github("nguyenngocbinh/vndirect")
pacman::p_load(knitr,
vndirect,
tidyverse,
magrittr,
hrbrthemes,
# Tidy Temporal Data Frames and Tools
tsibble,
# Feature Extraction and Statistics for Time Series
feasts,
fable,
fable.prophet)1 Dữ liệu
Sử dụng dữ liệu giá cổ phiếu nguồn vndirect. Trong phần này sẽ sử dụng chuỗi số liệu của 3 cổ phiếu TCB, LPB, BCC
# get data from vndirect source
rawdata <- vndirect::vnd_get_list_data(c('TCB', 'LPB', 'BCC'), 1000)Kiểm tra số lượng quan sát từng chuỗi
rawdata$ticker_name %>%
table() ## .
## BCC LPB TCB
## 1000 1000 9392 Tiền xử lý dữ liệu
Sử dụng giá đóng cửa hàng ngày để phân tích và dự báo. gói tsibble cho phép định dạng nhiều chuỗi trên cùng 1 bảng (tương tự như group_by trên dplyr). Bằng việc chọn key và index phù hợp ta sẽ định dạng được bảng dữ liệu gồm nhiều chuỗi thời gian
# convert to tsibble format
dt <- rawdata %>%
# select target variable
rename(value = adClose) %>%
# convert to tsibble format
as_tsibble(key = ticker_name, index = date) %>%
fill_gaps()Xử lý missing nếu có trong dữ liệu. Do fable yêu cầu fill_gaps, tức là điền các khoảng thời gian thiếu vào chuỗi. Vì vậy, sẽ có nhiều ngày không có dữ liệu (ví dụ, thứ 7, chủ nhật không có giao dịch). Phần dưới sẽ điền những giá trị missing này bằng giá trị ngày gần nhất.
fnc_na_vars <- function(x) {any(is.na(x))}
na_vars <- dt %>%
select_if(fnc_na_vars) %>%
names()
dt <- dt %>%
fill(all_of(na_vars), .direction = 'down')Chia dữ liệu thành nhiều tập train (fable không hỗ trợ chia train, test), sau đó sử dụng chuỗi dữ liệu ban đầu để đánh giá performance của mô hình.
# split data to cross validation
cv_dt <- dt %>%
stretch_tsibble(.init = 370, .step = 100) 3 Modelling
3.1 Các loại mô hình
Các loại mô hình được tham khảo từ Forecasting: Principles and Practice
3.1.1 Naive
Mô hình Naive gán tất cả các giá trị trong tương lai bằng giá trị quan sát cuối cùng.
3.1.2 Seasonal Naive
Mô hình Naive theo mùa được sử dụng cho dữ liệu theo mùa. Trong tình huống này, mỗi giá trị tương lai bằng giá trị cuối cùng của cùng một mùa.
3.1.3 Drift
Mô hình drift là một biến thể của mô hình Naive, cho phép dự báo tăng giảm theo thời gian.
3.1.4 Simple Exponential Smooth (SES)
Mô hình trơn theo cấp số nhân đơn giản được sử dụng khi không thể xác định được xu hướng rõ ràng hoặc mô hình theo mùa.
\[\begin{align*} \text{Forecast equation} && \hat{y}_{t+h|t} & = \ell_{t}\\ \text{Smoothing equation} && \ell_{t} & = \alpha y_{t} + (1 - \alpha)\ell_{t-1}, \end{align*}\]
Trong đó:
- \(\ell_{t}\): giá trị làm trơn tại thời gian t
- h: thời gian dự báo ngoài mẫu
- \(\alpha\): tham số làm trơn \(0\le\alpha\le1\)
3.1.5 Holt’s Linear
Mô hình Holt’s Linear là một phiên bản mở rộng của SES cho phép dự báo theo xu hướng.
\[\begin{align*} \text{Forecast equation}&& \hat{y}_{t+h|t} &= \ell_{t} + hb_{t} \\ \text{Level equation} && \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\ \text{Trend equation} && b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1}, \end{align*}\]
- \(b_t\): ước lượng xu hướng hay độ dốc tại thời điểm t
- \(\beta^*\): tham số làm trơn $ 0^* $
3.1.6 Damped Holt’s Linear
Là phiên bản mở rộng của Holt’s Linear. Mô hình Holt’s Linear cho ra chuỗi liên tục tăng hoặc liên tục giảm vô hạn trong tương lai vì vậy có thể dự báo quá mức (đặc biệt với chuỗi thời gian dài). Gardner & McKenzie (1985) đưa ra 1 tham số điều chỉnh xu hướng, cho phép xu hướng thay đổi theo thời gian với một tham số \(\phi\).
\[\begin{align*} \hat{y}_{t+h|t} &= \ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t} \\ \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\ b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)\phi b_{t-1}. \end{align*}\]
- \(0<\phi<1\), \(\phi=1\) ~ Holt’s Linear
3.1.7 Forecasting with Decomposition
Mô hình này được sử dụng để phân rã chuỗi thời gian theo mùa và theo xu hướng bằng cách sử dụng phương pháp lowess (Seasonal and Trend decomposition using Loess - STL). Sau đó, sử dụng ETS để dự báo dữ liệu được điều chỉnh theo mùa, và thành phần theo mùa có thể được thêm vào dữ liệu dự báo.
- Dạng mô hình cộng:
\[y_{t} = S_{t} + T_{t} + R_t\]
- Dạng mô hình nhân:
\[y_{t} = S_{t} \times T_{t} \times R_t\]
3.1.8 ARIMA
Mô hình ARIMA cung cấp một cách tiếp cận khác để dự báo chuỗi thời gian. Các mô hình ARIMA và ETS là hai cách tiếp cận được sử dụng rộng rãi nhất để dự báo chuỗi thời gian. Trong khi các mô hình ETS dựa trên đặc điểm về xu hướng và tính thời vụ trong dữ liệu, thì các mô hình ARIMA hướng tới tự tương quan trong dữ liệu. Ở đây, chúng tôi sẽ chọn mô hình ARIMA tự động bằng cách sử dụng AICC tối thiểu.
3.1.9 Timeseries Regression
Mô hình hồi quy, cho phép sử dụng biến độc lập khác.
\[ \begin{equation} y_t = \beta_{0} + \beta_{1} x_{1,t} + \beta_{2} x_{2,t} + \cdots + \beta_{k} x_{k,t} + \varepsilon_t, \tag{7.1} \end{equation} \]
3.1.10 Dynamic Harmonic Regression
Các phiên bản theo mùa của mô hình ARIMA và ETS được thiết kế cho các khoảng thời gian ngắn, chẳng hạn như 12 cho dữ liệu hàng tháng hoặc 4 cho dữ liệu hàng quý. Mô hình ETS() giới hạn thời vụ là khoảng thời gian tối đa là 24 để cho phép dữ liệu hàng giờ nhưng không cho phép dữ liệu có khoảng thời gian theo mùa lớn hơn. Hàm ARIMA () sẽ cho phép khoảng thời gian theo mùa lên đến m =350, nhưng trong thực tế, thường sẽ hết bộ nhớ bất cứ khi nào khoảng thời gian theo mùa nhiều hơn 200. Vì vậy, đối với chuỗi thời gian mà chia nhỏ đến phút hoặc giây thì 2 mô hình này không đáp ứng được.
Trong trường hợp này, sử dụng phương pháp hồi quy hài hòa trong đó mô hình theo mùa được mô hình hóa bằng cách sử dụng thuật ngữ Fourier với chuỗi thời gian ngắn hạn được xử lý bằng ARMA.
Đối với dữ liệu có nhiều hơn một khoảng thời gian theo mùa, có thể sử dụng Fourier với các mức khác nhau; độ mượt của mô hình theo mùa có thể được kiểm soát bởi K, bao gồm các cặp Fourier sin và cos - mô hình theo mùa mượt hơn đối với các giá trị nhỏ hơn của K. Chuỗi thời gian ngắn hạn dễ dàng được xử lý với mô hình ARMA đơn giản.
3.1.11 Ensemble
Mô hình tổng hợp chỉ đơn giản là sử dụng một số mô hình khác nhau cùng một lúc và tính toán giá trị trung bình của các dự báo kết quả. Ví dụ: ở đây, các mô hình ARIMA, SES & Decomposition được sử dụng cùng nhau để tính toán giá trị dự báo trung bình.
model_single_var_table <- cv_dt %>%
model(
## Model 1: Naive ----
naive_mod = NAIVE(log(value + 1)),
## Model 2: Snaive ----
snaive_mod = SNAIVE(log(value + 1)),
## Model 3: Drift ----
drift_mod = RW(log(value + 1) ~ drift()),
## Model 4: SES ----
ses_mod = ETS(log(value + 1) ~ error("A") + trend("N") + season("N"),
opt_crit = "mse"),
## Model 5: Holt's Linear ----
hl_mod = ETS(log(value + 1) ~ error("A") + trend("A") + season("N"),
opt_crit = "mse"),
## Model 6: Damped Holt's Linear ----
hldamp_mod = ETS(log(value + 1) ~ error("A") + trend("Ad") + season("N"),
opt_crit = "mse"),
## Model 7: STL decomposition with ETS ----
stl_ets_mod = decomposition_model(
STL(log(value + 1), ~ season(window = 7)),
ETS(season_adjust ~ season("N")),
NAIVE(season_week),
SNAIVE(season_year)
),
## Model 8: ARIMA ----
arima_mod = ARIMA(log(value + 1)),
## Model 9: Dynamic harmonic regression ----
dhr_mod = ARIMA(
log(value + 1) ~ PDQ(0, 0, 0) + fourier(K = 3)+ lag(basicPrice) + lag(high)
),
## Model 10: TSLM ----
tslm_mod = TSLM(
log(value + 1) ~ lag(basicPrice) + lag(high) + lag(low) + lag(nmVolume)
)
) %>%
## Model 99: Ensemble Model ----
mutate(ensemble_sm_mod = combination_ensemble(arima_mod, ses_mod, stl_ets_mod)) ## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning in sqrt(diag(best$var.coef)): NaNs produced## Warning: 12 errors (1 unique) encountered for stl_ets_mod
## [12] object 'season_year' not found4 Đánh giá khả năng dự báo
4.1 Training set accuracy
Do mô hình kết hợp không sử dụng cùng hàm accuracy(). Vì vậy, tôi viết lại hàm kết hợp các giá trị accuracy
comb_accuracy <- function(.model){
train_acc_single1 <- .model %>%
select(-ensemble_sm_mod) %>% # not include ensemble_sm_mod
accuracy()
train_acc_single2 <- .model %>%
select(ensemble_sm_mod) %>%
accuracy()
train_acc_single <- bind_rows(train_acc_single1, train_acc_single2)
return(train_acc_single)
}training_acc <- comb_accuracy(model_single_var_table) %>%
select(.id, ticker_name, .model, .type, RMSE, MAE, MAPE)
training_acc %>% head(15) %>% kable()| .id | ticker_name | .model | .type | RMSE | MAE | MAPE |
|---|---|---|---|---|---|---|
| 1 | BCC | naive_mod | Training | 0.1195488 | 0.0647371 | 1.063419 |
| 1 | BCC | snaive_mod | Training | 0.2839378 | 0.1958926 | 3.194357 |
| 1 | BCC | drift_mod | Training | 0.1195529 | 0.0652933 | 1.072899 |
| 1 | BCC | ses_mod | Training | 0.1192045 | 0.0665215 | 1.092302 |
| 1 | BCC | hl_mod | Training | 0.1192149 | 0.0667807 | 1.096848 |
| 1 | BCC | hldamp_mod | Training | 0.1189244 | 0.0667926 | 1.097196 |
| 1 | BCC | stl_ets_mod | Training | NaN | NaN | NaN |
| 1 | BCC | arima_mod | Training | 0.1192065 | 0.0665556 | 1.092881 |
| 1 | BCC | dhr_mod | Training | 0.1204791 | 0.0802047 | 1.320086 |
| 1 | BCC | tslm_mod | Training | 0.1306086 | 0.0885545 | 1.454254 |
| 1 | LPB | naive_mod | Training | 0.1643182 | 0.0915854 | 1.097041 |
| 1 | LPB | snaive_mod | Training | 0.4520146 | 0.3216777 | 3.806460 |
| 1 | LPB | drift_mod | Training | 0.1637109 | 0.0953367 | 1.145566 |
| 1 | LPB | ses_mod | Training | 0.1638137 | 0.0927775 | 1.111481 |
| 1 | LPB | hl_mod | Training | 0.1631510 | 0.0953328 | 1.144763 |
4.2 TSCV accuracy
Do mô hình Dynamic Harmonic Regression và Timeseries Regression không có dữ liệu trong tương lai, nên sẽ không forecast để tính accuracy
test_acc_single1 <- model_single_var_table %>%
select(-ensemble_sm_mod, -dhr_mod, -tslm_mod) %>% # not include ensemble_sm_mod, dhr_mod, tslm_mod
forecast(h = "3 months") %>%
accuracy(dt)## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 81 observations are missing between 2022-03-04 and 2022-05-23test_acc_single2 <- model_single_var_table %>%
select(ensemble_sm_mod) %>%
forecast(h = "3 months") %>%
accuracy(dt)## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 81 observations are missing between 2022-03-04 and 2022-05-23test_acc_single <- test_acc_single1 %>%
bind_rows(test_acc_single2) %>%
select(ticker_name, .model, .type, RMSE, MAE, MAPE)4.3 Mô hình tốt nhất theo MAPE
best_mod_training_acc <- training_acc %>%
group_by(ticker_name, .type) %>%
summarise(MAPE = min(MAPE, na.rm = T)) %>%
inner_join(training_acc)## `summarise()` has grouped output by 'ticker_name'. You can override using the
## `.groups` argument.
## Joining, by = c("ticker_name", ".type", "MAPE")best_mod_test_acc <- test_acc_single %>%
group_by(ticker_name, .type) %>%
summarise(MAPE = min(MAPE, na.rm = T)) %>%
inner_join(test_acc_single)## `summarise()` has grouped output by 'ticker_name'. You can override using the
## `.groups` argument.
## Joining, by = c("ticker_name", ".type", "MAPE")best_mod_acc <- best_mod_training_acc %>% bind_rows(best_mod_test_acc)
best_mod_acc %>% kable()| ticker_name | .type | MAPE | .id | .model | RMSE | MAE |
|---|---|---|---|---|---|---|
| BCC | Training | 1.0418665 | 4 | naive_mod | 0.1215965 | 0.0677534 |
| LPB | Training | 0.8977130 | 4 | naive_mod | 0.1313328 | 0.0682048 |
| TCB | Training | 0.8368715 | 3 | ses_mod | 0.4066599 | 0.2173702 |
| BCC | Test | 13.7929232 | NA | naive_mod | 2.9790262 | 1.7796963 |
| LPB | Test | 8.6506065 | NA | snaive_mod | 2.4714792 | 1.3628243 |
| TCB | Test | 8.7966750 | NA | snaive_mod | 4.0190035 | 2.8563180 |
Hơi bất ngờ khi các mô hình Naive lại có performance tốt nhất
5 Dự báo
5.1 Fit lại mô hình với toàn bộ chuỗi dữ liệu
final_model_single_var_table <- dt %>%
model(
## Model 1: Naive ----
naive_mod = NAIVE(log(value + 1)),
## Model 2: Snaive ----
snaive_mod = SNAIVE(log(value + 1)),
## Model 3: Drift ----
drift_mod = RW(log(value + 1) ~ drift()),
## Model 4: SES ----
ses_mod = ETS(log(value + 1) ~ error("A") + trend("N") + season("N"),
opt_crit = "mse"),
## Model 5: Holt's Linear ----
hl_mod = ETS(log(value + 1) ~ error("A") + trend("A") + season("N"),
opt_crit = "mse"),
## Model 6: Damped Holt's Linear ----
hldamp_mod = ETS(log(value + 1) ~ error("A") + trend("Ad") + season("N"),
opt_crit = "mse"),
## Model 7: STL decomposition with ETS ----
stl_ets_mod = decomposition_model(
STL(log(value + 1), ~ season(window = 7)),
ETS(season_adjust ~ season("N")),
NAIVE(season_week),
SNAIVE(season_year)
),
## Model 8: ARIMA ----
arima_mod = ARIMA(log(value + 1)),
## Model 9: Dynamic harmonic regression ----
dhr_mod = ARIMA(
log(value + 1) ~ PDQ(0, 0, 0) + fourier(K = 3)+ lag(basicPrice) + lag(high)
),
## Model 10: TSLM ----
tslm_mod = TSLM(
log(value + 1) ~ lag(basicPrice) + lag(high) + lag(low) + lag(nmVolume)
)
) %>%
## Model 99: Ensemble Model ----
mutate(ensemble_sm_mod = combination_ensemble(arima_mod, ses_mod, stl_ets_mod)) 5.2 Dự báo 3 tháng tiếp theo
final_fc <- final_model_single_var_table %>%
select(-dhr_mod, -tslm_mod) %>%
forecast(h = "3 months")final_fc %>%
filter(ticker_name =='BCC') %>%
autoplot(dt)+
theme_ipsum()+
labs(title = 'BCC')
final_fc %>%
filter(ticker_name =='LPB' & .model == 'snaive_mod') %>%
autoplot(dt)+
theme_ipsum()+
labs(title = 'LPB')
final_fc %>%
filter(ticker_name =='TCB' & .model == 'snaive_mod') %>%
autoplot(dt)+
theme_ipsum()+
labs(title = 'TCB')
5.3 Kết quả dự báo với mô hình có MAPE tốt nhất
best_fc <- best_mod_acc %>%
filter(.type == 'Test') %>%
inner_join(final_fc) ## Joining, by = c("ticker_name", ".model")