Các mô hình ngẫu nhiên

1. Chuẩn bị dữ liệu

setwd("F:/MHNN")
rm(list=ls())
library(PerformanceAnalytics)

## Loading required package: xts

## Loading required package: zoo

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

library(readxl)
library(tseries)

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

data <- read_excel("data.xlsx")
data<- na.omit(data)
KOS <- data$kospi
VNI <- data$vni

2. Đồ thị phân tán các hàm copula

#Đồ thị phân tán các hàm copula
library(copula)
library(MASS)
library(ggplot2)
library(VineCopula)

## 
## Attaching package: 'VineCopula'

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

library(ggplot2)
library(ggExtra)
library(VC2copula)

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

# Hàm vẽ mật độ copula
plot_copula_density <- function(copula_name, u, v) {
  df <- data.frame(u = u, v = v)
  p <- ggplot(df, aes(x = u, y = v)) +
    geom_point(alpha = 0.1, color = "blue") +  # Scatter plot layer
    geom_density_2d_filled(contour_var = "density", show.legend = FALSE) +
    labs(title = paste(copula_name, "Copula"), x = "u", y = "v") +
    theme_minimal()
  print(p)
}

# Số lượng mẫu
n_samples <- 1000

# Tạo mẫu ngẫu nhiên từ phân phối đồng đều
u <- runif(n_samples)
v <- runif(n_samples)

2.1 Gaussian copula

gaussian_cop <- normalCopula(0.5, dim = 2)
gaussian_samples <- rCopula(n_samples, gaussian_cop)
plot_copula_density("Gaussian", gaussian_samples[,1], gaussian_samples[,2])

2.2 Student-t Copula

student_t_cop <- tCopula(0.5, dim = 2, df = 4)
student_t_samples <- rCopula(n_samples, student_t_cop)
plot_copula_density("Student-t", student_t_samples[,1], student_t_samples[,2])

2.3 Clayton Copula

clayton_cop <- claytonCopula(2, dim = 2)
clayton_samples <- rCopula(n_samples, clayton_cop)
plot_copula_density("Clayton", clayton_samples[,1], clayton_samples[,2])

2.4 Gumbel Copula

gumbel_cop <- gumbelCopula(2, dim = 2)
gumbel_samples <- rCopula(n_samples, gumbel_cop)
plot_copula_density("Gumbel", gumbel_samples[,1], gumbel_samples[,2])

2.5 Frank Copula

frank_cop <- frankCopula(param = 5, dim = 2)
frank_samples <- rCopula(n_samples, frank_cop)
plot_copula_density("Frank", frank_samples[,1], frank_samples[,2])

2.6 Joe Copula

joe_cop <- joeCopula(2, dim = 2)
joe_samples <- rCopula(n_samples, joe_cop)
plot_copula_density("Joe", joe_samples[,1], joe_samples[,2])

2.7 Rotated Gumbel Copula (180 degrees)

rot_gumbel_samples <- 1 - gumbel_samples
plot_copula_density("Rot-Gumbel", rot_gumbel_samples[,1], rot_gumbel_samples[,2])

2.8 Rotated Clayton Copula (180 degrees)

rot_clayton_samples <- 1 - clayton_samples
plot_copula_density("Rot-Clayton", rot_clayton_samples[,1], rot_clayton_samples[,2])

2.9 BB1 Copula

bb1_cop <- BB1Copula(param = c(2, 3))
bb1_samples <- rCopula(n_samples, bb1_cop)
plot_copula_density("BB1", bb1_samples[,1], bb1_samples[,2])

2.10 BB6 Copula

bb6_cop <- BB6Copula(param = c(1.5, 1.5))
bb6_samples <- rCopula(n_samples, bb6_cop)
plot_copula_density("BB6", bb6_samples[,1], bb6_samples[,2])

2.11 BB7 Copula

bb7_cop <- BB7Copula(param = c(2, 3))
bb7_samples <- rCopula(n_samples, bb7_cop)
plot_copula_density("BB7", bb7_samples[,1], bb7_samples[,2])

2.12 BB8 Copula

bb8_cop <- BB8Copula(param = c(1.5, 0.5))
bb8_samples <- rCopula(n_samples, bb8_cop)
plot_copula_density("BB8", bb8_samples[,1], bb8_samples[,2])

