TIỂU LUẬN MÔ HÌNH NGẪU NHIÊN

1.ƯỚC LƯỢNG MÔ HÌNH COPULA

1.1. GÁN DỮ LIỆU

library(xlsx) 

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

library(PerformanceAnalytics)

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

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

#Nhập dữ liệu
data <- read.xlsx('C:/Users/admin/OneDrive - UFM/Documents/thnn1/THNN1/vniset.xlsx', sheetIndex = 1, header = T)
# Tính lợi suất logarit hàng ngày của cổ phiếu VNI và NYA
tsVNI <- diff(log(data$VNI), lag = 1)
tsSETI <- diff(log(data$SETI), lag = 1)
# Tạo một data frame mới chứa lợi suất logarit hàng ngày của cả hai cổ phiếu
rt <- data.frame(tsVNI = tsVNI, tsSETI = tsSETI)
summary(data)

##       Date                 VNI              SETI     
##  Min.   :2018-01-03   Min.   : 659.2   Min.   :1024  
##  1st Qu.:2019-07-03   1st Qu.: 962.9   1st Qu.:1528  
##  Median :2020-12-24   Median :1048.2   Median :1609  
##  Mean   :2020-12-30   Mean   :1089.7   Mean   :1573  
##  3rd Qu.:2022-07-04   3rd Qu.:1198.4   3rd Qu.:1654  
##  Max.   :2023-12-28   Max.   :1528.6   Max.   :1837

1.2. MA TRẬN HỆ SỐ

res <- cor(rt)
round <- round(res, 3)
library(corrplot)

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

## corrplot 0.92 loaded

