MHNN TT HNX SP500

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

1.1 GÁN DỮ LIỆU VÀO

#Sử dụng lệnh import file dữ liệu vào R
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(xlsx)
data <- read.xlsx("C:/Users/Admin/Desktop/Book1.xlsx", sheetName = 5, header = T)
LGSP <- diff(log(data$SP), lag = 1)
LGHNX <- diff(log(data$HNX), lag = 1)
mhnn <- data.frame(HNX = LGHNX, SP = LGSP)
summary(mhnn)

##       HNX                   SP            
##  Min.   :-0.1028979   Min.   :-0.1276573  
##  1st Qu.:-0.0049038   1st Qu.:-0.0038481  
##  Median : 0.0012719   Median : 0.0006183  
##  Mean   : 0.0005058   Mean   : 0.0003954  
##  3rd Qu.: 0.0072952   3rd Qu.: 0.0057139  
##  Max.   : 0.0539166   Max.   : 0.0896709

1.2 BIỂU ĐỒ NHIỄU

plot.ts(mhnn$HNX)

plot.ts(mhnn$SP)

1.3 MA TRẬN HỆ SỐ TƯƠNG QUAN

library(Hmisc)

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

res2<-rcorr(as.matrix(mhnn[,1:2]))
res2 

##      HNX   SP
## HNX 1.00 0.09
## SP  0.09 1.00
## 
## n= 2420 
## 
## 
## P
##     HNX SP
## HNX      0
## SP   0

# Hệ số tương quan Pearson
res <- cor(mhnn)
round(res, 2)

##      HNX   SP
## HNX 1.00 0.09
## SP  0.09 1.00

flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}
flattenCorrMatrix(res2$r, res2$P)

##   row column        cor           p
## 1 HNX     SP 0.09216803 5.58835e-06

symnum(res, abbr.colnames = FALSE)

##     HNX SP
## HNX 1     
## SP      1 
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1

library(corrplot)

## corrplot 0.92 loaded

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

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

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

#-Tương tự cho các cặp biến khác.
library(ggplot2)
library(ggcorrplot)
df <- dplyr::select_if(mhnn, is.numeric)
r <- cor(df, use="complete.obs")
ggcorrplot(r)

1.4 KIỂM ĐỊNH DỮ LIỆU

1.4.1 KIỂM ĐỊNH PHÂN PHỐI CHUẨN JARQUE - BERA

library(tseries)

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

# kiểm định cho chỉ số VNI 
Var1 <- mhnn$HNX
Var2 <- mhnn$SP
# Kiểm định phân phối chuẩn cho VNI 
result1 <-  jarque.bera.test(Var1)
print(result1)

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

result2 <-  jarque.bera.test(Var2)
print(result2)

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

1.4.2 KIỂM ĐỊNH TÍNH DỪNG: AUGMENTED DICKEY - FULLER

adf.test(Var1)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var1
## Dickey-Fuller = -11.495, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary

print(adf.test(Var1))

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var1
## Dickey-Fuller = -11.495, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary

adf.test(Var2)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var2
## Dickey-Fuller = -13.405, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary

print(adf.test(Var2))

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var2
## Dickey-Fuller = -13.405, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary

1.4.3 KIỂM ĐỊNH TƯƠNG QUAN CHUỖI: LJUNG - BOX

# Kiểm định tương quan chuỗi cho HNX 
library(stats)
result2 <-  Box.test(Var1, lag = 10, type = "Ljung-Box")
print(result2)

## 
##  Box-Ljung test
## 
## data:  Var1
## X-squared = 40.36, df = 10, p-value = 1.464e-05

result3 <-  Box.test(Var2, lag = 10, type = "Ljung-Box")
print(result3)

## 
##  Box-Ljung test
## 
## data:  Var2
## X-squared = 233.41, df = 10, p-value < 2.2e-16

1.5 MÔ HÌNH ARIMA

1.5.1 BIẾN HNX

library(tseries)
library(forecast)
modeld <- auto.arima(mhnn$HNX)
modeld

## Series: mhnn$HNX 
## ARIMA(2,0,1) with zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.9688  -0.0141  -0.9248
## s.e.  0.0597   0.0249   0.0561
## 
## sigma^2 = 0.0001742:  log likelihood = 7040.67
## AIC=-14073.34   AICc=-14073.33   BIC=-14050.18

1.5.2 BIẾN S&P500

library(tseries)
library(forecast)
modeld1 <- auto.arima(mhnn$SP)
modeld1

## Series: mhnn$SP 
## ARIMA(5,0,5) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2      ar3     ar4     ar5     ma1      ma2     ma3
##       -0.2912  0.2460  -0.3036  0.4421  0.7750  0.1763  -0.2227  0.3531
## s.e.   0.0488  0.0403   0.0314  0.0436  0.0409  0.0540   0.0348  0.0256
##           ma4      ma5   mean
##       -0.5520  -0.6815  4e-04
## s.e.   0.0398   0.0478  1e-04
## 
## sigma^2 = 0.000121:  log likelihood = 7485.09
## AIC=-14946.18   AICc=-14946.05   BIC=-14876.68

1.6 KIỂM ĐỊNH HIỆU ỨNG ARCH

1.6.1 BIẾN HNX

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

# Kiểm định hiệu ứng ARCH - LM  cho chỉ số chứng khoán 
arch_spec <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_HNX <- ugarchfit(spec = arch_spec, data = mhnn$HNX)

residuals <- residuals(arch_HNX)
n <- length(residuals)
x <- 1:n
# Tạo mô hình tuyến tính
arch_lm_model <- lm(residuals^2 ~ x)
# Kiểm định hiệu ứng ARCH-LM
arch2<- bptest(arch_lm_model)
# Hiển thị kết quả
arch2

## 
##  studentized Breusch-Pagan test
## 
## data:  arch_lm_model
## BP = 5.1787, df = 1, p-value = 0.02287

1.6.2 BIẾN S&P 500

library(lmtest)
library(fGarch)
library(rugarch)
# Kiểm định hiệu ứng ARCH - LM  cho chỉ số chứng khoán 
arch_spec1 <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_SP <- ugarchfit(spec = arch_spec1, data = mhnn$SP)

