CÁC MÔ HÌNH NGẪU NHIÊN

Nhập dữ liệu

library(PerformanceAnalytics)

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

## Loading required package: xts

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

## Loading required package: zoo

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

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'PerformanceAnalytics'

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

library(readxl)

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

library(xlsx)

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

data <- read_xlsx("D:/naaaaaa/MHNN/MSCI VNIndex.xlsx")
head(data,10)

## # A tibble: 10 × 3
##    Date                 MSCI   VNI
##    <dttm>              <dbl> <dbl>
##  1 2018-01-02 00:00:00 12768  996.
##  2 2018-01-03 00:00:00 12962 1006.
##  3 2018-01-04 00:00:00 13166 1020.
##  4 2018-01-05 00:00:00 13303 1013.
##  5 2018-01-08 00:00:00 13321 1023.
##  6 2018-01-09 00:00:00 13410 1034.
##  7 2018-01-10 00:00:00 13375 1038.
##  8 2018-01-11 00:00:00 13476 1048.
##  9 2018-01-12 00:00:00 13589 1050.
## 10 2018-01-16 00:00:00 13433 1063.

Chuyển đổi về tỷ suất lợi nhuận

MSCI <- diff(log(data$MSCI), lag = 1)
VNI <- diff(log(data$VNI), lag = 1)
mhnn <- data.frame( msci = MSCI, vni = VNI)

Tương quan

res <- cor(mhnn)
round(res, 2)

##      msci  vni
## msci 1.00 0.06
## vni  0.06 1.00

library(ggplot2)

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

ggplot(data.frame(price = mhnn$vni),  aes(sample = price)) +
  geom_qq() +
  geom_qq_line() +
  labs(title = "Q-Q Plot of VNI",
       x = "Theoretical Quantiles", 
       y = "Sample Quantiles") +  
   theme_bw() +
  theme_minimal()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  theme(aspect.ratio = 1, plot.title = element_text(hjust = 0.5))

library(ggplot2)
ggplot(data.frame(price = mhnn$msci),  aes(sample = price)) +
  geom_qq() +
  geom_qq_line() +
  labs(title = "Q-Q Plot of MSCI",
       x = "Theoretical Quantiles", 
       y = "Sample Quantiles") +  
     theme_bw() +
  theme_minimal()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+
  theme(aspect.ratio = 1, plot.title = element_text(hjust = 0.5))

cor(mhnn$vni,mhnn$msci, method = "pearson")

## [1] 0.06459518

cor(mhnn$vni,mhnn$msci, method = "spearman")

## [1] 0.03884516

cor(mhnn$vni,mhnn$msci, method = "kendall")

## [1] 0.0260835

library(corrplot)

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

## corrplot 0.92 loaded

res <- cor(mhnn)
rounded_res <- round(res, 2)
round(res, 2)

##      msci  vni
## msci 1.00 0.06
## vni  0.06 1.00

corrplot(rounded_res, method = "color")

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

library(ggplot2)
library(ggcorrplot)

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

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

# Tạo thống kê mô tả

library(tidyverse)

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

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

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first()  masks xts::first()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::last()   masks xts::last()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(psych)

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

## 
## Attaching package: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

summary <- summary(mhnn)
detailed <- describe(mhnn)
print(summary)

##       msci                vni            
##  Min.   :-0.144893   Min.   :-6.908e-02  
##  1st Qu.:-0.009900   1st Qu.:-4.802e-03  
##  Median : 0.001557   Median : 1.109e-03  
##  Mean   : 0.001025   Mean   : 8.705e-05  
##  3rd Qu.: 0.013023   3rd Qu.: 6.931e-03  
##  Max.   : 0.165811   Max.   : 4.860e-02

print(detailed)

##      vars    n mean   sd median trimmed  mad   min  max range  skew kurtosis se
## msci    1 1452    0 0.02      0       0 0.02 -0.14 0.17  0.31 -0.20     7.04  0
## vni     2 1452    0 0.01      0       0 0.01 -0.07 0.05  0.12 -0.94     3.64  0

plot.ts(mhnn$vni)

plot.ts(mhnn$msci)

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

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

# THỰC HÀNH CÁC KIỂM ĐỊNH

library(tseries)

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

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

library(xlsx)

kiểm định cho chỉ số MSCI

Var1 <- mhnn$msci

# Kiểm định phân phối chuẩn cho MSCI
result1 <-  jarque.bera.test(Var1)
print(result1)

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

Kiểm định hiệu ứng GARCH

library(rugarch)

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

## Loading required package: parallel

## 
## Attaching package: 'rugarch'

## The following object is masked from 'package:purrr':
## 
##     reduce

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

MSCIts <- ts(mhnn$msci)

kiểm định cho chỉ số VNI

Var2 <- mhnn$vni

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

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

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

Kiểm định tính dừng: Augmented Dickey–Fuller

Kiểm định tính dừng của MSCI

adf.test(Var1)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var1
## Dickey-Fuller = -10.738, Lag order = 11, 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 = -10.738, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

Kiểm định tính dừng của VNI

adf.test(Var2)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Var2
## Dickey-Fuller = -11.439, Lag order = 11, 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 = -11.439, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

kiểm định tự tương quan của MSCI

library(stats)
result3 <-  Box.test(Var1, lag = 2, type = "Ljung-Box")
print(result3)

## 
##  Box-Ljung test
## 
## data:  Var1
## X-squared = 24.594, df = 2, p-value = 4.565e-06

kiểm định tự tương quan của VNI

library(stats)
result4 <-  Box.test(Var2, lag = 2, type = "Ljung-Box")
print(result4)

## 
##  Box-Ljung test
## 
## data:  Var2
## X-squared = 5.7919, df = 2, p-value = 0.05525

Mô hình ARIMA

library(lmtest)

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

library(fGarch)

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

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

library(forecast)

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

modelD1 <- auto.arima(mhnn$msci)
modelD1

## Series: mhnn$msci 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1   mean
##       -0.1198  1e-03
## s.e.   0.0261  5e-04
## 
## sigma^2 = 0.000488:  log likelihood = 3476.53
## AIC=-6947.06   AICc=-6947.04   BIC=-6931.22

Trích xuất chuỗi phần dư của chuỗi lợi suất MSCI

library(forecast)
modelD<- auto.arima(mhnn$msci)
modelD

## Series: mhnn$msci 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1   mean
##       -0.1198  1e-03
## s.e.   0.0261  5e-04
## 
## sigma^2 = 0.000488:  log likelihood = 3476.53
## AIC=-6947.06   AICc=-6947.04   BIC=-6931.22

residuals <- residuals(modelD)
sigma <- sd(residuals)
msci.res <- residuals / sigma
fitdist(distribution = "sstd", msci.res, control = list())

