Pronóstico del precio del petróleo WTI y decisión de inversión
Trabajo final — Econometría Financiera 2026-2, ITM
Saray Hurtado Aguirre · Davinson Martinez Zapata · Wilson Alberto Orozco Valencia
Noviembre de 2026
FECHA_TOPE <- "2026-09" # último mes que usamos (septiembre de 2026)
PROP_TEST <- 0.20 # 20 % de los datos finales se guarda para el test (el enunciado pide entre 15 % y 25 %)
CAPITAL <- 100000 # capital de la inversión: USD 100.000
UMBRAL <- 0.05 # compramos solo si el precio mediano proyectado a 12 meses sube más de 5 %
H_DEC <- 12 # horizonte de la decisión y de la simulación: 12 meses (se entra en julio 2022; meses evaluados: agosto 2022 a julio 2023)
N_SIM <- 1000 # número de escenarios simulados (Se pide al menos 1.000)
SEMILLA <- 2026 # semilla fija: así la simulación da siempre lo mismo (reproducible)1. Resumen ejecutivo
Este trabajo estudia el precio spot del petróleo West Texas Intermediate (WTI), tomado de FRED en frecuencia mensual entre enero de 2006 y septiembre de 2026, para un total de 249 observaciones (U.S. Energy Information Administration, 2026b). La pregunta es concreta: ubicados en julio de 2022 con USD 100.000, ¿convenía comprar petróleo y mantenerlo doce meses, o dejar el capital en liquidez?
Para responderla se identificaron, estimaron y diagnosticaron diez modelos usando solo los datos hasta julio de 2022. El modelo elegido fue un ARIMA(0,1,1) sobre el logaritmo del precio, sin deriva, con \(\theta = 0{,}3194\). Se escogió porque fue el de menor AICc entre los modelos que cumplieron el diagnóstico y que se podían usar con la información disponible en la fecha de corte.
El modelo proyectó un precio mediano de 97,61 USD a doce meses, 3,95 % por debajo del precio de entrada de 101,62 USD. Como la regla definida exigía una subida mayor al 5 %, la decisión fue mantener liquidez. En los doce meses siguientes el precio cayó a 76,07 USD, de modo que comprar habría significado una pérdida de 25,14 %.
Este resultado debe leerse con cuidado. En el periodo de prueba de 50 meses el modelo superó a la caminata aleatoria por un margen pequeño, sus intervalos de predicción son muy anchos y la decisión habría sido la contraria si se usara la media del pronóstico en lugar de la mediana. Como los residuos mostraron volatilidad cambiante, se estimó además un modelo de la familia ARCH para la varianza, que describe ese comportamiento pero no modifica la decisión. Lo que el trabajo defiende es el procedimiento seguido y no la capacidad del modelo para anticipar el precio del petróleo.
2. Fuente y análisis de mercado
El West Texas Intermediate (WTI) es el crudo de referencia de Estados Unidos. Es un crudo liviano y de bajo contenido de azufre cuyo punto de entrega es Cushing, Oklahoma. No es un puerto: Cushing es un cruce de oleoductos, por eso su precio es sensible a la capacidad de almacenamiento en los tanques y de cuanto cudro se puede mover en los oleoductos (Fattouh, 2011).
Es un commodity y no una acción: no paga dividendos ni representa propiedad sobre una empresa, y su precio responde al balance físico entre oferta y demanda de petróleo.
Su principal mercado de negociación es el contrato de futuros CL de NYMEX, que pertenece a CME Group. Cada contrato equivale a 1.000 barriles, se cotiza en dólares y centavos por barril y se negocia de forma electrónica casi 23 horas al día, de domingo a viernes (CME Group, s. f.). Es el contrato de petróleo más líquido del mundo, y por eso el WTI sirve de referencia para fijar precios de otros crudos (Fattouh, 2011).
Es importante precisar que el Brent es el crudo de referencia mundial y el WTI es la referencia de estados unidos, es de resaltar que el brent es comercializado por vía marítima y el WTI es tierra adentro. Por esta razón los precios se separan y vienen en general de cambios en la capacidad de los oleoductos o en la producción de crudo de Estados Unidos (U.S. Energy Information Administration, 2017).
En este trabajo no se usa el futuro sino el precio spot que calcula la Administración de Información de Energía de Estados Unidos y que FRED publica como promedio mensual (U.S. Energy Information Administration, 2026b). Que sea un promedio importa para dos cosas que se retoman más adelante: suaviza los movimientos dentro del mes y no corresponde a un precio al que se pueda comprar en un día específico.
La frecuencia es mensual por dos razones. La primera es que una de las variables exógenas evaluadas, la producción industrial, solo se publica cada mes, y el enunciado pide considerar modelos con variables externas. La segunda es que la decisión de inversión es a un año, y para ese horizonte el ruido diario aporta poco. Con 249 meses hay más de veinte años de historia, cuatro veces el mínimo recomendado.
2.1 Evolución histórica
# Gráfico del precio del WTI en el tiempo: es la primera mirada a la serie (¿tiene tendencia?, ¿saltos?)
plot(datos_macro$WTI, main = "WTI: precio spot promedio mensual (USD/barril)")
La trayectoria general del precio WTI no gira alrededor de un nivel fijo y se alcanzan a apreciar periodos largos de subida y de bajada, se ven caídas fuertes alrededor del 2008, entre el 2014 y 2016 y en 2020 y subidas en el 2007, 2022 y 2026.
# Las cuatro series juntas, cada una con su propia escala, para comparar cómo se mueven
plot(datos_macro, multi.panel = TRUE, yaxis.same = FALSE, main = "Las cuatro series")
En el gráfico se logra identificar que en el 2008 y 2020 el WTI y la producción industrial caen juntos en las crisis, cuando la economía se frena, cae la demanda del petróleo y el dólar a su vez sube cuando el petróleo se cae, porque el petróleo se paga en dólares. En el 2022 la tasa de la FED sube con fuerza y el petróleo empieza a caer.
2.2 Estadística descriptiva y volatilidad
p <- as.numeric(datos_macro$WTI) # el precio como vector de números (sin fechas) para hacer cuentas
r <- diff(log(p)) # log-retorno: cuánto cambió el precio cada mes, en porcentaje aproximado
# Tabla de estadísticos básicos del precio y del retorno (toda la muestra, solo para describir)
tabla_descr <- data.frame(
serie = c("Precio", "Log-retorno"), # nombre de cada fila
media = c(mean(p), mean(r)), # promedio
mediana = c(median(p), median(r)), # valor del medio (no lo mueven los extremos)
desv_est = c(sd(p), sd(r)), # qué tanto se aleja normalmente del promedio
minimo = c(min(p), min(r)), # valor más bajo
maximo = c(max(p), max(r)), # valor más alto
asimetria = c(NA, mean((r - mean(r))^3) / sd(r)^3), # < 0: las caídas fuertes pesan más que las subidas
curtosis = c(NA, mean((r - mean(r))^4) / sd(r)^4) # > 3: colas pesadas (meses extremos más frecuentes que en la normal)
)
knitr::kable(tabla_descr, digits = 4, caption = "Estadísticos descriptivos (muestra completa)") # la muestra como tabla bonita| serie | media | mediana | desv_est | minimo | maximo | asimetria | curtosis |
|---|---|---|---|---|---|---|---|
| Precio | 72.7292 | 71.5800 | 21.0497 | 16.5500 | 133.8800 | NA | NA |
| Log-retorno | 0.0016 | 0.0131 | 0.1083 | -0.5681 | 0.5456 | -0.7505 | 10.6882 |
cat("Coeficiente de variación del precio:", round(sd(p) / mean(p) * 100, 1), "%\n") # desviación en % del promedio## Coeficiente de variación del precio: 28.9 %cat("Volatilidad anualizada del retorno:", round(sd(r) * sqrt(12) * 100, 1), "%\n") # volatilidad mensual llevada a un año (x raíz de 12)## Volatilidad anualizada del retorno: 37.5 %cat("Máximo:", max(p), "en", format(index(datos_macro)[which.max(p)]), "\n") # precio más alto y su fecha## Máximo: 133.88 en 2008-06-01cat("Mínimo:", min(p), "en", format(index(datos_macro)[which.min(p)]), "\n") # precio más bajo y su fecha## Mínimo: 16.55 en 2020-04-01# Volatilidad móvil: desviación del retorno en los últimos 12 meses, mes a mes, llevada a un año
vol_12m <- rollapply(r, width = 12, FUN = sd) * sqrt(12)
fechas_vol <- index(datos_macro)[-(1:12)] # fecha de cada ventana (el último mes de la ventana)
precio_vol <- p[-(1:12)] # precio en esas mismas fechas, para compararlos
# Gráfico: muestra los periodos de calma y de turbulencia del petróleo
plot(fechas_vol, vol_12m, type = "l", xlab = "", ylab = "Volatilidad anualizada",
main = "Volatilidad móvil de 12 meses")
abline(h = sd(r) * sqrt(12), col = "red", lty = 2) # línea roja: volatilidad promedio de toda la muestra
cat("Volatilidad mínima:", round(min(vol_12m) * 100, 1), "% en", format(fechas_vol[which.min(vol_12m)]), "\n") # el periodo más tranquilo## Volatilidad mínima: 13.4 % en 2014-07-01cat("Volatilidad máxima:", round(max(vol_12m) * 100, 1), "% en", format(fechas_vol[which.max(vol_12m)]), "\n") # el periodo más agitado## Volatilidad máxima: 107.7 % en 2021-01-01cat("Correlación volatilidad - precio:", round(cor(vol_12m, precio_vol), 3), "\n") # negativa: más volatilidad cuando el precio está bajo## Correlación volatilidad - precio: -0.494# Tabla con la fecha y el log-retorno de cada mes, para ordenarla
ret_df <- data.frame(fecha = index(datos_macro)[-1], log_retorno = round(r, 4))
head(ret_df[order(ret_df$log_retorno), ], 3) # los 3 peores meses (de menor a mayor retorno)## fecha log_retorno
## 171 2020-04-01 -0.5681
## 170 2020-03-01 -0.5483
## 35 2008-12-01 -0.3320head(ret_df[order(-ret_df$log_retorno), ], 3) # los 3 mejores meses (de mayor a menor)## fecha log_retorno
## 172 2020-05-01 0.5456
## 242 2026-03-01 0.3482
## 173 2020-06-01 0.2937# Precio del WTI en las fechas clave de los tres episodios: crisis 2008, caída 2014-2016 y pandemia 2020
datos_macro$WTI[as.Date(c("2008-06-01", "2008-12-01", "2014-06-01", "2016-02-01",
"2020-02-01", "2020-04-01", "2022-07-01"))]## WTI
## 2008-06-01 133.88
## 2008-12-01 41.12
## 2014-06-01 105.79
## 2016-02-01 30.32
## 2020-02-01 50.54
## 2020-04-01 16.55
## 2022-07-01 101.62En toda la muestra el precio promedió 72,73 USD, con una mediana muy cercana de 71,58 USD. El rango es amplio: desde 16,55 USD en abril de 2020 hasta 133,88 USD en junio de 2008. El coeficiente de variación es de 28,9 %, lo que confirma que es un activo de precio inestable.
El retorno logarítmico mensual tiene media prácticamente nula (0,16 %), en promedio de un mes a otro casi no sube ni baja. La desviación estándar es de 10,8 %: el precio se mueve hacia arriba o hacia abajo alrededor del 10 %, equivalente a una volatilidad anualizada de 37,5 %, que es el mismo movimiento de la desviación estándar pero anualizado. Su asimetría es negativa (−0,75): las caídas grandes son más frecuentes o más fuertes que las subidas grandes, y su curtosis es de 10,69, muy por encima del valor 3 de una distribución normal. Esto significa que los meses extremos son mucho más frecuentes de lo que predeciría una normal: un mes como abril de 2020 casi nunca debería pasar y en el periodo analizado pasó varias veces. Las colas pesadas son un rasgo habitual de los retornos financieros (Tsay, 2010). Para un inversionista, el riesgo del petróleo no está en la variación de un mes cualquiera sino en unos pocos meses muy malos.
La volatilidad tampoco es constante. Calculada en ventanas móviles de doce meses, va desde 13,4 % en la ventana que termina en julio de 2014 hasta 107,7 % en la que termina en enero de 2021. Esta última recoge los meses de marzo a mayo de 2020, que es donde ocurrió la turbulencia. Esta observación anticipa el resultado de la prueba ARCH de la sección 6. Además, la volatilidad y el nivel del precio tienen una correlación de −0,49: el mercado se vuelve más inestable cuando el precio está bajo.
Cinco hechos de mercado ayudan a entender los movimientos más grandes de la serie.
Crisis de 2008. El precio pasó de 133,88 USD en junio a 41,12 USD en diciembre, y diciembre de 2008 es el tercer peor mes de la serie (−0,332). Hamilton (2009) explica la subida previa por una demanda mundial fuerte frente a una producción que no crecía, y la caída posterior por el desplome de la actividad económica durante la crisis financiera.
Caída de 2014 a 2016. Entre junio de 2014 y febrero de 2016 el precio bajó de 105,79 a 30,32 USD. Baumeister y Kilian (2016) atribuyen la caída a una combinación de mayor oferta de crudo, acumulada desde antes de mediados de 2014, y un debilitamiento de la economía mundial, a lo que se sumó un cambio en las expectativas del mercado a finales de ese año.
Pandemia de 2020. El precio bajó de 50,54 USD en febrero a 16,55 USD en abril. Marzo y abril de 2020 son los dos peores meses de la serie (−0,548 y −0,568) y mayo de 2020 es el mejor (+0,546). El 20 de abril el contrato de futuros con vencimiento en mayo cerró en valores negativos por la falta de liquidez y de capacidad de almacenamiento (U.S. Energy Information Administration, 2020). Esto ocurrió en el futuro y no en el promedio mensual del spot, que nunca fue menor que 16,55 USD; por eso el logaritmo de la serie está definido en toda la muestra.
Invasión de Ucrania en 2022. El precio de entrada de este trabajo, 101,62 USD en julio de 2022, es alto frente al promedio histórico. En el primer semestre de 2022 la tensión geopolítica con Rusia, que terminó en la invasión de Ucrania el 24 de febrero, empujó los precios al alza en un momento en que los inventarios mundiales venían cayendo; en el segundo semestre los precios bajaron (U.S. Energy Information Administration, 2023). La fecha de corte queda entonces cerca de un máximo.
Conflicto en Medio Oriente en 2026. Marzo de 2026 es el segundo mejor mes de la muestra (+0,348): el promedio pasó de 64,51 a 91,38 USD. Este mes está dentro del periodo de prueba. El alza se debió al conflicto en Medio Oriente y a las interrupciones en el estrecho de Ormuz, aunque los inventarios altos de Estados Unidos limitaron la subida del WTI frente al Brent (U.S. Energy Information Administration, 2026a).
2.3 Variables exógenas e hipótesis de inversión
par(mfrow = c(3, 1), mar = c(3, 4, 2, 1)) # pone 3 gráficos uno debajo del otro, con márgenes pequeños
# Índice amplio del dólar: si el dólar se fortalece, el petróleo (que se paga en dólares) tiende a bajar
plot(index(datos_macro), as.numeric(datos_macro$DTWEXBGS), type = "l", xlab = "", ylab = "Índice",
main = "Índice amplio del dólar (DTWEXBGS)")
# Producción industrial: más producción = más demanda de energía = presión al alza en el petróleo
plot(index(datos_macro), as.numeric(datos_macro$INDPRO), type = "l", xlab = "", ylab = "Índice",
main = "Producción industrial de EE. UU. (INDPRO)")
# Tasa de la Fed: refleja el ciclo económico y el costo de financiar inventarios
plot(index(datos_macro), as.numeric(datos_macro$FEDFUNDS), type = "l", xlab = "", ylab = "%",
main = "Tasa de fondos federales (FEDFUNDS)")
par(mfrow = c(1, 1)) # vuelve a un gráfico por ventana
# Correlación de los cambios mensuales: ¿se mueven juntos? (solo descriptivo, no es causalidad)
round(cor(coredata(datos_diff)), 3)## dWTI dDTWEXBGS dINDPRO dFEDFUNDS
## dWTI 1.000 -0.421 0.436 0.256
## dDTWEXBGS -0.421 1.000 -0.159 -0.186
## dINDPRO 0.436 -0.159 1.000 0.280
## dFEDFUNDS 0.256 -0.186 0.280 1.000Se consideraron tres variables externas. El índice amplio del dólar nominal de la Reserva Federal (Board of Governors of the Federal Reserve System, 2026c), porque el crudo se cotiza en dólares y un dólar más fuerte lo encarece para el resto del mundo; se espera una relación negativa. La producción industrial de Estados Unidos (Board of Governors of the Federal Reserve System, 2026b), como medida de la demanda de energía; se espera una relación positiva. Y la tasa de fondos federales (Board of Governors of the Federal Reserve System, 2026a), como indicador del ciclo económico y de la política monetaria.
Las correlaciones entre los cambios mensuales van en la dirección esperada para las dos primeras: −0,42 con el dólar y +0,44 con la producción industrial. Con la tasa de la Fed la correlación es positiva (+0,26). Este signo no significa que subir tasas eleve el precio del petróleo. La explicación más razonable es que la Fed sube tasas cuando la economía y la inflación están fuertes, que es también cuando el petróleo tiende a subir: la tasa responde al mismo ciclo que mueve al crudo. Estas correlaciones se calcularon con toda la muestra y tienen un propósito solo descriptivo.
Hipótesis de inversión. Antes de estimar los modelos se fijó la siguiente regla: si el precio mediano que proyecta el modelo a doce meses supera en más de 5 % al precio de entrada, se compra con los USD 100.000 y se mantiene la posición un año; en caso contrario, el capital queda en liquidez.
Se usa la mediana y no la media porque, al modelar el logaritmo, la distribución del precio futuro es asimétrica: no puede bajar de cero pero no tiene techo. La media queda inflada por unos pocos escenarios muy altos, mientras que la mediana representa el caso típico, el que deja la mitad de los escenarios por debajo (Hyndman y Athanasopoulos, 2021).
El umbral del 5 % es la prima mínima que el equipo exige para asumir el riesgo de un activo con volatilidad anual cercana al 37 %. Es un criterio del equipo, fijado antes de estimar, y no una estimación del costo del dinero. En coherencia, la liquidez se evalúa con rendimiento de 0 %, como supuesto conservador.
El trabajo maneja dos horizontes distintos. El periodo de prueba es de 50 meses y sirve para evaluar la calidad del pronóstico. La decisión de inversión y la simulación son a 12 meses, que es un plazo en el que tiene sentido exigir un retorno anual.
3. Base de datos y tratamiento
if (file.exists("datos_procesados.csv")) { # ¿existe el archivo en la carpeta?
base <- read.csv("datos_procesados.csv") # lo lee
datos_macro <- xts(as.matrix(base[, c("WTI", "DTWEXBGS", "INDPRO", "FEDFUNDS")]),
order.by = as.Date(base$fecha)) # lo convierte en serie de tiempo con fechas
} else { # si no existe, descarga todo de FRED
getSymbols(c("MCOILWTICO", "DTWEXBGS", "INDPRO", "FEDFUNDS"), src = "FRED") # WTI, dólar, producción y tasa
DOLAR_m <- apply.monthly(na.omit(DTWEXBGS), mean) # el dólar viene diario: lo pasamos a promedio mensual
index(DOLAR_m) <- as.Date(cut(index(DOLAR_m), "month")) # le ponemos fecha del día 1, como las demás series
datos_macro <- na.omit(merge(MCOILWTICO, DOLAR_m, INDPRO, FEDFUNDS)) # une las 4 series y deja solo meses completos
colnames(datos_macro) <- c("WTI", "DTWEXBGS", "INDPRO", "FEDFUNDS") # nombres claros para cada columna
datos_macro <- datos_macro[paste0("/", FECHA_TOPE)] # corta en septiembre de 2026
write.csv(data.frame(fecha = index(datos_macro), coredata(datos_macro)),
"datos_procesados.csv", row.names = FALSE) # guarda la copia fija para la próxima vez
}# Revisiones que pide el enunciado: fechas, duplicados, faltantes y unidades
cat("Observaciones:", nrow(datos_macro), "\n") # cuántos meses hay (deben ser 249)## Observaciones: 249cat("Rango:", format(range(index(datos_macro))), "\n") # primera y última fecha## Rango: 2006-01-01 2026-09-01cat("Fechas duplicadas:", sum(duplicated(index(datos_macro))), "\n") # ningún mes debe repetirse## Fechas duplicadas: 0cat("Valores faltantes:", sum(is.na(datos_macro)), "\n") # no debe haber huecos## Valores faltantes: 0meses <- seq(start(datos_macro), end(datos_macro), by = "month") # la lista de todos los meses que DEBERÍAN estar
cat("Meses ausentes en la secuencia:", length(setdiff(meses, index(datos_macro))), "\n") # ¿falta alguno?## Meses ausentes en la secuencia: 0cat("Meses con WTI <= 0:", sum(datos_macro$WTI <= 0), "\n") # si hubiera precios negativos, no se podría sacar logaritmo## Meses con WTI <= 0: 0summary(coredata(datos_macro)) # mínimo, máximo, media y cuartiles de cada variable## WTI DTWEXBGS INDPRO FEDFUNDS
## Min. : 16.55 Min. : 86.32 Min. : 84.56 Min. :0.050
## 1st Qu.: 57.52 1st Qu.: 93.76 1st Qu.: 98.20 1st Qu.:0.120
## Median : 71.58 Median :110.85 Median :100.16 Median :0.650
## Mean : 72.73 Mean :106.54 Mean : 98.85 Mean :1.805
## 3rd Qu.: 87.90 3rd Qu.:117.10 3rd Qu.:101.44 3rd Qu.:3.750
## Max. :133.88 Max. :128.84 Max. :104.10 Max. :5.3303.1 Transformación
# Cambios mensuales de cada variable (para describirlas y compararlas)
dWTI <- diff(log(datos_macro$WTI)) # log-retorno del petróleo: cambio porcentual aproximado
dDTWEXBGS <- diff(log(datos_macro$DTWEXBGS)) # cambio porcentual del índice del dólar
dINDPRO <- diff(log(datos_macro$INDPRO)) # cambio porcentual de la producción industrial
dFEDFUNDS <- diff(datos_macro$FEDFUNDS) # la tasa ya está en %, así que se usa la diferencia simple
datos_diff <- na.omit(merge(dWTI, dDTWEXBGS, dINDPRO, dFEDFUNDS)) # une los cambios y quita el primer mes (no tiene cambio)
colnames(datos_diff) <- c("dWTI", "dDTWEXBGS", "dINDPRO", "dFEDFUNDS") # nombres de las columnas
# Objetos para los modelos: el precio como serie mensual y las exógenas en niveles.
# Arima() con d = 1 diferencia también las exógenas, así que el modelo trabaja con sus cambios.
y_ts <- ts(as.numeric(datos_macro$WTI), start = c(2006, 1), frequency = 12) # serie mensual (12 por año) desde enero 2006
X <- cbind(ldolar = log(as.numeric(datos_macro$DTWEXBGS)), # logaritmo del índice del dólar
lip = log(as.numeric(datos_macro$INDPRO)), # logaritmo de la producción industrial
ff = as.numeric(datos_macro$FEDFUNDS)) # tasa de la Fed (ya en %, sin logaritmo)Las cuatro series se descargaron de FRED y se guardaron en una copia fija, porque la producción industrial se revisa hacia atrás; con la copia, cualquier persona obtiene los mismos resultados.
| Variable | Código FRED | Fuente original | Unidad | Frecuencia original | Transformación |
|---|---|---|---|---|---|
| Precio spot del WTI | MCOILWTICO | U.S. Energy Information Administration (2026b) | USD por barril | Mensual (promedio) | Diferencia del logaritmo |
| Índice amplio del dólar | DTWEXBGS | Board of Governors (2026c) | Índice, enero de 2006 = 100 | Diaria, llevada a promedio mensual | Diferencia del logaritmo |
| Producción industrial | INDPRO | Board of Governors (2026b) | Índice, desestacionalizado | Mensual | Diferencia del logaritmo |
| Tasa de fondos federales | FEDFUNDS | Board of Governors (2026a) | Porcentaje | Mensual | Diferencia simple |
La ventana empieza en enero de 2006 porque ahí comienza el índice amplio del dólar, y termina en septiembre de 2026, el último mes con las cuatro series disponibles. La verificación no encontró problemas: 249 meses, sin fechas duplicadas, sin valores faltantes, sin meses ausentes en la secuencia y sin precios menores o iguales a cero.
Para el precio y los dos índices se usa la diferencia del logaritmo, que se lee como cambio porcentual aproximado. Para la tasa de interés, que ya está en porcentaje, se usa la diferencia simple en puntos porcentuales. La variable dependiente y las exógenas deben tener el mismo orden de integración: si se relacionan en niveles series que tienen tendencia estocástica, pueden aparecer relaciones que en realidad no existen, lo que se conoce como regresión espuria (Gujarati y Porter, 2010). Por eso el modelo trabaja con los cambios de todas las variables.
4. Fundamentación matemática
Operador de rezago. Para escribir los modelos de forma compacta se usa el operador de rezago \(B\), que desplaza la serie un periodo hacia atrás: \(B\,y_t = y_{t-1}\), y en general \(B^k y_t = y_{t-k}\). Con él, la primera diferencia se escribe \((1-B)\,y_t = y_t - y_{t-1}\). En algunos textos este operador se llama \(L\); aquí se usa \(B\), como en las diapositivas del curso (Rodríguez Guevara, 2026b).
Ruido blanco. Es la pieza básica de todos los modelos. Es una sucesión \(\varepsilon_t\) con \(E(\varepsilon_t)=0\), \(Var(\varepsilon_t)=\sigma^2\) y \(Cov(\varepsilon_t,\varepsilon_{t-k})=0\) para \(k \neq 0\). Representa la parte de la serie que no se puede predecir con su pasado: la noticia nueva de cada mes.
Modelo AR(p). En un proceso autorregresivo el valor de hoy depende de sus propios valores pasados:
\[y_t = c + \phi_1 y_{t-1} + \dots + \phi_p y_{t-p} + \varepsilon_t\]
Con el operador de rezago queda \(\Phi_p(B)\,y_t = c + \varepsilon_t\), donde \(\Phi_p(B) = 1 - \phi_1 B - \dots - \phi_p B^p\). El parámetro \(p\) es el número de rezagos. El proceso es estacionario si las raíces de \(\Phi_p(B)=0\) tienen módulo mayor que 1; en un AR(1) eso equivale a \(|\phi_1|<1\). Si la condición no se cumple, los choques no se desvanecen y la serie no vuelve a una media (Box et al., 2015).
Modelo MA(q). En un proceso de medias móviles el valor de hoy depende de los choques recientes:
\[y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}\]
Es decir, \(y_t = \mu + \Theta_q(B)\,\varepsilon_t\). El parámetro \(q\) indica cuántos periodos dura el efecto de un choque. Un MA siempre es estacionario; lo que se le exige es invertibilidad: que las raíces de \(\Theta_q(B)=0\) tengan módulo mayor que 1. La invertibilidad permite recuperar los choques a partir de los datos observados, y sin ella el modelo no sirve para pronosticar (Rodríguez Guevara, 2026b).
Modelo ARIMA(p,d,q). Muchas series financieras no son estacionarias en nivel pero sí lo son sus diferencias (Brooks, 2008). El parámetro \(d\) es el número de veces que hay que diferenciar. El modelo general combina las tres piezas:
\[\Phi_p(B)\,(1-B)^d\,y_t = c + \Theta_q(B)\,\varepsilon_t\]
En este trabajo \(y_t = \ln P_t\), el logaritmo del precio. Con \(d=1\), \((1-B)\ln P_t = \ln P_t - \ln P_{t-1}\), que es el retorno logarítmico mensual. Un ARIMA(p,1,q) sobre el logaritmo del precio es entonces un ARMA(p,q) sobre el retorno.
El modelo aplicado. El modelo seleccionado es un ARIMA(0,1,1) sin constante. Partiendo de la forma general con \(p=0\), \(d=1\) y \(q=1\):
\[(1-B)\ln P_t = (1+\theta B)\,\varepsilon_t \quad\Longrightarrow\quad \ln P_t - \ln P_{t-1} = \varepsilon_t + \theta\,\varepsilon_{t-1}\]
Con el valor estimado, el retorno de cada mes es el choque de ese mes más 0,3194 veces el choque del mes anterior. La raíz de \(1 + 0{,}3194B = 0\) es \(B = -3{,}13\), de módulo mayor que 1, así que el modelo es invertible. La autocorrelación teórica de primer orden de un MA(1) es \(\rho_1 = \theta/(1+\theta^2) = 0{,}29\), que coincide con la observada en los datos (0,294).
De esta ecuación sale la forma del pronóstico. Para el mes siguiente, el choque futuro tiene valor esperado cero y solo queda \(\theta\) por el último choque conocido: \(E_T(\Delta \ln P_{T+1}) = \theta\,\hat\varepsilon_T\). Para dos meses o más, todos los choques involucrados son futuros y el cambio esperado es cero. El pronóstico se ajusta una vez y después queda constante. Esta propiedad explica los resultados de las secciones 7 y 8.
Estructura estacional. Un modelo SARIMA\((p,d,q)(P,D,Q)_s\) agrega una segunda capa de dinámica con periodicidad \(s\):
\[\Phi_p(B)\,\Phi_P(B^s)\,(1-B)^d\,(1-B^s)^D\,y_t = \Theta_q(B)\,\Theta_Q(B^s)\,\varepsilon_t\]
Con datos mensuales \(s=12\). La dinámica regular compara un mes con el anterior; la estacional compara un mes con el mismo mes del año pasado. \(P\) y \(Q\) son los órdenes autorregresivo y de media móvil estacionales, y \(D\) es el número de diferencias estacionales, \((1-B^{12})\,y_t = y_t - y_{t-12}\). Como se muestra en la sección 5.5, en el WTI no hay evidencia de estacionalidad, por lo que el modelo final no incluye estos términos.
Variables exógenas. Un ARIMAX agrega variables explicativas \(X_t\). Siguiendo el planteamiento del curso, es una regresión cuyos errores siguen un proceso ARIMA (Rodríguez Guevara, 2026b):
\[y_t = \beta' X_t + n_t, \qquad n_t \sim ARIMA(p,d,q)\]
Cuando \(d=1\), al diferenciar ambos lados queda \(\Delta y_t = \beta' \Delta X_t + \Delta n_t\). Para el caso evaluado aquí:
\[\Delta \ln P_t = \beta_1\,\Delta \ln D_t + \beta_2\,\Delta \ln IP_t + \beta_3\,\Delta FF_t + \varepsilon_t + \theta\,\varepsilon_{t-1}\]
donde \(D_t\) es el índice del dólar, \(IP_t\) la producción industrial y \(FF_t\) la tasa de fondos federales. Las \(X_t\) pueden aportar información que el pasado del precio no contiene, como el estado de la demanda o el valor del dólar. Cada \(\beta\) mide cuánto cambia el retorno del WTI cuando cambia la variable en el mismo mes. Es una asociación estadística, no una relación causal: las tres variables responden a su vez al ciclo económico y al propio precio del petróleo. Además hay un problema práctico: para pronosticar \(P_{T+h}\) se necesitan los valores de \(X_{T+h}\), que en la fecha de corte no se conocen y deben suponerse o pronosticarse. El SARIMAX combina esta regresión con la estructura estacional anterior.
Varianza condicional: ARCH y GARCH. Los modelos anteriores describen la media del retorno y suponen que la varianza del error es constante. Cuando la volatilidad se agrupa en el tiempo ese supuesto no se cumple, y se agrega una segunda ecuación para la varianza. Siguiendo la notación del curso (Rodríguez Guevara, 2026a), el error se escribe como \(\varepsilon_t = \sqrt{h_t}\,z_t\), donde \(z_t\) es ruido blanco de varianza uno y \(h_t\) es la varianza condicional. En el modelo ARCH(1) de Engle (1982) esa varianza depende del choque al cuadrado del periodo anterior:
\[h_t = \omega + \alpha_1\,\varepsilon_{t-1}^2, \qquad \omega > 0,\; \alpha_1 \geq 0\]
Un choque grande hoy, positivo o negativo, eleva la varianza de mañana. Si \(\alpha_1 < 1\) existe una varianza de largo plazo igual a \(\omega/(1-\alpha_1)\). El modelo GARCH(1,1) de Bollerslev (1986) agrega la varianza del periodo anterior:
\[h_t = \omega + \alpha_1\,\varepsilon_{t-1}^2 + \beta_1\,h_{t-1}, \qquad \omega > 0,\; \alpha_1 \geq 0,\; \beta_1 \geq 0,\; \alpha_1 + \beta_1 < 1\]
La suma \(\alpha_1 + \beta_1\) se llama persistencia y mide qué tan lento se desvanece un episodio de alta volatilidad. Estos modelos no cambian la forma de la ecuación de la media, sino la medida de la incertidumbre alrededor del pronóstico. Sin embargo, como la media y la varianza se estiman en conjunto, el valor de los coeficientes de la media sí puede cambiar.
5. Identificación y calibración
5.1 Partición train / test
n <- length(y_ts) # total de meses: 249
n_test <- round(PROP_TEST * n) # 20 % de 249 = 49,8 -> 50 meses para el test
n_train <- n - n_test # 199 meses para entrenar (estimar) los modelos
H <- n_test # 50: horizonte del TEST, para evaluar el pronóstico (sección 7)
# Partición cronológica, como en las diapositivas (window): el pasado entrena, el futuro evalúa
train <- window(y_ts, end = time(y_ts)[n_train]) # enero 2006 a julio 2022
test <- window(y_ts, start = time(y_ts)[n_train + 1]) # agosto 2022 a septiembre 2026 (el modelo NO lo ve)
X_train <- X[1:n_train, ] # exógenas del mismo periodo de train
X_test <- X[(n_train + 1):n, ] # exógenas del test (solo para la discusión de leakage)
fechas <- index(datos_macro) # todas las fechas
T0 <- fechas[n_train] # fecha de corte: julio de 2022 (desde aquí "no sabemos" el futuro)
fechas_test <- fechas[(n_train + 1):n] # fechas del test
P0 <- as.numeric(tail(train, 1)) # precio de entrada: el último precio conocido en T0
cat("Train:", format(fechas[1]), "a", format(T0), "(", n_train, "meses )\n") # confirma el train## Train: 2006-01-01 a 2022-07-01 ( 199 meses )cat("Test: ", format(fechas_test[1]), "a", format(fechas[n]), "(", n_test, "meses )\n") # confirma el test## Test: 2022-08-01 a 2026-09-01 ( 50 meses )cat("Porcentaje de test:", round(n_test / n * 100, 1), "%\n") # debe estar entre 15 % y 25 %## Porcentaje de test: 20.1 %cat("Precio de entrada en T0:", P0, "USD/barril\n") # precio al que compraríamos## Precio de entrada en T0: 101.62 USD/barrilcat("Horizonte de la decisión:", H_DEC, "meses (", format(fechas_test[1]), "a", format(fechas_test[H_DEC]), ")\n") # agosto 2022 a julio 2023## Horizonte de la decisión: 12 meses ( 2022-08-01 a 2023-07-01 )# Gráfico de la partición: en un color lo que se usa para estimar y en otro lo que se guarda para evaluar
autoplot(y_ts, series = "Train") + autolayer(test, series = "Test") +
labs(title = "Partición cronológica", x = NULL, y = "USD/barril") + theme_minimal()
La muestra se dividió respetando el orden temporal. El entrenamiento va de enero de 2006 a julio de 2022 (199 meses) y la prueba de agosto de 2022 a septiembre de 2026 (50 meses, 20,1 %). La fecha de corte es julio de 2022 y el precio de entrada, 101,62 USD por barril.
Se reservó el 20 %, en la mitad del rango recomendado. Con 15 % la prueba tendría 37 meses, pocos para evaluar un pronóstico de largo plazo; con 25 % el entrenamiento perdería un año de historia. La partición se hizo antes de cualquier cálculo que dependa de los datos: el parámetro de Box-Cox, las pruebas de raíz unitaria, los correlogramas, los coeficientes y la selección del modelo se obtienen únicamente con el entrenamiento. Usar información de la prueba en cualquiera de esos pasos sería fuga de información y haría parecer el modelo mejor de lo que es (Hyndman y Athanasopoulos, 2021).
En el periodo de prueba el precio bajó desde el nivel de entrada hasta cerca de 70 USD a mediados de 2023, siguió cayendo hasta un mínimo de 57,97 USD en diciembre de 2025 y repuntó con fuerza desde marzo de 2026, para terminar en 97,31 USD. La fecha de corte quedó cerca de un máximo, lo que favorece en esta prueba a cualquier pronóstico que parta de un nivel más bajo.
5.2 Box-Cox y transformación
# Box-Cox (lo pide el enunciado): busca la transformación que estabiliza la varianza. lambda = 0 es logaritmo
cat("Lambda Box-Cox (Guerrero):", round(BoxCox.lambda(train, method = "guerrero"), 3), "\n") # método 1## Lambda Box-Cox (Guerrero): 0.752cat("Lambda Box-Cox (verosimilitud):", round(BoxCox.lambda(train, method = "loglik"), 3), "\n") # método 2## Lambda Box-Cox (verosimilitud): 0.55ly <- log(train) # logaritmo del precio (solo train): así los cambios se leen en porcentaje
dly <- diff(ly) # primera diferencia del logaritmo = log-retorno mensual (solo train)
# Control de sobrediferenciación: si diferenciamos otra vez y la varianza SUBE, una sola diferencia (d = 1) basta
data.frame(serie = c("1ra diferencia del log", "2da diferencia del log"),
varianza = c(var(dly), var(diff(dly))))## serie varianza
## 1 1ra diferencia del log 0.01294457
## 2 2da diferencia del log 0.01827132El parámetro de Box-Cox estimado con el entrenamiento fue 0,752 por el método de Guerrero y 0,55 por máxima verosimilitud. Ninguno está cerca de 0, que es el valor que corresponde al logaritmo; los dos apuntan a una transformación intermedia, parecida a la raíz cuadrada. Box-Cox, por sí solo, no respalda el logaritmo.
Aun así se trabaja con logaritmos por tres razones de orden financiero. La diferencia del logaritmo es el retorno compuesto continuo, que es la medida estándar de rendimiento (Tsay, 2010). Al devolver el pronóstico a la escala original el precio siempre resulta positivo. Y es la transformación que se usa en el curso para series de precios (Rodríguez Guevara, 2026b). Se reconoce que es una decisión por criterio financiero y de interpretación, no por resultado estadístico.
Para comprobar que una sola diferencia basta, se comparó la varianza de la primera diferencia del logaritmo (0,0129) con la de la segunda (0,0183). Al diferenciar otra vez la varianza sube, que es la señal de sobrediferenciación descrita en clase.
5.3 Estacionariedad
# acfinter() es la función del profesor: entrega ADF, KPSS y PP de una vez (y el ACF/PACF)
act_nivel <- acfinter(ly, lag = 36, ci.method = "white") # pruebas al logaritmo del precio
act_dif <- acfinter(dly, lag = 36, ci.method = "white") # pruebas al log-retorno
# ADF y PP: H0 = no estacionaria (queremos p < 0,05). KPSS: H0 = estacionaria (queremos p > 0,05)
cat("log(precio):\n"); print(act_nivel[["Stationary Test"]]) # esperamos: NO estacionaria## log(precio):## Statistic P_Value
## ADF -2.6253642 0.3142249
## KPSS-Level 0.8646651 0.0100000
## KPSS-Trend 0.2248383 0.0100000
## PP -15.9297003 0.2077812cat("\nlog-retorno:\n"); print(act_dif[["Stationary Test"]]) # esperamos: SÍ estacionaria##
## log-retorno:## Statistic P_Value
## ADF -6.21727765 0.01
## KPSS-Level 0.05893761 0.10
## KPSS-Trend 0.05184484 0.10
## PP -124.92309733 0.01# ndiffs() (diapositivas) dice cuántas diferencias necesita la serie para ser estacionaria
cat("\nndiffs(log precio) =", ndiffs(ly), "| ndiffs(log-retorno) =", ndiffs(dly), "\n") # 1 y 0 -> d = 1##
## ndiffs(log precio) = 1 | ndiffs(log-retorno) = 0Las pruebas se obtuvieron con la función acfinter del
paquete actfts (Rodríguez y Sierra, 2025). Se aplicaron tres pruebas con
hipótesis nulas opuestas. En ADF y Phillips-Perron la hipótesis nula es
que hay raíz unitaria; en KPSS es que la serie es estacionaria. Cuando
las tres llegan a la misma conclusión, la evidencia es más sólida
(Rodríguez Guevara, 2026b).
En el logaritmo del precio, ADF (p = 0,314) y PP (p = 0,208) no
rechazan la raíz unitaria, y KPSS rechaza la estacionariedad tanto en
nivel como en tendencia (p menor que 0,01). En el retorno ocurre lo
contrario: ADF y PP rechazan la raíz unitaria (p menor que 0,01) y KPSS
no rechaza la estacionariedad (p mayor que 0,10). Los valores 0,01 y
0,10 son los límites de las tablas que usa R, no p-valores exactos. La
conclusión es que el logaritmo del precio es integrado de orden uno,
\(d = 1\), y ndiffs lo
confirma.
Estacionariedad y estacionalidad de las exógenas
# Mismas pruebas para las exógenas (solo train): Y y X deben tener el mismo orden de integración
estac <- function(x) { # función que hace todas las pruebas a una serie
x <- as.numeric(x) # la serie en nivel
dx <- diff(x) # su primera diferencia
a <- acf(dx, lag.max = 36, plot = FALSE)$acf # autocorrelaciones de la diferencia (para ver estacionalidad)
c(ADF_nivel = adf.test(x)$p.value, KPSS_L_nivel = kpss.test(x, null = "Level")$p.value, # pruebas en nivel
KPSS_T_nivel = kpss.test(x, null = "Trend")$p.value, PP_nivel = pp.test(x)$p.value,
ADF_dif = adf.test(dx)$p.value, KPSS_L_dif = kpss.test(dx, null = "Level")$p.value, # pruebas en diferencia
KPSS_T_dif = kpss.test(dx, null = "Trend")$p.value, PP_dif = pp.test(dx)$p.value,
ndiffs = ndiffs(x), nsdiffs = nsdiffs(ts(x, frequency = 12)), # diferencias sugeridas: normal y estacional
ACF12_dif = a[13], ACF24_dif = a[25], ACF36_dif = a[37]) # ACF en los rezagos 12, 24 y 36
}
series_est <- list("log WTI" = as.numeric(ly), # las 4 series que comparamos
"log DTWEXBGS" = X_train[, "ldolar"],
"log INDPRO" = X_train[, "lip"],
"FEDFUNDS" = X_train[, "ff"])
tabla_est <- t(sapply(series_est, estac)) # aplica la función a cada serie y arma la tabla
knitr::kable(t(round(tabla_est, 3)),
caption = "p-valores de estacionariedad (nivel y primera diferencia), d, D y ACF estacional (train)")| log WTI | log DTWEXBGS | log INDPRO | FEDFUNDS | |
|---|---|---|---|---|
| ADF_nivel | 0.314 | 0.138 | 0.371 | 0.474 |
| KPSS_L_nivel | 0.010 | 0.010 | 0.010 | 0.010 |
| KPSS_T_nivel | 0.010 | 0.010 | 0.010 | 0.010 |
| PP_nivel | 0.208 | 0.516 | 0.463 | 0.926 |
| ADF_dif | 0.010 | 0.010 | 0.010 | 0.046 |
| KPSS_L_dif | 0.100 | 0.100 | 0.100 | 0.077 |
| KPSS_T_dif | 0.100 | 0.100 | 0.100 | 0.100 |
| PP_dif | 0.010 | 0.010 | 0.010 | 0.010 |
| ndiffs | 1.000 | 1.000 | 1.000 | 1.000 |
| nsdiffs | 0.000 | 0.000 | 0.000 | 0.000 |
| ACF12_dif | -0.017 | -0.068 | -0.039 | 0.071 |
| ACF24_dif | -0.065 | -0.046 | -0.074 | -0.136 |
| ACF36_dif | -0.044 | 0.078 | -0.051 | -0.020 |
Las tres variables externas se comportan igual que el WTI: no estacionarias en nivel y estacionarias en primera diferencia. En el dólar y la producción industrial el resultado es claro. En la tasa de la Fed es más justo: en diferencias, ADF da p = 0,046 y KPSS en nivel p = 0,077. Pasa al 5 %, con poco margen, lo que es coherente con una tasa que se mueve en ciclos largos de subidas y bajadas. Al ser todas I(1), se modelan en diferencias junto con el precio.
Ninguna muestra estacionalidad: nsdiffs es cero en las
cuatro series y ninguna autocorrelación en los rezagos 12, 24 y 36
supera la banda de ±0,139. En el caso de la producción industrial esto
era de esperar, porque FRED la publica desestacionalizada.
5.4 Identificación por ACF y PACF
tab_id <- as.data.frame(act_dif[["ACF-PACF Test"]]) # tabla de acfinter: ACF, PACF y Ljung-Box por rezago
banda <- 1.96 / sqrt(length(dly)) # banda al 95 %: 1,96 / raíz de 198 = 0,139
cat("Banda de significancia al 95 %: ±", round(banda, 4), "\n") # lo que pase de aquí NO es casualidad## Banda de significancia al 95 %: ± 0.1393round(head(tab_id, 12), 4) # los primeros 12 rezagos## lag acf pacf Box_Pierce Pv_Box Ljung_Box Pv_Ljung
## 1 1 0.2940 0.2940 17.1154 0.0000 17.3760 0.0000
## 2 2 -0.0290 -0.1264 17.2824 0.0002 17.5465 0.0002
## 3 3 -0.1196 -0.0810 20.1126 0.0002 20.4492 0.0001
## 4 4 -0.1014 -0.0472 22.1490 0.0002 22.5486 0.0002
## 5 5 -0.0452 -0.0152 22.5542 0.0004 22.9684 0.0003
## 6 6 -0.0362 -0.0420 22.8138 0.0009 23.2389 0.0007
## 7 7 0.0052 0.0127 22.8191 0.0018 23.2444 0.0015
## 8 8 -0.0376 -0.0630 23.0985 0.0032 23.5386 0.0027
## 9 9 -0.0511 -0.0346 23.6152 0.0050 24.0854 0.0042
## 10 10 -0.0330 -0.0183 23.8312 0.0081 24.3151 0.0068
## 11 11 -0.0350 -0.0391 24.0740 0.0124 24.5748 0.0105
## 12 12 -0.0166 -0.0160 24.1285 0.0195 24.6334 0.0167par(mfrow = c(3, 1), mar = c(4, 4, 2, 1)) # 3 gráficos en columna
# ACF: si se corta en el rezago q, sugiere un MA(q)
plot(tab_id$lag, tab_id$acf, type = "h", lwd = 2, xlab = "rezago", ylab = "ACF", main = "ACF del log-retorno (train)")
abline(h = 0); abline(h = c(-1, 1) * banda, lty = 2, col = "red") # línea en 0 y bandas rojas
# PACF: si se corta en el rezago p, sugiere un AR(p)
plot(tab_id$lag, tab_id$pacf, type = "h", lwd = 2, xlab = "rezago", ylab = "PACF", main = "PACF")
abline(h = 0); abline(h = c(-1, 1) * banda, lty = 2, col = "red")
# Ljung-Box: si p < 0,05, el retorno NO es ruido blanco -> hay estructura que modelar
plot(tab_id$lag, tab_id$Pv_Ljung, pch = 19, ylim = c(0, 1), xlab = "rezago", ylab = "p-valor",
main = "p-valor de Ljung-Box")
abline(h = 0.05, lty = 2, col = "red") # línea del 5 %
par(mfrow = c(1, 1)) # vuelve a un gráfico por ventanaLa banda de significancia al 95 % es \(\pm 1{,}96/\sqrt{198} = \pm 0{,}139\). En el retorno, solo el primer rezago la supera, con 0,294 tanto en la ACF como en la PACF. El segundo rezago de la PACF (−0,126) queda dentro. La prueba de Ljung-Box rechaza que el retorno sea ruido blanco, de modo que hay una estructura por modelar, aunque es corta.
En teoría, un MA(q) se reconoce porque la ACF se corta en el rezago \(q\) y un AR(p) porque la PACF se corta en el rezago \(p\) (Box et al., 2015). Con un solo rezago significativo en ambas no se puede distinguir un AR(1) de un MA(1). Por eso se proponen los dos, junto con el ARMA(1,1) que los combina. Se agregan como control tres modelos más grandes, ARMA(2,1), ARMA(1,2) y ARMA(2,2), para comprobar que los rezagos adicionales no aportan.
5.5 Estacionalidad
tab_id[c(12, 24, 36), c("lag", "acf", "pacf")] # ¿hay picos cada 12 meses? (compararlos con la banda de 0,139)## lag acf pacf
## 12 12 -0.01659012 -0.015950898
## 24 24 -0.06534672 0.001085151
## 36 36 -0.04375781 -0.051667337cat("nsdiffs(log precio) =", nsdiffs(ly), "\n") # diapositivas: diferencias estacionales necesarias (0 = no hay estacionalidad)## nsdiffs(log precio) = 0round(tapply(dly, cycle(dly), mean) * 100, 3) # retorno promedio de cada mes del año en % (como en las diapositivas)## 1 2 3 4 5 6 7 8 9 10 11 12
## -0.147 1.303 2.239 0.516 5.096 2.661 0.244 -2.691 -0.429 -0.752 -2.823 -3.225# Gráfico por mes: si las líneas horizontales (promedios) están a la misma altura, no hay estacionalidad
ggsubseriesplot(dly) + labs(title = "Log-retorno por mes del año (train)", y = "log-retorno") +
theme_minimal()
Las autocorrelaciones en los rezagos 12, 24 y 36 son −0,017, −0,065 y
−0,044, todas dentro de la banda. nsdiffs indica que no
hace falta ninguna diferencia estacional.
El retorno promedio por mes va de −3,2 % en diciembre a +5,1 % en mayo, una diferencia que a primera vista parece grande. No es evidencia de estacionalidad. Cada mes tiene solo 16 o 17 observaciones y el retorno mensual tiene una desviación cercana a 11 %, por lo que el error estándar de cada promedio ronda los 3 puntos porcentuales. Además, el promedio de mayo depende de un único dato: el rebote de mayo de 2020 (+54,6 %). Sin ese mes, mayo baja a cerca de +2 %. La prueba formal está en la sección siguiente: el coeficiente estacional del SARIMA no es significativo (p = 0,85). Tiene sentido económico. El WTI se negocia en un mercado global con inventarios y con futuros a distintos plazos, de manera que cualquier patrón de calendario conocido se incorpora al precio por anticipado. Se concluye que D=0 y que no corresponde incluir términos estacionales.
5.6 Modelos candidatos
# Todos los modelos se estiman sobre el precio con lambda = 0 (logaritmo) y d = 1
# Un ARIMA(p,1,q) sobre log(precio) es lo mismo que un ARMA(p,q) sobre el log-retorno.
# Cada modelo se describe con: order = (p, d, q), seas = parte estacional (P, D, Q) y x = exógenas que usa
specs <- list(
"ARIMA(1,1,0)" = list(order = c(1, 1, 0), seas = c(0, 0, 0), x = NULL), # AR(1): candidato por el PACF
"ARIMA(0,1,1)" = list(order = c(0, 1, 1), seas = c(0, 0, 0), x = NULL), # MA(1): candidato por el ACF
"ARIMA(1,1,1)" = list(order = c(1, 1, 1), seas = c(0, 0, 0), x = NULL), # ARMA(1,1): mezcla de los dos
"ARIMA(2,1,1)" = list(order = c(2, 1, 1), seas = c(0, 0, 0), x = NULL), # control: ¿aporta un rezago AR más?
"ARIMA(1,1,2)" = list(order = c(1, 1, 2), seas = c(0, 0, 0), x = NULL), # control: ¿aporta un rezago MA más?
"ARIMA(2,1,2)" = list(order = c(2, 1, 2), seas = c(0, 0, 0), x = NULL), # control: modelo grande
"SARIMA(0,1,1)(0,0,1)[12]" = list(order = c(0, 1, 1), seas = c(0, 0, 1), x = NULL), # prueba formal de estacionalidad
"ARIMAX(0,1,1) DTWEXBGS+INDPRO" =
list(order = c(0, 1, 1), seas = c(0, 0, 0), x = c("ldolar", "lip")), # MA(1) + dólar + producción
"ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS" =
list(order = c(0, 1, 1), seas = c(0, 0, 0), x = c("ldolar", "lip", "ff")), # MA(1) + las tres exógenas
"SARIMAX(0,1,1)(0,0,1)[12] DTWEXBGS+INDPRO" =
list(order = c(0, 1, 1), seas = c(0, 0, 1), x = c("ldolar", "lip")) # todo junto
)
# Función que estima un modelo SOLO con train. Si el método ML falla, intenta con CSS-ML (diapositivas)
ajustar <- function(sp, drift = TRUE) {
args <- list(y = train, order = sp$order, # serie y orden (p, d, q)
seasonal = list(order = sp$seas, period = 12), # parte estacional, periodo 12 meses
xreg = if (is.null(sp$x)) NULL else X_train[, sp$x, drop = FALSE], # exógenas, si el modelo las usa
include.drift = drift, lambda = 0, method = "ML") # drift = tendencia; lambda 0 = logaritmo
m <- tryCatch(do.call(Arima, args), error = function(e) NULL) # intenta estimar; si falla, devuelve NULL
if (is.null(m)) { # si falló con ML...
args$method <- "CSS-ML" # ...cambia de método
m <- tryCatch(do.call(Arima, args), error = function(e) NULL) # ...y vuelve a intentar
}
if (!is.null(m)) m$series <- "train" # nombre corto para que print() no muestre todos los datos
m # devuelve el modelo estimado
}
modelos <- lapply(specs, ajustar) # estima los 10 modelos
modelos <- modelos[!sapply(modelos, is.null)] # quita los que no se pudieron estimar (si hubo alguno)
cat("Modelos estimados:", length(modelos), "de", length(specs), "\n") # confirma cuántos salieron## Modelos estimados: 10 de 10# coeftest(): para cada coeficiente, Pr(>|z|) < 0,05 significa que sí aporta al modelo
for (nm in names(modelos)) {
cat("\n==========", nm, "==========\n") # título de cada modelo
print(coeftest(modelos[[nm]])) # tabla de coeficientes con su p-valor
}##
## ========== ARIMA(1,1,0) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.2948003 0.0679390 4.3392 0.0000143 ***
## drift 0.0018269 0.0109052 0.1675 0.867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== ARIMA(0,1,1) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3191657 0.0650390 4.9073 0.0000009234 ***
## drift 0.0019213 0.0101012 0.1902 0.8491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== ARIMA(1,1,1) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.0560985 0.1913454 0.2932 0.7694
## ma1 0.2706737 0.1808941 1.4963 0.1346
## drift 0.0018836 0.0103049 0.1828 0.8550
##
##
## ========== ARIMA(2,1,1) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.2611815 0.0667534 18.8931 < 0.00000000000000022 ***
## ar2 -0.3373292 0.0677942 -4.9758 0.0000006498 ***
## ma1 -0.9999997 0.0157966 -63.3047 < 0.00000000000000022 ***
## drift -0.0012481 0.0015556 -0.8023 0.4224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== ARIMA(1,1,2) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.3057372 0.6820474 -0.4483 0.6540
## ma1 0.6388047 0.6769126 0.9437 0.3453
## ma2 0.1390326 0.2225041 0.6249 0.5321
## drift 0.0018628 0.0104163 0.1788 0.8581
##
##
## ========== ARIMA(2,1,2) ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.1062068 0.1835582 6.0265 0.000000001676 ***
## ar2 -0.1898687 0.1777161 -1.0684 0.2853
## ma1 -0.8226693 0.1827679 -4.5012 0.000006758071 ***
## ma2 -0.1773298 0.1821446 -0.9736 0.3303
## drift -0.0011397 0.0016448 -0.6929 0.4884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== SARIMA(0,1,1)(0,0,1)[12] ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3192494 0.0650866 4.9050 0.0000009343 ***
## sma1 -0.0144913 0.0768617 -0.1885 0.8505
## drift 0.0019165 0.0099633 0.1924 0.8475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== ARIMAX(0,1,1) DTWEXBGS+INDPRO ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.1628555 0.0702200 2.3192 0.02038 *
## drift 0.0055652 0.0071148 0.7822 0.43410
## ldolar -3.7775059 0.4839595 -7.8054 0.00000000000000593 ***
## lip 3.1143355 0.4546713 6.8496 0.00000000000740354 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.1718816 0.0685938 2.5058 0.012218 *
## drift 0.0068264 0.0070580 0.9672 0.333448
## ldolar -3.5521012 0.4849436 -7.3248 0.0000000000002393 ***
## lip 2.7663677 0.4644392 5.9564 0.0000000025791636 ***
## ff 0.1151861 0.0431025 2.6724 0.007532 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ========== SARIMAX(0,1,1)(0,0,1)[12] DTWEXBGS+INDPRO ==========
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.1630617 0.0704419 2.3148 0.02062 *
## sma1 -0.0348603 0.0755036 -0.4617 0.64429
## drift 0.0055550 0.0068802 0.8074 0.41945
## ldolar -3.7986183 0.4860777 -7.8148 0.000000000000005503 ***
## lip 3.1119236 0.4543242 6.8496 0.000000000007407434 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Para cada modelo calculamos todo lo que necesitamos para elegir (solo con train)
# Raíz inversa más grande de la parte AR y MA: debe ser < 1 (estacionario e invertible).
# Si da casi 1, hay una raíz AR y una MA que se cancelan: el modelo tiene parámetros que sobran
raiz_max <- function(m) {
p <- m$arma[1]; q <- m$arma[2] # cuántos AR (p) y cuántos MA (q) tiene
r <- c(if (p > 0) abs(1 / polyroot(c(1, -coef(m)[paste0("ar", 1:p)]))), # raíces inversas de la parte AR
if (q > 0) abs(1 / polyroot(c(1, coef(m)[paste0("ma", 1:q)])))) # raíces inversas de la parte MA
if (length(r) == 0) 0 else max(r) # nos quedamos con la más grande
}
resumen_modelo <- function(m) {
ct <- coeftest(m) # coeficientes y sus p-valores
arma <- grepl("^(ar|ma|sar|sma)", rownames(ct)) # cuáles son coeficientes AR o MA (no drift ni exógenas)
ee <- suppressWarnings(sqrt(diag(m$var.coef))) # errores estándar (si salen NaN, la estimación falló)
data.frame(AIC = m$aic, AICc = m$aicc, BIC = m$bic, # criterios de información: menor es mejor
k = length(coef(m)), # número de parámetros
arma_signif = all(ct[arma, 4] < 0.05), # ¿TODOS los AR y MA son significativos?
ee_validos = all(is.finite(ee)), # ¿los errores estándar se calcularon bien?
raiz_max = raiz_max(m), # raíz inversa más grande
LjungBox_p = Box.test(residuals(m), lag = 24, type = "Ljung-Box",
fitdf = sum(m$arma[1:4]))$p.value, # ¿quedó autocorrelación en los residuos?
JarqueBera_p = jarque.bera.test(na.omit(residuals(m)))$p.value) # ¿los residuos son normales?
}
tabla <- do.call(rbind, lapply(modelos, resumen_modelo)) # una fila por modelo
tabla <- data.frame(modelo = names(modelos), tabla, row.names = NULL) # agrega el nombre del modelo
knitr::kable(tabla[order(tabla$AICc), ], digits = 4, row.names = FALSE,
caption = "Modelos candidatos (train), ordenados por AICc") # ordenada del mejor al peor AICc| modelo | AIC | AICc | BIC | k | arma_signif | ee_validos | raiz_max | LjungBox_p | JarqueBera_p |
|---|---|---|---|---|---|---|---|---|---|
| ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | -405.1090 | -404.6692 | -385.3794 | 5 | TRUE | TRUE | 0.1719 | 0.4506 | 0 |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO | -400.0096 | -399.6971 | -383.5683 | 4 | TRUE | TRUE | 0.1629 | 0.3698 | 0 |
| SARIMAX(0,1,1)(0,0,1)[12] DTWEXBGS+INDPRO | -398.2231 | -397.7833 | -378.4935 | 5 | FALSE | TRUE | 0.1631 | 0.3225 | 0 |
| ARIMA(2,1,1) | -316.1053 | -315.7928 | -299.6639 | 4 | TRUE | TRUE | 1.0000 | 0.8043 | 0 |
| ARIMA(2,1,2) | -315.0240 | -314.5843 | -295.2944 | 5 | FALSE | TRUE | 1.0000 | 0.8745 | 0 |
| ARIMA(0,1,1) | -313.8011 | -313.6774 | -303.9363 | 2 | TRUE | TRUE | 0.3192 | 0.8767 | 0 |
| ARIMA(1,1,1) | -311.8853 | -311.6780 | -298.7322 | 3 | FALSE | TRUE | 0.2707 | 0.8309 | 0 |
| ARIMA(1,1,0) | -311.7793 | -311.6556 | -301.9145 | 2 | TRUE | TRUE | 0.2948 | 0.6631 | 0 |
| SARIMA(0,1,1)(0,0,1)[12] | -311.8367 | -311.6294 | -298.6836 | 3 | FALSE | TRUE | 0.3192 | 0.8398 | 0 |
| ARIMA(1,1,2) | -310.1066 | -309.7941 | -293.6653 | 4 | FALSE | TRUE | 0.3729 | 0.7750 | 0 |
Todos los modelos se estimaron con la función Arima del
paquete forecast (Hyndman et al., 2026), con órdenes propuestos a mano a
partir del correlograma.
El AR(1) y el MA(1) tienen su coeficiente claramente significativo (0,295 y 0,319, con p menor que 0,001). En ambos la deriva no es significativa (p = 0,87 y 0,85), lo que indica que el retorno mensual promedio no se distingue de cero.
En el ARMA(1,1) los dos coeficientes pierden significancia (p = 0,77 y 0,13). Los dos términos intentan explicar la misma autocorrelación del primer rezago y se reparten el efecto, de modo que ninguno queda bien estimado. En el ARMA(1,2) ocurre lo mismo con tres coeficientes.
El ARIMA(2,1,1) merece atención. Todos sus coeficientes ARMA son significativos y su AICc es mejor que el del MA(1), pero el coeficiente de media móvil es prácticamente −1 (−0,99999): la raíz del polinomio MA está sobre el círculo unitario y el modelo no es invertible. Lo mismo pasa en el ARIMA(2,1,2). Un modelo no invertible no es válido para pronosticar, y por eso ambos se descartan a pesar de su AICc. Hay además una lectura de fondo: una raíz unitaria en la parte MA cancela la diferencia aplicada, lo que equivale a decir que el logaritmo del precio sería estacionario alrededor de una tendencia. Eso contradice las pruebas de la sección 5.3 y es la señal de sobrediferenciación que se menciona en clase (Rodríguez Guevara, 2026b).
El término estacional no es significativo ni en el SARIMA (p = 0,85) ni en el SARIMAX (p = 0,64), lo que confirma la sección 5.5.
En los modelos con exógenas, el dólar y la producción industrial son muy significativos y tienen los signos esperados. En el modelo con las tres variables, los coeficientes son −3,55 para el dólar y 2,77 para la producción industrial. Al estar en logaritmos se leen como elasticidades: un aumento de 1 % en el índice del dólar se asocia con una caída cercana a 3,5 % en el WTI en el mismo mes. La tasa de la Fed también es significativa (0,115, p = 0,008) y al incluirla el AICc mejora de −399,70 a −404,67. Son magnitudes grandes, y es probable que reflejen sobre todo los meses de 2008 y 2020 en que todas las variables se movieron a la vez. Son asociaciones contemporáneas y no efectos causales.
Todos los modelos tienen \(d = 1\) y la misma variable dependiente, así que sus criterios de información son comparables entre sí (Rodríguez Guevara y González Uribe, 2017).
5.7 Selección del modelo
# Regla fijada ANTES de mirar el test (el enunciado prohíbe elegir con el test):
# coeficientes AR/MA significativos + errores válidos + raíz inversa < 0,98 + Ljung-Box p > 0,05 -> menor AICc
validos <- subset(tabla, arma_signif & ee_validos & raiz_max < 0.98 & LjungBox_p > 0.05) # modelos que pasan todo
if (nrow(validos) == 0) validos <- tabla # respaldo: si ninguno cumpliera, se usa la tabla completa
knitr::kable(validos[order(validos$AICc), c("modelo", "AICc", "BIC", "LjungBox_p")], digits = 4,
row.names = FALSE, caption = "Modelos que cumplen la regla")| modelo | AICc | BIC | LjungBox_p |
|---|---|---|---|
| ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | -404.6692 | -385.3794 | 0.4506 |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO | -399.6971 | -383.5683 | 0.3698 |
| ARIMA(0,1,1) | -313.6774 | -303.9363 | 0.8767 |
| ARIMA(1,1,0) | -311.6556 | -301.9145 | 0.6631 |
con_x <- names(specs)[sapply(specs, function(s) !is.null(s$x))] # nombres de los modelos con exógenas
MODELO_REGLA <- validos$modelo[which.min(validos$AICc)] # el de menor AICc entre todos
univariados <- validos[!validos$modelo %in% con_x, ] # solo los que NO usan exógenas
MODELO_UNIV <- univariados$modelo[which.min(univariados$AICc)] # el mejor univariado
con_exog <- validos[validos$modelo %in% con_x, ] # solo los que SÍ usan exógenas
if (nrow(con_exog) == 0) con_exog <- tabla[tabla$modelo %in% con_x, ] # respaldo si ninguno cumpliera
MODELO_X <- con_exog$modelo[which.min(con_exog$AICc)] # el mejor con exógenas
# MODELO PARA LA DECISIÓN (con lo que se sabía en julio de 2022):
# las exógenas entran en el MISMO mes que el WTI; en T0 no se conocen sus valores futuros,
# e INDPRO además se publica con un mes de rezago. Por eso decidimos con el univariado.
MODELO_FINAL <- MODELO_UNIV
# MODELO_FINAL <- MODELO_X # alternativa: solo si el equipo lo justifica con argumentos de train
# Si el drift (tendencia) no es significativo, el modelo se vuelve a estimar sin él (parsimonia)
sin_drift_si_aplica <- function(nm) {
m <- modelos[[nm]] # toma el modelo
if ("drift" %in% names(coef(m)) && coeftest(m)["drift", 4] > 0.05) m <- ajustar(specs[[nm]], drift = FALSE) # p > 0,05 -> sin drift
m # devuelve el modelo final
}
fit <- sin_drift_si_aplica(MODELO_FINAL) # modelo para la decisión
usa_X <- !is.null(specs[[MODELO_FINAL]]$x) # ¿el modelo final usa exógenas? (aquí: no)
x_vars <- specs[[MODELO_FINAL]]$x # cuáles exógenas usa (aquí: ninguna)
fit_x <- sin_drift_si_aplica(MODELO_X) # mejor modelo con exógenas, solo para comparar
x_vars_x <- specs[[MODELO_X]]$x # sus exógenas
cat("Menor AICc de toda la tabla:", MODELO_REGLA, "\n") # resumen de la selección## Menor AICc de toda la tabla: ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDScat("Mejor univariado: ", MODELO_UNIV, "\n")## Mejor univariado: ARIMA(0,1,1)cat("Mejor con exógenas: ", MODELO_X, "\n")## Mejor con exógenas: ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDScat("MODELO PARA LA DECISIÓN: ", MODELO_FINAL, "\n\n")## MODELO PARA LA DECISIÓN: ARIMA(0,1,1)print(fit) # coeficientes, varianza y criterios del modelo final## Series: train
## ARIMA(0,1,1)
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ma1
## 0.3194
## s.e. 0.0650
##
## sigma^2 = 0.0117: log likelihood = 159.88
## AIC=-315.76 AICc=-315.7 BIC=-309.19coeftest(fit) # significancia de sus coeficientes##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.319389 0.065008 4.9131 0.0000008967 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("\n"); print(fit_x) # lo mismo para el modelo con exógenas## Series: train
## Regression with ARIMA(0,1,1) errors
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ma1 ldolar lip ff
## 0.1758 -3.5232 2.7805 0.1128
## s.e. 0.0681 0.4856 0.4653 0.0432
##
## sigma^2 = 0.007303: log likelihood = 208.09
## AIC=-406.18 AICc=-405.87 BIC=-389.74coeftest(fit_x)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.175847 0.068133 2.5809 0.009853 **
## ldolar -3.523206 0.485562 -7.2559 0.0000000000003989 ***
## lip 2.780544 0.465295 5.9759 0.0000000022885996 ***
## ff 0.112785 0.043201 2.6107 0.009035 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La regla de selección se definió antes de mirar la prueba y tiene dos pasos. Primero, solo compiten los modelos que cumplen el diagnóstico (coeficientes ARMA significativos, errores estándar válidos, invertibilidad y residuos sin autocorrelación según Ljung-Box) y que se pueden usar con la información disponible en la fecha de corte. Segundo, entre ellos se escoge el de menor AICc.
Cuatro modelos cumplen el diagnóstico: los dos ARIMAX, el ARIMA(0,1,1) y el ARIMA(1,1,0). El de menor AICc es el ARIMAX con las tres exógenas (−404,67). Sin embargo, no cumple la condición de disponibilidad. Sus variables entran en el mismo mes que el precio, así que para proyectar el WTI de agosto de 2022 habría que conocer el dólar, la producción industrial y la tasa de la Fed de agosto de 2022. En julio no se conocían. La producción industrial, además, se publica a mediados del mes siguiente y luego se revisa.
Este caso muestra por qué un menor criterio de información no basta para elegir. El AICc mide qué tan bien explica el modelo el pasado, y el ARIMAX explica muy bien el pasado porque usa información del mismo mes. Para decidir se necesita un modelo que funcione con lo que se sabe en el momento de decidir. Lo mismo pasó con el ARIMA(2,1,1): mejor AICc, pero no invertible.
El modelo de decisión es entonces el mejor univariado, el ARIMA(0,1,1), que supera al ARIMA(1,1,0) por dos puntos de AICc. La diferencia es pequeña y los dos modelos son casi equivalentes; se sigue la regla. Como su deriva no es significativa (p = 0,849), se reestima sin ella, porque no hay razón para imponer una tendencia fija al retorno de un commodity. El modelo final es:
\[\Delta \ln P_t = \varepsilon_t + 0{,}3194\,\varepsilon_{t-1}, \qquad \hat\sigma^2 = 0{,}0117, \qquad AICc = -315{,}70\]
El coeficiente tiene un error estándar de 0,065 y es significativo (p menor que 0,001) bajo el supuesto de varianza constante; la sección 6.1 revisa este resultado. El ARIMAX se conserva solo como punto de comparación. Esta elección no se modifica después de ver la prueba.
6. Diagnóstico de residuos
res <- na.omit(residuals(fit)) # residuos: lo que el modelo NO pudo explicar cada mes
checkresiduals(fit, lag = 24) # diapositivas: gráfico de residuos, su ACF, histograma y prueba de Ljung-Box
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)
## Q* = 15.474, df = 23, p-value = 0.8768
##
## Model df: 1. Total lags used: 24cat("Media de los residuos:", round(mean(res), 5), "\n\n") # debe ser casi 0 (el modelo no se equivoca hacia un lado)## Media de los residuos: 0.00155k <- sum(fit$arma[1:4]) # parámetros AR + MA del modelo (para fitdf)
# Ljung-Box a 6, 12 y 24 rezagos: p > 0,05 -> los residuos son ruido blanco (el modelo capturó la memoria)
for (L in c(6, 12, 24)) print(Box.test(res, lag = L, type = "Ljung-Box", fitdf = k))##
## Box-Ljung test
##
## data: res
## X-squared = 3.3916, df = 5, p-value = 0.6398
##
##
## Box-Ljung test
##
## data: res
## X-squared = 4.272, df = 11, p-value = 0.9613
##
##
## Box-Ljung test
##
## data: res
## X-squared = 15.474, df = 23, p-value = 0.8768# Correlograma de los residuos con acfinter (función del profesor)
tab_res <- as.data.frame(acfinter(as.numeric(res), lag = 36, ci.method = "white")[["ACF-PACF Test"]])
banda_r <- 1.96 / sqrt(length(res)) # banda al 95 %
par(mfrow = c(3, 1), mar = c(4, 4, 2, 1)) # 3 gráficos en columna
# Si ninguna barra pasa las bandas rojas, los residuos no tienen memoria
plot(tab_res$lag, tab_res$acf, type = "h", lwd = 2, xlab = "rezago", ylab = "ACF", main = "ACF de los residuos")
abline(h = 0); abline(h = c(-1, 1) * banda_r, lty = 2, col = "red")
plot(tab_res$lag, tab_res$pacf, type = "h", lwd = 2, xlab = "rezago", ylab = "PACF", main = "PACF de los residuos")
abline(h = 0); abline(h = c(-1, 1) * banda_r, lty = 2, col = "red")
# Todos los puntos por encima de la línea del 5 % = ruido blanco
plot(tab_res$lag, tab_res$Pv_Ljung, pch = 19, ylim = c(0, 1), xlab = "rezago", ylab = "p-valor",
main = "p-valor de Ljung-Box (residuos)")
abline(h = 0.05, lty = 2, col = "red")
par(mfrow = c(1, 1)) # vuelve a un gráfico por ventanapar(mfrow = c(1, 2)) # 2 gráficos lado a lado
hist(res, breaks = 30, freq = FALSE, main = "Histograma de residuos", xlab = "") # forma de los errores
curve(dnorm(x, mean(res), sd(res)), add = TRUE, col = "red", lwd = 2) # campana normal para comparar
qqnorm(res, main = "Gráfico Q-Q"); qqline(res, col = "red", lwd = 2) # si los puntos se salen de la línea en las puntas: colas pesadas
par(mfrow = c(1, 1))
jb <- jarque.bera.test(res); print(jb) # Jarque-Bera (la que pide el enunciado). H0: normal. p < 0,05 -> no son normales##
## Jarque Bera Test
##
## data: res
## X-squared = 786.67, df = 2, p-value < 0.00000000000000022print(shapiro.test(as.numeric(res))) # Shapiro-Wilk como complemento (mismo H0)##
## Shapiro-Wilk normality test
##
## data: as.numeric(res)
## W = 0.88507, p-value = 0.00000000003359# ¿La varianza cambia en el tiempo? (diapositivas de ARCH-GARCH). Solo se PRUEBA; no se estima un ARCH
# El efecto ARCH se prueba únicamente con ArchTest (prueba ARCH-LM de las diapositivas), con 12 rezagos
print(ArchTest(res, lags = 12)) # ARCH-LM de Engle. H0: no hay efecto ARCH. p < 0,05 -> hay efecto ARCH##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: res
## Chi-squared = 51.104, df = 12, p-value = 0.0000008928# Mismo diagnóstico para el mejor modelo con exógenas, para comparar
res_x <- na.omit(residuals(fit_x)) # sus residuos
cat("\nModelo con exógenas:", MODELO_X, "\n")##
## Modelo con exógenas: ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDScheckresiduals(fit_x, lag = 24) # gráficos + Ljung-Box
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(0,1,1) errors
## Q* = 23.217, df = 23, p-value = 0.4482
##
## Model df: 1. Total lags used: 24print(jarque.bera.test(res_x)) # normalidad##
## Jarque Bera Test
##
## data: res_x
## X-squared = 228.69, df = 2, p-value < 0.00000000000000022print(ArchTest(res_x, lags = 12)) # efecto ARCH##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: res_x
## Chi-squared = 7.1748, df = 12, p-value = 0.8459Media. La media de los residuos es 0,0016 en escala logarítmica, es decir, 0,16 % mensual. El modelo no se equivoca de forma sistemática hacia arriba ni hacia abajo.
Autocorrelación. La prueba de Ljung-Box no rechaza la ausencia de autocorrelación en ningún horizonte: p = 0,640 con 6 rezagos, 0,961 con 12 y 0,877 con 24. En el correlograma ninguna barra relevante supera las bandas. El MA(1) recogió toda la dependencia lineal del retorno.
Normalidad. Jarque-Bera rechaza la normalidad con claridad (estadístico de 786,7, p menor que 0,001) y Shapiro-Wilk coincide (W = 0,885). El histograma y el gráfico Q-Q muestran la causa: los puntos se separan de la línea en los dos extremos. Son las colas pesadas que ya se veían en los descriptivos.
Varianza. La prueba ARCH-LM de Engle (1982) con 12 rezagos, aplicada como en las diapositivas (Rodríguez Guevara, 2026a), rechaza la hipótesis de que no hay efecto ARCH (estadístico de 51,1, p menor que 0,001). La varianza no es constante: los meses de alta volatilidad tienden a venir seguidos, como en 2008 y 2020.
En el ARIMAX de comparación los residuos tampoco tienen autocorrelación (p = 0,448) ni son normales, pero la prueba ARCH no rechaza (p = 0,846). Una explicación posible es que la producción industrial y el dólar absorben los meses extremos de 2008 y 2020, que son los que generan la agrupación de volatilidad en el modelo univariado. Se plantea como hipótesis; no se comprobó.
Respuesta. Los residuos cumplen de forma parcial. Tienen media cero y no están autocorrelacionados, que es lo que se necesita para confiar en el pronóstico puntual. No son normales ni tienen varianza constante. Estos incumplimientos no invalidan el pronóstico puntual, pero hacen menos precisa la estimación del coeficiente y menos confiable su error estándar, y afectan todo lo que depende de la forma de la distribución (Brooks, 2008). Para el pronóstico, los intervalos de predicción suponen errores normales de varianza fija y por tanto no son confiables: no distinguen entre periodos agitados y periodos tranquilos. Para la simulación, generar choques normales subestimaría los escenarios extremos; por eso la sección 8 usa bootstrap de residuos.
Como la heterocedasticidad es marcada, en la sección 6.1 se modela la varianza y se vuelven a aplicar los criterios de modelación.
6.1 Modelo para la varianza: ARCH y GARCH
# Modelo para la varianza, con el procedimiento de las diapositivas de ARCH-GARCH (paquete rugarch).
# Todo se estima SOLO con train. La ecuación de la media es la del modelo final: MA(1) sin constante
r_train <- 100 * as.numeric(dly) # retorno logarítmico mensual en %: r_t = 100 * ln(P_t / P_(t-1)), como en las diapositivas
# Candidatos para la varianza. OJO: rugarch recibe garchOrder = c(rezagos ARCH, rezagos GARCH)
ordenes_var <- list("ARCH(1)" = c(1, 0), # h_t = omega + alpha1 * e_(t-1)^2
"GARCH(1,1)" = c(1, 1)) # h_t = omega + alpha1 * e_(t-1)^2 + beta1 * h_(t-1)
# Estima cada candidato: misma muestra, misma media y misma distribución (normal), para que sean comparables
ajustes_var <- lapply(ordenes_var, function(ord) {
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = ord), # ecuación de la varianza
mean.model = list(armaOrder = c(0, 1), include.mean = FALSE), # ecuación de la media: MA(1) sin constante
distribution.model = "norm") # innovaciones normales
tryCatch(ugarchfit(spec = spec, data = r_train, solver = "hybrid"), error = function(e) NULL) # si falla, devuelve NULL
})
# Solo se conservan los modelos que se pudieron estimar y que convergieron (convergence == 0, como piden las diapositivas)
ajustes_var <- ajustes_var[!sapply(ajustes_var, is.null)]
ajustes_var <- ajustes_var[sapply(ajustes_var, function(m) m@fit$convergence == 0)]
cat("Modelos de varianza estimados:", length(ajustes_var), "de", length(ordenes_var), "\n")## Modelos de varianza estimados: 2 de 2# Coeficientes de cada candidato con su error estándar y su p-valor (Pr(>|t|) < 0,05 = significativo)
for (nm in names(ajustes_var)) {
cat("\n==========", nm, "==========\n") # título de cada modelo
print(round(ajustes_var[[nm]]@fit$matcoef, 4)) # ma1 = media; omega, alpha1 y beta1 = varianza
}##
## ========== ARCH(1) ==========
## Estimate Std. Error t value Pr(>|t|)
## ma1 0.1157 0.0849 1.3635 0.1727
## omega 55.0405 8.2035 6.7094 0.0000
## alpha1 0.4885 0.1450 3.3687 0.0008
##
## ========== GARCH(1,1) ==========
## Estimate Std. Error t value Pr(>|t|)
## ma1 0.1475 0.0810 1.8214 0.0685
## omega 36.8947 12.7165 2.9013 0.0037
## alpha1 0.4631 0.1373 3.3742 0.0007
## beta1 0.2152 0.1558 1.3814 0.1672# Para cada candidato: criterios de información, condiciones de los parámetros y diagnóstico de residuos estandarizados
resumen_var <- function(m) {
z <- as.numeric(residuals(m, standardize = TRUE)) # residuos estandarizados: z_t = e_t / raíz de h_t
co <- m@fit$matcoef # tabla de coeficientes
es_var <- grepl("^(alpha|beta)", rownames(co)) # cuáles son coeficientes de la varianza
data.frame(AIC = infocriteria(m)[1], BIC = infocriteria(m)[2], # criterios de información: menor es mejor
persistencia = sum(co[es_var, 1]), # alpha1 + beta1: debe ser < 1 para que la varianza de largo plazo exista
var_signif = all(co[es_var, 4] < 0.05), # ¿TODOS los coeficientes de la varianza son significativos?
LjungBox_p = Box.test(z, lag = 12, type = "Ljung-Box", fitdf = 1)$p.value, # ¿quedó autocorrelación en la media?
ArchTest_p = ArchTest(z, lags = 12)$p.value) # ¿quedó efecto ARCH? Se quiere p > 0,05
}
tabla_var <- do.call(rbind, lapply(ajustes_var, resumen_var)) # una fila por modelo
tabla_var <- data.frame(modelo = names(ajustes_var), tabla_var, row.names = NULL) # agrega el nombre del modelo
knitr::kable(tabla_var[order(tabla_var$AIC), ], digits = 4, row.names = FALSE,
caption = "Modelos candidatos para la varianza (train), ordenados por AIC")| modelo | AIC | BIC | persistencia | var_signif | LjungBox_p | ArchTest_p |
|---|---|---|---|---|---|---|
| ARCH(1) | 7.3070 | 7.3568 | 0.4885 | TRUE | 0.309 | 0.9636 |
| GARCH(1,1) | 7.3085 | 7.3749 | 0.6783 | FALSE | 0.406 | 0.9907 |
# Misma lógica de la sección 5.7: primero deben cumplir el diagnóstico y después se escoge el de menor AIC
validos_var <- subset(tabla_var, var_signif & persistencia < 1 & LjungBox_p > 0.05 & ArchTest_p > 0.05) # modelos que pasan todo
if (nrow(validos_var) == 0) validos_var <- tabla_var # respaldo: si ninguno cumpliera, se usa la tabla completa
MODELO_VAR <- validos_var$modelo[which.min(validos_var$AIC)] # modelo elegido para la varianza
fit_var <- ajustes_var[[MODELO_VAR]] # su ajuste
cat("MODELO PARA LA VARIANZA:", MODELO_VAR, "\n")## MODELO PARA LA VARIANZA: ARCH(1)# Volatilidad condicional estimada mes a mes (raíz de h_t), en % mensual
vol_cond <- as.numeric(sigma(fit_var))
plot(fechas[2:n_train], vol_cond, type = "l", xlab = "", ylab = "% mensual",
main = paste("Volatilidad condicional estimada:", MODELO_VAR))
abline(h = sd(r_train), col = "red", lty = 2) # línea roja: desviación estándar del retorno en todo el train (volatilidad constante)
# Pronóstico de la volatilidad desde T0, como en las diapositivas (ugarchforecast)
fc_var <- ugarchforecast(fit_var, n.ahead = H_DEC) # 12 meses hacia adelante
vol_fc <- as.numeric(sigma(fc_var)) # volatilidad pronosticada (desviación estándar, % mensual)
co_var <- fit_var@fit$matcoef # coeficientes del modelo elegido
persist <- sum(co_var[grepl("^(alpha|beta)", rownames(co_var)), 1]) # alpha1 + beta1
vol_lp <- sqrt(co_var["omega", 1] / (1 - persist)) # volatilidad de largo plazo: raíz de omega / (1 - persistencia)
data.frame(h = c(1, 3, 6, 12), volatilidad = round(vol_fc[c(1, 3, 6, 12)], 2), largo_plazo = round(vol_lp, 2))## h volatilidad largo_plazo
## 1 1 11.55 10.37
## 2 3 10.67 10.37
## 3 6 10.41 10.37
## 4 12 10.37 10.37# Comprobación: ¿cambiaría la decisión si la media se pronostica con el coeficiente MA estimado junto con la varianza?
e_T_var <- tail(as.numeric(residuals(fit_var)), 1) # último residuo del modelo (julio 2022), en %
P_hat_var <- P0 * exp(co_var["ma1", 1] * e_T_var / 100) # precio proyectado: P0 * exp(theta * último residuo)
ret_var <- P_hat_var / P0 - 1 # retorno proyectado frente al precio de entrada
cat("Precio proyectado con la media de", MODELO_VAR, ":", round(P_hat_var, 2), "USD | Retorno:", round(ret_var * 100, 2),
"% | Decisión:", ifelse(ret_var > UMBRAL, "LARGO", "LIQUIDEZ"), "\n") # misma regla de la sección 8## Precio proyectado con la media de ARCH(1) : 100.14 USD | Retorno: -1.46 % | Decisión: LIQUIDEZComo la heterocedasticidad es marcada, se reaplicaron a la varianza los mismos criterios de modelación usados para la media, siguiendo el procedimiento de las diapositivas (Rodríguez Guevara, 2026a): se conserva la ecuación de la media, se proponen candidatos para la varianza, se estiman solo con el entrenamiento, se comparan con AIC y BIC, se revisan las condiciones de los parámetros y se diagnostican los residuos estandarizados. Los modelos se estimaron con el paquete rugarch (Galanos, 2026), con el retorno expresado en porcentaje.
Se propusieron dos candidatos, ambos con la media MA(1) sin constante del modelo final. El ARCH(1) de Engle (1982), en el que la varianza de un mes depende del choque al cuadrado del mes anterior, y el GARCH(1,1) de Bollerslev (1986), que agrega la varianza del mes anterior.
El ARCH(1) tiene el menor valor tanto en AIC (7,3070 frente a 7,3085) como en BIC (7,3568 frente a 7,3749). En el GARCH(1,1) el coeficiente de la varianza rezagada es 0,215 y no es significativo (p = 0,167): la varianza del mes anterior no aporta más de lo que ya explica el choque reciente. Aplicando la misma lógica de la sección 5.7 (coeficientes de la varianza significativos, persistencia menor que 1, residuos estandarizados sin autocorrelación y sin efecto ARCH, y entre los que cumplen el de menor AIC), el modelo elegido para la varianza es el ARCH(1):
\[h_t = 55{,}04 + 0{,}489\,\varepsilon_{t-1}^2\]
El coeficiente del choque al cuadrado, \(\alpha_1 = 0{,}489\), es significativo (p menor que 0,001): cerca de la mitad de un choque grande se traslada a la varianza del mes siguiente. Como es menor que 1, la varianza de largo plazo existe y equivale a una volatilidad de 10,4 % mensual, \(\sqrt{55{,}04/(1-0{,}489)}\).
El diagnóstico se hace sobre los residuos estandarizados, es decir, los residuos divididos por la volatilidad estimada de cada mes. La prueba ARCH-LM con 12 rezagos ya no rechaza la hipótesis de que no hay efecto ARCH (p = 0,964): el modelo recogió la agrupación de volatilidad. Los residuos estandarizados tampoco muestran autocorrelación según Ljung-Box (p = 0,309).
Hay un resultado que debe reportarse porque afecta la lectura del modelo de la media. Al estimar la media y la varianza en conjunto, el coeficiente de media móvil baja de 0,319 a 0,116 y deja de ser significativo (p = 0,173). Al modelar la varianza, los meses de volatilidad extrema pesan menos en la estimación y la autocorrelación del retorno se debilita. Esto indica que buena parte de la memoria que captura el MA(1) proviene de los episodios de 2008 y 2020, y que la significancia reportada en la sección 5.7 dependía del supuesto de varianza constante.
Este resultado no cambia la decisión de inversión. Con el coeficiente de 0,116, el precio proyectado sería de 100,14 USD, un retorno de −1,46 % frente al precio de entrada, también por debajo del umbral de 5 %. Con cualquiera de las dos estimaciones la regla indica liquidez.
El gráfico muestra que la volatilidad no es constante: se dispara en 2008 y en 2020 y permanece baja en los periodos tranquilos. Para agosto de 2022, el primer mes después de la fecha de corte, el modelo proyecta una volatilidad de 11,55 % mensual, por encima del promedio de largo plazo porque julio fue un mes de caída fuerte, y el pronóstico converge a 10,37 % a medida que crece el horizonte.
Implicaciones para el pronóstico y la simulación. Con varianza condicional, la incertidumbre de un pronóstico depende del momento en que se hace: después de un choque grande los intervalos deberían ser más anchos, y en periodos tranquilos más angostos. Los intervalos de la sección 7 no recogen esto porque suponen una varianza fija. En la simulación, los choques deberían generarse con la volatilidad proyectada mes a mes, de modo que los meses extremos aparezcan agrupados como ocurre en los datos.
Alcance. El modelo de varianza se usa para describir y diagnosticar la volatilidad. El modelo de decisión sigue siendo el ARIMA(0,1,1) seleccionado en la sección 5.7 con la regla fijada de antemano, y los intervalos de la sección 7 y la simulación de la sección 8 no se recalcularon con la varianza condicional. Esto se reconoce como limitación en la sección 9.4.
7. Train-test y capacidad predictiva
# Las exógenas del futuro NO se conocen en julio de 2022. Las proyectamos con una caminata
# aleatoria con deriva (siguen su tendencia promedio), usando SOLO el train: así no hay data leakage
X_esc <- sapply(colnames(X_train), function(v) as.numeric(rwf(X_train[, v], h = H, drift = TRUE)$mean))
# Tabla con el escenario supuesto para las exógenas: así quedan documentados los valores futuros que usa el pronóstico.
# Se imprimen los 12 meses de la decisión; los valores están en la escala del modelo (log del dólar, log de INDPRO, tasa en %)
knitr::kable(data.frame(fecha = fechas_test[1:H_DEC], round(X_esc[1:H_DEC, ], 4)),
caption = "Escenario de las exógenas para el pronóstico (caminata con deriva, solo con train)")| fecha | ldolar | lip | ff |
|---|---|---|---|
| 2022-08-01 | 4.8112 | 4.6175 | 1.6668 |
| 2022-09-01 | 4.8122 | 4.6176 | 1.6536 |
| 2022-10-01 | 4.8133 | 4.6178 | 1.6405 |
| 2022-11-01 | 4.8143 | 4.6179 | 1.6273 |
| 2022-12-01 | 4.8153 | 4.6181 | 1.6141 |
| 2023-01-01 | 4.8164 | 4.6182 | 1.6009 |
| 2023-02-01 | 4.8174 | 4.6184 | 1.5877 |
| 2023-03-01 | 4.8184 | 4.6186 | 1.5745 |
| 2023-04-01 | 4.8195 | 4.6187 | 1.5614 |
| 2023-05-01 | 4.8205 | 4.6189 | 1.5482 |
| 2023-06-01 | 4.8215 | 4.6190 | 1.5350 |
| 2023-07-01 | 4.8226 | 4.6192 | 1.5218 |
# Función que pronostica 50 meses; si el modelo usa exógenas, le pasa el escenario proyectado
pronosticar <- function(m, vars) {
if (is.null(vars)) forecast(m, h = H) else forecast(m, xreg = X_esc[, vars, drop = FALSE])
}
pronosticos <- lapply(names(modelos), function(nm) pronosticar(modelos[[nm]], specs[[nm]]$x)) # pronóstico de cada candidato
names(pronosticos) <- names(modelos) # con su nombre
pronosticos[[paste("FINAL:", MODELO_FINAL)]] <- pronosticar(fit, x_vars) # el modelo final (sin drift)
pronosticos[[paste("CON EXÓGENAS:", MODELO_X)]] <- pronosticar(fit_x, x_vars_x) # el de exógenas (sin drift)
pronosticos[["Benchmark: caminata aleatoria"]] <- naive(train, h = H) # el precio no cambia (último precio)
pronosticos[["Benchmark: caminata aleatoria con deriva"]] <- rwf(train, h = H, drift = TRUE) # sigue la tendencia promedio
# accuracy() (diapositivas): compara cada pronóstico con lo que pasó de verdad en el test
# MASE con la fórmula de las diapositivas: MAE del test dividido por el error del método ingenuo a un paso
# dentro del train, es decir, el promedio de |y_t - y_(t-1)|. (accuracy() usaría el ingenuo estacional de
# 12 meses por ser una serie mensual; por eso el MASE se calcula aquí directamente)
escala_mase <- mean(abs(diff(as.numeric(train)))) # denominador del MASE (USD): cambio absoluto promedio de un mes al siguiente en train
metricas <- do.call(rbind, lapply(names(pronosticos), function(nm) {
a <- accuracy(pronosticos[[nm]], test) # errores en train y en test
data.frame(modelo = nm, RMSE_train = a["Training set", "RMSE"], # error en train (qué tan bien ajusta)
RMSE = a["Test set", "RMSE"], MAE = a["Test set", "MAE"], # errores en test, en USD por barril
MAPE = a["Test set", "MAPE"], # error en %
MASE = a["Test set", "MAE"] / escala_mase) # MAE del test / ingenuo a un paso en train (fórmula de las diapositivas)
}))
# En la tabla, los candidatos se estimaron CON deriva; "FINAL" y "CON EXÓGENAS" son los reestimados SIN deriva (no era significativa)
knitr::kable(metricas[order(metricas$RMSE), ], digits = 3, row.names = FALSE,
caption = "Desempeño fuera de muestra sobre el PRECIO (USD), ordenado por RMSE del test. Candidatos con deriva; FINAL y CON EXÓGENAS sin deriva") # del mejor al peor| modelo | RMSE_train | RMSE | MAE | MAPE | MASE |
|---|---|---|---|---|---|
| ARIMA(2,1,2) | 6.043 | 16.107 | 11.813 | 14.378 | 2.422 |
| CON EXÓGENAS: ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | 5.298 | 16.558 | 15.411 | 21.135 | 3.159 |
| ARIMA(2,1,1) | 6.018 | 17.205 | 13.207 | 16.215 | 2.708 |
| FINAL: ARIMA(0,1,1) | 6.160 | 23.496 | 21.358 | 30.154 | 4.379 |
| ARIMA(1,1,0) | 6.120 | 26.970 | 24.622 | 34.693 | 5.048 |
| Benchmark: caminata aleatoria | 6.543 | 27.153 | 25.102 | 35.228 | 5.146 |
| ARIMA(1,1,2) | 6.132 | 27.246 | 24.889 | 35.059 | 5.103 |
| ARIMA(1,1,1) | 6.150 | 27.782 | 25.466 | 35.829 | 5.221 |
| SARIMA(0,1,1)(0,0,1)[12] | 6.170 | 27.859 | 25.559 | 35.949 | 5.240 |
| ARIMA(0,1,1) | 6.166 | 28.254 | 25.964 | 36.496 | 5.323 |
| SARIMAX(0,1,1)(0,0,1)[12] DTWEXBGS+INDPRO | 5.305 | 30.099 | 28.005 | 39.194 | 5.742 |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | 5.297 | 30.767 | 28.637 | 40.063 | 5.871 |
| Benchmark: caminata aleatoria con deriva | 6.540 | 31.732 | 29.735 | 41.507 | 6.096 |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO | 5.297 | 31.921 | 29.850 | 41.683 | 6.120 |
fc_final <- pronosticos[[paste("FINAL:", MODELO_FINAL)]] # pronóstico del modelo final
fc_x <- pronosticos[[paste("CON EXÓGENAS:", MODELO_X)]] # pronóstico del modelo con exógenas
# Gráfico: pronóstico con sus intervalos, contra lo que realmente pasó y contra la caminata aleatoria
autoplot(fc_final) +
autolayer(test, series = "Observado") +
autolayer(fc_x$mean, series = "Con exógenas") +
autolayer(pronosticos[["Benchmark: caminata aleatoria"]]$mean, series = "Caminata aleatoria") +
labs(title = paste("Pronóstico vs. observado:", MODELO_FINAL), x = NULL, y = "USD/barril") +
theme_minimal()
# Cobertura: % de meses reales que cayeron dentro de los intervalos del 80 % y del 95 %
cobertura <- function(fc) c(IC80 = mean(test >= fc$lower[, 1] & test <= fc$upper[, 1]),
IC95 = mean(test >= fc$lower[, 2] & test <= fc$upper[, 2])) * 100
round(rbind(final = cobertura(fc_final), con_exogenas = cobertura(fc_x)), 1)## IC80 IC95
## final 100 100
## con_exogenas 100 100# Qué tan ancho es el intervalo: a 12 meses (la decisión) y a 50 (el test). Más ancho = más incertidumbre
cat("IC 95 % del precio a", H_DEC, "meses:", round(fc_final$lower[H_DEC, 2], 2), "a",
round(fc_final$upper[H_DEC, 2], 2), "USD\n")## IC 95 % del precio a 12 meses: 37.69 a 252.77 USDcat("IC 95 % del precio a", H, "meses:", round(fc_final$lower[H, 2], 2), "a",
round(fc_final$upper[H, 2], 2), "USD\n")## IC 95 % del precio a 50 meses: 13.62 a 699.41 USD# SOLO PARA DISCUTIR EL DATA LEAKAGE: ¿qué pasaría si "hiciéramos trampa" usando las exógenas REALES del test?
# Cada modelo se compara consigo mismo: misma especificación (sin deriva si no es significativa),
# una vez con las X proyectadas solo con train y otra con las X reales del test
for (nm in names(modelos)[sapply(names(modelos), function(nm) !is.null(specs[[nm]]$x))]) { # solo modelos con exógenas
m_or <- sin_drift_si_aplica(nm) # el mismo criterio de deriva que el modelo de la tabla
vars_or <- specs[[nm]]$x # las exógenas que usa este modelo
fc_proy <- forecast(m_or, xreg = X_esc[, vars_or, drop = FALSE]) # pronóstico honesto: X proyectadas solo con train
fc_trampa <- forecast(m_or, xreg = X_test[, vars_or, drop = FALSE]) # pronóstico "con trampa": X reales del futuro
cat("\n", nm, "\n") # nombre del modelo
print(round(rbind(X_proyectadas = accuracy(fc_proy, test)["Test set", c("RMSE", "MAE", "MAPE")], # errores sin trampa
X_reales = accuracy(fc_trampa, test)["Test set", c("RMSE", "MAE", "MAPE")]), 3)) # errores con trampa
}##
## ARIMAX(0,1,1) DTWEXBGS+INDPRO
## RMSE MAE MAPE
## X_proyectadas 19.219 17.830 24.763
## X_reales 30.423 27.418 38.145
##
## ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS
## RMSE MAE MAPE
## X_proyectadas 16.558 15.411 21.135
## X_reales 66.584 63.912 86.782
##
## SARIMAX(0,1,1)(0,0,1)[12] DTWEXBGS+INDPRO
## RMSE MAE MAPE
## X_proyectadas 17.939 16.579 22.983
## X_reales 28.750 25.712 35.847# La tasa de la Fed subió muchísimo en el test: ayuda a explicar por qué la "trampa" sale tan mal
cat("\nFEDFUNDS: en T0 =", X_train[n_train, "ff"], "| primer mes del test =", X_test[1, "ff"], "| máximo del test =", max(X_test[, "ff"]),
"| último mes =", X_test[H, "ff"], "\n") # tasa de la Fed en T0 (julio 2022), al inicio del test, en su máximo y al final##
## FEDFUNDS: en T0 = 1.68 | primer mes del test = 2.33 | máximo del test = 5.33 | último mes = 3.75Las métricas se calculan sobre el precio en dólares, comparando el pronóstico a 50 meses con lo observado. RMSE y MAE están en USD por barril; el RMSE castiga más los errores grandes. El MAPE es el error porcentual promedio. El MASE no tiene unidades. En las cuatro, un valor menor es mejor.
El modelo final tuvo un RMSE de 23,50 USD y un MAE de 21,36 USD, con un MAPE de 30,2 %. Supera a la caminata aleatoria (RMSE 27,15) y a la caminata con deriva (31,73), que son los modelos de referencia. La ventaja debe entenderse bien. Tanto el modelo como la caminata aleatoria proyectan una línea horizontal; la del modelo está en 97,61 USD y la de la caminata en 101,62 USD. Como el precio promedio del periodo fue 76,54 USD, la línea más baja quedó más cerca. La diferencia de 3,66 USD viene del ajuste del primer mes y no de una capacidad de anticipar la trayectoria. Un error medio de 21 USD sobre precios entre 58 y 102 USD es grande.
El MASE se calculó con la fórmula de las diapositivas: el MAE de la prueba dividido por el error absoluto medio del método ingenuo a un paso dentro del entrenamiento, que es de 4,88 USD (Hyndman y Koehler, 2006; Rodríguez Guevara, 2026b). El modelo final obtiene 4,38 y la caminata aleatoria 5,15. Que sean mayores que 1 es esperable: se está comparando un error a 50 meses con el cambio típico de un solo mes. Lo que informa es la comparación entre modelos, y ahí el orden es el mismo que con el MAE.
El MAPE tiene una limitación conocida: se dispara cuando el valor observado se acerca a cero y trata distinto los errores por encima y por debajo (Hyndman y Koehler, 2006). Aquí no es un problema grave porque el precio mínimo de la prueba fue 57,97 USD.
Ajuste frente a generalización. El modelo que mejor ajusta en el entrenamiento es el ARIMAX: tiene el menor RMSE dentro de muestra (5,30 frente a 6,16 del modelo final) y el menor AICc. En la prueba su desempeño depende por completo de dos supuestos. Con deriva es de los peores (RMSE 30,77); sin deriva es el segundo mejor (16,56). Lo que cambia entre una versión y otra no es la información de las exógenas, sino la tendencia que se le impone al pronóstico y el escenario supuesto para las X. Ajustar bien el pasado no garantiza pronosticar bien.
Los dos modelos con menor error en la prueba, ARIMA(2,1,2) y ARIMA(2,1,1), son los que se descartaron por no ser invertibles. Su buen resultado tiene explicación: esos modelos implican que el precio vuelve hacia un nivel medio, y eso fue lo que ocurrió después de julio de 2022. No se rescatan. Hacerlo sería escoger el modelo con la prueba, y un único periodo favorable no corrige una especificación inválida.
¿Confirmó la prueba la elección? En parte. El modelo elegido superó a las referencias, pero no fue el mejor de la tabla. La elección se mantiene porque se hizo con los criterios del entrenamiento, y cambiarla ahora convertiría la prueba en un mecanismo de selección.
Intervalos. Todos los precios observados cayeron dentro de los intervalos del 80 % y del 95 %. No es un buen resultado: un intervalo del 80 % que acierta el 100 % de las veces es más ancho de lo que debería. El intervalo del 95 % a doce meses va de 37,69 a 252,77 USD y a cincuenta meses de 13,62 a 699,41 USD. Con esa amplitud no ayudan a decidir. Son tan anchos sobre todo porque, al tener raíz unitaria, la incertidumbre se acumula mes a mes. Además se calculan suponiendo errores normales de varianza constante, y la sección 6 rechazó ambos supuestos. El modelo de varianza de la sección 6.1 no se usó para recalcular estos intervalos.
Exógenas y fuga de información. Para pronosticar con el ARIMAX se necesitan los valores futuros de las X. En la fecha de corte no existen, así que se proyectaron con una caminata aleatoria con deriva calculada solo con el entrenamiento: cada variable sigue su tendencia promedio de 2006 a 2022. La tabla del escenario muestra los valores usados. Usar los valores reales de la prueba sería fuga de información.
Como ejercicio, se repitió el pronóstico de cada modelo con los valores reales de las exógenas. El resultado empeora en todos los casos. En el ARIMAX con las tres variables el RMSE pasa de 16,56 a 66,58 USD, y en el ARIMAX con dólar y producción industrial de 19,22 a 30,42 USD. La causa principal es la tasa de la Fed: subió de 1,68 % en julio de 2022 a 5,33 % en 2023 mientras el petróleo bajaba, lo contrario de la relación positiva estimada entre 2006 y 2022. Con un coeficiente de 0,115, esa subida de 3,65 puntos empuja el pronóstico cerca de 50 % hacia arriba. La relación no fue estable fuera de muestra, lo que confirma que los \(\beta\) miden una asociación histórica y no un efecto causal.
Forma del pronóstico. Como se dedujo en la sección 4, un MA(1) sin deriva solo usa el último choque. Su pronóstico se ajusta en el primer mes y queda plano. Por eso se parece tanto a la caminata aleatoria.
Alcance. Estas métricas evalúan un pronóstico a 50 pasos desde un único origen. No miden el desempeño de corto plazo, que es donde un MA(1) podría aportar.
8. Pronóstico e inversión
8.1 Regla de decisión
# Decisión a H_DEC = 12 meses (julio 2022 -> julio 2023)
P_hat <- as.numeric(fc_final$mean[H_DEC]) # precio pronosticado para el mes 12 (es la MEDIANA)
ret_esp <- P_hat / P0 - 1 # cuánto subiría o bajaría frente al precio de entrada
decision <- ifelse(ret_esp > UMBRAL, "LARGO", "LIQUIDEZ") # regla: si sube más de 5 % compramos; si no, liquidez
barriles <- ifelse(decision == "LARGO", CAPITAL / P0, 0) # barriles que compraríamos con los USD 100.000
cat("Fecha de corte:", format(T0), "| Precio de entrada:", P0, "USD/barril\n")## Fecha de corte: 2022-07-01 | Precio de entrada: 101.62 USD/barrilcat("Precio pronosticado (mediana) a", H_DEC, "meses:", round(P_hat, 2), "USD/barril\n")## Precio pronosticado (mediana) a 12 meses: 97.61 USD/barrilcat("Retorno mediano proyectado:", round(ret_esp * 100, 2), "% | Umbral:", UMBRAL * 100, "%\n")## Retorno mediano proyectado: -3.95 % | Umbral: 5 %cat("DECISIÓN:", decision, "| Barriles:", round(barriles, 2), "\n\n")## DECISIÓN: LIQUIDEZ | Barriles: 0# Con lambda = 0, forecast() entrega la MEDIANA del precio. La MEDIA se obtiene con biasadj = TRUE.
# La media siempre es mayor porque el precio no baja de 0 pero puede subir mucho (cola derecha larga)
media_H <- function(m, vars) {
f <- if (is.null(vars)) forecast(m, h = H, biasadj = TRUE) else
forecast(m, xreg = X_esc[, vars, drop = FALSE], biasadj = TRUE) # pronóstico de la media
as.numeric(f$mean[H_DEC]) # valor en el mes 12
}
# Sensibilidad: ¿qué habríamos decidido con otro modelo o con la media en vez de la mediana?
sens <- data.frame(modelo = c(MODELO_FINAL, MODELO_FINAL, MODELO_X, MODELO_X, "Caminata aleatoria con deriva"),
medida = c("mediana (regla)", "media", "mediana", "media", "punto"),
precio_12 = c(P_hat, media_H(fit, x_vars), fc_x$mean[H_DEC], media_H(fit_x, x_vars_x),
pronosticos[["Benchmark: caminata aleatoria con deriva"]]$mean[H_DEC]))
sens$retorno_pct <- (sens$precio_12 / P0 - 1) * 100 # retorno de cada alternativa
sens$decision <- ifelse(sens$retorno_pct / 100 > UMBRAL, "LARGO", "LIQUIDEZ") # la decisión que daría
knitr::kable(sens, digits = 2, caption = paste("Decisión según cada pronóstico a", H_DEC, "meses"))| modelo | medida | precio_12 | retorno_pct | decision |
|---|---|---|---|---|
| ARIMA(0,1,1) | mediana (regla) | 97.61 | -3.95 | LIQUIDEZ |
| ARIMA(0,1,1) | media | 109.11 | 7.37 | LARGO |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | mediana | 94.35 | -7.16 | LIQUIDEZ |
| ARIMAX(0,1,1) DTWEXBGS+INDPRO+FEDFUNDS | media | 99.93 | -1.66 | LIQUIDEZ |
| Caminata aleatoria con deriva | punto | 103.81 | 2.15 | LIQUIDEZ |
8.2 Evaluación financiera
precios_test <- as.numeric(test) # precios reales del test
precios_dec <- precios_test[1:H_DEC] # solo los 12 meses de la decisión
fechas_dec <- fechas_test[1:H_DEC] # sus fechas
valor_modelo <- if (decision == "LARGO") barriles * precios_dec else rep(CAPITAL, H_DEC) # valor de nuestra estrategia mes a mes
valor_bh <- CAPITAL / P0 * precios_dec # comprar y mantener: lo que habría pasado si comprábamos sí o sí
caida_max <- function(v) max(1 - v / cummax(v)) # caída máxima: la peor pérdida desde un máximo (drawdown)
# Tabla: cómo le fue a cada estrategia con los precios reales
eval_fin <- data.frame(
estrategia = c(paste("Regla del modelo:", decision), "Comprar y mantener", "Liquidez"),
valor_final = c(tail(valor_modelo, 1), tail(valor_bh, 1), CAPITAL), # plata al final
retorno_pct = (c(tail(valor_modelo, 1), tail(valor_bh, 1), CAPITAL) / CAPITAL - 1) * 100, # ganancia o pérdida en %
# La caída máxima se mide desde el capital inicial (USD 100.000 en T0); por eso se antepone CAPITAL a cada serie
caida_max_pct = c(caida_max(c(CAPITAL, valor_modelo)), caida_max(c(CAPITAL, valor_bh)), 0) * 100) # peor caída en %, incluyendo el punto de entrada
knitr::kable(eval_fin, digits = 2, caption = paste("Resultado real a", H_DEC, "meses (precio final:",
round(precios_dec[H_DEC], 2), "USD)"))| estrategia | valor_final | retorno_pct | caida_max_pct |
|---|---|---|---|
| Regla del modelo: LIQUIDEZ | 100000.00 | 0.00 | 0.00 |
| Comprar y mantener | 74857.31 | -25.14 | 30.87 |
| Liquidez | 100000.00 | 0.00 | 0.00 |
# Gráfico del valor del portafolio en los 12 meses
plot(fechas_dec, valor_bh, type = "l", lwd = 2, col = "grey40", xlab = "", ylab = "USD",
main = paste("Valor del portafolio en los", H_DEC, "meses de la decisión"), ylim = range(valor_bh, valor_modelo))
lines(fechas_dec, valor_modelo, col = "red", lwd = 2) # nuestra estrategia en rojo
abline(h = CAPITAL, lty = 2) # línea de los USD 100.000 iniciales
legend("topleft", c("Comprar y mantener", "Regla del modelo"), col = c("grey40", "red"), lwd = 2, bty = "n")
En la fecha de corte el precio de entrada es 101,62 USD por barril. El modelo proyecta un precio mediano de 97,61 USD a doce meses, lo que equivale a un retorno mediano de −3,95 %. Como no supera el umbral de 5 %, la regla indica mantener los USD 100.000 en liquidez. De haber comprado, la posición habría sido de 984 barriles.
Conviene entender de dónde sale el 97,61. Julio de 2022 fue un mes de caída: el residuo del modelo fue cercano a −12,6 %. El MA(1) traslada 0,3194 de ese choque al mes siguiente, y \(101{,}62 \times e^{0{,}3194 \times (-0{,}126)} = 97{,}61\). A partir de ahí el pronóstico no cambia, de modo que el valor a doce meses es el mismo que a un mes. La decisión se apoya en un solo dato, la caída de julio. El modelo no contiene información propia sobre lo que pasará dentro de un año.
Esto tiene una consecuencia sobre la regla. Para que la mediana proyectada supere el 5 % se necesita que el último mes haya tenido un choque positivo mayor que 15 %, algo poco frecuente en la muestra de entrenamiento. Con este modelo y este umbral, la regla recomienda liquidez casi siempre.
Mediana y media. Cuando se modela el logaritmo, el pronóstico devuelto a la escala original es la mediana del precio y no su media (Hyndman y Athanasopoulos, 2021). La media del mismo pronóstico es 109,11 USD, un 7,37 % por encima del precio de entrada, y con ella la regla habría indicado comprar. La decisión cambia según la medida. La mediana se fijó en la hipótesis por la razón expuesta en la sección 2.3: representa el escenario típico, mientras que la media queda cerca de 12 % por encima por el peso de unos pocos escenarios muy altos. Pero hay que reconocer que el resultado depende de esa elección.
Sensibilidad. Con el ARIMAX la decisión habría sido la misma con las dos medidas (−7,16 % con la mediana y −1,66 % con la media), y también con la caminata aleatoria con deriva (+2,15 %). De las cinco alternativas, solo la media del modelo final lleva a comprar. Tampoco cambia con el coeficiente de media móvil estimado junto con la varianza en la sección 6.1 (−1,46 %).
Resultado. Doce meses después, en julio de 2023, el precio promedio fue 76,07 USD. La estrategia de liquidez conservó los USD 100.000. Comprar y mantener habría terminado en USD 74.857, una pérdida de 25,14 %, y en el punto más bajo, junio de 2023, la caída desde la entrada habría sido de 30,87 %.
No se debe sobrevalorar este resultado. Es un solo periodo de doce meses que empieza cerca de un máximo de precios, y la regla habría recomendado liquidez en casi cualquier fecha. Que haya salido bien no demuestra que el modelo anticipara la caída. Lo que se puede defender es el procedimiento: una regla fijada antes de estimar y aplicada sin excepciones.
Supuestos. La liquidez se evaluó con rendimiento de 0 %, un supuesto conservador. Los precios de entrada y de salida son promedios mensuales, no precios negociables en un día. Y no se modela el instrumento: en la práctica la exposición se tomaría con futuros o con un fondo cotizado, que tienen costos no considerados aquí.
8.3 Simulación de escenarios
BOOTSTRAP <- jb$p.value < 0.05 # si Jarque-Bera rechaza la normalidad, usamos bootstrap (lo pide el enunciado)
cat("Jarque-Bera p =", signif(jb$p.value, 4), "->",
ifelse(BOOTSTRAP, "bootstrap de residuos", "innovaciones normales"), "\n")## Jarque-Bera p = 0 -> bootstrap de residuosset.seed(SEMILLA) # semilla fija: la simulación sale igual cada vez que se corre
# 1.000 caminos posibles del precio a 12 meses. Bootstrap = los choques se sacan al azar de los errores
# reales del modelo (en vez de inventarlos normales), así se conservan las colas pesadas
sims <- replicate(N_SIM, as.numeric(
simulate(fit, nsim = H_DEC, future = TRUE, bootstrap = BOOTSTRAP,
xreg = if (usa_X) X_esc[1:H_DEC, x_vars, drop = FALSE] else NULL)))
P_final <- sims[H_DEC, ] # precio del mes 12 en cada uno de los 1.000 caminos
valor_fin <- CAPITAL * P_final / P0 # cuánta plata tendríamos si hubiéramos comprado
# Resumen de la distribución: lo que pide el enunciado (media o mediana, percentiles y probabilidad de pérdida)
resumen_sim <- data.frame(
medida = c("Media", "Mediana (P50)", "Percentil 5", "Percentil 95", "Probabilidad de pérdida (%)"),
precio = c(mean(P_final), median(P_final), quantile(P_final, 0.05), quantile(P_final, 0.95), NA),
portafolio = c(mean(valor_fin), median(valor_fin), quantile(valor_fin, 0.05), quantile(valor_fin, 0.95),
mean(valor_fin < CAPITAL) * 100)) # % de caminos que terminan con menos de USD 100.000
knitr::kable(resumen_sim, digits = 2, caption = paste("Simulación de", N_SIM, "trayectorias a", H_DEC, "meses"))| medida | precio | portafolio |
|---|---|---|
| Media | 110.65 | 108882.16 |
| Mediana (P50) | 99.27 | 97692.34 |
| Percentil 5 | 47.28 | 46526.42 |
| Percentil 95 | 203.28 | 200037.16 |
| Probabilidad de pérdida (%) | NA | 52.40 |
banda_inf <- apply(sims, 1, quantile, 0.05) # percentil 5 de cada mes (escenario malo)
banda_sup <- apply(sims, 1, quantile, 0.95) # percentil 95 de cada mes (escenario bueno)
cat("Meses reales dentro de la banda 5 %-95 % simulada:",
round(mean(precios_dec >= banda_inf & precios_dec <= banda_sup) * 100, 1), "%\n") # ¿los precios reales cayeron en la banda?## Meses reales dentro de la banda 5 %-95 % simulada: 100 %# Gráfico de abanico: 200 de los 1.000 caminos, sus percentiles y el precio que realmente pasó
plot(fechas_dec, precios_dec, type = "n", xlab = "", ylab = "USD/barril",
ylim = range(sims[, 1:200], precios_dec),
main = paste(N_SIM, "trayectorias simuladas a", H_DEC, "meses (se muestran 200)")) # lienzo vacío
matlines(fechas_dec, sims[, 1:200], lty = 1, col = adjustcolor("steelblue", 0.15)) # caminos en azul claro
matlines(fechas_dec, t(apply(sims, 1, quantile, c(0.05, 0.5, 0.95))), lty = c(2, 1, 2), col = "black", lwd = 2) # percentiles 5/50/95
lines(fechas_dec, precios_dec, col = "red", lwd = 2) # precio real en rojo
abline(h = P0, lty = 3) # precio de entrada
legend("topleft", c("Percentiles 5/50/95", "Observado", "Precio de entrada"),
col = c("black", "red", "black"), lty = c(1, 1, 3), lwd = 2, bty = "n")
# Histograma del valor final: muestra todo el rango de posibles resultados de la inversión
hist(valor_fin, breaks = 50, col = "grey80", border = "white", xlab = "USD",
main = paste("Valor del portafolio a", H_DEC, "meses si se compra en T0"))
abline(v = CAPITAL, col = "red", lwd = 2, lty = 2) # rojo: los USD 100.000 iniciales
abline(v = quantile(valor_fin, c(0.05, 0.5, 0.95)), col = "steelblue", lwd = 2) # azul: percentiles 5, 50 y 95
Para medir el riesgo se simularon 1.000 trayectorias del precio a doce meses a partir del modelo final, con semilla fija. Como Jarque-Bera rechaza la normalidad de los residuos, los choques no se generan de una normal sino que se toman al azar de los residuos observados (bootstrap). Así se conservan las colas pesadas; con choques normales la simulación subestimaría la probabilidad de los escenarios extremos (Hyndman y Athanasopoulos, 2021).
Se simula la posición larga aunque la decisión fue no comprar. Con la decisión tomada el portafolio vale USD 100.000 en todos los escenarios; la simulación muestra la distribución de resultados de la alternativa descartada.
El precio simulado a doce meses tiene una mediana de 99,27 USD y una media de 110,65 USD. La diferencia refleja la asimetría hacia la derecha. El percentil 5 es 47,28 USD y el percentil 95 es 203,28 USD. En valor del portafolio: USD 46.526 en el escenario adverso, USD 97.692 en la mediana y USD 200.037 en el favorable. La probabilidad de terminar con menos del capital inicial es 52,4 %. La mediana simulada difiere un poco del pronóstico analítico (97,61) porque proviene de mil sorteos de residuos que no son simétricos.
La compra era, según el modelo, una apuesta con algo más de la mitad de probabilidad de pérdida, un escenario típico ligeramente negativo y un rango de resultados que va desde perder más de la mitad del capital hasta duplicarlo. El precio observado quedó dentro de la banda del 5 % al 95 % en todos los meses, pero eso se debe a lo ancha que es la banda.
Hay una limitación que afecta estos números. El bootstrap sortea los residuos de forma independiente y rompe la agrupación de volatilidad detectada por la prueba ARCH y modelada en la sección 6.1; la simulación no usa esa varianza condicional. En la realidad los meses extremos vienen juntos; en la simulación quedan dispersos. Estos escenarios son una herramienta para dimensionar el riesgo y no precios futuros garantizados.
9. Conclusiones y limitaciones
9.1 ¿Qué características de la serie fueron determinantes?
Cuatro rasgos definieron el modelo. El logaritmo del precio tiene raíz unitaria, lo que obliga a trabajar con el retorno (\(d = 1\)). El retorno tiene memoria muy corta, con un solo rezago significativo (0,294), lo que lleva a un modelo de orden uno. No hay estacionalidad, así que \(D = 0\) y no se incluyen términos estacionales. Y la varianza cambia en el tiempo, lo que no cambia la especificación de la media pero sí la precisión de su coeficiente y la confianza en sus intervalos; por eso se modeló con un ARCH(1). El precio se relaciona además con el dólar y la producción industrial, pero en el mismo mes, y por eso esa relación no se puede aprovechar para decidir.
9.2 ¿Cuál modelo fue seleccionado y por qué?
Se seleccionó un ARIMA(0,1,1) sobre el logaritmo del precio, sin deriva, con \(\theta = 0{,}3194\). Fue el de menor AICc entre los modelos que cumplieron el diagnóstico y que eran utilizables con la información de julio de 2022. El ARIMAX ajustaba mejor, pero requería conocer valores futuros de las exógenas; los modelos ARIMA(2,1,1) y (2,1,2) tenían mejor AICc, pero no eran invertibles. En la prueba el modelo superó por poco a la caminata aleatoria, sin ser el mejor de la tabla. Para la varianza se seleccionó un ARCH(1), con la misma lógica; al estimarlo, el coeficiente de media móvil bajó a 0,116 y dejó de ser significativo, pero la decisión de inversión es la misma con las dos estimaciones.
9.3 ¿Los residuos cumplen razonablemente las condiciones esperadas?
De forma parcial. No tienen autocorrelación y su media es cero, por lo que el pronóstico puntual es utilizable. No son normales y presentan efecto ARCH, por lo que los intervalos de predicción no son confiables y la simulación tuvo que hacerse con bootstrap. El efecto ARCH se modeló en la sección 6.1 con un ARCH(1), y sus residuos estandarizados ya no presentan ese efecto (p = 0,964).
9.4 Limitaciones y extensiones
- El modelo solo aporta información para el primer mes. Su pronóstico a doce meses es el de un mes, y la regla de decisión resulta poco sensible al modelo.
- La significancia del coeficiente de media móvil depende del supuesto de varianza constante: al modelar la varianza baja de 0,319 a 0,116 y deja de ser significativo. La memoria del retorno es más débil de lo que sugiere el ARIMA(0,1,1), y el modelo queda todavía más cerca de una caminata aleatoria.
- La decisión depende de usar la mediana y no la media del pronóstico.
- La serie es un promedio mensual de precios diarios. Working (1960) mostró que promediar una caminata aleatoria genera por sí solo una autocorrelación de primer orden cercana a 0,25 en los cambios, un valor parecido al 0,294 encontrado aquí. Parte de la estructura que captura el MA(1) puede venir de la forma en que está construido el dato. Además, los precios de entrada y salida no son negociables.
- Se detectó efecto ARCH y se estimó un modelo para la varianza, pero solo como diagnóstico: los intervalos de predicción y la simulación no se recalcularon con la varianza condicional, y la simulación trata los choques como independientes.
- Los modelos de varianza se estimaron con innovaciones normales y con 198 retornos mensuales, una muestra pequeña para este tipo de modelos; sus coeficientes deben leerse con cautela.
- Los intervalos de predicción suponen normalidad y varianza constante, y resultaron demasiado anchos para ser útiles.
- Box-Cox no respaldó el logaritmo; se usó por criterio financiero.
- Hay un solo origen de pronóstico y un solo periodo de evaluación financiera, que además empieza cerca de un máximo de precios.
- Las exógenas se proyectaron con un supuesto simple y la relación estimada con la tasa de la Fed no se mantuvo fuera de muestra.
- La liquidez se evaluó al 0 % y no se modelaron el instrumento ni los costos de transacción.
Como extensiones quedan usar la varianza condicional estimada en los intervalos y en la simulación, el uso de exógenas rezagadas que sí estén disponibles en la fecha de decisión, y repetir el ejercicio con precios de fin de mes.
La respuesta a la pregunta del trabajo es que, con la regla definida y el modelo seleccionado, en julio de 2022 no convenía comprar, y el resultado posterior fue coherente con esa decisión. La evidencia no permite afirmar que el modelo pueda anticipar el precio del petróleo a un año.
10. Referencias y reproducibilidad
Baumeister, C. y Kilian, L. (2016). Understanding the decline in the price of oil since June 2014. Journal of the Association of Environmental and Resource Economists, 3(1), 131-158.
Board of Governors of the Federal Reserve System. (2026a). Federal funds effective rate [FEDFUNDS]. FRED, Federal Reserve Bank of St. Louis. Recuperado en noviembre de 2026 de https://fred.stlouisfed.org/series/FEDFUNDS
Board of Governors of the Federal Reserve System. (2026b). Industrial production: Total index [INDPRO]. FRED, Federal Reserve Bank of St. Louis. Recuperado en noviembre de 2026 de https://fred.stlouisfed.org/series/INDPRO
Board of Governors of the Federal Reserve System. (2026c). Nominal broad U.S. dollar index [DTWEXBGS]. FRED, Federal Reserve Bank of St. Louis. Recuperado en noviembre de 2026 de https://fred.stlouisfed.org/series/DTWEXBGS
Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307-327.
Box, G. E. P., Jenkins, G. M., Reinsel, G. C. y Ljung, G. M. (2015). Time series analysis: Forecasting and control (5.ª ed.). Wiley.
Brooks, C. (2008). Introductory econometrics for finance (2.ª ed.). Cambridge University Press.
CME Group. (s. f.). Crude oil futures: Contract specs. Recuperado en noviembre de 2026 de https://www.cmegroup.com/trading/energy/crude-oil/light-sweet-crude_contractSpecs_futures.html
Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007.
Fattouh, B. (2011). An anatomy of the crude oil pricing system (OIES Paper WPM 40). Oxford Institute for Energy Studies.
Galanos, A. (2026). rugarch: Univariate GARCH models [paquete de R]. https://doi.org/10.32614/CRAN.package.rugarch
Gujarati, D. N. y Porter, D. C. (2010). Econometría (5.ª ed.). McGraw-Hill.
Hamilton, J. D. (2009). Causes and consequences of the oil shock of 2007-08. Brookings Papers on Economic Activity, primavera, 215-261.
Hyndman, R. J. y Athanasopoulos, G. (2021). Forecasting: Principles and practice (3.ª ed.). OTexts. https://otexts.com/fpp3/
Hyndman, R. J. y Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679-688.
Hyndman, R., Athanasopoulos, G., Bergmeir, C., Caceres, G., Chhay, L., O’Hara-Wild, M., Petropoulos, F., Razbash, S., Wang, E. y Yasmeen, F. (2026). forecast: Forecasting functions for time series and linear models (versión 9.0.2) [paquete de R]. https://doi.org/10.32614/CRAN.package.forecast
Rodríguez, D. y Sierra, S. (2025). actfts: Autocorrelation tools featured for time series (versión 0.3.0) [paquete de R]. https://doi.org/10.32614/CRAN.package.actfts
Rodríguez Guevara, D. E. (2026a). Modelos ARCH-GARCH [diapositivas de clase]. Instituto Tecnológico Metropolitano.
Rodríguez Guevara, D. E. (2026b). Modelos ARIMA [diapositivas de clase]. Instituto Tecnológico Metropolitano.
Rodríguez Guevara, D. E. y González Uribe, G. J. (2017). Principios de econometría. Fondo Editorial ITM.
Tsay, R. S. (2010). Analysis of financial time series (3.ª ed.). Wiley.
U.S. Energy Information Administration. (2017, 15 de noviembre). Transportation constraints and export costs widen the Brent-WTI crude oil price spread. Today in Energy. https://www.eia.gov/todayinenergy/detail.php?id=33752
U.S. Energy Information Administration. (2020, 27 de abril). Low liquidity and limited available storage pushed WTI crude oil futures prices below zero. Today in Energy. https://www.eia.gov/todayinenergy/detail.php?id=43495
U.S. Energy Information Administration. (2023, 4 de enero). Crude oil prices increased in first-half 2022 and declined in second-half 2022. Today in Energy. https://www.eia.gov/todayinenergy/detail.php?id=55079
U.S. Energy Information Administration. (2026a, 7 de abril). Crude oil and petroleum product prices increased sharply in the first quarter of 2026. Today in Energy. https://www.eia.gov/todayinenergy/detail.php?id=67424
U.S. Energy Information Administration. (2026b). Crude oil prices: West Texas Intermediate (WTI) – Cushing, Oklahoma [MCOILWTICO]. FRED, Federal Reserve Bank of St. Louis. Recuperado en noviembre de 2026 de https://fred.stlouisfed.org/series/MCOILWTICO
Working, H. (1960). Note on the correlation of first differences of averages in a random chain. Econometrica, 28(4), 916-918.
cat("Datos: datos_procesados.csv | Semilla de la simulación:", SEMILLA, "\n") # con estos dos, cualquiera obtiene lo mismo## Datos: datos_procesados.csv | Semilla de la simulación: 2026sessionInfo() # versión de R y de cada paquete (lo pide el enunciado para la reproducibilidad)## R version 4.6.2 (2026-10-29)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.3 lmtest_0.9-40 FinTS_0.4-9 tseries_0.10-63 forecast_9.0.2 quantmod_0.4.29 TTR_0.24.4 xts_0.14.3 zoo_1.9-1
## [10] actfts_0.3.0 rugarch_1.5-6
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.3 timeDate_4052.112 dplyr_1.2.1 farver_2.1.2
## [6] S7_0.2.2 fastmap_1.2.0 pracma_2.4.6 GeneralizedHyperbolic_0.8-7 digest_0.6.39
## [11] lifecycle_1.0.5 Rsolnp_2.0.1 magrittr_2.0.5 compiler_4.6.2 rlang_1.3.0
## [16] sass_0.4.10 tools_4.6.2 yaml_2.3.12 data.table_1.18.6.1 knitr_1.52
## [21] labeling_0.4.3 htmlwidgets_1.6.4 mclust_6.1.3 curl_8.1.0 RColorBrewer_1.1-3
## [26] KernSmooth_2.23-27 withr_3.0.3 purrr_1.2.2 numDeriv_2016.8-1.1 grid_4.6.2
## [31] colorspace_2.1-4 future_1.82.0 globals_0.20.0 scales_1.4.0 MASS_7.3-66
## [36] cli_3.6.6 mvtnorm_1.4-2 rmarkdown_2.32 generics_0.1.4 otel_0.3.0
## [41] rstudioapi_0.19.0 future.apply_1.21.0 httr_1.4.9 cachem_1.1.0 urca_1.3-4
## [46] vctrs_0.7.3 spd_2.0-1 Matrix_1.7-6 jsonlite_2.0.0 listenv_1.1.0
## [51] crosstalk_1.2.2 plotly_4.12.1 jquerylib_0.1.4 tidyr_1.3.2 glue_1.8.1
## [56] parallelly_1.50.0 nloptr_2.2.1 SkewHyperbolic_0.4-2 codetools_0.2-20 DistributionUtils_0.6-2
## [61] stringi_1.8.9 gtable_0.3.6 quadprog_1.5-8 tibble_3.3.1 pillar_1.11.1
## [66] htmltools_0.5.9 reactable_0.4.5 truncnorm_1.0-9 R6_2.6.1 ks_1.15.3
## [71] evaluate_1.0.5 lattice_0.23-1 openxlsx_4.2.9 fracdiff_1.5-4 bslib_0.12.0
## [76] Rcpp_1.1.2 zip_3.0.2 nlme_3.1-171 xfun_0.62 pkgconfig_2.0.3