residuals1 <- residuals(arch_SP)
n1 <- length(residuals1)
x1 <- 1:n
# Tạo mô hình tuyến tính
arch_lm_model1 <- lm(residuals1^2 ~ x)
arch1<- bptest(arch_lm_model1)
# Hiển thị kết quả
arch1

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

1.7 MÔ HÌNH GARCH

1.7.1 BIẾN HNX

library(rugarch)
HNXts<- ts(mhnn$HNX)
head(HNXts)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]  0.0005886682  0.0124279690  0.0094006043  0.0074573700  0.0165784881
## [6] -0.0056369935

library(lmtest)
HNXspec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1,1)), mean.model = list(armaOrder = c(2,1), include.mean =TRUE), distribution.model = 'norm')
print(HNXspec)

## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : gjrGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(2,0,1)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE

HNXfit <- ugarchfit(spec = HNXspec, HNXts)
print(HNXfit)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.000960    0.000240      3.9976 0.000064
## ar1     1.061405    0.019749     53.7459 0.000000
## ar2    -0.065973    0.020038     -3.2924 0.000993
## ma1    -0.987972    0.000029 -34216.0297 0.000000
## omega   0.000005    0.000001      7.4327 0.000000
## alpha1  0.127344    0.020647      6.1677 0.000000
## beta1   0.808513    0.011944     67.6927 0.000000
## gamma1  0.088521    0.026139      3.3866 0.000708
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.000960    0.000910  1.0555e+00 0.291203
## ar1     1.061405    0.082958  1.2795e+01 0.000000
## ar2    -0.065973    0.082155 -8.0303e-01 0.421956
## ma1    -0.987972    0.000052 -1.9073e+04 0.000000
## omega   0.000005    0.000002  2.3054e+00 0.021146
## alpha1  0.127344    0.042337  3.0078e+00 0.002631
## beta1   0.808513    0.018804  4.2996e+01 0.000000
## gamma1  0.088521    0.048519  1.8245e+00 0.068082
## 
## LogLikelihood : 7435.31 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.1383
## Bayes        -6.1191
## Shibata      -6.1383
## Hannan-Quinn -6.1313
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                        2.44 1.183e-01
## Lag[2*(p+q)+(p+q)-1][8]       8.41 1.164e-07
## Lag[4*(p+q)+(p+q)-1][14]     12.68 1.549e-02
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1683  0.6816
## Lag[2*(p+q)+(p+q)-1][5]    0.5930  0.9424
## Lag[4*(p+q)+(p+q)-1][9]    1.3950  0.9637
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     0.590 0.500 2.000  0.4424
## ARCH Lag[5]     0.609 1.440 1.667  0.8510
## ARCH Lag[7]     1.276 2.315 1.543  0.8649
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.9549
## Individual Statistics:              
## mu     0.16323
## ar1    0.07147
## ar2    0.06335
## ma1    0.19998
## omega  0.28210
## alpha1 0.40521
## beta1  0.35074
## gamma1 0.28002
## 
## 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           0.8960 0.37036    
## Negative Sign Bias  0.8647 0.38726    
## Positive Sign Bias  1.8518 0.06418   *
## Joint Effect        4.5986 0.20366    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     129.1    2.198e-18
## 2    30     142.3    6.845e-17
## 3    40     158.6    2.350e-16
## 4    50     178.7    1.184e-16
## 
## 
## Elapsed time : 1.013943

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

HNXst1<- ugarchfit(HNXst.spec,HNXts) 
print(HNXst1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.004898    0.001624      3.0158 0.002563
## ar1     1.035907    0.001115    928.8791 0.000000
## ar2    -0.037616    0.001002    -37.5513 0.000000
## ma1    -0.986804    0.000046 -21474.7043 0.000000
## omega   0.000005    0.000003      1.5876 0.112373
## alpha1  0.099148    0.034924      2.8389 0.004526
## beta1   0.826372    0.023249     35.5449 0.000000
## gamma1  0.094860    0.039799      2.3835 0.017150
## shape   4.526164    0.462573      9.7848 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.004898    0.001817  2.6957e+00 0.007024
## ar1     1.035907    0.001242  8.3383e+02 0.000000
## ar2    -0.037616    0.001053 -3.5713e+01 0.000000
## ma1    -0.986804    0.000057 -1.7222e+04 0.000000
## omega   0.000005    0.000011  4.8662e-01 0.626531
## alpha1  0.099148    0.095923  1.0336e+00 0.301317
## beta1   0.826372    0.055975  1.4763e+01 0.000000
## gamma1  0.094860    0.095281  9.9559e-01 0.319452
## shape   4.526164    0.780495  5.7991e+00 0.000000
## 
## LogLikelihood : 7553.457 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2351
## Bayes        -6.2135
## Shibata      -6.2351
## Hannan-Quinn -6.2273
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       5.061 2.447e-02
## Lag[2*(p+q)+(p+q)-1][8]      8.946 4.660e-09
## Lag[4*(p+q)+(p+q)-1][14]    12.691 1.540e-02
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.003579  0.9523
## Lag[2*(p+q)+(p+q)-1][5]  0.242857  0.9892
## Lag[4*(p+q)+(p+q)-1][9]  0.864121  0.9912
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2642 0.500 2.000  0.6072
## ARCH Lag[5]    0.3180 1.440 1.667  0.9354
## ARCH Lag[7]    0.8171 2.315 1.543  0.9414
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.3312
## Individual Statistics:              
## mu     0.27433
## ar1    0.59877
## ar2    0.61429
## ma1    0.06929
## omega  0.21719
## alpha1 0.51162
## beta1  0.44052
## gamma1 0.44778
## shape  0.16134
## 
## 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            0.883 0.3773    
## Negative Sign Bias   1.176 0.2396    
## Positive Sign Bias   1.567 0.1172    
## Joint Effect         4.254 0.2353    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     40.58    0.0027453
## 2    30     57.01    0.0014334
## 3    40     65.12    0.0054260
## 4    50     86.78    0.0007115
## 
## 
## Elapsed time : 1.739141

# Phân phối Student đối xứng (sstd)
HNXst1.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(2, 1), include.mean = TRUE), distribution.model = "sstd")

