In nguyenngocbinh/bml_stocks: Provide simple forecast on stocks price
knitr::opts_chunk$set(echo = FALSE,
warning = FALSE,
message = FALSE,
knitr.table.format = "html",
fig.align = "center")
# Ngày 11/11/2019 Sửa lỗi betas
# Ngày 12/1/2019: Sửa lỗi hậu kiểm
library(DT)
library(formattable)
library(fGarch)
library(FinTS)
# library(kableExtra)
# library(TSA) for acf
# library(tseries)
# library(robustbase)
# library(RcppRoll) # roll_mean
# library(zoo)
# library(TTR)
library(readxl)
#library(writexl)
library(rugarch)
library(PerformanceAnalytics)
library(tidyverse)
library(magrittr)
library(lubridate)
library(tidyquant)
rm(list = ls())
######### Input #############
stock_name <- "DHG" #
#############################
# Choose estimation period and validation period
############################## Input #############################
start_estimation_date <- as.Date("2012-12-01") #
end_estimation_date <- start_estimation_date %m+% years(5) - 1 #
start_validation_date <- start_estimation_date %m+% years(5) #
end_validation_date <- start_estimation_date %m+% years(6) -1 #
##################################################################
# format date to export
fm_start_estimation_date <- format(start_estimation_date, "%d/%m/%Y")
fm_end_estimation_date <- format(end_estimation_date, "%d/%m/%Y")
fm_start_validation_date <- format(start_validation_date, "%d/%m/%Y")
fm_end_validation_date <- format(end_validation_date, "%d/%m/%Y")
# http://www.unstarched.net/r-examples/rugarch/a-short-introduction-to-the-rugarch-package/
# legend.title = element_blank(),
# plot.title = element_text(hjust = 0.5),
# plot.subtitle = element_text(hjust = 0.5),
theme_nnb <- theme_tq() +
theme(plot.caption = element_text(hjust = 1),
legend.background = element_blank(),
panel.grid.minor.y = element_blank())
theme_set(theme_nnb)
Tóm lược
Tài liệu được soạn thảo với mục đích cung cấp thông tin chi tiết về kỹ thuật ước lượng và dự báo biến động giá cổ phiếu r stock_name. Tài liệu được cấu trúc theo 5 phần chính:
1. Diễn biến giá cổ phiếu r stock_name: phần này cung cấp cái nhìn tổng quát về diễn biến độ biến động giá cổ phiếu r stock_name.
2. Kiểm định giả thiết: kiểm định mức độ tự tương quan và phương sai sai số thay đổi của chuỗi lợi suất giá cổ phiếu r stock_name.
3. Ước lượng mô hình phù hợp: ước lượng mức giảm giá tối đa bằng phương pháp mô phỏng lịch sử và các mô hình GARCH, trong đó tập trung xác định mô hình GARCH phù hợp trong lớp các mô hình GARCH.
4. Hậu kiểm VaR: hậu kiểm giá trị VaR ước lượng được bằng phương pháp mô phỏng lịch sử và mô hình GARCH từ đó chọn ra phương pháp phù hợp.
5. Ước lượng và dự báo mức giảm giá tối đa: thực hiện ước lượng lại mức giảm giá tối đa bằng phương pháp phù hợp nhất đã lựa chọn được đối với khoảng dữ liệu mới nhất.
Các nội dung chính
-
Độ biến động và mức giảm giá tối đa của mặt hàng cổ phiếu r stock_name được ước lượng dựa trên dữ liệu giá cổ phiếu r stock_name hàng ngày – cophieu68.com trong khoảng thời gian từ r fm_start_estimation_date đến r fm_end_validation_date.
-
Dữ liệu được chia thành hai phần: khoảng thời gian từ r fm_start_estimation_date đến r fm_end_estimation_date được dùng để ước lượng độ biến động và mức giảm giá tối đa VaR95%, khoảng thời gian từ r fm_start_validation_date đến r fm_end_validation_date để hậu kiểm mức độ phù hợp của giá trị VaR95%. Sau khi xác định được phương pháp phù hợp để tính mức giảm giá tối đa, chúng tôi sẽ áp dụng phương pháp đó để ước lượng VaR95% trong khoảng thời gian từ r fm_start_estimation_date đến r fm_end_validation_date từ đó đưa ra dự báo về mức giảm giá tối đa của cổ phiếu r stock_name.
my_folder <- paste0("D:/stock/co_phieu/")
price_raw <- read_csv(paste0(my_folder,"/excel_",stock_name,".csv")) %>%
mutate(time_stamp = ymd(`<DTYYYYMMDD>`),
price = `<CloseFixed>`) %>%
na.omit()
# Change data type
data <- xts(x = price_raw$price, order.by = price_raw$time_stamp)
Diễn biến thực tế của giá cổ phiếu r stock_name
price_raw_filter <- price_raw %>%
filter(time_stamp >= end_validation_date - years(2))
price_raw_filter %>%
mutate(ma30 = rollmean(price, 30, align = "center", fill = NA),
ma7 = rollmean(price, 7, align = "center", fill = NA)) %>%
ggplot() +
geom_point(aes(x = time_stamp, y = price, color = "daily"), size = 1.5, alpha = 0.5) +
geom_line(aes(x = time_stamp, y = ma30, color = "ma30"), size = 1) +
geom_line(aes(x = time_stamp, y = ma7, color = "ma7"), size = 1)+
scale_color_tq()+
#scale_color_manual(values = c("#2c3e50", "#e31a1c", "#18BC9C"))+
scale_x_date(date_labels = "%b %Y") +
labs(x = NULL,
y = "Price",
title = paste("Market price of stock",stock_name),
subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)),
" to ", as.yearmon(end_validation_date)),
caption = "Author: Nguyen Ngoc Binh") +
theme(legend.position = c(0.75, 0.95),
legend.direction = "horizontal",
legend.title = element_blank(),
axis.text.x = element_text(angle = 30, hjust = 1))
Biểu đồ cho thấy giá cổ phiếu r stock_name nằm trong khoảng từ r min(price_raw_filter$price) đến r max(price_raw_filter$price) nghìn đồng. Để thấy rõ độ biến động của giá cổ phiếu r stock_name, chúng tôi sẽ quan sát biểu đồ chuỗi lợi suất (Return).
Tỷ suất lợi nhuận cổ phiếu r stock_name
# Calculate return
data_return <- Return.calculate(data, method="discrete") %>%
na.omit() %>%
set_names("Return") %>%
window(start = start_estimation_date, end = end_validation_date)
ggplot(aes(x = Index, y = Value), data = fortify(data_return, melt = TRUE)) +
geom_line(color = "skyblue3") +
labs(title = paste0("Daily return of stock ", stock_name %>% toupper()),
subtitle=paste0("From ", as.yearmon(start_estimation_date),
" to ", as.yearmon(end_validation_date)),
y = NULL,
x = NULL,
caption = "Author: Nguyen Ngoc Binh")+
scale_x_date(date_labels = "%b %Y") +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
- Độ biến động của giá cổ phiếu
r stock_name được thể hiện rõ hơn qua các giá trị đặc trưng của phân phối chuỗi lợi suất
data_return %>%
table.Stats() %>%
rownames_to_column() %>%
slice(c(1, 3, 6, 7, 9, 13, 15, 16)) %>%
set_names(c("Indicators", "Return")) %>%
datatable(rownames = FALSE, options = list(searching = FALSE))
- Làm sạch dữ liệu
Sử dụng phương pháp làm sạch dữ liệu (Clean data) bằng phương pháp “boudt” để quan sát rõ được độ biến động của chuỗi tỷ suất lợi nhuận giá cổ phiếu r stock_name với kỳ vọng khả năng xảy ra biến động lớn bất thường trong ngưỡng cho phép. Phương pháp này không làm thay đổi mức VaR mà ta muốn ước lượng và vẫn giữ nguyên được số lượng dữ liệu.
data_return_clean <- data_return %>% Return.clean(method = "boudt")
ggplot(aes(x = Index, y = Value), data = fortify(data_return_clean, melt = TRUE)) +
geom_line(color = "lightskyblue3") +
labs(title = paste0("Cleaned"," daily return of stock ", toupper(stock_name)),
subtitle=paste0("From ", as.yearmon(start_estimation_date),
" to ", as.yearmon(end_validation_date)),
y = NULL,
x = NULL,
caption = "Author: Nguyen Ngoc Binh")+
scale_x_date(date_labels = "%b %Y") +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
- Các giá trị đặc trưng của phân phối chuỗi lợi suất cổ phiếu
r stock_name sau khi làm sạch
data_return_clean %>%
table.Stats() %>%
rownames_to_column() %>%
slice(c(1, 3, 6, 7, 9, 13, 15, 16)) %>%
set_names(c("Indicators", "Cleaned Return")) %>%
datatable(rownames = FALSE, options = list(searching = FALSE))
Kiểm định giả thiết
Kiểm định mức độ tự tương quan
Thực hiện kiểm định về mức độ tự tương quan đối với bình phương sai số của chuỗi lợi suất, nhưng do trung bình của chuỗi lợi suất thường xấp xỉ bằng 0 nên để giảm tính phức tạp, chúng tôi thực hiện kiểm định mức độ tương quan đối với chuỗi bình phương giá trị lợi suất
b_acf <- function(dt){
p1 <- acf(dt, type = "correlation", lag.max = 30, plot = FALSE)
p2 <- acf(dt^2, type = "correlation", lag.max = 30, plot = FALSE)
p3 <- pacf(dt, lag.max = 30, plot = FALSE)
p4 <- pacf(dt^2, lag.max = 30, plot = FALSE)
series_name <- names(dt)
p <- data.frame(acf = p1$acf, lag = p1$lag, series_name = paste ("Acf -", series_name)) %>%
rbind(data.frame(acf = p2$acf, lag = p2$lag, series_name = paste ("Acf - Squared", series_name))) %>%
rbind(data.frame(acf = p3$acf, lag = p3$lag, series_name = paste ("Pacf -", series_name))) %>%
rbind(data.frame(acf = p4$acf, lag = p4$lag, series_name = paste ("Pacf - Squared", series_name)))
cutoff_upper <- 2/(length(dt))^0.5
cutoff_lower <- -2/(length(dt))^0.5
p %>% filter(lag != 0) %>%
ggplot(aes(x = lag, y = acf, color = series_name, group = series_name))+
geom_point(size = 1.5)+
geom_segment(aes(xend = lag, yend = 0), size = 1)+
geom_hline(yintercept = 0)+
geom_line(aes(y = cutoff_upper), color = "blue", linetype = 2)+
geom_line(aes(y = cutoff_lower), color = "blue", linetype = 2)+
facet_wrap(~ series_name, ncol = 2)+
expand_limits(y = c(-0.2, 0.2))+
scale_color_tq()+
theme_tq()+
labs(title = paste("ACF and PACF of", series_name),
subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)),
" to ", as.yearmon(end_validation_date)),
y = "Autocorrelation",
x = "Lags",
caption = "Author: Nguyen Ngoc Binh")+
theme(legend.position = "none"
)
}
b_acf(dt = data_return_clean)
Từ biểu đồ có thể thấy tồn tại hiện tượng tự tương quan trong chuỗi bình phương tỷ suất lợi nhuận của cổ phiếu r stock_name, điều này được thể hiện rõ hơn qua kiểm định Ljung-Box về mức độ tự tương quan.
Kiểm định Ljung-Box với cặp giả thiết:
- H0: Không tồn tại hiện tượng tự tương quan trong chuỗi bình phương lợi suất.
- H1: Tồn tại hiện tượng tự tương quan trong chuỗi bình phương lợi suất.
Kiểm định Ljung-Box là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Ljung-Box tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Ljung-Box với độ trễ lần lượt là 5, 10, 15 và độ tin cậy 95%, kết quả của kiểm định được thể hiện trong bảng dưới.
box_test <- lapply(c(5, 10, 15), function(i) {
z <- Box.test(data_return_clean^2, lag = i, type = "Ljung-Box")
data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>%
as.vector() %>%
set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>%
mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0",
TRUE ~ "Not reject H0")) %>%
return()
}) %>% bind_rows()
## Xuất bảng
box_test %>%
mutate_if(is.numeric, round, 4) %>%
datatable(options = list(searching = FALSE))
Cả ba kiểm định Ljung-Box đều cho kết quả tồn tại hiện tượng tự tương quan đối với chuỗi bình phương lợi suất, hay là tồn tại hiện tượng tự tương quan đối với độ biến động của chuỗi lợi suất. Điều này cho thấy độ biến động thời kỳ sau phụ thuộc vào độ biến động các thời kỳ trước đó (Conditional – time dependent).
Kiểm định phương sai của sai số thay đổi
Kiểm định phương sai của sai số thay đổi được thực hiện bằng kiểm định Engle’s ARCH.
Kiểm định Engle’s ARCH với cặp giả thiết:
- H0: Không tồn tại hiện tượng phương sai sai số thay đổi.
- H1: Tồn tại hiện tượng phương sai sai số thay đổi.
Kiểm định Engle’s ARCH là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Engle’s ARCH tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Engle’s ARCH với độ trễ lần lượt là 5, 10, 15 và độ tin cậy 95%, kết quả của kiểm định được thể hiện trong bảng
arch_test <- lapply(c(5, 10, 15), function(i) {
z <- FinTS::ArchTest(data_return_clean, lags = i)
data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>%
as.vector() %>%
set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>%
mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0",
TRUE ~ "Not reject H0")) %>%
return()
}) %>% bind_rows()
arch_test %>%
mutate_if(is.numeric, round, 4) %>%
datatable(options = list(searching = FALSE))
Cả 3 kiểm định Engle’s ARCH với độ trễ 5, 10 và 15 đều cho kết quả bác bỏ giả thiết không tồn tại hiện tượng phương sai sai số thay đổi. Điều này nghĩa là tồn tại hiện tượng phương sai sai số thay đổi (Conditional Heteroskedasticity) trong chuỗi tỷ suất lợi nhuận.
Ước lượng mô hình phù hợp
Phương pháp mô phỏng lịch sử
# ============================================================================
# Hàm mô phỏng lịch sử
# ============================================================================
f_hist_VaR <- function(i){
# Period
start.est.date <- start_estimation_date %m+% years(i)
end.est.date <- start.est.date %m+% years(2) - 1
start.val.date <- start.est.date %m+% years(2)
end.val.date <- start.est.date %m+% years(3) -1
period.est <- paste(format(start.est.date, format = "%d/%m/%Y"),
format(end.est.date, format = "%d/%m/%Y"), sep = " - ")
period.val <- paste(format(start.val.date, format = "%d/%m/%Y"),
format(end.val.date, format = "%d/%m/%Y"), sep = " - ")
# Set estimated data
estimation.return <- window(data_return_clean,
start = start.est.date,
end = end.est.date)
# Set validated data
validation.return <- window(data_return_clean,
start = start.val.date,
end = end.val.date)
# Length
est.length <- length(estimation.return)
val.length <- length(validation.return)
# VaR in percentage
VaR <- quantile(estimation.return, 0.05)
# No of times VaR > Validated period
n.over <- sum(validation.return < VaR)
# Critical value
critical.value <- val.length * 0.05 + qnorm(0.95) * sqrt(val.length * 0.05 * 0.95)
critical.value <- floor(critical.value)
data.frame(period.est, period.val, est.length, val.length, VaR, n.over, critical.value) %>%
mutate_if(is.factor, as.character) %>%
return()
}
hist_VaR <- map_dfr(0:3, f_hist_VaR)
# ---------------
# Export
# ---------------
hist_VaR %>%
mutate(VaR = VaR * 100) %>%
set_names(c("Estimated Period", "Validated Period", "No. of estimated observations",
"No. of validated observations", "VaR in percentage", "Test statistics",
"Critical value")) %>%
mutate_if(is.numeric, round, 2) %>%
datatable(options = list(searching = FALSE))
Phương pháp ước lượng độ biến động bằng mô hình GARCH
Ngoài phương pháp mô phỏng lịch sử, giá trị VaR95% có thể ước lượng được bằng việc áp dụng các mô hình GARCH . Để xác định được mô hình GARCH phù hợp cần lựa chọn được độ trễ (p,q), mô hình (GARCH, EGARCH, IGARCH, APARCH), phân phối của độ biến động chuẩn hóa.
Mẫu ước lượng
estimation_return <- data_return_clean %>%
window(start = start_estimation_date, end = end_estimation_date)
-
Khoảng thời gian r paste(fm_start_estimation_date, "-", fm_end_estimation_date)
-
Số quan sát r length(estimation_return)
Mẫu kiểm định
validation_return <- data_return_clean %>%
window(start = start_validation_date, end = end_validation_date)
-
Khoảng thời gian r paste(fm_start_estimation_date, "-", fm_end_estimation_date)
-
Số quan sát r length(validation_return)
Lựa chọn mô hình
Cách thức thực hiện như sau:
-
Bước 1: Ước lượng tất cả các mô hình với các cặp độ trễ, dạng mô hình, phân phối như bên dưới:
-
Độ trễ (p,q) = (1,1), (2,1), (1,2), (2,2)
-
Dạng mô hình: sGARCH, iGARCH, eGARCH, apARCH
-
Dạng phân phối: norm, snorm, std, sstd, ged, sged
-
Bước 2: Lựa chọn tất cả các mô hình có phân phối phù hợp:
-
Giá trị VaR ước lượng được rất khác nhau vì phân vị 5% của {et} khác nhau do các dạng phân phối khác nhau. Mục đích phần này là ước lượng được mức giảm giá tối đa với độ tin cậy cho trước (VaR) nên chúng tôi sẽ kiểm định mức độ phù hợp của mức VaR95% đối với từng phân phối.
-
Chúng tôi sử dụng phương pháp kiểm định nhị thức với độ tin cậy 95% để kiểm định mức độ phù hợp của mức VaR95%. Ở đây kiểm định sẽ được thực hiện với (mức VaR95%, 1 ngày) được xác định hàng ngày (mức VaR95% được xác định bằng độ lệch hàng ngày ước lượng được từ mô hình và mức phân vị 5% của phân phối {et}).
-
Chọn ra dạng mô hình, độ trễ, phân phối an toàn
-
Bước 3: Chọn mô hình cuối cùng có các tiêu chí AIC, BIC, SIC, HQIC (Akaike, Bayes, Shibata, Hannan-Quinn) nhỏ nhất:
-
Tính toán tất cả các chỉ số AIC, BIC, SIC, HQIC cho các mô hình
-
Tính hạng của từng mô hình theo các tiêu chí trên (tính dense_rank đối với từng chỉ số)
-
Tính tổng rank theo 4 tiêu chí
-
Chọn mô hình cuối cùng có tổng rank nhỏ nhất
f_backtest <- function(start,
end,
data = data_return_clean,
VaR.estimate = VaR_estimate,
alpha = 0.05,
conf.level = 0.99) {
start <- as.Date(start, origin = "1970-01-01")
end <- as.Date(end, origin = "1970-01-01")
data %>%
window(start = start, end = end) ->> validation.return
period.val <- paste(format(start, format = "%d/%m/%Y"),
format(end, format = "%d/%m/%Y"),
sep = " - ")
test.statistics <- sum(VaR.estimate > validation.return)
obs <- length(validation.return)
critical.value <- obs * alpha + qnorm(conf.level) * sqrt(obs * alpha * (1-alpha))
critical.value <- floor(critical.value)
conclusion <- case_when(test.statistics > critical.value ~ "Reject", TRUE ~ "Safety")
list(Period = period.val,
VaR.estimate = VaR.estimate,
obs = obs,
test.statistics = test.statistics,
critical.value = critical.value,
conclusion = conclusion) %>%
return()
}
ff <- function(p, q, models, distributions, data = estimation_return,
alpha = 0.05, conf.level = 0.99){
spec <- ugarchspec(variance.model = list(model = models, garchOrder = c(p, q)),
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
distribution.model = distributions)
modelfit <- ugarchfit(spec = spec, data = data)
#------------------------------------
# Kiểm định p-value các hệ số betas
#------------------------------------
betas <- modelfit@fit$robust.matcoef[,4]
check_betas <- case_when(all(betas[-1] <= 0.05) ~ "pass", TRUE ~ "reject")
#------------------------------------
# Kiểm định sự phù hợp của phân phối
#------------------------------------
VaR <- sigma(modelfit) * qdist(distribution = distributions,
p = 0.05,
mu = 0,
sigma = 1,
skew = coef(modelfit)["skew"],
shape = coef(modelfit)["shape"])
obs <- length(data)
critical.value <- obs * alpha + qnorm(conf.level) * sqrt(obs * alpha * (1-alpha))
critical.value <- floor(critical.value)
test.statistics <- sum(VaR > data)
check_distribution <- case_when(test.statistics > critical.value ~ "reject", TRUE ~ "pass")
#------------------------------------
# Hậu kiểm 1 năm và 6 tháng
#------------------------------------
VaR.estimate <- mean(tail(VaR, n = 30))
start.validation.date.6m <- start_validation_date
end.validation.date.6m <- start_validation_date %m+% months(6) -1
bt_1y <- f_backtest(start_validation_date,
end_validation_date,
data_return_clean,
VaR.estimate)$conclusion
bt_6m <- f_backtest(start.validation.date.6m,
end.validation.date.6m,
data_return_clean,
VaR.estimate)$conclusion
check_bt <- case_when(bt_1y == "Safety" & bt_6m == "Safety" ~ "pass", TRUE ~ "reject")
#------------------------------------
# Xuất kết quả
#------------------------------------
conclusion <- case_when(check_betas == "pass" &
check_distribution == "pass" &
check_bt == "pass" ~ "Safety",
TRUE ~ "Reject")
list(modelfit = modelfit,
VaR = VaR,
obs = obs,
critical.value = critical.value,
test.statistics = test.statistics,
conclusion = conclusion) %>%
return()
}
fit_function <- function(p, q, models, distributions, data = estimation_return,
alpha = 0.05, conf.level = 0.99){
out <- tryCatch(
{
message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions))
ff(p, q,
models,
distributions,
data = estimation_return,
alpha = 0.05,
conf.level = 0.99)
},
error = function(cond) {
message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions, "can't fit"))
message("Here's the original error message:")
message(cond)
return(NA)
},
warning = function(cond) {
message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions, "can't fit"))
message("Here's the original error message:")
message(cond)
return(NULL)
},
finally = {
message("Ok")
}
)
return(out)
}
## Run all models
garch_p <- 1:2
garch_q <- 1:2
garch_model <- c("sGARCH","iGARCH","eGARCH","apARCH")
garch_distribution <- c("norm","snorm","std","sstd","ged","sged")
garch_opt <- expand.grid(p = garch_p,
q = garch_q,
models = garch_model,
distributions = garch_distribution,
stringsAsFactors = FALSE)
garch_option <- garch_opt %>% as.list()
# fit all
all_results <- pmap(garch_option, fit_function)
# library(parallel)
# all_results2 <- parallel::mcmapply(FUN = fit_function, garch_opt$p, garch_opt$q, garch_opt$models, garch_opt$distributions)
Những mô hình có phân phối phù hợp
# Filter all results if conclusion = "Fail to Reject H0"
all_conlusions <- map(all_results, `[[`, "conclusion") == "Safety"
garch_opt %>%
cbind(all_conlusions) %>%
set_names(c("p", "q", "model name", "distribution", "is safety model")) %>%
datatable()
all_results <- all_results[all_conlusions]
garch_opt <- garch_opt %>% cbind(all_conlusions) %>% filter(all_conlusions == TRUE)
garch_lab <- garch_opt %>%
mutate(garch_labels = paste0(models, "(",p, ",", q, ")", " - ", distributions))
#
all_models <- all_results %>% map( `[[`, "modelfit")
# criteria
criteria <- all_models %>% map_dfc(likelihood) %>%
bind_rows(map_dfc(all_models, infocriteria)) %>%
t() %>%
as.data.frame() %>%
set_colnames(c("Log_Likelihood", "Akaike", "Bayes", "Shibata", "Hannan_Quinn")) %>%
bind_cols(garch_lab)
all_VaRs <- map(all_results, `[[`, "VaR")
- Mô hình có các tiêu chí Akaike, Bayes, Shibata, Hannan-Quinn tốt nhất
# Output
criteria_arrange <- criteria %>%
mutate(Akaike_rk = dense_rank(Akaike),
Bayes_rk = dense_rank(Bayes),
Shibata_rk = dense_rank(Shibata),
Hannan_Quinn_rk = dense_rank(Hannan_Quinn)) %>%
# Choose model by min dense_rank Criteria
mutate(min_rk = Akaike_rk + Bayes_rk + Shibata_rk + Hannan_Quinn_rk) %>%
arrange(min_rk) %>%
mutate_if(is.numeric, round, 4) %>%
mutate_at("Log_Likelihood", round, 1)
criteria_arrange %>%
select(garch_labels, Log_Likelihood, Akaike, Bayes, Shibata, Hannan_Quinn) %>%
formattable(
list(Log_Likelihood = color_tile("lightpink", "white"),
Akaike = color_tile("lightpink", "white"),
Bayes = color_tile("lightpink", "white"),
Shibata = color_tile("lightpink", "white"),
Hannan_Quinn = color_tile("lightpink", "white"))) %>%
as.datatable()
criteria_1 <- criteria_arrange %>% slice(1)
- Lựa chọn mô hình
r criteria_1$garch_labels để ước lượng và dự báo về độ biến động giá cổ phiếu r stock_name
garch_labels <- criteria_1$garch_labels
garch_model <- criteria_1$models
garch_order <- c(criteria_1$p, criteria_1$q)
garch_distribution <- criteria_1$distributions
garch_spec <- ugarchspec(variance.model = list(model = garch_model, garchOrder = garch_order),
mean.model = list(armaOrder=c(0,0)),
distribution.model = garch_distribution)
garch_fit <- ugarchfit(spec = garch_spec, data = estimation_return)
garch_VaR <- sigma(garch_fit) * qdist(distribution = garch_distribution,
p = 0.05,
mu = 0,
sigma = 1,
skew = coef(garch_fit)["skew"],
shape = coef(garch_fit)["shape"])
Kiểm định mức độ phù hợp của mô hình lựa chọn
estimation_return %>%
merge(garch_VaR) %>%
fortify(melt = FALSE) %>%
ggplot() +
geom_line(aes(x = Index, y = Return), color = "royalblue") +
geom_line(aes(x = Index, y = garch_VaR), color = "tomato") +
labs(title = paste("Return and Value at Risk of", stock_name),
subtitle=paste0("From ", as.yearmon(start_estimation_date),
" to ", as.yearmon(end_estimation_date)),
x = NULL,
y = NULL,
caption = "Author: Nguyen Ngoc Binh")
Sau khi ước lượng phương sai của sai số bằng mô hình r garch_labels, chúng tôi sẽ kiểm định mức độ tự tương quan và phương sai sai số thay đổi của chuỗi bình phương giá trị chuẩn hóa sai số
-
Kiểm định tự tương quan của sai số
-
Lược đồ ACF, PACF
# garch_residuals <- rugarch::residuals(garch_fit)
# names(garch_residuals) <- "Residuals"
# b_acf(garch_residuals)
# z <- residuals(garch_fit)/sigma(garch_fit) # = et
et <- rugarch::residuals(garch_fit, standardize = T)
names(et) <- "Standardized Residuals"
b_acf(et)
Từ biểu đồ trên có thể thấy đã không còn hiện tượng tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số (Squared Standardized Residuals)
- Kiểm định Ljung - Box
Sử dụng kiểm định Ljung-Box để kiểm định mức độ tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số. Kiểm định Ljung-Box với cặp giả thiết:
H0: Không tồn tại hiện tượng tự tương quan trong chuỗi bình phương giá trị chuẩn hóa sai số.
H1: Tồn tại hiện tượng tự tương quan trong chuỗi bình phương giá trị chuẩn hóa sai số.
Kiểm định Ljung-Box là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Ljung-Box tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Ljung-Box với độ trễ lần lượt là 5, 10, 15, kết quả của kiểm định được thể hiện trong bảng sau:
box_test_et <- lapply(c(5, 10, 15), function(i) {
z <- Box.test(et^2, lag = i, type = "Ljung-Box")
data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>%
as.vector() %>%
set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>%
mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0",
TRUE ~ "Not reject H0")) %>%
return()
}) %>% bind_rows()
## Xuất bảng
box_test_et %>%
mutate_if(is.numeric, round, 4) %>%
datatable()
Cả ba kiểm định Ljung-Box đều cho kết quả không tồn tại hiện tượng tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số
- Kiểm định phương sai sai số thay đổi
Tiến hành kiểm định Engle’s ARCH với độ tin cậy 95% cho chuỗi giá trị chuẩn hóa sai số, cặp giả thiết:
H0: Không tồn tại hiện tượng phương sai sai số thay đổi.
H1: Tồn tại hiện tượng phương sai sai số thay đổi.
Kiểm định Engle’s ARCH là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Engle’s ARCH tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Engle’s ARCH với độ trễ lần lượt là 5, 10, 15, kết quả của kiểm định được thể hiện trong bảng sau:
arch_test_et <- lapply(c(5,10,15),function(i) {
z <- FinTS::ArchTest(et, lags = i)
data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>%
as.vector() %>%
set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>%
mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0",
TRUE ~ "Not reject H0")) %>%
return()
}) %>% bind_rows()
arch_test_et %>%
mutate_if(is.numeric, round, 4) %>%
datatable()
Cả ba kiểm định Engle’s ARCH đều cho kết quả không tồn tại hiện tượng phương sai sai số thay đổi
- Kiểm định phân phối {et}
et %>%
fortify(melt = FALSE) %>%
set_names(c("Index", "resid")) %>%
ggplot(aes(sample = resid))+
stat_qq(alpha = 0.2, color = "#0055AA")+
geom_abline(color = "#C40003")+
labs(title = "QQ plot of Standardized residuals",
subtitle = paste(toupper(garch_distribution), "distribution -",
garch_model, "model"),
caption = "Author: Nguyen Ngoc Binh")
Một cách trực quan, có thể thấy phân phối của chuỗi sai số chuẩn hóa khá gần với phân phối r garch_distribution
- Kiểm định mức độ phù hợp của giá trị daily-VaR: Sử dụng kiểm định nhị thức với độ tin cậy 95%
Lớp các mô hình GARCH vẫn có khả năng dự báo nếu mức VaR xác định trong khoảng thời gian dự báo còn vượt qua được kiểm định mức độ phù hợp, trong đó mức VaR được xác định dựa trên các giá trị tham số của mô hình đã được xác định trong khoảng thời gian ước lượng. Trong thống kê, kiểm định như thế được gọi là kiểm định ngoài mẫu (out-of-sample). Chúng tôi thực hiện kiểm định ngoài mẫu đối với mô hình r garch_labels cũng cho kết quả phù hợp.
# Muc do phu hop cua Daily VaR - Kiem dinh VaR out-of-sample unclean
data_return %>%
window(start = start_validation_date, end = end_validation_date) ->> validation_return_unc
###################################### Choose 1 ###############################################
## Neu uoc luong tu du lieu khong lam sach #
# return_oos <- data_return #
#
## Neu uoc luong tu du lieu da lam sach #
estimation_return %>%
rbind(validation_return_unc) %>%
as.xts() ->> return_oos
################################################################################
test_spec <- garch_spec
test_fit <- ugarchfit(spec = test_spec,
data = return_oos,
out.sample = length(validation_return_unc))
test_fcst <- ugarchforecast(test_fit,
n.roll = length(validation_return_unc),
n.ahead = 1)
test_VaR <- sigma(test_fcst) * qdist(distribution = garch_distribution,
p = 0.05,
mu = 0,
sigma = 1,
skew = coef(test_fit)["skew"],
shape = coef(test_fit)["shape"])
test_VaR %>%
t() %>%
xts(order.by = ymd(rownames(.))) ->> test_VaR_backtest
no_obs <- length(validation_return_unc)
critical_value <- floor(no_obs*0.05 + qnorm(0.95) * sqrt(no_obs*0.95*0.05))
test_statistics <- sum(test_VaR_backtest > validation_return_unc)
VaR_out_of_sample <- c(no_obs, critical_value, test_statistics) %>%
set_names(c("No. of observations", "Critical value", "Test statistics"))
VaR_out_of_sample %>%
as.data.frame() %>%
datatable()
Mức độ phù hợp của giá trị daily - VaR95% được xác định từ mô hình r garch_labels còn có thể được theo dõi bằng biểu đồ sau:
validation_return_unc %>%
merge(test_VaR_backtest) %>%
fortify() %>%
ggplot() +
geom_line(aes(x = Index, y = Return), color = "dodgerblue")+
geom_line(aes(x = Index, y = `T.1`), color = "red") +
labs(x = NULL,
y = "Return",
title = paste(stock_name, "Price Volatility"),
subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)),
" to ", as.yearmon(end_validation_date)),
caption = "Author: Nguyen Ngoc Binh") +
theme(legend.title = element_blank(),
axis.text.x = element_text(angle = 30, hjust = 1))
Kiểm định ngoài mẫu đối với mô hình r garch_labels cho kết quả phù hợp.
- Kết luận:
Sau khi thực hiện kiểm định mức độ phù hợp của mô hình r garch_labels được lựa chọn, có thể thấy mô hình đã loại bỏ được hiện tượng tự tương quan đối với chuỗi bình phương sai số chuẩn hóa, loại bỏ được hiện tượng phương sai sai số thay đổi đối với chuỗi sai số chuẩn hóa và phân phối của chuỗi rất gần với phân phối r garch_distribution được lựa chọn. Thêm vào đó, tham số ước lượng được từ mô hình vượt qua được kiểm định ngoài mẫu về mức độ phù hợp của giá trị daily-VaR. Do đó, có thể kết luận mô hình r garch_labels phù hợp để ước lượng mức giảm giá tối đa (VaR95%).
Ước lượng và dự báo mức giảm giá tối đa từ mô hình r garch_labels
VaR tại những ngày khác nhau được ước tính khác nhau (do dữ liệu được cập nhật hằng ngày). Do đó, để giá trị (VaR95%, 1 ngày) có thể đại diện trong khoảng thời gian dài, chúng tôi sử dụng giá trị VaR95% trung bình trong một khoảng thời gian - n ngày quan sát cuối cùng.
Độ dài khoảng thời gian sử dụng là 30 ngày hoặc 60 ngày tùy thuộc vào độ dài dữ liệu và mức độ phù hợp của giá trị VaR trung bình trong khoảng thời gian đó .
VaR_estimate <- mean(tail(garch_VaR, n = 30))
Trong trường hợp với cổ phiếu r stock_name, giá trị VaR95% được xác định bằng trung bình của giá trị VaR95%,1 ngày ước lượng từ mô hình r garch_labels trong 30 ngày quan sát cuối cùng. Theo đó giá trị VaR95%,1 ngày được dùng để dự báo trong khoảng thời gian hậu kiểm xác định được là: r paste0(round(VaR_estimate * 100, 2),%)
Hậu kiểm giá trị VaR
################### Input ###########################################
start_validation_date_6m <- start_validation_date #
end_validation_date_6m <- start_validation_date %m+% months(6) -1 #
#####################################################################
garch_date_backtest <- list(start = c(start_validation_date, start_validation_date_6m),
end = c(end_validation_date, end_validation_date_6m))
garch_backtest <- pmap_dfr(garch_date_backtest, f_backtest)
garch_backtest %>%
mutate_at("VaR.estimate", function(x){paste0(round(x * 100, digits = 2), "%")}) %>%
datatable()
Hậu kiểm cho thấy ước lượng bằng mô hình r garch_labels cho ra kết quả an toàn
Nhìn chung cả hai phương pháp ước lượng bằng mô hình GARCH và phương pháp mô phỏng lịch sử đều cho kết quả hậu kiểm VaR an toàn. Tuy nhiên, có thể thấy chuỗi giá trị lợi suất tồn tại hiện tượng tự tương quan và phương sai sai số thay đổi, sau khi áp dụng mô hình r garch_labels đã triệt tiêu được hai hiện tượng này, do đó mô hình r garch_labels là mô hình phù hợp hơn so với phương pháp mô phỏng lịch sử để dự báo về mức giảm giá tối đa.
Ước lượng và dự báo mức giảm giá tối đa
Từ kết quả hậu kiểm ở trên, mô hình r garch_labels là mô hình phù hợp nhất để ước lượng và dự báo mức độ giảm giá tối đa với độ tin cậy 95% (VaR95%). Thực hiện lại ước lượng tham số cho mô hình r garch_labels trong khoảng thời gian từ r fm_start_estimation_date đến r fm_end_validation_date, sau đó lấy trung bình mức VaR95% trong 30 ngày quan sát gần nhất làm giá trị dự báo cho mức VaR95%.
Độ biến động trong 1 tháng, 3 tháng, 6 tháng tiếp theo được xác định từ độ biến động trong 1 ngày như sau:
- Độ biến động trong 1 tháng = Độ biến động trong 1 ngày x √21
- Độ biến động trong 3 tháng = Độ biến động trong 1 ngày x √63
- Độ biến động trong 6 tháng = Độ biến động trong 1 ngày x √126
# Re-estimated GARCH
spec <- ugarchspec(variance.model = list(model = garch_model, garchOrder = garch_order),
mean.model = list(armaOrder=c(0,0)),
distribution.model = garch_distribution)
fit <- ugarchfit(spec = spec, data = data_return_clean)
daily_VaR <- sigma(fit) * qdist(distribution = garch_distribution,
p = 0.05,
mu = 0,
sigma = 1,
skew = coef(fit)["skew"],
shape = coef(fit)["shape"])
daily_VaR %>%
tail(n = 30) %>%
mean() %>%
as.data.frame() %>%
set_names("VaR1day") %>%
mutate(VaR1month = VaR1day* sqrt(21),
VaR3month = VaR1day* sqrt(63),
VaR6month = VaR1day* sqrt(126)) %>%
set_names(c("VaR 1 day", "VaR 1 month next", "VaR 3 months next", "VaR 6 months next")) %>%
mutate_if(is.numeric, function(x){paste0(round(x * 100, digits = 2), "%")}) %>%
datatable(options = list(searching = FALSE))
Phương pháp ước tính độ biến động cho khoảng thời gian dài hơn 1 ngày từ độ biến động trong 1 ngày như trên theo công thức căn thời gian (square root of time formula) được áp dụng rộng rãi trong phân tích rủi ro thị trường (chẳng hạn, theo các hướng dẫn của Ủy ban Basel). Tuy nhiên, công thức được sử dụng như trên là công thức gần đúng và sai số sẽ càng lớn nếu khoảng thời gian ước tính càng dài. Đối với các kỳ dự báo dài (từ 1 tháng trở lên), phương pháp này sẽ đưa ra kết quả thận trọng hơn mức ước lượng chính xác, hay nói cách khác là đưa ra mức giảm giá tối đa tương ứng với từng mức độ tin cậy cao hơn mức giảm giá “thực” theo mô hình GARCH.
nguyenngocbinh/bml_stocks documentation built on March 20, 2022, 7:11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, knitr.table.format = "html", fig.align = "center") # Ngày 11/11/2019 Sửa lỗi betas # Ngày 12/1/2019: Sửa lỗi hậu kiểm
library(DT) library(formattable) library(fGarch) library(FinTS) # library(kableExtra) # library(TSA) for acf # library(tseries) # library(robustbase) # library(RcppRoll) # roll_mean # library(zoo) # library(TTR) library(readxl) #library(writexl) library(rugarch) library(PerformanceAnalytics) library(tidyverse) library(magrittr) library(lubridate) library(tidyquant)
rm(list = ls()) ######### Input ############# stock_name <- "DHG" # ############################# # Choose estimation period and validation period ############################## Input ############################# start_estimation_date <- as.Date("2012-12-01") # end_estimation_date <- start_estimation_date %m+% years(5) - 1 # start_validation_date <- start_estimation_date %m+% years(5) # end_validation_date <- start_estimation_date %m+% years(6) -1 # ################################################################## # format date to export fm_start_estimation_date <- format(start_estimation_date, "%d/%m/%Y") fm_end_estimation_date <- format(end_estimation_date, "%d/%m/%Y") fm_start_validation_date <- format(start_validation_date, "%d/%m/%Y") fm_end_validation_date <- format(end_validation_date, "%d/%m/%Y")
# http://www.unstarched.net/r-examples/rugarch/a-short-introduction-to-the-rugarch-package/
# legend.title = element_blank(), # plot.title = element_text(hjust = 0.5), # plot.subtitle = element_text(hjust = 0.5), theme_nnb <- theme_tq() + theme(plot.caption = element_text(hjust = 1), legend.background = element_blank(), panel.grid.minor.y = element_blank()) theme_set(theme_nnb)
Tóm lược
Tài liệu được soạn thảo với mục đích cung cấp thông tin chi tiết về kỹ thuật ước lượng và dự báo biến động giá cổ phiếu r stock_name. Tài liệu được cấu trúc theo 5 phần chính:
1. Diễn biến giá cổ phiếu r stock_name: phần này cung cấp cái nhìn tổng quát về diễn biến độ biến động giá cổ phiếu r stock_name.
2. Kiểm định giả thiết: kiểm định mức độ tự tương quan và phương sai sai số thay đổi của chuỗi lợi suất giá cổ phiếu r stock_name.
3. Ước lượng mô hình phù hợp: ước lượng mức giảm giá tối đa bằng phương pháp mô phỏng lịch sử và các mô hình GARCH, trong đó tập trung xác định mô hình GARCH phù hợp trong lớp các mô hình GARCH.
4. Hậu kiểm VaR: hậu kiểm giá trị VaR ước lượng được bằng phương pháp mô phỏng lịch sử và mô hình GARCH từ đó chọn ra phương pháp phù hợp.
5. Ước lượng và dự báo mức giảm giá tối đa: thực hiện ước lượng lại mức giảm giá tối đa bằng phương pháp phù hợp nhất đã lựa chọn được đối với khoảng dữ liệu mới nhất.
Các nội dung chính
-
Độ biến động và mức giảm giá tối đa của mặt hàng cổ phiếu
r stock_nameđược ước lượng dựa trên dữ liệu giá cổ phiếur stock_namehàng ngày – cophieu68.com trong khoảng thời gian từr fm_start_estimation_dateđếnr fm_end_validation_date. -
Dữ liệu được chia thành hai phần: khoảng thời gian từ
r fm_start_estimation_dateđếnr fm_end_estimation_dateđược dùng để ước lượng độ biến động và mức giảm giá tối đa VaR95%, khoảng thời gian từr fm_start_validation_dateđếnr fm_end_validation_dateđể hậu kiểm mức độ phù hợp của giá trị VaR95%. Sau khi xác định được phương pháp phù hợp để tính mức giảm giá tối đa, chúng tôi sẽ áp dụng phương pháp đó để ước lượng VaR95% trong khoảng thời gian từr fm_start_estimation_dateđếnr fm_end_validation_datetừ đó đưa ra dự báo về mức giảm giá tối đa của cổ phiếur stock_name.
my_folder <- paste0("D:/stock/co_phieu/") price_raw <- read_csv(paste0(my_folder,"/excel_",stock_name,".csv")) %>% mutate(time_stamp = ymd(`<DTYYYYMMDD>`), price = `<CloseFixed>`) %>% na.omit() # Change data type data <- xts(x = price_raw$price, order.by = price_raw$time_stamp)
Diễn biến thực tế của giá cổ phiếu r stock_name
price_raw_filter <- price_raw %>% filter(time_stamp >= end_validation_date - years(2)) price_raw_filter %>% mutate(ma30 = rollmean(price, 30, align = "center", fill = NA), ma7 = rollmean(price, 7, align = "center", fill = NA)) %>% ggplot() + geom_point(aes(x = time_stamp, y = price, color = "daily"), size = 1.5, alpha = 0.5) + geom_line(aes(x = time_stamp, y = ma30, color = "ma30"), size = 1) + geom_line(aes(x = time_stamp, y = ma7, color = "ma7"), size = 1)+ scale_color_tq()+ #scale_color_manual(values = c("#2c3e50", "#e31a1c", "#18BC9C"))+ scale_x_date(date_labels = "%b %Y") + labs(x = NULL, y = "Price", title = paste("Market price of stock",stock_name), subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)), " to ", as.yearmon(end_validation_date)), caption = "Author: Nguyen Ngoc Binh") + theme(legend.position = c(0.75, 0.95), legend.direction = "horizontal", legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1))
Biểu đồ cho thấy giá cổ phiếu r stock_name nằm trong khoảng từ r min(price_raw_filter$price) đến r max(price_raw_filter$price) nghìn đồng. Để thấy rõ độ biến động của giá cổ phiếu r stock_name, chúng tôi sẽ quan sát biểu đồ chuỗi lợi suất (Return).
Tỷ suất lợi nhuận cổ phiếu r stock_name
# Calculate return data_return <- Return.calculate(data, method="discrete") %>% na.omit() %>% set_names("Return") %>% window(start = start_estimation_date, end = end_validation_date) ggplot(aes(x = Index, y = Value), data = fortify(data_return, melt = TRUE)) + geom_line(color = "skyblue3") + labs(title = paste0("Daily return of stock ", stock_name %>% toupper()), subtitle=paste0("From ", as.yearmon(start_estimation_date), " to ", as.yearmon(end_validation_date)), y = NULL, x = NULL, caption = "Author: Nguyen Ngoc Binh")+ scale_x_date(date_labels = "%b %Y") + theme(axis.text.x = element_text(angle = 30, hjust = 1))
- Độ biến động của giá cổ phiếu
r stock_nameđược thể hiện rõ hơn qua các giá trị đặc trưng của phân phối chuỗi lợi suất
data_return %>% table.Stats() %>% rownames_to_column() %>% slice(c(1, 3, 6, 7, 9, 13, 15, 16)) %>% set_names(c("Indicators", "Return")) %>% datatable(rownames = FALSE, options = list(searching = FALSE))
- Làm sạch dữ liệu
Sử dụng phương pháp làm sạch dữ liệu (Clean data) bằng phương pháp “boudt” để quan sát rõ được độ biến động của chuỗi tỷ suất lợi nhuận giá cổ phiếu r stock_name với kỳ vọng khả năng xảy ra biến động lớn bất thường trong ngưỡng cho phép. Phương pháp này không làm thay đổi mức VaR mà ta muốn ước lượng và vẫn giữ nguyên được số lượng dữ liệu.
data_return_clean <- data_return %>% Return.clean(method = "boudt") ggplot(aes(x = Index, y = Value), data = fortify(data_return_clean, melt = TRUE)) + geom_line(color = "lightskyblue3") + labs(title = paste0("Cleaned"," daily return of stock ", toupper(stock_name)), subtitle=paste0("From ", as.yearmon(start_estimation_date), " to ", as.yearmon(end_validation_date)), y = NULL, x = NULL, caption = "Author: Nguyen Ngoc Binh")+ scale_x_date(date_labels = "%b %Y") + theme(axis.text.x = element_text(angle = 30, hjust = 1))
- Các giá trị đặc trưng của phân phối chuỗi lợi suất cổ phiếu
r stock_namesau khi làm sạch
data_return_clean %>% table.Stats() %>% rownames_to_column() %>% slice(c(1, 3, 6, 7, 9, 13, 15, 16)) %>% set_names(c("Indicators", "Cleaned Return")) %>% datatable(rownames = FALSE, options = list(searching = FALSE))
Kiểm định giả thiết
Kiểm định mức độ tự tương quan
Thực hiện kiểm định về mức độ tự tương quan đối với bình phương sai số của chuỗi lợi suất, nhưng do trung bình của chuỗi lợi suất thường xấp xỉ bằng 0 nên để giảm tính phức tạp, chúng tôi thực hiện kiểm định mức độ tương quan đối với chuỗi bình phương giá trị lợi suất
b_acf <- function(dt){ p1 <- acf(dt, type = "correlation", lag.max = 30, plot = FALSE) p2 <- acf(dt^2, type = "correlation", lag.max = 30, plot = FALSE) p3 <- pacf(dt, lag.max = 30, plot = FALSE) p4 <- pacf(dt^2, lag.max = 30, plot = FALSE) series_name <- names(dt) p <- data.frame(acf = p1$acf, lag = p1$lag, series_name = paste ("Acf -", series_name)) %>% rbind(data.frame(acf = p2$acf, lag = p2$lag, series_name = paste ("Acf - Squared", series_name))) %>% rbind(data.frame(acf = p3$acf, lag = p3$lag, series_name = paste ("Pacf -", series_name))) %>% rbind(data.frame(acf = p4$acf, lag = p4$lag, series_name = paste ("Pacf - Squared", series_name))) cutoff_upper <- 2/(length(dt))^0.5 cutoff_lower <- -2/(length(dt))^0.5 p %>% filter(lag != 0) %>% ggplot(aes(x = lag, y = acf, color = series_name, group = series_name))+ geom_point(size = 1.5)+ geom_segment(aes(xend = lag, yend = 0), size = 1)+ geom_hline(yintercept = 0)+ geom_line(aes(y = cutoff_upper), color = "blue", linetype = 2)+ geom_line(aes(y = cutoff_lower), color = "blue", linetype = 2)+ facet_wrap(~ series_name, ncol = 2)+ expand_limits(y = c(-0.2, 0.2))+ scale_color_tq()+ theme_tq()+ labs(title = paste("ACF and PACF of", series_name), subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)), " to ", as.yearmon(end_validation_date)), y = "Autocorrelation", x = "Lags", caption = "Author: Nguyen Ngoc Binh")+ theme(legend.position = "none" ) } b_acf(dt = data_return_clean)
Từ biểu đồ có thể thấy tồn tại hiện tượng tự tương quan trong chuỗi bình phương tỷ suất lợi nhuận của cổ phiếu r stock_name, điều này được thể hiện rõ hơn qua kiểm định Ljung-Box về mức độ tự tương quan.
Kiểm định Ljung-Box với cặp giả thiết:
- H0: Không tồn tại hiện tượng tự tương quan trong chuỗi bình phương lợi suất.
- H1: Tồn tại hiện tượng tự tương quan trong chuỗi bình phương lợi suất.
Kiểm định Ljung-Box là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Ljung-Box tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Ljung-Box với độ trễ lần lượt là 5, 10, 15 và độ tin cậy 95%, kết quả của kiểm định được thể hiện trong bảng dưới.
box_test <- lapply(c(5, 10, 15), function(i) { z <- Box.test(data_return_clean^2, lag = i, type = "Ljung-Box") data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>% as.vector() %>% set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>% mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0", TRUE ~ "Not reject H0")) %>% return() }) %>% bind_rows() ## Xuất bảng box_test %>% mutate_if(is.numeric, round, 4) %>% datatable(options = list(searching = FALSE))
Cả ba kiểm định Ljung-Box đều cho kết quả tồn tại hiện tượng tự tương quan đối với chuỗi bình phương lợi suất, hay là tồn tại hiện tượng tự tương quan đối với độ biến động của chuỗi lợi suất. Điều này cho thấy độ biến động thời kỳ sau phụ thuộc vào độ biến động các thời kỳ trước đó (Conditional – time dependent).
Kiểm định phương sai của sai số thay đổi
Kiểm định phương sai của sai số thay đổi được thực hiện bằng kiểm định Engle’s ARCH. Kiểm định Engle’s ARCH với cặp giả thiết:
- H0: Không tồn tại hiện tượng phương sai sai số thay đổi.
- H1: Tồn tại hiện tượng phương sai sai số thay đổi.
Kiểm định Engle’s ARCH là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Engle’s ARCH tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Engle’s ARCH với độ trễ lần lượt là 5, 10, 15 và độ tin cậy 95%, kết quả của kiểm định được thể hiện trong bảng
arch_test <- lapply(c(5, 10, 15), function(i) { z <- FinTS::ArchTest(data_return_clean, lags = i) data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>% as.vector() %>% set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>% mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0", TRUE ~ "Not reject H0")) %>% return() }) %>% bind_rows() arch_test %>% mutate_if(is.numeric, round, 4) %>% datatable(options = list(searching = FALSE))
Cả 3 kiểm định Engle’s ARCH với độ trễ 5, 10 và 15 đều cho kết quả bác bỏ giả thiết không tồn tại hiện tượng phương sai sai số thay đổi. Điều này nghĩa là tồn tại hiện tượng phương sai sai số thay đổi (Conditional Heteroskedasticity) trong chuỗi tỷ suất lợi nhuận.
Ước lượng mô hình phù hợp
Phương pháp mô phỏng lịch sử
# ============================================================================ # Hàm mô phỏng lịch sử # ============================================================================ f_hist_VaR <- function(i){ # Period start.est.date <- start_estimation_date %m+% years(i) end.est.date <- start.est.date %m+% years(2) - 1 start.val.date <- start.est.date %m+% years(2) end.val.date <- start.est.date %m+% years(3) -1 period.est <- paste(format(start.est.date, format = "%d/%m/%Y"), format(end.est.date, format = "%d/%m/%Y"), sep = " - ") period.val <- paste(format(start.val.date, format = "%d/%m/%Y"), format(end.val.date, format = "%d/%m/%Y"), sep = " - ") # Set estimated data estimation.return <- window(data_return_clean, start = start.est.date, end = end.est.date) # Set validated data validation.return <- window(data_return_clean, start = start.val.date, end = end.val.date) # Length est.length <- length(estimation.return) val.length <- length(validation.return) # VaR in percentage VaR <- quantile(estimation.return, 0.05) # No of times VaR > Validated period n.over <- sum(validation.return < VaR) # Critical value critical.value <- val.length * 0.05 + qnorm(0.95) * sqrt(val.length * 0.05 * 0.95) critical.value <- floor(critical.value) data.frame(period.est, period.val, est.length, val.length, VaR, n.over, critical.value) %>% mutate_if(is.factor, as.character) %>% return() } hist_VaR <- map_dfr(0:3, f_hist_VaR) # --------------- # Export # --------------- hist_VaR %>% mutate(VaR = VaR * 100) %>% set_names(c("Estimated Period", "Validated Period", "No. of estimated observations", "No. of validated observations", "VaR in percentage", "Test statistics", "Critical value")) %>% mutate_if(is.numeric, round, 2) %>% datatable(options = list(searching = FALSE))
Phương pháp ước lượng độ biến động bằng mô hình GARCH
Ngoài phương pháp mô phỏng lịch sử, giá trị VaR95% có thể ước lượng được bằng việc áp dụng các mô hình GARCH . Để xác định được mô hình GARCH phù hợp cần lựa chọn được độ trễ (p,q), mô hình (GARCH, EGARCH, IGARCH, APARCH), phân phối của độ biến động chuẩn hóa.
Mẫu ước lượng
estimation_return <- data_return_clean %>% window(start = start_estimation_date, end = end_estimation_date)
-
Khoảng thời gian
r paste(fm_start_estimation_date, "-", fm_end_estimation_date) -
Số quan sát
r length(estimation_return)
Mẫu kiểm định
validation_return <- data_return_clean %>% window(start = start_validation_date, end = end_validation_date)
-
Khoảng thời gian
r paste(fm_start_estimation_date, "-", fm_end_estimation_date) -
Số quan sát
r length(validation_return)
Lựa chọn mô hình
Cách thức thực hiện như sau:
-
Bước 1: Ước lượng tất cả các mô hình với các cặp độ trễ, dạng mô hình, phân phối như bên dưới:
-
Độ trễ (p,q) = (1,1), (2,1), (1,2), (2,2)
-
Dạng mô hình: sGARCH, iGARCH, eGARCH, apARCH
-
Dạng phân phối: norm, snorm, std, sstd, ged, sged
-
Bước 2: Lựa chọn tất cả các mô hình có phân phối phù hợp:
-
Giá trị VaR ước lượng được rất khác nhau vì phân vị 5% của {et} khác nhau do các dạng phân phối khác nhau. Mục đích phần này là ước lượng được mức giảm giá tối đa với độ tin cậy cho trước (VaR) nên chúng tôi sẽ kiểm định mức độ phù hợp của mức VaR95% đối với từng phân phối.
-
Chúng tôi sử dụng phương pháp kiểm định nhị thức với độ tin cậy 95% để kiểm định mức độ phù hợp của mức VaR95%. Ở đây kiểm định sẽ được thực hiện với (mức VaR95%, 1 ngày) được xác định hàng ngày (mức VaR95% được xác định bằng độ lệch hàng ngày ước lượng được từ mô hình và mức phân vị 5% của phân phối {et}).
-
Chọn ra dạng mô hình, độ trễ, phân phối an toàn
-
Bước 3: Chọn mô hình cuối cùng có các tiêu chí AIC, BIC, SIC, HQIC (Akaike, Bayes, Shibata, Hannan-Quinn) nhỏ nhất:
-
Tính toán tất cả các chỉ số AIC, BIC, SIC, HQIC cho các mô hình
-
Tính hạng của từng mô hình theo các tiêu chí trên (tính dense_rank đối với từng chỉ số)
-
Tính tổng rank theo 4 tiêu chí
-
Chọn mô hình cuối cùng có tổng rank nhỏ nhất
f_backtest <- function(start, end, data = data_return_clean, VaR.estimate = VaR_estimate, alpha = 0.05, conf.level = 0.99) { start <- as.Date(start, origin = "1970-01-01") end <- as.Date(end, origin = "1970-01-01") data %>% window(start = start, end = end) ->> validation.return period.val <- paste(format(start, format = "%d/%m/%Y"), format(end, format = "%d/%m/%Y"), sep = " - ") test.statistics <- sum(VaR.estimate > validation.return) obs <- length(validation.return) critical.value <- obs * alpha + qnorm(conf.level) * sqrt(obs * alpha * (1-alpha)) critical.value <- floor(critical.value) conclusion <- case_when(test.statistics > critical.value ~ "Reject", TRUE ~ "Safety") list(Period = period.val, VaR.estimate = VaR.estimate, obs = obs, test.statistics = test.statistics, critical.value = critical.value, conclusion = conclusion) %>% return() }
ff <- function(p, q, models, distributions, data = estimation_return, alpha = 0.05, conf.level = 0.99){ spec <- ugarchspec(variance.model = list(model = models, garchOrder = c(p, q)), mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), distribution.model = distributions) modelfit <- ugarchfit(spec = spec, data = data) #------------------------------------ # Kiểm định p-value các hệ số betas #------------------------------------ betas <- modelfit@fit$robust.matcoef[,4] check_betas <- case_when(all(betas[-1] <= 0.05) ~ "pass", TRUE ~ "reject") #------------------------------------ # Kiểm định sự phù hợp của phân phối #------------------------------------ VaR <- sigma(modelfit) * qdist(distribution = distributions, p = 0.05, mu = 0, sigma = 1, skew = coef(modelfit)["skew"], shape = coef(modelfit)["shape"]) obs <- length(data) critical.value <- obs * alpha + qnorm(conf.level) * sqrt(obs * alpha * (1-alpha)) critical.value <- floor(critical.value) test.statistics <- sum(VaR > data) check_distribution <- case_when(test.statistics > critical.value ~ "reject", TRUE ~ "pass") #------------------------------------ # Hậu kiểm 1 năm và 6 tháng #------------------------------------ VaR.estimate <- mean(tail(VaR, n = 30)) start.validation.date.6m <- start_validation_date end.validation.date.6m <- start_validation_date %m+% months(6) -1 bt_1y <- f_backtest(start_validation_date, end_validation_date, data_return_clean, VaR.estimate)$conclusion bt_6m <- f_backtest(start.validation.date.6m, end.validation.date.6m, data_return_clean, VaR.estimate)$conclusion check_bt <- case_when(bt_1y == "Safety" & bt_6m == "Safety" ~ "pass", TRUE ~ "reject") #------------------------------------ # Xuất kết quả #------------------------------------ conclusion <- case_when(check_betas == "pass" & check_distribution == "pass" & check_bt == "pass" ~ "Safety", TRUE ~ "Reject") list(modelfit = modelfit, VaR = VaR, obs = obs, critical.value = critical.value, test.statistics = test.statistics, conclusion = conclusion) %>% return() }
fit_function <- function(p, q, models, distributions, data = estimation_return, alpha = 0.05, conf.level = 0.99){ out <- tryCatch( { message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions)) ff(p, q, models, distributions, data = estimation_return, alpha = 0.05, conf.level = 0.99) }, error = function(cond) { message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions, "can't fit")) message("Here's the original error message:") message(cond) return(NA) }, warning = function(cond) { message(paste("model:", models, "p:", p, "q:", q, "distribution:", distributions, "can't fit")) message("Here's the original error message:") message(cond) return(NULL) }, finally = { message("Ok") } ) return(out) }
## Run all models garch_p <- 1:2 garch_q <- 1:2 garch_model <- c("sGARCH","iGARCH","eGARCH","apARCH") garch_distribution <- c("norm","snorm","std","sstd","ged","sged") garch_opt <- expand.grid(p = garch_p, q = garch_q, models = garch_model, distributions = garch_distribution, stringsAsFactors = FALSE) garch_option <- garch_opt %>% as.list() # fit all all_results <- pmap(garch_option, fit_function) # library(parallel) # all_results2 <- parallel::mcmapply(FUN = fit_function, garch_opt$p, garch_opt$q, garch_opt$models, garch_opt$distributions)
Những mô hình có phân phối phù hợp
# Filter all results if conclusion = "Fail to Reject H0" all_conlusions <- map(all_results, `[[`, "conclusion") == "Safety" garch_opt %>% cbind(all_conlusions) %>% set_names(c("p", "q", "model name", "distribution", "is safety model")) %>% datatable()
all_results <- all_results[all_conlusions] garch_opt <- garch_opt %>% cbind(all_conlusions) %>% filter(all_conlusions == TRUE) garch_lab <- garch_opt %>% mutate(garch_labels = paste0(models, "(",p, ",", q, ")", " - ", distributions))
# all_models <- all_results %>% map( `[[`, "modelfit") # criteria criteria <- all_models %>% map_dfc(likelihood) %>% bind_rows(map_dfc(all_models, infocriteria)) %>% t() %>% as.data.frame() %>% set_colnames(c("Log_Likelihood", "Akaike", "Bayes", "Shibata", "Hannan_Quinn")) %>% bind_cols(garch_lab)
all_VaRs <- map(all_results, `[[`, "VaR")
- Mô hình có các tiêu chí Akaike, Bayes, Shibata, Hannan-Quinn tốt nhất
# Output criteria_arrange <- criteria %>% mutate(Akaike_rk = dense_rank(Akaike), Bayes_rk = dense_rank(Bayes), Shibata_rk = dense_rank(Shibata), Hannan_Quinn_rk = dense_rank(Hannan_Quinn)) %>% # Choose model by min dense_rank Criteria mutate(min_rk = Akaike_rk + Bayes_rk + Shibata_rk + Hannan_Quinn_rk) %>% arrange(min_rk) %>% mutate_if(is.numeric, round, 4) %>% mutate_at("Log_Likelihood", round, 1) criteria_arrange %>% select(garch_labels, Log_Likelihood, Akaike, Bayes, Shibata, Hannan_Quinn) %>% formattable( list(Log_Likelihood = color_tile("lightpink", "white"), Akaike = color_tile("lightpink", "white"), Bayes = color_tile("lightpink", "white"), Shibata = color_tile("lightpink", "white"), Hannan_Quinn = color_tile("lightpink", "white"))) %>% as.datatable() criteria_1 <- criteria_arrange %>% slice(1)
- Lựa chọn mô hình
r criteria_1$garch_labelsđể ước lượng và dự báo về độ biến động giá cổ phiếur stock_name
garch_labels <- criteria_1$garch_labels garch_model <- criteria_1$models garch_order <- c(criteria_1$p, criteria_1$q) garch_distribution <- criteria_1$distributions garch_spec <- ugarchspec(variance.model = list(model = garch_model, garchOrder = garch_order), mean.model = list(armaOrder=c(0,0)), distribution.model = garch_distribution) garch_fit <- ugarchfit(spec = garch_spec, data = estimation_return) garch_VaR <- sigma(garch_fit) * qdist(distribution = garch_distribution, p = 0.05, mu = 0, sigma = 1, skew = coef(garch_fit)["skew"], shape = coef(garch_fit)["shape"])
Kiểm định mức độ phù hợp của mô hình lựa chọn
estimation_return %>% merge(garch_VaR) %>% fortify(melt = FALSE) %>% ggplot() + geom_line(aes(x = Index, y = Return), color = "royalblue") + geom_line(aes(x = Index, y = garch_VaR), color = "tomato") + labs(title = paste("Return and Value at Risk of", stock_name), subtitle=paste0("From ", as.yearmon(start_estimation_date), " to ", as.yearmon(end_estimation_date)), x = NULL, y = NULL, caption = "Author: Nguyen Ngoc Binh")
Sau khi ước lượng phương sai của sai số bằng mô hình r garch_labels, chúng tôi sẽ kiểm định mức độ tự tương quan và phương sai sai số thay đổi của chuỗi bình phương giá trị chuẩn hóa sai số
-
Kiểm định tự tương quan của sai số
-
Lược đồ ACF, PACF
# garch_residuals <- rugarch::residuals(garch_fit) # names(garch_residuals) <- "Residuals" # b_acf(garch_residuals) # z <- residuals(garch_fit)/sigma(garch_fit) # = et et <- rugarch::residuals(garch_fit, standardize = T) names(et) <- "Standardized Residuals" b_acf(et)
Từ biểu đồ trên có thể thấy đã không còn hiện tượng tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số (Squared Standardized Residuals)
- Kiểm định Ljung - Box
Sử dụng kiểm định Ljung-Box để kiểm định mức độ tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số. Kiểm định Ljung-Box với cặp giả thiết:
H0: Không tồn tại hiện tượng tự tương quan trong chuỗi bình phương giá trị chuẩn hóa sai số.
H1: Tồn tại hiện tượng tự tương quan trong chuỗi bình phương giá trị chuẩn hóa sai số.
Kiểm định Ljung-Box là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Ljung-Box tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Ljung-Box với độ trễ lần lượt là 5, 10, 15, kết quả của kiểm định được thể hiện trong bảng sau:
box_test_et <- lapply(c(5, 10, 15), function(i) { z <- Box.test(et^2, lag = i, type = "Ljung-Box") data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>% as.vector() %>% set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>% mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0", TRUE ~ "Not reject H0")) %>% return() }) %>% bind_rows() ## Xuất bảng box_test_et %>% mutate_if(is.numeric, round, 4) %>% datatable()
Cả ba kiểm định Ljung-Box đều cho kết quả không tồn tại hiện tượng tự tương quan đối với chuỗi bình phương giá trị chuẩn hóa sai số
- Kiểm định phương sai sai số thay đổi
Tiến hành kiểm định Engle’s ARCH với độ tin cậy 95% cho chuỗi giá trị chuẩn hóa sai số, cặp giả thiết:
H0: Không tồn tại hiện tượng phương sai sai số thay đổi.
H1: Tồn tại hiện tượng phương sai sai số thay đổi.
Kiểm định Engle’s ARCH là kiểm định một phía bên phải (nếu giá trị kiểm định – Test Statistics lớn hơn ngưỡng cho phép – Critical Value thì bác bỏ giả thiết H0). Giá trị kiểm định của kiểm định Engle’s ARCH tuân theo phân phối chi-squared với bậc tự do bằng độ trễ (lags) của kiểm định. Thực hiện kiểm định Engle’s ARCH với độ trễ lần lượt là 5, 10, 15, kết quả của kiểm định được thể hiện trong bảng sau:
arch_test_et <- lapply(c(5,10,15),function(i) { z <- FinTS::ArchTest(et, lags = i) data.frame(z$parameter, z$statistic, qchisq(.95, df = i), z$p.value) %>% as.vector() %>% set_names(c("Lag", "Test Statistic", "Critical Value", "p_value")) %>% mutate(Conclusion = case_when(p_value <= 0.05 ~ "Reject H0", TRUE ~ "Not reject H0")) %>% return() }) %>% bind_rows() arch_test_et %>% mutate_if(is.numeric, round, 4) %>% datatable()
Cả ba kiểm định Engle’s ARCH đều cho kết quả không tồn tại hiện tượng phương sai sai số thay đổi
- Kiểm định phân phối {et}
et %>% fortify(melt = FALSE) %>% set_names(c("Index", "resid")) %>% ggplot(aes(sample = resid))+ stat_qq(alpha = 0.2, color = "#0055AA")+ geom_abline(color = "#C40003")+ labs(title = "QQ plot of Standardized residuals", subtitle = paste(toupper(garch_distribution), "distribution -", garch_model, "model"), caption = "Author: Nguyen Ngoc Binh")
Một cách trực quan, có thể thấy phân phối của chuỗi sai số chuẩn hóa khá gần với phân phối r garch_distribution
- Kiểm định mức độ phù hợp của giá trị daily-VaR: Sử dụng kiểm định nhị thức với độ tin cậy 95%
Lớp các mô hình GARCH vẫn có khả năng dự báo nếu mức VaR xác định trong khoảng thời gian dự báo còn vượt qua được kiểm định mức độ phù hợp, trong đó mức VaR được xác định dựa trên các giá trị tham số của mô hình đã được xác định trong khoảng thời gian ước lượng. Trong thống kê, kiểm định như thế được gọi là kiểm định ngoài mẫu (out-of-sample). Chúng tôi thực hiện kiểm định ngoài mẫu đối với mô hình r garch_labels cũng cho kết quả phù hợp.
# Muc do phu hop cua Daily VaR - Kiem dinh VaR out-of-sample unclean data_return %>% window(start = start_validation_date, end = end_validation_date) ->> validation_return_unc ###################################### Choose 1 ############################################### ## Neu uoc luong tu du lieu khong lam sach # # return_oos <- data_return # # ## Neu uoc luong tu du lieu da lam sach # estimation_return %>% rbind(validation_return_unc) %>% as.xts() ->> return_oos ################################################################################ test_spec <- garch_spec test_fit <- ugarchfit(spec = test_spec, data = return_oos, out.sample = length(validation_return_unc)) test_fcst <- ugarchforecast(test_fit, n.roll = length(validation_return_unc), n.ahead = 1) test_VaR <- sigma(test_fcst) * qdist(distribution = garch_distribution, p = 0.05, mu = 0, sigma = 1, skew = coef(test_fit)["skew"], shape = coef(test_fit)["shape"]) test_VaR %>% t() %>% xts(order.by = ymd(rownames(.))) ->> test_VaR_backtest no_obs <- length(validation_return_unc) critical_value <- floor(no_obs*0.05 + qnorm(0.95) * sqrt(no_obs*0.95*0.05)) test_statistics <- sum(test_VaR_backtest > validation_return_unc) VaR_out_of_sample <- c(no_obs, critical_value, test_statistics) %>% set_names(c("No. of observations", "Critical value", "Test statistics")) VaR_out_of_sample %>% as.data.frame() %>% datatable()
Mức độ phù hợp của giá trị daily - VaR95% được xác định từ mô hình r garch_labels còn có thể được theo dõi bằng biểu đồ sau:
validation_return_unc %>% merge(test_VaR_backtest) %>% fortify() %>% ggplot() + geom_line(aes(x = Index, y = Return), color = "dodgerblue")+ geom_line(aes(x = Index, y = `T.1`), color = "red") + labs(x = NULL, y = "Return", title = paste(stock_name, "Price Volatility"), subtitle = paste0("From ", as.yearmon(end_validation_date - years(2)), " to ", as.yearmon(end_validation_date)), caption = "Author: Nguyen Ngoc Binh") + theme(legend.title = element_blank(), axis.text.x = element_text(angle = 30, hjust = 1))
Kiểm định ngoài mẫu đối với mô hình r garch_labels cho kết quả phù hợp.
- Kết luận:
Sau khi thực hiện kiểm định mức độ phù hợp của mô hình r garch_labels được lựa chọn, có thể thấy mô hình đã loại bỏ được hiện tượng tự tương quan đối với chuỗi bình phương sai số chuẩn hóa, loại bỏ được hiện tượng phương sai sai số thay đổi đối với chuỗi sai số chuẩn hóa và phân phối của chuỗi rất gần với phân phối r garch_distribution được lựa chọn. Thêm vào đó, tham số ước lượng được từ mô hình vượt qua được kiểm định ngoài mẫu về mức độ phù hợp của giá trị daily-VaR. Do đó, có thể kết luận mô hình r garch_labels phù hợp để ước lượng mức giảm giá tối đa (VaR95%).
Ước lượng và dự báo mức giảm giá tối đa từ mô hình r garch_labels
VaR tại những ngày khác nhau được ước tính khác nhau (do dữ liệu được cập nhật hằng ngày). Do đó, để giá trị (VaR95%, 1 ngày) có thể đại diện trong khoảng thời gian dài, chúng tôi sử dụng giá trị VaR95% trung bình trong một khoảng thời gian - n ngày quan sát cuối cùng.
Độ dài khoảng thời gian sử dụng là 30 ngày hoặc 60 ngày tùy thuộc vào độ dài dữ liệu và mức độ phù hợp của giá trị VaR trung bình trong khoảng thời gian đó .
VaR_estimate <- mean(tail(garch_VaR, n = 30))
Trong trường hợp với cổ phiếu r stock_name, giá trị VaR95% được xác định bằng trung bình của giá trị VaR95%,1 ngày ước lượng từ mô hình r garch_labels trong 30 ngày quan sát cuối cùng. Theo đó giá trị VaR95%,1 ngày được dùng để dự báo trong khoảng thời gian hậu kiểm xác định được là: r paste0(round(VaR_estimate * 100, 2),%)
Hậu kiểm giá trị VaR
################### Input ########################################### start_validation_date_6m <- start_validation_date # end_validation_date_6m <- start_validation_date %m+% months(6) -1 # ##################################################################### garch_date_backtest <- list(start = c(start_validation_date, start_validation_date_6m), end = c(end_validation_date, end_validation_date_6m)) garch_backtest <- pmap_dfr(garch_date_backtest, f_backtest) garch_backtest %>% mutate_at("VaR.estimate", function(x){paste0(round(x * 100, digits = 2), "%")}) %>% datatable()
Hậu kiểm cho thấy ước lượng bằng mô hình r garch_labels cho ra kết quả an toàn
Nhìn chung cả hai phương pháp ước lượng bằng mô hình GARCH và phương pháp mô phỏng lịch sử đều cho kết quả hậu kiểm VaR an toàn. Tuy nhiên, có thể thấy chuỗi giá trị lợi suất tồn tại hiện tượng tự tương quan và phương sai sai số thay đổi, sau khi áp dụng mô hình r garch_labels đã triệt tiêu được hai hiện tượng này, do đó mô hình r garch_labels là mô hình phù hợp hơn so với phương pháp mô phỏng lịch sử để dự báo về mức giảm giá tối đa.
Ước lượng và dự báo mức giảm giá tối đa
Từ kết quả hậu kiểm ở trên, mô hình r garch_labels là mô hình phù hợp nhất để ước lượng và dự báo mức độ giảm giá tối đa với độ tin cậy 95% (VaR95%). Thực hiện lại ước lượng tham số cho mô hình r garch_labels trong khoảng thời gian từ r fm_start_estimation_date đến r fm_end_validation_date, sau đó lấy trung bình mức VaR95% trong 30 ngày quan sát gần nhất làm giá trị dự báo cho mức VaR95%.
Độ biến động trong 1 tháng, 3 tháng, 6 tháng tiếp theo được xác định từ độ biến động trong 1 ngày như sau:
- Độ biến động trong 1 tháng = Độ biến động trong 1 ngày x √21
- Độ biến động trong 3 tháng = Độ biến động trong 1 ngày x √63
- Độ biến động trong 6 tháng = Độ biến động trong 1 ngày x √126
# Re-estimated GARCH spec <- ugarchspec(variance.model = list(model = garch_model, garchOrder = garch_order), mean.model = list(armaOrder=c(0,0)), distribution.model = garch_distribution) fit <- ugarchfit(spec = spec, data = data_return_clean) daily_VaR <- sigma(fit) * qdist(distribution = garch_distribution, p = 0.05, mu = 0, sigma = 1, skew = coef(fit)["skew"], shape = coef(fit)["shape"]) daily_VaR %>% tail(n = 30) %>% mean() %>% as.data.frame() %>% set_names("VaR1day") %>% mutate(VaR1month = VaR1day* sqrt(21), VaR3month = VaR1day* sqrt(63), VaR6month = VaR1day* sqrt(126)) %>% set_names(c("VaR 1 day", "VaR 1 month next", "VaR 3 months next", "VaR 6 months next")) %>% mutate_if(is.numeric, function(x){paste0(round(x * 100, digits = 2), "%")}) %>% datatable(options = list(searching = FALSE))
Phương pháp ước tính độ biến động cho khoảng thời gian dài hơn 1 ngày từ độ biến động trong 1 ngày như trên theo công thức căn thời gian (square root of time formula) được áp dụng rộng rãi trong phân tích rủi ro thị trường (chẳng hạn, theo các hướng dẫn của Ủy ban Basel). Tuy nhiên, công thức được sử dụng như trên là công thức gần đúng và sai số sẽ càng lớn nếu khoảng thời gian ước tính càng dài. Đối với các kỳ dự báo dài (từ 1 tháng trở lên), phương pháp này sẽ đưa ra kết quả thận trọng hơn mức ước lượng chính xác, hay nói cách khác là đưa ra mức giảm giá tối đa tương ứng với từng mức độ tin cậy cao hơn mức giảm giá “thực” theo mô hình GARCH.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.