Volatilidad y Value-at-Risk (VaR) con GARCH - Quant Risk Management Workshop

Función Generadora del Proceso de Volatilidad

Existen muchas formas para estimar la volatilidad. El principal problema al cual se enfrentan los Quant Risk Managers al momento de realizar estimaciones de la volatilidad: determinar en que forma deberían ponderar el pasado. Debido al fenómeno de clustering y la persistencia, el comportamiento de la volatilidad se mantiene oscilando entre niveles altos y bajos, pudiendo permanecer por tiempo indeterminado en cualquier estado.

La estimación del VaR histórico corre el riesgo de no funcionar por periodos prolongados, ya sea subestimando (estados turbulentos) o sobreestimando la volatilidad realizada (estados tranquilos). Este es el caso de la varianza incondicional: se otorga la misma ponderación a todas las observaciones (las observaciones de la semana de 25 al 29 de enero del 1968 aportan lo mismo a la estimación que las registradas en la semana anterior).
Tomando en cuenta la propiedad inherente de la “cola gorda” (eng: fat tail) en mercados financieros, la varianza histórica tiende a representar de manera desproporcional unas pocas observaciones extremas, pudiendo distorsionar la realidad representada en los modelos históricos.
Una solución factible de este problema consiste en definir la Función Generadora de Volatilidad (el modelo) con memoria limitada para fines de reducir la importancia de esas observaciones, y tomar en cuenta el resultado de las predicciones anteriores. Esto introduce el concepto de la varianza condicional en el tiempo y nivel.

Revisión teórica

EWMA

Exponentially Weighted Moving Average (EWMA) ha sido el predecesor de los modelos modernos de volatilidad realizada. Su principal innovación es condicionar la varianza al tiempo de ocurrencia. Su premisa principal: el pasado mas reciente tiene mayor importancia en determinar el futuro muy cercano.

En esencia: \[\begin{align*} volatilidad_{t} = \\ predicción_{t-1} * \lambda \\ + \\ noticia_{t-1,t} * (1-\lambda) \end{align*}\]

Formalmente: \[\begin{align*} \sigma^2 = \lambda \sigma_{t-1}^2 + (1-\lambda)R_{t-1,t}^2 \end{align*}\]

Donde, \(\sigma_{t-1}^2\) es la volatilidad estimada el día anterior, y \(R_{t-1,t}^2\) es el retorno (noticia) entre el último día y hoy. \(\lambda\) es in parámetro limitado al rango cerrado [0,1] que describe el decaimiento en la importancia como la función en el tiempo.

mas cercano a 0: mayor ponderación a observaciones recientes, mayor dependencia.
mas cercano a 1: menor ponderación a observaciones recientes, mayor aleatoriedad.

La principal ventaja de este tipo de modelos es que son faciles de implementar y computacionalmente convenientes. El modelo original de VaR (RiskMetrics) era un EWMA con \(\lambda\) = 0.94. El EWMA es un caso especial (y limitado) del GARCH.

GARCH

La familia de los modelos GARCH (Generalized AutoRegressive Conditional Heteroskedasticity), tal y como indica su nombre, permiten modelar los resultados de sus predicciones tomando en cuenta el comportamiento en el término del error (\(\epsilon\)) en el tiempo (término \(\lambda\), visto anteriormente), siguiendo un proceso autoregresivo con media móvil.

Existe una gran variedad de modelos GARCH: iGARCH, eGARCH, fGARCH, sGARCH, nGARCH, qGARCH, apGARCH, GJR-GARCH, T-GARCH, entre otros. Ver: Zivot, Bollerslev
El GARCH es menos restrictivo que el EWMA por lo que es mas preciso, sin embargo, dado que contiene una mayor cantidad de parámetros, es mas sensible al riesgo del modelo, dado que la optimización juega un rol mas importante.

\[\begin{align*} \sigma^2 = w + \alpha R_{t-1,t}^2 + \beta \sigma_{t-1,t}^2 \end{align*}\]

Donde, \(\alpha + \beta\) =< 1