HNXst2<- ugarchfit(HNXst1.spec,HNXts) 
print(HNXst2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.003982    0.001655      2.4061 0.016123
## ar1     1.028769    0.000920   1118.2078 0.000000
## ar2    -0.029247    0.000468    -62.5303 0.000000
## ma1    -0.988118    0.000022 -45281.6118 0.000000
## omega   0.000005    0.000003      1.9537 0.050737
## alpha1  0.097588    0.027058      3.6067 0.000310
## beta1   0.838084    0.018526     45.2375 0.000000
## gamma1  0.076212    0.032173      2.3688 0.017846
## skew    0.838197    0.023999     34.9263 0.000000
## shape   4.879363    0.496148      9.8345 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.003982    0.001894  2.1026e+00 0.035505
## ar1     1.028769    0.000924  1.1128e+03 0.000000
## ar2    -0.029247    0.000447 -6.5439e+01 0.000000
## ma1    -0.988118    0.000017 -5.8867e+04 0.000000
## omega   0.000005    0.000007  7.0261e-01 0.482297
## alpha1  0.097588    0.058644  1.6641e+00 0.096095
## beta1   0.838084    0.032632  2.5683e+01 0.000000
## gamma1  0.076212    0.063184  1.2062e+00 0.227744
## skew    0.838197    0.023484  3.5692e+01 0.000000
## shape   4.879363    0.530993  9.1891e+00 0.000000
## 
## LogLikelihood : 7572.526 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2500
## Bayes        -6.2261
## Shibata      -6.2501
## Hannan-Quinn -6.2413
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       7.333 0.0067699
## Lag[2*(p+q)+(p+q)-1][8]     12.315 0.0000000
## Lag[4*(p+q)+(p+q)-1][14]    16.204 0.0005566
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01043  0.9186
## Lag[2*(p+q)+(p+q)-1][5]   0.30819  0.9829
## Lag[4*(p+q)+(p+q)-1][9]   0.94441  0.9884
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2154 0.500 2.000  0.6425
## ARCH Lag[5]    0.2876 1.440 1.667  0.9435
## ARCH Lag[7]    0.7973 2.315 1.543  0.9441
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.1224
## Individual Statistics:              
## mu     0.04145
## ar1    0.49236
## ar2    0.50127
## ma1    0.07183
## omega  0.20119
## alpha1 0.57859
## beta1  0.51709
## gamma1 0.41417
## skew   0.08132
## shape  0.15396
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.29 2.54 3.05
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.2896 0.7721    
## Negative Sign Bias  1.2164 0.2239    
## Positive Sign Bias  1.2825 0.1998    
## Joint Effect        4.7646 0.1899    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     15.92       0.6628
## 2    30     25.00       0.6782
## 3    40     44.86       0.2396
## 4    50     43.10       0.7101
## 
## 
## Elapsed time : 2.120454

# Phân phối Generalized Error Distribution(ged)
HNXged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(2, 1), include.mean = TRUE), distribution.model = "ged")

HNXged1 <- ugarchfit(HNXged.spec,HNXts) 
print(HNXged1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : ged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.004386    0.000680     6.4491 0.000000
## ar1     1.022795    0.000841  1216.4387 0.000000
## ar2    -0.025098    0.001236   -20.3039 0.000000
## ma1    -0.985513    0.000147 -6694.2816 0.000000
## omega   0.000005    0.000006     0.8856 0.375834
## alpha1  0.109565    0.036058     3.0386 0.002377
## beta1   0.817235    0.029767    27.4547 0.000000
## gamma1  0.086624    0.041732     2.0757 0.037921
## shape   1.195800    0.033482    35.7144 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.004386    0.000317    13.82844  0.00000
## ar1     1.022795    0.005698   179.51039  0.00000
## ar2    -0.025098    0.003595    -6.98163  0.00000
## ma1    -0.985513    0.000115 -8563.25123  0.00000
## omega   0.000005    0.000034     0.15506  0.87677
## alpha1  0.109565    0.239309     0.45784  0.64707
## beta1   0.817235    0.195833     4.17312  0.00003
## gamma1  0.086624    0.288272     0.30049  0.76380
## shape   1.195800    0.171349     6.97872  0.00000
## 
## LogLikelihood : 7541.117 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2249
## Bayes        -6.2034
## Shibata      -6.2249
## Hannan-Quinn -6.2171
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       7.651 5.674e-03
## Lag[2*(p+q)+(p+q)-1][8]     11.473 1.110e-16
## Lag[4*(p+q)+(p+q)-1][14]    15.183 1.552e-03
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.05499  0.8146
## Lag[2*(p+q)+(p+q)-1][5]   0.30314  0.9834
## Lag[4*(p+q)+(p+q)-1][9]   1.00534  0.9860
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.3901 0.500 2.000  0.5322
## ARCH Lag[5]    0.4170 1.440 1.667  0.9078
## ARCH Lag[7]    0.9959 2.315 1.543  0.9142
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.356
## Individual Statistics:              
## mu     0.22956
## ar1    0.49405
## ar2    0.50180
## ma1    0.03618
## omega  0.23366
## alpha1 0.47005
## beta1  0.40022
## gamma1 0.41389
## shape  0.06803
## 
## 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           0.9206 0.35736    
## Negative Sign Bias  1.0319 0.30222    
## Positive Sign Bias  1.7802 0.07517   *
## Joint Effect        4.6636 0.19815    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     59.80    4.160e-06
## 2    30     82.17    5.544e-07
## 3    40     97.72    6.054e-07
## 4    50    100.62    1.996e-05
## 
## 
## Elapsed time : 2.16944

# Phân phối Generalized Error Distribution đối xứng ("sged")
HNXged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(2, 1), include.mean = TRUE), distribution.model = "sged")

