Análisis y pronóstico de la concentración de Ozono troposférico (O₃) en Cali mediante modelos de series de tiempo: metodología Box-Jenkins y ARIMA
Taller de Series de Tiempo — Gestión de Datos
Universidad del Valle (Cali) · Profesor: Orlando Joaqui Barandica
16 junio 2026
- 1 Introducción
- 2 Metodología
- 3 Descripción de la serie temporal
- 4 Resultados del modelo ARIMA
- 5 Pronóstico
- 6 Pruebas sobre los residuales
- 7 Gráfico final
- 8 Conclusiones
- 9 Bibliografía
1 Introducción
La calidad del aire constituye uno de los determinantes ambientales más relevantes para la salud pública en las grandes ciudades de América Latina, y Santiago de Cali no es la excepción. Como tercera ciudad más poblada de Colombia, ubicada en un valle interandino con condiciones meteorológicas particulares —alta radiación solar, baja velocidad del viento y temperaturas elevadas durante buena parte del año—, Cali presenta un escenario propicio para la formación de contaminantes secundarios como el ozono troposférico (O₃) (DAGMA, 2019; OMS, 2021).
A diferencia del ozono estratosférico, que protege la vida en la Tierra, el ozono a nivel del suelo es un contaminante perjudicial que se forma a partir de reacciones fotoquímicas entre óxidos de nitrógeno (NOₓ) y compuestos orgánicos volátiles (COV) en presencia de luz solar (EPA, 2023; OMS, 2021). Por esta razón, el O₃ no se emite directamente, sino que se genera en la atmósfera, alcanzando sus máximas concentraciones en las horas centrales del día, cuando la radiación solar es más intensa. La exposición prolongada a concentraciones elevadas de ozono se asocia con afecciones respiratorias, agravamiento del asma y reducción de la función pulmonar, lo que convierte su monitoreo y pronóstico en una herramienta fundamental para la gestión ambiental urbana (OMS, 2021; EPA, 2023).
El análisis de series de tiempo emerge como una metodología analítica poderosa para transformar el registro histórico de mediciones ambientales en horizontes de predicción, permitiendo identificar patrones diurnos, ciclos estacionales y tendencias que orientan la toma de decisiones en materia de política pública ambiental (Box & Jenkins, 1970; Hyndman & Athanasopoulos, 2018). En particular, cuando se aplica a contaminantes atmosféricos, esta metodología proporciona insumos valiosos para anticipar episodios críticos de contaminación, emitir alertas tempranas a la población vulnerable y evaluar la efectividad de las medidas de mitigación (DAGMA, 2019).
Este trabajo analiza la serie temporal de la concentración de ozono troposférico (O₃, en µg/m³) registrada por la red de monitoreo de calidad del aire de Cali durante el año 2018, con resolución horaria. A partir de la agregación de los datos a promedios diarios, se construye y valida un modelo ARIMA siguiendo la metodología de Box-Jenkins, con el objetivo de caracterizar la dinámica temporal de la serie, evaluar la calidad del ajuste mediante pruebas diagnósticas sobre los residuales y generar pronósticos de corto plazo con sus respectivas bandas de confianza. De esta manera, el trabajo busca cerrar la brecha entre la teoría del análisis de series temporales y su aplicación práctica a un problema ambiental concreto y de alta relevancia para la ciudad.
2 Metodología
2.1 Serie de tiempo
Una serie de tiempo es una secuencia de observaciones de una variable específica registradas en intervalos de tiempo regulares y ordenadas cronológicamente (Box & Jenkins, 1970). Estas observaciones pueden estar espaciadas en forma horaria, diaria, mensual o anual, dependiendo del fenómeno bajo estudio. En el contexto ambiental, las series de tiempo permiten analizar el comportamiento histórico de variables como la concentración de contaminantes, la temperatura, la humedad o la radiación solar, con el objetivo de identificar patrones subyacentes y realizar pronósticos sobre su evolución futura (Hyndman & Athanasopoulos, 2018).
El análisis de series temporales se fundamenta en el estudio de cuatro componentes que caracterizan el comportamiento de los datos a través del tiempo. La tendencia representa el patrón de largo plazo que muestra la dirección general de la serie, ya sea creciente, decreciente o constante. La estacionalidad corresponde a fluctuaciones periódicas que se repiten en intervalos fijos de tiempo —en el caso del ozono, un marcado ciclo diurno y posibles variaciones semanales asociadas a la actividad humana—. El componente cíclico refleja oscilaciones de período no fijo, vinculadas en este caso al régimen climático bimodal del valle del Cauca (dos temporadas secas y dos lluviosas al año). Finalmente, el componente aleatorio o residual representa la variabilidad que permanece después de remover los componentes anteriores.
2.2 Modelo ARIMA
El modelo ARIMA (AutoRegressive Integrated Moving Average) es una de las herramientas estadísticas más utilizadas y versátiles para el análisis y pronóstico de series temporales (Hyndman & Athanasopoulos, 2018). Combina tres componentes fundamentales que trabajan de manera conjunta para capturar la estructura de dependencia temporal de los datos. El componente Autorregresivo (AR) establece que el valor actual de la serie puede explicarse como una función lineal de sus valores pasados. El componente Integrado (I) se refiere al proceso de diferenciación aplicado a la serie para lograr estacionariedad, removiendo tendencias. El componente Media Móvil (MA) modela la dependencia entre una observación y los errores de pronóstico de observaciones pasadas.
La notación del modelo se expresa como ARIMA(p,d,q). El componente AR(p) se describe en la ecuación 1:
\[ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \ldots + \phi_p y_{t-p} + \varepsilon_t \tag{1} \]
donde \(\varepsilon_t\) es ruido blanco. El parámetro \(d\) indica el número de diferenciaciones necesarias para alcanzar la estacionariedad. La primera diferenciación (\(d=1\)) se muestra en la ecuación 2:
\[ y'_t = y_t - y_{t-1} \tag{2} \]
El componente MA(q) se presenta en la ecuación 3:
\[ y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \ldots + \theta_q \varepsilon_{t-q} \tag{3} \]
y el modelo ARMA(\(p, q\)) general se expresa en la ecuación 4:
\[ y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q} \tag{4} \]
Cuando la serie presenta estacionalidad de período \(s\), el modelo se extiende a su versión estacional SARIMA(p,d,q)(P,D,Q)\(_s\), que incorpora términos autorregresivos y de media móvil sobre los rezagos estacionales.
2.3 Metodología de Box-Jenkins
La estimación de modelos ARIMA se realiza tradicionalmente mediante la metodología de Box-Jenkins, un procedimiento sistemático e iterativo que consta de cuatro etapas fundamentales (Box & Jenkins, 1970).
2.3.1 Primera Etapa: Identificación del Modelo
En esta fase se debe asegurar que la serie sea estacionaria, es decir, que su media, varianza y autocorrelación se mantengan constantes a través del tiempo. Para verificar la estacionariedad se aplica la prueba de Dickey-Fuller Aumentada (ADF), cuya hipótesis nula establece que la serie posee una raíz unitaria y por tanto no es estacionaria:
- \(H_0\): La serie tiene raíz unitaria (no es estacionaria).
- \(H_1\): La serie no tiene raíz unitaria (es estacionaria).
Si la prueba ADF arroja un p-valor inferior al nivel de significancia (típicamente 0.05), se rechaza \(H_0\) y se concluye que la serie es estacionaria. Cuando la serie no es estacionaria, se aplica diferenciación. Los órdenes \(p\) y \(q\) se identifican mediante el análisis de las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF).
2.3.2 Segunda Etapa: Estimación de Parámetros
Una vez identificados los órdenes candidatos, se estiman los coeficientes del modelo mediante máxima verosimilitud. Se comparan varias especificaciones utilizando criterios de información como el AIC, el AICc y el BIC, prefiriéndose el modelo con menor valor, en equilibrio con la parsimonia.
2.3.3 Tercera Etapa: Diagnóstico y Validación del Modelo
Se verifica que los residuales se comporten como ruido blanco: media cero, varianza constante y ausencia de autocorrelación. Para ello se aplican la prueba de Ljung-Box (autocorrelación de los residuales), la prueba de Shapiro-Wilk (normalidad) y el análisis gráfico de la ACF/PACF de los residuales.
2.3.4 Cuarta Etapa: Predicción o Pronóstico
Validado el modelo, se generan pronósticos para horizontes futuros junto con sus intervalos de confianza, y se evalúa su desempeño predictivo sobre una ventana de prueba mediante métricas como MAE, RMSE y MAPE.
2.4 Descripción de las variables y fuente de los datos
La base de datos contiene 8.760 registros horarios correspondientes al año 2018 para una estación de monitoreo de la ciudad de Cali. El cuadro siguiente describe las variables disponibles; la variable de interés para este taller es la concentración de ozono troposférico (O₃).
# Tabla de descripción de variables
tabla_var <- data.frame(
Variable = c("Fecha & Hora", "O3", "Vel Viento", "Dir Viento",
"Temperatura", "Humedad", "Radiacion Solar", "Lluvia"),
Descripcion = c("Marca temporal de la medición (resolución horaria)",
"**Concentración de ozono troposférico (variable objetivo)**",
"Velocidad del viento",
"Dirección del viento",
"Temperatura del aire",
"Humedad relativa",
"Radiación solar global",
"Precipitación acumulada"),
Unidad = c("fecha-hora", "µg/m³", "m/s", "grados", "°C", "%", "Watt/m²", "mm"),
stringsAsFactors = FALSE
)
tabla_var %>%
gt() %>%
tab_header(
title = md("**Descripción de las variables disponibles**"),
subtitle = "Base de datos horaria de calidad del aire — Cali 2018"
) %>%
cols_label(
Variable = "Variable",
Descripcion = "Descripción",
Unidad = "Unidad"
) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_style(
style = list(cell_fill(color = verde$claro), cell_text(weight = "bold")),
locations = cells_body(columns = Variable)
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(columns = Variable, rows = Variable == "O3")
) %>%
fmt_markdown(columns = Descripcion) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(16),
heading.subtitle.font.size = px(13))| Descripción de las variables disponibles | ||
| Base de datos horaria de calidad del aire — Cali 2018 | ||
| Variable | Descripción | Unidad |
|---|---|---|
| Fecha & Hora | Marca temporal de la medición (resolución horaria) | fecha-hora |
| O3 | Concentración de ozono troposférico (variable objetivo) | µg/m³ |
| Vel Viento | Velocidad del viento | m/s |
| Dir Viento | Dirección del viento | grados |
| Temperatura | Temperatura del aire | °C |
| Humedad | Humedad relativa | % |
| Radiacion Solar | Radiación solar global | Watt/m² |
| Lluvia | Precipitación acumulada | mm |
El ozono troposférico (O₃) se selecciona como variable de interés por dos razones principales. En primer lugar, por su relevancia en salud pública: es un contaminante secundario que se forma fotoquímicamente y cuya exposición prolongada se asocia con afecciones respiratorias, por lo que la OMS establece un valor guía de 100 µg/m³ como media de 8 horas (OMS, 2021). En segundo lugar, por su dinámica temporal compleja y adecuación metodológica: al depender de la radiación solar y las condiciones meteorológicas, presenta marcados ciclos diurnos y estacionales (EPA, 2023), lo que lo convierte en un candidato ideal para el modelado con ARIMA y la metodología Box‑Jenkins. Las demás variables (temperatura, radiación, humedad, viento y lluvia) son factores explicativos, pero el ozono es el fenómeno de interés final por su impacto directo en la población y su estructura temporal rica en autocorrelación, especialmente relevante en el contexto del valle de Cali, donde las condiciones climáticas favorecen su formación (DAGMA, 2019).
3 Descripción de la serie temporal
3.1 Carga y preparación de los datos
A continuación se carga la base, se renombran las columnas a nombres manejables, se construye la marca temporal y se agrega la serie horaria a promedios diarios. La agregación diaria permite obtener una serie limpia de 365 observaciones para 2018, ideal para el modelado ARIMA, mientras que el ciclo diurno (intra-día) se analiza por separado en la sección descriptiva.
# --- Lectura del archivo Excel ---
ruta <- "Compartir2.xlsx"
crudo <- read_excel(ruta, sheet = 1)
# Renombramos las columnas
names(crudo) <- c("FechaHora", "O3", "VelViento", "DirViento",
"Temperatura", "Humedad", "Radiacion", "Lluvia")
# Aseguramos el tipo fecha-hora y extraemos la fecha (día) y la hora
crudo <- crudo %>%
mutate(
FechaHora = as.POSIXct(FechaHora),
Fecha = as.Date(FechaHora),
Hora = hour(FechaHora),
Mes = month(FechaHora, label = TRUE, abbr = TRUE),
DiaSemana = wday(FechaHora, label = TRUE, abbr = TRUE, week_start = 1)
) %>%
filter(Fecha >= as.Date("2018-01-01") & Fecha <= as.Date("2018-12-31"))
# --- Agregación a promedios DIARIOS---
diario <- crudo %>%
group_by(Fecha) %>%
summarise(
O3 = mean(O3, na.rm = TRUE),
Temperatura = mean(Temperatura, na.rm = TRUE),
Humedad = mean(Humedad, na.rm = TRUE),
Radiacion = mean(Radiacion, na.rm = TRUE),
VelViento = mean(VelViento, na.rm = TRUE),
Lluvia = sum(Lluvia, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(Fecha)
diario <- diario %>% mutate(across(where(is.numeric),
~ ifelse(is.nan(.x), NA, .x)))
# Construimos el objeto ts con periodicidad semanal (frecuencia 7) e
# interpolamos los pocos huecos para garantizar una serie continua.
serie_ts <- ts(diario$O3, frequency = 7)
serie_ts <- na.interp(serie_ts) # interpolación de valores faltantes
diario$O3 <- as.numeric(serie_ts) # serie diaria de O3 ya completa
datatable(
diario,
rownames = FALSE,
colnames = c("Fecha", "O₃ (µg/m³)", "Temperatura (°C)", "Humedad (%)",
"Radiación (W/m²)", "Vel. Viento (m/s)", "Lluvia (mm)"),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left;',
htmltools::strong('Serie diaria de variables ambientales — Cali 2018'),
htmltools::br(),
sprintf('Observaciones: %d días | Rango: %s a %s | O₃ medio: %.2f µg/m³',
nrow(diario),
format(min(diario$Fecha), "%d %b %Y"),
format(max(diario$Fecha), "%d %b %Y"),
mean(diario$O3))
),
options = list(pageLength = 7, dom = 'tip')
) %>%
formatRound(columns = 2:7, digits = 2)3.2 Estadísticas descriptivas
resumen_var <- function(x, etiqueta, unidad) {
data.frame(
Variable = etiqueta,
Unidad = unidad,
Media = round(mean(x, na.rm = TRUE), 2),
Mediana = round(median(x, na.rm = TRUE), 2),
Desv = round(sd(x, na.rm = TRUE), 2),
Minimo = round(min(x, na.rm = TRUE), 2),
Maximo = round(max(x, na.rm = TRUE), 2),
CV = round(sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100, 1),
stringsAsFactors = FALSE
)
}
tabla_desc <- bind_rows(
resumen_var(diario$O3, "Ozono (O3)", "ug/m3"),
resumen_var(diario$Temperatura, "Temperatura", "C"),
resumen_var(diario$Humedad, "Humedad relativa", "%"),
resumen_var(diario$Radiacion, "Radiacion solar", "W/m2"),
resumen_var(diario$VelViento, "Velocidad viento", "m/s"),
resumen_var(diario$Lluvia, "Lluvia (diaria)", "mm")
)
verde <- list(profundo = "#0B6E4F", medio = "#1B998B", claro = "#E8F5F1",
texto = "#FFFFFF")
tabla_desc %>%
gt() %>%
tab_header(
title = md("**Tabla 1. Estadísticas descriptivas de las variables (promedios diarios, 2018)**"),
subtitle = md("Estación de monitoreo de calidad del aire — Cali")
) %>%
cols_label(
Desv = "Desv. Est.", CV = "CV (%)"
) %>%
fmt_number(columns = c(Media, Mediana, Desv, Minimo, Maximo, CV), decimals = 2) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_style(
style = list(cell_fill(color = verde$claro), cell_text(weight = "bold")),
locations = cells_body(columns = Variable)
) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(16),
heading.subtitle.font.size = px(13))| Tabla 1. Estadísticas descriptivas de las variables (promedios diarios, 2018) | |||||||
| Estación de monitoreo de calidad del aire — Cali | |||||||
| Variable | Unidad | Media | Mediana | Desv. Est. | Minimo | Maximo | CV (%) |
|---|---|---|---|---|---|---|---|
| Ozono (O3) | ug/m3 | 29.59 | 29.40 | 7.55 | 12.70 | 58.68 | 25.50 |
| Temperatura | C | 28.28 | 28.43 | 1.21 | 24.34 | 30.93 | 4.30 |
| Humedad relativa | % | 57.66 | 57.06 | 5.53 | 42.50 | 73.20 | 9.60 |
| Radiacion solar | W/m2 | 108.24 | 104.14 | 38.00 | 11.17 | 253.80 | 35.10 |
| Velocidad viento | m/s | 2.51 | 2.48 | 0.31 | 1.80 | 3.77 | 12.50 |
| Lluvia (diaria) | mm | 6.59 | 0.00 | 16.37 | 0.00 | 108.96 | 248.50 |
La concentración media diaria de ozono en Cali durante 2018 fue de 29.6 µg/m³, con una mediana muy similar (29.4 µg/m³) y un valor máximo de 58.7 µg/m³, que se mantiene por debajo del umbral de referencia de la Organización Mundial de la Salud (OMS, 2021) de 100 µg/m³ como media de 8 horas. El coeficiente de variación del 25.5 % indica una variabilidad interdiaria moderada, reflejando la sensibilidad del O₃ a las condiciones meteorológicas cambiantes (EPA, 2023). La temperatura media (28.3 °C) y la radiación solar (108 W/m² en promedio, pero con alta dispersión) confirman un ambiente cálido y soleado, propicio para la fotoquímica, mientras que la humedad relativa (57.7 %) es consistente con el régimen bimodal seco‑húmedo del valle del Cauca (DAGMA, 2019).
La velocidad del viento es baja (2.5 m/s en promedio), lo que limita la dispersión de contaminantes y favorece la acumulación de precursores del ozono (EPA, 2023). La lluvia diaria presenta una mediana de 0 mm y un coeficiente de variación extremadamente alto (248.5 %), evidenciando que la precipitación se concentra en pocos días al año (DAGMA, 2019). Esta distribución sesgada implica que largos periodos secos —con alta radiación y escasa remoción húmeda— son los más críticos para episodios elevados de ozono (OMS, 2021). En conjunto, las variables meteorológicas descritas justifican la necesidad de modelos predictivos como el ARIMA presentado en este trabajo.
3.3 Evolución temporal de la serie
# Media móvil de 30 días para visualizar la tendencia de fondo
diario <- diario %>%
mutate(MM30 = zoo::rollmean(O3, k = 30, fill = NA, align = "center"))
p_serie <- ggplot(diario, aes(x = Fecha)) +
geom_line(aes(y = O3), color = "#1B998B", linewidth = 0.5, alpha = 0.8) +
geom_line(aes(y = MM30), color = "#0B6E4F", linewidth = 1.1, na.rm = TRUE) +
labs(
title = "Promedio diario de ozono troposférico (O₃) en Cali — 2018",
subtitle = "Serie diaria (línea clara) y media móvil de 30 días (línea oscura)",
x = NULL, y = expression(O[3]~"(µg/m"^3*")")
) +
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))
ggplotly(p_serie) %>% layout(hovermode = "x unified")La concentración diaria de ozono troposférico en Cali promedió 29.6 µg/m³ en 2018, con la mitad central de los datos situada entre 24.5 y 34.8 µg/m³. La serie, sin embargo, registró un episodio sobresaliente a finales de agosto y comienzos de septiembre, cuando el ozono escaló hasta su valor máximo anual de 58.7 µg/m³ el 27 de agosto. Esa cima, muy por encima de cualquier otra del año, coincidió con el periodo más seco y de mayor insolación del calendario climático caleño y puso de manifiesto la capacidad del valle para generar concentraciones puntuales elevadas bajo condiciones atmosféricas favorables.
Al suavizar las oscilaciones diarias con una media móvil de 30 días, se aprecia con nitidez la respiración bimodal del ozono, gobernada por el régimen de precipitaciones del suroccidente colombiano. Los meses de mayor concentración promedio fueron sept y ago, mientras que los más lluviosos may y nov mostraron los niveles más bajos. En estos últimos, la nubosidad persistente y el lavado atmosférico provocaron descensos diarios que llegaron hasta los 12.7 µg/m³ registrados en mayo. La alternancia seca‑húmeda, típica del valle del Cauca, imprime así una cadencia estacional que permite anticipar, con razonable certeza, los momentos de mayor acumulación de ozono.
Más allá de las modulaciones estacionales, la línea clara de la serie diaria expone una marcada erraticidad de corto plazo: entre dos días consecutivos el ozono llegó a cambiar hasta 28.6 µg/m³, reflejando la alta sensibilidad del contaminante a factores meteorológicos puntuales (nubosidad, ráfagas de viento, picos de tráfico). Esta variabilidad impone un reto al monitoreo y al pronóstico, pues episodios peligrosos pueden gestarse en cuestión de horas. La ventana crítica se extiende de julio a septiembre, cuando la combinación de tiempo seco, fuerte radiación y posible mayor emisión de precursores eleva no solo el nivel medio sino también la probabilidad de alcanzar picos que comprometan la salud respiratoria de la población caleña.
3.4 Estacionalidad: ciclo diurno y patrón mensual
El ozono es un contaminante fuertemente gobernado por la luz solar; por ello su comportamiento intra-día es uno de los rasgos más característicos de la serie. El gráfico siguiente, construido a partir de los datos horarios originales, muestra el ciclo diurno promedio del O₃.
perfil_horario <- crudo %>%
group_by(Hora) %>%
summarise(O3_medio = mean(O3, na.rm = TRUE),
O3_sd = sd(O3, na.rm = TRUE), .groups = "drop")
p_diurno <- ggplot(perfil_horario, aes(x = Hora, y = O3_medio)) +
geom_ribbon(aes(ymin = O3_medio - O3_sd, ymax = O3_medio + O3_sd),
fill = "#1B998B", alpha = 0.18) +
geom_line(color = "#0B6E4F", linewidth = 1.2) +
geom_point(color = "#0B6E4F", size = 1.8) +
scale_x_continuous(breaks = seq(0, 23, 2)) +
labs(title = "Ciclo diurno promedio del Ozono (O3) — Cali 2018",
subtitle = "Promedio horario ± 1 desviación estándar",
x = "Hora del día", y = expression(O[3]~"(µg/m"^3*")")) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))
p_diurno
A escala intradiaria (gráfico del ciclo diurno) se observa el característico perfil fotoquímico: las concentraciones son mínimas hacia las 6:00 horas, cuando la ausencia de radiación solar detiene la producción de ozono, y ascienden progresivamente hasta alcanzar su punto máximo alrededor de las 15:00 horas, coincidiendo con la mayor insolación del día. La banda sombreada (±1 desviación estándar) evidencia una variabilidad horaria considerable, especialmente en las horas centrales, lo que sugiere que factores locales como la nubosidad o la intensidad del tráfico pueden modificar el patrón diario.
# Cálculo de estadísticos por mes
estad_mes <- crudo %>%
group_by(Mes) %>%
summarise(
n = n(),
min = min(O3, na.rm = TRUE),
Q1 = quantile(O3, 0.25, na.rm = TRUE),
mediana = median(O3, na.rm = TRUE),
media = mean(O3, na.rm = TRUE),
Q3 = quantile(O3, 0.75, na.rm = TRUE),
max = max(O3, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
texto_caja = sprintf(
"<b>%s</b><br>Mín: %.1f | Q1: %.1f<br>Mediana: %.1f | Q3: %.1f<br>Máx: %.1f<br>Media: %.1f µg/m³",
Mes, min, Q1, mediana, Q3, max, media
)
)
crudo_label <- crudo %>% left_join(estad_mes %>% select(Mes, texto_caja), by = "Mes")
p_mes_fix <- ggplot(crudo_label, aes(x = Mes, y = O3, text = texto_caja)) +
geom_boxplot(
fill = "#1B998B",
alpha = 0.55,
outlier.alpha = 0.15,
outlier.size = 0.6,
color = "#0B6E4F"
) +
stat_summary(
aes(text = paste("Media:", round(..y.., 1), "µg/m³")),
fun = mean,
geom = "point",
color = "#C0392B",
size = 2.5,
alpha = 0.9
) +
labs(
title = "Distribución mensual del Ozono (O₃) — Cali 2018",
subtitle = "Diagrama de cajas horario por mes · punto rojo: media mensual",
x = NULL,
y = "O₃ (µg/m³)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", color = "#0B6E4F", size = 14),
plot.subtitle = element_text(color = "#2C3E50", size = 11),
panel.grid.major.x = element_blank(),
legend.position = "none"
)
ggplotly(p_mes_fix, tooltip = "text", dynamicTicks = TRUE) %>%
layout(
hovermode = "closest",
hoverlabel = list(bgcolor = "white", font = list(color = "#0B6E4F", size = 12)),
margin = list(t = 80),
yaxis = list(title = "O₃ (µg/m³)", zeroline = FALSE)
) %>%
config(displayModeBar = TRUE, displaylogo = FALSE)A escala mensual , la serie exhibe la huella del régimen climático bimodal del valle del Cauca. Los meses de menor precipitación como sept y ago presentan las medianas y medias más altas, mientras que los meses más lluviosos —may y nov— muestran concentraciones considerablemente menores. Esta alternancia confirma que la temporada de lluvias, con su mayor nubosidad y lavado atmosférico, suprime la formación fotoquímica del ozono, mientras que los periodos secos y soleados favorecen su acumulación.
4 Resultados del modelo ARIMA
4.1 Ventanas de entrenamiento y prueba
Siguiendo la metodología de Box‑Jenkins, dividimos la serie diaria de ozono en dos conjuntos: entrenamiento y prueba. El conjunto de entrenamiento incluye los primeros 344 días del año (desde el 1 de enero hasta el 9 de diciembre de 2018) y se utiliza para identificar el modelo, estimar sus parámetros y seleccionar la mejor especificación. El conjunto de prueba contiene los últimos 21 días (del 10 al 31 de diciembre de 2018), equivalentes a tres semanas completas, y se reserva exclusivamente para evaluar la capacidad predictiva del modelo final. Esta práctica es estándar en series de tiempo: se entrena con el pasado y se valida con el futuro más reciente, simulando cómo se comportaría el modelo en una aplicación real.
h <- 21 # horizonte de pronóstico (días)
n <- nrow(diario)
train_df <- diario[1:(n - h), ]
test_df <- diario[(n - h + 1):n, ]
ventana <- ts(train_df$O3, frequency = 7) # entrenamiento
ventana2 <- ts(test_df$O3, frequency = 7,
start = tsp(ventana)[2] + 1/7)
tabla_particion <- data.frame(
Conjunto = c("Entrenamiento", "Prueba"),
Observaciones = c(length(ventana), length(ventana2)),
Fecha_inicio = c(format(min(train_df$Fecha), "%d %b %Y"),
format(min(test_df$Fecha), "%d %b %Y")),
Fecha_fin = c(format(max(train_df$Fecha), "%d %b %Y"),
format(max(test_df$Fecha), "%d %b %Y"))
)
tabla_particion %>%
gt() %>%
tab_header(
title = md("**Ventanas de entrenamiento y prueba**"),
subtitle = md("Serie diaria de O<sub>3</sub> — Cali 2018")
) %>%
cols_label(
Observaciones = "Obs.",
Fecha_inicio = "Inicio",
Fecha_fin = "Fin"
) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_style(
style = list(cell_fill(color = verde$claro), cell_text(weight = "bold")),
locations = cells_body(columns = Conjunto)
) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(16),
heading.subtitle.font.size = px(13))| Ventanas de entrenamiento y prueba | |||
| Serie diaria de O3 — Cali 2018 | |||
| Conjunto | Obs. | Inicio | Fin |
|---|---|---|---|
| Entrenamiento | 344 | 01 ene. 2018 | 10 dic. 2018 |
| Prueba | 21 | 11 dic. 2018 | 31 dic. 2018 |
p_ventanas <- ggplot() +
geom_line(data = train_df,
aes(Fecha, O3, color = "Entrenamiento"),
linewidth = 0.6) +
geom_line(data = test_df,
aes(Fecha, O3, color = "Prueba"),
linewidth = 0.8) +
geom_vline(xintercept = as.numeric(max(train_df$Fecha)),
linetype = "dashed", color = "grey40", linewidth = 0.4) +
scale_color_manual(values = c("Entrenamiento" = "#1B998B",
"Prueba" = "#C0392B")) +
labs(title = "Partición de la serie diaria de O3",
x = NULL, y = expression(O[3]~"(µg/m"^3*")"),
color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", color = "#0B6E4F"))
ggplotly(p_ventanas, dynamicTicks = TRUE) %>%
layout(hovermode = "x unified",
legend = list(orientation = "h", y = 1.1))La elección de 21 días como horizonte de prueba obedece a que es un plazo razonable para pronósticos de calidad del aire, suficientemente largo para evaluar el desempeño, pero no tanto como para que la incertidumbre haga poco útiles las predicciones. Además, al tratarse de tres semanas completas, se cubren distintos patrones de días de semana y fines de semana, lo que permite verificar si el modelo captura variaciones asociadas a la actividad humana.
En el gráfico se observa la partición: la línea verde corresponde a los datos de entrenamiento (344 días) y la línea roja a los 21 días de prueba. La línea vertical punteada marca la frontera entre ambos conjuntos. Se aprecia que los valores de ozono en el periodo de prueba son relativamente bajos (la mayoría entre 15 y 36 µg/m³), mientras que en el entrenamiento hay mayor variabilidad, incluyendo el pico máximo del año cercano a 60 µg/m³. Esta diferencia es útil porque pone a prueba el modelo: si logra pronosticar bien un periodo con características distintas a las del entrenamiento.
4.2 Estacionariedad: prueba ADF y función de autocorrelación
Para verificar si la serie diaria de ozono es estacionaria, se aplicó la prueba de Dickey‑Fuller Aumentada (ADF). El estadístico obtenido fue de -3.10 y el p‑valor de 0.1115, que es superior al nivel de significancia del 5 % (0.05). Por lo tanto, no se rechaza la hipótesis nula de que la serie tiene una raíz unitaria, es decir, la serie no es estacionaria. Esto significa que la media y la varianza de los datos cambian a lo largo del tiempo, lo cual es común en variables ambientales con tendencias o ciclos.
# Prueba de Dickey-Fuller aumentada sobre la serie de entrenamiento
adf_original <- adf.test(ventana)
# --- Tabla con los resultados de la prueba ---
tabla_adf <- data.frame(
Prueba = "Dickey-Fuller Aumentada",
Estadístico = round(adf_original$statistic, 4),
`p-valor` = ifelse(adf_original$p.value < 0.001,
"< 0.001", round(adf_original$p.value, 4)),
`Rezagos` = adf_original$parameter,
check.names = FALSE
)
tabla_adf %>%
gt() %>%
tab_header(
title = md("**Tabla 2. Prueba de estacionariedad ADF**"),
subtitle = "H₀: raíz unitaria (serie no estacionaria). Si p < 0.05 se rechaza H₀."
) %>%
cols_label(
Estadístico = "Estad. ADF",
`p-valor` = "p‑valor",
`Rezagos` = "Lags"
) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(15),
heading.subtitle.font.size = px(12))| Tabla 2. Prueba de estacionariedad ADF | |||
| H₀: raíz unitaria (serie no estacionaria). Si p < 0.05 se rechaza H₀. | |||
| Prueba | Estad. ADF | p‑valor | Lags |
|---|---|---|---|
| Dickey-Fuller Aumentada | -3.1027 | 0.1115 | 6 |
# --- Gráficos ACF y PACF ---
p_acf <- ggAcf(ventana) + ggtitle("ACF — serie original") + theme_minimal(base_size = 11)
p_pacf <- ggPacf(ventana) + ggtitle("PACF — serie original") + theme_minimal(base_size = 11)
acf_plotly <- ggplotly(p_acf) %>% layout(hovermode = "x")
pacf_plotly <- ggplotly(p_pacf) %>% layout(hovermode = "x")
# Combinamos ambos en un solo gráfico
subplot(acf_plotly, pacf_plotly, nrows = 1, shareY = TRUE) %>%
layout(title = list(text = "Funciones de autocorrelación de la serie original",
font = list(color = "#0B6E4F", size = 14)),
margin = list(t = 50))Además de la prueba ADF, se examinaron las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie original. En la ACF se observa que las correlaciones disminuyen lentamente y se mantienen positivas durante varios rezagos, lo que confirma la falta de estacionariedad. Por su parte, la PACF muestra 4 picos significativos desde el rezago 1 hasta el rezago 4. Con base en estas evidencias, se decidió aplicar una diferenciación de primer orden (d = 1) para estabilizar la serie y hacerla estacionaria, tal como se muestra en la siguiente sección.
4.3 Diferenciación de la serie
Después de aplicar la diferenciación de primer orden (ΔO₃), se repitió la prueba de Dickey‑Fuller Aumentada. El estadístico obtenido fue de -9.94 y el p‑valor de 0.01, que es inferior al nivel de significancia del 5 %. Por lo tanto, se rechaza la hipótesis nula de raíz unitaria y se concluye que la serie diferenciada es estacionaria. Esto significa que ahora la media y la varianza se mantienen constantes a lo largo del tiempo, condición necesaria para identificar un modelo ARIMA.
# Primera diferencia
miserie <- diff(ventana)
# Prueba ADF sobre el vector numérico
adf_diff <- adf.test(na.omit(as.numeric(miserie)))
# Tabla de resultados ADF
tabla_adf_diff <- data.frame(
Prueba = "Dickey-Fuller Aumentada (diferenciada)",
Estadístico = round(adf_diff$statistic, 4),
`p-valor` = ifelse(adf_diff$p.value < 0.001,
"< 0.001", round(adf_diff$p.value, 4)),
Rezagos = adf_diff$parameter,
check.names = FALSE
)
tabla_adf_diff %>%
gt() %>%
tab_header(
title = md("**Tabla 3. Prueba de estacionariedad sobre ΔO₃**"),
subtitle = "H₀: raíz unitaria (no estacionaria). Se rechaza con p < 0.05."
) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(15),
heading.subtitle.font.size = px(12))| Tabla 3. Prueba de estacionariedad sobre ΔO₃ | |||
| H₀: raíz unitaria (no estacionaria). Se rechaza con p < 0.05. | |||
| Prueba | Estadístico | p-valor | Rezagos |
|---|---|---|---|
| Dickey-Fuller Aumentada (diferenciada) | -9.9416 | 0.01 | 6 |
# Gráfico: serie diferenciada + ACF
p_diff <- autoplot(miserie) +
ggtitle("Serie diferenciada (d = 1)") +
ylab("ΔO₃") +
theme_minimal(base_size = 11)
p_acf_diff <- ggAcf(miserie) +
ggtitle("ACF — serie diferenciada") +
theme_minimal(base_size = 11)
subplot(
ggplotly(p_diff) %>% layout(hovermode = "x"),
ggplotly(p_acf_diff) %>% layout(hovermode = "x"),
nrows = 1, shareY = FALSE
) %>%
layout(title = list(text = "Diferenciación de la serie O₃",
font = list(color = "#0B6E4F", size = 14)),
margin = list(t = 50))# ACF y PACF de la serie ya diferenciada
p_acf <- ggAcf(miserie) + ggtitle("ACF — ΔO₃") + theme_minimal(base_size = 11)
p_pacf <- ggPacf(miserie) + ggtitle("PACF — ΔO₃") + theme_minimal(base_size = 11)
acf_plotly <- ggplotly(p_acf) %>% layout(hovermode = "x")
pacf_plotly <- ggplotly(p_pacf) %>% layout(hovermode = "x")
subplot(acf_plotly, pacf_plotly, nrows = 1, shareY = TRUE) %>%
layout(title = list(text = "Funciones de autocorrelación de ΔO₃",
font = list(color = "#0B6E4F", size = 14)),
margin = list(t = 50))Las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie diferenciada también confirman la estacionariedad. En la ACF se observa un decaimiento rápido: las correlaciones son significativas solo en el primer rezago (lag 1) y luego caen dentro de las bandas de confianza. En la PACF también se ve un pico importante en el rezago 1 y después valores cercanos a cero. Este patrón ACF con corte en 1 y PACF con decaimiento o también corte en 1 es típico de un modelo de media móvil de orden 1 (MA(1)) o alternativamente de un modelo autorregresivo de orden 1 (AR(1)). Con esta información, se procede a estimar varios modelos candidatos, cuyos órdenes p y q se eligen con base en estos gráficos.
4.4 Identificación y estimación de modelos candidatos
Para seleccionar el mejor modelo, se combinaron dos enfoques. Primero, se utilizó la función auto.arima con una búsqueda exhaustiva (stepwise = FALSE, approximation = FALSE), la cual explora automáticamente múltiples combinaciones de órdenes autorregresivos (p, d, q) y estacionales (P, D, Q) y elige la que minimiza el AICc. En segundo lugar, con base en el análisis de las funciones ACF y PACF de la serie diferenciada (que sugerían órdenes bajos), se estimaron manualmente cinco modelos candidatos: ARIMA(1,1,1), ARIMA(2,1,1), ARIMA(1,1,2), ARIMA(2,1,2) y ARIMA(0,1,1). Todos los modelos incluyen una diferenciación de orden 1 (d=1), que ya había sido validada mediante la prueba ADF.
# Modelo sugerido automáticamente
modelo_auto <- auto.arima(ventana, seasonal = TRUE,
stepwise = FALSE, approximation = FALSE)
# Modelos candidatos estimados manualmente (metodología Box-Jenkins)
modelo1 <- Arima(ventana, order = c(1, 1, 1))
modelo2 <- Arima(ventana, order = c(2, 1, 1))
modelo3 <- Arima(ventana, order = c(1, 1, 2))
modelo4 <- Arima(ventana, order = c(2, 1, 2))
modelo5 <- Arima(ventana, order = c(0, 1, 1))
modelos <- list(Auto = modelo_auto, M1 = modelo1, M2 = modelo2,
M3 = modelo3, M4 = modelo4, M5 = modelo5)
# Etiqueta legible del orden ARIMA de cada modelo
etiqueta_arima <- function(m) {
ar <- arimaorder(m)
if (length(ar) == 7) {
sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[%d]",
ar[1], ar[2], ar[3], ar[4], ar[5], ar[6], ar[7])
} else {
sprintf("ARIMA(%d,%d,%d)", ar[1], ar[2], ar[3])
}
}
tabla_modelos <- data.frame(
Modelo = names(modelos),
Especificacion = sapply(modelos, etiqueta_arima),
AIC = round(sapply(modelos, AIC), 2),
AICc = round(sapply(modelos, function(m) m$aicc), 2),
BIC = round(sapply(modelos, BIC), 2),
row.names = NULL
)
tabla_modelos %>%
gt() %>%
tab_header(title = md("**Tabla 2. Comparación de modelos ARIMA candidatos**"),
subtitle = "Criterios de información (menor es mejor)") %>%
fmt_number(columns = c(AIC, AICc, BIC), decimals = 2) %>%
tab_style(style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())) %>%
tab_style(
style = list(cell_fill(color = "#E8F5F1"), cell_text(weight = "bold")),
locations = cells_body(rows = which.min(tabla_modelos$AICc))
) %>%
tab_options(table.font.size = px(13), heading.title.font.size = px(16))| Tabla 2. Comparación de modelos ARIMA candidatos | ||||
| Criterios de información (menor es mejor) | ||||
| Modelo | Especificacion | AIC | AICc | BIC |
|---|---|---|---|---|
| Auto | ARIMA(1,1,2) | 2,154.30 | 2,154.42 | 2,169.65 |
| M1 | ARIMA(1,1,1) | 2,156.43 | 2,156.50 | 2,167.94 |
| M2 | ARIMA(2,1,1) | 2,155.10 | 2,155.21 | 2,170.45 |
| M3 | ARIMA(1,1,2) | 2,154.30 | 2,154.42 | 2,169.65 |
| M4 | ARIMA(2,1,2) | 2,157.36 | 2,157.54 | 2,176.55 |
| M5 | ARIMA(0,1,1) | 2,156.33 | 2,156.37 | 2,164.01 |
# Seleccionamos el modelo con menor AICc
idx_mejor <- which.min(sapply(modelos, function(m) m$aicc))
modelo_sel <- modelos[[idx_mejor]]
nombre_sel <- etiqueta_arima(modelo_sel)La Tabla 2 compara los seis modelos según los criterios de información AIC, AICc y BIC. En todos los casos, el modelo con el valor más bajo de AICc (criterio recomendado para muestras pequeñas o moderadas) es el ARIMA(1,1,2), que aparece tanto en el modelo automático (Auto) como en el modelo manual M3. Sus valores AICc (2,154.42) son ligeramente inferiores a los del segundo mejor modelo, ARIMA(0,1,1) (2,156.37). La diferencia de aproximadamente 2 unidades en AICc indica que el ARIMA(1,1,2) es claramente preferible, pues ofrece un mejor equilibrio entre ajuste y complejidad. Por lo tanto, se selecciona este modelo para las siguientes etapas de diagnóstico y pronóstico.
4.5 Interpretación de los coeficientes
La Tabla 4 presenta los coeficientes estimados del modelo ARIMA(1,1,2) seleccionado. El modelo incluye un término autorregresivo de orden 1 (ar1 = 0.8027) y dos términos de media móvil de órdenes 1 y 2 (ma1 = -1.4296, ma2 = 0.4658). Todos los coeficientes son estadísticamente significativos, con p‑valores inferiores a 0.001. Esto significa que cada uno de estos parámetros aporta información relevante para explicar la dinámica temporal del ozono.
coefs <- coef(modelo_sel)
if (length(coefs) > 0) {
se <- sqrt(diag(modelo_sel$var.coef))
zval <- coefs / se
pval <- 2 * (1 - pnorm(abs(zval)))
tabla_coef <- data.frame(
Termino = names(coefs),
Coeficiente = round(coefs, 4),
Error_Est = round(se, 4),
z = round(zval, 3),
p_valor = ifelse(pval < 0.001, "< 0.001", sprintf("%.4f", pval)),
row.names = NULL
)
tabla_coef %>%
gt() %>%
tab_header(title = md(paste0("**Tabla 4. Coeficientes del modelo seleccionado — ",
nombre_sel, "**"))) %>%
cols_label(Error_Est = "Error Est.", p_valor = "p-valor") %>%
tab_style(style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())) %>%
tab_options(table.font.size = px(13), heading.title.font.size = px(15))
} else {
cat("El modelo seleccionado no contiene coeficientes estimados",
"(p. ej. un paseo aleatorio ARIMA(0,1,0)).\n")
}| Tabla 4. Coeficientes del modelo seleccionado — ARIMA(1,1,2) | ||||
| Termino | Coeficiente | Error Est. | z | p-valor |
|---|---|---|---|---|
| ar1 | 0.8027 | 0.1464 | 5.483 | < 0.001 |
| ma1 | -1.4296 | 0.1685 | -8.482 | < 0.001 |
| ma2 | 0.4658 | 0.1360 | 3.425 | < 0.001 |
El coeficiente autorregresivo (ar1 = 0.80) indica una fuerte persistencia: el nivel de ozono de un día depende positivamente del nivel del día anterior. Los coeficientes de media móvil (ma1 negativo y ma2 positivo) reflejan cómo los errores de pronóstico de los días recientes afectan el valor actual. En particular, un ma1 de -1.43 sugiere que un error grande en el día anterior se traduce en una corrección en sentido contrario al día siguiente, mientras que el ma2 positivo (0.47) indica que el error de hace dos días tiene un efecto más atenuado pero también significativo. En conjunto, estos coeficientes confirman que el modelo captura adecuadamente la estructura de correlación de la serie, tal como se validará en las pruebas diagnósticas sobre los residuales.
4.6 Desempeño predictivo sobre la ventana de prueba
La Tabla 5 muestra las métricas de error del modelo ARIMA(1,1,2) tanto en el conjunto de entrenamiento (344 días) como en el de prueba (21 días). En entrenamiento, el RMSE es de 5.51 µg/m³ y el MAE de 4.26 µg/m³, lo que indica que, en promedio, el modelo se equivoca por unos 4‑5 µg/m³ al predecir la concentración diaria de ozono. El MAPE es del 16.0 %, un error relativo aceptable para una variable ambiental con alta variabilidad diaria. Además, el MASE (0.656) es inferior a 1, lo que significa que el modelo pronostica mejor que un modelo ingenuo (que simplemente repite el valor del día anterior). El sesgo (ME) es prácticamente nulo (0.05), indicando que no hay sobreestimación ni subestimación sistemática en el entrenamiento.
# Pronóstico del modelo de entrenamiento sobre el horizonte de prueba
fc_test <- forecast(modelo_sel, h = h, level = c(80, 95))
# Métricas de exactitud comparando contra la ventana de prueba
acc <- accuracy(fc_test, ventana2)
acc_df <- as.data.frame(round(acc, 3))
acc_df <- cbind(Conjunto = rownames(acc_df), acc_df)
rownames(acc_df) <- NULL
acc_df %>%
gt() %>%
tab_header(title = md("**Tabla 4. Métricas de exactitud del modelo seleccionado**"),
subtitle = "Entrenamiento vs. prueba (últimos 21 días)") %>%
tab_style(style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())) %>%
tab_options(table.font.size = px(12), heading.title.font.size = px(15))| Tabla 4. Métricas de exactitud del modelo seleccionado | ||||||||
| Entrenamiento vs. prueba (últimos 21 días) | ||||||||
| Conjunto | ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil's U |
|---|---|---|---|---|---|---|---|---|
| Training set | 0.05 | 5.513 | 4.258 | -3.618 | 15.995 | 0.656 | -0.001 | NA |
| Test set | 2.95 | 6.598 | 5.488 | 5.519 | 19.307 | 0.846 | 0.466 | 0.825 |
En la ventana de prueba, las métricas se degradan ligeramente, como es normal al enfrentar datos no vistos. El RMSE sube a 6.60 µg/m³, el MAE a 5.49 µg/m³ y el MAPE al 19.3 %. El sesgo (ME) es de 2.95, lo que sugiere una leve subestimación del modelo en los últimos 21 días del año (el modelo predice valores un poco más bajos de lo que realmente ocurrieron). El MASE de 0.846 sigue siendo inferior a 1, confirmando que el modelo aún supera al ingenuo. El coeficiente de Theil (0.825) indica que el pronóstico es útil, aunque con margen de mejora. En conjunto, las métricas de entrenamiento y prueba son de magnitud comparable, lo que descarta un sobreajuste severo y valida la capacidad predictiva del modelo.
5 Pronóstico
Una vez validado el modelo ARIMA(1,1,2) con los datos de entrenamiento, se reestimó sobre la serie completa de 2018 (365 días) y se generó un pronóstico para los primeros 21 días de 2019 (del 1 al 21 de enero). La Tabla 5 presenta los valores puntuales pronosticados junto con sus intervalos de confianza al 80 % y al 95 %. El pronóstico para el 1 de enero es de 33.3 µg/m³ y desciende gradualmente hasta estabilizarse en 31.4 µg/m³ a partir del 14 de enero, manteniéndose prácticamente constante durante el resto del horizonte. Este comportamiento refleja la dinámica del modelo: la serie tiende a un nivel de equilibrio (media incondicional) después de unos días, lo cual es consistente con un proceso estacionario en diferencias con media no nula.
# Reestimamos el modelo seleccionado sobre TODA la serie
serie_full <- ts(diario$O3, frequency = 7)
orden_sel <- arimaorder(modelo_sel)
if (length(orden_sel) == 7) {
modelo_final <- Arima(serie_full,
order = orden_sel[1:3],
seasonal = list(order = orden_sel[4:6],
period = orden_sel[7]))
} else {
modelo_final <- Arima(serie_full, order = orden_sel[1:3])
}
fc_fut <- forecast(modelo_final, h = h, level = c(80, 95))
# Fechas futuras asociadas al pronóstico
fechas_fut <- max(diario$Fecha) + 1:h
tabla_fc <- data.frame(
Fecha = fechas_fut,
Pronostico = round(as.numeric(fc_fut$mean), 2),
Lo80 = round(as.numeric(fc_fut$lower[, 1]), 2),
Hi80 = round(as.numeric(fc_fut$upper[, 1]), 2),
Lo95 = round(as.numeric(fc_fut$lower[, 2]), 2),
Hi95 = round(as.numeric(fc_fut$upper[, 2]), 2)
)
datatable(tabla_fc,
rownames = FALSE,
colnames = c("Fecha", "Pronóstico (µg/m³)", "Inf. 80%", "Sup. 80%",
"Inf. 95%", "Sup. 95%"),
caption = "Tabla 5. Pronóstico de O3 a 21 días con bandas de confianza",
options = list(pageLength = 10, dom = "tip")) %>%
formatStyle(columns = "Pronostico", fontWeight = "bold")Las bandas de confianza se amplían ligeramente a medida que el horizonte de pronóstico aumenta, reflejando la incertidumbre creciente propia de cualquier predicción. Por ejemplo, para el 1 de enero, el intervalo del 95 % va de 22.4 a 44.1 µg/m³, mientras que para el 21 de enero se extiende de 15.5 a 47.4 µg/m³. Esta amplitud (aproximadamente ±12 µg/m³ alrededor del pronóstico puntual) indica que, aunque el modelo ofrece una estimación central útil, no se puede descartar que el ozono real se desvíe considerablemente, especialmente en días con condiciones meteorológicas extremas. No obstante, todos los valores pronosticados se mantienen dentro del rango histórico observado en 2018 (12‑59 µg/m³) y por debajo del umbral de alerta sanitaria de la OMS (100 µg/m³ como media de 8 horas). El pronóstico sugiere que durante las primeras tres semanas de 2019 se esperan concentraciones típicas de la temporada seca de comienzos de año.
autoplot(fc_fut, include = 90) +
labs(title = paste0("Pronóstico de O3 a 21 días — ", nombre_sel),
subtitle = "Últimos 90 días observados + pronóstico con bandas al 80% y 95%",
x = "Tiempo (semanas)", y = expression(O[3]~"(µg/m"^3*")")) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))
6 Pruebas sobre los residuales
La Tabla 6 presenta los resultados de dos pruebas fundamentales para validar el modelo ARIMA(1,1,2). La prueba de Ljung‑Box, aplicada con 14 rezagos, arroja un estadístico de 8.01 y un p‑valor de 0.7126, muy superior al nivel de significancia del 5 %. Por lo tanto, no se rechaza la hipótesis nula de que los residuales no están autocorrelacionados, es decir, se comportan como ruido blanco en cuanto a su estructura temporal. La prueba de Shapiro‑Wilk para normalidad da un p‑valor de 0.2111 (> 0.05), lo que también lleva a no rechazar la normalidad de los residuales. Ambas pruebas confirman que el modelo ha capturado adecuadamente la información relevante de la serie
Ljung-Box test
data: Residuals from ARIMA(1,1,2)
Q* = 8.0076, df = 11, p-value = 0.7126
Model df: 3. Total lags used: 14res <- residuals(modelo_final)
# Pruebas
lb <- Box.test(res, lag = 14, fitdf = length(coef(modelo_final)), type = "Ljung-Box")
sw <- shapiro.test(as.numeric(res))
# Tabla de pruebas
tabla_tests <- data.frame(
Prueba = c("Ljung-Box (lag = 14)", "Shapiro-Wilk"),
Estadistico = round(c(lb$statistic, sw$statistic), 4),
p_valor = c(ifelse(lb$p.value < 0.001, "< 0.001", sprintf("%.4f", lb$p.value)),
ifelse(sw$p.value < 0.001, "< 0.001", sprintf("%.4f", sw$p.value))),
Hipotesis = c("H0: residuales no autocorrelacionados (ruido blanco)",
"H0: residuales con distribución normal"),
Conclusion = c(ifelse(lb$p.value > 0.05,
"No se rechaza H0: sin autocorrelación",
"Se rechaza H0: autocorrelación remanente"),
ifelse(sw$p.value > 0.05,
"No se rechaza H0: normalidad",
"Se rechaza H0: no normalidad")),
row.names = NULL
)
tabla_tests %>%
gt() %>%
tab_header(title = md("**Tabla 6. Pruebas diagnósticas sobre los residuales**")) %>%
cols_label(p_valor = "p-valor") %>%
tab_style(style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())) %>%
tab_options(table.font.size = px(12), heading.title.font.size = px(15))| Tabla 6. Pruebas diagnósticas sobre los residuales | ||||
| Prueba | Estadistico | p-valor | Hipotesis | Conclusion |
|---|---|---|---|---|
| Ljung-Box (lag = 14) | 8.0076 | 0.7126 | H0: residuales no autocorrelacionados (ruido blanco) | No se rechaza H0: sin autocorrelación |
| Shapiro-Wilk | 0.9945 | 0.2111 | H0: residuales con distribución normal | No se rechaza H0: normalidad |
Los gráficos complementarios refuerzan estas conclusiones. En el panel superior izquierdo se observa que los residuales fluctúan alrededor de cero sin patrones evidentes (media cercana a cero y varianza aparentemente constante).
# Cálculo de estadísticos
estadisticos_res <- data.frame(
Estadístico = c("Media", "Desviación estándar", "Asimetría", "Curtosis"),
Valor = c(
round(mean(res, na.rm = TRUE), 4),
round(sd(res, na.rm = TRUE), 4),
round(skewness(res, na.rm = TRUE), 4),
round(kurtosis(res, na.rm = TRUE), 4)
),
Ideal = c("0", "—", "0", "3"),
stringsAsFactors = FALSE
)
estadisticos_res %>%
gt() %>%
tab_header(
title = md("**Estadísticos descriptivos de los residuales**"),
subtitle = "Ruido blanco ideal: media 0, asimetría 0, curtosis 3"
) %>%
cols_label(
Ideal = "Valor ideal"
) %>%
tab_style(
style = list(cell_fill(color = verde$profundo),
cell_text(color = verde$texto, weight = "bold")),
locations = cells_column_labels(everything())
) %>%
tab_options(table.font.size = px(13),
heading.title.font.size = px(15),
heading.subtitle.font.size = px(12))| Estadísticos descriptivos de los residuales | ||
| Ruido blanco ideal: media 0, asimetría 0, curtosis 3 | ||
| Estadístico | Valor | Valor ideal |
|---|---|---|
| Media | 0.1194 | 0 |
| Desviación estándar | 5.5236 | — |
| Asimetría | -0.1469 | 0 |
| Curtosis | 3.5279 | 3 |
Los estadísticos descriptivos de los residuales del modelo ARIMA(1,1,2) complementan las pruebas diagnósticas. La media es de 0.1194, muy cercana a cero, lo que indica que el modelo no presenta sesgo sistemático. La desviación estándar (5.52) es coherente con el RMSE del entrenamiento (5.51 µg/m³). La asimetría es ligeramente negativa (-0.1469), lo que sugiere una cola izquierda apenas más larga que la derecha, pero cercana a cero. La curtosis es de 3.53, ligeramente superior al valor ideal de 3 (distribución normal), lo que indica colas un poco más pesadas —algo común en series ambientales con episodios extremos ocasionales. En conjunto, estos valores confirman que los residuales se aproximan razonablemente a una distribución normal y a ruido blanco, validando el modelo.
# Gráficos ACF y PACF
p_acf_res <- ggAcf(res) + ggtitle("ACF de los residuales") + theme_minimal(base_size = 11)
p_pacf_res <- ggPacf(res) + ggtitle("PACF de los residuales") + theme_minimal(base_size = 11)
acf_res_plotly <- ggplotly(p_acf_res) %>% layout(hovermode = "x")
pacf_res_plotly <- ggplotly(p_pacf_res) %>% layout(hovermode = "x")
subplot(acf_res_plotly, pacf_res_plotly, nrows = 1, shareY = TRUE) %>%
layout(
title = list(text = "Autocorrelación de los residuales del modelo",
font = list(color = "#0B6E4F", size = 14)),
margin = list(t = 50)
)La función de autocorrelación (ACF) de los residuales muestra todos los rezagos dentro de las bandas de confianza, confirmando la ausencia de correlación serial. Finalmente, el histograma con la curva de densidad normal superpuesta revela una distribución simétrica y acampanada, consistente con la normalidad. En conjunto, estas evidencias indican que el modelo ARIMA(1,1,2) es adecuado para pronosticar la concentración de ozono, cumpliendo con los supuestos de la metodología Box‑Jenkins.
7 Gráfico final
El siguiente gráfico, listo para incluir en el informe, resume en una sola figura los elementos esenciales del análisis: la serie original, su tendencia de fondo, la componente estacional (extraída por descomposición STL) y el pronóstico con sus bandas de confianza.
# --- Descomposición STL para extraer tendencia y estacionalidad ---
stl_fit <- stl(serie_full, s.window = "periodic")
comp <- as.data.frame(stl_fit$time.series)
diario$Tendencia <- as.numeric(comp$trend)
diario$Estacional <- as.numeric(comp$seasonal)
# --- Data frame del pronóstico con fechas reales ---
df_fc <- data.frame(
Fecha = fechas_fut,
mean = as.numeric(fc_fut$mean),
lo80 = as.numeric(fc_fut$lower[, 1]), hi80 = as.numeric(fc_fut$upper[, 1]),
lo95 = as.numeric(fc_fut$lower[, 2]), hi95 = as.numeric(fc_fut$upper[, 2])
)
col_serie <- "#1B998B"; col_trend <- "#0B6E4F"
col_fc <- "#C0392B"; col_band95 <- "#E59866"; col_band80 <- "#D98880"
# Panel superior: serie + tendencia + pronóstico con bandas
p_main <- ggplot() +
geom_line(data = diario, aes(Fecha, O3, color = "Serie observada"),
linewidth = 0.45, alpha = 0.85) +
geom_line(data = diario, aes(Fecha, Tendencia, color = "Tendencia (STL)"),
linewidth = 1.2) +
geom_ribbon(data = df_fc, aes(Fecha, ymin = lo95, ymax = hi95),
fill = col_band95, alpha = 0.35) +
geom_ribbon(data = df_fc, aes(Fecha, ymin = lo80, ymax = hi80),
fill = col_band80, alpha = 0.45) +
geom_line(data = df_fc, aes(Fecha, mean, color = "Pronóstico"),
linewidth = 1.1) +
geom_vline(xintercept = as.numeric(max(diario$Fecha)),
linetype = "dashed", color = "grey40") +
scale_color_manual(values = c("Serie observada" = col_serie,
"Tendencia (STL)" = col_trend,
"Pronóstico" = col_fc)) +
scale_x_date(date_breaks = "2 months", date_labels = "%b %Y") +
labs(title = "Ozono troposférico (O3) en Cali: serie, tendencia y pronóstico",
subtitle = paste0("Modelo ", nombre_sel,
" · bandas de confianza al 80% y 95% · pronóstico a 21 días"),
x = NULL, y = expression(O[3]~"(µg/m"^3*")"), color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
plot.title = element_text(face = "bold", color = col_trend, size = 14),
panel.grid.minor = element_blank())
# Panel inferior: componente estacional (un ciclo semanal representativo)
p_season <- ggplot(diario[1:35, ], aes(Fecha, Estacional)) +
geom_hline(yintercept = 0, color = "grey60") +
geom_line(color = col_trend, linewidth = 0.9) +
geom_area(fill = col_serie, alpha = 0.20) +
labs(title = "Componente estacional (semanal) extraída por STL",
x = NULL, y = "Efecto estacional") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", color = col_trend, size = 11))
# Combinamos ambos paneles en una sola figura
g_final <- arrangeGrob(p_main, p_season, ncol = 1, heights = c(2.2, 1))
grid::grid.newpage()
grid::grid.draw(g_final)
8 Conclusiones
La serie horaria revela un marcado ciclo diurno de origen fotoquímico —mínimos en la madrugada y máximos hacia el mediodía—, mientras que los promedios diarios reflejan la modulación climática bimodal del valle del Cauca. Las concentraciones más elevadas se concentran en los meses secos (ene y nov), mientras que los periodos lluviosos (nov y ene) registran los valores más bajos. Esta caracterización descriptiva ya constituye un insumo valioso para la vigilancia de la calidad del aire en la ciudad.
La serie diaria requirió una diferenciación de primer orden para
alcanzar la estacionariedad (ADF: p < 0.01). La comparación
sistemática de modelos mediante AIC, AICc y BIC, junto con la sugerencia
de auto.arima, condujo al modelo
ARIMA(1,1,2) (AICc = 2154.42). Las pruebas de
diagnóstico sobre los residuales —Ljung‑Box (p = 0.7126) y Shapiro‑Wilk
(p = 0.2111)— confirmaron que el modelo captura la estructura de
dependencia temporal, comportándose los residuales aproximadamente como
ruido blanco en términos de autocorrelación, aunque con colas
ligeramente más pesadas que las de una normal, algo habitual en series
ambientales.
El modelo mostró una capacidad predictiva razonable sobre la ventana de prueba de 21 días, con métricas de error comparables entre entrenamiento y prueba (MAPE en prueba: 19.3 %), lo que descarta un sobreajuste importante. Los pronósticos a tres semanas, acompañados de bandas de confianza al 80 % y 95 %, ofrecen un horizonte de anticipación útil para la planeación de medidas preventivas. No obstante, la presencia de variabilidad extrema —cambios interdiarios de hasta 28.6 µg/m³— advierte que los episodios pico pueden superar las bandas estimadas por un modelo lineal univariado.
El modelo ARIMA empleado ignora variables meteorológicas
determinantes (radiación solar, temperatura, humedad, viento), por lo
que una extensión natural sería un modelo ARIMAX o de
regresión dinámica que incorpore dichos predictores. Además, la riqueza
de los datos horarios originales —con una fuerte estacionalidad diaria—
invita a explorar modelos de estacionalidad múltiple (p. ej.,
TBATS o SARIMA horario‑semanal) para aprovechar plenamente
la resolución temporal. En conjunto, el trabajo demuestra cómo las
técnicas de series de tiempo convierten un registro ambiental en
conocimiento accionable para la gestión de la calidad del aire en
Cali.
9 Bibliografía
Box, G. E. P., & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day.
Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal‑trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3‑73.
DAGMA. (2019). Informe de calidad del aire del municipio de Santiago de Cali. Departamento Administrativo de Gestión del Medio Ambiente, Alcaldía de Santiago de Cali.
EPA. (2023). Ground-level Ozone Basics. United States Environmental Protection Agency. https://www.epa.gov/ground-level-ozone-pollution
Hyndman, R. J., & Athanasopoulos, G. (2018). Forecasting: Principles and Practice (2.ª ed.). OTexts. https://otexts.com/fpp2/
Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1‑22.
Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297‑303.
Organización Mundial de la Salud (OMS). (2021). Directrices mundiales de la OMS sobre la calidad del aire: partículas (PM2.5 y PM10), ozono, dióxido de nitrógeno, dióxido de azufre y monóxido de carbono. OMS.
R Core Team. (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.R-project.org/
Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591‑611.
Sievert, C. (2020). Interactive web‑based data visualization with R, plotly, and shiny. CRC Press.