Untitled

CẤU TRÚC PHỤ THUỘC GIỮA THỊ TRƯỜNG CHỨNG KHOÁN VIỆT NAM VÀ THỊ TRƯỜNG CHỨNG KHOÁN MỸ: TIẾP CẬN BẰNG MÔ HÌNH COPULA

library(moments)
library(kableExtra)

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

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:kableExtra':
## 
##     group_rows

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

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

library(ggplot2)

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

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

library(VineCopula)

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

## 
## Attaching package: 'VineCopula'

## The following objects are masked from 'package:VC2copula':
## 
##     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

## The following object is masked from 'package:copula':
## 
##     pobs

library(gridGraphics)

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

## Loading required package: grid

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

## 
## ######################### 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: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

## 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 objects are masked from 'package:moments':
## 
##     kurtosis, skewness

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

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

library(ggcorrplot)

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

library(goftest)
library(VC2copula)
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

Các họ hàm cơ bản của Copula

Copula họ Elip

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")
  } 
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)
}
set.seed(123)
cop_nor <- normalCopula(param = 0.8, dim = 2)
random_nor <- rCopula(copula = cop_nor,n = 7600)
cop_std <- tCopula(param = 0.8, dim = 2, df = 1)
random_std <- rCopula(copula = cop_std,n = 7600)
#Vẽ biểu đồ 
png('nors.png') 
scatter_plot(random_nor, '#ADD8E6')
dev.off()

## png 
##   2

nors <- rasterGrob(readPNG('nors.png'), interpolate = TRUE) 
png('stds.png')
scatter_plot(random_std, '#F08080')
dev.off()

## png 
##   2

stds <- rasterGrob(readPNG('stds.png'), interpolate = TRUE)
nor_per <- persp_plot(cop_nor, 'norp.png', '#ADD8E6')
std_per <- persp_plot(cop_std, 'stdp.png', '#F08080')
legend <- legendGrob(
  labels = c("Gauss", "Student"), pch = 15,
  gp = gpar(col = c('#ADD8E6', '#F08080'), fill = c('#ADD8E6', '#F08080'))
)
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 : Biểu đồ phân tán và phối cảnh PDF của Copula họ Elip", 
                 gp = gpar(fontsize = 15, font = 2))
             )

## Copula họ Archimedean

set.seed(123)
cop_clay <- claytonCopula(param = 4, dim = 2)
random_clay <- rCopula(copula = cop_clay,n = 7600)
cop_gum <- gumbelCopula(param = 5, dim = 2)
random_gum <- rCopula(copula = cop_gum,n = 7600)
png('clays.png')
scatter_plot(random_clay,'#90EE90')
dev.off()

## png 
##   2

clays <- rasterGrob(readPNG('clays.png'), interpolate = TRUE)
png('gums.png')
scatter_plot(random_gum,'#FFA07A')
dev.off()

## png 
##   2

gums <- rasterGrob(readPNG('gums.png'), interpolate = TRUE)
clayp <- persp_plot(cop_clay,'clayp.png','#90EE90')
gump <- persp_plot(cop_gum,'gump.png','#FFA07A')
legend <- legendGrob(
  labels = c("Clayton", "Gumbel"), pch = 15,
  gp = gpar(col = c('#90EE90', '#FFA07A'), fill = c('#90EE90', '#FFA07A'))
)
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 : Biểu đồ phân tán và PDF của Copula Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
cop_surgum <- VC2copula::surGumbelCopula(param = 5)
random_surgum <- rCopula(copula = cop_surgum,n = 7600)
cop_surclay <- VC2copula::surClaytonCopula(param = 4)
random_surclay <- rCopula(copula = cop_surclay,n = 7600)
png('surclays.png')
scatter_plot(random_surclay,'#FFB6C1')
dev.off()

## png 
##   2

surclays <- rasterGrob(readPNG('surclays.png'), interpolate = TRUE)
png('surgums.png')
scatter_plot(random_surgum,'#E6E6FA')
dev.off()

## png 
##   2

surgums <- rasterGrob(readPNG('surgums.png'), interpolate = TRUE)
surclayp <- persp_plot(cop_surclay, "surclayp.png","#FFB6C1")
surgump <- persp_plot(cop_surgum, "surgump.png","#E6E6FA")
legend <- legendGrob(labels = c("Survival Clayton","Survival Gumbel"), pch = 15,
                    gp = gpar(col = c('#FFB6C1','#E6E6FA'), fill = c('#FFB6C1','#E6E6FA' )))
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 : Biểu đồ phân tán và PDF của Copula Sur Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
cop_frank <- frankCopula(param = 9.2)
random_frank <- rCopula(copula = cop_frank,n = 7600)
cop_joe <- joeCopula(param = 3)
random_joe <- rCopula(copula = cop_joe,n = 7600)
png('franks.png')
scatter_plot(random_frank,'#FFFACD')
dev.off()

