Mô hình ngẫu nhiên

Trích xuất dữ liệu

# Load các package cần thiết 
library(xlsx) 

## Warning: package 'xlsx' was built under R version 4.3.3

library(quantmod)

## Warning: package 'quantmod' was built under R version 4.3.3

## Loading required package: 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

## 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(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.3

## Warning: package 'ggplot2' was built under R version 4.3.3

## Warning: package 'tidyr' was built under R version 4.3.3

## Warning: package 'dplyr' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first()  masks xts::first()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::last()   masks xts::last()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(xts)
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

library(forecast)

## Warning: package 'forecast' 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:purrr':
## 
##     reduce
## 
## The following object is masked from 'package:stats':
## 
##     sigma

library(goftest)
library(VineCopula)

## Warning: package 'VineCopula' was built under R version 4.3.3

library(tseries)

## Warning: package 'tseries' was built under R version 4.3.3

library(FinTS)

## Warning: package 'FinTS' was built under R version 4.3.3

## 
## Attaching package: 'FinTS'
## 
## The following object is masked from 'package:forecast':
## 
##     Acf

library(ggplot2)
library(corrplot)

## Warning: package 'corrplot' was built under R version 4.3.3

## corrplot 0.92 loaded

library(ggcorrplot)

## Warning: package 'ggcorrplot' was built under R version 4.3.3

# Nhập dữ liệu 
data <- read.xlsx("C:/Users/ADMIN/Desktop/MHNN/dulieucd1.xlsx", sheetIndex = 1, header = T)

Tỷ suất lợi nhuận theo chuỗi data

# Tính tỷ suất lợi nhuận theo công thức log(t)-log(t-1)
lgvn<- diff(log(data$VNI), lag = 1)
lgtl<- diff(log(data$THAILAN), lag = 1)
mhnn <- data.frame(VNI= lgvn, TL = lgtl)

Thống kê mô tả

library(moments)

## 
## Attaching package: 'moments'

## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness

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

library(dplyr)
mhnn_df <- as.data.frame(mhnn)
a <- mhnn_df %>% summarise(Min = min(VNI),
                      Max = max(VNI),
                      Mean = mean(VNI),
                      StDev = sd(VNI),
                      Skewness = skewness(VNI),
                      Kurtosis = kurtosis(VNI))
b <- mhnn_df %>% summarise(Min = min(TL),
                      Max = max(TL),
                      Mean = mean(TL),
                      StDev = sd(TL),
                      Skewness = skewness(TL),
                      Kurtosis = kurtosis(TL))
