Đọc dữ liệu từ Excel

library(xlsx) 
## Warning: package 'xlsx' was built under R version 4.3.3
library(xts)
## Warning: package 'xts' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(quantmod)
## Warning: package 'quantmod' was built under R version 4.3.3
## Loading required package: TTR
## Warning: package 'TTR' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(PerformanceAnalytics)
## Warning: package 'PerformanceAnalytics' was built under R version 4.3.3
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
d1 <- read.xlsx("C:/Users/HP/Documents/VTH_K22/Dữ liệu Lịch sử VN Index.xlsx", sheetIndex = 1, header = T)
d1 <- d1[!is.na(d1$DATE), ]
d2 <- read.xlsx("C:/Users/HP/Documents/VTH_K22/Dữ liệu Lịch sử S&P 500 Buyback NTR.xlsx", sheetIndex = 1, header = T)

# Khai báo dữ liệu chuỗi thời gian, sắp xếp theo thời gian
VNI <- xts(d1[,-1], order.by = d1$DATE)
 SP500 <- xts(d2[,-1], order.by = d2$DATE)

#Gộp hai chuỗi giá dựa trên những ngày trùng nhau
d3 <- merge.xts(VNI,SP500,join = 'inner')

#Tính tỷ suất lợi nhuận theo công thức logt - logt-1
cr <- CalculateReturns(d3,method = 'log')
cr <- cr[-1]

Thống kê mô tả:

library(moments)
## 
## Attaching package: 'moments'
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
library(dplyr)
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.3.3
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
cr_df <- as.data.frame(cr)
a <- cr_df %>% summarise(Min = min(VNI),
                      Max = max(VNI),
                      Mean = mean(VNI),
                      StDev = sd(VNI),
                      Skewness = skewness(VNI),
                      Kurtosis = kurtosis(VNI))
b <- cr_df %>% summarise(Min = min(SP500),
                      Max = max(SP500),
                      Mean = mean(SP500),
                      StDev = sd(SP500),
                      Skewness = skewness(SP500),
                      Kurtosis = kurtosis(SP500))
m <- rbind(a,b)
rownames(m) <- c('VNI','SP500')
kable(m, format = 'pandoc', caption = 'Bảng 1: Thống kê mô tả chuỗi TSLN', table.attr = "style='width:100%;'") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed")): Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 1: Thống kê mô tả chuỗi TSLN
Min Max Mean StDev Skewness Kurtosis
VNI -0.0690762 0.0486002 0.0000871 0.0130956 -0.9422434 6.64908
SP500 -0.1460276 0.1147455 0.0003902 0.0154438 -0.8004484 17.05258

Bảng trên đây cho thấy kết quả thống kê mô tả của hai chuỗi tỷ suất sinh lợi của chỉ số VNINDEX đại diện cho TTCK Việt Nam và chỉ số SP500 đại diện cho TTCK Mỹ. Kết quả cho thấy, giá trị trung bình tỷ suất lợi nhuận của VNINDEX là 0.0000803 và của Mỹ là 0.000388. Chỉ số VNI có độ biến động (0.0109334) nhỏ hơn chỉ số Mỹ (0.000388), điều này chứng tỏ TTCK Mỹ giao động khá lớn và không ổn định so với TTCK Việt Nam, qua đó cho thấy rủi ro của TTCK Việt Nam thấp hơn. Hệ số Skewness của cả hai chỉ số đều nhỏ hơn 0 chứng tỏ rằng phân phối của cả hai chuỗi lợi suất đều lệch trái. Bên cạnh đó, các giá trị Kurtosis của cả hai đều lớn hơn rất nhiều so với mức 3, điều này chỉ ra rằng phân phối của cả hai chỉ số có đuôi dài và nhiều giá trị cực trị hơn so với phân phối chuẩn. Ta có thể trực quan hóa các mô tả cơ bản của hai chỉ số lợi suất thông qua biểu đồ hộp (Box-plot).

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
pivot_longer(cr_df, cols = everything(), names_to = "TSLN", values_to = "Value") %>% ggplot(aes(x = TSLN, y = Value)) +
  geom_boxplot(fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Hình 1: Biểu đồ hộp của TSLN VNI và S&P500",
       x = "Chỉ số",
       y = "Tỷ suất sinh lợi")

Biều đồ hiển thị bên trên cho ta thấy rằng ở cả hai chỉ số VNI thuộc TTCK Việt Nam và chỉ số S&P500 thuộc TTCK Mỹ đều có tỷ suất lợi nhuận chủ yếu tập trung quanh giá trị trung vị, với một số điểm ngoại lệ ở cả hai phía, điều này chỉ ra rằng có một số giá trị tỷ suất sinh lợi bất thường trong dữ liệu, nguyên nhân có thể đến từ cú sốc kinh tế của đại dịch Covid-19 vừa qua.

Kiểm định mô hình

