library(googlesheets4)
gs4_deauth()
df <- read_sheet("https://docs.google.com/spreadsheets/d/10nZIVWn3i5djs0tRZbI--8nntNbfbGoNvYKX8yI7Eeo/edit?gid=1109581498#gid=1109581498", sheet = "hose3")
df
# A tibble: 3,335 × 241
Date AAA AAM ABT ACB ACC ACL ADP AGR ANV
<dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2012-04-11 00:00:00 4582 9448 10913 3583 2760 5507 422 7418 1115
2 2012-04-12 00:00:00 4560 9557 10972 3681 2851 5446 461 7418 1115
3 2012-04-13 00:00:00 4476 9557 10972 3681 2851 5446 505 7339 1094
4 2012-04-16 00:00:00 4624 9557 11300 3681 2967 5538 477 7655 1115
5 2012-04-17 00:00:00 4624 9703 11420 3667 2993 5660 521 7971 1126
6 2012-04-18 00:00:00 4560 9849 11330 3611 3110 5568 532 8208 1115
7 2012-04-19 00:00:00 4348 9849 11390 3555 3226 5844 581 7813 1061
8 2012-04-20 00:00:00 4370 9885 11420 3569 3071 6119 581 7497 1105
9 2012-04-23 00:00:00 4370 9885 11420 3611 3006 6241 581 7339 1159
10 2012-04-24 00:00:00 4433 9995 11539 3611 3110 6241 526 7418 1207
# ℹ 3,325 more rows
# ℹ 231 more variables: APG <dbl>, ASM <dbl>, ASP <dbl>, BCE <dbl>, BIC <dbl>,
# BMC <dbl>, BMI <dbl>, BMP <dbl>, BRC <dbl>, BSI <dbl>, BTP <dbl>,
# BTT <dbl>, BVH <dbl>, C32 <dbl>, C47 <dbl>, CCI <dbl>, CCL <dbl>,
# CDC <dbl>, CHP <dbl>, CIG <dbl>, CII <dbl>, CLC <dbl>, CLW <dbl>,
# CMG <dbl>, CMV <dbl>, CMX <dbl>, CNG <dbl>, COM <dbl>, CSM <dbl>,
# CTD <dbl>, CTG <dbl>, CTI <dbl>, CTS <dbl>, CVT <dbl>, D2D <dbl>, …
vn <- df[, c("Date", "VNINDEX")]
vn
# A tibble: 3,335 × 2
Date VNINDEX
<dttm> <dbl>
1 2012-04-11 00:00:00 459.
2 2012-04-12 00:00:00 465.
3 2012-04-13 00:00:00 463.
4 2012-04-16 00:00:00 468.
5 2012-04-17 00:00:00 473.
6 2012-04-18 00:00:00 472.
7 2012-04-19 00:00:00 467.
8 2012-04-20 00:00:00 466.
9 2012-04-23 00:00:00 465.
10 2012-04-24 00:00:00 466.
# ℹ 3,325 more rows
str(vn)
tibble [3,335 × 2] (S3: tbl_df/tbl/data.frame)
$ Date : POSIXct[1:3335], format: "2012-04-11" "2012-04-12" ...
$ VNINDEX: num [1:3335] 459 465 463 468 473 ...
head(vn)
# A tibble: 6 × 2
Date VNINDEX
<dttm> <dbl>
1 2012-04-11 00:00:00 459.
2 2012-04-12 00:00:00 465.
3 2012-04-13 00:00:00 463.
4 2012-04-16 00:00:00 468.
5 2012-04-17 00:00:00 473.
6 2012-04-18 00:00:00 472.
tail(vn)
# A tibble: 6 × 2
Date VNINDEX
<dttm> <dbl>
1 2026-05-15 00:00:00 1922.
2 2026-05-18 00:00:00 1928.
3 2026-05-19 00:00:00 1913.
4 2026-05-20 00:00:00 1913.
5 2026-05-21 00:00:00 1897.
6 2026-05-22 00:00:00 1877.
summary(vn)
Date VNINDEX
Min. :2012-04-11 00:00:00 Min. : 375.3
1st Qu.:2015-09-16 12:00:00 1st Qu.: 594.0
Median :2019-04-23 00:00:00 Median : 957.1
Mean :2019-05-03 01:14:16 Mean : 941.2
3rd Qu.:2022-12-27 12:00:00 3rd Qu.:1228.3
Max. :2026-05-22 00:00:00 Max. :1927.9
4.1. Thống kê mô tả dữ liệu giá chứng khoán
library(fBasics)
basicStats(vn$VNINDEX) # cách 1
X..vn.VNINDEX
nobs 3.335000e+03
NAs 0.000000e+00
Minimum 3.752600e+02
Maximum 1.927940e+03
1. Quartile 5.940250e+02
3. Quartile 1.228350e+03
Mean 9.411805e+02
Median 9.571400e+02
Sum 3.138837e+06
SE Mean 6.318758e+00
LCL Mean 9.287914e+02
UCL Mean 9.535695e+02
Variance 1.331555e+05
Stdev 3.649048e+02
Skewness 4.207220e-01
Kurtosis -5.983190e-01
summary(vn$VNINDEX) # cách 2
Min. 1st Qu. Median Mean 3rd Qu. Max.
375.3 594.0 957.1 941.2 1228.3 1927.9
4.2. Thống kê mô tả dữ liệu return của giá chứng khoán
return_vn <- diff(log(vn$VNINDEX))
head(return_vn)
[1] 0.014112788 -0.005906590 0.012333898 0.009733368 -0.001439154
[6] -0.010817362
tail(return_vn)
[1] -0.0020067279 0.0032939030 -0.0078159773 0.0001568152 -0.0085772103
[6] -0.0104716883
basicStats(return_vn)
return_vn
nobs 3334.000000
NAs 0.000000
Minimum -0.069102
Maximum 0.077067
1. Quartile -0.004576
3. Quartile 0.006742
Mean 0.000423
Median 0.001194
Sum 1.409016
SE Mean 0.000208
LCL Mean 0.000015
UCL Mean 0.000830
Variance 0.000144
Stdev 0.012011
Skewness -0.753109
Kurtosis 4.949390
5.1. Vẽ đồ thị và hàm mật độ giá chứng khoán
ts.plot(vn$VNINDEX) # cách 1: đồ thị giá
plot(vn, type = 'l') # cách 2: đồ thị giá
hist(vn$VNINDEX)
d1 = density(vn$VNINDEX)
plot(d1, type = "l")
5.2. Vẽ đồ thị và hàm mật độ của return giá chứng khoán
plot(return_vn, type = 'l')
hist(return_vn)
d2 = density(return_vn)
plot(d2, type = "l")
library(tseries)
adf.test(vn$VNINDEX)
Augmented Dickey-Fuller Test
data: vn$VNINDEX
Dickey-Fuller = -2.265, Lag order = 14, p-value = 0.4661
alternative hypothesis: stationary
adf.test(return_vn)
Augmented Dickey-Fuller Test
data: return_vn
Dickey-Fuller = -14.706, Lag order = 14, p-value = 0.01
alternative hypothesis: stationary
kpss.test(return_vn) # H0: stationary
KPSS Test for Level Stationarity
data: return_vn
KPSS Level = 0.042477, Truncation lag parameter = 9, p-value = 0.1
normalTest(vn$VNINDEX, method = c("jb")) # H0: normal distribution
Title:
Jarque-Bera Normality Test
Test Results:
STATISTIC:
X-squared: 147.9806
P VALUE:
Asymptotic p Value: < 2.2e-16
normalTest(return_vn, method = c("jb")) # H0: normal distribution
Title:
Jarque-Bera Normality Test
Test Results:
STATISTIC:
X-squared: 3724.973
P VALUE:
Asymptotic p Value: < 2.2e-16
acf(return_vn)
pacf(return_vn)
Box.test(return_vn, lag = 12) # H0: no autocorrelation
Box-Pierce test
data: return_vn
X-squared = 37.152, df = 12, p-value = 0.000211
Box.test(return_vn, lag = 12, type = c("Ljung-Box")) # H0: no autocorrelation
Box-Ljung test
data: return_vn
X-squared = 37.24, df = 12, p-value = 0.0002042
# 8.1. Xây dựng mô hình AR
pacf(return_vn)
m1 <- ar(return_vn, method = c("mle"))
m1
Call:
ar(x = return_vn, method = c("mle"))
Coefficients:
1 2 3 4 5 6 7 8
0.0632 0.0336 -0.0038 -0.0162 0.0286 -0.0156 0.0057 0.0098
9 10 11 12
-0.0322 0.0340 0.0232 -0.0558
Order selected 12 sigma^2 estimated as 0.0001425
m1$aic
0 1 2 3 4 5 6 7
15.325157 4.306591 3.082695 5.078960 6.117447 5.341304 6.663063 8.624042
8 9 10 11 12
10.291214 9.409725 7.676903 8.383212 0.000000
# Ước lượng AR(p)
m21 <- arima(return_vn, order = c(1, 0, 0))
m21
Call:
arima(x = return_vn, order = c(1, 0, 0))
Coefficients:
ar1 intercept
0.0624 4e-04
s.e. 0.0173 2e-04
sigma^2 estimated as 0.0001437: log likelihood = 10019.08, aic = -20032.16
m22 <- arima(return_vn, order = c(12, 0, 0))
m22
Call:
arima(x = return_vn, order = c(12, 0, 0))
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
0.0632 0.0336 -0.0038 -0.0161 0.0286 -0.0157 0.0057 0.0098
s.e. 0.0173 0.0173 0.0173 0.0173 0.0173 0.0173 0.0173 0.0173
ar9 ar10 ar11 ar12 intercept
-0.0321 0.0339 0.0232 -0.0557 4e-04
s.e. 0.0173 0.0173 0.0173 0.0173 2e-04
sigma^2 estimated as 0.0001425: log likelihood = 10032.22, aic = -20036.43
AIC(m21, m22)
df AIC
m21 3 -20032.16
m22 14 -20036.43
# hàm tính p-value
coeftest <- function(model){
coef <- model$coef
se <- sqrt(diag(model$var.coef))
t <- coef / se
p <- 2 * (1 - pnorm(abs(t)))
data.frame(coef, se, t, p)
}
coeftest(m21)
coef se t p
ar1 0.062437853 0.0172868313 3.611874 0.0003039925
intercept 0.000422271 0.0002220688 1.901532 0.0572323850
coeftest(m22)
coef se t p
ar1 0.0631558558 0.0172937242 3.6519523 0.0002602542
ar2 0.0336015173 0.0173247513 1.9395094 0.0524393461
ar3 -0.0037952099 0.0173247080 -0.2190634 0.8266006396
ar4 -0.0161266241 0.0173153689 -0.9313474 0.3516738805
ar5 0.0286019207 0.0173129367 1.6520548 0.0985233738
ar6 -0.0156638994 0.0173230800 -0.9042214 0.3658780252
ar7 0.0056791668 0.0173218005 0.3278624 0.7430157176
ar8 0.0098392809 0.0173123573 0.5683386 0.5698050757
ar9 -0.0321138522 0.0173092381 -1.8553013 0.0635532103
ar10 0.0339422163 0.0173214056 1.9595532 0.0500480323
ar11 0.0231749045 0.0173209618 1.3379687 0.1809066464
ar12 -0.0557080012 0.0172915663 -3.2216862 0.0012743861
intercept 0.0004215796 0.0002240997 1.8812147 0.0599427210
# Kiểm định phần dư sau khi đã chọn bậc cho mô hình
plot(m21$residuals, type = 'l')
qqplot(m21$residuals, return_vn)
Box.test(m21$residuals, lag = 10)
Box-Pierce test
data: m21$residuals
X-squared = 14.949, df = 10, p-value = 0.1339
# 8.2. Xây dựng mô hình MA (chỉ là ví dụ, vì mô hình AR(1) phần dư đã hết tự tương quan)
acf(return_vn)
m3 <- arima(return_vn, order = c(1, 0, 1))
m3
Call:
arima(x = return_vn, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
0.3921 -0.3292 4e-04
s.e. 0.1940 0.1988 2e-04
sigma^2 estimated as 0.0001435: log likelihood = 10020.34, aic = -20032.67
coeftest(m3)
coef se t p
ar1 0.3920670065 0.1940188137 2.020768 0.04330379
ma1 -0.3291660365 0.1987908526 -1.655841 0.09775403
intercept 0.0004165141 0.0002296079 1.814023 0.06967414
# Kiểm định phần dư
plot(m3$residuals, type = 'l')
qqplot(m3$residuals, return_vn)
Box.test(m3$residuals, lag = 10)
Box-Pierce test
data: m3$residuals
X-squared = 12.578, df = 10, p-value = 0.2482
9.1. Kiểm định hiệu ứng ARCH (ARCH-LM test)
library(fGarch)
m4 <- garchFit(~ arma(1,0) + garch(1,0), data = return_vn, trace = FALSE)
summary(m4)
Title:
GARCH Modelling
Call:
garchFit(formula = ~arma(1, 0) + garch(1, 0), data = return_vn,
trace = FALSE)
Mean and Variance Equation:
data ~ arma(1, 0) + garch(1, 0)
<environment: 0x306baa4e0>
[data = return_vn]
Conditional Distribution:
norm
Coefficient(s):
mu ar1 omega alpha1
0.00054347 0.03762936 0.00010369 0.31379950
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
mu 5.435e-04 1.896e-04 2.866 0.00416 **
ar1 3.763e-02 2.271e-02 1.657 0.09753 .
omega 1.037e-04 3.413e-06 30.376 < 2e-16 ***
alpha1 3.138e-01 3.315e-02 9.465 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log Likelihood:
10148.65 normalized: 3.043986
Description:
Tue Jun 9 12:22:16 2026 by user:
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 3333.392168 0.0000000000
Shapiro-Wilk Test R W 0.937498 0.0000000000
Ljung-Box Test R Q(10) 24.669549 0.0060083299
Ljung-Box Test R Q(15) 39.082444 0.0006229258
Ljung-Box Test R Q(20) 41.750664 0.0029800275
Ljung-Box Test R^2 Q(10) 207.983937 0.0000000000
Ljung-Box Test R^2 Q(15) 267.981576 0.0000000000
Ljung-Box Test R^2 Q(20) 277.183072 0.0000000000
LM Arch Test R TR^2 165.246026 0.0000000000
Information Criterion Statistics:
AIC BIC SIC HQIC
-6.085573 -6.078240 -6.085576 -6.082949
9.2. Ước lượng ARIMA+GARCH
m5 <-garchFit(~ arma(1,0) + garch(1,1), data = return_vn, trace = FALSE)
summary(m5)
Title:
GARCH Modelling
Call:
garchFit(formula = ~arma(1, 0) + garch(1, 1), data = return_vn,
trace = FALSE)
Mean and Variance Equation:
data ~ arma(1, 0) + garch(1, 1)
<environment: 0x30588d4a0>
[data = return_vn]
Conditional Distribution:
norm
Coefficient(s):
mu ar1 omega alpha1 beta1
5.9490e-04 7.4157e-02 6.1102e-06 1.2476e-01 8.3884e-01
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
mu 5.949e-04 1.712e-04 3.475 0.000510 ***
ar1 7.416e-02 1.968e-02 3.768 0.000164 ***
omega 6.110e-06 1.050e-06 5.821 5.84e-09 ***
alpha1 1.248e-01 1.325e-02 9.418 < 2e-16 ***
beta1 8.388e-01 1.602e-02 52.367 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log Likelihood:
10340.51 normalized: 3.101533
Description:
Tue Jun 9 12:22:16 2026 by user:
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 4492.1725549 0.00000000
Shapiro-Wilk Test R W 0.9456598 0.00000000
Ljung-Box Test R Q(10) 15.2766879 0.12229521
Ljung-Box Test R Q(15) 26.4199217 0.03383384
Ljung-Box Test R Q(20) 28.0285591 0.10872504
Ljung-Box Test R^2 Q(10) 9.2064302 0.51263121
Ljung-Box Test R^2 Q(15) 10.1713349 0.80882626
Ljung-Box Test R^2 Q(20) 13.1578268 0.87051012
LM Arch Test R TR^2 10.0226301 0.61397524
Information Criterion Statistics:
AIC BIC SIC HQIC
-6.200066 -6.190900 -6.200070 -6.196786
predict(m5, 5)
meanForecast meanError standardDeviation
1 -0.0001816455 0.009490514 0.009490514
2 0.0005814312 0.009664193 0.009638532
3 0.0006380185 0.009805269 0.009779043
4 0.0006422149 0.009939188 0.009912555
5 0.0006425260 0.010066547 0.010039527
# Chạy hàm loss để tính VaR
nreturn_vn <- -return_vn
# 10.1. Ước tính VaR theo RISK METRICS
# Cách 1: tính thủ công
library(rugarch)
spec1 <- ugarchspec(mean.model = list(armaOrder=c(1,0)),
variance.model = list(model="iGARCH", garchOrder=c(1,1)))
m6 <- ugarchfit(spec=spec1,data = nreturn_vn)
Box.test(residuals(m6, standardize=TRUE), lag=10, type="Ljung")
Box-Ljung test
data: residuals(m6, standardize = TRUE)
X-squared = 14.848, df = 10, p-value = 0.1377
ugarchforecast(m6, n.ahead = 1)
*------------------------------------*
* GARCH Model Forecast *
*------------------------------------*
Model: iGARCH
Horizon: 1
Roll Steps: 0
Out of Sample: 0
0-roll forecast [T0=1979-02-17]:
Series Sigma
T+1 0.0001809 0.009661
# Tính VaR (Phương pháp của JP Morgan xem mean = 0):
VaR1 <- 0 + qnorm(0.95) * 0.009661
# Cách 2: tính từ RMeasure.R
setwd("~/Library/CloudStorage/OneDrive-Personal/0. Bai giang/QTRRTC/R")
source("RMeasure.R")
VaR2 <- RMeasure(0, 0.009661)
Risk Measures for selected probabilities:
prob VaR ES
[1,] 0.950 0.01589093 0.01992787
[2,] 0.990 0.02247485 0.02574863
[3,] 0.999 0.02985473 0.03252946
# 10.2. Ước tính VaR theo ECONOMICTRICS
m7 <- garchFit(~ arma(1,0) + garch(1,1), data = nreturn_vn, trace=F)
summary(m7)
Title:
GARCH Modelling
Call:
garchFit(formula = ~arma(1, 0) + garch(1, 1), data = nreturn_vn,
trace = F)
Mean and Variance Equation:
data ~ arma(1, 0) + garch(1, 1)
<environment: 0x30e5b9580>
[data = nreturn_vn]
Conditional Distribution:
norm
Coefficient(s):
mu ar1 omega alpha1 beta1
-5.9490e-04 7.4157e-02 6.1102e-06 1.2476e-01 8.3884e-01
Std. Errors:
based on Hessian
Error Analysis:
Estimate Std. Error t value Pr(>|t|)
mu -5.949e-04 1.712e-04 -3.475 0.000510 ***
ar1 7.416e-02 1.968e-02 3.768 0.000164 ***
omega 6.110e-06 1.050e-06 5.821 5.84e-09 ***
alpha1 1.248e-01 1.325e-02 9.418 < 2e-16 ***
beta1 8.388e-01 1.602e-02 52.367 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Log Likelihood:
10340.51 normalized: 3.101533
Description:
Tue Jun 9 12:22:16 2026 by user:
Standardised Residuals Tests:
Statistic p-Value
Jarque-Bera Test R Chi^2 4492.1724364 0.00000000
Shapiro-Wilk Test R W 0.9456598 0.00000000
Ljung-Box Test R Q(10) 15.2766878 0.12229521
Ljung-Box Test R Q(15) 26.4199217 0.03383384
Ljung-Box Test R Q(20) 28.0285592 0.10872504
Ljung-Box Test R^2 Q(10) 9.2064300 0.51263122
Ljung-Box Test R^2 Q(15) 10.1713347 0.80882627
Ljung-Box Test R^2 Q(20) 13.1578265 0.87051013
LM Arch Test R TR^2 10.0226299 0.61397526
Information Criterion Statistics:
AIC BIC SIC HQIC
-6.200066 -6.190900 -6.200070 -6.196786
predict(m7, 1)
meanForecast meanError standardDeviation
1 0.0001816456 0.009490514 0.009490514
# Tính VaR:
VaR3 <- 0.0001816456 + qnorm(0.95)*0.009490514
# 10.3. Ước lượng QUANTILE ESTIMATION
VaR4 <- quantile(nreturn_vn, 0.95)
# 10.4. Ước lượng MONTE CARLO SIMULATION
mean_vn <- mean(nreturn_vn)
sd_vn <- sd(nreturn_vn)
# Tiến hành mô phỏng Monte Carlo
set.seed(42)
sim1 <- rnorm(10000, mean = mean_vn, sd = sd_vn)
ts.plot(sim1)
VaR5 <- quantile(sim1, 0.95)
VaR1_100 = VaR1*100000
VaR1_100
[1] 1589.093
VaR3_100 = VaR3*100000
VaR3_100
[1] 1579.215
VaR4_100 = VaR4*100000
VaR4_100
95%
1893.13
VaR5_100 = VaR5*100000
VaR5_100
95%
1935.795