corrplot(res, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

chart.Correlation(rt, histogram=TRUE, pch=19)

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

library(ggplot2)
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(ggcorrplot)

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

df <- dplyr::select_if(rt, is.numeric)
r <- cor(df, use="complete.obs")
ggcorrplot(r)

Biến động của chuỗi tỷ suất lợi nhuận

plot.ts(rt$tsVNI)

Biến động của chuỗi tỷ suất lợi nhuận SETI

plot.ts(rt$tsSETI)

library(tidyr)
pivot_longer(rt, 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 7: Biểu đồ hộp của TSLN VNI và SETI",
       x = "Chỉ số",
       y = "Tỷ suất sinh lợi")

2. KIỂM ĐỊNH DỮ LIỆU

2.1. KIỂM ĐỊNH TÍNH DỪNG

library(tseries)

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

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

#Kiểm định tính dừng của chuỗi VNI
adf_vni <- adf.test(rt$tsVNI)

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

print(adf_vni)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  rt$tsVNI
## Dickey-Fuller = -11.512, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

#Kiểm định tính dừng của chuổi SETI
adf_seti <- adf.test(rt$tsSETI)

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

print(adf_seti)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  rt$tsSETI
## Dickey-Fuller = -8.9716, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

2.2. KIỂM ĐỊNH PHÂN PHỐI CHUẨN

# kiểm định phân phối chuẩn của VNI
library(tseries)
jq_vni <- jarque.bera.test(rt$tsVNI)
print(jq_vni)

## 
##  Jarque Bera Test
## 
## data:  rt$tsVNI
## X-squared = 941.92, df = 2, p-value < 2.2e-16

# kiểm định phân phối chuẩn của SETI
jq_seti <- jarque.bera.test(rt$tsSETI)
print(jq_seti)

## 
##  Jarque Bera Test
## 
## data:  rt$tsSETI
## X-squared = 36162, df = 2, p-value < 2.2e-16

3. MÔ HÌNH ARIMA

library(FinTS)

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

# Ước lượng và trích xuất phần dư từ mô hình ARMA tối ưu
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

arima_VNI <- autoarfima(rt$tsVNI,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

arima_seti <- autoarfima(rt$tsSETI,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

re_VNI <- arima_VNI$fit@fit$residuals
re_seti <- arima_seti$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_seti <- Box.test(re_seti,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_seti2 <- Box.test(re_seti^2,type = 'Ljung-Box', lag = 2)

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

#Trình bày kết quả
library(knitr)
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)
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_STI = c(adf_seti$p.value, jq_seti$p.value, lj_seti$p.value, lj_seti2$p.value, ar_seti$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_STI
ADF	0.0100000	0.0100000
J-B	0.0000000	0.0000000
Q(2)	0.9996102	0.9168313
Q(2)^2	0.0000000	0.0000000
ARCH(2)	0.0000000	0.0000000

chart.Correlation(rt, histogram=TRUE, pch=19)

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

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

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

4. MÔ HÌNH GARCH

4.1. SETI

# GJR-GARCH(11)SETI
seti.garch11n.spec <- ugarchspec(variance.model = list(model="gjrGARCH",garchOrder=c(1,1)),mean.model = list(armaOrder=c(1,0),include.mean=TRUE),distribution.model = "norm")
seti.garch11n.fit <- ugarchfit(spec = seti.garch11n.spec,  tsSETI)
seti.garch11t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "std")
 seti.garch11t.fit <- ugarchfit(spec =seti.garch11t.spec, tsSETI)
seti.garch11st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sstd")
 seti.garch11st.fit <- ugarchfit(spec = seti.garch11st.spec, tsSETI)
seti.garch11g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "ged")
seti.garch11g.fit <- ugarchfit(spec = seti.garch11g.spec, tsSETI)
seti.garch11sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sged")
seti.garch11sg.fit <- ugarchfit(spec = seti.garch11sg.spec, tsSETI)

#  GJR-GARCH(12)SETI
seti.garch12n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "norm")
seti.garch12n.fit <- ugarchfit(spec = seti.garch12n.spec, tsSETI)
seti.garch12t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "std")
seti.garch12t.fit <- ugarchfit(spec = seti.garch12t.spec, tsSETI)
seti.garch12st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sstd")
seti.garch12st.fit <- ugarchfit(spec = seti.garch12st.spec, tsSETI)
seti.garch12g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "ged")
seti.garch12g.fit <- ugarchfit(spec = seti.garch12g.spec, tsSETI)
seti.garch12sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sged")
seti.garch12sg.fit <- ugarchfit(spec = seti.garch12sg.spec, tsSETI)
# GJR-GARCH(21)SETI
seti.garch21n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "norm")
seti.garch21n.fit <- ugarchfit(spec = seti.garch21n.spec, tsSETI)
seti.garch21t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "std")
seti.garch21t.fit <- ugarchfit(spec = seti.garch21t.spec, tsSETI)
seti.garch21st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sstd")
seti.garch21st.fit <- ugarchfit(spec = seti.garch21st.spec, tsSETI)
seti.garch21g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "ged")
seti.garch21g.fit <- ugarchfit(spec = seti.garch21g.spec, tsSETI)
seti.garch21sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sged")
seti.garch21sg.fit <- ugarchfit(spec = seti.garch21sg.spec, tsSETI)
# GJR-GARCH(22)SETI
seti.garch22n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "norm")
seti.garch22n.fit <- ugarchfit(spec = seti.garch22n.spec, tsSETI)
seti.garch22t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "std")
seti.garch22t.fit <- ugarchfit(spec = seti.garch22t.spec, tsSETI)
seti.garch22st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sstd")
seti.garch22st.fit <- ugarchfit(spec = seti.garch22st.spec, tsSETI)
seti.garch22g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "ged")
seti.garch22g.fit <- ugarchfit(spec = seti.garch22g.spec, tsSETI)
seti.garch22sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sged")
seti.garch22sg.fit <- ugarchfit(spec = seti.garch22sg.spec, tsSETI)

4.2. VNi

# GJR-GARCH(11)VNI
vn.garch11n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "norm")
vn.garch11n.fit <- ugarchfit(spec = vn.garch11n.spec, tsVNI)
vn.garch11t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "std")
vn.garch11t.fit <- ugarchfit(spec = vn.garch11t.spec, tsVNI)
vn.garch11st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sstd")
vn.garch11st.fit <- ugarchfit(spec = vn.garch11st.spec, tsVNI)
vn.garch11g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "ged")
vn.garch11g.fit <- ugarchfit(spec = vn.garch11g.spec, tsVNI)
vn.garch11sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sged")
vn.garch11sg.fit <- ugarchfit(spec = vn.garch11sg.spec, tsVNI)