## png 
##   2

franks <- rasterGrob(readPNG('franks.png'), interpolate = TRUE)
png('joes.png')
scatter_plot(random_joe,'#FF6347')
dev.off()

## png 
##   2

joes <- rasterGrob(readPNG('joes.png'), interpolate = TRUE)
frankp <- persp_plot(cop_frank, "frankp.png","#FFFACD")
joep <- persp_plot(cop_joe, "joep.png","#FF6347")
legend <- legendGrob(labels = c("Franl","Joe"), pch = 15,
                    gp = gpar(col = c('#FFFACD','#FF6347'), fill = c('#FFFACD','#FF6347')))
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 : Biểu đồ phân tán và  PDF của Copula Frank và Joe", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
cop_bb1 <- VC2copula::BB1Copula(param = c(2,1.7))
random_bb1 <- rCopula(copula = cop_bb1,n = 7600)
cop_bb6 <- VC2copula::BB6Copula(param = c(2,4.5))
random_bb6 <- rCopula(copula = cop_bb6,n = 7600)
png('bb1s.png')
scatter_plot(random_bb1,'#6A5ACD')
dev.off()

## png 
##   2

bb1s <- rasterGrob(readPNG('bb1s.png'), interpolate = TRUE)
png('bb6s.png')
scatter_plot(random_bb6,'#DC143C')
dev.off()

## png 
##   2

bb6s <- rasterGrob(readPNG('bb6s.png'), interpolate = TRUE)
bb1p <- persp_plot(cop_bb1, "bb1p.png","#6A5ACD")
bb6p <- persp_plot(cop_bb6, "bb6p.png","#DC143C")
legend <- legendGrob(labels = c("BB1","BB6"), pch = 15,
                    gp = gpar(col = c('#6A5ACD','#DC143C'), fill = c('#6A5ACD','#DC143C')))
grid.arrange(bb1s, bb1p,bb6s, bb6p, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình: Biểu đồ phân tán và  PDF của Copula BB1 và BB6", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
cop_bb7 <- VC2copula::BB7Copula(param = c(2,4.5))
random_bb7 <- rCopula(copula = cop_bb7,n = 7600)
cop_bb8 <- VC2copula::BB8Copula(param = c(4,0.8))
random_bb8 <- rCopula(copula = cop_bb8,n = 7600)
png('bb7s.png')
scatter_plot(random_bb7,'#FFA07A')
dev.off()

## png 
##   2

bb7s <- rasterGrob(readPNG('bb7s.png'), interpolate = TRUE)
png('bb8s.png')
scatter_plot(random_bb8,'#C0C0C0')
dev.off()

## png 
##   2

bb8s <- rasterGrob(readPNG('bb8s.png'), interpolate = TRUE)
bb7p <- persp_plot(cop_bb7, "bb7p.png","#FFA07A")
bb8p <- persp_plot(cop_bb8, "bb8p.png","#C0C0C0")
legend <- legendGrob(labels = c("BB7","BB8"), pch = 15,
                    gp = gpar(col = c('#FFA07A','#C0C0C0'), fill = c('#FFA07A','#C0C0C0')))
grid.arrange(bb7s, bb7p,bb8s, bb8p, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình: Biểu đồ phân tán và  PDF của Copula BB7 và BB8", 
                 gp = gpar(fontsize = 15, font = 2))
             )

# Kết quả ## Thống kê mô tả

library(xlsx)

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

library(readxl)
data <- read_excel("D:/Hà/file excel/MOHINHNGAUNHIEN2.xlsx")
mydata <- data[,2:3]

data1 <- as.data.frame(data)
h <- data1 %>% summarise(Min = min(SP500),
                      Max = max(SP500),
                      Mean = mean(SP500),
                      StDev = sd(SP500),
                      Skewness = skewness(SP500),
                      Kurtosis = kurtosis(SP500))
k <- data1 %>% summarise(Min = min(FPT),
                      Max = max(FPT),
                      Mean = mean(FPT),
                      StDev = sd(FPT),
                      Skewness = skewness(FPT),
                      Kurtosis = kurtosis(FPT))
l <- rbind(h,k)
rownames(l) <- c('SP500','FPT')
kable(l, format = 'html', caption = 'Bảng : Thống kê mô tả chuỗi TSLN', table.attr = "style='width:100%;'") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bảng : Thống kê mô tả chuỗi TSLN

	Min	Max	Mean	StDev	Skewness	Kurtosis
