Datos y paquetes

Se usa el índice S&P 500 (precios de cierre diarios, más de 5 años).

# install.packages(c("quantmod","rugarch","moments","tseries"))   # solo la primera vez
library(quantmod); library(rugarch); library(moments); library(tseries)

precios <- Cl(getSymbols("^GSPC", from = "2019-01-01", to = "2026-11-01", auto.assign = FALSE))
# Alternativa con CSV:
# d <- read.csv("datos.csv"); precios <- xts::xts(d$Close, as.Date(d$Date))
cat("Observaciones de precios:", nrow(precios), "| Periodo:", format(start(precios)), "a", format(end(precios)))
## Observaciones de precios: 1969 | Periodo: 2019-01-02 a 2026-10-30

a) Retornos logarítmicos (%), estadísticas descriptivas y Jarque-Bera

r <- na.omit(100 * diff(log(precios)))
rn <- as.numeric(r)
plot(r, main = "Retornos logarítmicos diarios (%)", ylab = "r_t (%)")

desc <- data.frame(Media = mean(rn), Desv.Est = sd(rn), Asimetria = skewness(rn), Curtosis = kurtosis(rn))
round(desc, 4)
jb <- jarque.bera.test(rn); jb
## 
##  Jarque Bera Test
## 
## data:  rn
## X-squared = 20721, df = 2, p-value < 2.2e-16

Interpretación. La curtosis es 18.84 (la normal tiene 3) y la asimetría es -0.66. En la prueba de Jarque-Bera el valor p es <2e-16, por lo que se rechaza la normalidad: los retornos tienen colas pesadas.

b) ACF, Ljung-Box (10 rezagos) y ARCH-LM (5 rezagos)

par(mfrow = c(1, 2))
acf(rn,   lag.max = 30, main = "ACF de r_t")
acf(rn^2, lag.max = 30, main = "ACF de r_t^2")

par(mfrow = c(1, 1))
lb_r  <- Box.test(rn,   lag = 10, type = "Ljung-Box")
lb_r2 <- Box.test(rn^2, lag = 10, type = "Ljung-Box")

arch_lm <- function(x, q = 5) {          # LM = T * R^2 de la regresión auxiliar
  x2 <- x^2; X <- embed(x2, q + 1)
  R2 <- summary(lm(X[, 1] ~ X[, -1]))$r.squared
  LM <- nrow(X) * R2
  c(LM = LM, gl = q, p.valor = 1 - pchisq(LM, q))
}
lm5 <- arch_lm(rn - mean(rn), 5)

data.frame(Prueba = c("Ljung-Box r_t (10)", "Ljung-Box r_t^2 (10)", "ARCH-LM (5)"),
           Estadistico = c(lb_r$statistic, lb_r2$statistic, lm5["LM"]),
           p.valor = c(lb_r$p.value, lb_r2$p.value, lm5["p.valor"]))

Conclusión. Ljung-Box sobre r_t tiene p = <2e-16; sobre r_t^2, p = <2e-16; ARCH-LM, p = <2e-16. Se rechaza H0 en las pruebas sobre los cuadrados: existe autocorrelación en r_t^2 y efectos ARCH, lo que justifica un modelo GARCH.

c) Estimación de GARCH(1,1)-Normal, GARCH(1,1)-t y GJR-GARCH(1,1)-t