m <- rbind(a,b)
rownames(m) <- c('VNI','TL')
kable(m, format = 'pandoc', caption = 'Bảng: 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: Thống kê mô tả chuỗi TSLN

	Min	Max	Mean	StDev	Skewness	Kurtosis
VNI	-0.0648198	0.0620016	0.0004976	0.0109334	-0.7978217	8.318537
TL	-0.1140605	0.0811544	-0.0000457	0.0118501	-1.2811032	21.163499

Đồ thị Box-Plot

# Biểu đồ Box-plot
pivot_longer(mhnn_df, cols = everything(), names_to = "TSLN", values_to = "Value") %>% ggplot(aes(x = TSLN, y = Value)) +
  geom_boxplot(fill = "cyan", color = "black") +
  theme_minimal() +
  labs(title = "Hình: Biểu đồ hộp của TSLN VNI và TL",
       x = "Chỉ số",
       y = "Tỷ suất sinh lợi")

Các kiểm định cần thiết

library(tseries)
# Các kiểm định 
#Kiểm định tính dừng
adf_vni <- adf.test(mhnn$VNI)

## Warning in adf.test(mhnn$VNI): p-value smaller than printed p-value

adf_tl <- adf.test(mhnn$TL)

## Warning in adf.test(mhnn$TL): p-value smaller than printed p-value

#Kiểm định phân phối chuẩn
jq_vni <- jarque.bera.test(mhnn$VNI)
jq_tl <- jarque.bera.test(mhnn$TL)
# Ước lượng và trích xuất phần dư từ mô hình ARMA tối ưu
arima_VNI <- autoarfima(mhnn$VNI,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

arima_tl <- autoarfima(mhnn$TL,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

re_VNI <- arima_VNI$fit@fit$residuals
re_TL <- arima_tl$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_tl <- Box.test(re_TL,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_tl2 <- Box.test(re_TL^2,type = 'Ljung-Box', lag = 2)

#Kiểm định hiệu ứng ARCH
ar_vni <- ArchTest(re_VNI, lags = 2)
ar_tl <- ArchTest(re_TL, 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_TL = c(adf_tl$p.value, jq_tl$p.value, lj_tl$p.value, lj_tl2$p.value, ar_tl$p.value))
kable(test_result, 
  caption = "Bảng: 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: Kết quả các kiểm định

Test	P_value_VNI	P_value_TL
ADF	0.0100000	0.0100000
J-B	0.0000000	0.0000000
Q(2)	0.9982569	0.9543091
Q(2)^2	0.0000000	0.0000000
ARCH(2)	0.0000000	0.0000000

pearson <- cor(mhnn$VNI,mhnn$TL, method="pearson")
spearman <- cor(mhnn$VNI,mhnn$TL, method="spearman")
kendall <- cor(mhnn$VNI,mhnn$TL, method="kendall")

#Trình bày kết quả
relat <- data.frame('Tương quan' = 'VNI-STI',
                    Pearson = pearson,
                    spearman = spearman,
                    Kendall = kendall)
kable(relat, 
      col.names = c("Phương pháp", "Pearson", "Spearman","Kendall"),
      caption = "Bảng 4: 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 4: Kết quả hệ số tương quan

Phương pháp	Pearson	Spearman	Kendall
VNI-STI	0.3421649	0.2286002	0.1575133

Trực quan hóa hệ số tương quan

# Trực quan hóa hệ số tương quan 
par(mfrow = c(1,2))

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

library(PerformanceAnalytics)
chart.Correlation(mhnn, histogram=TRUE, pch=19)

## Warning in par(usr): argument 1 does not name a graphical parameter

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

par(mfrow = c(1,2))
a11 <- CalculateReturns(data, method = 'log')
ggplot(data, aes(x = Date, y= a11$VNI))+
  geom_line()+
  labs(title = "Hình: 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()`).

a22 <- CalculateReturns(data, method = 'log')

ggplot(data, aes(x = Date, y= a22$THAILAN))+
  geom_line()+
  labs(title = "Hình 10: Biến động theo ngày của TSLN Thái Lan",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()`).

Biểu đồ Scatter tỷ suất lợi nhuận

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

## `geom_smooth()` using formula = 'y ~ x'

Kiểm định cho chuối dữ liệu

# Kiểm định ARMA cho chỉ số VNI
library(forecast)
modeld1<- auto.arima(mhnn$VNI)
modeld1

## Series: mhnn$VNI 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2   mean
##       0.1001  0.0859  5e-04
## s.e.  0.0266  0.0266  3e-04
## 
## sigma^2 = 0.0001177:  log likelihood = 4376.02
## AIC=-8744.05   AICc=-8744.02   BIC=-8723.05

# Kiểm định ARMA cho chỉ số MSIC Thailand
library(forecast)
modeld2<- auto.arima(mhnn$TL)
modeld2

## Series: mhnn$TL 
## ARIMA(5,0,2) with zero mean 
## 
## Coefficients:
##           ar1      ar2     ar3     ar4     ar5     ma1     ma2
##       -1.3129  -0.8625  0.0666  0.0759  0.0684  1.2652  0.8382
## s.e.   0.0842   0.1036  0.0496  0.0481  0.0298  0.0795  0.0910
## 
## sigma^2 = 0.000137:  log likelihood = 4270.89
## AIC=-8525.78   AICc=-8525.67   BIC=-8483.77

#Garch(1,1)
#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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$VNI)

#Setup garch(1,2) and estimate 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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$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 = mhnn$VNI)

#Setup garch(1,1) and estimate for TL coef
#norm
spec_TL_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")

TL_garch11_norm <- ugarchfit(spec = spec_TL_ugarch11_norm, data = mhnn$TL)

#std
spec_TL_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")

TL_garch11_std <- ugarchfit(spec = spec_TL_ugarch11_std, data = mhnn$TL)

#sstd
spec_TL_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")

TL_garch11_sstd <- ugarchfit(spec = spec_TL_ugarch11_sstd, data = mhnn$TL)

#ged
spec_TL_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")

TL_garch11_ged <- ugarchfit(spec = spec_TL_ugarch11_ged, data = mhnn$TL)

#sged
spec_TL_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")

TL_garch11_sged <- ugarchfit(spec = spec_TL_ugarch11_sged, data = mhnn$TL)

#Setup garch(1,2) and estimate for TL coef
#norm
spec_TL_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")

TL_garch12_norm <- ugarchfit(spec = spec_TL_ugarch12_norm, data = mhnn$TL)

#std
spec_TL_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")

TL_garch12_std <- ugarchfit(spec = spec_TL_ugarch12_std, data = mhnn$TL)

#sstd
spec_TL_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")

TL_garch12_sstd <- ugarchfit(spec = spec_TL_ugarch12_sstd, data = mhnn$TL)

#ged
spec_TL_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")

TL_garch12_ged <- ugarchfit(spec = spec_TL_ugarch12_ged, data = mhnn$TL)

#sged
spec_TL_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")

TL_garch12_sged <- ugarchfit(spec = spec_TL_ugarch12_sged, data = mhnn$TL)

#Setup garch(2,1) and estimate for TL coef
#norm
spec_TL_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")

TL_garch21_norm <- ugarchfit(spec = spec_TL_ugarch21_norm, data = mhnn$TL)

#std
spec_TL_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")

TL_garch21_std <- ugarchfit(spec = spec_TL_ugarch21_std, data = mhnn$TL)

#sstd
spec_TL_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")

TL_garch21_sstd <- ugarchfit(spec = spec_TL_ugarch21_sstd, data = mhnn$TL)

#ged
spec_TL_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")

TL_garch21_ged <- ugarchfit(spec = spec_TL_ugarch21_ged, data = mhnn$TL)

#sged
spec_TL_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")

TL_garch21_sged <- ugarchfit(spec = spec_TL_ugarch21_sged, data = mhnn$TL)

#Setup garch(2,2) and estimate for TL coef
#norm
spec_TL_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")

TL_garch22_norm <- ugarchfit(spec = spec_TL_ugarch22_norm, data = mhnn$TL)

#std
spec_TL_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")

TL_garch22_std <- ugarchfit(spec = spec_TL_ugarch22_std, data = mhnn$TL)

#sstd
spec_TL_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")

TL_garch22_sstd <- ugarchfit(spec = spec_TL_ugarch22_sstd, data = mhnn$TL)

#ged
spec_TL_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")

TL_garch22_ged <- ugarchfit(spec = spec_TL_ugarch22_ged, data = mhnn$TL)

#sged
spec_TL_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")

TL_garch22_sged <- ugarchfit(spec = spec_TL_ugarch22_sged, data = mhnn$TL)

#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 sti coef
TL.model.list <- list(garch11n = TL_garch11_norm, garch11t = TL_garch11_std, garch11st = TL_garch11_sstd, 
                       garch11g = TL_garch11_ged, garch11sg = TL_garch11_sged, 
                       garch12n = TL_garch12_norm, garch12t = TL_garch12_std, garch12st = TL_garch12_sstd, 
                       garch12g = TL_garch12_ged, garch12sg = TL_garch12_sged,
                       garch21n = TL_garch21_norm, garch21t = TL_garch21_std, garch21st = TL_garch21_sstd, 
                       garch21g = TL_garch21_ged, garch21sg = TL_garch21_sged,
                       garch22n = TL_garch22_norm, garch22t = TL_garch22_std, garch22st = TL_garch22_sstd, 
                       garch22g = TL_garch22_ged, garch22sg = TL_garch22_sged
                       )

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

rownames(TL.info.mat) <- rownames(infocriteria(TL_garch11_norm))

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

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

Trình bày kết quả dạng mô hình phân phối biên

#Trình bày kết quả
mar_model <- data.frame(
  rate = c("VNI", "TL"),
  Format = c("ARMA(0,0,2)-GJR-GARCH(1,1)-Skewed Student t", "ARMA(5,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(0,0,2)-GJR-GARCH(1,1)-Skewed Student t
TL	ARMA(5,0,2)-GJR-GARCH(1,1)-Student t

Kết quả mô hình ARMA-GJR-GARCH của chỉ số VNI

# Trình bày các 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(TL_garch11_std)

kable(as.data.frame(fit1), 
  caption = "Bảng: Kết quả mô hình ARMA(0,0,2)-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: Kết quả mô hình ARMA(0,0,2)-GJR-Garch(1,1)-Skewed Student của biến VNI

	Estimate	Std. Error	Pr(>|z|)
mu	0.0004793	0.0003330	0.1500560
ar1	0.9368004	0.0446748	0.0000000
ma1	-0.9072217	0.0526405	0.0000000
omega	0.0000044	0.0000012	0.0001559
alpha1	0.0346699	0.0061189	0.0000000
beta1	0.8552800	0.0156780	0.0000000
gamma1	0.1377139	0.0323193	0.0000203
skew	0.8823255	0.0329038	0.0000000
shape	5.2500555	0.6718660	0.0000000

Kết quả mô hình ARMA-GJR-GARCH của chỉ số MSIC Thái Lan

# Trình bày các 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(TL_garch11_std)

kable(as.data.frame(fit2), 
  caption = "Bảng: Kết quả mô hình ARMA(5,0,2)-GJR-Garch(1,1)-Student của biến TL", 
  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: Kết quả mô hình ARMA(5,0,2)-GJR-Garch(1,1)-Student của biến TL

	Estimate	Std. Error	Pr(>|z|)
mu	0.0000508	0.0001929	0.7922457
ma1	-0.0044702	0.0256353	0.8615700
ma2	-0.0054035	0.0255649	0.8326045
omega	0.0000010	0.0000016	0.5286396
alpha1	0.0260499	0.0088047	0.0030899
beta1	0.9246161	0.0203925	0.0000000
gamma1	0.0852088	0.0050326	0.0000000
shape	5.3796013	0.5564032	0.0000000

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 TL
TL.res <- residuals(TL_garch11_std) / sigma(TL_garch11_std)  # Phần dư đã được chuẩn hóa

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

v <- pdist(distribution = "std", q = TL.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_tl <- ad.test(v, "punif")

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

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

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

kable(test, 
      caption = "Bảng: 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.

Bảng: 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.8290458	0.7483153	0.4891057
P_value TL	0.9608250	0.9951223	0.9992834

Trích xuất phần dư chuỗi

# Trích phần dư chuỗi
# Trích xuất phần dư chuỗi VNI
VNI.res <- residuals(VNI_garch11_sstd)/sigma(VNI_garch11_sstd)
fitdist(distribution = "sstd", VNI.res, control = list())

## $pars
##          mu       sigma        skew       shape 
## 0.004364585 1.000913590 0.884420004 5.221757263 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2068.772 1932.819 1932.819
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]       [,2]       [,3]      [,4]
## [1,] 1765.942250  260.09164 -499.98369  4.974985
## [2,]  260.091635 1569.46192  146.76419 46.686143
## [3,] -499.983689  146.76419 1174.26730 10.119976
## [4,]    4.974985   46.68614   10.11998  3.355333
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 115
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.14062 secs
## 
## $vscale
## [1] 1 1 1 1 1

s = pdist("sstd",VNI.res, mu =0.004364585, sigma =1.000913590 , skew =0.884420004, shape =5.221757263 )
head(s,10)

##  [1] 0.8298124 0.6382865 0.5567212 0.9962126 0.7551382 0.8550152 0.1605119
##  [8] 0.4542302 0.1947739 0.1826931

# Trích suất phần dư chuỗi TL
TL.res <- residuals(TL_garch11_std)/sigma(TL_garch11_std)
fitdist(distribution = "std", TL.res, control = list())

## $pars
##            mu         sigma         shape 
## -5.065177e-05  1.001621e+00  5.347285e+00 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 1951.862 1945.716 1945.716
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]        [,2]      [,3]
## [1,] 1794.103404    9.572421  2.241022
## [2,]    9.572421 1749.550633 47.318537
## [3,]    2.241022   47.318537  3.097759
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 112
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.1120169 secs
## 
## $vscale
## [1] 1 1 1 1

d = pdist("std",TL.res, mu =-5.065177e-05 , sigma = 1.001621e+00,  shape =5.347285e+00 )
head(d,10)

##  [1] 0.2095409 0.9485960 0.9386299 0.5360867 0.3214849 0.4160346 0.1273236
##  [8] 0.6414510 0.3163037 0.8632297

Biểu đồ các họ Coupula

library(knitr)

## Warning: package 'knitr' 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

## The following object is masked from 'package:lubridate':
## 
##     interval

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

#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, '#CCFF33')
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', '#CCFF33')

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

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,'#CC6633')
dev.off()

## png 
##   2

clays <- rasterGrob(readPNG('clays.png'), interpolate = TRUE)

png('gums.png')
scatter_plot(random_gum,'#FF6699')
dev.off()

## png 
##   2

gums <- rasterGrob(readPNG('gums.png'), interpolate = TRUE)
clayp <- persp_plot(cop_clay,'clayp.png','#CC6633')

gump <- persp_plot(cop_gum,'gump.png','#FF6699')


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

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))
             )

set.seed(123)
#Mô phỏng copula survival gumbel
cop_surgum <- VC2copula::surGumbelCopula(param = 5)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_surgum <- rCopula(copula = cop_surgum,n = 7600)
#Mô phỏng copula survival clayton  với p=4
cop_surclay <- VC2copula::surClaytonCopula(param = 4)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_surclay <- rCopula(copula = cop_surclay,n = 7600)
png('surclays.png')
scatter_plot(random_surclay,'#9966FF')
dev.off()

## png 
##   2

surclays <- rasterGrob(readPNG('surclays.png'), interpolate = TRUE)

png('surgums.png')
scatter_plot(random_surgum,'#FF3300')
dev.off()

## png 
##   2

surgums <- rasterGrob(readPNG('surgums.png'), interpolate = TRUE)
surclayp <- persp_plot(cop_surclay, "surclayp.png","#9966FF")
surgump <- persp_plot(cop_surgum, "surgump.png","#FF3300")

legend <- legendGrob(labels = c("Survival Clayton","Survival Gumbel"), pch = 15,
                    gp = gpar(col = c('#9966FF','#FF3300'), fill = c('#9966FF','#FF3300' )))

grid.arrange(surclays, surclayp,surgums, surgump, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 3: Biểu đồ phân tán và PDF của Copula Sur 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,'#FFCCFF')
dev.off()

## png 
##   2

franks <- rasterGrob(readPNG('franks.png'), interpolate = TRUE)

png('joes.png')
scatter_plot(random_joe,'#FFCC66')
dev.off()

## png 
##   2

joes <- rasterGrob(readPNG('joes.png'), interpolate = TRUE)
frankp <- persp_plot(cop_frank, "frankp.png","#FFCCFF")
joep <- persp_plot(cop_joe, "joep.png","#FFCC66")

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

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
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-TL'),
                    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 = "Bảng 9: 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.

Bảng 9: 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-TL	Survival Gumbel	1.177716	0	0.1986219	0	0.150899	-117.4945	-112.2438

Phối cảnh PDF theo chuỗi dữ liệu phân

persp(VC2copula::BB7Copula(param = c(est_opt_cop$par, est_opt_cop$par2)),dCopula, 
        xlab = 'u', ylab = 'v', col = 'lightblue', ltheta = 120,  
        ticktype = "detailed", cex.axis = 0.8,  main = 'Hình 13: Phối cảnh PDF của mô hình Copula Gumbel')

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
sample_size <- 500  # Kích thước mẫu
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

Phân phối mật độ xác suất của Copula Gumbel

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 ")
         ))

Đồ thị mật độ Gumbel Copula

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

df <- data.frame(u = s, k = d)

# Vẽ đồ thị mật độ
ggplot(df, aes(x = u, y = k)) +
  geom_point(alpha = 0.5) +
  geom_density_2d() +
  labs(title = "Đồ thị mật độ Gumbel Copula",
       x = "u",
       y = "k") +
  theme_minimal()

# Tải các thư viện cần thiết
library(copula)
library(MASS)

## 
## 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 (thay thế cd bằng dữ liệu thực của bạn)
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 = "turquoise")