# GJR-GARCH(12)VNI
vn.garch12n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "norm")
vn.garch12n.fit <- ugarchfit(spec = vn.garch12n.spec, tsVNI)
vn.garch12t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "std")
vn.garch12t.fit <- ugarchfit(spec = vn.garch12t.spec, tsVNI)
vn.garch12st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sstd")
vn.garch12st.fit <- ugarchfit(spec = vn.garch12st.spec, tsVNI)
vn.garch12g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "ged")
vn.garch12g.fit <- ugarchfit(spec = vn.garch12g.spec, tsVNI)
vn.garch12sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sged")
vn.garch12sg.fit <- ugarchfit(spec = vn.garch12sg.spec, tsVNI)

# GJR-GARCH(21)VNI
vn.garch21n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "norm")
vn.garch21n.fit <- ugarchfit(spec = vn.garch21n.spec, tsVNI)
vn.garch21t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "std")
vn.garch21t.fit <- ugarchfit(spec = vn.garch21t.spec, tsVNI)
vn.garch21st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sstd")
vn.garch21st.fit <- ugarchfit(spec = vn.garch21st.spec, tsVNI)

vn.garch21g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "ged")
vn.garch21g.fit <- ugarchfit(spec = vn.garch21g.spec, tsVNI)
vn.garch21sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sged")
vn.garch21sg.fit <- ugarchfit(spec = vn.garch21sg.spec, tsVNI)

# GJR-GARCH(22)VNI
vn.garch22n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "norm")
vn.garch22n.fit <- ugarchfit(spec = vn.garch22n.spec, tsVNI)
vn.garch22t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "std")
vn.garch22t.fit <- ugarchfit(spec = vn.garch22t.spec, tsVNI)
vn.garch22st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sstd")
vn.garch22st.fit <- ugarchfit(spec = vn.garch22st.spec, tsVNI)
vn.garch22g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model= list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "ged")
vn.garch22g.fit <- ugarchfit(spec = vn.garch22g.spec, tsVNI)
vn.garch22sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),mean.model = list(armaOrder = c(2, 2), include.mean = TRUE), distribution.model = "sged")
vn.garch22sg.fit <- ugarchfit(spec = vn.garch22sg.spec, tsVNI)

5. LỰA CHỌN MÔ HÌNH BIÊN PHÙ HỢP

5.1. LỰA CHỌN MÔ HÌNH BIÊN PHÙ HỢP NHẤT CHO CHUỖI VNI

vn.model.list <- list(garch11n = vn.garch11n.fit, garch11t = vn.garch11t.fit, garch11st = vn.garch11st.fit,garch11g = vn.garch11g.fit, garch11sg = vn.garch11sg.fit, garch12n = vn.garch12n.fit, garch12t = vn.garch12t.fit, garch12st = vn.garch12st.fit, garch12g =vn.garch12g.fit, garch12sg = vn.garch12sg.fit, garch21n = vn.garch21n.fit, garch21t = vn.garch21t.fit, garch21st = vn.garch21st.fit, garch21g =vn.garch21g.fit, garch21sg = vn.garch21sg.fit, garch22n = vn.garch22n.fit, garch22t = vn.garch22t.fit, garch22st = vn.garch22st.fit, garch22g =vn.garch22g.fit, garch22sg = vn.garch22sg.fit)
vn.info.mat <- sapply(vn.model.list, infocriteria)
rownames(vn.info.mat)<-rownames(infocriteria(vn.garch11n.fit))
vn.info.mat