3. Thống kê mô tả

max(KOS)

## [1] 0.08601244

max(VNI)

## [1] 0.04980051

min(KOS)

## [1] -0.08393665

min(VNI)

## [1] -0.06674444

mean(KOS)

## [1] 0.0001173994

mean(VNI)

## [1] 0.0001734899

sd(KOS)

## [1] 0.01181369

sd(VNI)

## [1] 0.0130814

library(e1071)

## 
## Attaching package: 'e1071'

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

skewness(KOS)

## [1] -0.0129739

skewness(VNI)

## [1] -0.8233062

kurtosis(KOS)

## [1] 6.655037

kurtosis(VNI)

## [1] 3.319956

summary_stats <- data.frame(
  Statistic = c("Max", "Min", "Mean", "Standard Deviation", "Skewness", "Kurtosis"),
  KOS = c(max(KOS), min(KOS), mean(KOS), sd(KOS), skewness(KOS), kurtosis(KOS)),
  VNI = c(max(VNI), min(VNI), mean(VNI), sd(VNI), skewness(VNI), kurtosis(VNI))
)
print(summary_stats)

##            Statistic           KOS           VNI
## 1                Max  0.0860124388  0.0498005130
## 2                Min -0.0839366516 -0.0667444425
## 3               Mean  0.0001173994  0.0001734899
## 4 Standard Deviation  0.0118136856  0.0130814021
## 5           Skewness -0.0129739045 -0.8233062135
## 6           Kurtosis  6.6550370957  3.3199557964

4. Thực hiện các kiểm định thống kê cho chuỗi VNINDEX và KOSPI

4.1 Kiểm định tính dừng

Chuỗi KOS

adf.test(KOS)

## Warning in adf.test(KOS): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  KOS
## Dickey-Fuller = -10.961, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

Chuỗi VNI

adf.test(VNI)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  VNI
## Dickey-Fuller = -11.01, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

4.2 Kiểm định phân phối chuẩn

Chuỗi KOS

jarque.bera.test(KOS)

## 
##  Jarque Bera Test
## 
## data:  KOS
## X-squared = 2662.6, df = 2, p-value < 2.2e-16

Chuỗi VNI

jarque.bera.test(VNI)

## 
##  Jarque Bera Test
## 
## data:  VNI
## X-squared = 826.13, df = 2, p-value < 2.2e-16

4.3 Kiểm định tương quan chuỗi

Chuỗi KOS