La especificación del modelo GARCH(p,q) controla por el rezago en el orden autoregresivo de la varianza \(\sigma^2\) (q) y el término de error \(\epsilon\) proveniente de \(R_{t-1,t}^2\) (p).

La especificación del modelo se ejecuta de la siguiente forma:

Se estima el modelo autorregresivo (AR), utilizando el retorno diario.
Se computan autocorrelaciones con variables rezagos.
La hipotésis nula: Inexistencia de errores GARCH. El rechazo indica que esos errores existen en la varianza condicional.

En la práctica, esto se hace computando el autocorrelograma y observando la significancia estadistica en los rezagos individuales.

library(PerformanceAnalytics)
library(quantmod)
library(rugarch)
library(car)
library(FinTS)
library(scales)
library(data.table)

# Descargar la data
symbols <- c("^GSPC")
getSymbols(symbols, from ="2000-01-03", to = "2021-12-31")

## [1] "^GSPC"

GSPC.ret <- CalculateReturns(Cl(GSPC), method="log")
GSPC.ret <- GSPC.ret[-1,]
colnames(GSPC.ret) <- "GSPC"

# Estadísticas Descriptiva
table.Stats(GSPC.ret)

##                      GSPC
## Observations    5534.0000
## NAs                0.0000
## Minimum           -0.1277
## Quartile 1        -0.0047
## Median             0.0006
## Arithmetic Mean    0.0002
## Geometric Mean     0.0001
## Quartile 3         0.0058
## Maximum            0.1096
## SE Mean            0.0002
## LCL Mean (0.95)   -0.0001
## UCL Mean (0.95)    0.0005
## Variance           0.0002
## Stdev              0.0124
## Skewness          -0.4004
## Kurtosis          11.0561

# Plot
dataToPlot <- cbind(GSPC.ret, GSPC.ret^2, abs(GSPC.ret))
colnames(dataToPlot) <- c("Retornos", "Retornos^2", "abs(Retornos)")
plot.zoo(dataToPlot, main=paste(symbols,"Retornos Diarios"), col="blue")

# Computar el rezago del retorno de ayer
GSPC.ret$yesterday.ret <- data.table::shift(GSPC.ret$GSPC, n=1, type = "lag")
GSPC.ret <- na.omit(GSPC.ret)
lm.autocorr <- lm(coredata(GSPC.ret$GSPC) ~ coredata(GSPC.ret$yesterday.ret))

  # visualizar
  # crear gráficos 2x2
  par(mfrow=c(2,2))

  # plot 1: Autocorrelación
  plot(coredata(GSPC.ret$yesterday.ret), coredata(GSPC.ret$GSPC), col=scales::alpha("black",0.6), pch=16, cex=0.8, main="Autocorrelación de retornos", xlab="retorno t-1", ylab="retorno t")
  abline(lm.autocorr, col="red")
  legend("topleft", legend=c(paste0("Pendiente: ",round(summary(lm.autocorr)$coefficients[2,1],3)," ( t-value ",round(summary(lm.autocorr)$coefficients[2,3],3),")"),
                             paste0("R^2: ",round(summary(lm.autocorr)$r.squared,3))), bty="n")
  box(col="grey")
  
  # plot 2: serie histórica residuales
  plot(lm.autocorr$residuals, type="l", main="Residuales")
  
  # plot 3: distribución residuales
  hist(lm.autocorr$residuals, breaks=64, main="PDF", col="grey", freq = FALSE, xlab="residuales")
  rug(lm.autocorr$residuals)
  box(col="grey")
  lines(density(lm.autocorr$residuals, bw=0.004), col="red")
  
  # plot 4: QQ-plot: revisión normalidad
  car::qqPlot(lm.autocorr$residuals, main="QQ-Plot")

## [1] 5080 2206

  par(mfrow=c(1,1)) 
  
  summary(lm.autocorr)