##               garch11n  garch11t garch11st  garch11g garch11sg  garch12n
## Akaike       -6.031198 -6.169380 -6.188726 -6.162834 -6.180806 -6.026379
## Bayes        -5.997582 -6.132029 -6.147640 -6.125483 -6.139720 -5.989028
## Shibata      -6.031279 -6.169481 -6.188847 -6.162934 -6.180927 -6.026479
## Hannan-Quinn -6.018634 -6.155420 -6.173370 -6.148874 -6.165450 -6.012419
##               garch12t garch12st  garch12g garch12sg  garch21n  garch21t
## Akaike       -6.168730 -6.188242 -6.162125 -6.182412 -6.030428 -6.170462
## Bayes        -6.127644 -6.143421 -6.121039 -6.137591 -5.989342 -6.125641
## Shibata      -6.168852 -6.188387 -6.162246 -6.182556 -6.030550 -6.170607
## Hannan-Quinn -6.153374 -6.171490 -6.146769 -6.165660 -6.015072 -6.153710
##              garch21st  garch21g garch21sg  garch22n  garch22t garch22st
## Akaike       -6.191820 -6.170359 -6.179975 -6.029005 -6.169039 -6.190396
## Bayes        -6.143263 -6.125538 -6.131418 -5.984184 -6.120483 -6.138105
## Shibata      -6.191989 -6.170503 -6.180144 -6.029149 -6.169208 -6.190592
## Hannan-Quinn -6.173671 -6.153607 -6.161826 -6.012253 -6.150891 -6.170852
##               garch22g garch22sg
## Akaike       -6.162231 -6.178551
## Bayes        -6.113675 -6.126260
## Shibata      -6.162400 -6.178747
## Hannan-Quinn -6.144083 -6.159007

vn.inds <- which(vn.info.mat == min(vn.info.mat), arr.ind=TRUE)
model.vn <- colnames(vn.info.mat)[vn.inds[,2]]
model.vn

## [1] "garch21st"

5.2. LỰA CHỌN MÔ HÌNH BIÊN PHÙ HỢP NHẤT CHO CHUỖI SETI

seti.model.list <- list(garch11n = seti.garch11n.fit, garch11t = seti.garch11t.fit,garch11st=seti.garch11st.fit, garch11g = seti.garch11g.fit, garch11sg = seti.garch11sg.fit,garch12n = seti.garch12n.fit, garch12t = seti.garch12t.fit, garch12st = seti.garch12st.fit,garch12g = seti.garch12g.fit, garch12sg = seti.garch12sg.fit,garch21n = seti.garch21n.fit, garch21t = seti.garch21t.fit, garch21st = seti.garch21st.fit,garch21g = seti.garch21g.fit, garch21sg = seti.garch21sg.fit, garch22n = seti.garch22n.fit, garch22t =seti.garch22t.fit, garch22st = seti.garch22st.fit, garch22g = seti.garch22g.fit, garch22sg = seti.garch22sg.fit)
seti.info.mat <- sapply(seti.model.list, infocriteria)

rownames(seti.info.mat) <- rownames(infocriteria(seti.garch11n.fit))
seti.info.mat

##               garch11n  garch11t garch11st  garch11g garch11sg  garch12n
## Akaike       -6.754805 -6.817557 -6.825277 -6.807757 -6.819077 -6.752890
## Bayes        -6.732395 -6.791411 -6.795397 -6.781611 -6.789197 -6.726744
## Shibata      -6.754842 -6.817606 -6.825342 -6.807806 -6.819142 -6.752939
## Hannan-Quinn -6.746429 -6.807785 -6.814109 -6.797985 -6.807909 -6.743118
##               garch12t garch12st  garch12g garch12sg  garch21n  garch21t
## Akaike       -6.815717 -6.823440 -6.805946 -6.817285 -6.753050 -6.815654
## Bayes        -6.785836 -6.789824 -6.776065 -6.783669 -6.723169 -6.782038
## Shibata      -6.815781 -6.823521 -6.806011 -6.817366 -6.753115 -6.815735
## Hannan-Quinn -6.804549 -6.810876 -6.794778 -6.804721 -6.741882 -6.803090
##              garch21st  garch21g garch21sg  garch22n  garch22t garch22st
## Akaike       -6.823044 -6.805745 -6.816632 -6.751627 -6.814339 -6.821684
## Bayes        -6.785693 -6.772129 -6.779282 -6.718012 -6.776988 -6.780598
## Shibata      -6.823144 -6.805826 -6.816733 -6.751709 -6.814439 -6.821806
## Hannan-Quinn -6.809083 -6.793181 -6.802672 -6.739063 -6.800379 -6.806328
##               garch22g garch22sg
## Akaike       -6.804368 -6.815222
## Bayes        -6.767017 -6.774136
## Shibata      -6.804468 -6.815343
## Hannan-Quinn -6.790408 -6.799865

seti.inds <- which(seti.info.mat == min(seti.info.mat), arr.ind = TRUE)
model.seti <- colnames(seti.info.mat)[seti.inds[, 2]]
model.seti