library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(rugarch)
## Warning: package 'rugarch' was built under R version 4.3.3
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
library(FinTS)
## Warning: package 'FinTS' was built under R version 4.3.3
#Kiểm định tính dừng
adf_vni <- adf.test(cr$VNI)
## Warning in adf.test(cr$VNI): p-value smaller than printed p-value
adf_SP500 <- adf.test(cr$SP500)
## Warning in adf.test(cr$SP500): p-value smaller than printed p-value
#Kiểm định phân phối chuẩn
jq_vni <- jarque.bera.test(cr$VNI)
jq_SP500 <- jarque.bera.test(cr$SP500)
# Ước lượng và trích xuất phần dư từ mô hình ARMA tối ưu
arima_VNI <- autoarfima(cr$VNI,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

arima_SP500 <- autoarfima(cr$SP500,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

re_VNI <- arima_VNI$fit@fit$residuals
re_SP500 <- arima_SP500$fit@fit$residuals
#Kiểm định tương quan chuỗi bậc 2 cho phần dư
lj_vni <- Box.test(re_VNI,type = 'Ljung-Box', lag = 2)
lj_SP500 <- Box.test(re_SP500,type = 'Ljung-Box', lag = 2)

#Kiểm định tương quan chuỗi bậc 2 cho phần dư bình phương
lj_vni2 <- Box.test(re_VNI^2,type = 'Ljung-Box', lag = 2)
lj_SP5002 <- Box.test(re_SP500^2,type = 'Ljung-Box', lag = 2)

#Kiểm định hiệu ứng ARCH
ar_vni <- ArchTest(re_VNI, lags = 2)
ar_SP500 <- ArchTest(re_SP500, lags = 2)

#Trình bày kết quả
test_result <- data.frame(Test = c("ADF","J-B","Q(2)","Q(2)^2", "ARCH(2)"),
                          P_value_VNI = c(adf_vni$p.value, jq_vni$p.value, lj_vni$p.value, lj_vni2$p.value, ar_vni$p.value),
                          P_value_SP500 = c(adf_SP500$p.value, jq_SP500$p.value, lj_SP500$p.value, lj_SP5002$p.value, ar_SP500$p.value))
kable(test_result, 
  caption = "Bảng 3: Kết quả các kiểm định", 
  label = 'Ghi chú: Q (2) and Q2 (2) lần lượt là kiểm định Ljung-Box Q2 cho tương quan chuỗi bậc 2 của phần dư và bình phương phần dư của lợi suất', 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 3: Kết quả các kiểm định
Test P_value_VNI P_value_SP500
ADF 0.0100000 0.0100000
J-B 0.0000000 0.0000000
Q(2) 0.9598288 0.2063587
Q(2)^2 0.0000000 0.0000000
ARCH(2) 0.0000000 0.0000000

Bảng trình bày các kết quả kiểm định của chuỗi tỷ suất sinh lợi và phần dư mô hình ARMA của nó. Kết quả cho thấy cả hai chuỗi tỷ suất đều là chuỗi dừng và không có phân phối chuẩn dựa trên kiểm định Augmented – Dickey – Fuller (ADF) và Jarque – Bera. Kiểm định Ljung-Box chứng tỏ phần dư của chuỗi VNI không có hiện tượng tương quan chuỗi ở cả chuỗi gốc và bình phương. Với phần dư chuỗi VNI và S&P500, kiểm định Ljung-Box cho chuỗi gốc (Q(2)) cho thấy không có hiện tượng tương quan chuỗi và kiểm định cho chuỗi bình phương (Q(2)^2) cho thấy hiện tượng tương quan chuỗi. Hiệu ứng ARCH được phát hiện ở phần dư của chuỗi S&P500 và chuỗi VNI cũng tồn tại hiệu ứng ARCH.

Phân tích hệ số tương quan

pearson <- cor(cr$VNI,cr$SP500, method="pearson")
spearman <- cor(cr$VNI,cr$SP500, method="spearman")
kendall <- cor(cr$VNI,cr$SP500, method="kendall")

#Trình bày kết quả
relat <- data.frame('Tương quan' = 'VNI-SP500',
                    Pearson = pearson[1,1],
                    spearman = spearman[1,1],
                    Kendall = kendall[1,1])
kable(relat, 
      col.names = c("Phương pháp", "Pearson", "Spearman","Kendall"),
      caption = "Bảng 2: Kết quả hệ số tương quan",
      format = 'pandoc',
      align = c("l", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 2: Kết quả hệ số tương quan
Phương pháp Pearson Spearman Kendall
VNI-SP500 0.1249314 0.0568061 0.0384619

Ta tiến hành trực quan hóa hệ số tương quan với phương pháp Pearson

library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.3.3
par(mfrow = c(1,2))

corr <- cor(cr)
ggcorrplot(corr, hc.order = TRUE,
   outline.col = "white",
   ggtheme = ggplot2::theme_gray,
   colors = c("lightblue", "white", "#E22229"),
   lab = TRUE, lab_col = 'white',title = 'Hình 2: Trực quan hóa hệ số tương quan với phương pháp Pearson')

chart.Correlation(cr, histogram=TRUE, pch=19)
## Warning in par(usr): argument 1 does not name a graphical parameter

Kết quả tính toán trên cho thấy hệ số tương quan của 3 phương pháp có mức chênh lệch nhau khá lớn, hệ số tương quan tuyến tính Pearson cho thấy mức độ tương quan lớn nhất với giá trị 0.1249314. Hai phương pháp còn lại là Spearman và Kendall cho kết quả có mức tương quan thấp hơn lần lượt là 0.0568061 và 0.0384619. Lưu ý rằng, hệ số tương quan tuyến tính Person giả định rằng các chuỗi lợi suất có phân phối chuẩn, nhưng thực tế thì lại không. Do đó, việc dùng hệ số tương quan để ước tính mối tương quan giữa 2 thị trường còn mang lại nhiều tranh cãi, và không thể giải thích được sự tương quan giữa các thị trường khi cả 2 biến động cực biên. Ngoài ra, hệ số tương quan hạng Spearman và Kendall không yêu cầu các chuỗi lợi suất có phân phối chuẩn. Hai hệ số tương quan này chưa phản ánh thông tin hoặc các cú sốc gây biến động thị trường, và chưa xem xét đến sự phụ thuộc đuôi của các thị trường. Do đó, để đánh giá khách quan hơn về mối quan hệ giữa hai thị trường, chúng ta sẽ tiến hành phương pháp Copula.

Phân tích sự biến động của chuỗi tỷ suất lợi nhuận

par(mfrow = c(1,2))
d11 <- CalculateReturns(d1, method = 'log')
ggplot(d1, aes(x = DATE, y= d11$VNI))+
  geom_line()+
  labs(title = "Hình 3: Biến động theo ngày của TSLN VNINDEX",x = "Ngày", y="Tỷ suất sinh lợi")+
  theme(plot.title = element_text(hjust = 0.5))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

d22 <- CalculateReturns(d2, method = 'log')

ggplot(d2, aes(x = DATE, y= d22$SP500))+
  geom_line()+
  labs(title = "Hình 4: Biến động theo ngày của TSLN SP500",x = "Ngày", y="Tỷ suất sinh lợi")+
  theme(plot.title = element_text(hjust = 0.5))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

Nhận xét chung: Hai hình ảnh hiển thị bên trên là 2 chuỗi biến động về tỷ suất sinh lợi ở chỉ số VNI và chỉ số S&P500 của Mỹ. Nhìn chung, chỉ số VNI có mức độ dao động thấp hơn chỉ số S&P500 của Mỹ, điều này chứng tỏ rằng rủi ro ở TTCK Việt Nam sẽ lơn hơn nhiều so với TTCK Mỹ. Theo tìm hiểu thì có lẽ biến động này diễn ra là do khoảng thời gian khó khăn khi diễn ra cú sốc kinh tế vì đại dịch Covid-19. Ở chỉ số VNINDEX, tỷ suất sinh lợi dao động mạnh nhất nằm trong khoảng thời gian từ năm 2021 đến khoảng cuối năm 2021 khi lợi nhuận giảm tối đa ở mức 0.05, đây cũng là khoảng thời gian mà dịch Covid-19 và cuộc xung đột diễn ra. Còn ở các giai đoạn còn lại thì sự dao động không quá lớn và dần có sự ổn định. Với chỉ số S&P500, nhìn chung mức độ dao động cao và cũng thấy rõ được đỉnh điểm của sự dao động bất ổn nhất là khoảng 2021 đến 2022. Điều này cho thấy có tồn tại sự ảnh hưởng của 2 cú sốc kinh tế đến thị trường chứng khoán ở cả Mỹ và Việt Nam. Để có thể dễ dàng quan sát và xác định rõ xu hướng biến động cùng với nhau giữa hai TTCK, chúng ta sẽ sử dụng biểu đồ Scatter với đường hồi quy để biểu diễn mối liên hệ tương quan giữa hai TTCK một cách trực quan.

cr %>% ggplot(aes(VNI, SP500)) +
              geom_point(col = '#BB444435',shape = TRUE) + 
              geom_smooth(method = 'lm',se = T, col = '#CCF247') + #thêm đường hồi quy với phương pháp hồi quy tuyến tính  
              labs(title = 'Hình 5: Biểu đồ Scatter TSLN của SP500 và VNI', x = 'VNI', y = 'SP500')
## `geom_smooth()` using formula = 'y ~ x'

Biểu đồ scatter plot và đường hồi quy tuyến tính giữa chỉ số VNI của thị trường chứng khoán Việt Nam và chỉ số S&P500 của thị trường chứng khoán Mỹ cung cấp nhiều thông tin có giá trị về mối quan hệ giữa hai thị trường này. Nhìn vào biểu đồ scatter plot, ta có thể thấy rằng các điểm dữ liệu tương đối tập trung xung quanh một đường thẳng. Điều này cho thấy có mối tương quan tuyến tính khá chặt chẽ giữa VNI và S&P500. Khi chỉ số VNI tăng lên, chỉ số S&P500 cũng có xu hướng tăng theo. Mối quan hệ này phản ánh sự liên kết và tương tác giữa hai thị trường chứng khoán này. Đường hồi quy tuyến tính thể hiện rõ hơn mối quan hệ này. Độ dốc của đường hồi quy cho biết mức độ thay đổi của S&P500 khi VNI thay đổi một đơn vị. Khoảng tin cậy 95% quanh đường hồi quy cũng cung cấp thông tin về độ chính xác của mối quan hệ tuyến tính này. Việc phân tích mối tương quan giữa VNI và S&P500 có ý nghĩa quan trọng đối với các nhà đầu tư và phân tích thị trường chứng khoán. Nó giúp họ hiểu rõ hơn về mối liên kết giữa hai thị trường, từ đó có thể đưa ra các chiến lược đầu tư và phòng ngừa rủi ro hiệu quả hơn. Đồng thời, phân tích này cũng cung cấp thông tin hữu ích cho việc so sánh, đánh giá và dự báo diễn biến của hai thị trường chứng khoán này.

Kết quả ước lượng cấu trúc phụ thuộc dựa trên phương pháp Copula có điều kiện

Xác định mô hình phân phối biên

#Setup garch(1,1) and eSP500mate for VNI coef
#Norm
spec_VNI_ugarch11_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch11_norm <- ugarchfit(spec = spec_VNI_ugarch11_norm, data = cr$VNI)

#Std
spec_VNI_ugarch11_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch11_std <- ugarchfit(spec = spec_VNI_ugarch11_std, data = cr$VNI)

#Sstd
spec_VNI_ugarch11_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch11_sstd <- ugarchfit(spec = spec_VNI_ugarch11_sstd, data = cr$VNI)

#ged
spec_VNI_ugarch11_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch11_ged <- ugarchfit(spec = spec_VNI_ugarch11_ged, data = cr$VNI)

#Sged
spec_VNI_ugarch11_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch11_sged <- ugarchfit(spec = spec_VNI_ugarch11_sged, data = cr$VNI)

#Setup garch(1,2) and eSP500mate for VNI coef
#norm
spec_VNI_ugarch12_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch12_norm <- ugarchfit(spec = spec_VNI_ugarch12_norm, data = cr$VNI)

#std
spec_VNI_ugarch12_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch12_std <- ugarchfit(spec = spec_VNI_ugarch12_std, data = cr$VNI)

#sstd
spec_VNI_ugarch12_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch12_sstd <- ugarchfit(spec = spec_VNI_ugarch12_sstd, data = cr$VNI)

#ged
spec_VNI_ugarch12_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch12_ged <- ugarchfit(spec = spec_VNI_ugarch12_ged, data = cr$VNI)

#sged
spec_VNI_ugarch12_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch12_sged <- ugarchfit(spec = spec_VNI_ugarch12_sged, data = cr$VNI)

#Setup garch(2,1) and estimate for VNI coef
#norm
spec_VNI_ugarch21_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch21_norm <- ugarchfit(spec = spec_VNI_ugarch21_norm, data = cr$VNI)

#std
spec_VNI_ugarch21_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch21_std <- ugarchfit(spec = spec_VNI_ugarch21_std, data = cr$VNI)

#sstd
spec_VNI_ugarch21_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch21_sstd <- ugarchfit(spec = spec_VNI_ugarch21_sstd, data = cr$VNI)

#ged
spec_VNI_ugarch21_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch21_ged <- ugarchfit(spec = spec_VNI_ugarch21_ged, data = cr$VNI)

#sged
spec_VNI_ugarch21_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch21_sged <- ugarchfit(spec = spec_VNI_ugarch21_sged, data = cr$VNI)

#Setup garch(2,2) and estimate for VNI coef
#norm
spec_VNI_ugarch22_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch22_norm <- ugarchfit(spec = spec_VNI_ugarch22_norm, data = cr$VNI)

#std
spec_VNI_ugarch22_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch22_std <- ugarchfit(spec = spec_VNI_ugarch22_std, data = cr$VNI)

#sstd
spec_VNI_ugarch22_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch22_sstd <- ugarchfit(spec = spec_VNI_ugarch22_sstd, data = cr$VNI)

#ged
spec_VNI_ugarch22_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch22_ged <- ugarchfit(spec = spec_VNI_ugarch22_ged, data = cr$VNI)

#sged
spec_VNI_ugarch22_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch22_sged <- ugarchfit(spec = spec_VNI_ugarch22_sged, data = cr$VNI)

#Setup garch(1,1) and estimate for SP500 coef
#norm
spec_SP500_ugarch11_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

SP500_garch11_norm <- ugarchfit(spec = spec_SP500_ugarch11_norm, data = cr$SP500)

#std
spec_SP500_ugarch11_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

SP500_garch11_std <- ugarchfit(spec = spec_SP500_ugarch11_std, data = cr$SP500)

#sstd
spec_SP500_ugarch11_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

SP500_garch11_sstd <- ugarchfit(spec = spec_SP500_ugarch11_sstd, data = cr$SP500)

#ged
spec_SP500_ugarch11_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

SP500_garch11_ged <- ugarchfit(spec = spec_SP500_ugarch11_ged, data = cr$SP500)

#sged
spec_SP500_ugarch11_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

SP500_garch11_sged <- ugarchfit(spec = spec_SP500_ugarch11_sged, data = cr$SP500)

#Setup garch(1,2) and eSP500mate for SP500 coef
#norm
spec_SP500_ugarch12_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

SP500_garch12_norm <- ugarchfit(spec = spec_SP500_ugarch12_norm, data = cr$SP500)

#std
spec_SP500_ugarch12_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

SP500_garch12_std <- ugarchfit(spec = spec_SP500_ugarch12_std, data = cr$SP500)

#sstd
spec_SP500_ugarch12_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

SP500_garch12_sstd <- ugarchfit(spec = spec_SP500_ugarch12_sstd, data = cr$SP500)

#ged
spec_SP500_ugarch12_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

SP500_garch12_ged <- ugarchfit(spec = spec_SP500_ugarch12_ged, data = cr$SP500)

#sged
spec_SP500_ugarch12_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

SP500_garch12_sged <- ugarchfit(spec = spec_SP500_ugarch12_sged, data = cr$SP500)

#Setup garch(2,1) and eSP500mate for SP500 coef
#norm
spec_SP500_ugarch21_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

SP500_garch21_norm <- ugarchfit(spec = spec_SP500_ugarch21_norm, data = cr$SP500)

#std
spec_SP500_ugarch21_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

SP500_garch21_std <- ugarchfit(spec = spec_SP500_ugarch21_std, data = cr$SP500)

#sstd
spec_SP500_ugarch21_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

SP500_garch21_sstd <- ugarchfit(spec = spec_SP500_ugarch21_sstd, data = cr$SP500)

#ged
spec_SP500_ugarch21_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

SP500_garch21_ged <- ugarchfit(spec = spec_SP500_ugarch21_ged, data = cr$SP500)

#sged
spec_SP500_ugarch21_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

SP500_garch21_sged <- ugarchfit(spec = spec_SP500_ugarch21_sged, data = cr$SP500)

#Setup garch(2,2) and eSP500mate for SP500 coef
#norm
spec_SP500_ugarch22_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

SP500_garch22_norm <- ugarchfit(spec = spec_SP500_ugarch22_norm, data = cr$SP500)

#std
spec_SP500_ugarch22_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

SP500_garch22_std <- ugarchfit(spec = spec_SP500_ugarch22_std, data = cr$SP500)

#sstd
spec_SP500_ugarch22_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

SP500_garch22_sstd <- ugarchfit(spec = spec_SP500_ugarch22_sstd, data = cr$SP500)

#ged
spec_SP500_ugarch22_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

SP500_garch22_ged <- ugarchfit(spec = spec_SP500_ugarch22_ged, data = cr$SP500)

#sged
spec_SP500_ugarch22_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

SP500_garch22_sged <- ugarchfit(spec = spec_SP500_ugarch22_sged, data = cr$SP500)

#Choose optimization model for vni coef
VNI.model.list <- list(garch11n = VNI_garch11_norm, garch11t = VNI_garch11_std, garch11st = VNI_garch11_sstd, 
                       garch11g = VNI_garch11_ged, garch11sg = VNI_garch11_sged, 
                       garch12n = VNI_garch12_norm, garch12t = VNI_garch12_std, garch12st = VNI_garch12_sstd, 
                       garch12g = VNI_garch12_ged, garch12sg = VNI_garch12_sged,
                       garch21n = VNI_garch21_norm, garch21t = VNI_garch21_std, garch21st = VNI_garch21_sstd, 
                       garch21g = VNI_garch21_ged, garch21sg = VNI_garch21_sged,
                       garch22n = VNI_garch22_norm, garch22t = VNI_garch22_std, garch22st = VNI_garch22_sstd, 
                       garch22g = VNI_garch22_ged, garch22sg = VNI_garch22_sged
                       ) 

VNI.info.mat <- sapply(VNI.model.list, infocriteria)

rownames(VNI.info.mat) <- rownames(infocriteria(VNI_garch11_norm))

VNI.inds <- which(VNI.info.mat == min(VNI.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.VNI <- colnames(VNI.info.mat)[VNI.inds[,2]]


#Choose optimization model for SP500 coef
SP500.model.list <- list(garch11n = SP500_garch11_norm, garch11t = SP500_garch11_std, garch11st = SP500_garch11_sstd, 
                       garch11g = SP500_garch11_ged, garch11sg = SP500_garch11_sged, 
                       garch12n = SP500_garch12_norm, garch12t = SP500_garch12_std, garch12st = SP500_garch12_sstd, 
                       garch12g = SP500_garch12_ged, garch12sg = SP500_garch12_sged,
                       garch21n = SP500_garch21_norm, garch21t = SP500_garch21_std, garch21st = SP500_garch21_sstd, 
                       garch21g = SP500_garch21_ged, garch21sg = SP500_garch21_sged,
                       garch22n = SP500_garch22_norm, garch22t = SP500_garch22_std, garch22st = SP500_garch22_sstd, 
                       garch22g = SP500_garch22_ged, garch22sg = SP500_garch22_sged
                       )

SP500.info.mat <- sapply(SP500.model.list, infocriteria)

rownames(SP500.info.mat) <- rownames(infocriteria(SP500_garch11_norm))

SP500.inds <- which(SP500.info.mat == min(SP500.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.SP500 <- colnames(SP500.info.mat)[SP500.inds[,2]]
#Trình bày kết quả
mar_model <- data.frame(
  rate = c("VNI", "SP500"),
  Format = c("ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t", "ARMA(0,2)-GJR-GARCH(1,1)-Student t"))

# Render the table
kable(mar_model, col.names = c("Tỷ suất sinh lợi", "Dạng mô hình phân phối biên"), 
      caption = "Bảng 5: Mô hình phân phối biên tối ưu", format = "pandoc", 
      table.attr = "style='width:100%;'") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed")): Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 5: Mô hình phân phối biên tối ưu
Tỷ suất sinh lợi Dạng mô hình phân phối biên
VNI ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t
SP500 ARMA(0,2)-GJR-GARCH(1,1)-Student t

Dựa trên tiêu chuẩn giá trị AIC nhỏ nhất, kết quả thu được dạng mô hình phân phối biên cho tỷ suất sinh lợi của VNI và SP500 lần lượt là mô hình ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t và ARMA(0,2)-GJR-GARCH(1,1)-Student t. Điều này một lần nữa xác nhận rằng, hàm phân phối biên của tất cả các chuỗi lợi suất chứng khoán được xem xét không tuân theo phân phối chuẩn.

#Trình bày kết quả

extract_garch_results <- function(fit) {
  coef <- coef(fit)
  se <- fit@fit$se.coef
  pvalues <- 2 * (1 - pnorm(abs(fit@fit$tval))) # tính giá trị p từ giá trị t
  results <- cbind(coef, se, pvalues)
  colnames(results) <- c("Estimate", "Std. Error", "Pr(>|z|)")
  return(results)
}

fit1 <- extract_garch_results(VNI_garch11_sstd)
fit2 <- extract_garch_results(SP500_garch11_std)

kable(as.data.frame(fit1), 
  caption = "Bảng 3: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Skewed Student của biến VNI", 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 3: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Skewed Student của biến VNI
Estimate Std. Error Pr(>|z|)
mu 0.0002241 0.0003846 0.5601985
ar1 0.9330998 0.0455384 0.0000000
ma1 -0.9041205 0.0551232 0.0000000
omega 0.0000108 0.0000004 0.0000000
alpha1 0.0285809 0.0133323 0.0320540
beta1 0.7931812 0.0197323 0.0000000
gamma1 0.2349827 0.0489172 0.0000016
skew 0.8220094 0.0305851 0.0000000
shape 4.1728728 0.4137612 0.0000000

Bảng 5 kết quả mô hình ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student cho tỷ suất sinh lợi của chỉ số VNI cho thấy một số đặc điểm quan trọng về cấu trúc và động thái của dữ liệu. Tham số AR(1) có giá trị dương gần 1 và có ý nghĩa thống kê cao, cho thấy sự tự duy trì mạnh mẽ của lợi nhuận qua các thời kỳ. Ngược lại, tham số MA(1) có giá trị âm và cũng rất có ý nghĩa, điều này ngụ ý rằng các cú sốc quá khứ có ảnh hưởng tiêu cực mạnh mẽ đến lợi nhuận hiện tại. Trong phần mô hình GARCH, tham số omega rất nhỏ nhưng có ý nghĩa, biểu thị mức độ biến động cơ bản thấp. Tham số beta1 cao và có ý nghĩa, cho thấy sự duy trì mạnh mẽ của biến động qua thời gian. Tham số alpha1 không có ý nghĩa thống kê, nghĩa là các cú sốc quá khứ không ảnh hưởng nhiều đến biến động hiện tại. Tuy nhiên, tham số gamma1 lại có ý nghĩa, cho thấy hiệu ứng đòn bẩy, nơi mà các cú sốc tiêu cực có ảnh hưởng mạnh hơn đến biến động so với các cú sốc tích cực có cùng độ lớn. Cuối cùng, tham số skew và shape đều có ý nghĩa thống kê, chỉ ra rằng phân phối sai số có độ lệch và đuôi dày hơn so với phân phối chuẩn. Cụ thể, độ lệch dương cho thấy phân phối bị lệch về phía phải, trong khi giá trị shape cho thấy sự hiện diện của các đuôi nặng, ám chỉ rằng các giá trị cực đoan xảy ra thường xuyên hơn so với phân phối chuẩn.

kable(as.data.frame(fit2), 
  caption = "Bảng 4: Kết quả mô hình ARMA(0,2)-GJR-Garch(1,1)-Student của biến SP500", 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 4: Kết quả mô hình ARMA(0,2)-GJR-Garch(1,1)-Student của biến SP500
Estimate Std. Error Pr(>|z|)
mu 0.0005908 0.0002492 0.0177649
ma1 -0.0060688 0.0278055 0.8272264
ma2 0.0072306 0.0282396 0.7979164
omega 0.0000039 0.0000022 0.0801615
alpha1 0.0413252 0.0212955 0.0523111
beta1 0.8339081 0.0265387 0.0000000
gamma1 0.2195926 0.0434293 0.0000004
shape 7.8950122 1.5426984 0.0000003

Với mô hình ARMA(0,2)-GJR-GARCH(1,1)-Student cho tỷ suất sinh lợi của chỉ số SP500 cho thấy tham số MA(1) và MA(2) đều không có ý nghĩa thống kê, nghĩa là các cú sốc quá khứ gần đây không có ảnh hưởng đáng kể đến lợi nhuận hiện tại. Đối với mô hình GARCH, tham số omega rất nhỏ nhưng có ý nghĩa thống kê, biểu thị mức độ biến động cơ bản thấp. Tương tự như trên, tham số alpha1 của mô hình SP500 cũng không có ý nghĩa thống kê, cho thấy các cú sốc quá khứ không ảnh hưởng nhiều đến biến động hiện tại. Ngược lại, tham số beta1 cao và có ý nghĩa, cho thấy sự duy trì mạnh mẽ của biến động qua thời gian. Tham số gamma1 có ý nghĩa thống kê, chỉ ra hiệu ứng đòn bẩy, nơi mà các cú sốc tiêu cực có ảnh hưởng mạnh hơn đến biến động so với các cú sốc tích cực có cùng độ lớn. Bên cạnh đó, tham số shape cũng chỉ ra rằng phân phối sai số có đuôi dày hơn so với phân phối chuẩn, ám chỉ rằng các giá trị cực đoan xảy ra thường xuyên hơn. Tóm lại, kết quả mô hình cho thấy rằng mặc dù các thành phần trung bình và cú sốc quá khứ gần đây không có ảnh hưởng đáng kể đến lợi nhuận hiện tại, nhưng mức độ biến động có tính duy trì cao và có hiệu ứng đòn bẩy đáng kể. Phân phối sai số cũng cho thấy đuôi dày, ngụ ý rằng các biến động lớn xảy ra thường xuyên hơn so với phân phối chuẩn.

Kiểm định tính phù hợp của mô hình phân phối biên

# Tính tích phân xác suất cho phần dư của VNI
VNI.res <- residuals(VNI_garch11_sstd) / sigma(VNI_garch11_sstd)  # Phần dư đã được chuẩn hóa

a <- fitdist(distribution = "sstd", VNI.res, control = list())  # Khớp phân phối cho phần dư này

u <- pdist(distribution = "sstd", q = VNI.res, mu = a$pars['mu'], sigma = a$pars['sigma'],
           skew = a$pars['skew'], shape = a$pars['shape'])  # Tính tích phân xác suất

# Tính tích phân xác suất cho phần dư của SP500
SP500.res <- residuals(SP500_garch11_std) / sigma(SP500_garch11_std)  # Phần dư đã được chuẩn hóa

b <- fitdist(distribution = "std", SP500.res, control = list())  # Khớp phân phối cho phần dư này

v <- pdist(distribution = "std", q = SP500.res, mu = b$pars['mu'], sigma = b$pars['sigma'],
           skew = b$pars['skew'], shape = b$pars['shape'])  # Tính tích phân xác suất

# Kiểm định với Anderson-Darling
library(goftest)
ad_vni <- ad.test(u, "punif")
ad_SP500 <- ad.test(v, "punif")

# Kiểm định Cramer-von Mises
cvm_vni <- cvm.test(u, "punif")
cvm_SP500 <- cvm.test(v, "punif")

# Kiểm định Kolmogorov-Smirnov
ks_vni <- ks.test(u, "punif")
ks_SP500 <- ks.test(v, "punif")

# Trình bày kết quả
test <- data.frame(test = c('P_value VNI', 'P_value SP500'),
                   AD = c(ad_vni$p.value, ad_SP500$p.value),
                   CVM = c(cvm_vni$p.value, cvm_SP500$p.value),
                   KS = c(ks_vni$p.value, ks_SP500$p.value))

kable(test, 
      caption = "Kết quả các kiểm định cho mô hình phân phối biên", 
      col.names = c("Các kiểm định", "Anderson-Darling", "Cramer-von Mises", "Kolmogorov-Smirnov"), 
      format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Kết quả các kiểm định cho mô hình phân phối biên
Các kiểm định Anderson-Darling Cramer-von Mises Kolmogorov-Smirnov
P_value VNI 0.4654821 0.3664421 0.1200949
P_value SP500 0.1912726 0.4778018 0.5389404

Bảng 7 trình bày kết quả các kiểm định phân phối cho mô hình của hai biến VNI và S&P500, với ba loại kiểm định: Anderson-Darling, Cramer-von Mises và Kolmogorov-Smirnov. Các giá trị P_value cao và lớn hơn mức ý nghĩa từ cả ba kiểm định cho thấy rằng không có bằng chứng thống kê để bác bỏ giả thuyết rằng dữ liệu của các chuỗi tỷ suất VNI và S&P500 tuân theo phân phối được giả định trong mô hình. Cụ thể, giá trị P_value của kiểm định Anderson-Darling là 0.4654821 cho VNI và 0.1912726 cho S&P500, giá trị P_value của kiểm định Cramer-von Mises là 0.3664421 cho VNI và 0.4778018 cho S&P500, và giá trị P_value của kiểm định Kolmogorov-Smirnov là 0.1200949 cho VNI và 0.5389404 cho S&P500. Những kết quả này củng cố kết luận rằng mô hình Copula sẽ là phù hợp với dữ liệu của hai biến này.

Ước lượng tham số mô hình copula

Sau khi lựa chọn được mô hình biên phù hợp, em sẽ ước lượng mô hình Copula dựa trên phương pháp MLE, mô hình Copula tối ưu được lựa chọn dựa trên chỉ số AIC và BIC nhỏ nhất với hàm BiCopSelect nằm trong package VineCopula.

library(VineCopula)
## Warning: package 'VineCopula' was built under R version 4.3.3
library(copula)
## Warning: package 'copula' was built under R version 4.3.3
## 
## Attaching package: 'copula'
## The following object is masked from 'package:VineCopula':
## 
##     pobs
#Ước lượng
gau <- BiCopEst(u,v,family = 1, method = "mle", se = T, max.df = 10)
stu <- BiCopEst(u,v,family = 2, method = "mle", se = T, max.df = 10)
#Trình bày kết quả

est <- data.frame(mqh = c('VNI-SP500' ,'VNI-SP500'),
                    Copula = c('Gauss', 'Student'),
                    Thamso = c(gau$par, stu$par),
                  Thamso2 = c(gau$par2, stu$par2),
                   duoi = c(gau$taildep$lower,stu$taildep$lower))

kable(est, 
  caption = "Bảng 7: Kết quả ước lượng Copula họ Elip", col.names = c("Chỉ số","Copula", "Par","Par2", "Phụ thuộc đuôi"), 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Bảng 7: Kết quả ước lượng Copula họ Elip
Chỉ số Copula Par Par2 Phụ thuộc đuôi
VNI-SP500 Gauss 0.1023965 0 0.0000000
VNI-SP500 Student 0.0911169 10 0.0115105

Kết quả ước lượng copula họ Elip trong điều kiện thị trường bình thường cho mối quan hệ giữa hai biến VNI và ST&P500, bao gồm hai loại copula: Gaussian và Student. Kết quả từ copula Gaussian cho thấy một mức độ tương quan dương giữa hai biến với giá trị Par là 0.1023965, nhưng không thể mô hình hóa phụ thuộc đuôi, như thể hiện qua giá trị phụ thuộc đuôi là 0. Ngược lại, copula Student không chỉ cho thấy mức độ tương quan dương tương tự với giá trị Par là 0.0911169, mà còn có khả năng mô hình hóa các mối quan hệ phi tuyến và phụ thuộc đuôi. Tham số Par2 của copula Student là 10, cho thấy độ tự do của phân phối Student t, và giá trị phụ thuộc đuôi là 0.0115105, biểu thị mức độ phụ thuộc đuôi dương giữa VNI và S&P500. Điều này ngụ ý rằng có một xác suất nhất định rằng hai biến sẽ cùng trải qua các biến động cực đoan theo cùng một hướng. Cụ thể, nếu TTCK Mỹ biến động tăng (giảm) giá thì khả năng TTCK Việt Nam tăng (giảm) theo là khoảng 1.15%.

Biểu đồ các họ Coupula

library(knitr)
library(copula)
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 4.3.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.3.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(VC2copula)
## Warning: package 'VC2copula' was built under R version 4.3.3
## 
## Attaching package: 'VC2copula'
## The following objects are masked from 'package:VineCopula':
## 
##     BB1Copula, BB6Copula, BB7Copula, BB8Copula, copulaFromFamilyIndex,
##     joeBiCopula, r270BB1Copula, r270BB6Copula, r270BB7Copula,
##     r270BB8Copula, r270ClaytonCopula, r270GumbelCopula,
##     r270JoeBiCopula, r270TawnT1Copula, r270TawnT2Copula, r90BB1Copula,
##     r90BB6Copula, r90BB7Copula, r90BB8Copula, r90ClaytonCopula,
##     r90GumbelCopula, r90JoeBiCopula, r90TawnT1Copula, r90TawnT2Copula,
##     surBB1Copula, surBB6Copula, surBB7Copula, surBB8Copula,
##     surClaytonCopula, surGumbelCopula, surJoeBiCopula, surTawnT1Copula,
##     surTawnT2Copula, tawnT1Copula, tawnT2Copula, vineCopula
library(VineCopula)
library(gridGraphics)
## Warning: package 'gridGraphics' was built under R version 4.3.3
## Loading required package: grid
library(png)
#Tạo hàm vẽ biểu đồ
scatter_plot <- function(random_data, cl) {
  ggplot(data.frame(random_data), aes(random_data[,1], random_data[,2])) +
  geom_point(alpha = 0.5, col = cl) +
  theme_minimal() +
  labs(x = "u", y = "v")
  } #for scatter

persp_plot <- function(copula_obj, file_name, cl) {
  png(file_name)
  persp(copula_obj, dCopula, 
        xlab = 'u', ylab = 'v', col = cl, ltheta = 120,  
        ticktype = "detailed", cex.axis = 0.8)
  dev.off()
  rasterGrob(readPNG(file_name),interpolate = TRUE)
}# for pdf

set.seed(123)
#Mô phỏng copula gauss với p=0.8
cop_nor <- normalCopula(param = 0.8, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_nor <- rCopula(copula = cop_nor,n = 7600)
#Mô phỏng copula với p=0.8
cop_std <- tCopula(param = 0.8, dim = 2, df = 1)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_std <- rCopula(copula = cop_std,n = 7600)

Copula Gauss và Student

library(ggplot2)
library(plotly)
#Vẽ biểu đồ 
png('nors.png') #Lưu ảnh biểu đồ
scatter_plot(random_nor, '#99FFCC')

dev.off()
## png 
##   2
nors <- rasterGrob(readPNG('nors.png'), interpolate = TRUE) #Chuyển ảnh sang dạng Grob

png('stds.png')
scatter_plot(random_std, '#CC5333')
dev.off()
## png 
##   2
stds <- rasterGrob(readPNG('stds.png'), interpolate = TRUE)
nor_per <- persp_plot(cop_nor, 'norp.png', '#99FFCC')
std_per <- persp_plot(cop_std, 'stdp.png', '#CC5333')

legend <- legendGrob(
  labels = c("Gauss", "Student"), pch = 15,
  gp = gpar(col = c('#99FFCC', '#CC5333'), fill = c('#99FFCC', '#CC5333'))
)

grid.arrange(nors, nor_per, stds, std_per, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 1: Biểu đồ phân tán và phối cảnh PDF của Copula họ Elip", 
                 gp = gpar(fontsize = 15, font = 2))
             )

Copula Clayton và Gumbel

set.seed(123)
#Mô phỏng copula clayton với p=4
cop_clay <- claytonCopula(param = 4, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_clay <- rCopula(copula = cop_clay,n = 7600)
#Mô phỏng copula gumbel với p=5
cop_gum <- gumbelCopula(param = 5, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_gum <- rCopula(copula = cop_gum,n = 7600)
#Vẽ biểu đồ 
png('clays.png')
scatter_plot(random_clay,'gray')
dev.off()
## png 
##   2
clays <- rasterGrob(readPNG('clays.png'), interpolate = TRUE)

png('gums.png')
scatter_plot(random_gum,'red')
dev.off()
## png 
##   2
gums <- rasterGrob(readPNG('gums.png'), interpolate = TRUE)
clayp <- persp_plot(cop_clay,'clayp.png','gray')

gump <- persp_plot(cop_gum,'gump.png','red')


legend <- legendGrob(
  labels = c("Clayton", "Gumbel"), pch = 15,
  gp = gpar(col = c('gray', 'red'), fill = c('gray', 'red'))
)

grid.arrange(clays, clayp, gums, gump, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 2: Biểu đồ phân tán và PDF của Copula Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

Copula Franl và Joe

set.seed(123)
#Mô phỏng copula Frank
cop_frank <- frankCopula(param = 9.2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_frank <- rCopula(copula = cop_frank,n = 7600)
#Mô phỏng copula survival clayton  với p=4
cop_joe <- joeCopula(param = 3)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_joe <- rCopula(copula = cop_joe,n = 7600)
png('franks.png')
scatter_plot(random_frank,'darkblue')
dev.off()
## png 
##   2
franks <- rasterGrob(readPNG('franks.png'), interpolate = TRUE)

png('joes.png')
scatter_plot(random_joe,'pink')
dev.off()
## png 
##   2
joes <- rasterGrob(readPNG('joes.png'), interpolate = TRUE)
frankp <- persp_plot(cop_frank, "frankp.png",'darkblue')
joep <- persp_plot(cop_joe, "joep.png",'pink')

legend <- legendGrob(labels = c("Franl","Joe"), pch = 15,
                    gp = gpar(col = c('darkblue','pink'), fill = c('darkblue','pink')))

grid.arrange(franks, frankp,joes, joep, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 4: Biểu đồ phân tán và  PDF của Copula Franl và Joe", 
                 gp = gpar(fontsize = 15, font = 2))
             )

Lựa chọn và ước lượng

#Lựa chọn và ước lượng
opt_cop <- BiCopSelect(u,v,selectioncrit = "AIC",method = "mle")
est_opt_cop <- BiCopEst(u,v,family = opt_cop$family,method = "mle",max.df = 30)

#Trình bày kết quả
est_opt_cop_dt <- data.frame(mqh = c('VNI-SP500'),
                    Copula = est_opt_cop$familyname,
                    Thamso = est_opt_cop$par,
                  Thamso2 = est_opt_cop$par2,
                   duoiduoi = est_opt_cop$taildep$lower,
                  duoitren = est_opt_cop$taildep$upper,
                  tau = est_opt_cop$tau,
                  aic = est_opt_cop$AIC,
                  bic = est_opt_cop$BIC)

kable(est_opt_cop_dt, 
  caption = "Kết quả ước lượng mô hình Copula tối ưu", col.names = c("Chỉ số","Copula", "Par","Par2", "Đuôi dưới", "Đuôi trên","Hệ số Kendall", "AIC","BIC"), 
  format = 'pandoc') %>% kable_styling(
bootstrap_options = c("striped", "hover", "condensed"), 
full_width = F)
## Warning in kable_styling(., bootstrap_options = c("striped", "hover",
## "condensed"), : Please specify format in kable. kableExtra can customize either
## HTML or LaTeX outputs. See https://haozhu233.github.io/kableExtra/ for details.
Kết quả ước lượng mô hình Copula tối ưu
Chỉ số Copula Par Par2 Đuôi dưới Đuôi trên Hệ số Kendall AIC BIC
VNI-SP500 Survival Gumbel 1.058816 0 0.0755428 0 0.0555484 -25.41188 -20.13187 
persp(VC2copula::BB7Copula(param = c(est_opt_cop$par, est_opt_cop$par2)),dCopula, 
        xlab = 'u', ylab = 'v', col = 'yellow', ltheta = 120,  
        ticktype = "detailed", cex.axis = 0.8,  main = 'Phối cảnh PDF của mô hình Copula Gumbel')

Đáng chú ý, phân bố xác suất không đối xứng, có phần nghiêng về phía trái, điều này có thể cho thấy rằng các biến động của chỉ số MSCI có xu hướng lớn hơn so với biến động của chỉ số VNI. Điều này có thể được giải thích bởi sự ảnh hưởng của các yếu tố vĩ mô và địa chính trị toàn cầu đối với thị trường chứng khoán quốc tế (thể hiện qua MSCI) so với thị trường chứng khoán Việt Nam (thể hiện qua VNI). Ngoài ra, ta cũng có thể quan sát thấy rằng phân bố xác suất có vùng đuôi dày, điều này có thể cho thấy rằng xác suất xảy ra các biến động cùng chiều lớn (tăng hoặc giảm đồng thời) giữa hai chỉ số này là khá cao. Điều này phù hợp với đặc điểm của thị trường tài chính, khi các biến động lớn thường xảy ra đồng thời trên nhiều thị trường.

library(plotly)
library(ggplot2)
set.seed(123)

# Tham số cho copula Gumbel (ví dụ với theta = 2)
theta_gumbel <- 2

# Tạo copula Gumbel
cop_gumbel <- gumbelCopula(param = theta_gumbel, dim = 2)

# Mô phỏng dữ liệu từ copula Gumbel

# Kích thước mẫu
sample_size <- 500  
sample_gumbel <- rCopula(sample_size, cop_gumbel)

# Hiển thị một số mẫu của copula Gumbel
head(sample_gumbel)
##            [,1]       [,2]
## [1,] 0.08097408 0.04203784
## [2,] 0.87624894 0.98004472
## [3,] 0.45315700 0.48634119
## [4,] 0.95903035 0.81314647
## [5,] 0.71034376 0.88160340
## [6,] 0.82655551 0.38809806
# [1,] 0.08097408 0.04203784
# [2,] 0.87624894 0.98004472
# [3,] 0.45315700 0.48634119
# [4,] 0.95903035 0.81314647
# [5,] 0.71034376 0.88160340
# [6,] 0.82655551 0.38809806
set.seed(123)

# Tham số cho copula Gumbel
theta_gumbel <- 2

# Tạo copula Gumbel
cop_gumbel <- gumbelCopula(param = theta_gumbel, dim = 2)

# Mô phỏng dữ liệu từ copula Gumbel
sample_size <- 500
sample_gumbel <- rCopula(sample_size, cop_gumbel)

# Tạo lưới điểm để tính toán mật độ
x_seq <- seq(0, 1, length.out = 30)
y_seq <- seq(0, 1, length.out = 30)
grid <- expand.grid(x = x_seq, y = y_seq)

# Tính toán mật độ của copula Gumbel cho lưới điểm
pdf_values <- apply(grid, 1, function(p) {
  dCopula(p, cop_gumbel)
})
pdf_matrix <- matrix(pdf_values, nrow = length(x_seq), ncol = length(y_seq))

# Vẽ biểu đồ 3D của phân phối mật độ xác suất
plot_ly(z = pdf_matrix, x = x_seq, y = y_seq, type = "surface") %>%
  layout(title = "Phân phối mật độ xác suất của Copula Gumbel",
         scene = list(
           xaxis = list(title = "X"),
           yaxis = list(title = "Y"),
           zaxis = list(title = "Tỷ trọng ")
         ))

Hình ảnh này thể hiện một mô hình copula Gumbel, được sử dụng để phân tích mối quan hệ phụ thuộc giữa hai biến ngẫu nhiên. Cụ thể, chúng ta đang xem xét sự liên kết giữa chỉ số chứng khoán VNI của Việt Nam và chỉ số S&P500 của Mỹ. Copula Gumbel thuộc nhóm copula Archimedean và đặc biệt phù hợp để mô hình hóa mối quan hệ phụ thuộc có đuôi dày, tức là những trường hợp phụ thuộc mạnh mẽ giữa các biến. Trong hình, phân phối mật độ xác suất của copula Gumbel cho thấy sự tập trung rõ rệt ở góc trên bên trái, điều này chỉ ra sự phụ thuộc mạnh mẽ giữa hai chỉ số. Nói cách khác, khi chỉ số VNI của Việt Nam tăng, chỉ số S&P500 của Mỹ cũng có xu hướng tăng theo, và ngược lại. Mô hình copula Gumbel cho phép chúng ta phân tích mối quan hệ này với độ chính xác cao hơn so với các phương pháp truyền thống như hệ số tương quan. Hiểu rõ cấu trúc phụ thuộc giữa các chỉ số chứng khoán từ các quốc gia khác nhau giúp các nhà đầu tư và quản lý rủi ro đưa ra các quyết định đầu tư và chiến lược phòng ngừa rủi ro hiệu quả hơn. Mô hình copula Gumbel là một công cụ quan trọng trong việc phân tích và mô hình hóa các mối quan hệ đa biến trong tài chính.

# Tải các thư viện cần thiết
library(copula)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked _by_ '.GlobalEnv':
## 
##     SP500
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
library(scatterplot3d)

# Định nghĩa tham số cho copula Gumbel
theta <- 1.18

# Tạo đối tượng copula Gumbel
gumbel_copula <- gumbelCopula(param = theta)

# Dữ liệu ví dụ cho kde2d

set.seed(123)
cd <- cbind(rnorm(1000), rnorm(1000))

# Thực hiện ước lượng mật độ kernel
kde2d_result <- kde2d(cd[,1], cd[,2], n = 50)

# Chuyển đổi kết quả kde2d thành định dạng phù hợp với scatterplot3d
x <- kde2d_result$x
y <- kde2d_result$y
z <- kde2d_result$z

# Tạo lưới điểm cho scatterplot3d
grid <- expand.grid(x = x, y = y)

# Làm phẳng z cho scatterplot3d
z_flat <- as.vector(t(z))

# Đảm bảo grid và z_flat có cùng độ dài
if (length(z_flat) != nrow(grid)) {
  stop("Độ dài của z_flat không khớp với số lượng điểm trong grid")
}

# Vẽ đồ thị mật độ sử dụng scatterplot3d
scatterplot3d(grid$x, grid$y, z_flat, type = "h", angle = 30, 
              main = "Biểu đồ Mật độ của Copula Gumbel",
              xlab = "Rainfall Depth", 
              ylab = "Shot Duration", 
              zlab = "Phân phối hai biến",
              pch = 20, color = "darkgray")

Phân tích mối quan hệ giữa các biến trong thực tế có thể gặp khó khăn khi các biến không tuân theo phân phối chuẩn đơn giản. Trong những tình huống này, mô hình Copula có thể cung cấp một phương pháp linh hoạt và hiệu quả để mô hình hóa các liên kết phức tạp giữa các biến. Đoạn mã R dưới đây minh họa cách sử dụng thư viện copula và MASS để thiết lập và vẽ biểu đồ mật độ dựa trên mô hình Copula Gumbel. Copula Gumbel là một công cụ thống kê mạnh mẽ, cho phép mô hình hóa mối quan hệ phi tuyến và phụ thuộc mạnh mẽ giữa các biến. Ví dụ này sử dụng dữ liệu mô phỏng ngẫu nhiên để tạo ra một tập dữ liệu mẫu, nhưng để phân tích mối quan hệ giữa các chỉ số thực tế như VNI và S&P500, cần sử dụng dữ liệu thực tế. Biểu đồ mật độ 3D cho thấy cấu trúc phân phối phức tạp của dữ liệu, điều này phản ánh rằng mối quan hệ giữa hai biến không phải là tuyến tính đơn giản. Sự xuất hiện của các điểm dữ liệu ở các cạnh của biểu đồ cho thấy sự phụ thuộc mạnh mẽ giữa hai biến trong những trường hợp cụ thể. Để đánh giá chính xác độ phù hợp của mô hình Copula Gumbel, cần thực hiện các bước bổ sung như ước lượng tham số, kiểm tra độ phù hợp của mô hình và so sánh với các mô hình khác. Tuy nhiên, biểu đồ mật độ 3D đã cung cấp cái nhìn trực quan về cấu trúc phân phối phức tạp và hỗ trợ việc sử dụng mô hình Copula Gumbel trong việc mô hình hóa mối quan hệ giữa các biến.