sp <- function(modelo, dist) ugarchspec(
  variance.model = list(model = modelo, garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
  distribution.model = dist)

fits <- list(
  "GARCH(1,1)-Normal" = ugarchfit(sp("sGARCH",   "norm"), rn),
  "GARCH(1,1)-t"      = ugarchfit(sp("sGARCH",   "std"),  rn),
  "GJR-GARCH(1,1)-t"  = ugarchfit(sp("gjrGARCH", "std"),  rn)
)
n <- length(rn)

for (nm in names(fits)) { cat("\n====", nm, "====\n"); print(round(fits[[nm]]@fit$robust.matcoef, 4)) }
## 
## ==== GARCH(1,1)-Normal ====
##         Estimate  Std. Error  t value Pr(>|t|)
## mu        0.0838      0.0171   4.8949        0
## omega     0.0497      0.0123   4.0559        0
## alpha1    0.1511      0.0285   5.2950        0
## beta1     0.8079      0.0288  28.0434        0
## 
## ==== GARCH(1,1)-t ====
##         Estimate  Std. Error  t value Pr(>|t|)
## mu        0.1144      0.0145   7.8868        0
## omega     0.0373      0.0085   4.3602        0
## alpha1    0.1519      0.0254   5.9876        0
## beta1     0.8231      0.0227  36.2631        0
## shape     6.1041      0.8196   7.4472        0
## 
## ==== GJR-GARCH(1,1)-t ====
##         Estimate  Std. Error  t value Pr(>|t|)
## mu        0.0833      0.0156   5.3469        0
## omega     0.0348      0.0083   4.2142        0
## alpha1    0.0000      0.0191   0.0000        1
## beta1     0.8473      0.0268  31.6372        0
## gamma1    0.2382      0.0446   5.3433        0
## shape     6.4484      0.9464   6.8132        0
tabla <- do.call(rbind, lapply(names(fits), function(nm) {
  f <- fits[[nm]]; ic <- infocriteria(f)
  data.frame(Modelo = nm, LogLik = likelihood(f), Parametros = length(coef(f)),
             AIC = ic["Akaike", ] * n, BIC = ic["Bayes", ] * n)
}))
rownames(tabla) <- NULL
tabla[, -1] <- round(tabla[, -1], 4); tabla
elegido <- tabla$Modelo[which.min(tabla$BIC)]

(Los errores estándar mostrados son los robustos de cuasi-máxima verosimilitud.) El modelo elegido por menor BIC es GJR-GARCH(1,1)-t.

d) Diagnóstico de los residuos estandarizados del modelo elegido

fit <- fits[[elegido]]
z <- as.numeric(residuals(fit, standardize = TRUE))

lb_z  <- Box.test(z,   lag = 10, type = "Ljung-Box")
lb_z2 <- Box.test(z^2, lag = 10, type = "Ljung-Box")
lm_z  <- arch_lm(z, 5)

data.frame(Prueba = c("Ljung-Box z_t (10)", "Ljung-Box z_t^2 (10)", "ARCH-LM z_t (5)"),
           Estadistico = c(lb_z$statistic, lb_z2$statistic, lm_z["LM"]),
           p.valor = c(lb_z$p.value, lb_z2$p.value, lm_z["p.valor"]))
dist <- fit@model$modeldesc$distribution
if (dist == "std") {
  nu <- coef(fit)["shape"]
  qqplot(qdist("std", ppoints(length(z)), shape = nu), z, main = "QQ-plot de z_t vs t de Student",
         xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales")
} else {
  qqnorm(z, main = "QQ-plot de z_t vs Normal")
}
abline(0, 1, col = "red")

Interpretación. Los p-valores sobre z_t^2 (0.759) y ARCH-LM (0.551) no rechazan H0: el modelo captura adecuadamente la heterocedasticidad condicional. En el QQ-plot, los puntos cercanos a la recta indican buen ajuste de la distribución; desviaciones en los extremos señalan colas mal capturadas.

e) Volatilidad condicional, pronóstico a 10 días y VaR al 99 %

plot(sigma(fit), main = paste("Volatilidad condicional estimada -", elegido), ylab = "sigma_t (%)")

fc <- ugarchforecast(fit, n.ahead = 10)
pron <- data.frame(Dia = 1:10, Volatilidad_pct = round(as.numeric(sigma(fc)), 4))
pron
plot(pron$Dia, pron$Volatilidad_pct, type = "b", xlab = "Horizonte (días)", ylab = "sigma (%)",
     main = "Pronóstico de volatilidad a 10 días")

mu     <- coef(fit)["mu"]
sigma1 <- as.numeric(sigma(fc))[1]
q01    <- if (dist == "std") qdist("std", 0.01, shape = coef(fit)["shape"]) else qnorm(0.01)
VaR99  <- -(mu + sigma1 * q01)
c(mu = unname(mu), sigma_t1 = sigma1, q_0.01 = unname(q01), VaR_99 = unname(VaR99))
##          mu    sigma_t1      q_0.01      VaR_99 
##  0.08334164  1.20260241 -2.55054599  2.98395110

Resultado. La volatilidad pronosticada para el día siguiente es 1.2026 % y el VaR a un día al 99 % es 2.984 %: con 99 % de confianza, la pérdida diaria no debería superar ese porcentaje del valor del portafolio. El pronóstico a 10 días muestra la reversión de la volatilidad hacia su nivel de largo plazo.