## [1] "garch11st"

6. THAM SỐ ƯỚC LƯỢNG MÔ HÌNH BIÊN PHÙ HỢP NHẤT:

vn.garch21st.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(2,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000306    0.000276   1.110100 0.266956
## ar1    -0.965442    0.039883 -24.206619 0.000000
## ar2    -0.884502    0.080041 -11.050586 0.000000
## ma1     0.998065    0.031621  31.563440 0.000000
## ma2     0.916532    0.064962  14.108768 0.000000
## omega   0.000006    0.000001   5.070666 0.000000
## alpha1  0.000000    0.034772   0.000006 0.999995
## alpha2  0.043830    0.038462   1.139555 0.254472
## beta1   0.870317    0.016139  53.927637 0.000000
## gamma1  0.367808    0.101054   3.639718 0.000273
## gamma2 -0.253131    0.092194  -2.745621 0.006040
## skew    0.825169    0.031515  26.183465 0.000000
## shape   4.064768    0.414307   9.811014 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000306    0.000310   0.986309 0.323982
## ar1    -0.965442    0.029769 -32.430823 0.000000
## ar2    -0.884502    0.122135  -7.242022 0.000000
## ma1     0.998065    0.022092  45.178360 0.000000
## ma2     0.916532    0.098481   9.306641 0.000000
## omega   0.000006    0.000002   3.381772 0.000720
## alpha1  0.000000    0.031700   0.000007 0.999994
## alpha2  0.043830    0.038150   1.148872 0.250609
## beta1   0.870317    0.013005  66.919206 0.000000
## gamma1  0.367808    0.085493   4.302218 0.000017
## gamma2 -0.253131    0.087239  -2.901583 0.003713
## skew    0.825169    0.033205  24.850473 0.000000
## shape   4.064768    0.395831  10.268954 0.000000
## 
## LogLikelihood : 4362.753 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.1918
## Bayes        -6.1433
## Shibata      -6.1920
## Hannan-Quinn -6.1737
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       3.240 0.0718632
## Lag[2*(p+q)+(p+q)-1][11]     8.209 0.0004455
## Lag[4*(p+q)+(p+q)-1][19]    12.091 0.1839062
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       2.075  0.1498
## Lag[2*(p+q)+(p+q)-1][8]      2.968  0.6914
## Lag[4*(p+q)+(p+q)-1][14]     5.126  0.7537
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]   0.06077 0.500 2.000  0.8053
## ARCH Lag[6]   0.11087 1.461 1.711  0.9864
## ARCH Lag[8]   0.62063 2.368 1.583  0.9712
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  6.5797
## Individual Statistics:              
## mu     0.26388
## ar1    0.06189
## ar2    0.19050
## ma1    0.08022
## ma2    0.16516
## omega  0.67789
## alpha1 0.06682
## alpha2 0.05680
## beta1  0.07981
## gamma1 0.07785
## gamma2 0.07532
## skew   0.13648
## shape  0.08689
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.89 3.15 3.69
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias            1.324 0.18584    
## Negative Sign Bias   1.720 0.08557   *
## Positive Sign Bias   1.110 0.26705    
## Joint Effect         6.508 0.08935   *
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     39.28     0.004057
## 2    30     55.18     0.002362
## 3    40     58.43     0.023464
## 4    50     69.06     0.030978
## 
## 
## Elapsed time : 1.711117