HNXged2 <- ugarchfit(HNXged.spec,HNXts) 
print(HNXged2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(2,0,1)
## Distribution : sged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.002512    0.001510  1.6637e+00 0.096182
## ar1     1.002393    0.003497  2.8664e+02 0.000000
## ar2    -0.003231    0.001754 -1.8418e+00 0.065500
## ma1    -0.989055    0.000016 -6.0223e+04 0.000000
## omega   0.000005    0.000010  4.9118e-01 0.623299
## alpha1  0.105777    0.075584  1.3995e+00 0.161673
## beta1   0.832424    0.045685  1.8221e+01 0.000000
## gamma1  0.064690    0.088781  7.2865e-01 0.466213
## skew    0.836304    0.020728  4.0347e+01 0.000000
## shape   1.249604    0.054034  2.3126e+01 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.002512    0.013597  1.8478e-01 0.853402
## ar1     1.002393    0.035224  2.8458e+01 0.000000
## ar2    -0.003231    0.016679 -1.9373e-01 0.846384
## ma1    -0.989055    0.000052 -1.9100e+04 0.000000
## omega   0.000005    0.000108  4.3950e-02 0.964944
## alpha1  0.105777    0.818265  1.2927e-01 0.897144
## beta1   0.832424    0.482040  1.7269e+00 0.084190
## gamma1  0.064690    0.953524  6.7843e-02 0.945910
## skew    0.836304    0.089403  9.3544e+00 0.000000
## shape   1.249604    0.305509  4.0902e+00 0.000043
## 
## LogLikelihood : 7563.64 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2427
## Bayes        -6.2187
## Shibata      -6.2427
## Hannan-Quinn -6.2340
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       15.29 9.219e-05
## Lag[2*(p+q)+(p+q)-1][8]      21.39 0.000e+00
## Lag[4*(p+q)+(p+q)-1][14]     25.41 1.256e-08
## d.o.f=3
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0000297  0.9957
## Lag[2*(p+q)+(p+q)-1][5] 0.2263125  0.9905
## Lag[4*(p+q)+(p+q)-1][9] 0.9031302  0.9899
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.2755 0.500 2.000  0.5997
## ARCH Lag[5]    0.3099 1.440 1.667  0.9376
## ARCH Lag[7]    0.8835 2.315 1.543  0.9317
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  2.1986
## Individual Statistics:              
## mu     0.07972
## ar1    0.40272
## ar2    0.41043
## ma1    0.05368
## omega  0.15038
## alpha1 0.52469
## beta1  0.44476
## gamma1 0.37598
## skew   0.15637
## shape  0.07454
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.29 2.54 3.05
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.08975 0.9285    
## Negative Sign Bias 1.03758 0.2996    
## Positive Sign Bias 1.37261 0.1700    
## Joint Effect       5.24611 0.1546    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     21.79      0.29508
## 2    30     32.64      0.29265
## 3    40     44.60      0.24818
## 4    50     63.88      0.07501
## 
## 
## Elapsed time : 3.987283

1.7.2 LỰA CHỌN MÔ HÌNH TỐI ƯU CHO BIẾN HNX

#Tạo danh sách mô hình lựa chọn
HNX.model.list <- list(
  garch11n = HNXfit,
  garch11t = HNXst1,
  garch11st =HNXst2,
  garch11g = HNXged1,
  garch11sg = HNXged2
)
#Tính toán lựa chọn mô hình
HNX.info.mat <- sapply(HNX.model.list, infocriteria)
rownames(HNX.info.mat) <- rownames(infocriteria(HNXfit))
HNX.info.mat

##               garch11n  garch11t garch11st  garch11g garch11sg
## Akaike       -6.138273 -6.235089 -6.250022 -6.224890 -6.242678
## Bayes        -6.119128 -6.213550 -6.226090 -6.203351 -6.218746
## Shibata      -6.138295 -6.235116 -6.250055 -6.224918 -6.242712
## Hannan-Quinn -6.131311 -6.227256 -6.241319 -6.217058 -6.233975

1.8 TRÍCH CHUỖI LỢI NHUẬN CHO MÔ HÌNH KIỂM ĐỊNH BIÊN CỦA BIẾN HNX

# Trích xuất chuỗi phần dư v của chuỗi lợi suất HNX
HNX.res <- residuals(HNXst2)/sigma(HNXst2)
fitdist(distribution = "sstd", HNX.res, control = list())

## $pars
##          mu       sigma        skew       shape 
## -0.01325317  1.00244760  0.83176853  4.89568261 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 3564.440 3286.508 3286.508
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]     [,2]        [,3]       [,4]
## [1,]  3438.22380  879.103 -1108.92060  31.466855
## [2,]   879.10303 3095.334   191.44304 118.715972
## [3,] -1108.92060  191.443  2288.57404  11.453314
## [4,]    31.46685  118.716    11.45331   8.684299
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 122
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.217237 secs
## 
## $vscale
## [1] 1 1 1 1 1

HNX.res

##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 0001-01-01                                 -0.25680144
## 0002-01-01                                  0.63920032
## 0003-01-01                                  0.39562074
## 0004-01-01                                  0.25910014
## 0005-01-01                                  1.10030522
## 0006-01-01                                 -0.94033093
## 0007-01-01                                  0.00734706
## 0008-01-01                                  0.44375807
## 0009-01-01                                  0.94990269
## 0010-01-01                                  0.03370242
##        ...                                            
## 2411-01-01                                 -0.07545673
## 2412-01-01                                 -0.59054231
## 2413-01-01                                  0.75245480
## 2414-01-01                                  0.40593442
## 2415-01-01                                  0.14164373
## 2416-01-01                                 -0.12426240
## 2417-01-01                                  1.63998099
## 2418-01-01                                 -0.40115354
## 2419-01-01                                  0.38060994
## 2420-01-01                                 -0.20625148

s01 = pdist("sstd",HNX.res, mu =-0.01325317    
      
, sigma =  1.00244760, skew =  0.83176853 , shape = 4.89568261  )
head(s01,10)

##  [1] 0.3529759 0.7709597 0.6625621 0.5952479 0.9054319 0.1439761 0.4714384
##  [8] 0.6855289 0.8720351 0.4840772

1.8.1 BIẾN S&P 500

library(rugarch)
SPts<- ts(mhnn$SP)
head(SPts,10)

## Time Series:
## Start = 1 
## End = 10 
## Frequency = 1 
##  [1] -0.0003275646 -0.0025148994  0.0060578132 -0.0002176634  0.0003264773
##  [6]  0.0023366401 -0.0126722255  0.0107707228  0.0051528332 -0.0013534366

library(lmtest)
SPspec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1,1)), mean.model = list(armaOrder = c(5,5), include.mean =TRUE), distribution.model = 'norm')
print(SPspec)

## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : gjrGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(5,0,5)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE

SPfit <- ugarchfit(spec = SPspec, SPts)
print(SPfit)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(5,0,5)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.000430    0.000135  3.1818e+00 0.001464
## ar1    -0.242172    0.018275 -1.3251e+01 0.000000
## ar2    -0.004167    0.017817 -2.3385e-01 0.815100
## ar3    -0.431464    0.014085 -3.0633e+01 0.000000
## ar4     0.299702    0.018360  1.6324e+01 0.000000
## ar5     0.597360    0.015527  3.8472e+01 0.000000
## ma1     0.190603    0.001315  1.4491e+02 0.000000
## ma2     0.045483    0.000427  1.0647e+02 0.000000
## ma3     0.433967    0.001155  3.7589e+02 0.000000
## ma4    -0.319816    0.001059 -3.0189e+02 0.000000
## ma5    -0.558203    0.000044 -1.2691e+04 0.000000
## omega   0.000004    0.000000  7.4721e+00 0.000000
## alpha1  0.037375    0.008048  4.6440e+00 0.000003
## beta1   0.797203    0.013671  5.8313e+01 0.000000
## gamma1  0.274259    0.032233  8.5086e+00 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu      0.000430    0.000212     2.03132 0.042223
## ar1    -0.242172    0.018737   -12.92465 0.000000
## ar2    -0.004167    0.020766    -0.20064 0.840977
## ar3    -0.431464    0.017683   -24.39929 0.000000
## ar4     0.299702    0.021878    13.69903 0.000000
## ar5     0.597360    0.017699    33.75035 0.000000
## ma1     0.190603    0.002581    73.84063 0.000000
## ma2     0.045483    0.003055    14.88952 0.000000
## ma3     0.433967    0.000962   451.15736 0.000000
## ma4    -0.319816    0.003213   -99.54129 0.000000
## ma5    -0.558203    0.000061 -9212.98784 0.000000
## omega   0.000004    0.000001     2.95785 0.003098
## alpha1  0.037375    0.027848     1.34211 0.179561
## beta1   0.797203    0.016481    48.37063 0.000000
## gamma1  0.274259    0.063573     4.31409 0.000016
## 
## LogLikelihood : 8104.44 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.6855
## Bayes        -6.6496
## Shibata      -6.6856
## Hannan-Quinn -6.6724
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.7752  0.3786
## Lag[2*(p+q)+(p+q)-1][29]   13.3378  0.9983
## Lag[4*(p+q)+(p+q)-1][49]   21.1384  0.8539
## d.o.f=10
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.4098  0.5221
## Lag[2*(p+q)+(p+q)-1][5]    3.3198  0.3517
## Lag[4*(p+q)+(p+q)-1][9]    6.6265  0.2325
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.7866 0.500 2.000 0.37512
## ARCH Lag[5]    5.9167 1.440 1.667 0.06284
## ARCH Lag[7]    7.2107 2.315 1.543 0.07820
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  9.6398
## Individual Statistics:              
## mu     0.60056
## ar1    0.17328
## ar2    0.10198
## ar3    0.24179
## ar4    0.01637
## ar5    0.15427
## ma1    0.10391
## ma2    0.08301
## ma3    0.16353
## ma4    0.02271
## ma5    0.12733
## omega  0.43304
## alpha1 0.68219
## beta1  1.06403
## gamma1 0.65318
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.26 3.54 4.07
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.8697 0.004144 ***
## Negative Sign Bias  1.4643 0.143229    
## Positive Sign Bias  0.1616 0.871647    
## Joint Effect       11.6338 0.008749 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     90.45    2.765e-11
## 2    30    105.63    1.206e-10
## 3    40    126.91    3.095e-11
## 4    50    145.70    1.560e-11
## 
## 
## Elapsed time : 2.498694

SPst.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(5, 5), include.mean = TRUE), distribution.model = "std")

SPst1<- ugarchfit(SPst.spec,SPts) 
print(SPst1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(5,0,5)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000650    0.000065  10.008708 0.000000
## ar1    -0.246391    0.063687  -3.868768 0.000109
## ar2     0.007840    0.053711   0.145963 0.883951
## ar3     0.001240    0.073854   0.016796 0.986599
## ar4     0.237122    0.111961   2.117903 0.034183
## ar5     0.860141    0.096775   8.888075 0.000000
## ma1     0.221619    0.051366   4.314521 0.000016
## ma2    -0.003144    0.053297  -0.058998 0.952954
## ma3    -0.011450    0.067127  -0.170577 0.864556
## ma4    -0.237273    0.104417  -2.272363 0.023065
## ma5    -0.884834    0.076880 -11.509251 0.000000
## omega   0.000003    0.000003   0.983720 0.325253
## alpha1  0.014666    0.037130   0.394987 0.692853
## beta1   0.814944    0.004057 200.870275 0.000000
## gamma1  0.296378    0.059128   5.012470 0.000001
## shape   5.959861    0.605796   9.838071 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000650    0.001071  0.606912 0.543909
## ar1    -0.246391    0.241622 -1.019737 0.307853
## ar2     0.007840    0.079607  0.098481 0.921551
## ar3     0.001240    0.329686  0.003763 0.996998
## ar4     0.237122    0.340990  0.695392 0.486809
## ar5     0.860141    0.410464  2.095536 0.036123
## ma1     0.221619    0.176283  1.257178 0.208689
## ma2    -0.003144    0.106416 -0.029548 0.976428
## ma3    -0.011450    0.303339 -0.037748 0.969889
## ma4    -0.237273    0.314726 -0.753904 0.450907
## ma5    -0.884834    0.365515 -2.420791 0.015487
## omega   0.000003    0.000024  0.115963 0.907682
## alpha1  0.014666    0.356977  0.041084 0.967229
## beta1   0.814944    0.133748  6.093114 0.000000
## gamma1  0.296378    0.617350  0.480082 0.631169
## shape   5.959861    3.111056  1.915704 0.055403
## 
## LogLikelihood : 8175.866 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.7437
## Bayes        -6.7054
## Shibata      -6.7438
## Hannan-Quinn -6.7298
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.2282  0.6329
## Lag[2*(p+q)+(p+q)-1][29]   12.0371  1.0000
## Lag[4*(p+q)+(p+q)-1][49]   19.8043  0.9287
## d.o.f=10
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1163  0.7331
## Lag[2*(p+q)+(p+q)-1][5]    1.6468  0.7038
## Lag[4*(p+q)+(p+q)-1][9]    3.9968  0.5894
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1326 0.500 2.000  0.7157
## ARCH Lag[5]    3.0931 1.440 1.667  0.2766
## ARCH Lag[7]    4.1039 2.315 1.543  0.3318
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  5.8107
## Individual Statistics:              
## mu     1.34898
## ar1    0.02209
## ar2    0.10292
## ar3    0.01735
## ar4    0.10700
## ar5    0.09737
## ma1    0.03192
## ma2    0.11273
## ma3    0.01909
## ma4    0.06227
## ma5    0.12022
## omega  3.49546
## alpha1 1.72109
## beta1  2.17515
## gamma1 1.49973
## shape  1.31101
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.46 3.75 4.3
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                     t-value     prob sig
## Sign Bias           2.83283 0.004652 ***
## Negative Sign Bias  1.88309 0.059809   *
## Positive Sign Bias  0.02461 0.980367    
## Joint Effect       10.64863 0.013786  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     49.79    1.410e-04
## 2    30     73.15    1.111e-05
## 3    40     81.02    8.993e-05
## 4    50    106.78    3.549e-06
## 
## 
## Elapsed time : 3.393515

# Phân phối Student đối xứng (sstd)
SPst1.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(5, 5), include.mean = TRUE), distribution.model = "sstd")