## $pars
##           mu        sigma         skew        shape 
## -0.004740383  1.006250197  0.911740028  3.863497521 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2043.823 1929.298 1929.298
## 
## $lagrange
## [1] 0
## 
## $hessian
##            [,1]       [,2]        [,3]       [,4]
## [1,] 2145.21221  288.19039 -778.518835  20.558567
## [2,]  288.19039 1676.64173   40.967984 146.018156
## [3,] -778.51884   40.96798 1136.963504   3.338239
## [4,]   20.55857  146.01816    3.338239  18.930815
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 103
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.2223241 secs
## 
## $vscale
## [1] 1 1 1 1 1

d <- msci.res

Kiểm định hiệu ứng ARCH - LM cho chỉ số chứng khoán MSCI

arch_spec1 <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_MSCI <- ugarchfit(spec = arch_spec1, data = mhnn$msci)

residuals <- residuals(arch_MSCI)

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

arch1 <- bptest(arch_lm_model)

Hiển thị kết quả

arch1

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

library(aTSA)

## 
## Attaching package: 'aTSA'

## The following object is masked from 'package:forecast':
## 
##     forecast

## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test

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

aTSA::arch.test(arima(mhnn$msci,order = c(1,0,0)))

## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order  PQ p.value
## [1,]     4 263       0
## [2,]     8 534       0
## [3,]    12 635       0
## [4,]    16 673       0
## [5,]    20 684       0
## [6,]    24 692       0
## Lagrange-Multiplier test: 
##      order   LM p.value
## [1,]     4 1466       0
## [2,]     8  444       0
## [3,]    12  283       0
## [4,]    16  208       0
## [5,]    20  164       0
## [6,]    24  135       0

library(forecast)
modelD01 <- auto.arima(mhnn$vni)
modelD01

## Series: mhnn$vni 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1      ma1
##       0.7040  -0.6582
## s.e.  0.1845   0.1961
## 
## sigma^2 = 0.0001709:  log likelihood = 4238.31
## AIC=-8470.62   AICc=-8470.61   BIC=-8454.78

Kiểm định hiệu ứng ARCH: ARCH-LM

library(FinTS)

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

## 
## Attaching package: 'FinTS'

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

Kiểm định hiệu ứng ARCH - LM cho chỉ số chứng khoán VNI

arch_spec2 <- ugarchspec(variance.model = list(model = "sGARCH"))
arch_VNI <- ugarchfit(spec = arch_spec2, data = mhnn$vni)

residuals1 <- residuals(arch_VNI)

n <- length(residuals1)
x <- 1:n

Tạo mô hình tuyến tính

arch_lm_model1 <- lm(residuals1^2 ~ x)

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

arch2 <- bptest(arch_lm_model1)

Hiển thị kết quả

arch2

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

library(aTSA)
aTSA::arch.test(arima(mhnn$vni,order = c(1,0,1)))

## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order  PQ p.value
## [1,]     4 181       0
## [2,]     8 289       0
## [3,]    12 371       0
## [4,]    16 399       0
## [5,]    20 410       0
## [6,]    24 411       0
## Lagrange-Multiplier test: 
##      order   LM p.value
## [1,]     4 1034       0
## [2,]     8  452       0
## [3,]    12  286       0
## [4,]    16  207       0
## [5,]    20  163       0
## [6,]    24  131       0

## Trích xuất chuỗi phần dư của chuỗi lợi suất VNI

library(forecast)
modeld2<- auto.arima(mhnn$msci)
modeld2

## Series: mhnn$msci 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1   mean
##       -0.1198  1e-03
## s.e.   0.0261  5e-04
## 
## sigma^2 = 0.000488:  log likelihood = 3476.53
## AIC=-6947.06   AICc=-6947.04   BIC=-6931.22

residuals <- residuals(modeld2)
sigma <- sd(residuals)
vn.res <- residuals / sigma
fitdist(distribution = "sstd", vn.res, control = list())

## $pars
##           mu        sigma         skew        shape 
## -0.004740383  1.006250197  0.911740028  3.863497521 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2043.823 1929.298 1929.298
## 
## $lagrange
## [1] 0
## 
## $hessian
##            [,1]       [,2]        [,3]       [,4]
## [1,] 2145.21221  288.19039 -778.518835  20.558567
## [2,]  288.19039 1676.64173   40.967984 146.018156
## [3,] -778.51884   40.96798 1136.963504   3.338239
## [4,]   20.55857  146.01816    3.338239  18.930815
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 103
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.11941 secs
## 
## $vscale
## [1] 1 1 1 1 1

VNIts <- ts(mhnn$vni)

Ước lượng mô hình GARCH(1,1) theo phân phối chuẩn cho chỉ số chứng khoán VNI

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

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

Mô hình Grach(1,1) với phân phối Student’s t

VNIspec1 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                                mean.model = list(armaOrder = c(1, 1)),
                                distribution.model = "std")

VNIfit <- ugarchfit(spec = VNIspec1, data = VNI)

VNIfit <- ugarchfit(spec = VNIspec, VNIts)
print(VNIfit)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000059    0.000348  0.16912 0.865705
## ar1     0.665587    0.178519  3.72838 0.000193
## ma1    -0.592629    0.189085 -3.13420 0.001723
## omega   0.000009    0.000000 74.51158 0.000000
## alpha1  0.021004    0.007444  2.82156 0.004779
## beta1   0.841500    0.010338 81.39982 0.000000
## gamma1  0.143951    0.025904  5.55717 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000059    0.000361  0.16301 0.870510
## ar1     0.665587    0.186705  3.56490 0.000364
## ma1    -0.592629    0.190662 -3.10827 0.001882
## omega   0.000009    0.000000 44.23816 0.000000
## alpha1  0.021004    0.010766  1.95100 0.051057
## beta1   0.841500    0.018299 45.98559 0.000000
## gamma1  0.143951    0.040611  3.54466 0.000393
## 
## LogLikelihood : 4396.643 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.0463
## Bayes        -6.0209
## Shibata      -6.0464
## Hannan-Quinn -6.0368
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1982  0.6562
## Lag[2*(p+q)+(p+q)-1][5]    0.2283  1.0000
## Lag[4*(p+q)+(p+q)-1][9]    0.5791  1.0000
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.09492  0.7580
## Lag[2*(p+q)+(p+q)-1][5]   0.47759  0.9611
## Lag[4*(p+q)+(p+q)-1][9]   1.15770  0.9788
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.5091 0.500 2.000  0.4755
## ARCH Lag[5]    0.6713 1.440 1.667  0.8322
## ARCH Lag[7]    0.8416 2.315 1.543  0.9379
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  30.3367
## Individual Statistics:              
## mu     0.16944
## ar1    0.18449
## ma1    0.18840
## omega  5.63503
## alpha1 0.10173
## beta1  0.08032
## gamma1 0.05789
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.4812 0.6305    
## Negative Sign Bias  0.1016 0.9191    
## Positive Sign Bias  1.2987 0.1942    
## Joint Effect        4.4686 0.2151    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     116.0    6.207e-16
## 2    30     138.5    3.165e-16
## 3    40     154.2    1.274e-15
## 4    50     170.1    2.799e-15
## 
## 
## Elapsed time : 0.8583422

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