seti.garch11st.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu     -0.000366    0.000217 -1.68456 0.092074
## ar1     0.024742    0.027082  0.91359 0.360934
## omega   0.000003    0.000004  0.74441 0.456631
## alpha1  0.021362    0.022226  0.96111 0.336496
## beta1   0.882705    0.023965 36.83272 0.000000
## gamma1  0.119083    0.029759  4.00153 0.000063
## skew    0.875611    0.033034 26.50642 0.000000
## shape   5.903650    0.932408  6.33162 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu     -0.000366    0.000495 -0.73940 0.459664
## ar1     0.024742    0.048525  0.50989 0.610132
## omega   0.000003    0.000027  0.10686 0.914899
## alpha1  0.021362    0.129252  0.16527 0.868729
## beta1   0.882705    0.128528  6.86783 0.000000
## gamma1  0.119083    0.093003  1.28042 0.200398
## skew    0.875611    0.056889 15.39164 0.000000
## shape   5.903650    1.702178  3.46829 0.000524
## 
## LogLikelihood : 4802.757 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.8253
## Bayes        -6.7954
## Shibata      -6.8253
## Hannan-Quinn -6.8141
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.309  0.2526
## Lag[2*(p+q)+(p+q)-1][2]     1.374  0.5097
## Lag[4*(p+q)+(p+q)-1][5]     2.790  0.4860
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1137  0.7359
## Lag[2*(p+q)+(p+q)-1][5]    3.9049  0.2662
## Lag[4*(p+q)+(p+q)-1][9]    5.3572  0.3789
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     4.737 0.500 2.000 0.02952
## ARCH Lag[5]     5.697 1.440 1.667 0.07074
## ARCH Lag[7]     5.887 2.315 1.543 0.14900
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  7.5132
## Individual Statistics:              
## mu     0.09630
## ar1    0.04243
## omega  0.85619
## alpha1 0.27775
## beta1  0.14873
## gamma1 0.06361
## skew   0.06598
## shape  0.04847
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.3769 0.1687    
## Negative Sign Bias  0.9885 0.3231    
## Positive Sign Bias  0.3875 0.6985    
## Joint Effect        2.1717 0.5376    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     11.24       0.9154
## 2    30     25.41       0.6567
## 3    40     39.19       0.4616
## 4    50     48.13       0.5083
## 
## 
## Elapsed time : 0.4570529

7. KIỂM ĐỊNH SỰ PHÙ HỢP CỦA MÔ HÌNH BIÊN

seti.res <- residuals(seti.garch11st.fit)/sigma(seti.garch11st.fit)
fitdist(distribution = "sstd", seti.res, control = list())

## $pars
##         mu      sigma       skew      shape 
## 0.01189854 0.99227744 0.88124855 6.00644104 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2068.117 1932.468 1932.468
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]       [,2]        [,3]      [,4]
## [1,] 1761.381361  323.99940 -481.508214  5.008668
## [2,]  323.999399 2011.57727  130.597322 38.268395
## [3,] -481.508214  130.59732 1195.814434  3.645681
## [4,]    5.008668   38.26839    3.645681  1.953110
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 114
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.07715392 secs
## 
## $vscale
## [1] 1 1 1 1 1

u <- pdist(distribution = "sstd", q = seti.res, mu = -0.03294783 , sigma = 1.00098561 , skew =  0.88663219 , shape = 6.27208051 )
## Trích xuất chuỗi phần dư v của chuỗi lợi suất VN
vn.res <- residuals(vn.garch21st.fit)/sigma(vn.garch21st.fit)
fitdist(distribution = "sstd", vn.res, control = list())

## $pars
##           mu        sigma         skew        shape 
## -0.005743176  1.006734640  0.822052716  4.022462887 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2041.359 1861.450 1861.450
## 
## $lagrange
## [1] 0
## 
## $hessian
##           [,1]      [,2]       [,3]      [,4]
## [1,] 2124.8203  538.3958 -697.97625  30.42440
## [2,]  538.3958 1708.6974  205.35459 129.16418
## [3,] -697.9762  205.3546 1323.90914  27.78374
## [4,]   30.4244  129.1642   27.78374  14.58298
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 111
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.07035494 secs
## 
## $vscale
## [1] 1 1 1 1 1

v <- pdist("sstd",vn.res, mu = 0.01683462, sigma = 0.99427949, skew = 0.83629797, shape =  3.99205291)

8. CÁC KIỂM ĐỊNH SỰ PHÙ HỢP CỦA MÔ HÌNH BIÊN

# Kiem dinh Anderson-Darling
library(nortest)
ad.test(v)

## 
##  Anderson-Darling normality test
## 
## data:  v
## A = 14.959, p-value < 2.2e-16

ad.test(v)

## 
##  Anderson-Darling normality test
## 
## data:  v
## A = 14.959, p-value < 2.2e-16

9. KIỂM ĐỊNH CRAMER-VON MISES

# Null hypothesis: uniform distribution
cvm.test(u)

## Warning in cvm.test(u): p-value is smaller than 7.37e-10, cannot be computed
## more accurately

