Proyecto Final Riesgos - Box Jenkins
Juan Vicente Olave
r Sys.Date()
## Warning: package 'quantmod' was built under R version 4.5.2## Cargando paquete requerido: xts## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Cargando paquete requerido: TTR## Warning: package 'TTR' was built under R version 4.5.2## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Warning: package 'tidyverse' was built under R version 4.5.2## Warning: package 'ggplot2' was built under R version 4.5.2## Warning: package 'tidyr' was built under R version 4.5.2## Warning: package 'readr' was built under R version 4.5.2## Warning: package 'purrr' was built under R version 4.5.2## Warning: package 'forcats' was built under R version 4.5.2## Warning: package 'lubridate' was built under R version 4.5.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks xts::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks xts::last()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Warning: package 'plotly' was built under R version 4.5.2##
## Adjuntando el paquete: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout## Warning: package 'PerformanceAnalytics' was built under R version 4.5.2##
## Adjuntando el paquete: 'PerformanceAnalytics'
##
## The following object is masked from 'package:graphics':
##
## legendINTRODUCCION
NVIDIA
1. Serie de precios de NVIDIA
# 1) Descargar precios ajustados NVDA
ticker_nvda <- "NVDA"
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
# Descargar desde Yahoo y extraer Precio Ajustado
getSymbols(ticker_nvda, src = "yahoo", from = start_date, to = end_date, auto.assign = TRUE)## [1] "NVDA"prices_nvda <- Ad(get(ticker_nvda))
colnames(prices_nvda) <- "NVDA"
prices_nvda <- na.omit(prices_nvda)
# Resumen
cat("NVDA - observaciones:", nrow(prices_nvda), "\n")## NVDA - observaciones: 2009## NVDA
## 2017-11-06 5.179376
## 2017-11-07 5.238674
## 2017-11-08 5.167763
## 2017-11-09 5.072887
## 2017-11-10 5.340220
## 2017-11-13 5.253497## NVDA
## 2025-10-27 191.49
## 2025-10-28 201.03
## 2025-10-29 207.04
## 2025-10-30 202.89
## 2025-10-31 202.49
## 2025-11-03 206.882. Series de rendimientos financieros y rendimientos financieros al cuadrado
# 2) Cálculo de rendimientos logarÃtmicos diarios y rendimientos al cuadrado
library(quantmod)
r_nvda <- dailyReturn(prices_nvda, type = "log")
colnames(r_nvda) <- "r_NVDA"
r_nvda_sq <- r_nvda^2
colnames(r_nvda_sq) <- "r_NVDA_sq"
nvda_data <- merge(prices_nvda, r_nvda, r_nvda_sq, join = "inner")
# Resumen
summary(nvda_data)## Index NVDA r_NVDA r_NVDA_sq
## Min. :2017-11-06 Min. : 3.151 Min. :-0.207711 Min. :0.0000000
## 1st Qu.:2019-11-05 1st Qu.: 6.231 1st Qu.:-0.015097 1st Qu.:0.0000606
## Median :2021-11-02 Median : 16.508 Median : 0.002810 Median :0.0003078
## Mean :2021-11-03 Mean : 40.905 Mean : 0.001835 Mean :0.0010440
## 3rd Qu.:2023-11-01 3rd Qu.: 46.804 3rd Qu.: 0.019562 3rd Qu.:0.0010562
## Max. :2025-11-03 Max. :207.040 Max. : 0.218088 Max. :0.04756233. Graficas de las 3 series: precios, rendimientos financieros y rendimientos financieros al cuadrado
library(ggplot2)
library(zoo)
library(tibble)
# Preparar data.frame para ggplot
df_nvda <- fortify.zoo(nvda_data, name = "date") %>% as_tibble()
# Colores
col_price <- "#2c7fb8"
col_ret <- "#238b45"
col_ret2 <- "#ae017e"
# Gráfico precio
p_price <- ggplot(df_nvda, aes(x = date, y = NVDA)) +
geom_line(color = col_price) +
labs(title = "NVDA — Precio ajustado", x = NULL, y = "USD") +
theme_minimal() + theme(plot.title = element_text(face = "bold"))
# Gráfico rendimientos
p_ret <- ggplot(df_nvda, aes(x = date, y = r_NVDA)) +
geom_line(color = col_ret) +
labs(title = "NVDA — Rendimientos logarÃtmicos diarios", x = NULL, y = "r_t") +
theme_minimal() + theme(plot.title = element_text(face = "bold"))
# Gráfico rendimientos^2
p_ret_sq <- ggplot(df_nvda, aes(x = date, y = r_NVDA_sq)) +
geom_line(color = col_ret2) +
labs(title = "NVDA — Rendimientos al cuadrado", x = NULL, y = "r_t^2") +
theme_minimal() + theme(plot.title = element_text(face = "bold"))
# Mostrar los gráficos (aparecen uno debajo del otro)
p_price
4. Estimacion del correlograma para cada una de las series del numeral 3
# 4) Correlogramas: ACF y PACF para precio, r_t y r_t^2
graphics::par(mfrow = c(3,2), mar = c(4,4,2,1))
# Precio (niveles)
acf(coredata(nvda_data$NVDA), main = "ACF - NVDA Precio")
pacf(coredata(nvda_data$NVDA), main = "PACF - NVDA Precio")
# Rendimientos
acf(coredata(nvda_data$r_NVDA), main = "ACF - NVDA Rendimientos")
pacf(coredata(nvda_data$r_NVDA), main = "PACF - NVDA Rendimientos")
# Rendimientos^2
acf(coredata(nvda_data$r_NVDA_sq), main = "ACF - NVDA Rendimientos^2")
pacf(coredata(nvda_data$r_NVDA_sq), main = "PACF - NVDA Rendimientos^2")
5.Revisar la normalidad en cada una de las series del numeral 3
library(tseries) # jarque.bera.test
library(moments) # skewness, kurtosis
library(ggplot2)
library(gridExtra)
library(ggpubr)
# Extraer series como vectores numéricos
x_price <- coredata(nvda_data$NVDA)
x_ret <- coredata(nvda_data$r_NVDA)
x_ret2 <- coredata(nvda_data$r_NVDA_sq)
# Función auxiliar para gráficos y tests
check_normality <- function(x, name, trim_outliers = FALSE) {
# opcional: quitar NA (ya deberÃan estar omitidos) y posibles outliers si se desea
x <- na.omit(as.numeric(x))
if (trim_outliers) {
q <- quantile(x, probs = c(0.01, 0.99))
x <- x[x >= q[1] & x <= q[2]]
}
# EstadÃsticas descriptivas
n <- length(x)
m <- mean(x)
s <- sd(x)
sk <- skewness(x)
kt <- kurtosis(x) # Fisher: kurtosis normal = 3 -> aquà devuelve 3 para normal?
# moments::kurtosis por defecto devuelve la curtosis de Fisher (normal = 3)
cat("===", name, "===\n")
cat("N =", n, "\n")
cat("Media =", signif(m,6), " SD =", signif(s,6), "\n")
cat("Skewness =", signif(sk,6), " Kurtosis =", signif(kt,6), "\n\n")
# Pruebas de normalidad
# Shapiro-Wilk: recomendado para n <= 5000 (si n muy grande, siempre rechaza H0)
if (n <= 5000) {
sh <- shapiro.test(x)
cat("Shapiro-Wilk: W =", signif(sh$statistic,6), " p =", signif(sh$p.value,6), "\n")
} else {
cat("Shapiro-Wilk: n >", 5000, " (no recomendable). Se omite.\n")
}
# Jarque-Bera (adecuado para series financieras)
jb <- tryCatch(jarque.bera.test(x), error = function(e) NULL)
if (!is.null(jb)) {
cat("Jarque-Bera: X2 =", signif(jb$statistic,6), " p =", signif(jb$p.value,6), "\n")
}
cat("\n")
# Gráficos: histograma + densidad normal y QQ-plot
df <- data.frame(x = x)
p_hist <- ggplot(df, aes(x = x)) +
geom_histogram(aes(y = ..density..), bins = 40, fill = "#2c7fb8", alpha = 0.6, color = "white") +
geom_density(color = "#238b45", size = 1) +
stat_function(fun = dnorm, args = list(mean = m, sd = s), color = "red", linetype = "dashed", size = 1) +
labs(title = paste(name, "- Histograma y densidad (normal en rojo)"), x = NULL, y = "Densidad") +
theme_minimal()
p_qq <- ggqqplot(df$x, title = paste(name, "- QQ-plot"), ggtheme = theme_minimal()) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red")
list(summary = list(n = n, mean = m, sd = s, skew = sk, kurt = kt),
shapiro = if (n <= 5000) sh else NULL,
jarque = jb,
plots = list(hist = p_hist, qq = p_qq))
}
# Aplicar a cada serie
res_price <- check_normality(x_price, "NVDA - Precio (niveles)", trim_outliers = FALSE)## === NVDA - Precio (niveles) ===
## N = 2009
## Media = 40.9045 SD = 50.5336
## Skewness = 1.48733 Kurtosis = 3.9
##
## Shapiro-Wilk: W = 0.723406 p = 1.55688e-49
## Jarque-Bera: X2 = 808.504 p = 0## === NVDA - Rendimientos log (r_t) ===
## N = 2009
## Media = 0.00183547 SD = 0.032267
## Skewness = -0.197361 Kurtosis = 7.74228
##
## Shapiro-Wilk: W = 0.955945 p = 3.55495e-24
## Jarque-Bera: X2 = 1895.58 p = 0## === NVDA - Rendimientos^2 (r_t^2) ===
## N = 2009
## Media = 0.00104401 SD = 0.00269638
## Skewness = 9.5069 Kurtosis = 128.744
##
## Shapiro-Wilk: W = 0.342519 p = 9.28302e-65
## Jarque-Bera: X2 = 1353830 p = 0# Mostrar gráficos en panel (2 columnas)
grid.arrange(res_price$plots$hist, res_price$plots$qq,
res_ret$plots$hist, res_ret$plots$qq,
res_ret2$plots$hist, res_ret2$plots$qq,
ncol = 2)
6. Estimar la estacionareidad de cada una de las 3 series del numeral 3
library(tseries) # adf.test, pp.test
library(urca) # ur.kpss
library(dplyr)
library(ggplot2)
library(gridExtra)
library(forecast)
# Extraer series (ya definidas anteriormente)
x_price <- coredata(nvda_data$NVDA)
x_ret <- coredata(nvda_data$r_NVDA)
x_ret2 <- coredata(nvda_data$r_NVDA_sq)
# Función auxiliar para pruebas de estacionariedad
test_stationarity <- function(x, name, maxlag = NULL) {
x <- na.omit(as.numeric(x))
n <- length(x)
# ADF: adf.test from tseries (null: unit root -> non-stationary)
adf <- tryCatch(adf.test(x, alternative = "stationary", k = maxlag), error = function(e) NULL)
# PP: Phillips-Perron (null: unit root -> non-stationary)
pp <- tryCatch(pp.test(x, alternative = "stationary"), error = function(e) NULL)
# KPSS from urca package (null: stationary)
# usamos tipo "tau" (trend) y "mu" (level) para comparar. Aqui se usa level ("mu")
kpss_mu <- tryCatch(ur.kpss(x, type = "mu", lags = "short"), error = function(e) NULL)
kpss_tau <- tryCatch(ur.kpss(x, type = "tau", lags = "short"), error = function(e) NULL)
cat("====", name, "====\n")
if (!is.null(adf)) {
cat("ADF: stat =", signif(adf$statistic,6), " p =", signif(adf$p.value,6), "(H0: unit root -> non-stationary)\n")
} else cat("ADF: error\n")
if (!is.null(pp)) {
cat("PP: stat =", signif(pp$statistic,6), " p =", signif(pp$p.value,6), "(H0: unit root -> non-stationary)\n")
} else cat("PP: error\n")
if (!is.null(kpss_mu)) {
cat("KPSS (level) LM stat =", signif(kpss_mu@teststat,6),
" critvals:", paste(names(kpss_mu@cval), signif(kpss_mu@cval,4), collapse = ", "),
"(H0: stationary)\n")
} else cat("KPSS (level): error\n")
if (!is.null(kpss_tau)) {
cat("KPSS (trend) LM stat =", signif(kpss_tau@teststat,6),
" critvals:", paste(names(kpss_tau@cval), signif(kpss_tau@cval,4), collapse = ", "), "\n")
} else cat("KPSS (trend): error\n")
cat("\n")
# Gráficos de la serie y de la ACF (útiles para ver comportamiento)
df <- data.frame(date = index(nvda_data), y = x)
p_ts <- ggplot(df, aes(x = date, y = y)) + geom_line(color = "#2c7fb8") +
labs(title = paste(name, "- Serie temporal"), x = NULL, y = NULL) + theme_minimal()
p_acf <- ggAcf(x, lag.max = 40) + ggtitle(paste(name, "- ACF (40 lags)")) + theme_minimal()
list(adf = adf, pp = pp, kpss_mu = kpss_mu, kpss_tau = kpss_tau,
plots = list(ts = p_ts, acf = p_acf))
}
# Aplicar pruebas a las 3 series
res_price <- test_stationarity(x_price, "NVDA - Precio (niveles)")## ==== NVDA - Precio (niveles) ====
## ADF: error
## PP: stat = 0.855745 p = 0.99 (H0: unit root -> non-stationary)
## KPSS (level) LM stat = 15.9319 critvals: 0.347, 0.463, 0.574, 0.739 (H0: stationary)
## KPSS (trend) LM stat = 4.34961 critvals: 0.119, 0.146, 0.176, 0.216## ==== NVDA - Rendimientos log (r_t) ====
## ADF: error
## PP: stat = -2185.32 p = 0.01 (H0: unit root -> non-stationary)
## KPSS (level) LM stat = 0.218962 critvals: 0.347, 0.463, 0.574, 0.739 (H0: stationary)
## KPSS (trend) LM stat = 0.0564293 critvals: 0.119, 0.146, 0.176, 0.216## ==== NVDA - Rendimientos^2 (r_t^2) ====
## ADF: error
## PP: stat = -2036.78 p = 0.01 (H0: unit root -> non-stationary)
## KPSS (level) LM stat = 0.136937 critvals: 0.347, 0.463, 0.574, 0.739 (H0: stationary)
## KPSS (trend) LM stat = 0.095654 critvals: 0.119, 0.146, 0.176, 0.216# Mostrar gráficos
grid.arrange(res_price$plots$ts, res_price$plots$acf,
res_ret$plots$ts, res_ret$plots$acf,
res_ret2$plots$ts, res_ret2$plots$acf,
ncol = 2)
# Resumen tabular de resultados (p/ informe)
summ_row <- function(name, res) {
adf_stat <- if (!is.null(res$adf)) res$adf$statistic else NA
adf_p <- if (!is.null(res$adf)) res$adf$p.value else NA
pp_stat <- if (!is.null(res$pp)) res$pp$statistic else NA
pp_p <- if (!is.null(res$pp)) res$pp$p.value else NA
kpss_mu_stat <- if (!is.null(res$kpss_mu)) as.numeric(res$kpss_mu@teststat) else NA
kpss_tau_stat<- if (!is.null(res$kpss_tau)) as.numeric(res$kpss_tau@teststat) else NA
data.frame(serie = name,
adf_stat = adf_stat, adf_p = adf_p,
pp_stat = pp_stat, pp_p = pp_p,
kpss_mu = kpss_mu_stat, kpss_tau = kpss_tau_stat)
}
tab_results <- bind_rows(
summ_row("Precio (niveles)", res_price),
summ_row("Rendimientos (r_t)", res_ret),
summ_row("Rendimientos^2 (r_t^2)", res_ret2)
)
print(tab_results)## serie adf_stat adf_p pp_stat
## Dickey-Fuller Z(alpha)...1 Precio (niveles) NA NA 0.8557445
## Dickey-Fuller Z(alpha)...2 Rendimientos (r_t) NA NA -2185.3174035
## Dickey-Fuller Z(alpha)...3 Rendimientos^2 (r_t^2) NA NA -2036.7783960
## pp_p kpss_mu kpss_tau
## Dickey-Fuller Z(alpha)...1 0.99 15.9319395 4.34960937
## Dickey-Fuller Z(alpha)...2 0.01 0.2189617 0.05642927
## Dickey-Fuller Z(alpha)...3 0.01 0.1369369 0.095654027. Seleccion de series que cumplen con la propiedad de estacionareidad
Elijo trabajar con los rendimientos continuamente compuestos (logarÃtmicos) y no con los rendimientos al cuadrado porque los log‑rendimientos son la variable relevante para analizar la rentabilidad: convierten la serie de precios no estacionaria en una serie estacionaria adecuada para modelización, eliminan la dependencia de la escala y permiten la composición aditiva de retornos (la suma de log‑retornos iguala el log‑cambio total), mientras que los rendimientos al cuadrado sólo reflejan la varianza o volatilidad instantánea y no la dirección ni el nivel esperado del retorno
8. Con la serie seleccionada () estimacion de los modelos AR(p), MA(q) y ARMA(p,q)
pkgs <- c("quantmod", "forecast", "tseries", "FinTS", "rugarch", "xts", "zoo", "dplyr", "knitr")
install_if_missing <- function(p) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
invisible(lapply(pkgs, install_if_missing))
library(quantmod); library(forecast); library(tseries); library(FinTS); library(rugarch)
library(xts); library(zoo); library(dplyr); library(knitr)
ticker <- "NVDA"
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
maxP <- 5
maxQ <- 5
include.mean <- TRUE
lag_test <- 20
investment_ref <- 100000
lambda_rm <- 0.94
niveles_conf <- c(0.90, 0.95, 0.99)
res_folder <- "results_arma_garch_var"
if (!dir.exists(res_folder)) dir.create(res_folder, recursive = TRUE)
getSymbols(ticker, src = "yahoo", from = start_date, to = end_date, auto.assign = TRUE, warnings = FALSE)## [1] "NVDA"prices <- tryCatch(Cl(get(ticker)), error = function(e) stop("Error descargando ", ticker))
write.zoo(prices, file.path(res_folder, paste0(ticker, "_prices.csv")), sep = ",")
rendimientos_xts <- na.omit(diff(log(prices)))
colnames(rendimientos_xts) <- paste0("R_", ticker)
serie_xts <- rendimientos_xts[,1]
serie <- as.numeric(na.omit(coredata(serie_xts)))
n <- length(serie)
if (n < 100) warning("Pocos datos (<100). Considera ampliar rango de fechas.")
cat("Observaciones de rendimientos:", n, "\n")## Observaciones de rendimientos: 2008fmt <- function(x, digits = 6) ifelse(is.na(x), NA, formatC(as.numeric(x), format = "f", digits = digits))
plot_acf_eviews <- function(x, main = "", lag.max = 20) {
x_core <- as.numeric(na.omit(x))
acf_obj <- acf(x_core, lag.max = lag.max, plot = FALSE)
acf_vals <- acf_obj$acf[-1]
lags <- seq_along(acf_vals)
ci <- 1.96 / sqrt(length(x_core))
ylim_low <- min(acf_vals, -ci) - 0.02 * diff(range(c(acf_vals, -ci, ci)))
ylim_high <- max(acf_vals, ci) + 0.02 * diff(range(c(acf_vals, -ci, ci)))
plot(NA, xlim = c(0.5, length(acf_vals) + 0.5), ylim = c(ylim_low, ylim_high),
xlab = "Rezagos", ylab = "ACF", main = main, xaxt = "n", bty = "n")
axis(1, at = lags)
abline(h = ci, lty = 2, col = "blue"); abline(h = -ci, lty = 2, col = "blue"); abline(h = 0)
ancho <- 0.4
for (i in lags) {
top <- acf_vals[i]
if (top >= 0) rect(i - ancho, 0, i + ancho, top, col = "grey60", border = "grey40")
else rect(i - ancho, top, i + ancho, 0, col = "grey60", border = "grey40")
}
}
fit_safe <- function(p, q, x) {
tryCatch({
fit <- arima(x, order = c(p, 0, q), include.mean = include.mean, method = "ML")
list(fit = fit, AIC = AIC(fit))
}, error = function(e) list(fit = NULL, AIC = Inf))
}
best_ar <- list(AIC = Inf)
for (p in 1:maxP) {
r <- fit_safe(p, 0, serie)
if (r$AIC < best_ar$AIC) best_ar <- c(list(p = p), r)
}
best_ar_p <- best_ar$p; model_ar <- best_ar$fit
best_ma <- list(AIC = Inf)
for (q in 1:maxQ) {
r <- fit_safe(0, q, serie)
if (r$AIC < best_ma$AIC) best_ma <- c(list(q = q), r)
}
best_ma_q <- best_ma$q; model_ma <- best_ma$fit
best_arma <- list(AIC = Inf)
for (p in 1:maxP) for (q in 1:maxQ) {
r <- fit_safe(p, q, serie)
if (r$AIC < best_arma$AIC) best_arma <- c(list(p = p, q = q), r)
}
best_arma_p <- best_arma$p; best_arma_q <- best_arma$q; model_arma <- best_arma$fit
aic_family <- data.frame(
Modelo = c(paste0("AR(",best_ar_p,")"), paste0("MA(",best_ma_q,")"), paste0("ARMA(",best_arma_p,",",best_arma_q,")")),
AIC = c(AIC(model_ar), AIC(model_ma), AIC(model_arma)),
stringsAsFactors = FALSE
)
write.csv(aic_family, file.path(res_folder, paste0(ticker, "_aic_family.csv")), row.names = FALSE)
kable(aic_family, caption = "AIC - Mejor AR/MA/ARMA por familia")| Modelo | AIC |
|---|---|
| AR(2) | -8099.532 |
| MA(3) | -8099.780 |
| ARMA(4,3) | -8105.891 |
9. En cada uno de los modelos estimados en el numeral 8, revisar la presencia de posibles problemas de autocorrelación y heterocedasticidad
run_model_diagnostics <- function(model, model_name, lag_test = 20) {
res <- residuals(model)
lb <- Box.test(res, lag = lag_test, type = "Ljung-Box")
arch <- tryCatch(FinTS::ArchTest(res, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
par(mfrow = c(2,1))
plot_acf_eviews(res, main = paste("ACF residuales -", model_name), lag.max = lag_test)
plot_acf_eviews(res^2, main = paste("ACF residuales^2 -", model_name), lag.max = lag_test)
par(mfrow = c(1,1))
data.frame(
Prueba = c("Ljung-Box (res.)", "ARCH-LM (res.)"),
Estadistico = c(round(as.numeric(lb$statistic),4), round(as.numeric(arch$statistic),4)),
Valor_p = c(lb$p.value, arch$p.value),
Decision = c(ifelse(lb$p.value < 0.05, "Rechazar H0 (autocorrel.)", "No rechazar H0"),
ifelse(arch$p.value < 0.05, "Rechazar H0 (heteroced.)", "No rechazar H0")),
stringsAsFactors = FALSE
)
}
diag_ar <- run_model_diagnostics(model_ar, paste0("AR(",best_ar_p,")"), lag_test); diag_ar
diag_arma <- run_model_diagnostics(model_arma, paste0("ARMA(",best_arma_p,",",best_arma_q,")"), lag_test); diag_arma
10. Si encuentre problemas mencionados en el numeral 9, proceder a estimar los mismos modelos pero bajo metodologÃa de estimación GARCH
## --- Configuración de selección ---
prefer_garch <- TRUE # si TRUE, preferir modelos GARCH cuando convergen y mejoran diagnósticos
garch_order <- c(1,1) # orden GARCH por defecto
dist_garch <- "std" # distribución en GARCH ("norm","std",...)
## --- 10bis: Estimar modelos GARCH para AR/MA/ARMA si es necesario o solicitado ---
need_garch <- any(c(diag_ar$Valor_p[2], diag_ma$Valor_p[2], diag_arma$Valor_p[2]) < 0.05, na.rm = TRUE) || prefer_garch
fit_ar_garch <- fit_ma_garch <- fit_arma_garch <- NULL
conv_ar_garch <- conv_ma_garch <- conv_arma_garch <- NA
if (need_garch) {
# Especificaciones
spec_ar <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_ar_p, 0), include.mean = include.mean),
distribution.model = dist_garch)
spec_ma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(0, best_ma_q), include.mean = include.mean),
distribution.model = dist_garch)
spec_arma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_arma_p, best_arma_q), include.mean = include.mean),
distribution.model = dist_garch)
# Ajustes con tryCatch para no interrumpir flujo
fit_ar_garch <- tryCatch(ugarchfit(spec_ar, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("AR-GARCH error: ", e$message); NULL })
fit_ma_garch <- tryCatch(ugarchfit(spec_ma, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("MA-GARCH error: ", e$message); NULL })
fit_arma_garch <- tryCatch(ugarchfit(spec_arma, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("ARMA-GARCH error: ", e$message); NULL })
# Convergencia (0 = OK)
conv_ar_garch <- if (!is.null(fit_ar_garch)) tryCatch(fit_ar_garch@fit$convergence, error = function(e) NA) else NA
conv_ma_garch <- if (!is.null(fit_ma_garch)) tryCatch(fit_ma_garch@fit$convergence, error = function(e) NA) else NA
conv_arma_garch <- if (!is.null(fit_arma_garch)) tryCatch(fit_arma_garch@fit$convergence, error = function(e) NA) else NA
}
conv_df <- data.frame(
Modelo = c("AR-GARCH", "MA-GARCH", "ARMA-GARCH"),
Converge = c(conv_ar_garch, conv_ma_garch, conv_arma_garch),
stringsAsFactors = FALSE
)
write.csv(conv_df, file.path(res_folder, paste0(ticker, "_garch_convergence.csv")), row.names = FALSE)
print(conv_df) # o kable(conv_df) en Rmd## Modelo Converge
## 1 AR-GARCH 0
## 2 MA-GARCH 0
## 3 ARMA-GARCH 011. Seleccionar el modelo de mejor ajuste. Usar el criterio AIC
## --- 11: Comparación ampliada por AIC + diagnósticos para favorecer GARCH si corresponde ---
# Construir tabla de comparación base con AIC de modelos lineales (si existen)
comparacion <- data.frame(Model = character(), AIC = numeric(), Converge = logical(),
ARCH_p = numeric(), LB_p_resid = numeric(), LB_p_resid2 = numeric(),
IsGARCH = logical(), stringsAsFactors = FALSE)
# helper para tests diagnósticos sobre un objeto (puede ser arima o uGARCHfit)
diagnosticar_modelo <- function(obj) {
out <- list(arch_p = NA, lb_p = NA, lb_p2 = NA, converge = NA)
if (is.null(obj)) return(out)
if (inherits(obj, "uGARCHfit")) {
out$converge <- tryCatch(obj@fit$convergence == 0, error = function(e) NA)
res_std <- tryCatch(residuals(obj, standardize = TRUE), error = function(e) residuals(obj))
res_plain <- tryCatch(residuals(obj), error = function(e) res_std)
resq <- res_std^2
out$lb_p <- tryCatch(Box.test(as.numeric(res_std), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$lb_p2 <- tryCatch(Box.test(as.numeric(resq), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_std), lags = lag_test)$p.value, error = function(e) NA)
} else {
# arima / lm style objects -> residuals()
out$converge <- TRUE
res_plain <- tryCatch(residuals(obj), error = function(e) NULL)
if (is.null(res_plain)) return(out)
resq <- res_plain^2
out$lb_p <- tryCatch(Box.test(as.numeric(res_plain), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$lb_p2 <- tryCatch(Box.test(as.numeric(resq), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_plain), lags = lag_test)$p.value, error = function(e) NA)
}
return(out)
}
# Añadir modelos lineales si existen (model_ar, model_ma, model_arma)
if (exists("model_ar") && !is.null(model_ar)) {
di <- diagnosticar_modelo(model_ar)
comparacion <- rbind(comparacion, data.frame(Model = paste0("AR(",best_ar_p,")"), AIC = AIC(model_ar),
Converge = di$converge, ARCH_p = di$arch_p,
LB_p_resid = di$lb_p, LB_p_resid2 = di$lb_p2, IsGARCH = FALSE))
}
if (exists("model_ma") && !is.null(model_ma)) {
di <- diagnosticar_modelo(model_ma)
comparacion <- rbind(comparacion, data.frame(Model = paste0("MA(",best_ma_q,")"), AIC = AIC(model_ma),
Converge = di$converge, ARCH_p = di$arch_p,
LB_p_resid = di$lb_p, LB_p_resid2 = di$lb_p2, IsGARCH = FALSE))
}
if (exists("model_arma") && !is.null(model_arma)) {
di <- diagnosticar_modelo(model_arma)
comparacion <- rbind(comparacion, data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")"), AIC = AIC(model_arma),
Converge = di$converge, ARCH_p = di$arch_p,
LB_p_resid = di$lb_p, LB_p_resid2 = di$lb_p2, IsGARCH = FALSE))
}
# Añadir modelos GARCH estimados (si existen)
add_garch_entry <- function(fit_obj, label) {
if (is.null(fit_obj)) return(NULL)
aicv <- tryCatch(infocriteria(fit_obj)["Akaike", 1], error = function(e) NA)
di <- diagnosticar_modelo(fit_obj)
return(data.frame(Model = label, AIC = aicv, Converge = di$converge, ARCH_p = di$arch_p,
LB_p_resid = di$lb_p, LB_p_resid2 = di$lb_p2, IsGARCH = TRUE, stringsAsFactors = FALSE))
}
if (!is.null(fit_ar_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_ar_garch, paste0("AR(",best_ar_p,")-GARCH(1,1)")))
if (!is.null(fit_ma_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_ma_garch, paste0("MA(",best_ma_q,")-GARCH(1,1)")))
if (!is.null(fit_arma_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_arma_garch, paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)")))
# Ordenar por AIC
comparacion <- comparacion %>% arrange(AIC)
write.csv(comparacion, file.path(res_folder, paste0(ticker, "_comparacion_aic_diag.csv")), row.names = FALSE)
print(comparacion) # o kable(comparacion)## Model AIC Converge ARCH_p LB_p_resid
## Chi-squared2 ARMA(4,3) -8105.890690 TRUE 2.990871e-11 0.9550833
## Chi-squared1 MA(3) -8099.780487 TRUE 4.794662e-13 0.2650999
## Chi-squared AR(2) -8099.532097 TRUE 7.023041e-13 0.1579667
## Chi-squared4 MA(3)-GARCH(1,1) -4.238223 TRUE 9.979106e-01 0.8491478
## Chi-squared3 AR(2)-GARCH(1,1) -4.238166 TRUE 9.980078e-01 0.7581683
## Chi-squared5 ARMA(4,3)-GARCH(1,1) -4.236978 TRUE 9.973380e-01 0.8700331
## LB_p_resid2 IsGARCH
## Chi-squared2 0.0000000 FALSE
## Chi-squared1 0.0000000 FALSE
## Chi-squared 0.0000000 FALSE
## Chi-squared4 0.9982469 TRUE
## Chi-squared3 0.9983352 TRUE
## Chi-squared5 0.9977780 TRUE# Selección final: si prefer_garch TRUE, buscar el primer modelo GARCH que converge y mejora diagnóstico
mejor_modelo <- NA
if (prefer_garch) {
# obtener modelos GARCH válidos (converge == TRUE o 0)
garch_candidates <- comparacion %>% filter(IsGARCH == TRUE & (Converge == TRUE | Converge == 0))
if (nrow(garch_candidates) > 0) {
# preferir aquél con menor AIC entre candidatos
mejor_modelo <- garch_candidates$Model[which.min(garch_candidates$AIC)]
mensaje_pref <- paste0("Se prefirió modelo GARCH por prefer_garch: ", mejor_modelo)
message(mensaje_pref)
} else {
# si no hay GARCH válidos, seleccionar por AIC global
mejor_modelo <- comparacion$Model[1]
message("No hay GARCH válidos convergidos; se selecciona mejor por AIC global: ", mejor_modelo)
}
} else {
# no preferir GARCH: elegir mejor por AIC
mejor_modelo <- comparacion$Model[1]
message("Selección por AIC: ", mejor_modelo)
}12.Al mejor modelo, aplicar los test de heterocedasticidad y a través del correlograma de los residuales revisar si persiste el problema de autocorrelación.
objs <- list(model_ar, model_ma, model_arma, fit_ar_garch, fit_ma_garch, fit_arma_garch)
names(objs) <- c(paste0("AR(", best_ar_p, ")"),
paste0("MA(", best_ma_q, ")"),
paste0("ARMA(", best_arma_p, ",", best_arma_q, ")"),
paste0("AR(", best_ar_p, ")-GARCH(1,1)"),
paste0("MA(", best_ma_q, ")-GARCH(1,1)"),
paste0("ARMA(", best_arma_p, ",", best_arma_q, ")-GARCH(1,1)"))
modelo_final <- objs[[mejor_modelo]]
if (is.null(modelo_final)) stop("No se encontró el objeto del modelo seleccionado: ", mejor_modelo)
if (inherits(modelo_final, "uGARCHfit")) {
residuales <- residuals(modelo_final, standardize = TRUE)
residuales_cuad <- residuales^2
lb_res <- Box.test(residuales, lag = lag_test, type = "Ljung-Box")
lb_res2 <- Box.test(residuales_cuad, lag = lag_test, type = "Ljung-Box")
arch_res <- tryCatch(FinTS::ArchTest(residuales, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
} else {
residuales <- residuals(modelo_final)
residuales_cuad <- residuales^2
lb_res <- Box.test(residuales, lag = lag_test, type = "Ljung-Box")
lb_res2 <- Box.test(residuales_cuad, lag = lag_test, type = "Ljung-Box")
arch_res <- tryCatch(FinTS::ArchTest(residuales, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
}
diag_final <- data.frame(Test = c("LjungBox resid.", "LjungBox resid.^2", "ARCH-LM"),
Stat = c(lb_res$statistic, lb_res2$statistic, arch_res$statistic),
P_value = c(lb_res$p.value, lb_res2$p.value, arch_res$p.value))
write.csv(diag_final, file.path(res_folder, paste0(ticker, "_diag_final_modelo.csv")), row.names = FALSE)
kable(diag_final, caption = paste("Diagnóstico final -", mejor_modelo))| Test | Stat | P_value |
|---|---|---|
| LjungBox resid. | 13.62169 | 0.8491478 |
| LjungBox resid.^2 | 6.39554 | 0.9982469 |
| ARCH-LM | 6.55452 | 0.9979106 |
par(mfrow = c(2,1))
plot_acf_eviews(residuales, main = paste("ACF residuales -", mejor_modelo), lag.max = lag_test)
plot_acf_eviews(residuales_cuad, main = paste("ACF residuales^2 -", mejor_modelo), lag.max = lag_test)
13.Con la ecuación del modelo seleccionado, calcular la volatilidad y posteriormente el VaR
# Mapear el nombre a objeto real
objs <- list()
if (exists("model_ar")) objs[[paste0("AR(",best_ar_p,")")]] <- model_ar
if (exists("model_ma")) objs[[paste0("MA(",best_ma_q,")")]] <- model_ma
if (exists("model_arma")) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")")]] <- model_arma
if (!is.null(fit_ar_garch)) objs[[paste0("AR(",best_ar_p,")-GARCH(1,1)")]] <- fit_ar_garch
if (!is.null(fit_ma_garch)) objs[[paste0("MA(",best_ma_q,")-GARCH(1,1)")]] <- fit_ma_garch
if (!is.null(fit_arma_garch)) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)")]] <- fit_arma_garch
if (!mejor_modelo %in% names(objs)) stop("No se encontró el objeto asociado al modelo seleccionado: ", mejor_modelo)
modelo_final <- objs[[mejor_modelo]]
# Obtener residuos (estandarizados si es GARCH)
if (inherits(modelo_final, "uGARCHfit")) {
residuales <- tryCatch(residuals(modelo_final, standardize = TRUE), error = function(e) residuals(modelo_final))
} else {
residuales <- residuals(modelo_final)
}
residuales_cuad <- residuales^2
# Tests
lb_res <- tryCatch(Box.test(as.numeric(residuales), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
lb_res2 <- tryCatch(Box.test(as.numeric(residuales_cuad), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
arch_res <- tryCatch(FinTS::ArchTest(as.numeric(residuales), lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
diag_final <- data.frame(Test = c("LjungBox resid.", "LjungBox resid.^2", "ARCH-LM"),
Stat = c(lb_res$statistic, lb_res2$statistic, arch_res$statistic),
P_value = c(lb_res$p.value, lb_res2$p.value, arch_res$p.value),
stringsAsFactors = FALSE)
write.csv(diag_final, file.path(res_folder, paste0(ticker, "_diag_final_modelo.csv")), row.names = FALSE)
print(diag_final) # o kable(diag_final)## Test Stat P_value
## 1 LjungBox resid. 13.62169 0.8491478
## 2 LjungBox resid.^2 6.39554 0.9982469
## 3 ARCH-LM 6.55452 0.9979106# Graficar ACF residuales y residuales^2
par(mfrow = c(2,1))
plot_acf_eviews(residuales, main = paste("ACF residuales -", mejor_modelo), lag.max = lag_test)
plot_acf_eviews(residuales_cuad, main = paste("ACF residuales^2 -", mejor_modelo), lag.max = lag_test)
## $fit
##
## Call:
## arima(x = x, order = c(p, 0, q), include.mean = include.mean, method = "ML")
##
## Coefficients:## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 ma4
## -0.3952 0.0277 -0.6041 -0.2138 0.5426 0.3325 -0.0159 0.5804 0.148
## s.e. NaN NaN NaN NaN NaN NaN NaN NaN NaN
## ma5 intercept
## -0.5285 0.0018
## s.e. NaN 0.0007
##
## sigma^2 estimated as 0.001023: log likelihood = 4063.3, aic = -8102.6
##
## $AIC
## [1] -8102.604## --- 13: Volatilidad y VaR según el modelo seleccionado ---
# Calcular volatilidad condicional y VaR; si el modelo es GARCH usamos sigma(); si no, usamos histórico simple
if (inherits(modelo_final, "uGARCHfit")) {
vol_cond <- tryCatch(sigma(modelo_final), error = function(e) NULL)
# Alinear fechas: asumir longitud vol_cond <= length(serie)
if (is.null(vol_cond)) stop("No se pudo extraer sigma() del modelo GARCH seleccionado.")
# extraer última sigma
ultima_sigma <- tail(as.numeric(vol_cond), 1)
dist_model <- modelo_final@model$modeldesc$distribution
shape_param <- if (dist_model == "std" && "shape" %in% names(coef(modelo_final))) coef(modelo_final)["shape"] else NULL
calcular_var_garch <- function(sigma_val, nivel_conf) {
alpha <- 1 - nivel_conf
if (!is.null(shape_param) && dist_model == "std") {
# qt gives quantile for central t; usamos qstudent adaptado: q = qt(alpha, df)
q <- qt(alpha, df = shape_param)
} else {
q <- qnorm(alpha)
}
return(as.numeric(sigma_val * q))
}
var_garch <- sapply(niveles_conf, function(p) calcular_var_garch(ultima_sigma, p))
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_garch,
VaR_percent = var_garch * 100,
Perdida = round(abs(var_garch) * investment_ref, 2),
stringsAsFactors = FALSE)
} else {
sd_hist <- sd(serie, na.rm = TRUE)
var_simple <- sd_hist * qnorm(1 - niveles_conf)
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_simple,
VaR_percent = var_simple * 100,
Perdida = round(abs(var_simple) * investment_ref, 2),
stringsAsFactors = FALSE)
}
write.csv(var_garch_df, file.path(res_folder, paste0(ticker, "_var_modelo.csv")), row.names = FALSE)
print(var_garch_df) # o kable(var_garch_df)## Nivel VaR_decimal VaR_percent Perdida
## 1 90% -0.03779698 -3.779698 3779.70
## 2 95% -0.05124874 -5.124874 5124.87
## 3 99% -0.08394774 -8.394774 8394.7714.Comparar esta medida de VaR con el VaR riskmetrics.
## --- 14: VaR RiskMetrics y comparación ---
rend_ts <- serie
n <- length(rend_ts)
var_rm <- numeric(n); var_rm[1] <- var(rend_ts, na.rm = TRUE)
for (t in 2:n) var_rm[t] <- lambda_rm * var_rm[t-1] + (1 - lambda_rm) * (rend_ts[t-1]^2)
sigma_rm <- sqrt(var_rm)
ultima_sigma_rm <- tail(sigma_rm, 1)
var_rm_vals <- sapply(niveles_conf, function(p) ultima_sigma_rm * qnorm(1 - p))
var_rm_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_RiskMetrics = var_rm_vals,
VaR_RM_percent = var_rm_vals * 100,
Perdida_RM = round(abs(var_rm_vals) * investment_ref, 2),
stringsAsFactors = FALSE)
write.csv(var_rm_df, file.path(res_folder, paste0(ticker, "_var_riskmetrics.csv")), row.names = FALSE)
print(var_rm_df) # o kable(var_rm_df)## Nivel VaR_RiskMetrics VaR_RM_percent Perdida_RM
## 1 90% -0.02883942 -2.883942 2883.94
## 2 95% -0.03701499 -3.701499 3701.50
## 3 99% -0.05235101 -5.235101 5235.10comparativa <- data.frame(Nivel = var_garch_df$Nivel,
VaR_Model_percent = var_garch_df$VaR_percent,
Perdida_Model = var_garch_df$Perdida,
VaR_RM_percent = var_rm_df$VaR_RM_percent,
Perdida_RM = var_rm_df$Perdida_RM,
stringsAsFactors = FALSE)
write.csv(comparativa, file.path(res_folder, paste0(ticker, "_comparativa_var.csv")), row.names = FALSE)
print(comparativa) # o kable(comparativa)## Nivel VaR_Model_percent Perdida_Model VaR_RM_percent Perdida_RM
## 1 90% -3.779698 3779.70 -2.883942 2883.94
## 2 95% -5.124874 5124.87 -3.701499 3701.50
## 3 99% -8.394774 8394.77 -5.235101 5235.10Vanguard Total Stock Market Index Fund
1. Serie de precios de Vanguard Total Stock Market Index Fund
pkgs <- c("quantmod","tidyquant","forecast","tseries","FinTS","rugarch","xts","zoo","dplyr",
"ggplot2","gridExtra","ggpubr","knitr","moments","scales","urca")
install_if_missing <- function(p) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
invisible(lapply(pkgs, install_if_missing))
library(quantmod); library(tidyquant); library(forecast); library(tseries); library(FinTS); library(rugarch)
library(xts); library(zoo); library(dplyr); library(ggplot2); library(gridExtra); library(ggpubr)
library(knitr); library(moments); library(scales); library(urca)
ticker <- "VTI"
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
maxP <- 5; maxQ <- 5; include.mean <- TRUE; lag_test <- 20
investment_ref <- 100000; lambda_rm <- 0.94
niveles_conf <- c(0.90, 0.95, 0.99)
res_folder <- "results_arma_garch_var_vti"
if (!dir.exists(res_folder)) dir.create(res_folder, recursive = TRUE)
prices_vti <- NULL
ok <- FALSE
# 1) Intento con quantmod (Yahoo)
try({
suppressWarnings(getSymbols(ticker, src = "yahoo", from = start_date, to = end_date, auto.assign = TRUE))
env_names <- ls(.GlobalEnv)
cand <- grep(paste0("\\^?", gsub("\\^","",ticker)), env_names, value = TRUE)
if (length(cand) > 0) {
symname <- cand[1]
tmp <- try(Ad(get(symname)), silent = TRUE)
if (!inherits(tmp, "try-error") && !is.null(tmp) && nrow(tmp) > 0) {
prices_vti <- na.omit(tmp); colnames(prices_vti) <- "VTI"; ok <- TRUE
message("Datos obtenidos con quantmod::getSymbols (Yahoo).")
}
}
}, silent = TRUE)
# 2) Fallback tidyquant
if (!ok) {
message("Fallback: tidyquant::tq_get() ...")
tq_try <- try(tq_get(ticker, get = "stock.prices", from = start_date, to = end_date), silent = TRUE)
if (!inherits(tq_try, "try-error") && !is.null(tq_try) && nrow(tq_try) > 0) {
prices_vti <- xts(tq_try$adjusted, order.by = as.Date(tq_try$date)); colnames(prices_vti) <- "VTI"; ok <- TRUE
message("Datos obtenidos con tidyquant::tq_get().")
}
}
# 3) Fallback CSV local (si existe)
if (!ok && file.exists("VTI_prices_local.csv")) {
message("Cargando CSV local 'VTI_prices_local.csv' ...")
tmp <- try(read.csv("VTI_prices_local.csv", stringsAsFactors = FALSE), silent = TRUE)
if (!inherits(tmp, "try-error")) {
nm <- names(tmp)
date_col <- nm[which(tolower(nm) == "date")[1]]
adj_col <- nm[which(tolower(nm) %in% c("adjusted","adj","close","close_adj","adjusted_close"))[1]]
if (!is.na(date_col) && !is.na(adj_col)) {
prices_vti <- xts(tmp[[adj_col]], order.by = as.Date(tmp[[date_col]])); colnames(prices_vti) <- "VTI"; prices_vti <- na.omit(prices_vti); ok <- TRUE
message("CSV local cargado.")
} else message("CSV local no tiene columnas esperadas (date + adjusted).")
} else message("Error leyendo CSV local.")
}
if (!ok) stop("No se pudo obtener datos para VTI. Comprueba conexión o proporciona 'VTI_prices_local.csv'.")
cat("VTI - observaciones (precios):", nrow(prices_vti), "\n")## VTI - observaciones (precios): 2009## VTI
## 2017-11-06 117.0599
## 2017-11-07 116.8576
## 2017-11-08 117.0511
## 2017-11-09 116.5761
## 2017-11-10 116.5937
## 2017-11-13 116.6553## VTI
## 2025-10-27 337.43
## 2025-10-28 337.95
## 2025-10-29 337.71
## 2025-10-30 334.08
## 2025-10-31 335.42
## 2025-11-03 335.802. Rendimientos logarÃtmicos y rendimientos^2
r_vti <- dailyReturn(prices_vti, type = "log"); colnames(r_vti) <- "r_VTI"
r_vti_sq <- r_vti^2; colnames(r_vti_sq) <- "r_VTI_sq"
vti_data <- merge(prices_vti, r_vti, r_vti_sq, join = "inner")
summary(vti_data)## Index VTI r_VTI r_VTI_sq
## Min. :2017-11-06 Min. :102.8 Min. :-0.1208225 Min. :0.000e+00
## 1st Qu.:2019-11-05 1st Qu.:137.7 1st Qu.:-0.0044231 1st Qu.:5.435e-06
## Median :2021-11-02 Median :192.9 Median : 0.0008549 Median :3.237e-05
## Mean :2021-11-03 Mean :194.4 Mean : 0.0005246 Mean :1.564e-04
## 3rd Qu.:2023-11-01 3rd Qu.:225.7 3rd Qu.: 0.0067622 3rd Qu.:1.155e-04
## Max. :2025-11-03 Max. :338.0 Max. : 0.0966333 Max. :1.460e-023. Gráficas — precio, rendimientos y rendimientos^2
df_vti <- fortify.zoo(vti_data, name = "date") %>% as.data.frame()
col_price <- "#238b45" # teal/verde oscuro
col_ret <- "#6a51a3" # púrpura
col_ret2 <- "#ef8a62" # naranja suave
alpha_shade <- 0.18
p_price <- ggplot(df_vti, aes(x = date, y = VTI)) +
geom_area(fill = col_price, alpha = alpha_shade) +
geom_line(color = col_price, size = 0.6) +
labs(title = "VTI — Precio ajustado", x = NULL, y = "Nivel") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
p_ret <- ggplot(df_vti, aes(x = date, y = r_VTI)) +
geom_hline(yintercept = 0, color = "grey70", linetype = "dashed") +
geom_line(color = col_ret, size = 0.5) +
geom_smooth(method = "loess", color = "black", se = FALSE, size = 0.4) +
labs(title = "VTI — Rendimientos logarÃtmicos diarios", x = NULL, y = "r_t") +
theme_light(base_size = 12) + theme(plot.title = element_text(face = "bold"))
p_ret2 <- ggplot(df_vti, aes(x = date, y = r_VTI_sq)) +
geom_col(fill = col_ret2, alpha = 0.75) +
labs(title = "VTI — Rendimientos al cuadrado (volatilidad instantánea)", x = NULL, y = "r_t^2") +
theme_classic(base_size = 12) + theme(plot.title = element_text(face = "bold"))
grid.arrange(p_price, p_ret, p_ret2, ncol = 1)
4. Correlogramas (ACF/PACF) para precio, r_t y r_t^2 (estilo simple)
plot_acf_eviews <- function(x, main = "", lag.max = 36, col = "#2c3e50") {
x_core <- as.numeric(na.omit(x))
acf_obj <- acf(x_core, lag.max = lag.max, plot = FALSE)
acf_vals <- acf_obj$acf[-1]; lags <- seq_along(acf_vals); ci <- 1.96/sqrt(length(x_core))
plot(0, type = "n", xlim = c(0.5, length(acf_vals)+0.5), ylim = c(min(acf_vals, -ci), max(acf_vals, ci)),
xlab = "Rezagos", ylab = "ACF", main = main, xaxt = "n")
axis(1, at = lags); abline(h = ci, col = "blue", lty = 2); abline(h = -ci, col = "blue", lty = 2); abline(h = 0)
for (i in lags) { v <- acf_vals[i]; rect(i-0.35, 0, i+0.35, v, col = adjustcolor(col, alpha.f = 0.6), border = NA) }
}
par(mfrow = c(3,2), mar = c(4,4,2,1))
plot_acf_eviews(vti_data$VTI, "ACF - VTI Precio", col = col_price)
pacf(coredata(vti_data$VTI), main = "PACF - VTI Precio", col = col_price)
plot_acf_eviews(vti_data$r_VTI, "ACF - VTI Rendimientos", col = col_ret)
pacf(coredata(vti_data$r_VTI), main = "PACF - VTI Rendimientos", col = col_ret)
plot_acf_eviews(vti_data$r_VTI_sq, "ACF - VTI Rendimientos^2", col = col_ret2)
pacf(coredata(vti_data$r_VTI_sq), main = "PACF - VTI Rendimientos^2", col = col_ret2)
5.Revisar la normalidad en cada una de las series del numeral 3 (histogramas + Jarque-Bera + QQ)
check_normality_plot <- function(x, name, fill="#d9f0d3", col_line="#238b45") {
x <- na.omit(as.numeric(x)); m <- mean(x); s <- sd(x)
p_hist <- ggplot(data.frame(x = x), aes(x = x)) +
geom_histogram(aes(y = ..density..), bins = 45, fill = fill, color = "white") +
geom_density(color = col_line, size = 0.8) +
stat_function(fun = dnorm, args = list(mean = m, sd = s), color = "red", linetype = "dashed") +
labs(title = paste(name, "- Histograma")) + theme_minimal()
p_qq <- ggqqplot(x, title = paste(name, "- QQ-plot"), ggtheme = theme_minimal()) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
list(hist = p_hist, qq = p_qq)
}
p_price <- check_normality_plot(coredata(vti_data$VTI), "VTI - Precio", fill="#edf8fb", col_line=col_price)
p_ret <- check_normality_plot(coredata(vti_data$r_VTI), "VTI - Rendimientos", fill="#f2e6f7", col_line=col_ret)
p_ret2 <- check_normality_plot(coredata(vti_data$r_VTI_sq), "VTI - Rendimientos^2", fill="#fff5eb", col_line=col_ret2)
grid.arrange(p_price$hist, p_price$qq, p_ret$hist, p_ret$qq, p_ret2$hist, p_ret2$qq, ncol = 2)
##
## Jarque Bera Test
##
## data: coredata(vti_data$r_VTI)
## X-squared = 16071, df = 2, p-value < 2.2e-166.Estimar la estacionariedad de cada una de las 3 series del numeral 3
# Usamos adf.test (Dickey-Fuller aumentado)
adf_price <- tryCatch(adf.test(coredata(vti_data$VTI)), error = function(e) NULL)
adf_ret <- tryCatch(adf.test(coredata(vti_data$r_VTI)), error = function(e) NULL)
adf_ret2 <- tryCatch(adf.test(coredata(vti_data$r_VTI_sq)), error = function(e) NULL)
cat("ADF - Precio: p-value =", if(!is.null(adf_price)) adf_price$p.value else NA, "\n")## ADF - Precio: p-value = 0.6504599## ADF - Rendimientos: p-value = 0.01## ADF - Rendimientos^2: p-value = 0.017. Seleccionarlas series que cumplan con la propiedad de estacionariedad
La prueba ADF muestra que la serie de precios no es estacionaria (p = 0.6505), mientras que tanto los rendimientos logarÃtmicos como sus cuadrados son estacionarios (p ≈ 0.01). Por tanto, procederemos a modelar y analizar los rendimientos logarÃtmicos diarios de VTI (r_t), ya que cumplen la condición de estacionariedad requerida para estimar modelos AR/MA/ARMA y modelos de volatilidad GARCH.
8. Con la serie seleccionada estimar los siguientes modelos: AR(p), MA(q), ARMA (p,q)
serie <- as.numeric(na.omit(coredata(vti_data$r_VTI)))
fit_safe <- function(p, q, x) {
tryCatch({
fit <- arima(x, order = c(p, 0, q), include.mean = include.mean, method = "ML")
list(fit = fit, AIC = AIC(fit))
}, error = function(e) list(fit = NULL, AIC = Inf))
}
best_ar <- list(AIC = Inf)
for (p in 1:maxP) { r <- fit_safe(p, 0, serie); if (r$AIC < best_ar$AIC) best_ar <- c(list(p = p), r) }
best_ma <- list(AIC = Inf)
for (q in 1:maxQ) { r <- fit_safe(0, q, serie); if (r$AIC < best_ma$AIC) best_ma <- c(list(q = q), r) }
best_arma <- list(AIC = Inf)
for (p in 1:maxP) for (q in 1:maxQ) { r <- fit_safe(p, q, serie); if (r$AIC < best_arma$AIC) best_arma <- c(list(p = p, q = q), r) }
model_ar <- best_ar$fit; best_ar_p <- best_ar$p
model_ma <- best_ma$fit; best_ma_q <- best_ma$q
model_arma <- best_arma$fit; best_arma_p <- best_arma$p; best_arma_q <- best_arma$q
aic_family <- data.frame(Modelo = c(paste0("AR(",best_ar_p,")"), paste0("MA(",best_ma_q,")"), paste0("ARMA(",best_arma_p,",",best_arma_q,")")),
AIC = c(AIC(model_ar), AIC(model_ma), AIC(model_arma)))
write.csv(aic_family, file.path(res_folder, "VTI_aic_family.csv"), row.names = FALSE)
kable(aic_family, caption = "AIC - Mejor AR/MA/ARMA por familia")| Modelo | AIC |
|---|---|
| AR(5) | -11955.14 |
| MA(4) | -11949.96 |
| ARMA(4,4) | -12025.21 |
9. En cada uno de los modelos estimados en el numeral 8, revisar la presencia de posibles problemas de autocorrelación y heterocedasticidad
run_model_diagnostics <- function(model, model_name, lag_test = 20) {
res <- residuals(model)
lb <- Box.test(res, lag = lag_test, type = "Ljung-Box")
arch <- tryCatch(FinTS::ArchTest(res, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
par(mfrow = c(2,1))
plot_acf_eviews(res, main = paste("ACF residuales -", model_name), col = "#2c3e50")
plot_acf_eviews(res^2, main = paste("ACF residuales^2 -", model_name), col = col_ret2)
par(mfrow = c(1,1))
data.frame(Model = model_name,
LjungBox_stat = round(as.numeric(lb$statistic),4), LjungBox_p = lb$p.value,
ARCH_stat = round(as.numeric(arch$statistic),4), ARCH_p = arch$p.value)
}
diag_ar <- run_model_diagnostics(model_ar, paste0("AR(",best_ar_p,")"))
diag_all <- bind_rows(diag_ar, diag_ma, diag_arma)
write.csv(diag_all, file.path(res_folder, "VTI_diag_models.csv"), row.names = FALSE)
kable(diag_all, caption = "Diagnóstico modelos")| Model | LjungBox_stat | LjungBox_p | ARCH_stat | ARCH_p | |
|---|---|---|---|---|---|
| Chi-squared…1 | AR(5) | 85.8694 | 0.0000000 | 541.4718 | 0 |
| Chi-squared…2 | MA(4) | 98.5683 | 0.0000000 | 537.2959 | 0 |
| Chi-squared…3 | ARMA(4,4) | 7.0505 | 0.9965161 | 511.9803 | 0 |
10. Si encuentre problemas mencionados en el numeral 9, proceder a estimar los mismos modelos pero bajo metodologÃa de estimación GARCH
need_garch <- any(diag_all$ARCH_p < 0.05, na.rm = TRUE)
fit_ar_garch <- fit_ma_garch <- fit_arma_garch <- NULL
if (need_garch) {
spec_ar <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(best_ar_p, 0), include.mean = include.mean),
distribution.model = "std")
spec_ma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0, best_ma_q), include.mean = include.mean),
distribution.model = "std")
spec_arma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = c(best_arma_p, best_arma_q), include.mean = include.mean),
distribution.model = "std")
fit_ar_garch <- tryCatch(ugarchfit(spec_ar, data = serie, solver = "hybrid"), error = function(e) NULL)
fit_ma_garch <- tryCatch(ugarchfit(spec_ma, data = serie, solver = "hybrid"), error = function(e) NULL)
fit_arma_garch <- tryCatch(ugarchfit(spec_arma, data = serie, solver = "hybrid"), error = function(e) NULL)
}
conv_df <- data.frame(Modelo = c("AR-GARCH","MA-GARCH","ARMA-GARCH"),
Converge = c(if(!is.null(fit_ar_garch)) fit_ar_garch@fit$convergence else NA,
if(!is.null(fit_ma_garch)) fit_ma_garch@fit$convergence else NA,
if(!is.null(fit_arma_garch)) fit_arma_garch@fit$convergence else NA))
write.csv(conv_df, file.path(res_folder, "VTI_garch_conv.csv"), row.names = FALSE)
kable(conv_df, caption = "Convergencia GARCH")| Modelo | Converge |
|---|---|
| AR-GARCH | 0 |
| MA-GARCH | 0 |
| ARMA-GARCH | 0 |
11. Seleccionar el modelo de mejor ajuste. Usar el criterio AIC
comparacion <- data.frame(Model = character(), AIC = numeric(), stringsAsFactors = FALSE)
comparacion <- rbind(comparacion, data.frame(Model = paste0("AR(",best_ar_p,")"), AIC = AIC(model_ar)))
comparacion <- rbind(comparacion, data.frame(Model = paste0("MA(",best_ma_q,")"), AIC = AIC(model_ma)))
comparacion <- rbind(comparacion, data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")"), AIC = AIC(model_arma)))
if (!is.null(fit_ar_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("AR(",best_ar_p,")-GARCH(1,1)"), AIC = infocriteria(fit_ar_garch)["Akaike",1]))
if (!is.null(fit_ma_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("MA(",best_ma_q,")-GARCH(1,1)"), AIC = infocriteria(fit_ma_garch)["Akaike",1]))
if (!is.null(fit_arma_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)"), AIC = infocriteria(fit_arma_garch)["Akaike",1]))
comparacion <- comparacion %>% arrange(AIC)
write.csv(comparacion, file.path(res_folder, "VTI_comparacion_aic.csv"), row.names = FALSE)
kable(comparacion, caption = "Comparación por AIC (menor = mejor)")| Model | AIC |
|---|---|
| ARMA(4,4) | -12025.211630 |
| AR(5) | -11955.137650 |
| MA(4) | -11949.958192 |
| ARMA(4,4)-GARCH(1,1) | -6.477592 |
| MA(4)-GARCH(1,1) | -6.473416 |
| AR(5)-GARCH(1,1) | -6.472454 |
## Mejor modelo por AIC: ARMA(4,4)12. Al mejor modelo, aplicar los test de heterocedasticidad y a través del correlograma de los residuales revisar si persiste el problema de autocorrelación.
objs <- list(model_ar, model_ma, model_arma, fit_ar_garch, fit_ma_garch, fit_arma_garch)
names(objs) <- c(paste0("AR(", best_ar_p, ")"),
paste0("MA(", best_ma_q, ")"),
paste0("ARMA(", best_arma_p, ",", best_arma_q, ")"),
paste0("AR(", best_ar_p, ")-GARCH(1,1)"),
paste0("MA(", best_ma_q, ")-GARCH(1,1)"),
paste0("ARMA(", best_arma_p, ",", best_arma_q, ")-GARCH(1,1)"))
modelo_final <- objs[[mejor_modelo]]
if (is.null(modelo_final)) stop("No se encontró el objeto del modelo seleccionado: ", mejor_modelo)
if (inherits(modelo_final, "uGARCHfit")) {
residuales <- residuals(modelo_final, standardize = TRUE); residuales_cuad <- residuales^2
} else {
residuales <- residuals(modelo_final); residuales_cuad <- residuales^2
}
lb_res <- Box.test(residuales, lag = lag_test, type = "Ljung-Box")
lb_res2 <- Box.test(residuales_cuad, lag = lag_test, type = "Ljung-Box")
arch_res <- tryCatch(FinTS::ArchTest(residuales, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
diag_final <- data.frame(Test = c("LjungBox resid.", "LjungBox resid.^2", "ARCH-LM"),
Stat = c(lb_res$statistic, lb_res2$statistic, arch_res$statistic),
P_value = c(lb_res$p.value, lb_res2$p.value, arch_res$p.value))
write.csv(diag_final, file.path(res_folder, "VTI_diag_final.csv"), row.names = FALSE)
kable(diag_final, caption = paste("Diagnóstico final -", mejor_modelo))| Test | Stat | P_value |
|---|---|---|
| LjungBox resid. | 7.050502 | 0.9965161 |
| LjungBox resid.^2 | 1724.163897 | 0.0000000 |
| ARCH-LM | 511.980323 | 0.0000000 |
par(mfrow = c(2,1))
plot_acf_eviews(residuales, main = paste("ACF residuales -", mejor_modelo), col = "#2c3e50")
plot_acf_eviews(residuales_cuad, main = paste("ACF residuales^2 -", mejor_modelo), col = col_ret2)
13-14. Con la ecuación del modelo seleccionado, calcular la volatilidad y posteriormente el VaR, y comparaciòn con Riskmetrics
if (inherits(modelo_final, "uGARCHfit")) {
vol_cond <- sigma(modelo_final)
idx <- index(r_vti); vol_xts <- xts(vol_cond, order.by = idx[(length(idx)-length(vol_cond)+1):length(idx)])
dist_model <- modelo_final@model$modeldesc$distribution
shape_param <- if (dist_model == "std") coef(modelo_final)["shape"] else NULL
calcular_var_garch <- function(sigma, nivel_conf) {
alpha <- 1 - nivel_conf
if (dist_model == "std" && !is.null(shape_param)) q <- qt(alpha, df = shape_param) else q <- qnorm(alpha)
return(sigma * q)
}
ultima_sigma <- tail(vol_cond, 1)
var_garch <- sapply(niveles_conf, function(p) calcular_var_garch(ultima_sigma, p))
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_garch,
VaR_percent = var_garch*100,
Perdida = round(abs(var_garch)*investment_ref, 2))
} else {
sd_hist <- sd(serie, na.rm = TRUE)
var_simple <- sd_hist * qnorm(1 - niveles_conf)
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_simple,
VaR_percent = var_simple*100,
Perdida = round(abs(var_simple)*investment_ref, 2))
}
write.csv(var_garch_df, file.path(res_folder, "VTI_var_modelo.csv"), row.names = FALSE)
kable(var_garch_df, caption = paste("VaR según modelo seleccionado -", mejor_modelo))| Nivel | VaR_decimal | VaR_percent | Perdida |
|---|---|---|---|
| 90% | -0.0160156 | -1.601562 | 1601.56 |
| 95% | -0.0205558 | -2.055583 | 2055.58 |
| 99% | -0.0290725 | -2.907250 | 2907.25 |
# RiskMetrics
rend_ts <- serie; n <- length(rend_ts)
var_rm <- numeric(n); var_rm[1] <- var(rend_ts, na.rm = TRUE)
for (t in 2:n) var_rm[t] <- lambda_rm * var_rm[t-1] + (1 - lambda_rm) * (rend_ts[t-1]^2)
sigma_rm <- sqrt(var_rm); ultima_sigma_rm <- tail(sigma_rm,1)
var_rm_vals <- sapply(niveles_conf, function(p) ultima_sigma_rm * qnorm(1 - p))
var_rm_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_RiskMetrics = var_rm_vals,
VaR_RM_percent = var_rm_vals*100,
Perdida_RM = round(abs(var_rm_vals)*investment_ref, 2))
write.csv(var_rm_df, file.path(res_folder, "VTI_var_riskmetrics.csv"), row.names = FALSE)
kable(var_rm_df, caption = "VaR RiskMetrics")| Nivel | VaR_RiskMetrics | VaR_RM_percent | Perdida_RM |
|---|---|---|---|
| 90% | -0.0099124 | -0.9912411 | 991.24 |
| 95% | -0.0127224 | -1.2722441 | 1272.24 |
| 99% | -0.0179936 | -1.7993591 | 1799.36 |
comparativa <- data.frame(Nivel = var_garch_df$Nivel,
VaR_Model_percent = var_garch_df$VaR_percent,
Perdida_Model = var_garch_df$Perdida,
VaR_RM_percent = var_rm_df$VaR_RM_percent,
Perdida_RM = var_rm_df$Perdida_RM)
write.csv(comparativa, file.path(res_folder, "VTI_comparativa_var.csv"), row.names = FALSE)
kable(comparativa, caption = "Comparativa VaR: Modelo vs RiskMetrics")| Nivel | VaR_Model_percent | Perdida_Model | VaR_RM_percent | Perdida_RM |
|---|---|---|---|---|
| 90% | -1.601562 | 1601.56 | -0.9912411 | 991.24 |
| 95% | -2.055583 | 2055.58 | -1.2722441 | 1272.24 |
| 99% | -2.907250 | 2907.25 | -1.7993591 | 1799.36 |
Tipo de cambio USD/MXN (USDMXN=X)
1. Serie de precios de Tipo de cambio USD/MXN
pkgs <- c("quantmod","tidyquant","forecast","tseries","FinTS","rugarch","xts","zoo","dplyr",
"ggplot2","gridExtra","ggpubr","knitr","moments","scales","urca")
install_if_missing <- function(p) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
invisible(lapply(pkgs, install_if_missing))
library(quantmod); library(tidyquant); library(forecast); library(tseries); library(FinTS); library(rugarch)
library(xts); library(zoo); library(dplyr); library(ggplot2); library(gridExtra); library(ggpubr)
library(knitr); library(moments); library(scales); library(urca)
ticker <- "USDMXN=X"
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
maxP <- 5; maxQ <- 5; include.mean <- TRUE; lag_test <- 20
investment_ref <- 100000; lambda_rm <- 0.94
niveles_conf <- c(0.90, 0.95, 0.99)
res_folder <- "results_arma_garch_var_usdmxn"
if (!dir.exists(res_folder)) dir.create(res_folder, recursive = TRUE)
prices_fx <- NULL
ok <- FALSE
# 1) Intento con quantmod (Yahoo)
try({
suppressWarnings(getSymbols(ticker, src = "yahoo", from = start_date, to = end_date, auto.assign = TRUE))
env_names <- ls(.GlobalEnv)
# Buscar nombre cargado (Yahoo suele crear "USDMXN=X" como nombre complicado, buscamos GLOB)
cand <- grep("USDMXN|USDMXN=X", env_names, value = TRUE)
if (length(cand) == 0) {
# alternativa: intentar con versión sin '='
cand <- grep("USDMXN", env_names, value = TRUE)
}
if (length(cand) > 0) {
symname <- cand[1]
tmp <- try(Ad(get(symname)), silent = TRUE)
if (!inherits(tmp, "try-error") && !is.null(tmp) && nrow(tmp) > 0) {
prices_fx <- na.omit(tmp); colnames(prices_fx) <- "USDMXN"; ok <- TRUE
message("Datos obtenidos con quantmod::getSymbols (Yahoo).")
}
}
}, silent = TRUE)
# 2) Fallback tidyquant
if (!ok) {
message("Fallback: tidyquant::tq_get() ...")
tq_try <- try(tq_get("USDMXN=X", get = "stock.prices", from = start_date, to = end_date), silent = TRUE)
if (!inherits(tq_try, "try-error") && !is.null(tq_try) && nrow(tq_try) > 0) {
prices_fx <- xts(tq_try$adjusted, order.by = as.Date(tq_try$date)); colnames(prices_fx) <- "USDMXN"; ok <- TRUE
message("Datos obtenidos con tidyquant::tq_get().")
}
}
# 3) Fallback CSV local
if (!ok && file.exists("USDMXN_prices_local.csv")) {
message("Cargando CSV local 'USDMXN_prices_local.csv' ...")
tmp <- try(read.csv("USDMXN_prices_local.csv", stringsAsFactors = FALSE), silent = TRUE)
if (!inherits(tmp, "try-error")) {
nm <- names(tmp)
date_col <- nm[which(tolower(nm) == "date")[1]]
adj_col <- nm[which(tolower(nm) %in% c("adjusted","adj","close","close_adj","adjusted_close","price"))[1]]
if (!is.na(date_col) && !is.na(adj_col)) {
prices_fx <- xts(tmp[[adj_col]], order.by = as.Date(tmp[[date_col]])); colnames(prices_fx) <- "USDMXN"; prices_fx <- na.omit(prices_fx); ok <- TRUE
message("CSV local cargado.")
} else message("CSV local no tiene columnas esperadas (date + adjusted).")
} else message("Error leyendo CSV local.")
}
if (!ok) stop("No se pudo obtener datos para USDMXN=X. Comprueba conexión o proporciona 'USDMXN_prices_local.csv'.")
cat("USDMXN - observaciones (precios):", nrow(prices_fx), "\n")## USDMXN - observaciones (precios): 2083## USDMXN
## 2017-11-06 19.04645
## 2017-11-07 19.13350
## 2017-11-08 19.09040
## 2017-11-09 19.02555
## 2017-11-10 19.10148
## 2017-11-13 19.10700## USDMXN
## 2025-10-28 18.38933
## 2025-10-29 18.42526
## 2025-10-30 18.46570
## 2025-10-31 18.52540
## 2025-11-03 18.55408
## 2025-11-04 18.491672. Rendimientos logarÃtmicos y rendimientos^2 (tipo de cambio)
# Para tipo de cambio, usamos retornos logarÃtmicos (log del tipo de cambio)
r_fx <- dailyReturn(prices_fx, type = "log"); colnames(r_fx) <- "r_USDMXN"
r_fx_sq <- r_fx^2; colnames(r_fx_sq) <- "r_USDMXN_sq"
fx_data <- merge(prices_fx, r_fx, r_fx_sq, join = "inner")
summary(fx_data)## Index USDMXN r_USDMXN r_USDMXN_sq
## Min. :2017-11-06 Min. :16.31 Min. :-3.775e-02 Min. :0.000e+00
## 1st Qu.:2019-11-05 1st Qu.:18.66 1st Qu.:-4.181e-03 1st Qu.:3.527e-06
## Median :2021-11-03 Median :19.45 Median :-5.591e-04 Median :1.612e-05
## Mean :2021-11-02 Mean :19.46 Mean :-1.419e-05 Mean :6.043e-05
## 3rd Qu.:2023-11-01 3rd Qu.:20.22 3rd Qu.: 3.639e-03 3rd Qu.:5.644e-05
## Max. :2025-11-04 Max. :25.34 Max. : 6.084e-02 Max. :3.702e-033. Gráficas — paletas distintas
df_fx <- fortify.zoo(fx_data, name = "date") %>% as.data.frame()
col_price <- "#0b3d91" # azul oscuro
col_ret <- "#ff7f50" # coral
col_ret2 <- "#2ca25f" # verde
alpha_shade <- 0.16
p_price <- ggplot(df_fx, aes(x = date, y = USDMXN)) +
geom_area(fill = col_price, alpha = alpha_shade) +
geom_line(color = col_price, size = 0.6) +
labs(title = "Tipo de cambio USD/MXN (USDMXN) — Nivel", x = NULL, y = "MXN por USD") +
theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold"))
p_ret <- ggplot(df_fx, aes(x = date, y = r_USDMXN)) +
geom_hline(yintercept = 0, color = "grey70", linetype = "dashed") +
geom_line(color = col_ret, size = 0.5) +
geom_smooth(method = "loess", color = "black", se = FALSE, size = 0.4) +
labs(title = "USD/MXN — Rendimientos logarÃtmicos diarios", x = NULL, y = "r_t") +
theme_light(base_size = 12) + theme(plot.title = element_text(face = "bold"))
p_ret2 <- ggplot(df_fx, aes(x = date, y = r_USDMXN_sq)) +
geom_col(fill = col_ret2, alpha = 0.75) +
labs(title = "USD/MXN — Rendimientos al cuadrado", x = NULL, y = "r_t^2") +
theme_classic(base_size = 12) + theme(plot.title = element_text(face = "bold"))
grid.arrange(p_price, p_ret, p_ret2, ncol = 1)
4. Correlogramas ACF / PACF
plot_acf_eviews <- function(x, main = "", lag.max = 36, col = "#2c3e50") {
x_core <- as.numeric(na.omit(x))
acf_obj <- acf(x_core, lag.max = lag.max, plot = FALSE)
acf_vals <- acf_obj$acf[-1]; lags <- seq_along(acf_vals); ci <- 1.96/sqrt(length(x_core))
plot(0, type = "n", xlim = c(0.5, length(acf_vals)+0.5), ylim = c(min(acf_vals, -ci), max(acf_vals, ci)),
xlab = "Rezagos", ylab = "ACF", main = main, xaxt = "n")
axis(1, at = lags); abline(h = ci, col = "blue", lty = 2); abline(h = -ci, col = "blue", lty = 2); abline(h = 0)
for (i in lags) { v <- acf_vals[i]; rect(i-0.35, 0, i+0.35, v, col = adjustcolor(col, alpha.f = 0.6), border = NA) }
}
par(mfrow = c(3,2), mar = c(4,4,2,1))
plot_acf_eviews(fx_data$USDMXN, "ACF - Precio (USDMXN)", col = col_price)
pacf(coredata(fx_data$USDMXN), main = "PACF - Precio (USDMXN)", col = col_price)
plot_acf_eviews(fx_data$r_USDMXN, "ACF - Rendimientos", col = col_ret)
pacf(coredata(fx_data$r_USDMXN), main = "PACF - Rendimientos", col = col_ret)
plot_acf_eviews(fx_data$r_USDMXN_sq, "ACF - Rendimientos^2", col = col_ret2)
pacf(coredata(fx_data$r_USDMXN_sq), main = "PACF - Rendimientos^2", col = col_ret2)
5. Normalidad: histogramas + QQ + Jarque‑Bera
check_normality_plot <- function(x, name, fill="#c6dbef", col_line="#0b3d91") {
x <- na.omit(as.numeric(x)); m <- mean(x); s <- sd(x)
p_hist <- ggplot(data.frame(x = x), aes(x = x)) +
geom_histogram(aes(y = ..density..), bins = 45, fill = fill, color = "white") +
geom_density(color = col_line, size = 0.8) +
stat_function(fun = dnorm, args = list(mean = m, sd = s), color = "red", linetype = "dashed") +
labs(title = paste(name, "- Histograma")) + theme_minimal()
p_qq <- ggqqplot(x, title = paste(name, "- QQ-plot"), ggtheme = theme_minimal()) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
list(hist = p_hist, qq = p_qq)
}
p1 <- check_normality_plot(coredata(fx_data$USDMXN), "USDMXN - Nivel")
p2 <- check_normality_plot(coredata(fx_data$r_USDMXN), "USDMXN - Rendimientos")
p3 <- check_normality_plot(coredata(fx_data$r_USDMXN_sq), "USDMXN - Rendimientos^2")
grid.arrange(p1$hist, p1$qq, p2$hist, p2$qq, p3$hist, p3$qq, ncol = 2)
##
## Jarque Bera Test
##
## data: coredata(fx_data$r_USDMXN)
## X-squared = 2211.9, df = 2, p-value < 2.2e-166. Estacionariedad — sólo Dickey–Fuller (ADF)
adf_price <- tryCatch(adf.test(coredata(fx_data$USDMXN)), error = function(e) NULL)
adf_ret <- tryCatch(adf.test(coredata(fx_data$r_USDMXN)), error = function(e) NULL)
adf_ret2 <- tryCatch(adf.test(coredata(fx_data$r_USDMXN_sq)), error = function(e) NULL)
cat("ADF - Precio: p-value =", if(!is.null(adf_price)) adf_price$p.value else NA, "\n")## ADF - Precio: p-value = 0.1388196## ADF - Rendimientos: p-value = 0.01## ADF - Rendimientos^2: p-value = 0.017. Seleccionar las series que cumplan con la propiedad de estacionariedad
Trabajamos con rendimientos logarÃtmicos simples (r_t = ln(P_t / P_{t-1})) porque convierten la serie de precios no estacionaria en una serie (aproximadamente) estacionaria, apropiada para modelar con AR/MA/ARMA y GARCH. Los rendimientos logarÃtmicos tienen varias propiedades útiles: son aditivos a lo largo del tiempo (la suma de rendimientos log en periodos sucesivos da el rendimiento logarÃtmico acumulado), facilitan la interpretación porcentual (r * 100 ≈ % cambio cuando los cambios son pequeños) y su distribución suele acomodar mejor los supuestos de modelos de series temporales (por ejemplo, homogeneidad en la escala y menor heterocedasticidad relativa). Además, la ADF mostró que la serie de rendimientos y sus cuadrados son estacionarias (p ≈ 0.01), confirmando que trabajar con r_t es la opción correcta para estimar modelos lineales y de volatilidad.
8. Estimación AR/MA/ARMA sobre rendimientos
serie <- as.numeric(na.omit(coredata(fx_data$r_USDMXN)))
fit_safe <- function(p, q, x) {
tryCatch({
fit <- arima(x, order = c(p, 0, q), include.mean = include.mean, method = "ML")
list(fit = fit, AIC = AIC(fit))
}, error = function(e) list(fit = NULL, AIC = Inf))
}
best_ar <- list(AIC = Inf)
for (p in 1:maxP) { r <- fit_safe(p, 0, serie); if (r$AIC < best_ar$AIC) best_ar <- c(list(p = p), r) }
best_ma <- list(AIC = Inf)
for (q in 1:maxQ) { r <- fit_safe(0, q, serie); if (r$AIC < best_ma$AIC) best_ma <- c(list(q = q), r) }
best_arma <- list(AIC = Inf)
for (p in 1:maxP) for (q in 1:maxQ) { r <- fit_safe(p, q, serie); if (r$AIC < best_arma$AIC) best_arma <- c(list(p = p, q = q), r) }
model_ar <- best_ar$fit; best_ar_p <- best_ar$p
model_ma <- best_ma$fit; best_ma_q <- best_ma$q
model_arma <- best_arma$fit; best_arma_p <- best_arma$p; best_arma_q <- best_arma$q
aic_family <- data.frame(Modelo = c(paste0("AR(",best_ar_p,")"), paste0("MA(",best_ma_q,")"), paste0("ARMA(",best_arma_p,",",best_arma_q,")")),
AIC = c(AIC(model_ar), AIC(model_ma), AIC(model_arma)))
write.csv(aic_family, file.path(res_folder, "USDMXN_aic_family.csv"), row.names = FALSE)
kable(aic_family, caption = "AIC - Mejor AR/MA/ARMA por familia")| Modelo | AIC |
|---|---|
| AR(1) | -14317.06 |
| MA(1) | -14317.06 |
| ARMA(3,2) | -14324.51 |
9. Diagnóstico: Ljung‑Box y ARCH‑LM
run_model_diagnostics <- function(model, model_name, lag_test = 20) {
res <- residuals(model)
lb <- Box.test(res, lag = lag_test, type = "Ljung-Box")
arch <- tryCatch(FinTS::ArchTest(res, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
par(mfrow = c(2,1))
plot_acf_eviews(res, main = paste("ACF residuales -", model_name), col = "#2c3e50")
plot_acf_eviews(res^2, main = paste("ACF residuales^2 -", model_name), col = col_ret2)
par(mfrow = c(1,1))
data.frame(Model = model_name,
LjungBox_stat = round(as.numeric(lb$statistic),4), LjungBox_p = lb$p.value,
ARCH_stat = round(as.numeric(arch$statistic),4), ARCH_p = arch$p.value)
}
diag_ar <- run_model_diagnostics(model_ar, paste0("AR(",best_ar_p,")"))
diag_all <- bind_rows(diag_ar, diag_ma, diag_arma)
write.csv(diag_all, file.path(res_folder, "USDMXN_diag_models.csv"), row.names = FALSE)
kable(diag_all, caption = "Diagnóstico modelos")| Model | LjungBox_stat | LjungBox_p | ARCH_stat | ARCH_p | |
|---|---|---|---|---|---|
| Chi-squared…1 | AR(1) | 27.9722 | 0.1100601 | 345.3548 | 0 |
| Chi-squared…2 | MA(1) | 27.9719 | 0.1100662 | 345.3572 | 0 |
| Chi-squared…3 | ARMA(3,2) | 11.7523 | 0.9243378 | 332.8155 | 0 |
10. Estimación GARCH(1,1) si procede
# ----- Parámetros de control -----
prefer_garch <- TRUE # si TRUE, preferir modelos GARCH cuando convergen y cumplen diagnósticos
require_lb2_non_sig <- TRUE # si TRUE, exigir LB(resid.^2) p > 0.05 para aceptar un GARCH candidato
garch_order <- c(1,1)
dist_garch <- "std" # "norm" o "std"
outfile_prefix <- "VTI" # prefijo para archivos (VTI)
# ----- 10: Estimar modelos GARCH si procede -----
need_garch <- any(diag_all$ARCH_p < 0.05, na.rm = TRUE) || prefer_garch
fit_ar_garch <- fit_ma_garch <- fit_arma_garch <- NULL
conv_ar_garch <- conv_ma_garch <- conv_arma_garch <- NA
if (need_garch) {
spec_ar <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_ar_p, 0), include.mean = include.mean),
distribution.model = dist_garch)
spec_ma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(0, best_ma_q), include.mean = include.mean),
distribution.model = dist_garch)
spec_arma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_arma_p, best_arma_q), include.mean = include.mean),
distribution.model = dist_garch)
fit_ar_garch <- tryCatch(ugarchfit(spec_ar, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("AR-GARCH error: ", e$message); NULL })
fit_ma_garch <- tryCatch(ugarchfit(spec_ma, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("MA-GARCH error: ", e$message); NULL })
fit_arma_garch <- tryCatch(ugarchfit(spec_arma, data = serie, solver = "hybrid", silent = TRUE), error = function(e) { message("ARMA-GARCH error: ", e$message); NULL })
conv_ar_garch <- if (!is.null(fit_ar_garch)) tryCatch(fit_ar_garch@fit$convergence, error = function(e) NA) else NA
conv_ma_garch <- if (!is.null(fit_ma_garch)) tryCatch(fit_ma_garch@fit$convergence, error = function(e) NA) else NA
conv_arma_garch <- if (!is.null(fit_arma_garch)) tryCatch(fit_arma_garch@fit$convergence, error = function(e) NA) else NA
}
conv_df <- data.frame(Modelo = c("AR-GARCH","MA-GARCH","ARMA-GARCH"),
Converge = c(conv_ar_garch, conv_ma_garch, conv_arma_garch),
stringsAsFactors = FALSE)
write.csv(conv_df, file.path(res_folder, paste0(outfile_prefix, "_garch_conv.csv")), row.names = FALSE)
print(conv_df) # o kable(conv_df) en Rmd ## Modelo Converge
## 1 AR-GARCH 0
## 2 MA-GARCH 0
## 3 ARMA-GARCH 011. Selección por AIC entre modelos
# ----- 11: Comparación por AIC + diagnósticos -----
library(dplyr)
comparacion <- data.frame(Model = character(), AIC = numeric(), Converge = logical(),
ARCH_p = numeric(), LB_p_resid = numeric(), LB_p_resid2 = numeric(),
IsGARCH = logical(), stringsAsFactors = FALSE)
# helper diagnóstico
diagnosticar_modelo <- function(obj) {
out <- list(converge = NA, arch_p = NA, lb_p = NA, lb_p2 = NA)
if (is.null(obj)) return(out)
if (inherits(obj, "uGARCHfit")) {
out$converge <- tryCatch(obj@fit$convergence == 0, error = function(e) NA)
res_std <- tryCatch(residuals(obj, standardize = TRUE), error = function(e) residuals(obj))
res_plain <- tryCatch(residuals(obj), error = function(e) res_std)
out$lb_p <- tryCatch(Box.test(as.numeric(res_std), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$lb_p2 <- tryCatch(Box.test(as.numeric(res_std^2), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_std), lags = lag_test)$p.value, error = function(e) NA)
} else {
out$converge <- TRUE
res_plain <- tryCatch(residuals(obj), error = function(e) NULL)
if (is.null(res_plain)) return(out)
out$lb_p <- tryCatch(Box.test(as.numeric(res_plain), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$lb_p2 <- tryCatch(Box.test(as.numeric(res_plain^2), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_plain), lags = lag_test)$p.value, error = function(e) NA)
}
return(out)
}
# añadir modelos AR/MA/ARMA (si existen)
if (exists("model_ar") && !is.null(model_ar)) {
di <- diagnosticar_modelo(model_ar)
comparacion <- rbind(comparacion,
data.frame(Model = paste0("AR(",best_ar_p,")"),
AIC = AIC(model_ar),
Converge = di$converge,
ARCH_p = di$arch_p,
LB_p_resid = di$lb_p,
LB_p_resid2 = di$lb_p2,
IsGARCH = FALSE,
stringsAsFactors = FALSE))
}
if (exists("model_ma") && !is.null(model_ma)) {
di <- diagnosticar_modelo(model_ma)
comparacion <- rbind(comparacion,
data.frame(Model = paste0("MA(",best_ma_q,")"),
AIC = AIC(model_ma),
Converge = di$converge,
ARCH_p = di$arch_p,
LB_p_resid = di$lb_p,
LB_p_resid2 = di$lb_p2,
IsGARCH = FALSE,
stringsAsFactors = FALSE))
}
if (exists("model_arma") && !is.null(model_arma)) {
di <- diagnosticar_modelo(model_arma)
comparacion <- rbind(comparacion,
data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")"),
AIC = AIC(model_arma),
Converge = di$converge,
ARCH_p = di$arch_p,
LB_p_resid = di$lb_p,
LB_p_resid2 = di$lb_p2,
IsGARCH = FALSE,
stringsAsFactors = FALSE))
}
# añadir entradas GARCH si existen
add_garch_entry <- function(fit_obj, label) {
if (is.null(fit_obj)) return(NULL)
aicv <- tryCatch(infocriteria(fit_obj)["Akaike",1], error = function(e) NA)
di <- diagnosticar_modelo(fit_obj)
data.frame(Model = label, AIC = aicv, Converge = di$converge,
ARCH_p = di$arch_p, LB_p_resid = di$lb_p, LB_p_resid2 = di$lb_p2,
IsGARCH = TRUE, stringsAsFactors = FALSE)
}
if (!is.null(fit_ar_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_ar_garch, paste0("AR(",best_ar_p,")-GARCH(1,1)")))
if (!is.null(fit_ma_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_ma_garch, paste0("MA(",best_ma_q,")-GARCH(1,1)")))
if (!is.null(fit_arma_garch)) comparacion <- rbind(comparacion, add_garch_entry(fit_arma_garch, paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)")))
comparacion <- comparacion %>% arrange(AIC)
write.csv(comparacion, file.path(res_folder, paste0(outfile_prefix, "_comparacion_aic_diag.csv")), row.names = FALSE)
print(comparacion) # o kable(comparacion) ## Model AIC Converge ARCH_p
## Chi-squared2 ARMA(3,2) -14324.511869 TRUE 1.530456e-58
## Chi-squared1 MA(1) -14317.056077 TRUE 4.028252e-61
## Chi-squared AR(1) -14317.056025 TRUE 4.032857e-61
## Chi-squared5 ARMA(3,2)-GARCH(1,1) -7.160773 TRUE 8.194907e-01
## Chi-squared4 MA(1)-GARCH(1,1) -7.158108 TRUE 8.378771e-01
## Chi-squared3 AR(1)-GARCH(1,1) -7.158098 TRUE 8.377194e-01
## LB_p_resid LB_p_resid2 IsGARCH
## Chi-squared2 0.9243378 0.0000000 FALSE
## Chi-squared1 0.1100662 0.0000000 FALSE
## Chi-squared 0.1100601 0.0000000 FALSE
## Chi-squared5 0.7848813 0.8159798 TRUE
## Chi-squared4 0.6187248 0.8296671 TRUE
## Chi-squared3 0.6193327 0.8295042 TRUE# Selección final
mejor_modelo <- NA
if (prefer_garch) {
# candidatos GARCH que convergieron
garch_cands <- comparacion %>% filter(IsGARCH == TRUE & (Converge == TRUE | Converge == 0))
if (nrow(garch_cands) > 0) {
if (require_lb2_non_sig) {
garch_cands <- garch_cands %>% filter(!is.na(LB_p_resid2) & LB_p_resid2 > 0.05)
}
if (nrow(garch_cands) > 0) {
mejor_modelo <- garch_cands$Model[which.min(garch_cands$AIC)]
message("Se seleccionó modelo GARCH candidato: ", mejor_modelo)
} else {
mejor_modelo <- comparacion$Model[1]
message("No hay GARCH candidatos que cumplan requisitos diagnósticos; se selecciona por AIC global: ", mejor_modelo)
}
} else {
mejor_modelo <- comparacion$Model[1]
message("No hay modelos GARCH convergidos; se selecciona por AIC global: ", mejor_modelo)
}
} else {
mejor_modelo <- comparacion$Model[1]
message("Selección por AIC: ", mejor_modelo)
} 12. Diagnóstico final del mejor modelo
# construir lista de objetos
objs <- list()
if (exists("model_ar")) objs[[paste0("AR(",best_ar_p,")")]] <- model_ar
if (exists("model_ma")) objs[[paste0("MA(",best_ma_q,")")]] <- model_ma
if (exists("model_arma")) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")")]] <- model_arma
if (!is.null(fit_ar_garch)) objs[[paste0("AR(",best_ar_p,")-GARCH(1,1)")]] <- fit_ar_garch
if (!is.null(fit_ma_garch)) objs[[paste0("MA(",best_ma_q,")-GARCH(1,1)")]] <- fit_ma_garch
if (!is.null(fit_arma_garch)) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)")]] <- fit_arma_garch
if (!mejor_modelo %in% names(objs)) stop("No se encontró objeto para el modelo seleccionado: ", mejor_modelo)
modelo_final <- objs[[mejor_modelo]]
# residuos (estandarizados si GARCH)
if (inherits(modelo_final, "uGARCHfit")) {
residuales <- tryCatch(residuals(modelo_final, standardize = TRUE), error = function(e) residuals(modelo_final))
} else {
residuales <- residuals(modelo_final)
}
residuales_cuad <- residuales^2
lb_res <- tryCatch(Box.test(as.numeric(residuales), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
lb_res2 <- tryCatch(Box.test(as.numeric(residuales_cuad), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
arch_res <- tryCatch(FinTS::ArchTest(as.numeric(residuales), lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
diag_final <- data.frame(Test = c("LjungBox resid.", "LjungBox resid.^2", "ARCH-LM"),
Stat = c(lb_res$statistic, lb_res2$statistic, arch_res$statistic),
P_value = c(lb_res$p.value, lb_res2$p.value, arch_res$p.value),
stringsAsFactors = FALSE)
write.csv(diag_final, file.path(res_folder, paste0(outfile_prefix, "_diag_final.csv")), row.names = FALSE)
print(diag_final) # o kable(diag_final)## Test Stat P_value
## 1 LjungBox resid. 14.85078 0.7848813
## 2 LjungBox resid.^2 14.28077 0.8159798
## 3 ARCH-LM 14.21383 0.8194907# Graficar ACF (residuos y residuos^2)
par(mfrow = c(2,1))
plot_acf_eviews(residuales, main = paste("ACF residuales -", mejor_modelo), col = "#2c3e50")
plot_acf_eviews(residuales_cuad, main = paste("ACF residuales^2 -", mejor_modelo), col = col_ret2)
13-14. Volatilidad condicional y VaR (modelo vs RiskMetrics)
# ----- 13: Volatilidad condicional y VaR del modelo seleccionado -----
if (inherits(modelo_final, "uGARCHfit")) {
vol_cond <- sigma(modelo_final)
# alinear Ãndices con la serie original (serie_xts/nombres) si dispones de Ãndices; aquà tomo tail
ultima_sigma <- tail(as.numeric(vol_cond), 1)
dist_model <- modelo_final@model$modeldesc$distribution
shape_param <- if (dist_model == "std" && "shape" %in% names(coef(modelo_final))) coef(modelo_final)["shape"] else NULL
calcular_var_garch <- function(sigma_val, nivel_conf) {
alpha <- 1 - nivel_conf
if (!is.null(shape_param) && dist_model == "std") q <- qt(alpha, df = shape_param) else q <- qnorm(alpha)
sigma_val * q
}
var_garch <- sapply(niveles_conf, function(p) calcular_var_garch(ultima_sigma, p))
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_garch,
VaR_percent = var_garch * 100,
Perdida = round(abs(var_garch) * investment_ref, 2),
stringsAsFactors = FALSE)
} else {
sd_hist <- sd(serie, na.rm = TRUE)
var_simple <- sd_hist * qnorm(1 - niveles_conf)
var_garch_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_decimal = var_simple,
VaR_percent = var_simple * 100,
Perdida = round(abs(var_simple) * investment_ref, 2),
stringsAsFactors = FALSE)
}
write.csv(var_garch_df, file.path(res_folder, paste0(outfile_prefix, "_var_modelo.csv")), row.names = FALSE)
print(var_garch_df) # o kable(var_garch_df)## Nivel VaR_decimal VaR_percent Perdida
## 1 90% -0.006531400 -0.6531400 653.14
## 2 95% -0.008820447 -0.8820447 882.04
## 3 99% -0.014288749 -1.4288749 1428.87# ----- 14: VaR RiskMetrics y comparación -----
rend_ts <- serie
n <- length(rend_ts)
var_rm <- numeric(n); var_rm[1] <- var(rend_ts, na.rm = TRUE)
for (t in 2:n) var_rm[t] <- lambda_rm * var_rm[t-1] + (1 - lambda_rm) * (rend_ts[t-1]^2)
sigma_rm <- sqrt(var_rm); ultima_sigma_rm <- tail(sigma_rm, 1)
var_rm_vals <- sapply(niveles_conf, function(p) ultima_sigma_rm * qnorm(1 - p))
var_rm_df <- data.frame(Nivel = paste0(niveles_conf*100, "%"),
VaR_RiskMetrics = var_rm_vals,
VaR_RM_percent = var_rm_vals * 100,
Perdida_RM = round(abs(var_rm_vals) * investment_ref, 2),
stringsAsFactors = FALSE)
write.csv(var_rm_df, file.path(res_folder, paste0(outfile_prefix, "_var_riskmetrics.csv")), row.names = FALSE)
print(var_rm_df) # o kable(var_rm_df)## Nivel VaR_RiskMetrics VaR_RM_percent Perdida_RM
## 1 90% -0.003810370 -0.3810370 381.04
## 2 95% -0.004890557 -0.4890557 489.06
## 3 99% -0.006916808 -0.6916808 691.68comparativa <- data.frame(Nivel = var_garch_df$Nivel,
VaR_Model_percent = var_garch_df$VaR_percent,
Perdida_Model = var_garch_df$Perdida,
VaR_RM_percent = var_rm_df$VaR_RM_percent,
Perdida_RM = var_rm_df$Perdida_RM,
stringsAsFactors = FALSE)
write.csv(comparativa, file.path(res_folder, paste0(outfile_prefix, "_comparativa_var.csv")), row.names = FALSE)
print(comparativa) # o kable(comparativa)## Nivel VaR_Model_percent Perdida_Model VaR_RM_percent Perdida_RM
## 1 90% -0.6531400 653.14 -0.3810370 381.04
## 2 95% -0.8820447 882.04 -0.4890557 489.06
## 3 99% -1.4288749 1428.87 -0.6916808 691.68Rendimiento Treasury 10y (^TNX) — ADF
1. Serie de rendimientos de Rendimiento Treasury 10y (^TNX) — ADF
pkgs <- c("quantmod","tidyquant","forecast","tseries","FinTS","rugarch","xts","zoo","dplyr",
"ggplot2","gridExtra","ggpubr","knitr","moments","scales","urca")
install_if_missing <- function(p) if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
invisible(lapply(pkgs, install_if_missing))
library(quantmod); library(tidyquant); library(forecast); library(tseries); library(FinTS); library(rugarch)
library(xts); library(zoo); library(dplyr); library(ggplot2); library(gridExtra); library(ggpubr)
library(knitr); library(moments); library(scales); library(urca)
ticker <- "^TNX"
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
maxP <- 5; maxQ <- 5; include.mean <- TRUE; lag_test <- 20
investment_ref <- 100000 # valor nominal para convertir VaR a pérdidas monetarias
Duration_mod <- 8 # duración modificada por defecto (ajusta si quieres)
lambda_rm <- 0.94
niveles_conf <- c(0.90, 0.95, 0.99)
res_folder <- "results_arma_garch_var_tnx"
if (!dir.exists(res_folder)) dir.create(res_folder, recursive = TRUE)
prices_tnx <- NULL
ok <- FALSE
# 1) quantmod (Yahoo)
try({
suppressWarnings(getSymbols(ticker, src = "yahoo", from = start_date, to = end_date, auto.assign = TRUE))
env_names <- ls(.GlobalEnv)
cand <- grep("TNX|\\^TNX", env_names, value = TRUE)
if (length(cand) > 0) {
symname <- cand[1]
tmp <- try(Ad(get(symname)), silent = TRUE)
if (!inherits(tmp, "try-error") && !is.null(tmp) && nrow(tmp) > 0) {
prices_tnx <- na.omit(tmp); colnames(prices_tnx) <- "TNX"; ok <- TRUE
message("Datos obtenidos con quantmod::getSymbols (Yahoo).")
}
}
}, silent = TRUE)
# 2) tidyquant (fallback)
if (!ok) {
message("Fallback: tidyquant::tq_get() ...")
tq_try <- try(tq_get("^TNX", get = "stock.prices", from = start_date, to = end_date), silent = TRUE)
if (!inherits(tq_try, "try-error") && !is.null(tq_try) && nrow(tq_try) > 0) {
prices_tnx <- xts(tq_try$adjusted, order.by = as.Date(tq_try$date)); colnames(prices_tnx) <- "TNX"; ok <- TRUE
message("Datos obtenidos con tidyquant::tq_get().")
}
}
# 3) CSV local
if (!ok && file.exists("TNX_prices_local.csv")) {
message("Cargando CSV local 'TNX_prices_local.csv' ...")
tmp <- try(read.csv("TNX_prices_local.csv", stringsAsFactors = FALSE), silent = TRUE)
if (!inherits(tmp, "try-error")) {
nm <- names(tmp)
date_col <- nm[which(tolower(nm) == "date")[1]]
adj_col <- nm[which(tolower(nm) %in% c("adjusted","adj","close","close_adj","adjusted_close","price","value"))[1]]
if (!is.na(date_col) && !is.na(adj_col)) {
prices_tnx <- xts(tmp[[adj_col]], order.by = as.Date(tmp[[date_col]])); colnames(prices_tnx) <- "TNX"; prices_tnx <- na.omit(prices_tnx); ok <- TRUE
message("CSV local cargado.")
} else message("CSV local no tiene columnas esperadas (date + adjusted).")
} else message("Error leyendo CSV local.")
}
if (!ok) stop("No se pudo obtener datos para ^TNX. Comprueba conexión o proporciona 'TNX_prices_local.csv'.")
cat("^TNX - observaciones (yield):", nrow(prices_tnx), "\n")## ^TNX - observaciones (yield): 2009## TNX
## 2017-11-06 2.320
## 2017-11-07 2.307
## 2017-11-08 2.325
## 2017-11-09 2.331
## 2017-11-10 2.400
## 2017-11-13 2.400## TNX
## 2025-10-27 3.997
## 2025-10-28 3.983
## 2025-10-29 4.058
## 2025-10-30 4.093
## 2025-10-31 4.101
## 2025-11-03 4.1062. Transformación: primera diferencia de la yield (Δy_t) y Δy_t^2
# TNX viene expresado en % (por ejemplo 3.5 = 3.5%)
# Convertimos a decimal (opcional) o trabajamos en puntos porcentuales. Aquà trabajo en unidades de porcentaje como están.
y_t <- prices_tnx$TNX
dy <- diff(y_t) # Δy_t = y_t - y_{t-1}
dy <- na.omit(dy)
dy_sq <- dy^2
tnx_data <- merge(y_t, dy, dy_sq, join = "inner")
colnames(tnx_data) <- c("TNX", "dTNX", "dTNX_sq")
summary(tnx_data)## Index TNX dTNX dTNX_sq
## Min. :2017-11-07 Min. :0.499 Min. :-0.3220000 Min. :0.000000
## 1st Qu.:2019-11-05 1st Qu.:1.659 1st Qu.:-0.0320001 1st Qu.:0.000289
## Median :2021-11-02 Median :2.876 Median : 0.0000000 Median :0.001089
## Mean :2021-11-04 Mean :2.811 Mean : 0.0008894 Mean :0.003253
## 3rd Qu.:2023-11-01 3rd Qu.:4.026 3rd Qu.: 0.0340000 3rd Qu.:0.003249
## Max. :2025-11-03 Max. :4.988 Max. : 0.2690000 Max. :0.1036843. Gráficas (yield, cambios, cambios^2)
df_tnx <- fortify.zoo(tnx_data, name = "date") %>% as.data.frame()
col_price <- "#2b8cbe"
col_ret <- "#f03b20"
col_ret2 <- "#7fc97f"
alpha_shade <- 0.16
p_price <- ggplot(df_tnx, aes(x = date, y = TNX)) +
geom_line(color = col_price, size = 0.6) +
labs(title = "Yield 10y (^TNX) — Nivel (%)", x = NULL, y = "Yield (%)") + theme_minimal()
p_dy <- ggplot(df_tnx, aes(x = date, y = dTNX)) +
geom_hline(yintercept = 0, color = "grey70", linetype = "dashed") +
geom_line(color = col_ret, size = 0.5) +
labs(title = "ΔYield (y_t - y_{t-1}) — Cambios diarios", x = NULL, y = "ΔYield (pp)") + theme_light()
p_dy2 <- ggplot(df_tnx, aes(x = date, y = dTNX_sq)) +
geom_col(fill = col_ret2, alpha = 0.75) +
labs(title = "ΔYield^2 — Indicador de volatilidad", x = NULL, y = "ΔYield^2") + theme_classic()
grid.arrange(p_price, p_dy, p_dy2, ncol = 1)
4. Correlogramas ACF / PACF (nivel y cambios)
plot_acf_eviews <- function(x, main = "", lag.max = 36, col = "#2c3e50") {
x_core <- as.numeric(na.omit(x))
acf_obj <- acf(x_core, lag.max = lag.max, plot = FALSE)
acf_vals <- acf_obj$acf[-1]; lags <- seq_along(acf_vals); ci <- 1.96/sqrt(length(x_core))
plot(0, type = "n", xlim = c(0.5, length(acf_vals)+0.5), ylim = c(min(acf_vals, -ci), max(acf_vals, ci)),
xlab = "Rezagos", ylab = "ACF", main = main, xaxt = "n")
axis(1, at = lags); abline(h = ci, col = "blue", lty = 2); abline(h = -ci, col = "blue", lty = 2); abline(h = 0)
for (i in lags) { v <- acf_vals[i]; rect(i-0.35, 0, i+0.35, v, col = adjustcolor(col, alpha.f = 0.6), border = NA) }
}
par(mfrow = c(3,2), mar = c(4,4,2,1))
plot_acf_eviews(tnx_data$TNX, "ACF - TNX Nivel", col = col_price)
pacf(coredata(tnx_data$TNX), main = "PACF - TNX Nivel", col = col_price)
plot_acf_eviews(tnx_data$dTNX, "ACF - ΔYield", col = col_ret)
pacf(coredata(tnx_data$dTNX), main = "PACF - ΔYield", col = col_ret)
plot_acf_eviews(tnx_data$dTNX_sq, "ACF - ΔYield^2", col = col_ret2)
pacf(coredata(tnx_data$dTNX_sq), main = "PACF - ΔYield^2", col = col_ret2)
5. Normalidad y Jarque‑Bera sobre ΔYield
x <- coredata(tnx_data$dTNX)
p_hist <- ggplot(data.frame(x = x), aes(x = x)) +
geom_histogram(aes(y = ..density..), bins = 45, fill = "#edf8fb", color = "white") +
geom_density(color = "#2b8cbe", size = 0.8) +
labs(title = "ΔYield - Histograma y densidad") + theme_minimal()
p_qq <- ggqqplot(x, title = "ΔYield - QQ-plot", ggtheme = theme_minimal()) + geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
grid.arrange(p_hist, p_qq, ncol = 2)
##
## Jarque Bera Test
##
## data: x
## X-squared = 350.29, df = 2, p-value < 2.2e-166. Estacionariedad — sólo Dickey–Fuller (ADF)
adf_level <- tryCatch(adf.test(coredata(tnx_data$TNX)), error = function(e) NULL)
adf_dy <- tryCatch(adf.test(coredata(tnx_data$dTNX)), error = function(e) NULL)
adf_dy2 <- tryCatch(adf.test(coredata(tnx_data$dTNX_sq)), error = function(e) NULL)
cat("ADF - Nivel TNX: p-value =", if(!is.null(adf_level)) adf_level$p.value else NA, "\n")## ADF - Nivel TNX: p-value = 0.8032471## ADF - ΔYield: p-value = 0.01## ADF - ΔYield^2: p-value = 0.017. Seleccionar las series que cumplan con la propiedad de estacionariedad
Dado que la prueba ADF muestra que el nivel del yield (^TNX) no es estacionario (p = 0.803) pero su primera diferencia ΔYield y ΔYield^2 sà lo son (p ≈ 0.01), modelamos ΔYield (cambios absolutos en puntos porcentuales) con AR/MA/ARMA; si además el test ARCH sobre los residuos resulta significativo, estimamos un AR/MA/ARMA–GARCH(1,1) (distribución t si hay colas pesadas) — se selecciona el modelo final por AIC y se convierten los VaR en pérdidas monetarias usando duración modificada (Loss ≈ −Duration_mod × Δy_decimal × Notional).
8. Estimación AR/MA/ARMA sobre ΔYield
serie <- as.numeric(na.omit(coredata(tnx_data$dTNX)))
fit_safe <- function(p, q, x) {
tryCatch({
fit <- arima(x, order = c(p, 0, q), include.mean = include.mean, method = "ML")
list(fit = fit, AIC = AIC(fit))
}, error = function(e) list(fit = NULL, AIC = Inf))
}
best_ar <- list(AIC = Inf)
for (p in 1:maxP) { r <- fit_safe(p, 0, serie); if (r$AIC < best_ar$AIC) best_ar <- c(list(p = p), r) }
best_ma <- list(AIC = Inf)
for (q in 1:maxQ) { r <- fit_safe(0, q, serie); if (r$AIC < best_ma$AIC) best_ma <- c(list(q = q), r) }
best_arma <- list(AIC = Inf)
for (p in 1:maxP) for (q in 1:maxQ) { r <- fit_safe(p, q, serie); if (r$AIC < best_arma$AIC) best_arma <- c(list(p = p, q = q), r) }
model_ar <- best_ar$fit; best_ar_p <- best_ar$p
model_ma <- best_ma$fit; best_ma_q <- best_ma$q
model_arma <- best_arma$fit; best_arma_p <- best_arma$p; best_arma_q <- best_arma$q
aic_family <- data.frame(Modelo = c(paste0("AR(",best_ar_p,")"), paste0("MA(",best_ma_q,")"), paste0("ARMA(",best_arma_p,",",best_arma_q,")")),
AIC = c(AIC(model_ar), AIC(model_ma), AIC(model_arma)))
write.csv(aic_family, file.path(res_folder, "TNX_aic_family.csv"), row.names = FALSE)
kable(aic_family, caption = "AIC - Mejor AR/MA/ARMA por familia")| Modelo | AIC |
|---|---|
| AR(2) | -5803.033 |
| MA(2) | -5802.977 |
| ARMA(2,1) | -5801.730 |
9. Diagnóstico (Ljung‑Box y ARCH‑LM)
run_model_diagnostics <- function(model, model_name, lag_test = 20) {
res <- residuals(model)
lb <- Box.test(res, lag = lag_test, type = "Ljung-Box")
arch <- tryCatch(FinTS::ArchTest(res, lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
par(mfrow = c(2,1))
plot_acf_eviews(res, main = paste("ACF residuales -", model_name), col = "#2c3e50")
plot_acf_eviews(res^2, main = paste("ACF residuales^2 -", model_name), col = col_ret2)
par(mfrow = c(1,1))
data.frame(Model = model_name,
LjungBox_stat = round(as.numeric(lb$statistic),4), LjungBox_p = lb$p.value,
ARCH_stat = round(as.numeric(arch$statistic),4), ARCH_p = arch$p.value)
}
diag_ar <- run_model_diagnostics(model_ar, paste0("AR(",best_ar_p,")"))
diag_all <- bind_rows(diag_ar, diag_ma, diag_arma)
write.csv(diag_all, file.path(res_folder, "TNX_diag_models.csv"), row.names = FALSE)
kable(diag_all, caption = "Diagnóstico modelos")| Model | LjungBox_stat | LjungBox_p | ARCH_stat | ARCH_p | |
|---|---|---|---|---|---|
| Chi-squared…1 | AR(2) | 16.0185 | 0.7154777 | 279.8127 | 0 |
| Chi-squared…2 | MA(2) | 16.1609 | 0.7065908 | 279.8566 | 0 |
| Chi-squared…3 | ARMA(2,1) | 15.3306 | 0.7571906 | 280.0855 | 0 |
10. Estimación GARCH(1,1) si procede
# ----- Parámetros de control -----
prefer_garch <- TRUE # si TRUE, preferir modelos GARCH cuando convergen y cumplen diagnósticos
require_lb2_non_sig <- TRUE # si TRUE, exigir LB(resid.^2) p > 0.05 para aceptar un GARCH candidato
garch_order <- c(1,1)
dist_garch <- "std" # "norm" o "std"
outfile_prefix <- "TNX" # prefijo para archivos (TNX)
# ----- 10: Estimar modelos GARCH si procede -----
need_garch <- any(diag_all$ARCH_p < 0.05, na.rm = TRUE) || prefer_garch
fit_ar_garch <- fit_ma_garch <- fit_arma_garch <- NULL
conv_ar_garch <- conv_ma_garch <- conv_arma_garch <- NA
if (need_garch) {
spec_ar <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_ar_p, 0), include.mean = include.mean),
distribution.model = dist_garch)
spec_ma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(0, best_ma_q), include.mean = include.mean),
distribution.model = dist_garch)
spec_arma <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = garch_order),
mean.model = list(armaOrder = c(best_arma_p, best_arma_q), include.mean = include.mean),
distribution.model = dist_garch)
fit_ar_garch <- tryCatch(ugarchfit(spec_ar, data = serie, solver = "hybrid", silent = TRUE),
error = function(e) { message("AR-GARCH error: ", e$message); NULL })
fit_ma_garch <- tryCatch(ugarchfit(spec_ma, data = serie, solver = "hybrid", silent = TRUE),
error = function(e) { message("MA-GARCH error: ", e$message); NULL })
fit_arma_garch <- tryCatch(ugarchfit(spec_arma, data = serie, solver = "hybrid", silent = TRUE),
error = function(e) { message("ARMA-GARCH error: ", e$message); NULL })
conv_ar_garch <- if (!is.null(fit_ar_garch)) tryCatch(fit_ar_garch@fit$convergence, error = function(e) NA) else NA
conv_ma_garch <- if (!is.null(fit_ma_garch)) tryCatch(fit_ma_garch@fit$convergence, error = function(e) NA) else NA
conv_arma_garch <- if (!is.null(fit_arma_garch)) tryCatch(fit_arma_garch@fit$convergence, error = function(e) NA) else NA
}
conv_df <- data.frame(Modelo = c("AR-GARCH","MA-GARCH","ARMA-GARCH"),
Converge = c(conv_ar_garch, conv_ma_garch, conv_arma_garch),
stringsAsFactors = FALSE)
write.csv(conv_df, file.path(res_folder, paste0(outfile_prefix, "_garch_conv.csv")), row.names = FALSE)
print(conv_df) # o kable(conv_df) ## Modelo Converge
## 1 AR-GARCH 0
## 2 MA-GARCH 0
## 3 ARMA-GARCH 011. Selección por AIC entre modelos
comparacion <- data.frame(Model = character(), AIC = numeric(), stringsAsFactors = FALSE)
comparacion <- rbind(comparacion, data.frame(Model = paste0("AR(",best_ar_p,")"), AIC = AIC(model_ar)))
comparacion <- rbind(comparacion, data.frame(Model = paste0("MA(",best_ma_q,")"), AIC = AIC(model_ma)))
comparacion <- rbind(comparacion, data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")"), AIC = AIC(model_arma)))
if (!is.null(fit_ar_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("AR(",best_ar_p,")-GARCH(1,1)"), AIC = infocriteria(fit_ar_garch)["Akaike",1]))
if (!is.null(fit_ma_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("MA(",best_ma_q,")-GARCH(1,1)"), AIC = infocriteria(fit_ma_garch)["Akaike",1]))
if (!is.null(fit_arma_garch)) comparacion <- rbind(comparacion, data.frame(Model = paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)"), AIC = infocriteria(fit_arma_garch)["Akaike",1]))
comparacion <- comparacion %>% arrange(AIC)
write.csv(comparacion, file.path(res_folder, "TNX_comparacion_aic.csv"), row.names = FALSE)
kable(comparacion, caption = "Comparación por AIC (menor = mejor)")| Model | AIC |
|---|---|
| AR(2) | -5803.033264 |
| MA(2) | -5802.976588 |
| ARMA(2,1) | -5801.730009 |
| MA(2)-GARCH(1,1) | -3.077811 |
| AR(2)-GARCH(1,1) | -3.077770 |
| ARMA(2,1)-GARCH(1,1) | -3.077766 |
## Mejor modelo por AIC: AR(2)12. Diagnóstico final del mejor modelo
# construir lista de objetos
objs <- list()
if (exists("model_ar")) objs[[paste0("AR(",best_ar_p,")")]] <- model_ar
if (exists("model_ma")) objs[[paste0("MA(",best_ma_q,")")]] <- model_ma
if (exists("model_arma")) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")")]] <- model_arma
if (!is.null(fit_ar_garch)) objs[[paste0("AR(",best_ar_p,")-GARCH(1,1)")]] <- fit_ar_garch
if (!is.null(fit_ma_garch)) objs[[paste0("MA(",best_ma_q,")-GARCH(1,1)")]] <- fit_ma_garch
if (!is.null(fit_arma_garch)) objs[[paste0("ARMA(",best_arma_p,",",best_arma_q,")-GARCH(1,1)")]] <- fit_arma_garch
if (!mejor_modelo %in% names(objs)) stop("No se encontró objeto para el modelo seleccionado: ", mejor_modelo)
modelo_final <- objs[[mejor_modelo]]
# residuos (estandarizados si GARCH)
if (inherits(modelo_final, "uGARCHfit")) {
residuales <- tryCatch(residuals(modelo_final, standardize = TRUE), error = function(e) residuals(modelo_final))
} else {
residuales <- residuals(modelo_final)
}
residuales_cuad <- residuales^2
lb_res <- tryCatch(Box.test(as.numeric(residuales), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
lb_res2 <- tryCatch(Box.test(as.numeric(residuales_cuad), lag = lag_test, type = "Ljung-Box"), error = function(e) list(statistic = NA, p.value = NA))
arch_res <- tryCatch(FinTS::ArchTest(as.numeric(residuales), lags = lag_test), error = function(e) list(statistic = NA, p.value = NA))
diag_final <- data.frame(Test = c("LjungBox resid.", "LjungBox resid.^2", "ARCH-LM"),
Stat = c(lb_res$statistic, lb_res2$statistic, arch_res$statistic),
P_value = c(lb_res$p.value, lb_res2$p.value, arch_res$p.value),
stringsAsFactors = FALSE)
write.csv(diag_final, file.path(res_folder, paste0(outfile_prefix, "_diag_final.csv")), row.names = FALSE)
print(diag_final) # o kable(diag_final)## Test Stat P_value
## 1 LjungBox resid. 16.01847 7.154777e-01
## 2 LjungBox resid.^2 763.50702 0.000000e+00
## 3 ARCH-LM 279.81267 1.049243e-47# Graficar ACF (residuos y residuos^2)
par(mfrow = c(2,1))
plot_acf_eviews(residuales, main = paste("ACF residuales -", mejor_modelo), col = "#2c3e50")
plot_acf_eviews(residuales_cuad, main = paste("ACF residuales^2 -", mejor_modelo), col = col_ret2)
13-14. Volatilidad condicional y VaR (en unidades de yield y conversión a pérdidas monetarias)
if (inherits(modelo_final, "uGARCHfit")) {
vol_cond <- sigma(modelo_final)
if (is.null(vol_cond)) stop("No se pudo extraer sigma() del modelo GARCH seleccionado.")
ultima_sigma <- tail(as.numeric(vol_cond), 1)
dist_model <- modelo_final@model$modeldesc$distribution
shape_param <- if (dist_model == "std" && "shape" %in% names(coef(modelo_final))) coef(modelo_final)["shape"] else NULL
calcular_var_model <- function(sigma_val, nivel_conf) {
alpha <- 1 - nivel_conf
if (!is.null(shape_param) && dist_model == "std") q <- qt(alpha, df = shape_param) else q <- qnorm(alpha)
sigma_val * q
}
var_model <- sapply(niveles_conf, function(p) calcular_var_model(ultima_sigma, p))
} else {
sd_hist <- sd(serie, na.rm = TRUE)
ultima_sigma <- sd_hist
var_model <- sd_hist * qnorm(1 - niveles_conf)
}
var_model_df <- data.frame(
Nivel = paste0(niveles_conf*100, "%"),
VaR_yield = var_model, # VaR in yield units (percentage points)
VaR_yield_abs = abs(var_model)
)
# ----- RiskMetrics (EWMA) sobre ΔYield -----
rend_ts <- serie; n <- length(rend_ts)
var_rm <- numeric(n); var_rm[1] <- var(rend_ts, na.rm = TRUE)
for (t in 2:n) var_rm[t] <- lambda_rm * var_rm[t-1] + (1 - lambda_rm) * (rend_ts[t-1]^2)
sigma_rm <- sqrt(var_rm); ultima_sigma_rm <- tail(sigma_rm,1)
var_rm_vals <- sapply(niveles_conf, function(p) ultima_sigma_rm * qnorm(1 - p))
var_rm_df <- data.frame(
Nivel = paste0(niveles_conf*100, "%"),
VaR_yield_RM = var_rm_vals,
VaR_yield_RM_abs = abs(var_rm_vals)
)
# ----- Conversión a pérdidas monetarias usando Duration_mod -----
to_monetary_loss <- function(dy_pp) {
# dy_pp: change in yield in same units as serie (percent). Convert to decimal fraction:
dy_dec <- dy_pp / 100
loss <- - Duration_mod * dy_dec * investment_ref
return(loss)
}
var_model_df$Perdida_USD <- round(to_monetary_loss(var_model_df$VaR_yield), 2)
var_rm_df$Perdida_USD <- round(to_monetary_loss(var_rm_df$VaR_yield_RM), 2)
write.csv(var_model_df, file.path(res_folder, paste0(outfile_prefix, "_var_modelo.csv")), row.names = FALSE)
write.csv(var_rm_df, file.path(res_folder, paste0(outfile_prefix, "_var_riskmetrics.csv")), row.names = FALSE)
# Mostrar tablas
print(var_model_df) # o kable(var_model_df)## Nivel VaR_yield VaR_yield_abs Perdida_USD
## 1 90% -0.07309888 0.07309888 584.79
## 2 95% -0.09382140 0.09382140 750.57
## 3 99% -0.13269340 0.13269340 1061.55## Nivel VaR_yield_RM VaR_yield_RM_abs Perdida_USD
## 1 90% -0.04742579 0.04742579 379.41
## 2 95% -0.06087034 0.06087034 486.96
## 3 99% -0.08609009 0.08609009 688.72comparativa <- data.frame(
Nivel = var_model_df$Nivel,
VaR_Model_yield_pp = var_model_df$VaR_yield,
Perdida_Model_USD = var_model_df$Perdida_USD,
VaR_RM_yield_pp = var_rm_df$VaR_yield_RM,
Perdida_RM_USD = var_rm_df$Perdida_USD,
stringsAsFactors = FALSE
)
write.csv(comparativa, file.path(res_folder, paste0(outfile_prefix, "_comparativa_var.csv")), row.names = FALSE)
print(comparativa) # o kable(comparativa)## Nivel VaR_Model_yield_pp Perdida_Model_USD VaR_RM_yield_pp Perdida_RM_USD
## 1 90% -0.07309888 584.79 -0.04742579 379.41
## 2 95% -0.09382140 750.57 -0.06087034 486.96
## 3 99% -0.13269340 1061.55 -0.08609009 688.72Resumen del trabajo
Tabla Resumen del trabajo
## --------- Resumen del trabajo (Mejor modelo debe ser GARCH si es posible) ----------
library(quantmod); library(tidyquant); library(forecast); library(tseries)
library(FinTS); library(rugarch); library(dplyr); library(tibble); library(knitr)
# Parámetros generales (ajusta si hace falta)
start_date <- as.Date("2017-11-04")
end_date <- as.Date("2025-11-04")
lambda_rm <- 0.94
niveles_conf <- c(0.90, 0.95, 0.99)
maxP <- 5; maxQ <- 5
include.mean <- TRUE
lag_test <- 20
notional <- 100000 # monto de referencia para pérdidas
duration_tnx <- 8 # duración modificada para el 10y
activos <- tribble(
~Activo, ~Ticker, ~Tipo,
"NVIDIA (NVDA)", "NVDA", "equity",
"Vanguard Total Market", "VTI", "equity",
"USD/MXN spot", "USDMXN=X", "fx",
"UST 10y yield (^TNX)", "^TNX", "yield"
)
# Funciones auxiliares (safe_get, búsqueda ARMA, diagnóstico, ajuste GARCH y VaR)
safe_get <- function(ticker, from, to) {
# Intenta quantmod y fallback tidyquant; devuelve xts (adjusted price) o error
tryCatch({
suppressWarnings(getSymbols(ticker, src = "yahoo", from = from, to = to, auto.assign = FALSE))
}, error = function(e) {
tq <- tryCatch(tq_get(ticker, get = "stock.prices", from = from, to = to), error = function(e) NULL)
if (is.null(tq) || nrow(tq) == 0) stop("No se pudo descargar ", ticker)
xts(tq$adjusted, order.by = as.Date(tq$date))
})
}
select_best_arma <- function(x, maxP = 5, maxQ = 5, include.mean = TRUE) {
best_AIC <- Inf; best_fit <- NULL; best_p <- NA; best_q <- NA
for (p in 0:maxP) for (q in 0:maxQ) {
if (p == 0 && q == 0) next
fit <- tryCatch(arima(x, order = c(p,0,q), include.mean = include.mean, method = "ML"), error = function(e) NULL)
if (!is.null(fit)) {
a <- AIC(fit)
if (a < best_AIC) { best_AIC <- a; best_fit <- fit; best_p <- p; best_q <- q }
}
}
list(fit = best_fit, p = best_p, q = best_q, AIC = best_AIC)
}
diagnosticar_res <- function(obj, is_garch = FALSE, lag_test = 20) {
# Retorna lista con p-values: lb_res, lb_res2, arch_p
out <- list(lb_res_p = NA_real_, lb_res2_p = NA_real_, arch_p = NA_real_)
if (is.null(obj)) return(out)
if (is_garch) {
res_std <- tryCatch(residuals(obj, standardize = TRUE), error = function(e) residuals(obj))
out$lb_res_p <- tryCatch(Box.test(as.numeric(res_std), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA_real_)
out$lb_res2_p <- tryCatch(Box.test(as.numeric(res_std^2), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA_real_)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_std), lags = lag_test)$p.value, error = function(e) NA_real_)
} else {
res_plain <- tryCatch(residuals(obj), error = function(e) NULL)
if (is.null(res_plain)) return(out)
out$lb_res_p <- tryCatch(Box.test(as.numeric(res_plain), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA_real_)
out$lb_res2_p <- tryCatch(Box.test(as.numeric(res_plain^2), lag = lag_test, type = "Ljung-Box")$p.value, error = function(e) NA_real_)
out$arch_p <- tryCatch(FinTS::ArchTest(as.numeric(res_plain), lags = lag_test)$p.value, error = function(e) NA_real_)
}
out
}
try_fit_garch <- function(x, arma_order, include.mean = TRUE, dist = "std") {
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = arma_order, include.mean = include.mean),
distribution.model = dist)
fit <- tryCatch(ugarchfit(spec, data = x, solver = "hybrid", silent = TRUE), error = function(e) { message("GARCH error: ", e$message); NULL })
fit
}
calc_var_values <- function(x, fit_garch, niveles_conf, lambda_rm) {
n <- length(x)
# var model
if (!is.null(fit_garch) && inherits(fit_garch, "uGARCHfit")) {
svec <- tryCatch(sigma(fit_garch), error = function(e) NULL)
if (is.null(svec)) {
sigma_last <- sd(x, na.rm = TRUE)
var_model <- sapply(niveles_conf, function(p) sigma_last * qnorm(1 - p))
} else {
sigma_last <- tail(as.numeric(svec), 1)
dist_model <- fit_garch@model$modeldesc$distribution
coefs <- coef(fit_garch)
shape <- if ("shape" %in% names(coefs)) coefs["shape"] else NA
qfun <- function(alpha) {
if (!is.na(shape) && dist_model == "std") qt(alpha, df = shape) else qnorm(alpha)
}
var_model <- sapply(niveles_conf, function(p) sigma_last * qfun(1 - p))
}
} else {
sigma_last <- sd(x, na.rm = TRUE)
var_model <- sapply(niveles_conf, function(p) sigma_last * qnorm(1 - p))
}
# var riskmetrics
var_ewma <- numeric(n); var_ewma[1] <- var(x, na.rm = TRUE)
if (n >= 2) for (t in 2:n) var_ewma[t] <- lambda_rm * var_ewma[t-1] + (1 - lambda_rm) * x[t-1]^2
sigma_rm_last <- sqrt(tail(var_ewma, 1))
var_rm <- sapply(niveles_conf, function(p) sigma_rm_last * qnorm(1 - p))
list(var_model = var_model, var_rm = var_rm, sigma_last = sigma_last, sigma_rm_last = sigma_rm_last)
}
# Loop principal: para cada activo seleccionamos ARMA, estimamos 3 GARCH y elegimos mejor GARCH
resumen_list <- vector("list", nrow(activos))
for (i in seq_len(nrow(activos))) {
act <- activos$Activo[i]; tick <- activos$Ticker[i]; tipo <- activos$Tipo[i]
cat("Procesando:", act, "(", tick, ")\n")
px <- tryCatch(safe_get(tick, start_date, end_date), error = function(e) { cat("Error descarga:", e$message, "\n"); NULL })
if (is.null(px)) next
# construir serie a modelar
if (tipo == "yield") {
serie <- diff(as.numeric(px)) # Δyield (pp)
serie_label <- "Δyield (pp)"
} else {
serie <- as.numeric(dailyReturn(px, type = "log"))
serie_label <- "log-ret"
}
serie <- na.omit(serie); nobs <- length(serie)
if (nobs < 30) { cat("Datos insuficientes para", tick, "\n"); next }
# ADF
adf_level <- tryCatch(adf.test(as.numeric(px)), error = function(e) NULL)
adf_model <- tryCatch(adf.test(serie), error = function(e) NULL)
# Selección ARMA
sel <- select_best_arma(serie, maxP = maxP, maxQ = maxQ, include.mean = include.mean)
arma_fit <- sel$fit; arma_ord <- c(sel$p, sel$q); arma_name <- paste0("ARMA(", sel$p, ",", sel$q, ")")
# Diagnóstico residuos ARMA
diag_arma <- diagnosticar_res(arma_fit, is_garch = FALSE, lag_test = lag_test)
arch_pval <- diag_arma$arch_p
# Estimar los 3 GARCH: AR-GARCH, MA-GARCH, ARMA-GARCH (siempre intentamos, no solo si ARCH significativo)
fit_ar_garch <- try_fit_garch(serie, arma_order = c(sel$p, 0), include.mean = include.mean)
fit_ma_garch <- try_fit_garch(serie, arma_order = c(0, sel$q), include.mean = include.mean)
fit_arma_garch <- try_fit_garch(serie, arma_order = c(sel$p, sel$q), include.mean = include.mean)
# Recopilar info de convergencia/AIC/diagnósticos
entries <- list(
list(label = paste0("AR(", sel$p, ")-GARCH(1,1)"), fit = fit_ar_garch),
list(label = paste0("MA(", sel$q, ")-GARCH(1,1)"), fit = fit_ma_garch),
list(label = paste0("ARMA(", sel$p, ",", sel$q, ")-GARCH(1,1)"), fit = fit_arma_garch)
)
garch_info <- lapply(entries, function(en) {
fit <- en$fit
if (is.null(fit)) return(data.frame(Model = en$label, Converge = FALSE, AIC = NA_real_, LB2_p = NA_real_, ARCH_p = NA_real_, stringsAsFactors = FALSE))
converge <- tryCatch(fit@fit$convergence == 0, error = function(e) NA)
aicv <- tryCatch(infocriteria(fit)["Akaike",1], error = function(e) NA_real_)
diagv <- diagnosticar_res(fit, is_garch = TRUE, lag_test = lag_test)
data.frame(Model = en$label, Converge = isTRUE(converge), AIC = aicv, LB2_p = diagv$lb_res2_p, ARCH_p = diagv$arch_p, stringsAsFactors = FALSE)
})
garch_tab <- bind_rows(garch_info)
# Selección del mejor modelo: prioridad 1) Mejor GARCH que converge y pasa diagnósticos (LB2_p>0.05 & ARCH_p>0.05)
# 2) Si ninguno, mejor GARCH convergido por AIC. 3) Si ninguno convergió, seleccionar ARMA (y marcar GARCH NA).
best_garch <- NULL
garch_valid <- garch_tab %>% filter(Converge == TRUE & !is.na(LB2_p) & !is.na(ARCH_p) & LB2_p > 0.05 & ARCH_p > 0.05)
if (nrow(garch_valid) > 0) {
best_garch <- garch_valid %>% arrange(AIC) %>% slice(1)
} else {
# intentar cualquier GARCH convergido
garch_conv <- garch_tab %>% filter(Converge == TRUE & !is.na(AIC))
if (nrow(garch_conv) > 0) best_garch <- garch_conv %>% arrange(AIC) %>% slice(1)
}
if (!is.null(best_garch)) {
mejor_modelo <- best_garch$Model[1]
# recuperar objeto fit correspondiente (uso if/else seguro en lugar de switch)
if (mejor_modelo == paste0("AR(", sel$p, ")-GARCH(1,1)")) {
modelo_final <- fit_ar_garch
} else if (mejor_modelo == paste0("MA(", sel$q, ")-GARCH(1,1)")) {
modelo_final <- fit_ma_garch
} else if (mejor_modelo == paste0("ARMA(", sel$p, ",", sel$q, ")-GARCH(1,1)")) {
modelo_final <- fit_arma_garch
} else {
modelo_final <- NULL
}
garch_selected_flag <- TRUE
# obtener parámetros
coefs <- tryCatch(coef(modelo_final), error = function(e) NULL)
omega <- if (!is.null(coefs) && "omega" %in% names(coefs)) coefs["omega"] else NA_real_
alpha1 <- if (!is.null(coefs) && any(grepl("^alpha", names(coefs)))) coefs[grep("^alpha", names(coefs))[1]] else NA_real_
beta1 <- if (!is.null(coefs) && any(grepl("^beta", names(coefs)))) coefs[grep("^beta", names(coefs))[1]] else NA_real_
alpha_beta <- as.numeric(alpha1) + as.numeric(beta1)
# diagnostics for selected
diag_sel <- diagnosticar_res(modelo_final, is_garch = TRUE, lag_test = lag_test)
} else {
# no hubo GARCH convergido -> fallback a ARMA
modelo_final <- arma_fit
mejor_modelo <- arma_name
garch_selected_flag <- FALSE
omega <- alpha1 <- beta1 <- alpha_beta <- NA_real_
diag_sel <- diagnosticar_res(modelo_final, is_garch = FALSE, lag_test = lag_test)
}
# Calcular VaR y pérdidas
varvals <- calc_var_values(serie, if (garch_selected_flag) modelo_final else NULL, niveles_conf, lambda_rm)
var_model <- varvals$var_model; var_rm <- varvals$var_rm
if (tipo == "yield") {
loss_model <- - (var_model / 100) * duration_tnx * notional
loss_rm <- - (var_rm / 100) * duration_tnx * notional
} else {
loss_model <- - var_model * notional
loss_rm <- - var_rm * notional
}
resumen_list[[i]] <- tibble(
Activo = act, Ticker = tick, Tipo = tipo, Serie_modelada = serie_label,
ADF_p_level = if (!is.null(adf_level)) round(adf_level$p.value, 6) else NA_real_,
ADF_p_serie = if (!is.null(adf_model)) round(adf_model$p.value, 6) else NA_real_,
Mejor_modelo = mejor_modelo,
Es_GARCH = garch_selected_flag,
AIC_mejor = if (!is.null(best_garch)) round(best_garch$AIC, 4) else round(sel$AIC, 4),
ARCH_p_ARMA = round(arch_pval, 6),
Seleccion_GARCH_table = list(garch_tab),
omega = round(as.numeric(omega), 6),
alpha1 = round(as.numeric(alpha1), 6),
beta1 = round(as.numeric(beta1), 6),
alpha_plus_beta = round(as.numeric(alpha_beta), 6),
LB2_p_selected = round(diag_sel$lb_res2_p, 6),
ARCH_p_selected = round(diag_sel$arch_p, 6),
sigma_last = round(varvals$sigma_last, 6),
VaR_model_90 = round(var_model[1], 6), VaR_model_95 = round(var_model[2], 6), VaR_model_99 = round(var_model[3], 6),
VaR_RM_90 = round(var_rm[1], 6), VaR_RM_95 = round(var_rm[2], 6), VaR_RM_99 = round(var_rm[3], 6),
Loss_model_90 = round(loss_model[1], 2), Loss_model_95 = round(loss_model[2], 2), Loss_model_99 = round(loss_model[3], 2),
Loss_RM_90 = round(loss_rm[1], 2), Loss_RM_95 = round(loss_rm[2], 2), Loss_RM_99 = round(loss_rm[3], 2)
)
} ## Procesando: NVIDIA (NVDA) ( NVDA )## Procesando: Vanguard Total Market ( VTI )## Procesando: USD/MXN spot ( USDMXN=X )## Procesando: UST 10y yield (^TNX) ( ^TNX )tabla_resumen <- bind_rows(resumen_list)
# Mostrar tabla (resumen) con kable
kable(tabla_resumen, caption = "Resumen comparativo (Mejor modelo prioriza GARCH)", digits = 4) | Activo | Ticker | Tipo | Serie_modelada | ADF_p_level | ADF_p_serie | Mejor_modelo | Es_GARCH | AIC_mejor | ARCH_p_ARMA | Seleccion_GARCH_table | omega | alpha1 | beta1 | alpha_plus_beta | LB2_p_selected | ARCH_p_selected | sigma_last | VaR_model_90 | VaR_model_95 | VaR_model_99 | VaR_RM_90 | VaR_RM_95 | VaR_RM_99 | Loss_model_90 | Loss_model_95 | Loss_model_99 | Loss_RM_90 | Loss_RM_95 | Loss_RM_99 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NVIDIA (NVDA) | NVDA | equity | log-ret | 0.01 | 0.01 | AR(4)-GARCH(1,1) | TRUE | -4.2388 | 0 | AR(4)-GARCH(1,1) , MA(5)-GARCH(1,1) , ARMA(4,5)-GARCH(1,1), TRUE , TRUE , TRUE , -4.23878709379274 , -4.23788233160719 , -4.23782640483954 , 0.998238200587823 , 0.998189204352688 , 0.997550035645691 , 0.998115484385657 , 0.998086380137436 , 0.997475406159233 | 0 | 0.0955 | 0.8790 | 0.9745 | 0.9982 | 0.9981 | 0.0259 | -0.0377 | -0.0512 | -0.0840 | -0.0288 | -0.0370 | -0.0524 | 3772.15 | 5118.82 | 8403.86 | 2883.94 | 3701.50 | 5235.10 |
| Vanguard Total Market | VTI | equity | log-ret | 0.01 | 0.01 | ARMA(4,4)-GARCH(1,1) | TRUE | -6.4703 | 0 | AR(4)-GARCH(1,1) , MA(4)-GARCH(1,1) , ARMA(4,4)-GARCH(1,1), TRUE , TRUE , TRUE , -6.46515706890259 , -6.46555369127165 , -6.4703424249007 , 0.799698327552294 , 0.799563861239973 , 0.896060388823766 , 0.792446700277363 , 0.792417768274921 , 0.892546947085759 | 0 | 0.1704 | 0.8167 | 0.9871 | 0.8961 | 0.8925 | 0.0082 | -0.0118 | -0.0159 | -0.0258 | -0.0099 | -0.0127 | -0.0180 | 1176.67 | 1589.66 | 2577.90 | 990.48 | 1271.27 | 1797.98 |
| USD/MXN spot | USDMXN=X | fx | log-ret | NA | 0.01 | ARMA(2,3)-GARCH(1,1) | TRUE | -7.1604 | 0 | AR(2)-GARCH(1,1) , MA(3)-GARCH(1,1) , ARMA(2,3)-GARCH(1,1), TRUE , TRUE , TRUE , -7.15708155824615 , -7.15627774245626 , -7.16038750977678 , 0.81434949883795 , 0.828026724071456 , 0.867090193182517 , 0.82469769875898 , 0.83518910234124 , 0.874334185149089 | 0 | 0.1770 | 0.7731 | 0.9502 | 0.8671 | 0.8743 | 0.0045 | -0.0065 | -0.0087 | -0.0141 | -0.0038 | -0.0049 | -0.0069 | 647.79 | 873.41 | 1408.65 | 381.04 | 489.06 | 691.68 |
| UST 10y yield (^TNX) | ^TNX | yield | Δyield (pp) | NA | 0.01 | MA(4)-GARCH(1,1) | TRUE | -6.4478 | 1 | AR(4)-GARCH(1,1) , MA(4)-GARCH(1,1) , ARMA(4,4)-GARCH(1,1), FALSE , TRUE , TRUE , NA , -6.44783277666115 , -3.01448462360774 , NA , 1 , 0 , NA , 1 , 0 | 0 | 0.3292 | 0.6560 | 0.9852 | 1.0000 | 1.0000 | 0.0342 | -0.0495 | -0.0670 | -0.1093 | -0.0474 | -0.0608 | -0.0860 | 396.01 | 536.12 | 874.44 | 379.02 | 486.47 | 688.02 |
Grà ficas Resumen del trabajo
# Chunk completo: descargar, modelar y graficar volatilidades GARCH vs RiskMetrics
# Requiere internet. Ajusta start_date / res_root / lambda_rm si quieres.
# Paquetes necesarios (instala si falta)
needed <- c("quantmod","tidyquant","xts","zoo","dplyr","ggplot2","tidyr","rugarch","FinTS","tseries")
inst <- needed[!(needed %in% installed.packages()[,"Package"])]
if (length(inst) > 0) {
message("Instalando paquetes faltantes: ", paste(inst, collapse = ", "))
install.packages(inst, dependencies = TRUE)
}
invisible(lapply(needed, require, character.only = TRUE))
# CONFIG
assets <- list(
NVDA = list(ticker = "NVDA", type = "equity"),
VTI = list(ticker = "VTI", type = "equity"),
USDMXN = list(ticker = "USDMXN=X", type = "fx"),
TNX = list(ticker = "^TNX", type = "yield")
)
start_date <- as.Date("2017-11-04")
end_date <- Sys.Date()
res_root <- "results_all_assets_recalc"
if (!dir.exists(res_root)) dir.create(res_root, recursive = TRUE)
lambda_rm <- 0.94
# Helpers
safe_get <- function(ticker, from, to) {
tryCatch({
suppressWarnings(quantmod::getSymbols(ticker, src = "yahoo", from = from, to = to, auto.assign = FALSE))
}, error = function(e) {
# fallback tidyquant
tq <- tryCatch(tidyquant::tq_get(ticker, from = from, to = to), error = function(e) NULL)
if (is.null(tq) || nrow(tq) == 0) stop(paste("No se pudo descargar", ticker))
xts(tq$adjusted, order.by = as.Date(tq$date))
})
}
calc_ewma_sigma <- function(x, lambda = 0.94) {
x <- as.numeric(na.omit(x)); n <- length(x)
if (n == 0) return(numeric(0))
var_rm <- numeric(n)
var_rm[1] <- var(x, na.rm = TRUE)
if (is.na(var_rm[1]) || var_rm[1] == 0) var_rm[1] <- mean(x^2, na.rm = TRUE)
if (n >= 2) {
for (t in 2:n) var_rm[t] <- lambda * var_rm[t-1] + (1 - lambda) * (x[t-1]^2)
}
sqrt(var_rm)
}
# Procesamiento por activo
vol_plots <- list()
for (asset in names(assets)) {
info <- assets[[asset]]
folder <- file.path(res_root, asset)
if (!dir.exists(folder)) dir.create(folder, recursive = TRUE)
cat("\n--- Procesando:", asset, "(", info$ticker, ") ---\n")
# Descargar
dat_raw <- tryCatch(safe_get(info$ticker, start_date, end_date), error = function(e) { cat("ERROR descarga:", e$message, "\n"); NULL })
if (is.null(dat_raw)) next
if (!inherits(dat_raw, "xts")) {
cat(" Descarga no es xts; intentando coerción...\n")
dat_raw <- tryCatch(as.xts(dat_raw), error = function(e) NULL)
if (is.null(dat_raw)) { cat(" No se pudo convertir descarga a xts. Saltando.\n"); next }
}
# Seleccionar serie de precios/valor relevante (Adjusted si existe)
datp <- dat_raw
if (NCOL(dat_raw) > 1) {
adj_idx <- grep("adjusted|Adjusted|Adj", colnames(dat_raw))
if (length(adj_idx) >= 1) datp <- dat_raw[, adj_idx[1]] else datp <- dat_raw[,1]
}
colnames(datp) <- asset
# Construir serie de análisis: log-returns para equity/fx, diff para yield
serie_xts <- NULL
if (info$type %in% c("equity", "fx")) {
prices <- as.numeric(datp)
dates_prices <- index(datp)
ok <- which(!is.na(prices))
if (length(ok) < 2) { cat(" Datos insuficientes para retornos. Saltando.\n"); next }
prices <- prices[ok]; dates_prices <- dates_prices[ok]
logp <- log(prices)
ret <- diff(logp)
dates_ret <- dates_prices[-1]
# Alinear longitudes defensivamente
minlen <- min(length(ret), length(dates_ret))
ret <- ret[seq_len(minlen)]; dates_ret <- dates_ret[seq_len(minlen)]
# eliminar NAs
good <- which(!is.na(ret))
if (length(good) == 0) { cat(" Todos los retornos son NA. Saltando.\n"); next }
ret <- ret[good]; dates_ret <- as.Date(dates_ret[good])
serie_xts <- tryCatch(xts(as.numeric(ret), order.by = dates_ret), error = function(e) NULL)
if (is.null(serie_xts)) { cat(" No se pudo construir xts para retornos. Saltando.\n"); next }
} else {
# yield: diferencias simples
y <- as.numeric(datp); dates_y <- index(datp)
ok <- which(!is.na(y))
if (length(ok) < 2) { cat(" Datos insuficientes para yield. Saltando.\n"); next }
y <- y[ok]; dates_y <- dates_y[ok]
dy <- diff(y); dates_dy <- dates_y[-1]
minlen <- min(length(dy), length(dates_dy))
dy <- dy[seq_len(minlen)]; dates_dy <- dates_dy[seq_len(minlen)]
good <- which(!is.na(dy))
if (length(good) == 0) { cat(" Todas las diferencias de yield son NA. Saltando.\n"); next }
dy <- dy[good]; dates_dy <- as.Date(dates_dy[good])
serie_xts <- tryCatch(xts(as.numeric(dy), order.by = dates_dy), error = function(e) NULL)
if (is.null(serie_xts)) { cat(" No se pudo construir xts para Δyield. Saltando.\n"); next }
}
cat(" Serie construida. Observaciones:", NROW(serie_xts), "\n")
# Guardar serie
write.csv(data.frame(Date = as.character(index(serie_xts)), serie = coredata(serie_xts)),
file = file.path(folder, paste0(asset, "_serie.csv")), row.names = FALSE)
# Selección ARMA por AIC (0..5)
bestAIC <- Inf; best_fit <- NULL; best_pq <- c(NA, NA)
xdata <- coredata(serie_xts)
for (p in 0:5) for (q in 0:5) {
if (p == 0 && q == 0) next
fit_try <- tryCatch(arima(xdata, order = c(p,0,q), include.mean = TRUE, method = "ML"), error = function(e) NULL)
if (!is.null(fit_try)) {
aicv <- AIC(fit_try)
if (aicv < bestAIC) { bestAIC <- aicv; best_fit <- fit_try; best_pq <- c(p,q) }
}
}
cat(" Mejor ARMA:", paste0("ARMA(", best_pq[1], ",", best_pq[2], ")"), "AIC=", round(bestAIC, 2), "\n")
# Test ARCH en residuales
arch_test <- NULL; arch_p <- NA
if (!is.null(best_fit)) {
resid <- residuals(best_fit)
arch_test <- tryCatch(FinTS::ArchTest(resid, lags = 20), error = function(e) NULL)
if (!is.null(arch_test)) arch_p <- arch_test$p.value
}
cat(" ARCH p-value:", arch_p, "\n")
# Ajuste GARCH(1,1) si ARCH significativo
fit_garch <- NULL; sigma_garch_xts <- NULL; omega <- alpha1 <- beta1 <- NA
if (!is.null(arch_test) && arch_test$p.value < 0.05) {
cat(" Ajustando GARCH(1,1) con mean ARMA:", paste0("(", best_pq[1], ",", best_pq[2], ")"), "\n")
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1)),
mean.model = list(armaOrder = best_pq, include.mean = TRUE),
distribution.model = "std")
fitg <- tryCatch(ugarchfit(spec, data = xdata, solver = "hybrid", silent = TRUE), error = function(e) { cat(" Error GARCH:", e$message, "\n"); NULL })
if (!is.null(fitg) && inherits(fitg, "uGARCHfit")) {
fit_garch <- fitg
cof <- tryCatch(coef(fitg), error = function(e) NULL)
if (!is.null(cof)) {
omega <- ifelse("omega" %in% names(cof), as.numeric(cof["omega"]), NA)
an <- names(cof)[grepl("alpha", names(cof), ignore.case = TRUE)]
bn <- names(cof)[grepl("beta", names(cof), ignore.case = TRUE)]
if (length(an) > 0) alpha1 <- as.numeric(cof[an[1]])
if (length(bn) > 0) beta1 <- as.numeric(cof[bn[1]])
}
svec <- tryCatch(sigma(fitg), error = function(e) NULL)
if (!is.null(svec)) {
# alinear defensivamente con serie_xts
if (length(svec) == NROW(serie_xts)) {
sigma_garch_xts <- xts(as.numeric(svec), order.by = index(serie_xts))
} else if (length(svec) < NROW(serie_xts)) {
sigma_garch_xts <- xts(as.numeric(svec), order.by = tail(index(serie_xts), length(svec)))
cat(" Nota: sigma GARCH más corta que la serie; alineada por cola\n")
} else {
sigma_garch_xts <- xts(as.numeric(svec[seq_len(NROW(serie_xts))]), order.by = index(serie_xts))
cat(" Nota: sigma GARCH más larga que la serie; truncada\n")
}
} else {
cat(" No se extrajo sigma desde fit_garch\n")
}
} else {
cat(" GARCH no convergió o falló\n")
}
} else {
cat(" No se ajustó GARCH (ARCH no significativo o error)\n")
}
# Calcular sigma RiskMetrics (EWMA)
sigma_rm_vec <- calc_ewma_sigma(coredata(serie_xts), lambda = lambda_rm)
# Alinear longitudes defensivamente
if (length(sigma_rm_vec) != NROW(serie_xts)) {
minlen <- min(length(sigma_rm_vec), NROW(serie_xts))
sigma_rm_xts <- xts(tail(sigma_rm_vec, minlen), order.by = tail(index(serie_xts), minlen))
if (!is.null(sigma_garch_xts)) sigma_garch_xts <- last(sigma_garch_xts, minlen)
cat(" Ajustada longitud RM a", minlen, "\n")
} else {
sigma_rm_xts <- xts(sigma_rm_vec, order.by = index(serie_xts))
}
# Si no hay sigma_garch crear serie NA alineada a RM
if (is.null(sigma_garch_xts)) sigma_garch_xts <- xts(rep(NA_real_, NROW(sigma_rm_xts)), order.by = index(sigma_rm_xts))
# Alinear por intersección o cola
common_dates <- intersect(index(sigma_garch_xts), index(sigma_rm_xts))
if (length(common_dates) >= 10) {
sigma_garch_xts <- sigma_garch_xts[common_dates]; sigma_rm_xts <- sigma_rm_xts[common_dates]
} else {
nmin <- min(NROW(sigma_garch_xts), NROW(sigma_rm_xts))
sigma_garch_xts <- last(sigma_garch_xts, nmin); sigma_rm_xts <- last(sigma_rm_xts, nmin)
}
# Guardar resultados: CSV + RDS (cada activo)
df_plot <- data.frame(Date = as.Date(index(sigma_rm_xts)),
Sigma_RM = as.numeric(sigma_rm_xts),
Sigma_GARCH = as.numeric(sigma_garch_xts))
write.csv(df_plot, file = file.path(folder, paste0(asset, "_vols.csv")), row.names = FALSE)
results_obj <- list(
serie = serie_xts,
best_fit = best_fit,
arch_test = arch_test,
fit_garch = fit_garch,
sigma_garch = sigma_garch_xts,
sigma_rm = sigma_rm_xts,
omega = omega, alpha1 = alpha1, beta1 = beta1
)
saveRDS(results_obj, file = file.path(folder, paste0(asset, "_results.rds")))
# Plot: GARCH vs RM
df_long <- tidyr::pivot_longer(df_plot, cols = c("Sigma_RM", "Sigma_GARCH"),
names_to = "Method", values_to = "Sigma")
p <- ggplot(df_long, aes(x = Date, y = Sigma, color = Method)) +
geom_line(na.rm = TRUE, size = 0.7) +
scale_color_manual(values = c("Sigma_GARCH" = "#D55E00", "Sigma_RM" = "#0072B2"),
labels = c("GARCH sigma_t", "RiskMetrics EWMA")) +
labs(title = paste("Volatility: GARCH vs RiskMetrics —", asset),
subtitle = paste0("lambda_RM = ", lambda_rm),
x = "", y = "Sigma (std dev)") +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
# Guardar PNG local y en carpeta global de comparaciones
png_local <- file.path(folder, paste0(asset, "_vol_comparison.png"))
tryCatch(ggsave(png_local, p, width = 10, height = 4.5, dpi = 300), error = function(e) cat("Error guardando PNG local:", e$message, "\n"))
out_plot_dir <- file.path(res_root, "volatility_comparisons")
if (!dir.exists(out_plot_dir)) dir.create(out_plot_dir, recursive = TRUE)
png_global <- file.path(out_plot_dir, paste0(asset, "_vol_comparison.png"))
tryCatch(ggsave(png_global, p, width = 10, height = 4.5, dpi = 300), error = function(e) cat("Error guardando PNG global:", e$message, "\n"))
cat(" Guardado: ", png_local, " and ", png_global, "\n")
vol_plots[[asset]] <- p
}##
## --- Procesando: NVDA ( NVDA ) ---
## Serie construida. Observaciones: 2018## Mejor ARMA: ARMA(4,4) AIC= -8153.36
## ARCH p-value: 4.367456e-11
## Ajustando GARCH(1,1) con mean ARMA: (4,4)## Guardado: results_all_assets_recalc/NVDA/NVDA_vol_comparison.png and results_all_assets_recalc/volatility_comparisons/NVDA_vol_comparison.png
##
## --- Procesando: VTI ( VTI ) ---
## Serie construida. Observaciones: 2018## Mejor ARMA: ARMA(4,4) AIC= -12082.41
## ARCH p-value: 3.205867e-96
## Ajustando GARCH(1,1) con mean ARMA: (4,4)## Guardado: results_all_assets_recalc/VTI/VTI_vol_comparison.png and results_all_assets_recalc/volatility_comparisons/VTI_vol_comparison.png
##
## --- Procesando: USDMXN ( USDMXN=X ) ---
## Serie construida. Observaciones: 2092## Mejor ARMA: ARMA(4,2) AIC= -14391.66
## ARCH p-value: 1.898094e-57
## Ajustando GARCH(1,1) con mean ARMA: (4,2)## Guardado: results_all_assets_recalc/USDMXN/USDMXN_vol_comparison.png and results_all_assets_recalc/volatility_comparisons/USDMXN_vol_comparison.png
##
## --- Procesando: TNX ( ^TNX ) ---
## Serie construida. Observaciones: 2018## Mejor ARMA: ARMA(2,0) AIC= -5837.08
## ARCH p-value: 4.695068e-48
## Ajustando GARCH(1,1) con mean ARMA: (2,0)## Guardado: results_all_assets_recalc/TNX/TNX_vol_comparison.png and results_all_assets_recalc/volatility_comparisons/TNX_vol_comparison.pngcat("\nProceso completado. Revisa las PNG en:", file.path(res_root, "volatility_comparisons"), "\n")##
## Proceso completado. Revisa las PNG en: results_all_assets_recalc/volatility_comparisons# Imprimir plots en el documento si el motor lo permite
if (length(vol_plots) > 0) {
for (asset in names(vol_plots)) print(vol_plots[[asset]])
}