SPst2<- ugarchfit(SPst1.spec,SPts) 
print(SPst2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(5,0,5)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000483    0.000142     3.3999 0.000674
## ar1    -0.792969    0.022697   -34.9371 0.000000
## ar2    -0.236382    0.032566    -7.2586 0.000000
## ar3    -0.712699    0.027852   -25.5886 0.000000
## ar4    -0.103392    0.041612    -2.4847 0.012967
## ar5     0.601040    0.023189    25.9196 0.000000
## ma1     0.741171    0.001075   689.5274 0.000000
## ma2     0.178287    0.010375    17.1839 0.000000
## ma3     0.671166    0.000144  4665.7053 0.000000
## ma4     0.033318    0.009590     3.4744 0.000512
## ma5    -0.644882    0.000518 -1244.1045 0.000000
## omega   0.000002    0.000002     1.2281 0.219411
## alpha1  0.026449    0.020825     1.2701 0.204063
## beta1   0.830222    0.023261    35.6915 0.000000
## gamma1  0.254339    0.063293     4.0184 0.000059
## skew    0.820993    0.029754    27.5925 0.000000
## shape   6.991662    0.938353     7.4510 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000483    0.000473    1.02166 0.306940
## ar1    -0.792969    0.119874   -6.61503 0.000000
## ar2    -0.236382    0.163722   -1.44380 0.148796
## ar3    -0.712699    0.146917   -4.85102 0.000001
## ar4    -0.103392    0.220448   -0.46901 0.639065
## ar5     0.601040    0.123251    4.87655 0.000001
## ma1     0.741171    0.005603  132.27179 0.000000
## ma2     0.178287    0.015659   11.38585 0.000000
## ma3     0.671166    0.002209  303.79036 0.000000
## ma4     0.033318    0.044319    0.75179 0.452179
## ma5    -0.644882    0.001638 -393.69511 0.000000
## omega   0.000002    0.000013    0.18425 0.853815
## alpha1  0.026449    0.093254    0.28363 0.776697
## beta1   0.830222    0.118484    7.00702 0.000000
## gamma1  0.254339    0.349564    0.72759 0.466867
## skew    0.820993    0.087011    9.43553 0.000000
## shape   6.991662    1.053917    6.63398 0.000000
## 
## LogLikelihood : 8194.278 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.7581
## Bayes        -6.7174
## Shibata      -6.7582
## Hannan-Quinn -6.7433
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.081  0.2986
## Lag[2*(p+q)+(p+q)-1][29]    14.622  0.7376
## Lag[4*(p+q)+(p+q)-1][49]    22.356  0.7555
## d.o.f=10
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.03904  0.8434
## Lag[2*(p+q)+(p+q)-1][5]   2.15623  0.5820
## Lag[4*(p+q)+(p+q)-1][9]   4.54420  0.4993
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.6584 0.500 2.000  0.4171
## ARCH Lag[5]    3.9188 1.440 1.667  0.1814
## ARCH Lag[7]    4.8361 2.315 1.543  0.2416
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  14.4228
## Individual Statistics:              
## mu     1.04312
## ar1    0.07470
## ar2    0.05699
## ar3    0.04042
## ar4    0.02611
## ar5    0.03707
## ma1    0.03910
## ma2    0.06551
## ma3    0.02707
## ma4    0.04894
## ma5    0.07728
## omega  7.80071
## alpha1 1.68267
## beta1  2.01136
## gamma1 1.39095
## skew   0.29193
## shape  0.99897
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.64 3.95 4.51
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.9035 0.003724 ***
## Negative Sign Bias  1.6468 0.099735   *
## Positive Sign Bias  0.4215 0.673433    
## Joint Effect       12.9605 0.004723 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     26.07      0.12836
## 2    30     49.99      0.00905
## 3    40     50.78      0.09806
## 4    50     71.74      0.01878
## 
## 
## Elapsed time : 6.0145

# Phân phối Generalized Error Distribution(ged)
SPged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(5, 5), include.mean = TRUE), distribution.model = "ged")