SP500	-0.1396504	0.0951132	0.0002559	0.0134855	-2.0896071	19.822935
FPT	-0.0293705	0.0742481	-0.0004911	0.0081411	0.7639885	7.347041

Phân tích biến động

ggplot(data1, aes(x = Date, y = SP500)) +
  geom_line(color = "black") +
  labs(title = "Biến động tỷ suất lợi nhuận của chỉ số S&P 500",
       x = "Ngày",
       y = "Tỷ suất lợi nhuận") +
  theme_minimal()

ggplot(data1, aes(x = Date, y = FPT)) +
  geom_line(color = "black") +
  labs(title = "Biến động tỷ suất lợi nhuận của chỉ số FPT",
       x = "Ngày",
       y = "Tỷ suất lợi nhuận") +
  theme_minimal()

# Tạo một dataframe mới để vẽ đồ thị ghép
data_plot <- data1 %>%
  pivot_longer(cols = c(SP500, FPT), names_to = "Index", values_to = "Return")

# Vẽ đồ thị ghép
ggplot(data_plot, aes(x = Date, y = Return)) +
  geom_line() +
  facet_wrap(~ Index, ncol = 2) +  
  labs(title = "Biến động tỷ suất lợi nhuận của các chỉ số",
       x = "Ngày",
       y = "Tỷ suất lợi nhuận") +
  theme_minimal()

## Kiểm định thống kê

#Kiểm định tính dừng
adf_SP500 <-adf.test(data$SP500)

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

print(adf.test(data$SP500))

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

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

adf_FPT <- adf.test(data$FPT)

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

print(adf.test(data$FPT))

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

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

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