## 
## Call:
## lm(formula = coredata(GSPC.ret$GSPC) ~ coredata(GSPC.ret$yesterday.ret))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.117897 -0.004927  0.000563  0.005477  0.107994 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       0.0002462  0.0001654   1.488    0.137    
## coredata(GSPC.ret$yesterday.ret) -0.1126219  0.0133484  -8.437   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0123 on 5531 degrees of freedom
## Multiple R-squared:  0.01271,    Adjusted R-squared:  0.01253 
## F-statistic: 71.19 on 1 and 5531 DF,  p-value: < 2.2e-16

# plot autocorrelaciones de retornos, retornos^2 and abs(retornos)
par(mfrow=c(3,1))
  acf(GSPC.ret$GSPC, main="^GSPC Retornos")
  acf(GSPC.ret$GSPC^2, main="^GSPC Retornos^2s")
  acf(abs(GSPC.ret$GSPC), main="^GSPC Abs Retornos")

par(mfrow=c(1,1)) 

# Testeando por efectos ARCH/GARCH en retornos
Box.test(coredata(GSPC.ret$GSPC^2), type="Ljung-Box", lag = 12)

## 
##  Box-Ljung test
## 
## data:  coredata(GSPC.ret$GSPC^2)
## X-squared = 6403.1, df = 12, p-value < 2.2e-16

# ArchTest() del FinTS package para Engle's LM test
ArchTest(GSPC.ret$GSPC)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  GSPC.ret$GSPC
## Chi-squared = 1660.3, df = 12, p-value < 2.2e-16

Modelo GARCH

Realizamos un fit de prueba, GARCH(1,1) Normal.

# Construir el modelo GARCH
spec = ugarchspec() # Default: GARCH(1,1)

# GJR GARCH acomoda a la asimetria en choques y saltos presentes en la volatilidad
# la distribución Standard Skewed
# spec = ugarchspec(variance.model=list(model="gjrGARCH", garchOrder=c(1,1)), 
#                    mean.model=list(armaOrder=c(0,0), include.mean=TRUE),  
#                    distribution.model="sstd")
fit = ugarchfit(data = GSPC.ret$GSPC, spec = spec)
fit

## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(1,0,1)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000637    0.000087   7.2982 0.000000
## ar1     0.743216    0.117393   6.3310 0.000000
## ma1    -0.793121    0.106879  -7.4207 0.000000
## omega   0.000002    0.000001   3.3060 0.000946
## alpha1  0.124888    0.009795  12.7497 0.000000
## beta1   0.857968    0.010203  84.0919 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000637    0.000086  7.43496 0.000000
## ar1     0.743216    0.167512  4.43679 0.000009
## ma1    -0.793121    0.154207 -5.14322 0.000000
## omega   0.000002    0.000004  0.57284 0.566754
## alpha1  0.124888    0.027744  4.50139 0.000007
## beta1   0.857968    0.040639 21.11217 0.000000
## 
## LogLikelihood : 17915.46 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.4737
## Bayes        -6.4665
## Shibata      -6.4737
## Hannan-Quinn -6.4712
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.004771 0.94493
## Lag[2*(p+q)+(p+q)-1][5]  3.963250 0.07301
## Lag[4*(p+q)+(p+q)-1][9]  6.405415 0.19619
## d.o.f=2
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.8028  0.3703
## Lag[2*(p+q)+(p+q)-1][5]    5.5692  0.1134
## Lag[4*(p+q)+(p+q)-1][9]    7.1832  0.1843
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.05131 0.500 2.000  0.8208
## ARCH Lag[5]   2.12082 1.440 1.667  0.4451
## ARCH Lag[7]   2.71000 2.315 1.543  0.5699
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  23.3593
## Individual Statistics:              
## mu     0.57382
## ar1    0.04204
## ma1    0.03823
## omega  2.13718
## alpha1 0.36827
## beta1  0.89771
## 
## 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           3.6299 2.861e-04 ***
## Negative Sign Bias  0.8648 3.872e-01    
## Positive Sign Bias  2.6156 8.932e-03 ***
## Joint Effect       42.3943 3.309e-09 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     205.7    2.529e-33
## 2    30     234.7    4.611e-34
## 3    40     281.3    1.884e-38
## 4    50     293.5    1.421e-36
## 
## 
## Elapsed time : 0.5754659