SPged1 <- ugarchfit(SPged.spec,SPts) 
print(SPged1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(5,0,5)
## Distribution : ged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000600    0.000074    8.1231 0.000000
## ar1    -0.914328    0.005952 -153.6236 0.000000
## ar2     0.525931    0.010244   51.3385 0.000000
## ar3     0.697135    0.015598   44.6949 0.000000
## ar4     0.286209    0.018364   15.5853 0.000000
## ar5     0.293406    0.018622   15.7557 0.000000
## ma1     0.861033    0.004389  196.1778 0.000000
## ma2    -0.549901    0.006755  -81.4121 0.000000
## ma3    -0.643758    0.006251 -102.9848 0.000000
## ma4    -0.278068    0.009392  -29.6074 0.000000
## ma5    -0.319431    0.002206 -144.8266 0.000000
## omega   0.000003    0.000001    3.7642 0.000167
## alpha1  0.023897    0.002388   10.0093 0.000000
## beta1   0.804327    0.015263   52.6976 0.000000
## gamma1  0.289071    0.032846    8.8009 0.000000
## shape   1.325046    0.049646   26.6901 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000600    0.000090   6.70233 0.000000
## ar1    -0.914328    0.014212 -64.33676 0.000000
## ar2     0.525931    0.005211 100.93534 0.000000
## ar3     0.697135    0.022481  31.00959 0.000000
## ar4     0.286209    0.019686  14.53853 0.000000
## ar5     0.293406    0.024449  12.00065 0.000000
## ma1     0.861033    0.002189 393.41283 0.000000
## ma2    -0.549901    0.015506 -35.46367 0.000000
## ma3    -0.643758    0.018170 -35.42984 0.000000
## ma4    -0.278068    0.009815 -28.33145 0.000000
## ma5    -0.319431    0.024047 -13.28386 0.000000
## omega   0.000003    0.000002   1.50648 0.131943
## alpha1  0.023897    0.038604   0.61904 0.535889
## beta1   0.804327    0.019296  41.68293 0.000000
## gamma1  0.289071    0.067532   4.28053 0.000019
## shape   1.325046    0.065828  20.12888 0.000000
## 
## LogLikelihood : 8165.074 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.7348
## Bayes        -6.6965
## Shibata      -6.7349
## Hannan-Quinn -6.7208
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.7203  0.3960
## Lag[2*(p+q)+(p+q)-1][29]   12.2880  1.0000
## Lag[4*(p+q)+(p+q)-1][49]   20.0781  0.9162
## d.o.f=10
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1092  0.7411
## Lag[2*(p+q)+(p+q)-1][5]    1.9318  0.6346
## Lag[4*(p+q)+(p+q)-1][9]    4.1490  0.5638
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     0.502 0.500 2.000  0.4786
## ARCH Lag[5]     3.135 1.440 1.667  0.2707
## ARCH Lag[7]     4.039 2.315 1.543  0.3410
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  5.3617
## Individual Statistics:              
## mu     1.08000
## ar1    0.04679
## ar2    0.04462
## ar3    0.10448
## ar4    0.15988
## ar5    0.25989
## ma1    0.05195
## ma2    0.06785
## ma3    0.12831
## ma4    0.23102
## ma5    0.31626
## omega  0.50014
## alpha1 1.31227
## beta1  1.79044
## gamma1 1.25189
## shape  0.75364
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.46 3.75 4.3
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           2.7152 0.006671 ***
## Negative Sign Bias  1.6485 0.099389   *
## Positive Sign Bias  0.2332 0.815616    
## Joint Effect       10.6783 0.013599  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     58.78    6.032e-06
## 2    30     73.87    8.806e-06
## 3    40     82.88    5.326e-05
## 4    50     87.73    5.653e-04
## 
## 
## Elapsed time : 7.487664

# Phân phối Generalized Error Distribution đối xứng ("sged")
SPged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(5, 5), include.mean = TRUE), distribution.model = "sged")

SPged2 <- ugarchfit(SPged.spec,SPts) 
print(SPged2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(5,0,5)
## Distribution : sged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.000385    0.000087    4.43284 0.000009
## ar1    -0.803468    0.009250  -86.85753 0.000000
## ar2    -0.229046    0.005003  -45.78145 0.000000
## ar3    -0.668655    0.003848 -173.74445 0.000000
## ar4    -0.073686    0.003989  -18.47275 0.000000
## ar5     0.587097    0.008773   66.92301 0.000000
## ma1     0.745114    0.009443   78.90761 0.000000
## ma2     0.164712    0.011484   14.34241 0.000000
## ma3     0.634208    0.012656   50.11018 0.000000
## ma4     0.007851    0.010076    0.77916 0.435887
## ma5    -0.628531    0.012122  -51.85175 0.000000
## omega   0.000003    0.000002    1.55023 0.121087
## alpha1  0.030118    0.017925    1.68018 0.092923
## beta1   0.823979    0.013937   59.12112 0.000000
## gamma1  0.252694    0.017662   14.30687 0.000000
## skew    0.834577    0.022154   37.67119 0.000000
## shape   1.416800    0.053606   26.42976 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000385    0.000475   0.80983  0.41803
## ar1    -0.803468    0.013210 -60.82285  0.00000
## ar2    -0.229046    0.018591 -12.32035  0.00000
## ar3    -0.668655    0.028468 -23.48774  0.00000
## ar4    -0.073686    0.061447  -1.19919  0.23045
## ar5     0.587097    0.031113  18.87009  0.00000
## ma1     0.745114    0.013909  53.57121  0.00000
## ma2     0.164712    0.027993   5.88413  0.00000
## ma3     0.634208    0.047284  13.41286  0.00000
## ma4     0.007851    0.075958   0.10336  0.91768
## ma5    -0.628531    0.042468 -14.80028  0.00000
## omega   0.000003    0.000010   0.25694  0.79722
## alpha1  0.030118    0.141293   0.21316  0.83120
## beta1   0.823979    0.044130  18.67153  0.00000
## gamma1  0.252694    0.232945   1.08478  0.27802
## skew    0.834577    0.028863  28.91488  0.00000
## shape   1.416800    0.113748  12.45565  0.00000
## 
## LogLikelihood : 8182.412 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.7483
## Bayes        -6.7076
## Shibata      -6.7484
## Hannan-Quinn -6.7335
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       1.784  0.1816
## Lag[2*(p+q)+(p+q)-1][29]    14.461  0.8203
## Lag[4*(p+q)+(p+q)-1][49]    21.574  0.8219
## d.o.f=10
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1277  0.7208
## Lag[2*(p+q)+(p+q)-1][5]    2.0984  0.5954
## Lag[4*(p+q)+(p+q)-1][9]    4.4244  0.5185
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.7228 0.500 2.000  0.3952
## ARCH Lag[5]    3.5818 1.440 1.667  0.2158
## ARCH Lag[7]    4.5388 2.315 1.543  0.2755
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  4.7869
## Individual Statistics:              
## mu     1.04613
## ar1    0.07399
## ar2    0.04847
## ar3    0.06077
## ar4    0.02946
## ar5    0.11173
## ma1    0.06273
## ma2    0.05584
## ma3    0.11032
## ma4    0.04733
## ma5    0.13139
## omega  1.74734
## alpha1 1.33876
## beta1  1.76027
## gamma1 1.26039
## skew   0.31069
## shape  0.66434
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.64 3.95 4.51
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias            3.277 0.001063 ***
## Negative Sign Bias   1.778 0.075548   *
## Positive Sign Bias   0.393 0.694384    
## Joint Effect        15.964 0.001153 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     38.07     0.005821
## 2    30     47.98     0.014756
## 3    40     65.32     0.005182
## 4    50     76.57     0.007111
## 
## 
## Elapsed time : 8.674628

1.8.2 LỰA CHỌN MÔ HÌNH TỐI ƯU CHO BIẾN S&P 500

#Tạo danh sách mô hình lựa chọn
SP.model.list <- list(
  garchn = SPfit,
  garcht = SPst1,
  garchst = SPst2,
  garchg = SPged1,
  garchsg = SPged2
)
#Tính toán lựa chọn mô hình
SP.info.mat <- sapply(SP.model.list, infocriteria)
rownames(SP.info.mat) <- rownames(infocriteria(SPfit))
SP.info.mat

##                 garchn    garcht   garchst    garchg   garchsg
## Akaike       -6.685488 -6.743691 -6.758081 -6.734772 -6.748275
## Bayes        -6.649590 -6.705400 -6.717397 -6.696481 -6.707590
## Shibata      -6.685564 -6.743777 -6.758179 -6.734859 -6.748372
## Hannan-Quinn -6.672433 -6.729766 -6.743286 -6.720847 -6.733480

1.9 TRÍCH SUẤT CHUỖI PHẦN DƯ CHO MÔ HÌNH KIỂM ĐỊNH BIÊN CỦA BIẾN S&P 500

# Trích xuất chuỗi phần dư v của chuỗi lợi suất S&P
SP.res <- residuals(SPst2)/sigma(SPst2)
fitdist(distribution = "sstd", SP.res, control = list())

## $pars
##           mu        sigma         skew        shape 
## -0.006756084  0.999503039  0.818894211  7.054790485 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 3631.927 3351.086 3351.086
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]       [,2]      [,3]      [,4]
## [1,] 3180.083442  793.36941 -716.4760  6.509737
## [2,]  793.369413 3337.82312  442.7873 39.985434
## [3,] -716.476046  442.78727 1805.9856 12.849401
## [4,]    6.509737   39.98543   12.8494  1.644984
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 131
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.2309201 secs
## 
## $vscale
## [1] 1 1 1 1 1