arima_FPT <- autoarfima(data$FPT,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

re_SP500 <- arima_SP500$fit@fit$residuals
re_FPT <- arima_FPT$fit@fit$residuals
#Kiểm định tương quan chuỗi bậc 2 cho phần dư
lj_SP500 <- Box.test(re_SP500,type = 'Ljung-Box', lag = 2)
lj_FPT <- Box.test(re_FPT,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_SP5002 <- Box.test(re_SP500^2,type = 'Ljung-Box', lag = 2)
lj_FPT2 <- Box.test(re_FPT^2,type = 'Ljung-Box', lag = 2)

#Kiểm định hiệu ứng ARCH
ar_SP500 <- ArchTest(re_SP500, lags = 2)
ar_FPT <- ArchTest(re_FPT, 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_SP500 = c(adf_SP500$p.value, jq_SP500$p.value, lj_SP500$p.value, lj_SP5002$p.value, ar_SP500$p.value),
                          P_value_FPT = c(adf_FPT$p.value, jq_FPT$p.value, lj_FPT$p.value, lj_FPT2$p.value, ar_FPT$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 = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng : Kết quả các kiểm định

Test	P_value_SP500	P_value_FPT
ADF	0.0100000	0.0100000
J-B	0.0000000	0.0000000
Q(2)	0.9035411	0.6245544
Q(2)^2	0.0000000	0.0003582
ARCH(2)	0.0000000	0.0006392

#Trình bày kết quả
test_result1 <- data.frame(Test = c("ADF","J-B","Q(2)","Q(2)^2", "ARCH(2)"),
                          statistic_SP500 = c(adf_SP500$statistic, jq_SP500$statistic, lj_SP500$statistic, lj_SP5002$statistic, ar_SP500$statistic),
                          statistic_FPT = c(adf_FPT$statistic, jq_FPT$statistic, lj_FPT$statistic, lj_FPT2$statistic, ar_FPT$statistic))
kable(test_result1, 
  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 = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng : Kết quả các kiểm định

Test	statistic_SP500	statistic_FPT
ADF	-10.2240005	-11.6331367
J-B	24830.0812910	3406.9808223
Q(2)	0.2028674	0.9414337
Q(2)^2	179.5306021	15.8686400
ARCH(2)	149.2965324	14.7106659

Hệ số tương quan

pearson <- cor(data$SP500,data$FPT, method="pearson")
spearman <- cor(data$SP500,data$FPT, method="spearman")
kendall <- cor(data$SP500,data$FPT, method="kendall")

#Trình bày kết quả
relat <- data.frame('Tương quan' = 'SP500-FPT',
                    Pearson = pearson,
                    Spearman = spearman,
                    Kendall = kendall)
kable(relat, 
      col.names = c("Phương pháp", "Pearson", "Spearman","Kendall"),
      caption = "Bảng : Kết quả hệ số tương quan",
      format = 'html',
      align = c("l", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Bảng : Kết quả hệ số tương quan

Phương pháp	Pearson	Spearman	Kendall
SP500-FPT	-0.0648805	-0.0254924	-0.0168918

par(mfrow = c(1,2))
corr <- cor(mydata)
ggcorrplot(corr, hc.order = TRUE,
   ggtheme = ggplot2::theme_gray,  
   colors = c("#6D9EC1", "white", "#FFB6C1"),  
   lab = TRUE, 
   lab_col = 'black',  
   title = 'Hình: Trực quan hóa hệ số tương quan với phương pháp Pearson')

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

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

## Ước lượng mô hình phân phối biên cho mỗi chuỗi lợi suất ### GJR-GARCH(11)SP500

sp500.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")
sp500.garch11n.fit <- ugarchfit(spec = sp500.garch11n.spec, data = data$SP500)
sp500.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") 
sp500.garch11t.fit <- ugarchfit(spec = sp500.garch11t.spec, data = data$SP500)
sp500.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") 
sp500.garch11st.fit <- ugarchfit(spec = sp500.garch11st.spec, data = data$SP500)
sp500.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") 
sp500.garch11g.fit <- ugarchfit(spec = sp500.garch11g.spec, data = data$SP500)
sp500.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") 
sp500.garch11sg.fit <- ugarchfit(spec = sp500.garch11sg.spec, data = data$SP500)

GJR-GARCH(12)SP500

sp500.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")
sp500.garch12n.fit <- ugarchfit(spec = sp500.garch12n.spec, data = data$SP500)
sp500.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") 
sp500.garch12t.fit <- ugarchfit(spec = sp500.garch12t.spec, data = data$SP500)
sp500.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") 
sp500.garch12st.fit <- ugarchfit(spec = sp500.garch12st.spec, data = data$SP500)
sp500.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") 
sp500.garch12g.fit <- ugarchfit(spec = sp500.garch12g.spec, data = data$SP500)
sp500.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") 
sp500.garch12sg.fit <- ugarchfit(spec = sp500.garch12sg.spec, data = data$SP500)

GJR-GARCH(21)SP500

sp500.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")
sp500.garch21n.fit <- ugarchfit(spec = sp500.garch21n.spec, data = data$SP500)
sp500.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") 
sp500.garch21t.fit <- ugarchfit(spec = sp500.garch21t.spec, data = data$SP500)
sp500.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") 
sp500.garch21st.fit <- ugarchfit(spec = sp500.garch21st.spec, data = data$SP500)
sp500.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") 
sp500.garch21g.fit <- ugarchfit(spec = sp500.garch21g.spec, data = data$SP500)
sp500.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") 
sp500.garch21sg.fit <- ugarchfit(spec = sp500.garch21sg.spec, data = data$SP500)

GJR-GARCH(22)SP500

sp500.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")
sp500.garch22n.fit <- ugarchfit(spec = sp500.garch22n.spec, data = data$SP500)
sp500.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") 
sp500.garch22t.fit <- ugarchfit(spec = sp500.garch22t.spec, data = data$SP500)
sp500.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") 
sp500.garch22st.fit <- ugarchfit(spec = sp500.garch22st.spec, data = data$SP500)
sp500.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") 
sp500.garch22g.fit <- ugarchfit(spec = sp500.garch22g.spec, data = data$SP500)
sp500.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") 
sp500.garch22sg.fit <- ugarchfit(spec = sp500.garch22sg.spec, data = data$SP500)

GJR-GARCH(11)FPT

fpt.garch11n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "norm")
fpt.garch11n.fit <- ugarchfit(spec = fpt.garch11n.spec, data = data$FPT)
fpt.garch11t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "std") 
fpt.garch11t.fit <- ugarchfit(spec = fpt.garch11t.spec, data = data$FPT)
fpt.garch11st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sstd") 
fpt.garch11st.fit <- ugarchfit(spec = fpt.garch11st.spec, data = data$FPT)
fpt.garch11g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "ged") 
fpt.garch11g.fit <- ugarchfit(spec = fpt.garch11g.spec, data = data$FPT)
fpt.garch11sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sged") 
fpt.garch11sg.fit <- ugarchfit(spec = fpt.garch11sg.spec, data = data$FPT)

GJR-GARCH(12)FPT

fpt.garch12n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "norm")
fpt.garch12n.fit <- ugarchfit(spec = fpt.garch12n.spec, data = data$FPT)
fpt.garch12t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "std") 
fpt.garch12t.fit <- ugarchfit(spec = fpt.garch12t.spec, data = data$FPT)
fpt.garch12st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sstd") 
fpt.garch12st.fit <- ugarchfit(spec = fpt.garch12st.spec, data = data$FPT)
fpt.garch12g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "ged") 
fpt.garch12g.fit <- ugarchfit(spec = fpt.garch12g.spec, data = data$FPT)
fpt.garch12sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sged") 
fpt.garch12sg.fit <- ugarchfit(spec = fpt.garch12sg.spec, data = data$FPT)

GJR-GARCH(21)FPT

fpt.garch21n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "norm")
fpt.garch21n.fit <- ugarchfit(spec = fpt.garch21n.spec, data = data$FPT)
fpt.garch21t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "std") 
fpt.garch21t.fit <- ugarchfit(spec = fpt.garch21t.spec, data = data$FPT)
fpt.garch21st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sstd") 
fpt.garch21st.fit <- ugarchfit(spec = fpt.garch21st.spec, data = data$FPT)
fpt.garch21g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "ged") 
fpt.garch21g.fit <- ugarchfit(spec = fpt.garch21g.spec, data = data$FPT)
fpt.garch21sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sged") 
fpt.garch21sg.fit <- ugarchfit(spec = fpt.garch21sg.spec, data = data$FPT)