VNIst1<- ugarchfit(VNIst.spec,VNIts) 
print(VNIst1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001042    0.000273   3.8105 0.000139
## ar1     0.728204    0.193922   3.7551 0.000173
## ma1    -0.676872    0.207848  -3.2566 0.001128
## omega   0.000014    0.000001  16.5463 0.000000
## alpha1  0.024729    0.012229   2.0221 0.043164
## beta1   0.754007    0.023856  31.6062 0.000000
## gamma1  0.291132    0.056710   5.1337 0.000000
## shape   3.770174    0.329116  11.4555 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001042    0.000275   3.7953 0.000147
## ar1     0.728204    0.224779   3.2397 0.001197
## ma1    -0.676872    0.239818  -2.8224 0.004766
## omega   0.000014    0.000001  11.2544 0.000000
## alpha1  0.024729    0.020558   1.2029 0.229022
## beta1   0.754007    0.024884  30.3007 0.000000
## gamma1  0.291132    0.055795   5.2179 0.000000
## shape   3.770174    0.335222  11.2468 0.000000
## 
## LogLikelihood : 4503.07 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.1916
## Bayes        -6.1625
## Shibata      -6.1916
## Hannan-Quinn -6.1807
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.07435  0.7851
## Lag[2*(p+q)+(p+q)-1][5]   0.13838  1.0000
## Lag[4*(p+q)+(p+q)-1][9]   0.45892  1.0000
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.154  0.2828
## Lag[2*(p+q)+(p+q)-1][5]     2.181  0.5763
## Lag[4*(p+q)+(p+q)-1][9]     2.976  0.7630
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.140 0.500 2.000  0.2856
## ARCH Lag[5]     1.343 1.440 1.667  0.6344
## ARCH Lag[7]     1.453 2.315 1.543  0.8307
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  33.0803
## Individual Statistics:             
## mu     0.2411
## ar1    0.2369
## ma1    0.2261
## omega  6.1196
## alpha1 0.2173
## beta1  0.2134
## gamma1 0.1798
## shape  0.2129
## 
## 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.1668 0.86758    
## Negative Sign Bias  1.0370 0.29989    
## Positive Sign Bias  1.8019 0.07176   *
## Joint Effect        4.5142 0.21103    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     46.82    3.795e-04
## 2    30     66.84    8.094e-05
## 3    40     75.60    3.963e-04
## 4    50     91.53    2.201e-04
## 
## 
## Elapsed time : 1.096038

Phân phối Student đối xứng (sstd)

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

VNIst2<- ugarchfit(VNIst1.spec,VNIts) 
print(VNIst2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000225    0.000384   0.58568 0.558089
## ar1     0.933222    0.045284  20.60822 0.000000
## ma1    -0.904387    0.054817 -16.49822 0.000000
## omega   0.000011    0.000000  27.40696 0.000000
## alpha1  0.028816    0.013440   2.14406 0.032028
## beta1   0.790288    0.019874  39.76541 0.000000
## gamma1  0.238731    0.049480   4.82480 0.000001
## skew    0.821607    0.030720  26.74491 0.000000
## shape   4.181782    0.413777  10.10637 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000225    0.000391   0.57593 0.564661
## ar1     0.933222    0.050438  18.50226 0.000000
## ma1    -0.904387    0.059184 -15.28083 0.000000
## omega   0.000011    0.000001  20.29122 0.000000
## alpha1  0.028816    0.014698   1.96053 0.049934
## beta1   0.790288    0.018640  42.39653 0.000000
## gamma1  0.238731    0.050069   4.76809 0.000002
## skew    0.821607    0.032279  25.45295 0.000000
## shape   4.181782    0.394834  10.59125 0.000000
## 
## LogLikelihood : 4517.301 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2098
## Bayes        -6.1770
## Shibata      -6.2099
## Hannan-Quinn -6.1976
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.400  0.2368
## Lag[2*(p+q)+(p+q)-1][5]     1.649  0.9939
## Lag[4*(p+q)+(p+q)-1][9]     2.153  0.9771
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.9564  0.3281
## Lag[2*(p+q)+(p+q)-1][5]    1.9853  0.6219
## Lag[4*(p+q)+(p+q)-1][9]    2.8636  0.7813
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.172 0.500 2.000  0.2791
## ARCH Lag[5]     1.406 1.440 1.667  0.6173
## ARCH Lag[7]     1.508 2.315 1.543  0.8198
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  46.3348
## Individual Statistics:               
## mu      0.18787
## ar1     0.18382
## ma1     0.19837
## omega  10.22014
## alpha1  0.16101
## beta1   0.16674
## gamma1  0.11629
## skew    0.09423
## shape   0.14523
## 
## 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.7617 0.4463    
## Negative Sign Bias  1.2376 0.2161    
## Positive Sign Bias  1.4315 0.1525    
## Joint Effect        5.1617 0.1603    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     33.81    0.0193347
## 2    30     53.62    0.0035771
## 3    40     65.96    0.0044625
## 4    50     85.88    0.0008821
## 
## 
## Elapsed time : 1.446039

Phân phối Generalized Error Distribution(ged)

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

VNIged1 <- ugarchfit(VNIged.spec,VNIts) 
print(VNIged1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : ged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001068    0.000189   5.6611 0.000000
## ar1     0.900464    0.014575  61.7794 0.000000
## ma1    -0.875741    0.015767 -55.5413 0.000000
## omega   0.000012    0.000000  35.3305 0.000000
## alpha1  0.025019    0.013332   1.8767 0.060562
## beta1   0.787964    0.019197  41.0455 0.000000
## gamma1  0.203358    0.044491   4.5708 0.000005
## shape   1.055056    0.046279  22.7979 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.001068    0.000132    8.0847 0.000000
## ar1     0.900464    0.007099  126.8491 0.000000
## ma1    -0.875741    0.006682 -131.0578 0.000000
## omega   0.000012    0.000000   28.8490 0.000000
## alpha1  0.025019    0.014840    1.6859 0.091809
## beta1   0.787964    0.020616   38.2218 0.000000
## gamma1  0.203358    0.045173    4.5018 0.000007
## shape   1.055056    0.050514   20.8866 0.000000
## 
## LogLikelihood : 4497.014 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.1832
## Bayes        -6.1541
## Shibata      -6.1833
## Hannan-Quinn -6.1724
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.186  0.2761
## Lag[2*(p+q)+(p+q)-1][5]     1.414  0.9990
## Lag[4*(p+q)+(p+q)-1][9]     1.879  0.9890
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.6882  0.4068
## Lag[2*(p+q)+(p+q)-1][5]    1.3995  0.7646
## Lag[4*(p+q)+(p+q)-1][9]    2.0769  0.8956
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.8937 0.500 2.000  0.3445
## ARCH Lag[5]    1.0191 1.440 1.667  0.7273
## ARCH Lag[7]    1.0909 2.315 1.543  0.8983
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  38.9469
## Individual Statistics:              
## mu      0.2255
## ar1     0.2120
## ma1     0.2307
## omega  10.1467
## alpha1  0.1961
## beta1   0.1442
## gamma1  0.1428
## shape   0.1357
## 
## 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.05976 0.95236    
## Negative Sign Bias 0.58014 0.56191    
## Positive Sign Bias 1.68435 0.09233   *
## Joint Effect       3.76763 0.28767    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     52.79    5.029e-05
## 2    30     73.87    8.806e-06
## 3    40     82.33    6.230e-05
## 4    50    112.46    6.755e-07
## 
## 
## Elapsed time : 1.613254

Phân phối Generalized Error Distribution đối xứng (“sged”)

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

VNIged2 <- ugarchfit(VNIged.spec,VNIts) 
print(VNIged2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : sged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu     -0.000028    0.000289  -0.096336  0.92325
## ar1     0.931718    0.010541  88.389658  0.00000
## ma1    -0.908576    0.012376 -73.414387  0.00000
## omega   0.000011    0.000000  21.739708  0.00000
## alpha1  0.028727    0.002638  10.889689  0.00000
## beta1   0.799823    0.012210  65.505060  0.00000
## gamma1  0.191675    0.036138   5.303921  0.00000
## skew    0.862093    0.007572 113.859909  0.00000
## shape   1.116545    0.038127  29.284826  0.00000
## 
## Robust Standard Errors:
##         Estimate  Std. Error     t value Pr(>|t|)
## mu     -0.000028    0.000339   -0.082027 0.934625
## ar1     0.931718    0.005183  179.772331 0.000000
## ma1    -0.908576    0.006566 -138.368163 0.000000
## omega   0.000011    0.000001   19.418417 0.000000
## alpha1  0.028727    0.009642    2.979279 0.002889
## beta1   0.799823    0.015671   51.037923 0.000000
## gamma1  0.191675    0.031432    6.098091 0.000000
## skew    0.862093    0.020069   42.957164 0.000000
## shape   1.116545    0.048368   23.084523 0.000000
## 
## LogLikelihood : 4512.593 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2033
## Bayes        -6.1706
## Shibata      -6.2034
## Hannan-Quinn -6.1911
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.884  0.1699
## Lag[2*(p+q)+(p+q)-1][5]     2.394  0.8317
## Lag[4*(p+q)+(p+q)-1][9]     2.829  0.9149
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      0.685  0.4079
## Lag[2*(p+q)+(p+q)-1][5]     1.524  0.7339
## Lag[4*(p+q)+(p+q)-1][9]     2.340  0.8608
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     1.015 0.500 2.000  0.3136
## ARCH Lag[5]     1.209 1.440 1.667  0.6721
## ARCH Lag[7]     1.291 2.315 1.543  0.8622
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  42.378
## Individual Statistics:               
## mu      0.13428
## ar1     0.19021
## ma1     0.22836
## omega  10.88679
## alpha1  0.17153
## beta1   0.14555
## gamma1  0.12971
## skew    0.05078
## shape   0.10032
## 
## 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.112 0.2661    
## Negative Sign Bias   1.057 0.2909    
## Positive Sign Bias   1.166 0.2439    
## Joint Effect         5.305 0.1508    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     36.04    0.0104247
## 2    30     59.03    0.0008135
## 3    40     60.78    0.0143217
## 4    50     82.78    0.0018201
## 
## 
## Elapsed time : 2.699843

#Tạo danh sách mô hình lựa chọn
VNI.model.list <- list(
  garch11n = VNIfit,
  garch11t = VNIst1,
  garch11st = VNIst2,
  garch11g = VNIged1,
  garch11sg = VNIged2
)
#Tính toán lựa chọn mô hình
VNI.info.mat <- sapply(VNI.model.list, infocriteria)
rownames(VNI.info.mat) <- rownames(infocriteria(VNIfit))
VNI.info.mat

##               garch11n  garch11t garch11st  garch11g garch11sg
## Akaike       -6.046340 -6.191557 -6.209781 -6.183215 -6.203296
## Bayes        -6.020882 -6.162462 -6.177050 -6.154120 -6.170565
## Shibata      -6.046386 -6.191617 -6.209857 -6.183275 -6.203373
## Hannan-Quinn -6.036841 -6.180700 -6.197568 -6.172359 -6.191083

Trích xuất chuỗi phần dư v của chuỗi lợi suất VNI

VNI.res <- residuals(VNIst2)/sigma(VNIst2)
fitdist(distribution = "sstd", VNI.res, control = list())

## $pars
##           mu        sigma         skew        shape 
## -0.002559858  1.005232150  0.820400209  4.140198895 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2116.226 1930.639 1930.639
## 
## $lagrange
## [1] 0
## 
## $hessian
##            [,1]      [,2]       [,3]      [,4]
## [1,] 2127.21606  589.1352 -562.86284  36.14824
## [2,]  589.13516 1882.8160  309.16160 130.24625
## [3,] -562.86284  309.1616 1378.41303  39.22691
## [4,]   36.14824  130.2462   39.22691  14.00550
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 101
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.5711279 secs
## 
## $vscale
## [1] 1 1 1 1 1

VNI.res

##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 0001-01-01                                  0.73892507
## 0002-01-01                                  1.09800802
## 0003-01-01                                 -0.67931003
## 0004-01-01                                  0.82038760
## 0005-01-01                                  0.87314305
## 0006-01-01                                  0.31932004
## 0007-01-01                                  0.86904768
## 0008-01-01                                  0.05749605
## 0009-01-01                                  1.20994187
## 0010-01-01                                 -3.20784459
##        ...                                            
## 1443-01-01                                 -0.70741848
## 1444-01-01                                 -0.91253404
## 1445-01-01                                  0.40341958
## 1446-01-01                                  0.41410202
## 1447-01-01                                  0.15953903
## 1448-01-01                                  0.05878007
## 1449-01-01                                  1.93800406
## 1450-01-01                                 -0.08795017
## 1451-01-01                                  0.64947130
## 1452-01-01                                  0.02652322

s = pdist("sstd",VNI.res, mu =-0.002559858       
, sigma =  1.005232150, skew = 0.820400209, shape = 4.140198895 )
head(s,10)

##  [1] 0.815652756 0.911745423 0.192954719 0.843198259 0.859122515 0.622573778
##  [7] 0.857938881 0.485465114 0.930198347 0.008260853

Ước lượng mô hình GARCH(1,1) theo phân phối chuẩn cho chỉ số chứng khoán MSCI

MSCIts <- ts(mhnn$msci)

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

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

MSCIfit <- ugarchfit(spec = MSCIspec, MSCIts)
print(MSCIfit)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001251    0.000486   2.5723 0.010103
## ar1    -0.033425    0.028760  -1.1622 0.245156
## omega   0.000023    0.000008   2.9955 0.002740
## alpha1  0.043057    0.021902   1.9659 0.049315
## beta1   0.864722    0.033675  25.6788 0.000000
## gamma1  0.077558    0.028334   2.7373 0.006195
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001251    0.000534  2.34244 0.019158
## ar1    -0.033425    0.029093 -1.14891 0.250594
## omega   0.000023    0.000018  1.29641 0.194835
## alpha1  0.043057    0.047627  0.90404 0.365975
## beta1   0.864722    0.077026 11.22644 0.000000
## gamma1  0.077558    0.041007  1.89134 0.058579
## 
## LogLikelihood : 3618.984 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -4.9766
## Bayes        -4.9547
## Shibata      -4.9766
## Hannan-Quinn -4.9684
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1540  0.6947
## Lag[2*(p+q)+(p+q)-1][2]    0.1583  0.9995
## Lag[4*(p+q)+(p+q)-1][5]    0.4857  0.9942
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01238  0.9114
## Lag[2*(p+q)+(p+q)-1][5]   0.19020  0.9933
## Lag[4*(p+q)+(p+q)-1][9]   0.34310  0.9996
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.03071 0.500 2.000  0.8609
## ARCH Lag[5]   0.04162 1.440 1.667  0.9962
## ARCH Lag[7]   0.17008 2.315 1.543  0.9979
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.3968
## Individual Statistics:              
## mu     0.05221
## ar1    0.26978
## omega  0.39055
## alpha1 0.14334
## beta1  0.32709
## gamma1 0.11767
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.9600 0.3372    
## Negative Sign Bias  0.1877 0.8512    
## Positive Sign Bias  0.7856 0.4322    
## Joint Effect        1.0954 0.7782    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     57.31    1.022e-05
## 2    30     66.14    1.003e-04
## 3    40     85.19    2.734e-05
## 4    50    104.47    6.832e-06
## 
## 
## Elapsed time : 3.093076

Ước lượng mô hình GARCH(1,1) theo phân phối chuẩn cho chỉ số chứng khoán MSCI

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

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

MSCIfit <- ugarchfit(spec = MSCIspec, MSCIts)
print(MSCIfit)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001251    0.000486   2.5723 0.010103
## ar1    -0.033425    0.028760  -1.1622 0.245156
## omega   0.000023    0.000008   2.9955 0.002740
## alpha1  0.043057    0.021902   1.9659 0.049315
## beta1   0.864722    0.033675  25.6788 0.000000
## gamma1  0.077558    0.028334   2.7373 0.006195
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001251    0.000534  2.34244 0.019158
## ar1    -0.033425    0.029093 -1.14891 0.250594
## omega   0.000023    0.000018  1.29641 0.194835
## alpha1  0.043057    0.047627  0.90404 0.365975
## beta1   0.864722    0.077026 11.22644 0.000000
## gamma1  0.077558    0.041007  1.89134 0.058579
## 
## LogLikelihood : 3618.984 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -4.9766
## Bayes        -4.9547
## Shibata      -4.9766
## Hannan-Quinn -4.9684
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1540  0.6947
## Lag[2*(p+q)+(p+q)-1][2]    0.1583  0.9995
## Lag[4*(p+q)+(p+q)-1][5]    0.4857  0.9942
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01238  0.9114
## Lag[2*(p+q)+(p+q)-1][5]   0.19020  0.9933
## Lag[4*(p+q)+(p+q)-1][9]   0.34310  0.9996
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.03071 0.500 2.000  0.8609
## ARCH Lag[5]   0.04162 1.440 1.667  0.9962
## ARCH Lag[7]   0.17008 2.315 1.543  0.9979
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.3968
## Individual Statistics:              
## mu     0.05221
## ar1    0.26978
## omega  0.39055
## alpha1 0.14334
## beta1  0.32709
## gamma1 0.11767
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.9600 0.3372    
## Negative Sign Bias  0.1877 0.8512    
## Positive Sign Bias  0.7856 0.4322    
## Joint Effect        1.0954 0.7782    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     57.31    1.022e-05
## 2    30     66.14    1.003e-04
## 3    40     85.19    2.734e-05
## 4    50    104.47    6.832e-06
## 
## 
## Elapsed time : 1.141091

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

MSCIst1<- ugarchfit(MSCIst.spec,MSCIts) 
print(MSCIst1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001760    0.000408   4.3127 0.000016
## ar1    -0.031348    0.027043  -1.1592 0.246385
## omega   0.000012    0.000002   5.3070 0.000000
## alpha1  0.064207    0.001433  44.7955 0.000000
## beta1   0.871536    0.016944  51.4378 0.000000
## gamma1  0.077802    0.031604   2.4618 0.013825
## shape   5.687384    0.700993   8.1133 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001760    0.000459   3.8340 0.000126
## ar1    -0.031348    0.027519  -1.1392 0.254637
## omega   0.000012    0.000004   2.9663 0.003014
## alpha1  0.064207    0.029733   2.1595 0.030813
## beta1   0.871536    0.020196  43.1532 0.000000
## gamma1  0.077802    0.040628   1.9150 0.055498
## shape   5.687384    1.165937   4.8780 0.000001
## 
## LogLikelihood : 3695.836 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.0810
## Bayes        -5.0556
## Shibata      -5.0811
## Hannan-Quinn -5.0715
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1808  0.6707
## Lag[2*(p+q)+(p+q)-1][2]    0.1831  0.9992
## Lag[4*(p+q)+(p+q)-1][5]    0.5072  0.9933
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.002565  0.9596
## Lag[2*(p+q)+(p+q)-1][5]  0.403854  0.9715
## Lag[4*(p+q)+(p+q)-1][9]  0.736034  0.9947
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1925 0.500 2.000  0.6608
## ARCH Lag[5]    0.2365 1.440 1.667  0.9566
## ARCH Lag[7]    0.5146 2.315 1.543  0.9770
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  13.7595
## Individual Statistics:             
## mu     0.3158
## ar1    0.3625
## omega  2.5302
## alpha1 0.1771
## beta1  0.2211
## gamma1 0.1380
## shape  0.2453
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.9122 0.3618    
## Negative Sign Bias  0.6728 0.5012    
## Positive Sign Bias  0.2111 0.8328    
## Joint Effect        0.9500 0.8133    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     33.73      0.01977
## 2    30     41.10      0.06751
## 3    40     47.72      0.15943
## 4    50     57.37      0.19275
## 
## 
## Elapsed time : 0.860446

Phân phối Student đối xứng (sstd)

MSCIst1.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sstd")

MSCIst2<- ugarchfit(MSCIst1.spec,MSCIts) 
print(MSCIst2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001331    0.000433   3.0741 0.002111
## ar1    -0.028116    0.026961  -1.0428 0.297038
## omega   0.000012    0.000002   6.7180 0.000000
## alpha1  0.054953    0.005237  10.4923 0.000000
## beta1   0.877619    0.015573  56.3560 0.000000
## gamma1  0.084102    0.030204   2.7845 0.005361
## skew    0.904623    0.034083  26.5414 0.000000
## shape   5.786893    0.754279   7.6721 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001331    0.000471   2.8298 0.004658
## ar1    -0.028116    0.027585  -1.0192 0.308092
## omega   0.000012    0.000003   4.0558 0.000050
## alpha1  0.054953    0.022889   2.4008 0.016359
## beta1   0.877619    0.017800  49.3057 0.000000
## gamma1  0.084102    0.038305   2.1956 0.028121
## skew    0.904623    0.035812  25.2600 0.000000
## shape   5.786893    1.133010   5.1075 0.000000
## 
## LogLikelihood : 3699.197 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.0843
## Bayes        -5.0552
## Shibata      -5.0844
## Hannan-Quinn -5.0734
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1319  0.7164
## Lag[2*(p+q)+(p+q)-1][2]    0.1327  0.9998
## Lag[4*(p+q)+(p+q)-1][5]    0.4151  0.9966
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.001059  0.9740
## Lag[2*(p+q)+(p+q)-1][5]  0.367025  0.9761
## Lag[4*(p+q)+(p+q)-1][9]  0.687421  0.9957
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1611 0.500 2.000  0.6881
## ARCH Lag[5]    0.1931 1.440 1.667  0.9671
## ARCH Lag[7]    0.4730 2.315 1.543  0.9808
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  16.0314
## Individual Statistics:             
## mu     0.2740
## ar1    0.3694
## omega  2.9337
## alpha1 0.1785
## beta1  0.2207
## gamma1 0.1414
## skew   1.0348
## shape  0.2297
## 
## 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.9693 0.3326    
## Negative Sign Bias  0.6705 0.5027    
## Positive Sign Bias  0.2827 0.7774    
## Joint Effect        1.0194 0.7965    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     26.40       0.1194
## 2    30     28.41       0.4959
## 3    40     45.74       0.2125
## 4    50     53.65       0.3007
## 
## 
## Elapsed time : 0.9777579

Phân phối Generalized Error Distribution(ged)

MSCIged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "ged")

MSCIged1 <- ugarchfit(MSCIged.spec,MSCIts) 
print(MSCIged1)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : ged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001632    0.000361  4.51671 0.000006
## ar1    -0.020531    0.021592 -0.95086 0.341675
## omega   0.000015    0.000004  3.45919 0.000542
## alpha1  0.048863    0.006442  7.58553 0.000000
## beta1   0.876356    0.014990 58.46146 0.000000
## gamma1  0.084309    0.031384  2.68636 0.007224
## shape   1.269396    0.040578 31.28250 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001632    0.000412   3.9593 0.000075
## ar1    -0.020531    0.019682  -1.0431 0.296880
## omega   0.000015    0.000009   1.5825 0.113527
## alpha1  0.048863    0.041125   1.1881 0.234780
## beta1   0.876356    0.031840  27.5240 0.000000
## gamma1  0.084309    0.045574   1.8499 0.064324
## shape   1.269396    0.127480   9.9576 0.000000
## 
## LogLikelihood : 3679.318 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.0583
## Bayes        -5.0328
## Shibata      -5.0583
## Hannan-Quinn -5.0488
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.003030  0.9561
## Lag[2*(p+q)+(p+q)-1][2]  0.006255  1.0000
## Lag[4*(p+q)+(p+q)-1][5]  0.313970  0.9987
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  9.180e-06  0.9976
## Lag[2*(p+q)+(p+q)-1][5] 3.095e-01  0.9827
## Lag[4*(p+q)+(p+q)-1][9] 5.842e-01  0.9975
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1238 0.500 2.000  0.7249
## ARCH Lag[5]    0.1404 1.440 1.667  0.9788
## ARCH Lag[7]    0.3803 2.315 1.543  0.9879
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  11.6052
## Individual Statistics:             
## mu     0.3949
## ar1    0.3748
## omega  0.4654
## alpha1 0.1536
## beta1  0.2889
## gamma1 0.1305
## shape  0.2556
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.8891 0.3741    
## Negative Sign Bias  0.5057 0.6131    
## Positive Sign Bias  0.4208 0.6739    
## Joint Effect        0.7984 0.8498    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     45.19    0.0006446
## 2    30     50.64    0.0076899
## 3    40     64.69    0.0059925
## 4    50     75.34    0.0091735
## 
## 
## Elapsed time : 0.949688

Phân phối Generalized Error Distribution đối xứng (“sged”)

MSCIged.spec <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(1, 0), include.mean = TRUE), distribution.model = "sged")

MSCIged2 <- ugarchfit(MSCIged.spec,MSCIts) 
print(MSCIged2)

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(1,0,0)
## Distribution : sged 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001205    0.000410  2.94195 0.003262
## ar1    -0.015048    0.024490 -0.61445 0.538916
## omega   0.000014    0.000003  5.42654 0.000000
## alpha1  0.042036    0.010068  4.17515 0.000030
## beta1   0.880951    0.014872 59.23535 0.000000
## gamma1  0.089897    0.030476  2.94979 0.003180
## skew    0.909639    0.029966 30.35607 0.000000
## shape   1.292804    0.049740 25.99103 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.001205    0.000448  2.68859 0.007176
## ar1    -0.015048    0.025326 -0.59417 0.552399
## omega   0.000014    0.000005  3.04505 0.002326
## alpha1  0.042036    0.027761  1.51422 0.129970
## beta1   0.880951    0.025517 34.52409 0.000000
## gamma1  0.089897    0.042032  2.13880 0.032452
## skew    0.909639    0.039896 22.80042 0.000000
## shape   1.292804    0.113724 11.36787 0.000000
## 
## LogLikelihood : 3682.503 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.0613
## Bayes        -5.0322
## Shibata      -5.0614
## Hannan-Quinn -5.0504
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.007459  0.9312
## Lag[2*(p+q)+(p+q)-1][2]  0.008783  1.0000
## Lag[4*(p+q)+(p+q)-1][5]  0.278051  0.9991
## d.o.f=1
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0002052  0.9886
## Lag[2*(p+q)+(p+q)-1][5] 0.2777698  0.9860
## Lag[4*(p+q)+(p+q)-1][9] 0.5416014  0.9980
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.0976 0.500 2.000  0.7547
## ARCH Lag[5]    0.1084 1.440 1.667  0.9853
## ARCH Lag[7]    0.3484 2.315 1.543  0.9900
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  14.3236
## Individual Statistics:             
## mu     0.3230
## ar1    0.4046
## omega  0.6435
## alpha1 0.1578
## beta1  0.2817
## gamma1 0.1315
## skew   0.7389
## shape  0.2361
## 
## 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.9151 0.3603    
## Negative Sign Bias  0.5115 0.6091    
## Positive Sign Bias  0.4574 0.6474    
## Joint Effect        0.8428 0.8392    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     35.66      0.01162
## 2    30     45.36      0.02717
## 3    40     52.96      0.06718
## 4    50     53.23      0.31458
## 
## 
## Elapsed time : 1.609371

#Tạo danh sách mô hình lựa chọn
MSCI.model.list <- list(
  garch11n = MSCIfit,
  garch11t = MSCIst1,
  garch11st = MSCIst2,
  garch11g = MSCIged1,
  garch11sg = MSCIged2
)
#Tính toán lựa chọn mô hình
MSCI.info.mat <- sapply(MSCI.model.list, infocriteria)
rownames(MSCI.info.mat) <- rownames(infocriteria(MSCIfit))
MSCI.info.mat

##               garch11n  garch11t garch11st  garch11g garch11sg
## Akaike       -4.976562 -5.081041 -5.084293 -5.058290 -5.061299
## Bayes        -4.954741 -5.055583 -5.055198 -5.032832 -5.032204
## Shibata      -4.976596 -5.081087 -5.084353 -5.058336 -5.061359
## Hannan-Quinn -4.968420 -5.071541 -5.073437 -5.048790 -5.050442

Trích xuất chuỗi phần dư v của chuỗi lợi suất MSCI

MSCI.res <- residuals(MSCIst2)/sigma(MSCIst2)
fitdist(distribution = "sstd", MSCI.res, control = list())

## $pars
##          mu       sigma        skew       shape 
## -0.01360653  0.99830939  0.89839192  5.84893811 
## 
## $convergence
## [1] 0
## 
## $values
## [1] 2126.516 2005.335 2005.335
## 
## $lagrange
## [1] 0
## 
## $hessian
##             [,1]       [,2]         [,3]       [,4]
## [1,] 1844.872079  210.72229 -446.4081373  2.5316736
## [2,]  210.722291 1928.26847   50.7439268 40.9417852
## [3,] -446.408137   50.74393 1014.0620518  0.2719823
## [4,]    2.531674   40.94179    0.2719823  2.4184551
## 
## $ineqx0
## NULL
## 
## $nfuneval
## [1] 114
## 
## $outer.iter
## [1] 2
## 
## $elapsed
## Time difference of 0.1076548 secs
## 
## $vscale
## [1] 1 1 1 1 1

s1 = pdist("sstd",MSCI.res, mu = -0.01360653, sigma =  0.99830939, skew = 0.89839192, shape = 5.84893811 )
head(s1,10)

##  [1] 0.7611569 0.7856273 0.6965685 0.4914036 0.6185254 0.3895934 0.6498857
##  [8] 0.6856351 0.1859523 0.7417784

#test with Anderson-Darling
library(goftest)
ad_vni <- ad.test(s, "punif")
ad_vni

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

ad_msci <- ad.test(s1, "punif")
ad_msci

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

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

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

cvm_msci <- cvm.test(s1, "punif")
cvm_msci

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

# Kiem dinh ks-test
ks_vni <- ks.test(s, "punif")
ks_vni

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  s
## D = 0.031534, p-value = 0.1114
## alternative hypothesis: two-sided

ks_msci <- ks.test(s1, "punif")
ks_msci

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  s1
## D = 0.014005, p-value = 0.9383
## alternative hypothesis: two-sided

library(VineCopula)

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

BiCopSelect(s,s1, familyset= c(1:10), selectioncrit="AIC",indeptest = FALSE, level = 0.05)

## Bivariate copula: t (par = 0.05, par2 = 10.22, tau = 0.03)

Stu <- BiCopEst(s, s1, family = 2, method = "mle", se = T, max.df = 10)
summary(Stu)

## Family
## ------ 
## No:    2
## Name:  t
## 
## Parameter(s)
## ------------
## par:  0.05  (SE = 0.03)
## par2: 10  (SE = NA)
## Dependence measures
## -------------------
## Kendall's tau:    0.03 (empirical = 0.03, p value = 0.06)
## Upper TD:         0.01 
## Lower TD:         0.01 
## 
## Fit statistics
## --------------
## logLik:  8.67 
## AIC:    -13.33 
## BIC:    -2.77

Tạo copula

library("copula")

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

## 
## Attaching package: 'copula'

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

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

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

##  [1] 0.80247763 0.91190640 0.18719890 0.83757743 0.85960083 0.63799036
##  [7] 0.85684790 0.49346180 0.92842395 0.01238816

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

##  [1] 0.7598073 0.7790778 0.6896077 0.5037853 0.6180317 0.3950447 0.6441844
##  [8] 0.6758431 0.1892636 0.7391604

# Tạo copula
copula_model <- normalCopula(dim = 2)

# Ước lượng copula từ dữ liệu chuẩn hóa
fit_copula <- fitCopula(copula_model, cbind(u,k), method = "ml")

# Xem kết quả ước lượng
summary(fit_copula)

## Call: fitCopula(copula_model, data = cbind(u, k), ... = pairlist(method = "ml"))
## Fit based on "maximum likelihood" and 1452 2-dimensional observations.
## Normal copula, dim. d = 2 
##       Estimate Std. Error
## rho.1   0.0534      0.026
## The maximized loglikelihood is 2.037 
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient 
##        6        6

library(VineCopula)
library(copula)
mycop <- normalCopula(c(-0.1),dim=2,dispstr="ex")
mymvd <- mvdc(copula=mycop,margins =c("norm","norm"),paramMargins=list(list(mean=0,sd=1),list(mean=1,2)))
mymvd

## Multivariate Distribution Copula based ("mvdc")
##  @ copula:
## Normal copula, dim. d = 2 
## Dimension:  2 
## Parameters:
##   rho.1   = -0.1
##  @ margins:
## [1] "norm" "norm"
##    with 2 (not identical)  margins; with parameters (@ paramMargins) 
## List of 2
##  $ :List of 2
##   ..$ mean: num 0
##   ..$ sd  : num 1
##  $ :List of 2
##   ..$ mean: num 1
##   ..$ sd  : num 2

# Tạo dữ liệu mẫu
library(VineCopula)
set.seed(123)
sample_size <- 1452
CPL <- rMvdc(sample_size, mymvd)

# Vẽ biểu đồ copula 2D
library(VineCopula)
plot(mymvd@copula, CPL, n = 1452)

#Vẽ biểu đồ copula 3D
library(rgl)

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

CPL1 <- plot3d(CPL[,1], CPL[,2], dMvdc(CPL, mymvd), 
       col = "darkblue", 
       size = 1,
       xlab = "Dimension 1",
       ylab = "Dimension 2", 
       zlab = "Density")
CPL1

Vẽ scatter-plot với đường phụ thuộc copula

library(copula)
# Tính toán các giá trị phân phối tích lũy
t1 <- pobs(u)
v1 <- pobs(k)
plot(t1, v1, pch = 20, col = "pink")

plot(Stu, add = TRUE , main = "VNI-MSCI: Copula Student-t")

# Vẽ đồ thị contour

library(VineCopula)
plot(Stu, type = "contour", n= 1452, main = "Đồ thị contour của Copula Student-t")

library(VineCopula)

# Tạo đồ thị contour
plot(Stu, type = "contour", n= 500, main = "Đồ thị contour của Copula Student-t")

# Vẽ điểm tâm đỏ ở giữa
points(0, 0, pch = 16, col = "red", cex = 1.5)

Đồ thị mật độ của copula Student

student.cop <- tCopula(param = 0.05, dim = 2)
library(rgl)
persp(student.cop, dCopula, xlab = "t1", ylab = "v1", zlab = "Mật độ",main = "Đồ thị mật độ của copula Student", col="lightblue", expand = 0.45)

# Vẽ biểu đồ 3D
library(scatterplot3d)
scatterplot3d(CPL[,1], CPL[,2], pch = 16, color = "blue",
             type = "h", highlight.3d = TRUE,
             xlab = "U1", ylab = "U2", zlab = "Probability")

## Warning in scatterplot3d(CPL[, 1], CPL[, 2], pch = 16, color = "blue", type =
## "h", : color is ignored when highlight.3d = TRUE

library('rvinecopulib')

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

tcop <- bicop_dist('t', parameters=c(0.05, 10.22))
print(tcop)

## Bivariate copula ('bicop_dist'): family = t, rotation = 0, parameters = 0.05, 10.22, var_types = c,c

plot(tcop)

library(knitr)
library(copula)
library(ggplot2)
library(plotly)

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

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

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

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

library(gridExtra)

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

## 
## Attaching package: 'gridExtra'

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

library(VC2copula)

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

## 
## Attaching package: 'VC2copula'

## The following objects are masked from 'package:VineCopula':
## 
##     BB1Copula, BB6Copula, BB7Copula, BB8Copula, copulaFromFamilyIndex,
##     joeBiCopula, r270BB1Copula, r270BB6Copula, r270BB7Copula,
##     r270BB8Copula, r270ClaytonCopula, r270GumbelCopula,
##     r270JoeBiCopula, r270TawnT1Copula, r270TawnT2Copula, r90BB1Copula,
##     r90BB6Copula, r90BB7Copula, r90BB8Copula, r90ClaytonCopula,
##     r90GumbelCopula, r90JoeBiCopula, r90TawnT1Copula, r90TawnT2Copula,
##     surBB1Copula, surBB6Copula, surBB7Copula, surBB8Copula,
##     surClaytonCopula, surGumbelCopula, surJoeBiCopula, surTawnT1Copula,
##     surTawnT2Copula, tawnT1Copula, tawnT2Copula, vineCopula

library(VineCopula)
library(gridGraphics)

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

## Loading required package: grid

library(png)
#Tạo hàm vẽ biểu đồ
scatter_plot <- function(random_data, cl) {
  ggplot(data.frame(random_data), aes(random_data[,1], random_data[,2])) +
  geom_point(alpha = 0.5, col = cl) +
  theme_minimal() +
  labs(x = "u", y = "k")
  } #for scatter

persp_plot <- function(copula_obj, file_name, cl) {
  png(file_name)
  persp(copula_obj, dCopula, 
        xlab = 'u', ylab = 'k', col = cl, ltheta = 120,  
        ticktype = "detailed", cex.axis = 0.8)
  dev.off()
  rasterGrob(readPNG(file_name),interpolate = TRUE)
}# for pdf

set.seed(123)
#Mô phỏng copula gauss với p=0.8
cop_nor <- normalCopula(param = 0.8, dim = 2)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_nor <- rCopula(copula = cop_nor,n = 1452)
#Mô phỏng copula với p=0.8
cop_std <- tCopula(param = 0.8, dim = 2, df = 1)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_std <- rCopula(copula = cop_std,n = 1452)

#Vẽ biểu đồ 
png('nors.png') #Lưu ảnh biểu đồ
scatter_plot(random_nor, '#D3D3D3')
dev.off()

## png 
##   2

nors <- rasterGrob(readPNG('nors.png'), interpolate = TRUE) #Chuyển ảnh sang dạng Grob

png('stds.png')
scatter_plot(random_std, '#FFB6C1')
dev.off()

## png 
##   2

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

nor_per <- persp_plot(cop_nor, 'norp.png', '#D3D3D3')
std_per <- persp_plot(cop_std, 'stdp.png', '#FFB6C1')

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

grid.arrange(nors, nor_per, stds, std_per, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 1: Biểu đồ phân tán và phối cảnh PDF của Copula họ Elip", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula clayton với p=4
cop_clay <- claytonCopula(param = 4, dim = 2)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_clay <- rCopula(copula = cop_clay,n = 1452)
#Mô phỏng copula gumbel với p=5
cop_gum <- gumbelCopula(param = 5, dim = 2)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_gum <- rCopula(copula = cop_gum,n = 1452)

#Vẽ biểu đồ 

png('clays.png')
scatter_plot(random_clay,'#CC4488')
dev.off()

## png 
##   2

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

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

## png 
##   2

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

clayp <- persp_plot(cop_clay,'clayp.png','#CC4488')

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


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

grid.arrange(clays, clayp, gums, gump, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 2: Biểu đồ phân tán và PDF của Copula Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2)))

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

png('surclays.png')
scatter_plot(random_surclay,'#7766FF')
dev.off()

## png 
##   2

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

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

## png 
##   2

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

surclayp <- persp_plot(cop_surclay, "surclayp.png","#7766FF")
surgump <- persp_plot(cop_surgum, "surgump.png","#FF1199")

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

grid.arrange(surclays, surclayp,surgums, surgump, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 3: Biểu đồ phân tán và PDF của Copula Sur Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula Frank
cop_frank <- frankCopula(param = 9.2)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_frank <- rCopula(copula = cop_frank,n = 1452)
#Mô phỏng copula survival clayton  với p=4
cop_joe <- joeCopula(param = 3)
#Mô phỏng 1452 quan sát ngẫu nhiên dựa trên copula có sẵn
random_joe <- rCopula(copula = cop_joe,n = 1452)

png('franks.png')
scatter_plot(random_frank,'#EEAAEE')
dev.off()

## png 
##   2

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

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

## png 
##   2

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

frankp <- persp_plot(cop_frank, "frankp.png","#EEAAEE")
joep <- persp_plot(cop_joe, "joep.png","#FFCC44")

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

grid.arrange(franks, frankp,joes, joep, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 4: Biểu đồ phân tán và  PDF của Copula Frank và Joe", 
                 gp = gpar(fontsize = 15, font = 2)))