SP.res

##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 0001-01-01                                 -0.07211003
## 0002-01-01                                 -0.26669392
## 0003-01-01                                  0.49592905
## 0004-01-01                                 -0.06233329
## 0005-01-01                                 -0.01392689
## 0006-01-01                                  0.17290521
## 0007-01-01                                 -1.33741718
## 0008-01-01                                  0.86527841
## 0009-01-01                                  0.47429832
## 0010-01-01                                 -0.18068765
##        ...                                            
## 2411-01-01                                 -0.01857413
## 2412-01-01                                  0.81715495
## 2413-01-01                                  1.32899038
## 2414-01-01                                 -3.15740983
## 2415-01-01                                  1.08372380
## 2416-01-01                                  0.11441942
## 2417-01-01                                  0.49587936
## 2418-01-01                                  0.17303175
## 2419-01-01                                  0.07801875
## 2420-01-01                                 -0.50970764

s02 = pdist("sstd",SP.res, mu = -0.006756084     , sigma =  0.999503039 , skew = 0.818894211, shape =7.054790485  )
head(s02,10)

##  [1] 0.43617463 0.35652452 0.69246692 0.44037196 0.46138512 0.54530998
##  [7] 0.08826068 0.83223260 0.68296447 0.39075421

1.9.1 KIỂM ĐỊNH SỰ PHÙ HỢP MÔ HÌNH BIÊN

1.9.1.1 KIỂM ĐỊNH CHO BIẾN HNX

# Kiểm định Anderson_Darling
library(nortest)
ad.test(s01)

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

# Kiểm định Cramer-von Mises
cvm.test(s01)

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

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

# Kiểm định Kolmogorov-Smirnov
ks.test(s01, y = "punif")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  s01
## D = 0.013052, p-value = 0.8041
## alternative hypothesis: two-sided

1.9.1.2 KIỂM ĐỊNH CHO BIẾN S&P 500

# Kiểm định Anderson_Darling
library(nortest)
ad.test(s02)

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

# Kiểm định Cramer-von Mises
cvm.test(s02)

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

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

# Kiểm định Kolmogorov-Smirnov
ks.test(s02, y = "punif")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  s02
## D = 0.022151, p-value = 0.1859
## alternative hypothesis: two-sided

1.10 MÔ HÌNH COPULA

1.10.1 ƯỚC LƯỢNG HỌ MÔ HÌNH COPULA

library(VineCopula)
library("copula")

## 
## Attaching package: 'copula'

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

library("scatterplot3d")
#Chuyển đôi dữ liệu phân phối đều
u <- pobs(s01)
head(u,10)

##  [1] 0.3519207 0.7686906 0.6670797 0.6009913 0.9070632 0.1420900 0.4779017
##  [8] 0.6906237 0.8719537 0.4907063

k <- pobs(s02)
head(k,10)

##  [1] 0.44609665 0.35315985 0.69434118 0.45187939 0.47666254 0.56381660
##  [7] 0.09417596 0.82280050 0.68566708 0.39487815

# ƯỚC LƯỢNG THAM SỐ MÔ HÌNH COPULA
BiCopSelect(u,k, familyset= c(1:10,13,16:20), selectioncrit="AIC",indeptest = FALSE, level = 0.05)

## Bivariate copula: Survival Joe (par = 1.07, tau = 0.04)

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

library(copula)
# 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, k)
fit <- fitCopula(joeCopula(), data = tnt, method = "ml")  
# Maximum Likelihood Estimation (MLE)
theta <- coef(fit)
print(theta)

##    alpha 
## 1.027689

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

##     alpha 
## 0.0157639

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

1.10.3 ĐỒ THỊ PHÂN TÁN

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

plot(k, main = "Đồ thị phân tán của biến S&P500")

1.10.4 ĐỒ THỊ MẬT ĐỘ

# Tính toán hàm mật độ và vẽ đồ thị
library(MASS)
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")

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

1.10.6 ƯỚ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.09394604

print(coef(fit_gumbel))

##    alpha 
## 1.030923

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

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

# Create histograms of marginal distributions
p_hist_X1 <- ggplot(tnt, aes(x = s01)) +
  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 = s02)) +
  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(s01,s02, 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(s01,s02)),method="ml")
coef(fit)

##    alpha 
## 1.030332

a$par

## [1] 1.045522

a$par2

## [1] 0

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