Volatilidad y Value-at-Risk (VaR) con GARCH: selección y prueba de modelo - Quant Risk Management Workshop
S.B.
Previamente
- Volatilidad y Value-at-Risk (VaR) con GARCH
- Modelos de Volatilidad: Introducción rápida
- R-configuración
Anteriormente, fue definido y calibrado el modelo básico GARCH(1,1) con choques Normales. Se introdujo el concepto de simulación MC, pero existen miles de modelos posibles. ¿Como podemos tener la certeza que hemos elegido apropiadamente?
- ¿Qué tan aproximada a la volatilidad empírica se encuentra la función generadora (PGF) modelada?
- Sobre la distribución de errores, ¿es la mas adecuada?
- ¿Como reacciona a las noticias positivas vs. negativas.
- ¿De que forma acomoda la característica de saltos?
** En esta clase: Selección automatizada del modelo y sus paramétros. Test de modelo. **
Pre-procesamiento
# la función: https://rpubs.com/whitenoise1/969634
# la misma se puede pegar dentro del documento
source("~/My Drive/R/R codes/_github/functions/GARCH_msperlin/GARCH_msperlin_functions.r")
# leer desde el directorio donde se encuentra seteada la sesión de R.
# getwd() para ver. setwd() para asignar.
# source("GARCH_msperlin_functions.r")
# Si desea correr el analísis sobre una compañia publica.
# symbol = "^GSPC"
# ticker = getSymbols(symbol, src = 'yahoo', from = '2010-01-01', auto.assign = FALSE)
# Descarga la serie diaria GOBIXDR y realiza el pre-procesamiento.
symbol = "GOBIX"
ticker <- read.csv(url("https://www.bvrd.com.do/indice/Data/GobixDataIRP.csv"))
ticker <- ticker[,1:2]
colnames(ticker)[1] <- "fecha"
ticker$fecha <- lubridate::mdy(ticker[,1])
ticker <- as.xts(ticker[,-1], order.by=ticker$fecha)
colnames(ticker)[1] <- "Close.GOBIX"
# IDEA:
# 1. Divide la muestra en dos conjuntos.
# a. Primero, para calibrar el modelo mas apropiado (2014-2021).
# b. Segundo, realizar la prueba "Out-of-Sample".
# 2. Realizar una predicción para el año 2022 completo.
# 3. Interpretar el performance del GOBIX, sobre la proyección de precio construido a partir de la banda de volatilidad, estimada con el proceso GARCH.
# Prueba Out of Sample
H1.oos <- ticker["2022::"]
# Traininig: In the Sample
ticker <- ticker["2014::2021"]
ticker$daily.vol <- diff(log(Cl(ticker)))
ticker <- na.omit(ticker)Encontrar el modelo mas apropiado, sus paramétros.
tab_out <- do_arch_test(x = ticker$daily.vol, max_lag = 5)
max_lag_AR <- 1 # m - parametro AR
max_lag_MA <- 1 # n - parametro MA
max_lag_ARCH <- 1 # p - parametro ARCH
max_lag_GARCH <- 1 # q - parametro GARCH
dist_to_use <- c('norm','std','sstd','jsu') # ver rugarch::ugarchspecs
models_to_estimate <- c('sGARCH', 'eGARCH', 'gjrGARCH') # see rugarch::rugarchspec help for more
out <- find_best_arch_model(x = ticker$daily.vol,
type_models = models_to_estimate,
dist_to_use = dist_to_use,
max_lag_AR = max_lag_AR,
max_lag_MA = max_lag_MA,
max_lag_ARCH = max_lag_ARCH,
max_lag_GARCH = max_lag_GARCH)
# tabla resumen con resultados
tab_out <- out$tab_out
print(tab_out)## # A tibble: 48 × 9
## lag_ar lag_ma lag_arch lag_garch AIC BIC type_model type_dist model_name
## <int> <int> <int> <int> <dbl> <dbl> <chr> <chr> <chr>
## 1 0 0 1 1 -9.37 -9.36 sGARCH norm ARMA(0,0)+…
## 2 0 1 1 1 -9.37 -9.36 sGARCH norm ARMA(0,1)+…
## 3 1 0 1 1 -9.37 -9.36 sGARCH norm ARMA(1,0)+…
## 4 1 1 1 1 NA NA sGARCH norm ARMA(1,1)+…
## 5 0 0 1 1 NA NA sGARCH std ARMA(0,0)+…
## 6 0 1 1 1 NA NA sGARCH std ARMA(0,1)+…
## 7 1 0 1 1 NA NA sGARCH std ARMA(1,0)+…
## 8 1 1 1 1 -9.63 -9.61 sGARCH std ARMA(1,1)+…
## 9 0 0 1 1 NA NA sGARCH sstd ARMA(0,0)+…
## 10 0 1 1 1 -9.63 -9.61 sGARCH sstd ARMA(0,1)+…
## # … with 38 more rows
## # ℹ Use `print(n = ...)` to see more rowsmodels_names <- unique(tab_out$model_name)
best_models <- c(tab_out$model_name[which.min(tab_out$AIC)],
tab_out$model_name[which.min(tab_out$BIC)])
print(best_models)## [1] "ARMA(0,0)+eGARCH(1,1) jsu" "ARMA(0,0)+eGARCH(1,1) jsu"# Define el modelo que minimiza el BIC
best_spec = ugarchspec(variance.model = list(model = out$best_bic$type_model,
garchOrder = c(out$best_bic$lag_arch,
out$best_bic$lag_garch)),
mean.model = list(armaOrder = c(out$best_bic$lag_ar,
out$best_bic$lag_ma)),
distribution = out$best_bic$type_dist)El modelo mas apropiado es:
- Según el criterio Akaike (BIC): ARMA(0,0)+eGARCH(1,1) jsu.
- Según el criterio Bayesiano (BIC): ARMA(0,0)+eGARCH(1,1) jsu.
Fit
Procedemos con el fit sobre el mejor modelo.
# Implementa el modelo definido
my_best_garch <- ugarchfit(spec = best_spec,
data = ticker$daily.vol)
my_best_garch##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : jsu
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000053 0.000045 1.1690 0.242415
## omega -0.678756 0.062880 -10.7944 0.000000
## alpha1 0.048373 0.020974 2.3064 0.021090
## beta1 0.944006 0.005170 182.6038 0.000000
## gamma1 0.233898 0.035801 6.5333 0.000000
## skew 0.090578 0.043882 2.0642 0.039003
## shape 1.047690 0.049035 21.3663 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000053 0.000055 0.96823 0.332930
## omega -0.678756 0.032993 -20.57265 0.000000
## alpha1 0.048373 0.024189 1.99975 0.045527
## beta1 0.944006 0.002674 352.99885 0.000000
## gamma1 0.233898 0.038641 6.05309 0.000000
## skew 0.090578 0.048661 1.86142 0.062684
## shape 1.047690 0.070244 14.91510 0.000000
##
## LogLikelihood : 9683.306
##
## Information Criteria
## ------------------------------------
##
## Akaike -9.6426
## Bayes -9.6230
## Shibata -9.6426
## Hannan-Quinn -9.6354
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.08407 0.771852
## Lag[2*(p+q)+(p+q)-1][2] 0.22324 0.840065
## Lag[4*(p+q)+(p+q)-1][5] 12.35397 0.002248
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.02022 0.8869
## Lag[2*(p+q)+(p+q)-1][5] 1.21439 0.8095
## Lag[4*(p+q)+(p+q)-1][9] 2.50723 0.8367
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.005356 0.500 2.000 0.9417
## ARCH Lag[5] 2.313865 1.440 1.667 0.4059
## ARCH Lag[7] 2.848208 2.315 1.543 0.5427
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 5.6987
## Individual Statistics:
## mu 0.3298
## omega 1.3961
## alpha1 0.1186
## beta1 1.3529
## gamma1 0.1967
## skew 0.4464
## shape 4.3677
##
## 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.6668 0.5050
## Negative Sign Bias 0.2122 0.8319
## Positive Sign Bias 0.7453 0.4562
## Joint Effect 0.6842 0.8769
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 52.30 5.961e-05
## 2 30 68.23 5.272e-05
## 3 40 76.69 2.965e-04
## 4 50 117.94 1.297e-07
##
##
## Elapsed time : 0.7364221plot(my_best_garch, which="all")##
## please wait...calculating quantiles...
# model params
infocriteria(my_best_garch) ##
## Akaike -9.642558
## Bayes -9.623011
## Shibata -9.642582
## Hannan-Quinn -9.635382# coefficients
coef(my_best_garch) ## mu omega alpha1 beta1 gamma1
## 5.284391e-05 -6.787563e-01 4.837295e-02 9.440064e-01 2.338979e-01
## skew shape
## 9.057834e-02 1.047690e+00# unconditional mean in mean equation
uncmean(my_best_garch) ## [1] 5.284391e-05# unconditional variance: omega/(alpha1 + beta1)
uncvariance(my_best_garch) ## [1] 5.438392e-06# persistence: alpha1 + beta1
persistence(my_best_garch) ## [1] 0.9440064# half-life: no. días donde la volatilidad revierte 50%
halflife(my_best_garch) ## [1] 12.02914# predicción 10 días hacia adelante
garch.forecast = ugarchforecast(my_best_garch, n.ahead = 10)
garch.forecast##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: eGARCH
## Horizon: 10
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=2021-12-31]:
## Series Sigma
## T+1 5.284e-05 0.003838
## T+2 5.284e-05 0.003732
## T+3 5.284e-05 0.003635
## T+4 5.284e-05 0.003546
## T+5 5.284e-05 0.003464
## T+6 5.284e-05 0.003388
## T+7 5.284e-05 0.003318
## T+8 5.284e-05 0.003253
## T+9 5.284e-05 0.003193
## T+10 5.284e-05 0.003137# visualizar el forecast
par(mfrow=c(2,1))
plot(garch.forecast, which=1)
plot(garch.forecast, which=3)
par(mfrow=c(1,1))Algunos elementos importantes de este ejercicio:
La predicción de los próximos 10 días indica:
\[\begin{align*} E_{retorno_t+1}: 0.0053\\ E_{retorno_t+2}: 0.0053\\ \\ \sigma_{diario_t+1}: 0.384 \\ \sigma_{anualizado_t+1}: 6.093 \\ \end{align*}\]
La volatilidad incondicional diaria es de 0.002. Traducido a la volatilida anual y expresada en términos porcentuales, esto representa 3.702%.
La persistencia persistencia del modelo es de 0.944 < 1, indicando que el modelo es estacionario, lo cual proporciona convergencia y estabilidad al proceso modelado.
La velocidad a la cual el modelo reviere a la media indica que el 50% de la subida se normaliza dentro de 12 días.
Simulación Monte-Carlo
- Con el mejor modelo identificado se procede a realizar 5,000 simulaciones.
days.ahead = 365
n.sim <- 5000
garch.sim <- matrix(nrow = days.ahead, ncol=n.sim)
set.seed(123)
for(i in 1:n.sim){
p.sim = ugarchsim(my_best_garch, n.sim=days.ahead, startMethod="sample")
garch.sim[,i] <- p.sim@simulation$seriesSim
}
garch.sim <- as.data.frame(garch.sim)
# Determinar los percentiles correspondientes
garch.sim$Q25 <- NA
garch.sim$Q05 <- NA
garch.sim$Q01 <- NA
garch.sim$Q75 <- NA
garch.sim$Q95 <- NA
garch.sim$Q99 <- NA
garch.sim$Q01 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.01)}, MARGIN = 1)
garch.sim$Q05 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.05)}, MARGIN = 1)
garch.sim$Q25 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.25)}, MARGIN = 1)
garch.sim$Q75 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.75)}, MARGIN = 1)
garch.sim$Q95 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.95)}, MARGIN = 1)
garch.sim$Q99 <- apply(garch.sim[,2:(ncol(garch.sim)-6)], FUN = function(x){quantile(na.omit(x),0.99)}, MARGIN = 1)Estimación probabilidades en el precio futuro
Procedemos con el fit sobre el mejor modelo.
# PRICE LEVEL FORECAST ----
# Mapeo de predicciones de volatilidad por vía de la simulación en el precio
sim.dates <- seq(lubridate::ymd(last(index(ticker))+1), last(index(ticker)) + (days.ahead), by = "days")
df <- cbind(as.data.frame(sim.dates), garch.sim)
df <- as.xts(df[,-1], order.by=df$sim.dates)
df <- as.numeric(tail(Cl(ticker),1)) * cumprod(1 + df)
df.combined <- rbind(Cl(ticker),Cl(H1.oos))
df.combined <- merge(df.combined,df[,1:(ncol(df)-1)]) # do not merge the empty column
last.day <- index(tail(Cl(ticker),1))
df.combined$Q01 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-1)], FUN = function(x){quantile(na.omit(x),0.01)}, MARGIN = 1))
df.combined$Q05 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-2)], FUN = function(x){quantile(na.omit(x),0.05)}, MARGIN = 1))
df.combined$Q25 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-3)], FUN = function(x){quantile(na.omit(x),0.25)}, MARGIN = 1))
df.combined$Q75 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-4)], FUN = function(x){quantile(na.omit(x),0.75)}, MARGIN = 1))
df.combined$Q95 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-5)], FUN = function(x){quantile(na.omit(x),0.95)}, MARGIN = 1))
df.combined$Q99 <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-6)], FUN = function(x){quantile(na.omit(x),0.99)}, MARGIN = 1))
df.combined$mean <- as.numeric(apply(df.combined[,2:(ncol(df.combined)-7)], FUN = function(x){mean(na.omit(x))}, MARGIN = 1))
tail(df.combined[,c('Q01','Q05','Q25','mean','Q75','Q95','Q99')],15)## Q01 Q05 Q25 mean Q75 Q95 Q99
## 2022-12-17 114.3657 118.6377 124.2506 128.4556 132.2650 139.3950 146.5314
## 2022-12-18 114.2925 118.6773 124.2613 128.4543 132.2649 139.3784 146.5258
## 2022-12-19 114.4803 118.6622 124.2964 128.4593 132.2951 139.3428 146.4400
## 2022-12-20 114.3552 118.6145 124.2959 128.4747 132.3064 139.4625 146.4942
## 2022-12-21 114.5616 118.7044 124.3096 128.4791 132.2381 139.5193 146.4996
## 2022-12-22 114.5373 118.6696 124.2887 128.4845 132.2475 139.5070 146.6776
## 2022-12-23 114.3914 118.4832 124.2768 128.4867 132.2663 139.5013 146.5134
## 2022-12-24 114.3482 118.5356 124.2663 128.5009 132.2862 139.4876 146.4236
## 2022-12-25 114.3427 118.5275 124.2577 128.5077 132.4178 139.5625 146.8116
## 2022-12-26 114.3470 118.5529 124.2707 128.5162 132.3807 139.5432 146.5604
## 2022-12-27 114.3672 118.6014 124.2335 128.5185 132.4411 139.5403 146.7051
## 2022-12-28 114.3707 118.6528 124.2318 128.5226 132.4361 139.5939 146.7157
## 2022-12-29 114.3971 118.6476 124.2471 128.5325 132.3995 139.6563 146.7794
## 2022-12-30 114.3645 118.6450 124.2687 128.5460 132.4519 139.6417 146.7196
## 2022-12-31 114.5162 118.6522 124.2766 128.5593 132.4473 139.6673 146.9332#Visualizacion
par(mfrow=c(1,1))
garch.sim.plot <- plot(Cl(df.combined["2019::"]), ylim=c(min(na.omit(df.combined["2019::",1]))*0.6,max(na.omit(df.combined["2020::",1]))*1.3),
main=paste0("Simulated ",
out$best_bic$type_model, "(", out$best_bic$lag_arch, "," ,out$best_bic$lag_garch, ") + ARMA(",out$best_bic$lag_ar, "," ,out$best_bic$lag_ma,") + ", out$best_bic$type_dist, " distribution"), grid.col=NA)
garch.sim.plot <- lines(df.combined[,2:(ncol(df.combined)-6)], col=alpha("grey",0.5), on=1, lty=2, lwd=0.75)
garch.sim.plot <- lines(df.combined[,'Q01'], col="red", on=1, lty=1, lwd=1)
garch.sim.plot <- lines(df.combined[,'Q05'], col="red", on=1, lty=2, lwd=1.5)
garch.sim.plot <- lines(df.combined[,'Q25'], col="red", on=1, lty=3, lwd=1.5)
garch.sim.plot <- lines(df.combined[,'mean'], col="black", on=1, lty=2, lwd=1.5)
garch.sim.plot <- lines(df.combined[,'Q75'], col="blue", on=1, lty=3, lwd=1.5)
garch.sim.plot <- lines(df.combined[,'Q95'], col="blue", on=1, lty=2, lwd=1.5)
garch.sim.plot <- lines(df.combined[,'Q99'], col="blue", on=1, lty=1, lwd=1)
garch.sim.plot <- lines(Cl(df.combined["2022::"]), col="black", on=1, lwd=1.25)
garch.sim.plot <- points(tail(df.combined[,'mean'],1), col="red", pch=16)
garch.sim.plot
Validación de Modelo
# FUNCIONES
# funcion auxiliar para leer el nombre
convert.to.char <- function(v1) { deparse(substitute(v1)) }
# sumariza el resultado del backtest
rolling.var.backtest <- function(rolling.garch.object, symbol, return.series, n){
require(xts)
require(dplyr)
require(lubridate)
garch.object.name <- convert.to.char(rolling.garch.object)
VaR <- paste0(garch.object.name,"@forecast$VaR")
VaR <- eval(parse(text = VaR))
#forecast.var <- paste0(garch.object.name,"@forecast$VaR")
realized.var.name <- paste0(garch.object.name,"@forecast$VaR$realized")
#assign("realized.var", realized.var.name)
realized.var <- eval(parse(text = realized.var.name))
var05.breach <- which(rolling.garch.object@forecast$VaR$realized < as.numeric(quantile(return.series, 0.05)))
var01.breach <- which(rolling.garch.object@forecast$VaR$realized < as.numeric(quantile(return.series, 0.01)))
var95.breach <- which(rolling.garch.object@forecast$VaR$realized > as.numeric(quantile(return.series, 0.95)))
var99.breach <- which(rolling.garch.object@forecast$VaR$realized > as.numeric(quantile(return.series, 0.99)))
rows <- index(rolling.garch.object@forecast$VaR[var95.breach,'realized'])
forecast.VaR <- as.data.frame(rolling.garch.object@forecast$VaR)
forecast.VaR$var01.rank <- ntile(forecast.VaR$`alpha(1%)`,100)
forecast.VaR$var99.rank <- ntile(forecast.VaR$`alpha(99%)`,100)
forecast.VaR$var05.breach <- ifelse(forecast.VaR$realized < rolling.garch.object@forecast$VaR$`alpha(5%)`,1,0)
forecast.VaR$var01.breach <- ifelse(forecast.VaR$realized < rolling.garch.object@forecast$VaR$`alpha(1%)`,1,0)
forecast.VaR$var95.breach <- ifelse(forecast.VaR$realized > rolling.garch.object@forecast$VaR$`alpha(95%)`,1,0)
forecast.VaR$var99.breach <- ifelse(forecast.VaR$realized > rolling.garch.object@forecast$VaR$`alpha(99%)`,1,0)
forecast.VaR$nextday.ret <- lead(forecast.VaR$realized)
var.breaches <- list()
# breaches are identified in cols: 8-11
for(i in 1:4){
breach.days <- forecast.VaR[forecast.VaR[,i+7] == 1,]
breach.return <- forecast.VaR[forecast.VaR[,i+7] == 1,'nextday.ret']
breach.days$return <- cumprod(1 + na.omit(breach.return))
forecast.VaR <- as.xts(tail(forecast.VaR,n))
forecast.density <- as.data.frame(rolling.garch.object@forecast$density)
forecast.density <- as.xts(tail(forecast.density,n))
forecast.density$price <- as.numeric(tail(return.series,n))
# see consecutive breaches
var.breaches[[i]] <- as.data.frame(breach.days)
var.breaches[[i]]$date <- ymd(rownames(var.breaches[[i]]))
var.breaches[[i]]$days.between <- NA
#for(j in 2:nrow(var.breaches[[i]])){
# var.breaches[[i]]$days.between[j] <- var.breaches[[i]]$date[j] - var.breaches[[i]]$date[j-1]
#}
}
# upside breaches
upside.breaking.VaR <- forecast.VaR[forecast.VaR$realized > 0,]
upside.breaking.VaR$var99.breach.levels <- upside.breaking.VaR$realized - upside.breaking.VaR$`alpha(99%)`
mean.upside.breaking <- mean(upside.breaking.VaR[upside.breaking.VaR$var99.breach.levels > 0,'var99.breach.levels'])
# downside breaches
downside.breaking.VaR <- forecast.VaR[forecast.VaR$realized < 0,]
downside.breaking.VaR$var01.breach.levels <- downside.breaking.VaR$realized - downside.breaking.VaR$`alpha(1%)`
mean.downside.breaking <- mean(downside.breaking.VaR[downside.breaking.VaR$var01.breach.levels < 0,'var01.breach.levels'])
VaR.backtest <- list(forecast.VaR,
var.breaches,
upside.breaking.VaR,
downside.breaking.VaR)
return(VaR.backtest)
}Recalibración y Prueba
# CALIBRACION CORRIDA Y VALIDACION ----
# Idea: tal vez el modelo debería ser recalibrado de vez en cuando
# Ejecuta la predicción de 1-día hacía adelante (rolling GARCH) con recalibración cada 30 días
garch.roll = ugarchroll(best_spec, ticker$daily.vol, n.ahead=1,
forecast.length = 1000, solver = "hybrid",
refit.every=30, refit.window="moving", VaR.alpha=c(0.01, 0.05, 0.95, 0.99))
# Prueba de modelo
report(garch.roll, type="VaR")## VaR Backtest Report
## ===========================================
## Model: eGARCH-jsu
## Backtest Length: 1000
## Data:
##
## ==========================================
## alpha: 1%
## Expected Exceed: 10
## Actual VaR Exceed: 16
## Actual %: 1.6%
##
## Unconditional Coverage (Kupiec)
## Null-Hypothesis: Correct Exceedances
## LR.uc Statistic: 3.077
## LR.uc Critical: 3.841
## LR.uc p-value: 0.079
## Reject Null: NO
##
## Conditional Coverage (Christoffersen)
## Null-Hypothesis: Correct Exceedances and
## Independence of Failures
## LR.cc Statistic: 8.212
## LR.cc Critical: 5.991
## LR.cc p-value: 0.016
## Reject Null: YESreport(garch.roll, type="fpm")##
## GARCH Roll Mean Forecast Performance Measures
## ---------------------------------------------
## Model : eGARCH
## No.Refits : 34
## No.Forecasts: 1000
##
## Stats
## MSE 7.148e-06
## MAE 1.693e-03
## DAC 5.220e-01a <- rolling.var.backtest(rolling.garch.object=garch.roll, symbol=symbol, return.series=ticker$daily.vol, n=1000)
# static VaR. LT
asset1.ret <- ticker$daily.vol
var05.breach <- which(garch.roll@forecast$VaR$realized < as.numeric(quantile(asset1.ret, 0.05)))
var01.breach <- which(garch.roll@forecast$VaR$realized < as.numeric(quantile(asset1.ret, 0.01)))
var95.breach <- which(garch.roll@forecast$VaR$realized > as.numeric(quantile(asset1.ret, 0.95)))
var99.breach <- which(garch.roll@forecast$VaR$realized > as.numeric(quantile(asset1.ret, 0.99)))
rows <- index(garch.roll@forecast$VaR[var95.breach,'realized'])
forecast.VaR <- as.data.frame(garch.roll@forecast$VaR)
forecast.VaR$var05.breach <- ifelse(forecast.VaR$realized < garch.roll@forecast$VaR$`alpha(5%)`,1,0)
forecast.VaR$var01.breach <- ifelse(forecast.VaR$realized < garch.roll@forecast$VaR$`alpha(1%)`,1,0)
forecast.VaR$var95.breach <- ifelse(forecast.VaR$realized > garch.roll@forecast$VaR$`alpha(95%)`,1,0)
forecast.VaR$var99.breach <- ifelse(forecast.VaR$realized > garch.roll@forecast$VaR$`alpha(99%)`,1,0)
forecast.VaR$nextday.ret <- lead(forecast.VaR$realized)
breach.days <- forecast.VaR[forecast.VaR$var95.breach == 1,]
breach.return <- forecast.VaR[forecast.VaR$var95.breach == 1,'nextday.ret']
breach.return <- cumprod(1 + na.omit(breach.return))
forecast.VaR <- as.xts(forecast.VaR)
forecast.density <- as.data.frame(garch.roll@forecast$density)
forecast.density <- as.xts(forecast.density)
forecast.density$price <- as.numeric(tail(asset1.ret,1000))
# see consecutive breaches
forecast.VaR <- a[[1]]
forecast.VaR1 <- as.data.frame(a[[1]][a[[1]]$var95.breach==1,])
forecast.VaR1$date <- ymd(rownames(forecast.VaR1))
forecast.VaR1$days.between <- NA
for(i in 2:nrow(forecast.VaR1)){
forecast.VaR1$days.between[i] <- forecast.VaR1$date[i] - forecast.VaR1$date[i-1]
}
# Plot 1. VaR Only
par(mfrow=c(1,1), mar=c(3,3,3,3))
rolling.plot <- plot(forecast.VaR$realized, ylim=c(min(forecast.VaR$`alpha(1%)` -0.01), max(forecast.VaR$`alpha(99%)` +0.01)), main=paste(symbol,"t+1 Rolling eGARCH VaR Estimation w/ 30-Day Refit"), type="h", col="black", lwd=0.50)
rolling.plot <- lines(forecast.VaR$`alpha(5%)`, on=1, col="orange", lty=1, lwd=0.75)
rolling.plot <- lines(forecast.VaR$`alpha(1%)`, on=1, col="red", lty=1, lwd=1)
rolling.plot <- lines(forecast.VaR$`alpha(95%)`, on=1, col="orange", lty=1, lwd=0.75)
rolling.plot <- lines(forecast.VaR$`alpha(99%)`, on=1, col="red", lty=1, lwd=1)
rolling.plot <- points(forecast.VaR[forecast.VaR$var99.breach==1,5], col="red", cex=1.25, pch=16)
rolling.plot <- points(forecast.VaR[forecast.VaR$var01.breach==1,1], col="red", cex=1.25, pch=16)
rolling.plot <- addLegend("bottomleft", legend.names=c(paste("05pctile: Current", round(last(forecast.VaR$`alpha(5%)`)*100,2),"%"),
paste("01pctile: Current", round(last(forecast.VaR$`alpha(1%)`)*100,2),"%")), lty=c(1,1), col=c("orange","red"), bty="n", y.intersp = 0.75)
rolling.plot <- addLegend("topleft", legend.names=c(paste("99pctile: Current", round(last(forecast.VaR$`alpha(99%)`)*100,2),"%"),
paste("95pctile: Current", round(last(forecast.VaR$`alpha(95%)`)*100,2),"%")), lty=c(1,1), col=c("red","orange"), bty="n", y.intersp = 0.75)
plot(rolling.plot)
# Plot 2. VaR + Sigma + Mu
rolling.plot <- plot(forecast.VaR$realized, ylim=c(min(forecast.VaR$`alpha(1%)` -0.01), max(forecast.VaR$`alpha(99%)` +0.01)), main=paste(symbol,"t+1 Rolling eGARCH VaR Estimation w/ 30-Day Refit"), type="h", col="black", lwd=0.50)
forecast.density.plot <- lines(forecast.VaR$`alpha(5%)`, on=1, col="orange", lty=1, lwd=0.75)
forecast.density.plot <- lines(forecast.VaR$`alpha(1%)`, on=1, col="red", lty=1, lwd=1)
forecast.density.plot <- lines(forecast.VaR$`alpha(95%)`, on=1, col="orange", lty=1, lwd=0.75)
forecast.density.plot <- lines(forecast.VaR$`alpha(99%)`, on=1, col="red", lty=1, lwd=1)
forecast.density.plot <- points(forecast.VaR[forecast.VaR$var99.breach==1,5], col="red", cex=1.25, pch=16)
forecast.density.plot <- points(forecast.VaR[forecast.VaR$var01.breach==1,1], col="red", cex=1.25, pch=16)
forecast.density.plot <- addLegend("bottomleft", legend.names=c(paste("05pctile: Current", round(last(forecast.VaR$`alpha(5%)`)*100,2),"%"),
paste("01pctile: Current", round(last(forecast.VaR$`alpha(1%)`)*100,2),"%")), lty=c(1,1), col=c("orange","red"), bty="n", y.intersp = 0.75)
forecast.density.plot <- addLegend("topleft", legend.names=c(paste("99pctile: Current", round(last(forecast.VaR$`alpha(99%)`)*100,2),"%"),
paste("95pctile: Current", round(last(forecast.VaR$`alpha(95%)`)*100,2),"%")), lty=c(1,1), col=c("red","orange"), bty="n", y.intersp = 0.75)
forecast.density.plot <- lines(abs(forecast.density[,6]), on=NA, main="Sigma", col="grey", lwd=0.5)
forecast.density.plot <- lines(forecast.density[,2], on=2, col="red", lwd=2)
forecast.density.plot <- lines(forecast.density[,1], on=NA, main="Mu", col="black", lwd=1)
plot(forecast.density.plot)
Conclusiones
El mejor modelo es ARMA(0,0)+eGARCH(1,1) jsu. Este permite la asimetría en choques. Pasa la prueba de modelo incondicional de Kupiec, pero no logra determinar que los excesos son independientes (Christoffersen). ¿La dependencia en los excesos registrados en rebases de la banda de volatilidad, es producto de mala suerte o señala que el modelo sea inadecuado?