## 
##  Cramer-von Mises normality test
## 
## data:  u
## W = 1.9606, p-value = 7.37e-10

cvm.test(v)

## Warning in cvm.test(v): p-value is smaller than 7.37e-10, cannot be computed
## more accurately

## 
##  Cramer-von Mises normality test
## 
## data:  v
## W = 1.9153, p-value = 7.37e-10

10. KIỂM ĐỊNH KS-TEST

#Null hypothesis: uniform distribution
ks.test(u, v)

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  u and v
## D = 0.059786, p-value = 0.01318
## alternative hypothesis: two-sided

ks.test(v, "punif")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  v
## D = 0.033662, p-value = 0.08282
## alternative hypothesis: two-sided

11. ƯỚC LƯỢNG THAM SỐ MÔ HÌNH COPULA

library(VineCopula)

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

BiCopSelect(u, v, familyset= NA, selectioncrit="AIC",indeptest = FALSE, level = 0.05) 

## Bivariate copula: Survival Gumbel (par = 1.17, tau = 0.15)

Stu <- BiCopEst(u, v, family = 2, method = "mle", se = T, max.df = 10)

summary(Stu)

## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  0.23  (SE = 0.03)
## par2: 10  (SE = NA)
## Dependence measures
## -------------------
## Kendall's tau:    0.15 (empirical = 0.14, p value < 0.01)
## Upper TD:         0.02 
## Lower TD:         0.02 
## 
## Fit statistics
## --------------
## logLik:  44.02 
## AIC:    -84.03 
## BIC:    -73.53

12. TÍNH TOÁN THAM SỐ ĐỂ TẠO BIỂU ĐỒ COPULA SURVIVAL JOE

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

# Thiết lập Joe copula với tham số par = 1.07
theta <- 1.07
joe_cop <- joeCopula(param = theta)
# Kiểm tra Kendall's tau của Joe copula
tau <- tau(joe_cop)
print(tau)

## [1] 0.03879337

tnt <- cbind(u, v)
fit <- fitCopula(joeCopula(), data = tnt, method = "ml")  

## Warning in fitCopula.ml(copula, u = data, method = method, start = start, :
## possible convergence problem: optim() gave code=52

# Maximum Likelihood Estimation (MLE)
theta <- coef(fit)
print(theta)

##    alpha 
## 1.121689

joe_cop <- joeCopula(param = theta)
tau <- tau(joe_cop)
print(tau)

##      alpha 
## 0.06532958

par.joeh<-BiCopTau2Par(6, .04, check.taus = TRUE)
par.joel<-BiCopTau2Par(36, -.04, check.taus = TRUE)
obj.joeh <- BiCop(family = 6, par = par.joeh)
obj.joel <- BiCop(family = 36, par = par.joel)
sim.joeh<-BiCopSim(1409,obj.joeh)
sim.joel<-BiCopSim(1409,obj.joel)
plot(sim.joeh, xlab=expression("u"[1]),ylab=expression("u"[2]),xlim=c(0,1),ylim=c(0,1))

13. ĐỒ THỊ PHÂN TÁN

# Vẽ đồ thị phân tán
plot(u, main = "Đồ thị phân tán của biến VNI")

plot(v, main = "Đồ thị phân tán của biến SETI")

14. ĐỒ THỊ MẬT ĐỘ

# Tính toán hàm mật độ và vẽ đồ thị
library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

kde2d_plot <- kde2d(tnt[,1], tnt[,2], n = 50)
persp(kde2d_plot, phi = 50, theta = coef(fit), main = "Đồ thị mật độ của copula JOE")

15. BẢN ĐỒ NHIỆT KERNEL

# Vẽ bản đồ nhiệt của hàm mật độ kernel
image(kde2d_plot, main = "Bản đồ nhiệt của hàm mật độ kernel")
# Vẽ các đường đồng mức của hàm mật độ kernel
contour(kde2d_plot, add = TRUE)

17. ƯỚC LƯỢNG SO SÁNH VỚI CÁC MÔ HÌNH COPULA KHÁC

# Ước lượng và so sánh với một số mô hình copula khác
fit_clayton <- fitCopula(claytonCopula(), data = tnt, method = "ml")
fit_gumbel <- fitCopula(gumbelCopula(), data = tnt, method = "ml")
# In các tham số ước lượng
print(coef(fit_clayton))