GJR-GARCH(22)FPT

fpt.garch22n.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "norm") 
fpt.garch22n.fit <- ugarchfit(spec = fpt.garch22n.spec, data = data$FPT)
fpt.garch22t.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "std") 
fpt.garch22t.fit <- ugarchfit(spec = fpt.garch22t.spec, data = data$FPT)
fpt.garch22st.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sstd") 
fpt.garch22st.fit <- ugarchfit(spec = fpt.garch22st.spec, data = data$FPT)
fpt.garch22g.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), mean.model 
= list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "ged") 
fpt.garch22g.fit <- ugarchfit(spec = fpt.garch22g.spec, data =data$FPT)
fpt.garch22sg.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE), distribution.model = "sged") 
fpt.garch22sg.fit <- ugarchfit(spec = fpt.garch22sg.spec, data = data$FPT)

Lựa chọn mô hình biên phù hợp nhất cho chuỗi FPT

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

##       garch11n  garch11t garch11st  garch11g garch11sg  garch12n garch12st
## [1,] -6.973104 -7.043321 -7.042381 -7.041838 -7.040476 -6.971775 -7.041155
## [2,] -6.947646 -7.014226 -7.009650 -7.012744 -7.007745 -6.942681 -7.004786
## [3,] -6.973150 -7.043381 -7.042457 -7.041899 -7.040553 -6.971836 -7.041249
## [4,] -6.963605 -7.032464 -7.030168 -7.030982 -7.028263 -6.960919 -7.027584
##       garch12g garch12sg  garch21n  garch21t garch21st  garch21g garch21sg
## [1,] -7.040507 -7.039150 -6.970619 -7.040701 -7.039718 -7.039136 -7.037771
## [2,] -7.007775 -7.002781 -6.937887 -7.004332 -6.999713 -7.002767 -6.997766
## [3,] -7.040583 -7.039244 -6.970695 -7.040795 -7.039831 -7.039230 -7.037884
## [4,] -7.028293 -7.025579 -6.958405 -7.027130 -7.024790 -7.025565 -7.022843
##       garch22n  garch22t garch22st  garch22g garch22sg
## [1,] -6.969479 -7.039405 -7.039687 -7.037769 -7.036407
## [2,] -6.933110 -6.999400 -6.996045 -6.997763 -6.992765
## [3,] -6.969573 -7.039519 -7.039822 -7.037882 -7.036542
## [4,] -6.955908 -7.024478 -7.023402 -7.022841 -7.020122

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

## [1] "garch11t"

Lựa chọn mô hình biên phù hợp nhất cho chuỗi SP500

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

##       garch11n  garch11t garch11st  garch11g garch11sg  garch12n  garch12t
## [1,] -6.375299 -6.416092 -6.421371 -6.402325 -6.414656 -6.371281 -6.414714
## [2,] -6.342567 -6.379723 -6.381365 -6.365957 -6.374651 -6.334912 -6.374709
## [3,] -6.375375 -6.416186 -6.421484 -6.402419 -6.414770 -6.371375 -6.414828
## [4,] -6.363085 -6.402521 -6.406443 -6.388755 -6.399729 -6.357710 -6.399787
##      garch12st  garch12g garch12sg  garch21n  garch21t garch21st  garch21g
## [1,] -6.419993 -6.408829 -6.409725 -6.375534 -6.412928 -6.420412 -6.405582
## [2,] -6.376351 -6.368823 -6.366083 -6.335529 -6.369286 -6.373133 -6.361940
## [3,] -6.420128 -6.408942 -6.409860 -6.375648 -6.413063 -6.420570 -6.405717
## [4,] -6.403709 -6.393901 -6.393440 -6.360606 -6.396643 -6.402770 -6.389297
##      garch21sg  garch22n  garch22t garch22st  garch22g garch22sg
## [1,] -6.414454 -6.374457 -6.413498 -6.419035 -6.404149 -6.409190
## [2,] -6.367175 -6.330814 -6.366219 -6.368119 -6.356870 -6.358275
## [3,] -6.414612 -6.374592 -6.413657 -6.419218 -6.404308 -6.409374
## [4,] -6.396812 -6.358172 -6.395856 -6.400036 -6.386508 -6.390192

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