# coef(fit)
plot(fit, which="all")

## 
## please wait...calculating quantiles...

# coeficientes
#coef(fit)
# media incondicional
uncmean(fit)

## [1] 0.0006370267

# varianza incondicional: omega/(alpha1 + beta1)
uncvariance(fit)

## [1] 0.0001388314

# persistencia: alpha1 + beta1
persistence(fit)

## [1] 0.9828551

# half-life: tiempo (días) que toma recuperar la mitad del incremento en la varianza
halflife(fit)

## [1] 40.08113

# plot
dataToPlot <- cbind(residuals(fit), sigma(fit))
colnames(dataToPlot) <- c("Residual", "Volatilidad_condicional")
plot.zoo(dataToPlot, main=paste(symbols,"Retornos Diarios"), col="blue")

\[\begin{align*} Resultados: \\ Volatilidad - largo.plazo: 0.0123809\\ w: 2.380254\times 10^{-6}\\ \alpha: 0.1248875\\ \beta: 0.8579675\\ \\ Volatilidad.incondicional_{diaria}: 0.0117827\\ Volatilidad.incondicional_{anualizada}: 0.1870441\\ \\ \sigma^2 = w + \alpha R_{t-1,t}^2 + \beta \sigma_{t-1,t}^2 \\ \\ Sustituyendo:\\ Escenario_1:\\ R_{t+1}: Retorno +5.00 pct \\ Nueva.volatilidad_{condicional}: 0.0211214\\ \\ Escenario_2:\\ R_{t+1}: Retorno +0.00 pct \\ Nueva.volatilidad_{condicional}: 0.0115713 \end{align*}\]

Simulación Básica GARCH

Procedemos con una simulación a partir del GARCH(1,1) Normal.

garch11.sim = ugarchsim(fit, n.sim=nrow(GSPC.ret), rseed=111)

# plotear retornos reales y simulados
par(mfrow=c(2,1))
  plot(GSPC.ret$GSPC, main=paste(symbols,"Retornos realizados"), grid.col=NA, lwd=1)
  plot(as.xts(garch11.sim@simulation$seriesSim, 
              order.by=index(GSPC.ret)),
       main="GARCH(1,1) Retornos Simulados", grid.col=NA, lwd=1)

par(mfrow=c(1,1))

precio.actual = as.numeric(head(GSPC$GSPC.Adjusted[1]))
precio.futuro.sim1 <- precio.actual*(cumprod(1+garch11.sim@simulation$seriesSim))
GSPC$precio.futuro.sim1[3:nrow(GSPC)] <- precio.futuro.sim1
GSPC <- na.omit(GSPC)

plot.sim <- plot(GSPC$precio.futuro.sim1["2000"], type="l", ylim=c(min(GSPC$precio.futuro.sim1["2000"],GSPC$GSPC.Close["2000"]),
       max(GSPC$precio.futuro.sim1["2000"],GSPC$GSPC.Close["2000"])),main="Futuro alterno simulado - data sintética", grid.col=NA)
plot.sim <- lines(GSPC$GSPC.Close, col="red", on=1)
plot.sim <- addLegend("topleft", legend.names=c("Realizado","Simulado"), col=c("red","black"), lwd=2)
plot.sim

Conclusiones

En la práctica, el modelo GARCH(1,1) con innovaciones Normales no se adapta adecuadamente. Importantes limitaciones:

La especificación de la función generadora modelada, ¿es la mas aproximada a la volatilidad empírica?
Sobre la distribución de los errores, ¿es la mas adecuada?
No permite la asimetría.El mercado no reacciona de la misma manera a las noticias positivas y negativas.
El modelo debería acomodar la posibilidad de saltos.

Conclusión: Es necesario realizar un fit iterativo que pruebe distintos modelos, en distintas configuraciones, incluyendo otras distribuciones de probabilidad. Entonces, elegir el modelo que minimiza el error.