##     alpha 
## 0.3341853

print(coef(fit_gumbel))

##    alpha 
## 1.132574

# Load necessary libraries
library(copula)
library(MASS)

# Define the parameter for the Joe copula
theta <- 1.07  # Tham số của bạn

# Create the Joe copula object
joe_copula <- joeCopula(param = theta)

# Generate a sample from the Joe copula
set.seed(123)
n <- 1409
sample <- rCopula(n, joe_copula)

# Perform Kernel Density Estimation
kde2d_plot <- kde2d(tnt[,1], tnt[,2], n = 50)

# Plot the density
persp(kde2d_plot, phi = 50, theta = 50, col = "lightblue",
      xlab = "U", ylab = "V", zlab = "Density", 
      main = "Density Plot of Joe Copula")

# Load necessary libraries
library(copula)
library(MASS)
library(scatterplot3d)

# Define the parameter for the Joe copula
theta <- 1.07

# Create the Joe copula object
joe_copula <- joeCopula(param = theta)


# Perform Kernel Density Estimation
kde2d_result <- kde2d(tnt[,1], tnt[,2], n = 50)

# Convert kde2d result to a format suitable for scatterplot3d
x <- kde2d_result$x
y <- kde2d_result$y
z <- kde2d_result$z

# Create a grid of points for scatterplot3d
grid <- expand.grid(x = x, y = y)

# Flatten z for scatterplot3d
z_flat <- as.vector(z)

# Plot the density using scatterplot3d
scatterplot3d(grid$x, grid$y, z_flat, type = "h", angle = 30, 
              main = "Density Plot of Joe Copula",
              xlab = "Rainfall Depth", 
              ylab = "Shot Duration", 
              zlab = "Bivariate Distribution",
              pch = 20, color = "turquoise")

# Load necessary libraries
library(copula)
library(MASS)

# Define the parameter for the Joe copula
theta <- 1.07

# Create the Joe copula object
joe_copula <- joeCopula(param = theta)

# Perform Kernel Density Estimation
kde2d_result <- kde2d(tnt[,1], tnt[,2], n = 1409)

# Extract the results
x <- kde2d_result$x
y <- kde2d_result$y
z <- kde2d_result$z

# Plot the PDF using persp
persp(x, y, z, phi = 30, theta = 1.07, col = "red", 
      xlab = "Discharge", ylab = "Duration", zlab = "Density", 
      main = "PDF of Joe Copula")

# Plot the contour plot
contour(x, y, z, xlab = "Discharge", ylab = "Duration", 
        main = "Contour Plot of Joe Copula")

# Load the required libraries
library(copula)
library(ggplot2)
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

# Define the Joe copula
joe_copula <- joeCopula(param = 1.07)

# Create histograms of marginal distributions
p_hist_X1 <- ggplot(tnt, aes(x = tsVNI)) +
  geom_histogram(bins = 30, fill = "skyblue", color = "black") +
  ggtitle("Phân Phối Marginal của X1") +
  xlab("X1") +
  ylab("Số lượng") +
  theme_minimal()

p_hist_X2 <- ggplot(tnt, aes(x = tsSETI)) +
  geom_histogram(bins = 30, fill = "salmon", color = "black") +
  ggtitle("Phân Phối Marginal của X2") +
  xlab("X2") +
  ylab("Số lượng") +
  theme_minimal()

# Arrange histograms side by side
grid.arrange(p_hist_X1, p_hist_X2, ncol = 2)

library(copula)
library(VineCopula)
a <- BiCopSelect(u,v, familyset = c(1:10), selectioncrit = "AIC", indeptest = FALSE, level = 0.05)
sur.joe <- joeCopula(param= 1.07)
fit <- fitCopula(sur.joe, data= as.matrix(cbind(u,v)),method="ml")

## Warning in fitCopula.ml(copula, u = data, method = method, start = start, :
## possible convergence problem: optim() gave code=52

coef(fit)

##    alpha 
## 1.121689

a$par

## [1] 1.170376

a$par2

## [1] 0

persp(joeCopula(param = a$par), dCopula, col= "pink", border= "black")