## [1] "garch11st"

Ước lượng mô hình phân phối biên cho mỗi chuỗi lợi suất

mar_model <- data.frame(
  rate = c("S&P500", "FPT"),
  Format = c("ARMA(2,2)-GJR-GARCH(1,1)-Skewed Student t", "ARMA(1,1)-GJR-GARCH(1,1)-Student t"))#
kable(mar_model, col.names = c("Chuỗi lợi suất", "Dạng mô hình phân phối biên"), 
      caption = "Bảng : Mô hình phân phối biên tối ưu", format = "html", 
      table.attr = "style='width:100%;'") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bảng : Mô hình phân phối biên tối ưu

Chuỗi lợi suất	Dạng mô hình phân phối biên
S&P500	ARMA(2,2)-GJR-GARCH(1,1)-Skewed Student t
FPT	ARMA(1,1)-GJR-GARCH(1,1)-Student t

Kết quả ước lượng mô hình phân phối biên

fpt.garch11t.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu     -0.000568    0.000154 -3.67431 0.000238
## ar1     0.665156    0.447029  1.48795 0.136765
## ma1    -0.684297    0.435857 -1.57001 0.116414
## omega   0.000003    0.000004  0.81552 0.414772
## alpha1  0.094745    0.036997  2.56090 0.010440
## beta1   0.841118    0.034573 24.32883 0.000000
## gamma1  0.034457    0.039482  0.87271 0.382821
## shape   4.662810    0.658883  7.07684 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu     -0.000568    0.000204 -2.77861 0.005459
## ar1     0.665156    0.336614  1.97602 0.048153
## ma1    -0.684297    0.331087 -2.06682 0.038751
## omega   0.000003    0.000028  0.12104 0.903655
## alpha1  0.094745    0.168867  0.56107 0.574753
## beta1   0.841118    0.191187  4.39944 0.000011
## gamma1  0.034457    0.063978  0.53858 0.590180
## shape   4.662810    2.089916  2.23110 0.025675
## 
## LogLikelihood : 5121.451 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -7.0433
## Bayes        -7.0142
## Shibata      -7.0434
## Hannan-Quinn -7.0325
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.0470  0.8284
## Lag[2*(p+q)+(p+q)-1][5]    0.5191  1.0000
## Lag[4*(p+q)+(p+q)-1][9]    1.5195  0.9968
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1416  0.7067
## Lag[2*(p+q)+(p+q)-1][5]    0.9781  0.8646
## Lag[4*(p+q)+(p+q)-1][9]    1.8375  0.9235
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.7008 0.500 2.000  0.4025
## ARCH Lag[5]    1.3183 1.440 1.667  0.6414
## ARCH Lag[7]    1.6008 2.315 1.543  0.8009
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.8242
## Individual Statistics:              
## mu     0.06889
## ar1    0.03789
## ma1    0.04027
## omega  0.86524
## alpha1 0.23278
## beta1  0.32221
## gamma1 0.26251
## shape  0.80935
## 
## 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.3591 0.1743    
## Negative Sign Bias  1.3746 0.1695    
## Positive Sign Bias  0.5715 0.5678    
## Joint Effect        2.5035 0.4747    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     24.91    0.1633725
## 2    30     55.23    0.0023280
## 3    40     72.13    0.0009804
## 4    50     84.02    0.0013669
## 
## 
## Elapsed time : 1.036947

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(fpt.garch11t.fit)
kable(as.data.frame(fit1), 
  caption = "Bảng: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Student của biến FPT", 
  format = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Student của biến FPT

	Estimate	Std. Error	Pr(>|z|)
mu	-0.0005676	0.0001545	0.0002385
ar1	0.6651562	0.4470294	0.1367648
ma1	-0.6842973	0.4358565	0.1164137
omega	0.0000033	0.0000041	0.4147723
alpha1	0.0947452	0.0369968	0.0104401
beta1	0.8411180	0.0345729	0.0000000
gamma1	0.0344567	0.0394825	0.3828214
shape	4.6628104	0.6588831	0.0000000

sp500.garch11st.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001259    0.000226   5.5618        0
## ar1    -1.582113    0.034454 -45.9196        0
## ar2    -0.789746    0.049455 -15.9691        0
## ma1     1.586436    0.030680  51.7096        0
## ma2     0.804360    0.053222  15.1133        0
## omega   0.000012    0.000001  14.6627        0
## alpha1  0.652706    0.078109   8.3563        0
## beta1   0.616491    0.024308  25.3613        0
## gamma1 -0.568012    0.084430  -6.7276        0
## skew    0.875081    0.032978  26.5353        0
## shape   7.306058    0.995423   7.3397        0
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.001259    0.000201    6.2592        0
## ar1    -1.582113    0.014232 -111.1649        0
## ar2    -0.789746    0.036618  -21.5673        0
## ma1     1.586436    0.010365  153.0538        0
## ma2     0.804360    0.041434   19.4128        0
## omega   0.000012    0.000001   10.7891        0
## alpha1  0.652706    0.086901    7.5110        0
## beta1   0.616491    0.026287   23.4520        0
## gamma1 -0.568012    0.093420   -6.0802        0
## skew    0.875081    0.037001   23.6503        0
## shape   7.306058    1.159339    6.3019        0
## 
## LogLikelihood : 4672.915 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.4214
## Bayes        -6.3814
## Shibata      -6.4215
## Hannan-Quinn -6.4064
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.971  0.1603
## Lag[2*(p+q)+(p+q)-1][11]     4.846  0.9783
## Lag[4*(p+q)+(p+q)-1][19]     8.155  0.7723
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.09802  0.7542
## Lag[2*(p+q)+(p+q)-1][5]   1.50071  0.7397
## Lag[4*(p+q)+(p+q)-1][9]   3.72640  0.6355
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.07938 0.500 2.000  0.7781
## ARCH Lag[5]   1.13094 1.440 1.667  0.6945
## ARCH Lag[7]   2.44625 2.315 1.543  0.6235
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  30.9132
## Individual Statistics:              
## mu     0.16434
## ar1    0.06785
## ar2    0.13976
## ma1    0.05498
## ma2    0.12427
## omega  9.91334
## alpha1 0.20138
## beta1  0.09690
## gamma1 0.13067
## skew   0.13891
## shape  0.10211
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.49 2.75 3.27
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value    prob sig
## Sign Bias           1.4263 0.15398    
## Negative Sign Bias  0.8575 0.39130    
## Positive Sign Bias  1.7059 0.08824   *
## Joint Effect        3.7165 0.29374    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     17.26       0.5725
## 2    30     24.03       0.7273
## 3    40     41.99       0.3424
## 4    50     58.54       0.1651
## 
## 
## Elapsed time : 3.660176

fit2 <- extract_garch_results(sp500.garch11st.fit)
kable(as.data.frame(fit2), 
  caption = "Bảng : Kết quả mô hình ARMA(2,2)-GJR-Garch(1,1)-Skewed Student của biến SP500", 
  format = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng : Kết quả mô hình ARMA(2,2)-GJR-Garch(1,1)-Skewed Student của biến SP500

	Estimate	Std. Error	Pr(>|z|)
mu	0.0012589	0.0002263	0
ar1	-1.5821125	0.0344540	0
ar2	-0.7897456	0.0494548	0
ma1	1.5864362	0.0306797	0
ma2	0.8043598	0.0532221	0
omega	0.0000122	0.0000008	0
alpha1	0.6527060	0.0781093	0
beta1	0.6164907	0.0243084	0
gamma1	-0.5680120	0.0844296	0
skew	0.8750805	0.0329780	0
shape	7.3060579	0.9954227	0

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

sp500.res <- residuals(sp500.garch11st.fit)/sigma(sp500.garch11st.fit) 
a <- fitdist(distribution = "sstd", sp500.res, control = list()) 
u <- pdist(distribution = "sstd", q = sp500.res, mu = a$pars['mu'] , sigma = a$pars['sigma'],
           skew = a$pars['skew'], shape = a$pars['shape'])
fpt.res <- residuals(fpt.garch11t.fit)/sigma(fpt.garch11t.fit)
b <- fitdist(distribution = "std", fpt.res, control = list())
v <- pdist(distribution = "std", q = fpt.res, mu = b$pars['mu'] , sigma = b$pars['sigma'],
           skew = b$pars['skew'], shape = b$pars['shape'])
#test with Anderson-Darling
ad_sp500 <- ad.test(u, "punif")
ad_fpt <- ad.test(v, "punif")
#Kiểm định Cramer-von Mises
cvm_sp500 <- cvm.test(u, "punif")
cvm_fpt <- cvm.test(v, "punif")
# Kiem dinh ks-test
ks_sp500 <- ks.test(u, "punif")
ks_fpt <- ks.test(v, "punif")
#Trình bày kết quả
test <- data.frame(test = c('P_value SP500' ,'P_value FPT'),
                    AD = c(ad_sp500$p.value, ad_fpt$p.value),
                    CVM = c(cvm_sp500$p.value, cvm_fpt$p.value),
                   KS = c(ks_sp500$p.value,ks_fpt$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 = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

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 SP500	0.9130240	0.8257465	0.8343235
P_value FPT	0.8222822	0.8516722	0.7608929

Kết quả ước lượng tham số Copula

x <- data$SP500
y <- data$FPT
# Chuyển đổi dữ liệu sang phân phối chuẩn (Uniform Margins)
u1 <- pobs(u)
u2 <- pobs(v)
# Khởi tạo danh sách các copula bằng mã số
copula_codes <- list(
  "gauss" = 1, "t" = 2, "clayton" = 3, "gumbel" = 4, "frank" = 5, "joe" = 6,
  "surclayton" = 13, "surgumbel" = 14, "surjoe" = 16,
  "surbb1" = 23, "surbb6" = 26, "surbb7" = 27, "surbb8" = 28,
  "bb1" = 7, "bb6" = 8, "bb7" = 9, "bb8" = 10
)
# Hàm tính các hệ số và chỉ số cho mỗi copula
results <- sapply(copula_codes, function(copula_code) {
  fit <- BiCopSelect(u1, u2, familyset = copula_code)
  
# Lấy thông tin từ mô hình copula
tau <- fit$tau
aic <- fit$AIC
bic <- fit$BIC
  
# Tính toán các hệ số phụ thuộc đuôi (tail dependence coefficients)
lambda_u <- BiCopPar2TailDep(fit$family, fit$par, fit$par2)["upper"]
lambda_l <- BiCopPar2TailDep(fit$family, fit$par, fit$par2)["lower"]
return(c(tau, aic, bic, lambda_u, lambda_l))
})

# Đặt tên cho các kết quả
results <- t(results)
colnames(results) <- c("Kendall's τ", "AIC", "BIC", "λ_upper", "λ_lower")

# In ra kết quả
print(results)

##            Kendall's τ AIC       BIC       λ_upper      λ_lower     
## gauss      -0.02993768 -1.155976 4.124721  0            0           
## t          -0.02857813 0.1115954 10.67299  2.029418e-06 2.029418e-06
## clayton    -0.03343148 -4.345226 0.9354708 0            0           
## gumbel     -0.03628659 -7.538298 -2.257601 0            0           
## frank      -0.0227935  0.3189155 5.599613  0            0           
## joe        -0.03053553 -8.594599 -3.313902 0            0           
## surclayton -0.03343148 -4.345226 0.9354708 0            0           
## surgumbel  -0.03628659 -7.538298 -2.257601 0            0           
## surjoe     -0.03053553 -8.594599 -3.313902 0            0           
## surbb1     -0.03343148 -4.345226 0.9354708 0            0           
## surbb6     -0.03053553 -8.594599 -3.313902 0            0           
## surbb7     -0.03649457 -5.494494 5.066901  0            0           
## surbb8     -0.03078011 -6.581545 3.97985   0            0           
## bb1        -0.03649457 -5.494494 5.066901  0            0           
## bb6        -0.03078011 -6.581545 3.97985   0            0           
## bb7        -0.03088114 -6.588571 3.972823  0            0           
## bb8        -0.03054195 -6.594599 3.966795  0            0

library(copula)
library(png)
library(grid)
library(rgl)

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

library(VineCopula) 
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('SP500-FPT'),
                    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 : 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 = 'html') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng : 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
SP500-FPT	Rotated Joe 90 degrees	-1.05585	0	0	0	-0.0312287	-8.212424	-2.931727

install.packages("VC2copula")

## Warning: package 'VC2copula' is in use and will not be installed

library(VC2copula)
copula <- VC2copula::r90JoeBiCopula(param = est_opt_cop$par)
# Sử dụng hàm persp để vẽ đồ thị 3D
persp(copula, dCopula, xlab = "u", ylab = "v",
      main = "Đồ thị PDF của mô hình Rotated Joe 90 degrees", ,
      col = "lightpink")