library(tibble)
library(httr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
library(xml2)
library(rvest)
library(quantmod)
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(tidyquant)
## Loading required package: lubridate
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
## Loading required package: PerformanceAnalytics
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(ggplot2)
library(tseries)
library(moments)
##
## Attaching package: 'moments'
## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewness
library(AdequacyModel)
library(rugarch)
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
library(lmtest)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(forecast)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(GAS)
##
## Attaching package: 'GAS'
## The following objects are masked from 'package:rugarch':
##
## convergence, pit, residuals
## The following object is masked from 'package:PerformanceAnalytics':
##
## ES
## The following object is masked from 'package:stats':
##
## residuals
library(knitr)
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 object is masked from 'package:TTR':
##
## volatility
library(purrr)
##
## Attaching package: 'purrr'
## The following object is masked from 'package:car':
##
## some
## The following object is masked from 'package:rugarch':
##
## reduce
library(quarks)
## ********************************************************************************
## Welcome to the package 'quarks'!
## ********************************************************************************
##
## Please report any possible errors and bugs to sebastian.letmathe@uni-paderborn.de.
## Thank you.
##
## ********************************************************************************
#set data
VN30_Index<- read.csv("D:/Jupyter/Untitled Folder/new_d.csv")
VN30_Index <- VN30_Index %>%
mutate(symbol = "^VN30")
VN30_Index$date <- as.Date(VN30_Index$date)
VN30_Index %>%
ggplot(aes(date, Closed_price, color = symbol))+
geom_line(show.legend = FALSE)+
labs(x = "Year" , y = "INDEX", title = "Figure 1: VN30 Index Trend")
VN30_Index %>%
group_by(symbol)%>%
mutate(Log_Price_Close = log(Closed_price)) %>%
ungroup() %>%
select(symbol,date,Log_Price_Close) -> VN30_Index_log
VN30_Index_Modified<-VN30_Index[,c(1,5)]
Date = as.Date(VN30_Index_Modified$date, format = "%Y%m%d")
VN30_Index_Modified<-VN30_Index_Modified[order(Date),]
VN30_Index_Modified1<-xts(x = VN30_Index_Modified$Closed_price, order.by = Date)
names(VN30_Index_Modified1)<- "VN30"
VN30.logret<- diff(log(VN30_Index_Modified1), lag = 1)
VN30.logret[1,]<- 0
hist(VN30.logret[1:1960] ,breaks=50,col="deeppink1",type="l")
## Warning in plot.window(xlim, ylim, "", ...): graphical parameter "type" is
## obsolete
## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## graphical parameter "type" is obsolete
## Warning in axis(1, ...): graphical parameter "type" is obsolete
## Warning in axis(2, ...): graphical parameter "type" is obsolete
lines(x = density(x = VN30.logret[1:1960]), col = "blue")
plot(VN30.logret[1:1960], col = "deeppink1", lwd = 1)
ind.summary<- c("Average" = mean(VN30.logret[1:1960]) , "Standard Deviation" = sd(VN30.logret), "Skewness" = skewness(VN30.logret), "Kurtosis" = kurtosis(VN30.logret))
ind.summary %>%
kbl(caption = "Table 2: Data indicators summary", escape = TRUE)%>%
kable_classic(full_width = FALSE, html_font = "cambria")
Table 2: Data indicators summary
|
|
x
|
|
Average
|
0.0004452
|
|
Standard Deviation
|
0.0709336
|
|
Skewness.VN30
|
-0.1993389
|
|
Kurtosis.VN30
|
10.4149608
|
acf(VN30.logret[1:1960])
pacf(VN30.logret[1:1960])
adf.test(VN30.logret[1:1960])
## Warning in adf.test(VN30.logret[1:1960]): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: VN30.logret[1:1960]
## Dickey-Fuller = -18.44, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(VN30.logret[1:1960], null="Trend",lshort=FALSE)
## Warning in kpss.test(VN30.logret[1:1960], null = "Trend", lshort = FALSE):
## p-value greater than printed p-value
##
## KPSS Test for Trend Stationarity
##
## data: VN30.logret[1:1960]
## KPSS Trend = 0.038056, Truncation lag parameter = 25, p-value = 0.1
model.arima = auto.arima(VN30.logret[1:1960] , max.order = c(10,10,10) , stationary = TRUE , trace = T , ic = 'aicc')
##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2) with non-zero mean : -5863.833
## ARIMA(0,0,0) with non-zero mean : -5236.384
## ARIMA(1,0,0) with non-zero mean : -5659.806
## ARIMA(0,0,1) with non-zero mean : -5825.391
## ARIMA(0,0,0) with zero mean : -5238.292
## ARIMA(1,0,2) with non-zero mean : -5877.706
## ARIMA(0,0,2) with non-zero mean : -5828.298
## ARIMA(1,0,1) with non-zero mean : -5832.243
## ARIMA(1,0,3) with non-zero mean : -5881.408
## ARIMA(0,0,3) with non-zero mean : -5849.028
## ARIMA(2,0,3) with non-zero mean : -5866.024
## ARIMA(1,0,4) with non-zero mean : -5880.984
## ARIMA(0,0,4) with non-zero mean : -5857.6
## ARIMA(2,0,4) with non-zero mean : -5869.195
## ARIMA(1,0,3) with zero mean : -5879.599
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(1,0,3) with non-zero mean : -5882.142
##
## Best model: ARIMA(1,0,3) with non-zero mean
model.arima
## Series: VN30.logret[1:1960]
## ARIMA(1,0,3) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 ma3 mean
## 0.6968 -1.3230 0.4386 -0.0592 4e-04
## s.e. 0.0542 0.0577 0.0483 0.0243 2e-04
##
## sigma^2 = 0.0029: log likelihood = 2947.09
## AIC=-5882.18 AICc=-5882.14 BIC=-5848.7
model.arima$residuals %>% ggtsdisplay(plot.type = 'hist' , lag.max = 14)
ar.res = model.arima$residuals
Box.test(model.arima$residuals , lag = 7 , fitdf = 2 , type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: model.arima$residuals
## X-squared = 10.46, df = 5, p-value = 0.06321
adf.test(model.arima$residuals)
## Warning in adf.test(model.arima$residuals): p-value smaller than printed
## p-value
##
## Augmented Dickey-Fuller Test
##
## data: model.arima$residuals
## Dickey-Fuller = -13.278, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(model.arima$residuals, null="Trend",lshort=FALSE)
## Warning in kpss.test(model.arima$residuals, null = "Trend", lshort = FALSE):
## p-value greater than printed p-value
##
## KPSS Test for Trend Stationarity
##
## data: model.arima$residuals
## KPSS Trend = 0.11239, Truncation lag parameter = 25, p-value = 0.1
jarque.bera.test(model.arima$residuals)
##
## Jarque Bera Test
##
## data: model.arima$residuals
## X-squared = 3626.8, df = 2, p-value < 2.2e-16
pacf(model.arima$residuals)
tsdisplay(ar.res^2 , main = 'Squared Residuals')
model.spec = ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1, 1)),mean.model = list(armaOrder = c(0, 0)))
model.fit = ugarchfit(spec = model.spec , data = ar.res , solver = 'solnp')
model.fit@fit$matcoef
## Estimate Std. Error t value Pr(>|t|)
## mu -9.271511e-04 9.131351e-04 -1.015349 3.099393e-01
## omega 2.480356e-05 5.424337e-06 4.572644 4.816073e-06
## alpha1 5.276220e-02 6.637215e-03 7.949449 1.776357e-15
## beta1 9.392346e-01 6.602173e-03 142.261438 0.000000e+00
quantile(VN30.logret[1:1960] , 0.05)
## 5%
## -0.1000222
qplot(VN30.logret[1:1960] , geom = 'histogram') + geom_histogram(fill = 'lightblue' , bins = 30) +geom_histogram(aes(VN30.logret[1:1960][VN30.logret[1:1960] < quantile(VN30.logret[1:1960] , 0.05)]) , fill = 'red' , bins = 30) +labs(x = 'Daily Log Returns')
## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Don't know how to automatically pick scale for object of type <xts/zoo>.
## Defaulting to continuous.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
qqnorm(VN30.logret[1:1960])
qqline(VN30.logret[1:1960])
jarque.bera.test(ar.res^2)
##
## Jarque Bera Test
##
## data: ar.res^2
## X-squared = 473996, df = 2, p-value < 2.2e-16
jarque.bera.test(ar.res)
##
## Jarque Bera Test
##
## data: ar.res
## X-squared = 3626.8, df = 2, p-value < 2.2e-16
jarque.bera.test(VN30.logret[1:1960])
##
## Jarque Bera Test
##
## data: VN30.logret[1:1960]
## X-squared = 4174.7, df = 2, p-value < 2.2e-16
#create density plots
curve(dt(x, df=7), from=-4, to=4, col='blue')
curve(dt(x, df=30), from=-4, to=4, col='green', add=TRUE)
curve(dnorm(x, mean =0, sd=1), from=-4, to=4, col='red', add=TRUE)
#add legend
legend(-4, .3, legend=c("df=2", "df=30","normal"),
col=c("blue", "green","red"), lty=1, cex=1.2)
p2_1 = qplot(VN30.logret[1:1960] , geom = 'density') + geom_density(fill = 'blue' , alpha = 0.4) +
geom_density(aes(rnorm(200000 , 0 , sd(VN30.logret[1:1960]))) , fill = 'red' , alpha = 0.25) +
labs(x = '')
p2_2 = qplot(VN30.logret[1:1960] , geom = 'density') + geom_density(fill = 'blue' , alpha = 0.4) +
geom_density(aes(rnorm(200000 , 0 , sd(VN30.logret[1:1960]))) , fill = 'red' , alpha = 0.25) +
coord_cartesian(xlim = c(-0.07 , -0.02) , ylim = c(0 , 10)) +
geom_vline(xintercept = c(qnorm(p = c(0.01 ,0.025, 0.05, 0.1) , mean = mean(VN30.logret[1:1960]) , sd = sd(VN30.logret[1:1960]))) ,
color = c('darkgreen' , 'green', 'blue','red') , size = 1) + labs(x = 'Daily Log Returns')
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
grid.arrange(p2_1 , p2_2 , ncol = 1)
## Don't know how to automatically pick scale for object of type <xts/zoo>.
## Defaulting to continuous.
## Don't know how to automatically pick scale for object of type <xts/zoo>.
## Defaulting to continuous.
fitdist(distribution = 'std' , x = VN30.logret[1:1960])$pars
## mu sigma shape
## 0.00110026 0.31295155 2.01000001
shape <- fitdist(distribution = 'std' , x = VN30.logret[1:1960])$pars['shape']
cat("For a = 0.1 the quantile value of normal distribution is: ",
qnorm(p=0.1),"\n" ,
"For a = 0.1 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = shape , p = 0.1) , "\n" , "\n" ,
"For a = 0.025 the quantile value of normal distribution is: ",
qnorm(p=0.025),"\n" ,
"For a = 0.025 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = shape , p = 0.025) , "\n" , "\n" ,
"For a = 0.05 the quantile value of normal distribution is: " ,
qnorm(p = 0.05) , "\n" ,
"For a = 0.05 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = shape , p = 0.05) , "\n" , "\n" ,
'For a = 0.01 the quantile value of normal distribution is: ' ,
qnorm(p = 0.01) , "\n" ,
"For a = 0.01 the quantile value of t-distribution is: " ,
qdist(distribution = 'std' , shape = shape , p = 0.01) , sep = "")
## For a = 0.1 the quantile value of normal distribution is: -1.281552
## For a = 0.1 the quantile value of t-distribution is: -0.1327103
##
## For a = 0.025 the quantile value of normal distribution is: -1.959964
## For a = 0.025 the quantile value of t-distribution is: -0.3020435
##
## For a = 0.05 the quantile value of normal distribution is: -1.644854
## For a = 0.05 the quantile value of t-distribution is: -0.2052629
##
## For a = 0.01 the quantile value of normal distribution is: -2.326348
## For a = 0.01 the quantile value of t-distribution is: -0.4879308
# CL 99%
qbinom(0.95,494,0.01,lower.tail=TRUE,log.p=FALSE)
## [1] 9
#32
qbinom(0.9999,494,0.01,lower.tail=TRUE,log.p=FALSE)
## [1] 15
# CL 97.5%
qbinom(0.95,494,0.025,lower.tail=TRUE,log.p=FALSE)
## [1] 18
qbinom(0.9999,494,0.025,lower.tail=TRUE,log.p=FALSE)
## [1] 27
# CL 95%
qbinom(0.95,494,0.05,lower.tail=TRUE,log.p=FALSE)
## [1] 33
qbinom(0.9999,494,0.05,lower.tail=TRUE,log.p=FALSE)
## [1] 44
# CL 90%
qbinom(0.95,494,0.1,lower.tail=TRUE,log.p=FALSE)
## [1] 61
qbinom(0.9999,494,0.1,lower.tail=TRUE,log.p=FALSE)
## [1] 76
binom_sum <- function(size, prob, lb, ub, col = 4, lwd = 1, ...) {
x <- 0:size
if (missing(lb)) {
lb <- min(x)
}
if (missing(ub)) {
ub <- max(x)
}
plot(dbinom(x, size = size, prob = prob), type = "h", lwd = lwd, ...)
if(lb == min(x) & ub == max(x)) {
color <- col
} else {
color <- rep(1, length(x))
color[(lb + 1):ub ] <- col
}
lines(dbinom(x, size = size, prob = prob), type = "h",
col = color, lwd = lwd, ...)
}
binom_sum(size = 494, prob = 0.01, ub =5, lwd = 1,
ylab = "P(X = x)", xlab = "Number of successes")
binom_sum(size = 494, prob = 0.01, ub =10, lwd = 1,
ylab = "P(X = x)", xlab = "Number of successes")
model.roll = ugarchroll(spec = model.spec , data = VN30.logret , n.start = 1960 , ,refit.window = 'moving')
VaR90_td = model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=shape, p=0.1)
qplot(y = VaR90_td , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR90_td)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with GARCH approach: ' , (sum(VN30.logret[1961:2454] < VaR90_td)) , sep = '')
## Number of exceptions with GARCH approach: 135
VaRTest(alpha=0.1,actual=VN30.logret[1961:2454],VaR=VaR90_td,conf.level = 0.95)
## $expected.exceed
## [1] 49
##
## $actual.exceed
## [1] 135
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 117.8914
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 118.1979
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist <- function(alpha, actual, var)
{
n00 <- 0
n01 <- 0
n10 <- 0
n11 <- 0
if (actual[1]<var[1])
{
n01 <- n01+1
}
if (actual[1]>var[1])
{
n00 <- n00+1
}
for (i in 2:length(var))
{
if (actual[i]<var[i] & actual[i-1] < var[i-1])
{
n11 <- n11+1
}
if (actual[i]<var[i] & actual[i-1] > var[i-1])
{
n01 <- n01+1
}
if (actual[i]>var[i] & actual[i-1] > var[i-1])
{
n00 <- n00+1
}
if (actual[i]>var[i] & actual[i-1] < var[i-1])
{
n10 <- n10+1
}
}
pi0 <- n01/(n00+n01)
pi1 <- n11/(n10+n11)
pi <- (n01+n11)/(n10+n11+n00+n01)
tuso <- (1-pi)^(n00+n10)*pi^(n01+n11)
mauso <- (1-pi0)^n00*pi0^n01*(1-pi1)^n10*pi1^n11
z <- -2*log(tuso/mauso)
if (z < qchisq(p=1-alpha, df=1, ncp = 0, lower.tail = TRUE, log.p = FALSE))
{
cat("Accept H0 with LRstat is: ",z)
}
if (z > qchisq(p=1-alpha, df=1, ncp = 0, lower.tail = TRUE, log.p = FALSE))
{
cat("Reject H0 with LRstat is: ",z)
}
}
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR90_td)
## Accept H0 with LRstat is: 0.3572901
VaRDurTest_Haas <- function(alpha, actual, var,conf)
{
vec <- c(0)
stat <- c()
for (i in 1:length(var))
{
if (actual[i]<var[i])
{
vec <- cbind(vec,i)
}
}
for (j in 1:length(vec))
{
stat[j] <- -2*log((alpha*(1-alpha)^(vec[j+1]-vec[j]-1))/
((1/(vec[j+1]-vec[j]))*(1-(1/(vec[j+1]-vec[j])))^(vec[j+1]-vec[j]-1)))
}
z <- sum(na.omit(stat))
if (z < qchisq(p=conf, df=length(vec)-1, ncp = 0, lower.tail = TRUE, log.p = FALSE))
{
cat("Accept H0 with LRstat is: ",z)
}
if (z > qchisq(p=conf, df=length(vec)-1, ncp = 0, lower.tail = TRUE, log.p = FALSE))
{
cat("Reject H0 with LRstat is: ",z)
}
}
VaRDurTest_Haas(alpha=0.1,actual=as.numeric(VN30.logret[1961:2454]),var=VaR90_td,conf=0.95)
## Reject H0 with LRstat is: 262.6247
my_function <- function(x)
{
qt(p=x,df=shape)
}
ES90_td = (1/(1-0.1))*model.roll@forecast$density[,'Sigma']*integrate(my_function,lower = 0.0000001,upper = 0.1)$value*(1/(shape-1))*(shape+(qt(p=0.1,df=shape))^2)
qplot(y = VaR90_td , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES90_td,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR90_td)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
# Acerbi/Szekely - Test 1 (2014)
as1 = function(x,var,es){
# Inputs:
# x | time series
# var | VaR prediction
# es | ES prediction
# Output:
# z | test statistic
z = 0
n = length(x)
nt= 0
for (i in 1:n)
{
if (x[i] < as.numeric(var[i]))
{
nt = nt + 1
z = z + x[i]/as.numeric(-es[i])
}
}
z = z/nt + 1
return(z)
}
z1_90 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR90_td,es=ES90_td)
z1_90
## [1] 0.545887
# Acerbi/Szekely - Test 2 (2014)
as2 = function(x,var,es,alpha){
# Inputs:
# x | time series
# var | VaR prediction
# es | ES prediction
# alpha | confidence level
# Output:
# z | test statistic
n = length(x)
z = 0
for (i in 1:n){
if (x[i] < var[i])
{
z = z + (x[i])/(n*alpha*-es[i])
}
}
z = z+1
return(z)
}
z2_90 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR90_td,es=ES90_td,alpha=0.05)
z2_90
## [1] -1.481994
forc <- c(as.numeric(VN30.logret[1:1960]),model.roll@forecast$density[,'Sigma'])*rdist(distribution='std', shape=shape, n=490)
## Warning in c(as.numeric(VN30.logret[1:1960]), model.roll@forecast$density[, :
## longer object length is not a multiple of shorter object length
DoGARCH_VaR <- function(y, probability, portfolio_value = 1, WE)
{
# Number of observations
n <- length(y)
# Initialize empty VaR vector
VaR <- rep(0, n)
# Do a loop for the forecast
for (i in 1:(n-WE)){
# Subset the dataset to the estimation window
window <- y[i:(i+WE-1)]
op <- ceiling(length(window)*probability)
# Allocate the VaR forecast in the vector
VaR[i+WE] <- sort(as.numeric(window))[op]
}
# Return the VaR vector
return(VaR)
}
VaR_90_his <- DoGARCH_VaR(y=forc, probability=0.1, WE=1960)
VaR_90_his <- VaR_90_his[1961:2454]
qplot(y = VaR_90_his , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_90_his)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with His approach: ' , (sum(VN30.logret[1961:2454] < VaR_90_his)) , sep = '')
## Number of exceptions with His approach: 168
VaRTest(alpha=0.1,actual=VN30.logret[1:494],VaR=VaR_90_his,conf.level = 0.95)
## $expected.exceed
## [1] 49
##
## $actual.exceed
## [1] 168
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 208.9675
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 221.3252
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRTest(alpha=0.1,actual=VN30.logret[1961:2454],VaR=VaR_90_his,conf.level = 0.95)
## $expected.exceed
## [1] 49
##
## $actual.exceed
## [1] 168
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 208.9675
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 210.2328
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_90_his)
## Accept H0 with LRstat is: 1.365219
VaRDurTest_Haas(alpha=0.1,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_90_his,conf=0.95)
## Reject H0 with LRstat is: 371.071
DoGARCH_ES <- function(y, probability, portfolio_value = 1, WE)
{
# Number of observations
n <- length(y)
# Initialize empty VaR vector
ES <- rep(0, n)
# Do a loop for the forecast
for (i in 1:(n-WE)){
# Subset the dataset to the estimation window
window <- y[i:(i+WE-1)]
op <- ceiling(length(window)*probability)
# Allocate the VaR forecast in the vector
ES[i+WE] <- mean(sort(window)[1:op])
}
# Return the VaR vector
return(ES)
}
ES_90_his <- DoGARCH_ES(y=forc, probability=0.1, WE=1960)
ES_90_his <- ES_90_his[1961:2454]
qplot(y = VaR_90_his , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES_90_his,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_90_his)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_his_90 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR_90_his,es=ES_90_his)
z1_his_90
## [1] -2.423387
z2_his_90 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR_90_his,es=ES_90_his,alpha=0.05)
z2_his_90
## [1] -22.28457
model.roll = ugarchroll(spec = model.spec , data = VN30.logret , n.start = 1960 , ,refit.window = 'moving')
VaR95_td = model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=shape, p=0.05)
qplot(y = VaR95_td , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR95_td)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with GARCH approach: ' , (sum(VN30.logret[1961:2454] < VaR95_td)) , sep = '')
## Number of exceptions with GARCH approach: 120
VaRTest(alpha=0.05,actual=VN30.logret[1961:2454],VaR=VaR95_td,conf.level = 0.95)
## $expected.exceed
## [1] 24
##
## $actual.exceed
## [1] 120
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 209.5794
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 209.6112
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR95_td)
## Accept H0 with LRstat is: 0.04975041
VaRDurTest_Haas(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR95_td,conf=0.95)
## Reject H0 with LRstat is: 344.1655
ES95_td = (1/(1-0.05))*model.roll@forecast$density[,'Sigma']*integrate(my_function,lower = 0.0000001,upper = 0.05)$value*(1/(shape-1))*(shape+(qt(p=0.05,df=shape))^2)
qplot(y = VaR95_td , x = 1:494, geom = 'line') +
geom_line(aes(x=1:494,y=ES95_td,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR95_td)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
VN30.logret1 <- diff(log(VN30_Index_Modified$Closed_price),lag=1)
z1_95 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR95_td,es=ES95_td)
z1_95
## [1] 0.6123557
z2_95 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR95_td,es=ES95_td,alpha=0.05)
z2_95
## [1] -0.883292
VaR_95_his <- DoGARCH_VaR(y=forc, probability=0.05, WE=1960)
VaR_95_his <- VaR_95_his[1961:2454]
qplot(y = VaR_95_his , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_95_his)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with His approach: ' , (sum(VN30.logret[1961:2454] < VaR_95_his)) , sep = '')
## Number of exceptions with His approach: 127
VaRTest(alpha=0.05,actual=VN30.logret[1961:2454],VaR=VaR_95_his,conf.level = 0.95)
## $expected.exceed
## [1] 24
##
## $actual.exceed
## [1] 127
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 235.4193
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 237.821
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest(alpha=0.05,actual=VN30.logret[1961:2454],VaR=VaR_95_his,conf.level = 0.95)
## $b
## [1] 1.427238
##
## $uLL
## [1] -286.3017
##
## $rLL
## [1] -298.148
##
## $LRp
## [1] 1.130242e-06
##
## $H0
## [1] "Duration Between Exceedances have no memory (Weibull b=1 = Exponential)"
##
## $Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_95_his)
## Accept H0 with LRstat is: 2.354077
VaRDurTest_Haas(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_95_his,conf=0.95)
## Reject H0 with LRstat is: 353.9289
ES_95_his <- DoGARCH_ES(y=forc, probability=0.05, WE=1960)
ES_95_his <- ES_95_his[1961:2454]
qplot(y = VaR_95_his , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES_95_his,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_95_his)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_his_95 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR_95_his,es=ES_95_his)
z1_his_95
## [1] -1.676681
z2_his_95 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR_95_his,es=ES_95_his,alpha=0.05)
z2_his_95
## [1] -12.76269
model.roll = ugarchroll(spec = model.spec , data = VN30.logret , n.start = 1960 , ,refit.window = 'moving')
VaR975_td = model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=shape, p=0.025)
qplot(y = VaR975_td , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR975_td)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with GARCH approach: ' , (sum(VN30.logret[1961:2454] < VaR975_td)) , sep = '')
## Number of exceptions with GARCH approach: 101
VaRTest(alpha=0.025,actual=VN30.logret[1961:2454],VaR=VaR975_td,conf.level = 0.95)
## $expected.exceed
## [1] 12
##
## $actual.exceed
## [1] 101
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 264.6171
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 264.6354
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR975_td)
## Accept H0 with LRstat is: 0.01535075
VaRDurTest_Haas(alpha=0.025,actual=as.numeric(VN30.logret[1961:2454]),var=VaR975_td,conf=0.95)
## Reject H0 with LRstat is: 376.6435
ES975_td = (1/(1-0.025))*model.roll@forecast$density[,'Sigma']*integrate(my_function,lower = 0.0000001,upper = 0.025)$value*(1/(shape-1))*(shape+(qt(p=0.025,df=shape))^2)
qplot(y = VaR975_td , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES975_td,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR975_td)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_975 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR975_td,es=ES975_td)
z1_975
## [1] 0.6690893
z2_975 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR975_td,es=ES975_td,alpha=0.05)
z2_975
## [1] -0.3531164
VaR_975_his <- DoGARCH_VaR(y=forc, probability=0.025, WE=1960)
VaR_975_his <- VaR_975_his[1961:2454]
qplot(y = VaR_975_his , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_975_his)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with His approach: ' , (sum(VN30.logret[1961:2454] < VaR_975_his)) , sep = '')
## Number of exceptions with His approach: 104
VaRTest(alpha=0.025,actual=VN30.logret[1961:2454],VaR=VaR_975_his,conf.level = 0.95)
## $expected.exceed
## [1] 12
##
## $actual.exceed
## [1] 104
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 278.5575
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 278.781
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_975_his)
## Accept H0 with LRstat is: 0.212415
VaRDurTest_Haas(alpha=0.025,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_975_his,conf=0.95)
## Reject H0 with LRstat is: 387.5602
ES_975_his <- DoGARCH_ES(y=forc, probability=0.025, WE=1960)
ES_975_his <- ES_975_his[1961:2454]
qplot(y = VaR_975_his , x = 1:494, geom = 'line') +
geom_line(aes(x=1:494,y=ES_975_his,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_975_his)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_his_975 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR_975_his,es=ES_975_his)
z1_his_975
## [1] -1.018027
z2_his_975 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR_975_his,es=ES_975_his,alpha=0.05)
z2_his_975
## [1] -7.496957
model.roll = ugarchroll(spec = model.spec , data = VN30.logret , n.start = 1960 , ,refit.window = 'moving')
VaR99_td = model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=shape, p=0.01)
qplot(y = VaR99_td , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR99_td)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with GARCH approach: ' , (sum(VN30.logret[1961:2454] < VaR99_td)) , sep = '')
## Number of exceptions with GARCH approach: 81
VaRTest(alpha=0.01,actual=VN30.logret[1961:2454],VaR=VaR99_td,conf.level = 0.95)
## $expected.exceed
## [1] 4
##
## $actual.exceed
## [1] 81
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 313.5025
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 313.8663
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR99_td)
## Accept H0 with LRstat is: 0.3745313
VaRDurTest_Haas(alpha=0.01,actual=as.numeric(VN30.logret[1961:2454]),var=VaR99_td,conf=0.95)
## Reject H0 with LRstat is: 413.1763
alpha=0.01
actual=as.numeric(VN30.logret[1961:2454])
var=VaR99_td
vec <- c(0)
stat <- c()
for (i in 1:length(var))
{
if (actual[i]<var[i])
{
vec <- cbind(vec,i)
}
}
for (j in 1:length(vec))
{
stat[j] <- -2*log((alpha*(1-alpha)^(vec[j+1]-vec[j]-1))/
((1/(vec[j+1]-vec[j]))*(1-(1/(vec[j+1]-vec[j])))^(vec[j+1]-vec[j]-1)))
}
z <- sum(na.omit(stat))
ES99_td = (1/(1-0.01))*model.roll@forecast$density[,'Sigma']*integrate(my_function,lower = 0.0000001,upper = 0.01)$value*(1/(shape-1))*(shape+(qt(p=0.01,df=shape))^2)
qplot(y = VaR99_td , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES99_td,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR99_td)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_99 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR99_td,es=ES99_td)
z1_99
## [1] 0.7443275
z2_99<- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR99_td,es=ES99_td,alpha=0.05)
z2_99
## [1] 0.1615597
VaR_99_his <- DoGARCH_VaR(y=forc, probability=0.01, WE=1960)
VaR_99_his <- VaR_99_his[1961:2454]
qplot(y = VaR_99_his , x = 1:494 , geom = 'line') +
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_99_his)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
cat('Number of exceptions with His approach: ' , (sum(VN30.logret[1961:2454] < VaR_99_his)) , sep = '')
## Number of exceptions with His approach: 73
VaRTest(alpha=0.01,actual=VN30.logret[1961:2454],VaR=VaR_99_his,conf.level = 0.95)
## $expected.exceed
## [1] 4
##
## $actual.exceed
## [1] 73
##
## $uc.H0
## [1] "Correct Exceedances"
##
## $uc.LRstat
## [1] 267.0161
##
## $uc.critical
## [1] 3.841459
##
## $uc.LRp
## [1] 0
##
## $uc.Decision
## [1] "Reject H0"
##
## $cc.H0
## [1] "Correct Exceedances & Independent"
##
## $cc.LRstat
## [1] 267.0733
##
## $cc.critical
## [1] 5.991465
##
## $cc.LRp
## [1] 0
##
## $cc.Decision
## [1] "Reject H0"
VaRDurTest_Chist(alpha=0.05,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_99_his)
## Accept H0 with LRstat is: 0.05357649
VaRDurTest_Haas(alpha=0.01,actual=as.numeric(VN30.logret[1961:2454]),var=VaR_99_his,conf=0.95)
## Reject H0 with LRstat is: 348.5173
ES_99_his <- DoGARCH_ES(y=forc, probability=0.01, WE=1960)
ES_99_his <- ES_99_his[1961:2454]
qplot(y = VaR_99_his , x = 1:494 , geom = 'line') +
geom_line(aes(x=1:494,y=ES_99_his,color='green'))+
geom_point(aes(x = 1:494 , y = VN30.logret[1961:2454] , color = as.factor(VN30.logret[1961:2454] < VaR_99_his)) , size = 2) + scale_color_manual(values = c('gray' , 'green','red')) +
labs(y = 'Daily Log Returns' , x = 'Test set Observation') + theme_light() +
theme(legend.position = 'none')
z1_his_99 <- as1(x=as.numeric(VN30.logret[1961:2454]),var=VaR_99_his,es=ES_99_his)
z1_his_99
## [1] -0.5495941
z2_his_99 <- as2(x=as.numeric(VN30.logret[1961:2454]),var=VaR_99_his,es=ES_99_his,alpha=0.01)
z2_his_99
## [1] -21.89886