library(stats)
Box.test(KOS, lag = 2, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  KOS
## X-squared = 12.422, df = 2, p-value = 0.002007

Chuỗi VNI

Box.test(VNI, lag = 2, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  VNI
## X-squared = 8.0437, df = 2, p-value = 0.01792

4.4 Kiểm định hiệu ứng ARCH

Chuỗi KOS

library(lmtest)
library(fGarch)

## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")

## 
## Attaching package: 'fGarch'

## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     ES, VaR

library(rugarch)

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:stats':
## 
##     sigma

arch_spec1 <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_KOS <- ugarchfit(spec = arch_spec1, data = KOS)
residuals1 <- residuals(arch_KOS)
n1 <- length(residuals1)
x1 <- 1:n1
arch_lm_model1 <- lm(residuals1^2 ~ x1)
arch1 <- bptest(arch_lm_model1)
arch1

## 
##  studentized Breusch-Pagan test
## 
## data:  arch_lm_model1
## BP = 0.74466, df = 1, p-value = 0.3882

Chuỗi VNI

arch_spec2 <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_VNI <- ugarchfit(spec = arch_spec2, data = VNI)
residuals2 <- residuals(arch_VNI)
n2 <- length(residuals2)
x2 <- 1:n2
arch_lm_model2 <- lm(residuals2^2 ~ x2)
arch2 <- bptest(arch_lm_model2)
arch2

## 
##  studentized Breusch-Pagan test
## 
## data:  arch_lm_model2
## BP = 0.022695, df = 1, p-value = 0.8803

5. Hệ số tương quan

pearson_corr <- cor(KOS, VNI, method = "pearson")
pearson_corr

## [1] 0.3016879

spearman_corr <- cor(KOS, VNI, method = "spearman")
spearman_corr

## [1] 0.2784454

kendall_corr <- cor(KOS, VNI, method = "kendall")
kendall_corr

## [1] 0.1915764

correlation_stats <- data.frame(
  Method = c("Pearson", "Spearman", "Kendall"),
  Correlation = c(pearson_corr, spearman_corr, kendall_corr)
)
print(correlation_stats)

##     Method Correlation
## 1  Pearson   0.3016879
## 2 Spearman   0.2784454
## 3  Kendall   0.1915764

5.1 Kiểm định Pearson

pearson_test <- cor.test(KOS, VNI, method = "pearson")
print(pearson_test)

## 
##  Pearson's product-moment correlation
## 
## data:  KOS and VNI
## t = 11.987, df = 1435, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2539381 0.3479707
## sample estimates:
##       cor 
## 0.3016879

5.2 Kiểm định Spearman

spearman_test <- cor.test(KOS, VNI, method = "spearman")
print(spearman_test)

## 
##  Spearman's rank correlation rho
## 
## data:  KOS and VNI
## S = 356851904, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2784454

5.3 Kiểm định Kendall

kendall_test <- cor.test(KOS, VNI, method = "kendall")
print(kendall_test)

## 
##  Kendall's rank correlation tau
## 
## data:  KOS and VNI
## z = 10.88, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.1915764

5.4 Phân phối và hệ số tương quan Pearson giữa chỉ số VNINDEX và KOSPI

library(ggplot2)
library(patchwork)

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:MASS':
## 
##     area

# Biểu đồ histogram của KOS
hist_btc <- ggplot(data.frame(KOS), aes(x = KOS)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black") +
  geom_density(color = "blue") +
  labs(title = "KOS") +
  theme_minimal()

# Biểu đồ histogram của VNI
hist_eth <- ggplot(data.frame(VNI), aes(x = VNI)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black") +
  geom_density(color = "blue") +
  labs(title = "VNI") +
  theme_minimal()

# Biểu đồ scatter thể hiện mối tương quan giữa KOS và VNI
scatter_plot <- ggplot(data.frame(KOS, VNI), aes(x = KOS, y = VNI)) +
  geom_point(color = "black") +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(x = "KOS", y = "VNI") +
  theme_minimal()

# Hệ số tương quan Pearson
corr_text <- ggplot() +
  geom_text(aes(0.5, 0.5, label = sprintf("%.2f", pearson_corr)), size = 20) +
  theme_void()

# Kết hợp các biểu đồ lại với nhau bằng patchwork
final_plot <- (hist_btc + corr_text) / (scatter_plot + hist_eth)

# Hiển thị đồ thị
print(final_plot)

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

6. Biến động tỷ suất sinh lợi 2018-2023

Chuỗi KOS

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 object is masked from 'package:MASS':
## 
##     select

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

data %>% ggplot(aes(x = Date, y = kospi)) +
  geom_line(color = "black", group = 1) +
  labs(title = "KOS", x = "Ngày", y = "Độ biến động") +
  theme_minimal()

Chuỗi VNI

data %>% ggplot(aes(x = Date, y = vni)) +
  geom_line(color = "black", group = 1) +
  labs(title = "VNI", x = "Ngày", y = "Độ biến động") +
  theme_minimal()

7. Kết quả ước lượng mức độ phản ứng của TTCK Việt Nam đối với TTCK Hàn Quốc

7.1 Mô hình ARMA

Chuỗi KOS

autoarfima(KOS,ar.max = 2, ma.max = 2, criterion = "AIC", method = "full")

## $fit
## 
## *----------------------------------*
## *          ARFIMA Model Fit        *
## *----------------------------------*
## Mean Model   : ARFIMA(0,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##        Estimate  Std. Error  t value Pr(>|t|)
## ma1    0.000000          NA       NA       NA
## ma2    0.096449    0.026557   3.6318 0.000281
## sigma  0.011757    0.000219  53.6097 0.000000
## 
## Robust Standard Errors:
##        Estimate  Std. Error  t value Pr(>|t|)
## ma1    0.000000          NA       NA       NA
## ma2    0.096449    0.053752   1.7943 0.072759
## sigma  0.011757    0.000446  26.3837 0.000000
## 
## LogLikelihood : 4345.981 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.0459
## Bayes        -6.0386
## Shibata      -6.0459
## Hannan-Quinn -6.0432
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2496  0.6173
## Lag[2*(p+q)+(p+q)-1][5]    1.3343  0.9995
## Lag[4*(p+q)+(p+q)-1][9]    3.2166  0.8551
## 
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      397.9       0
## Lag[2*(p+q)+(p+q)-1][2]     499.7       0
## Lag[4*(p+q)+(p+q)-1][5]     723.9       0
## 
## 
## ARCH LM Tests
## ------------------------------------
##              Statistic DoF P-Value
## ARCH Lag[2]      416.2   2       0
## ARCH Lag[5]      443.9   5       0
## ARCH Lag[10]     477.5  10       0
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.0492
## Individual Statistics:            
## ma2   0.2210
## sigma 0.8608
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          0.61 0.749 1.07
## Individual Statistic:     0.35 0.47 0.75
## 
## 
## Elapsed time : 0.02326703 
## 
## 
## $rank.matrix
##    ar1 ar2 ma1 ma2 im arf       AIC converged
## 1    0   0   0   1  0   0 -6.045903         1
## 2    1   1   1   1  0   0 -6.045869         1
## 3    0   1   0   0  0   0 -6.045483         1
## 4    0   0   1   1  0   0 -6.044771         1
## 5    1   0   0   1  0   0 -6.044683         1
## 6    0   1   0   1  0   0 -6.044624         1
## 7    0   0   0   1  1   0 -6.044595         1
## 8    1   1   1   1  1   0 -6.044563         1
## 9    0   1   1   0  0   0 -6.044264         1
## 10   1   1   0   0  0   0 -6.044195         1
## 11   0   1   0   0  1   0 -6.044174         1
## 12   1   0   1   1  0   0 -6.044071         1
## 13   0   1   1   1  0   0 -6.043576         1
## 14   1   1   0   1  0   0 -6.043473         1
## 15   0   0   1   1  1   0 -6.043466         1
## 16   1   1   1   0  0   0 -6.043421         1
## 17   1   0   0   1  1   0 -6.043377         1
## 18   0   1   0   1  1   0 -6.043317         1
## 19   0   1   1   0  1   0 -6.042958         1
## 20   1   1   0   0  1   0 -6.042889         1
## 21   1   0   1   1  1   0 -6.042760         1
## 22   0   1   1   1  1   0 -6.042273         1
## 23   1   1   0   1  1   0 -6.042169         1
## 24   1   1   1   0  1   0 -6.042111         1
## 25   1   0   0   0  0   0 -6.037059         1
## 26   0   0   1   0  0   0 -6.037039         1
## 27   0   0   0   0  1   0 -6.037029         1
## 28   1   0   1   0  0   0 -6.037000         1
## 29   1   0   0   0  1   0 -6.035768         1
## 30   0   0   1   0  1   0 -6.035747         1
## 31   1   0   1   0  1   0 -6.035709         1

Chuỗi VNI

autoarfima(VNI,ar.max = 2, ma.max = 2, criterion = "AIC", method = "full")

## $fit
## 
## *----------------------------------*
## *          ARFIMA Model Fit        *
## *----------------------------------*
## Mean Model   : ARFIMA(0,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##        Estimate  Std. Error  t value Pr(>|t|)
## ma1    0.052910    0.026359   2.0073 0.044715
## ma2    0.050629    0.026331   1.9228 0.054505
## sigma  0.013046    0.000243  53.5916 0.000000
## 
## Robust Standard Errors:
##        Estimate  Std. Error  t value Pr(>|t|)
## ma1    0.052910    0.030638   1.7269 0.084178
## ma2    0.050629    0.030415   1.6646 0.095989
## sigma  0.013046    0.000685  19.0394 0.000000
## 
## LogLikelihood : 4196.843 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.8369
## Bayes        -5.8259
## Shibata      -5.8370
## Hannan-Quinn -5.8328
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0000666  0.9935
## Lag[2*(p+q)+(p+q)-1][5] 0.2273164  1.0000
## Lag[4*(p+q)+(p+q)-1][9] 1.3817937  0.9982
## 
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      57.75 2.975e-14
## Lag[2*(p+q)+(p+q)-1][2]     89.16 0.000e+00
## Lag[4*(p+q)+(p+q)-1][5]    158.60 0.000e+00
## 
## 
## ARCH LM Tests
## ------------------------------------
##              Statistic DoF P-Value
## ARCH Lag[2]      100.1   2       0
## ARCH Lag[5]      135.8   5       0
## ARCH Lag[10]     173.2  10       0
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  0.7041
## Individual Statistics:             
## ma1   0.06383
## ma2   0.31161
## sigma 0.33962
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          0.846 1.01 1.35
## Individual Statistic:     0.35 0.47 0.75
## 
## 
## Elapsed time : 0.02058506 
## 
## 
## $rank.matrix
##    ar1 ar2 ma1 ma2 im arf       AIC converged
## 1    0   0   1   1  0   0 -5.836942         1
## 2    1   0   0   1  0   0 -5.836927         1
## 3    0   1   1   0  0   0 -5.836887         1
## 4    1   1   0   0  0   0 -5.836875         1
## 5    1   0   1   0  0   0 -5.836546         1
## 6    1   0   0   0  0   0 -5.836039         1
## 7    0   0   1   0  0   0 -5.835765         1
## 8    0   0   1   1  1   0 -5.835698         1
## 9    1   0   0   1  1   0 -5.835683         1
## 10   0   1   1   0  1   0 -5.835644         1
## 11   1   1   0   0  1   0 -5.835631         1
## 12   0   1   0   0  0   0 -5.835559         1
## 13   0   0   0   1  0   0 -5.835532         1
## 14   1   0   1   1  0   0 -5.835530         1
## 15   1   1   1   0  0   0 -5.835496         1
## 16   0   1   1   1  0   0 -5.835495         1
## 17   1   1   0   1  0   0 -5.835485         1
## 18   1   0   1   0  1   0 -5.835294         1
## 19   1   0   0   0  1   0 -5.834806         1
## 20   0   0   1   0  1   0 -5.834535         1
## 21   0   1   0   0  1   0 -5.834329         1
## 22   0   0   0   1  1   0 -5.834302         1
## 23   1   0   1   1  1   0 -5.834284         1
## 24   1   1   1   0  1   0 -5.834253         1
## 25   0   1   1   1  1   0 -5.834252         1
## 26   1   1   0   1  1   0 -5.834240         1
## 27   0   1   0   1  0   0 -5.834231         1
## 28   1   1   1   1  0   0 -5.834104         1
## 29   0   0   0   0  1   0 -5.833163         1
## 30   0   1   0   1  1   0 -5.832997         1
## 31   1   1   1   1  1   0 -5.832860         1

7.2 Mô hình GJR-GARCH

Chuỗi KOS

Chuỗi VNI

7.3 Lựa chọn mô hình GJR-GARCH

Chuỗi KOS

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

##               garch11n  garch11t garch11st  garch11g garch11sg  garch12n
## Akaike       -6.291504 -6.323824 -6.331609 -6.316950 -6.325407 -6.290256
## Bayes        -6.265831 -6.294483 -6.298601 -6.287610 -6.292399 -6.260916
## Shibata      -6.291551 -6.323885 -6.331687 -6.317012 -6.325485 -6.290318
## Hannan-Quinn -6.281920 -6.312870 -6.319286 -6.305996 -6.313084 -6.279302
##               garch12t garch12st  garch12g garch12sg  garch21n  garch21t
## Akaike       -6.322434 -6.330220 -6.315568 -6.324057 -6.288844 -6.321170
## Bayes        -6.289425 -6.293544 -6.282560 -6.287381 -6.255836 -6.284495
## Shibata      -6.322511 -6.330315 -6.315646 -6.324153 -6.288922 -6.321266
## Hannan-Quinn -6.310110 -6.316527 -6.303245 -6.310365 -6.276521 -6.307478
##              garch21st  garch21g garch21sg  garch22n  garch22t garch22st
## Akaike       -6.329042 -6.314241 -6.322842 -6.287452 -6.319778 -6.327651
## Bayes        -6.288699 -6.277565 -6.282499 -6.250776 -6.279435 -6.283640
## Shibata      -6.329158 -6.314337 -6.322958 -6.287548 -6.319894 -6.327788
## Hannan-Quinn -6.313981 -6.300548 -6.307780 -6.273760 -6.304717 -6.311220
##               garch22g garch22sg
## Akaike       -6.312849 -6.321450
## Bayes        -6.272505 -6.277439
## Shibata      -6.312965 -6.321588
## Hannan-Quinn -6.297787 -6.305019

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

## [1] "garch11st"

Chuỗi VNI

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

##               garch11n  garch11t garch11st  garch11g garch11sg  garch12n
## Akaike       -6.059219 -6.185005 -6.199770 -6.174255 -6.192352 -6.060219
## Bayes        -6.033546 -6.155665 -6.166762 -6.144914 -6.159344 -6.030878
## Shibata      -6.059266 -6.185067 -6.199848 -6.174316 -6.192430 -6.060280
## Hannan-Quinn -6.049634 -6.174051 -6.187447 -6.163301 -6.180029 -6.049265
##               garch12t garch12st  garch12g garch12sg  garch21n  garch21t
## Akaike       -6.184758 -6.199960 -6.173935 -6.192264 -6.066962 -6.187025
## Bayes        -6.151750 -6.163284 -6.140927 -6.155588 -6.033954 -6.150349
## Shibata      -6.184836 -6.200056 -6.174013 -6.192360 -6.067040 -6.187121
## Hannan-Quinn -6.172435 -6.186268 -6.161612 -6.178571 -6.054639 -6.173332
##              garch21st  garch21g garch21sg  garch22n  garch22t garch22st
## Akaike       -6.202897 -6.176795 -6.195293 -6.065570 -6.185633 -6.201505
## Bayes        -6.162553 -6.140119 -6.154950 -6.028894 -6.145290 -6.157494
## Shibata      -6.203013 -6.176891 -6.195409 -6.065666 -6.185749 -6.201643
## Hannan-Quinn -6.187835 -6.163103 -6.180231 -6.051877 -6.170571 -6.185074
##               garch22g garch22sg
## Akaike       -6.175403 -6.193901
## Bayes        -6.135060 -6.149890
## Shibata      -6.175519 -6.194039
## Hannan-Quinn -6.160342 -6.177470

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

## [1] "garch21st"

7.4 Tham số ước lượng biên phù hợp nhất

Chuỗi KOS

kospi.garch11st.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(0,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu     -0.000023    0.000268 -0.084684 0.932513
## ma1     0.014277    0.027301  0.522948 0.601010
## ma2     0.045430    0.027284  1.665084 0.095896
## omega   0.000010    0.000000 25.193167 0.000000
## alpha1  0.034466    0.011737  2.936563 0.003319
## beta1   0.781227    0.019155 40.784042 0.000000
## gamma1  0.213892    0.045070  4.745720 0.000002
## skew    0.872063    0.033063 26.375637 0.000000
## shape   8.067648    1.432226  5.632945 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu     -0.000023    0.000271 -0.083632 0.933349
## ma1     0.014277    0.024638  0.579465 0.562276
## ma2     0.045430    0.026864  1.691108 0.090816
## omega   0.000010    0.000000 20.953991 0.000000
## alpha1  0.034466    0.010385  3.318971 0.000903
## beta1   0.781227    0.019403 40.263410 0.000000
## gamma1  0.213892    0.049154  4.351475 0.000014
## skew    0.872063    0.032832 26.561194 0.000000
## shape   8.067648    1.513908  5.329020 0.000000
## 
## LogLikelihood : 4558.261 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.3316
## Bayes        -6.2986
## Shibata      -6.3317
## Hannan-Quinn -6.3193
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.3942  0.5301
## Lag[2*(p+q)+(p+q)-1][5]    0.9130  1.0000
## Lag[4*(p+q)+(p+q)-1][9]    1.4437  0.9977
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01919  0.8898
## Lag[2*(p+q)+(p+q)-1][5]   1.27178  0.7957
## Lag[4*(p+q)+(p+q)-1][9]   2.20343  0.8794
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.101 0.500 2.000  0.2940
## ARCH Lag[5]     1.315 1.440 1.667  0.6424
## ARCH Lag[7]     1.537 2.315 1.543  0.8140
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  55.8794
## Individual Statistics:              
## mu     0.25573
## ma1    0.04001
## ma2    0.82376
## omega  9.79476
## alpha1 0.50782
## beta1  0.65012
## gamma1 0.64914
## skew   0.40266
## shape  0.61560
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.1 2.32 2.82
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.6062 0.1084    
## Negative Sign Bias  1.3593 0.1743    
## Positive Sign Bias  0.4555 0.6488    
## Joint Effect        2.9720 0.3960    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     12.59       0.8589
## 2    30     21.10       0.8555
## 3    40     45.28       0.2263
## 4    50     44.73       0.6467
## 
## 
## Elapsed time : 1.471031

Chuỗi VNI

vni.garch21st.fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(2,1)
## Mean Model   : ARFIMA(0,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000287    0.000300  0.95667 0.338733
## ma1     0.067107    0.028610  2.34563 0.018995
## ma2     0.050789    0.025134  2.02068 0.043312
## omega   0.000007    0.000000 17.59629 0.000000
## alpha1  0.000000    0.039145  0.00001 0.999992
## alpha2  0.025486    0.044947  0.56701 0.570707
## beta1   0.842559    0.016913 49.81858 0.000000
## gamma1  0.456838    0.113063  4.04057 0.000053
## gamma2 -0.276530    0.100532 -2.75066 0.005948
## skew    0.829943    0.031942 25.98293 0.000000
## shape   4.549656    0.515973  8.81762 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000287    0.000317  0.905178 0.365371
## ma1     0.067107    0.028389  2.363845 0.018086
## ma2     0.050789    0.025570  1.986249 0.047006
## omega   0.000007    0.000001 13.970793 0.000000
## alpha1  0.000000    0.034296  0.000011 0.999991
## alpha2  0.025486    0.039087  0.652014 0.514392
## beta1   0.842559    0.013410 62.832352 0.000000
## gamma1  0.456838    0.091607  4.986929 0.000001
## gamma2 -0.276530    0.082665 -3.345181 0.000822
## skew    0.829943    0.032720 25.364893 0.000000
## shape   4.549656    0.465581  9.771987 0.000000
## 
## LogLikelihood : 4467.781 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2029
## Bayes        -6.1626
## Shibata      -6.2030
## Hannan-Quinn -6.1878
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.5115  0.4745
## Lag[2*(p+q)+(p+q)-1][5]    0.9088  1.0000
## Lag[4*(p+q)+(p+q)-1][9]    1.7250  0.9933
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       2.945 0.08614
## Lag[2*(p+q)+(p+q)-1][8]      4.339 0.45060
## Lag[4*(p+q)+(p+q)-1][14]     6.801 0.53168
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]     1.613 0.500 2.000  0.2041
## ARCH Lag[6]     1.699 1.461 1.711  0.5607
## ARCH Lag[8]     1.884 2.368 1.583  0.7649
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  28.8981
## Individual Statistics:              
## mu     0.28856
## ma1    0.07365
## ma2    0.36455
## omega  4.28047
## alpha1 0.17612
## alpha2 0.11579
## beta1  0.12785
## gamma1 0.09438
## gamma2 0.08851
## skew   0.10807
## shape  0.16310
## 
## 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            2.557 0.010673  **
## Negative Sign Bias   2.599 0.009447 ***
## Positive Sign Bias   0.768 0.442626    
## Joint Effect        13.284 0.004061 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     20.75      0.35104
## 2    30     33.29      0.26611
## 3    40     43.33      0.29166
## 4    50     64.01      0.07352
## 
## 
## Elapsed time : 1.969304

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

Chuỗi KOS

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

## $pars
##          mu       sigma        skew       shape 
## 0.001069737 0.998058659 0.872681026 8.121847871 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2160.462 2006.472 2006.472
## 
## $lagrange
## [1] 0
## 
## $hessian
##              [,1]       [,2]        [,3]       [,4]
## [1,] 1916.6084007  420.72389 -493.547666  0.5744737
## [2,]  420.7238922 1934.84568   83.463352 14.1786825
## [3,] -493.5476656   83.46335 1135.649847  1.2064234
## [4,]    0.5744737   14.17868    1.206423  0.5550982
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 127
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.165642 secs
## 
## $vscale
## [1] 1 1 1 1 1

u <-pdist(distribution = "sstd", q = kospi.res, mu = 0.001069737 , sigma = 0.998058659, skew = 0.872681026, shape = 8.121847871)

Chuỗi VNI

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

## $pars
##           mu        sigma         skew        shape 
## -0.002596318  1.003719372  0.828651094  4.515788329 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2104.142 1935.011 1935.011
## 
## $lagrange
## [1] 0
## 
## $hessian
##            [,1]      [,2]       [,3]      [,4]
## [1,] 1983.16121  468.1805 -602.74472 18.048769
## [2,]  468.18046 1676.7981  213.00139 88.269000
## [3,] -602.74472  213.0014 1248.49163 20.592853
## [4,]   18.04877   88.2690   20.59285  7.868569
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 101
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.1344161 secs
## 
## $vscale
## [1] 1 1 1 1 1

v <-pdist(distribution = "sstd", q = vni.res, mu = -0.002596318 , sigma = 1.003719372 , skew = 0.828651094 , shape = 4.515788329)

8. Kiểm định sự phù hợp của mô hình biên

8.1 Kiểm định Anderson-Darling

library(tseries)
library(rugarch)
library(goftest)
ad.test(u, "punif")

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  u
## An = 0.27056, p-value = 0.9585

ad.test(v, "punif")

## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  v
## An = 0.50372, p-value = 0.743

8.2 Kiểm định Cramer-von Mises

cvm.test(u, "punif")

## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  u
## omega2 = 0.037476, p-value = 0.9458

cvm.test(v, "punif")

## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  v
## omega2 = 0.092068, p-value = 0.6252

8.3 Kiểm định ks-test

ks.test(u, "punif")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.014552, p-value = 0.9212
## alternative hypothesis: two-sided

ks.test(v, "punif")

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

9. Mô hình copula

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

library(VineCopula)
cop <- BiCopSelect(u, v, familyset= 1:9, selectioncrit="AIC",indeptest = FALSE, level = 0.05)
summary(cop)

## Family
## ------ 
## No:    14
## Name:  Survival Gumbel
## 
## Parameter(s)
## ------------
## par:  1.2
## 
## Dependence measures
## -------------------
## Kendall's tau:    0.16 (empirical = 0.18, p value < 0.01)
## Upper TD:         0 
## Lower TD:         0.22 
## 
## Fit statistics
## --------------
## logLik:  62.28 
## AIC:    -122.57 
## BIC:    -117.3

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

## Family
## ------ 
## No:    4
## Name:  Gumbel
## 
## Parameter(s)
## ------------
## par:  1.17  (SE = 0.02)
## 
## Dependence measures
## -------------------
## Kendall's tau:    0.14 (empirical = 0.18, p value < 0.01)
## Upper TD:         0.19 
## Lower TD:         0 
## 
## Fit statistics
## --------------
## logLik:  36.98 
## AIC:    -71.97 
## BIC:    -66.7

9.2 Trực quan hóa mô hình copula

# Tạo lưới giá trị cho u và v
u_grid <- seq(0.01, 0.99, length.out = 100)
v_grid <- seq(0.01, 0.99, length.out = 100)

# Tạo ma trận chứa các cặp (u,v)
grid <- expand.grid(u = u_grid, v = v_grid)

# Tính toán mật độ copula
density <- BiCopPDF(grid$u, grid$v, Stu)

# Chuyển đổi mật độ sang dạng ma trận
density_matrix <- matrix(density, nrow = 100, ncol = 100)
library(ggplot2)
library(reshape2)
# Chuyển dữ liệu sang định dạng cho ggplot2
density_data <- melt(density_matrix)
colnames(density_data) <- c("u", "v", "density")
density_data$u <- u_grid[density_data$u]
density_data$v <- v_grid[density_data$v]

# Tạo biểu đồ đồng mức với ggplot2
ggplot(density_data, aes(x = u, y = v, z = density)) +
  geom_tile(aes(fill = density)) +
  geom_contour(color = "white") +
  scale_fill_viridis_c() +  # Sử dụng thang màu viridis
  labs(x = "U", y = "V", title = "Biểu đồ đồng mức của hàm copula") +
  theme_minimal()

# Biểu đồ bề mặt 3D
persp(u_grid, v_grid, density_matrix, theta = 30, phi = 30, expand = 0.5, col = "lightblue", xlab = "U", ylab = "V", zlab = "Density", main = "Biểu đồ bề mặt 3D của hàm copula")