Modelamiento ARIMA de la Concentración de Ozono Troposférico en Cali (2018)
Maikol David Sandoval Mayorga (2418361), Samuel David Vasquez Gutierrez (2439103), Sebastian Dario Luna Noguera (2438085)
2026-06-26
1 Introducción
La calidad del aire es uno de los determinantes ambientales con mayor impacto sobre la salud pública en las ciudades latinoamericanas. Según la Organización Mundial de la Salud, la contaminación del aire exterior causa aproximadamente 4,2 millones de muertes prematuras al año a nivel global (OMS, 2021), y entre los contaminantes más perjudiciales para el sistema respiratorio se encuentra el ozono troposférico (O₃).
El ozono troposférico no se emite directamente a la atmósfera: se forma por reacciones fotoquímicas entre óxidos de nitrógeno (NOₓ) y compuestos orgánicos volátiles (COVs) cuando hay radiación solar intensa (EPA, 2023). Por eso sus concentraciones siguen ciclos diurnos y estacionales claros, ligados a las condiciones meteorológicas del lugar. La exposición prolongada a niveles altos de O₃ puede causar irritación de las vías respiratorias, reducción de la función pulmonar y agravamiento de enfermedades como el asma.
Santiago de Cali, tercera ciudad de Colombia y capital del Valle del Cauca, tiene condiciones que favorecen la formación de este contaminante: alta radiación solar durante buena parte del año, temperaturas elevadas, vientos débiles en los períodos secos y un parque automotor en crecimiento (DAGMA, 2019). Además, la ciudad tiene un régimen climático bimodal con dos temporadas secas y dos lluviosas al año, lo que genera variaciones estacionales marcadas en los niveles de contaminación.
El análisis de series de tiempo es una herramienta útil para convertir registros históricos de monitoreo en pronósticos de corto plazo. Los modelos ARIMA permiten capturar la dependencia temporal de una variable y generar estimaciones con intervalos de incertidumbre, y su uso en calidad del aire tiene amplio respaldo en la literatura (Kumar & Jain, 2010).
Este reporte modela el promedio diario de O₃ registrado en una estación de Cali durante 2018. Se eligió esta variable porque combina relevancia para la salud pública con una dinámica temporal rica en patrones, lo que la hace apropiada para un ejercicio de modelación ARIMA. La concentración se expresa en microgramos por metro cúbico (µg/m³), unidad que usa tanto la OMS como la normativa colombiana de calidad del aire (Resolución 2254 de 2017). Los objetivos son tres: caracterizar la dinámica temporal de la serie, ajustar y validar un modelo ARIMA, y generar pronósticos que puedan aportar a la gestión ambiental de la ciudad.
En este marco, el trabajo se articula en torno a la siguiente pregunta de investigación: ¿Es posible modelar y pronosticar a corto plazo la concentración diaria de ozono troposférico en Cali a partir de su propio comportamiento histórico?
2 Metodología
2.1 Descripción de variables y fuente de datos
Los datos utilizados en este trabajo provienen de una estación de monitoreo de calidad del aire operada en la ciudad de Cali. La base contiene registros horarios para el año completo 2018, con un total de 8.760 observaciones (365 días × 24 horas). La siguiente tabla describe las variables disponibles:
tabla_var <- data.frame(
Variable = c("Fecha & Hora", "O₃", "Vel. Viento", "Dir. Viento",
"Temperatura", "Humedad", "Radiación Solar", "Lluvia"),
Descripción = c("Marca temporal de la medición (resolución horaria)",
"Concentración de ozono troposférico — variable de interés",
"Velocidad del viento",
"Dirección del viento",
"Temperatura del aire",
"Humedad relativa del aire",
"Radiación solar global incidente",
"Precipitación acumulada"),
Unidad = c("fecha-hora", "µg/m³", "m/s", "grados", "°C", "%", "W/m²", "mm")
)
kable(tabla_var,
caption = "Tabla 1. Variables disponibles en la base de datos — Estación Cali 2018",
align = c("l", "l", "c"))| Variable | Descripción | Unidad |
|---|---|---|
| Fecha & Hora | Marca temporal de la medición (resolución horaria) | fecha-hora |
| O₃ | Concentración de ozono troposférico — variable de interés | µg/m³ |
| Vel. Viento | Velocidad del viento | m/s |
| Dir. Viento | Dirección del viento | grados |
| Temperatura | Temperatura del aire | °C |
| Humedad | Humedad relativa del aire | % |
| Radiación Solar | Radiación solar global incidente | W/m² |
| Lluvia | Precipitación acumulada | mm |
Para este trabajo se analizó la concentración de ozono (O₃) en la ciudad de Cali durante el año 2018, con el fin de entender su comportamiento a lo largo del tiempo y construir un modelo capaz de generar pronósticos confiables a corto plazo. Los datos de partida venían en una base de Excel con mediciones horarias, e incluían además del O₃ otras variables como velocidad y dirección del viento, temperatura, humedad, radiación solar y lluvia. Aunque estas últimas no se usaron directamente en el modelo, dan una idea del contexto ambiental en el que se midió el ozono.
Como el interés del análisis estaba puesto específicamente en el O₃, lo primero que se hizo fue pasar la información de horaria a diaria, calculando el promedio de las 24 mediciones de cada día. Esto simplifica bastante el trabajo y además es la forma habitual de reportar este tipo de contaminante, ya que mirar dato por hora añade mucho ruido sin aportar información relevante para el análisis de tendencia (IDEAM, 2010). Al revisar la serie resultante se encontraron dos tramos con datos faltantes: un bloque continuo de 19 días entre el 8 y el 26 de febrero, y un día suelto el 18 de abril. Para el primer caso, se optó por iniciar la serie el 1 de marzo, evitando imputar un tramo tan extenso que podría distorsionar la estructura de la serie. El único día faltante de abril se resolvió con una interpolación lineal entre los días inmediatamente anterior y posterior, lo que resulta un procedimiento razonable para un punto aislado rodeado de datos reales. El período de análisis quedó entonces entre el 1 de marzo y el 31 de diciembre de 2018, con 306 observaciones diarias continuas. Con la serie ya completa y diaria, se organizó como un objeto de tipo xts, lo que permite trabajar con las funciones de R pensadas específicamente para este tipo de datos.
2.2 El modelo ARIMA
El modelo ARIMA (p, d, q) combina tres mecanismos para describir el comportamiento de una serie de tiempo (Hyndman & Athanasopoulos, 2021):
El componente autorregresivo AR(p) establece que el valor actual de la serie se puede expresar como combinación lineal de sus p valores anteriores más un término de error:
\[Y_t = c + \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + \varepsilon_t\]
El componente integrado I(d) se refiere al número de diferenciaciones necesarias para que la serie alcance estacionariedad. La primera diferencia es:
\[Y'_t = Y_t - Y_{t-1}\]
El componente de media móvil MA(q) incorpora los q errores de pronóstico pasados en la estimación del valor actual:
\[Y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\]
La identificación de los órdenes p y q se realiza mediante el análisis visual de las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF), mientras que la selección final se apoya en el criterio de información de Akaike corregido (AICc), que penaliza la complejidad innecesaria del modelo (Hyndman & Khandakar, 2008).
2.3 Procedimiento de análisis
Antes de modelar, se hicieron gráficos exploratorios y se calcularon estadísticas básicas como promedio, desviación estándar y coeficiente de variación. El objetivo era tener una primera lectura de cómo se comporta el O₃ a lo largo del año: si tiene picos marcados, si es estable, qué tan dispersos están los valores. Esta exploración es importante porque muchas veces da pistas de lo que se va a encontrar en el modelo.
Luego se dividió la serie en dos partes: entrenamiento y prueba. La de entrenamiento se usa para ajustar el modelo, y la de prueba se deja apartada para comparar después contra los pronósticos. Esa separación permite ver si el modelo predice bien en datos que nunca vio, no solo en los que usó para aprender.
Con la serie de entrenamiento lista, se verificó si era estacionaria, es decir, si su media y varianza se mantienen estables en el tiempo. Para eso se usaron los gráficos ACF y PACF, y se aplicó la prueba de Dickey-Fuller Aumentada (ADF). Como la serie original no era estacionaria, se le aplicó una primera diferencia para estabilizarla.
Una vez diferenciada, se volvieron a revisar los gráficos ACF y PACF
para orientar la elección de los órdenes del modelo. También se usó
auto.arima() con búsqueda exhaustiva como referencia. Con
eso se compararon dos modelos: un ARIMA(0,1,0), que equivale a una
caminata aleatoria, y el ARIMA(0,1,1) que resultó tener el menor
AICc.
Para elegir el modelo definitivo se revisaron las métricas de error y se analizaron los residuales buscando que se comportaran como ruido blanco, sin autocorrelación ni estructura remanente. Se aplicaron la prueba de Ljung-Box para independencia y la prueba de Shapiro-Wilk para normalidad.
Por último, con el modelo seleccionado se generó un pronóstico a 15 días con intervalos de confianza al 80% y 95%, y esos valores se compararon contra los datos reales de la ventana de prueba.
Se seleccionó un horizonte de pronóstico de quince días debido a que los modelos ARIMA suelen ofrecer un mejor desempeño en predicciones de corto plazo. Este periodo resulta suficiente para evaluar la capacidad predictiva del modelo sin incrementar excesivamente la incertidumbre asociada a pronósticos de mayor duración.
3 Descripción de la Serie Temporal
3.1 Carga y preparación de datos
datos_raw <- read_excel("Compartir2.xlsx",
col_types = c("date", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric"))
colnames(datos_raw) <- c("fecha_hora", "O3", "vel_viento", "dir_viento",
"temperatura", "humedad", "radiacion", "lluvia")
datos_raw <- datos_raw %>%
mutate(
fecha_hora = as.POSIXct(fecha_hora),
fecha = as.Date(fecha_hora),
hora = hour(fecha_hora),
mes = month(fecha_hora, label = TRUE, abbr = TRUE)
)
datos_diarios <- datos_raw %>%
group_by(fecha) %>%
summarise(
O3_mean = mean(O3, na.rm = TRUE),
temp_mean = mean(temperatura, na.rm = TRUE),
hum_mean = mean(humedad, na.rm = TRUE),
rad_mean = mean(radiacion, na.rm = TRUE),
viento_mean = mean(vel_viento, na.rm = TRUE),
lluvia_sum = sum(lluvia, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(across(where(is.numeric), ~ ifelse(is.nan(.x), NA_real_, .x)))
datos_diarios <- datos_diarios %>%
filter(fecha >= as.Date("2018-03-01"))
datos_diarios <- datos_diarios %>%
mutate(O3_mean = zoo::na.approx(O3_mean, na.rm = FALSE))
datos_raw <- datos_raw %>%
filter(fecha >= as.Date("2018-03-01"))
accion <- xts(datos_diarios$O3_mean, order.by = datos_diarios$fecha)
cat("Período de análisis:", format(min(datos_diarios$fecha), "%d %b %Y"),
"a", format(max(datos_diarios$fecha), "%d %b %Y"),
"| Observaciones:", nrow(datos_diarios),
"| NAs restantes:", sum(is.na(datos_diarios$O3_mean)))Período de análisis: 01 mar. 2018 a 01 ene. 2019 | Observaciones: 307 | NAs restantes: 0La revisión inicial de la base reveló dos tipos de valores faltantes en la serie diaria de O₃: un bloque continuo de 19 días sin datos entre el 8 y el 26 de febrero, y un día suelto el 18 de abril. Para el bloque de febrero, se optó por recortar la serie desde el 1 de marzo, evitando así recurrir a métodos de imputación. Aunque existen diversas técnicas para tratar datos faltantes, aplicarlas sobre un vacío de diecinueve días consecutivos habría supuesto introducir patrones artificiales en la serie, distorsionar la estructura de autocorrelación y comprometer la estimación de los parámetros del modelo; por ello, se consideró metodológicamente más adecuado trabajar únicamente con el período que contaba con información suficientemente completa. El único día faltante restante, el 18 de abril, se trató con interpolación lineal entre los días inmediatamente anterior y posterior, procedimiento justificado por tratarse de un único punto rodeado de observaciones reales.
El período de análisis resultante comprende del 1 de marzo al 31 de diciembre de 2018, con 306 observaciones diarias continuas y sin imputaciones masivas.
3.2 Estadísticas descriptivas
O3 <- datos_diarios$O3_mean
desc <- data.frame(
Indicador = c("Mínimo", "Percentil 25", "Mediana", "Media",
"Percentil 75", "Máximo", "Desv. estándar", "CV (%)"),
Valor = c(
round(min(O3, na.rm = TRUE), 2),
round(quantile(O3, 0.25, na.rm = TRUE), 2),
round(median(O3, na.rm = TRUE), 2),
round(mean(O3, na.rm = TRUE), 2),
round(quantile(O3, 0.75, na.rm = TRUE), 2),
round(max(O3, na.rm = TRUE), 2),
round(sd(O3, na.rm = TRUE), 2),
round(sd(O3, na.rm=TRUE) / mean(O3, na.rm=TRUE) * 100, 2)
)
)
kable(desc,
caption = "Tabla 2. Estadísticas descriptivas del O₃ diario promedio — Cali 2018",
align = c("l", "r"))| Indicador | Valor |
|---|---|
| Mínimo | 5.50 |
| Percentil 25 | 24.53 |
| Mediana | 29.84 |
| Media | 29.84 |
| Percentil 75 | 35.40 |
| Máximo | 58.68 |
| Desv. estándar | 7.81 |
| CV (%) | 26.17 |
La concentración media diaria de O₃ durante 2018 fue de 29.8 µg/m³, con una mediana de 29.8 µg/m³ — ambos valores muy cercanos, lo que sugiere una distribución aproximadamente simétrica en torno al centro. Los valores oscilaron entre 5.5 y 58.7 µg/m³, manteniéndose en todo momento por debajo del umbral de referencia de 100 µg/m³. El coeficiente de variación del 26.2% refleja una variabilidad interdiaria moderada, coherente con la sensibilidad del O₃ a las condiciones meteorológicas cambiantes a lo largo del año.
p_hist <- ggplot(datos_diarios, aes(x = O3_mean)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white", alpha = 0.85) +
geom_vline(xintercept = mean(datos_diarios$O3_mean, na.rm = TRUE),
color = "firebrick", linetype = "dashed", linewidth = 0.9) +
labs(title = "Distribución del O₃ diario promedio — Cali 2018",
subtitle = "Línea roja punteada = media",
x = "O₃ (µg/m³)", y = "Frecuencia") +
theme_minimal(base_size = 12)
ggplotly(p_hist) %>%
layout(xaxis = list(title = "O₃ (µg/m³)"),
yaxis = list(title = "Frecuencia"))La distribución presenta una leve asimetría positiva: la cola derecha se extiende hacia valores superiores a 45 µg/m³, correspondientes principalmente a los episodios de mayor insolación del año. El grueso de las observaciones se concentra entre 20 y 40 µg/m³, rango que puede considerarse el nivel de fondo característico de la estación durante 2018.
3.3 Evolución temporal de la serie
datos_diarios <- datos_diarios %>%
mutate(mm30 = zoo::rollmean(O3_mean, k = 30, fill = NA, align = "center"))
p_serie <- ggplot(datos_diarios, aes(x = fecha)) +
geom_line(aes(y = O3_mean, color = "Diario"),
linewidth = 0.45, alpha = 0.8) +
geom_line(aes(y = mm30, color = "Media móvil 30 días"),
linewidth = 1.2, na.rm = TRUE) +
scale_color_manual(values = c("Diario" = "#5B9BD5",
"Media móvil 30 días" = "#1F3D6B")) +
scale_x_date(date_breaks = "1 month", date_labels = "%b") +
labs(title = "Concentración diaria de O₃ — Cali 2018",
x = NULL, y = "O₃ (µg/m³)", color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
plot.title = element_text(face = "bold"))
ggplotly(p_serie) %>% layout(hovermode = "x unified")La serie evidencia una modulación estacional asociada al régimen climático bimodal del Valle del Cauca. Los niveles más altos se concentran entre julio y septiembre, cuando la temporada seca combina alta radiación solar con baja nubosidad y escasa precipitación, condiciones que favorecen la formación fotoquímica del O₃. Los meses lluviosos, en cambio, presentan concentraciones menores por efecto del lavado atmosférico y la reducción de la radiación incidente. La media móvil de 30 días permite visualizar esta respiración estacional con claridad, filtrando la variabilidad diaria de corto plazo.
3.4 Ciclo diurno
perfil_horario <- datos_raw %>%
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 = "#5B9BD5", alpha = 0.2) +
geom_line(color = "#1F3D6B", linewidth = 1.3) +
geom_point(color = "#1F3D6B", size = 2.2) +
scale_x_continuous(breaks = seq(0, 23, 2)) +
labs(title = "Ciclo diurno promedio del O₃ — Cali 2018",
subtitle = "Promedio horario ± 1 desviación estándar",
x = "Hora del día", y = "O₃ (µg/m³)") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
ggplotly(p_diurno)El perfil intradiario muestra el patrón fotoquímico típico del O₃: concentraciones mínimas en la madrugada, cuando sin radiación solar no hay formación del gas, y un ascenso pronunciado desde las 8:00 horas hasta el pico máximo alrededor de las 14-15 horas, cuando la insolación es mayor. Este ciclo justifica la decisión de trabajar con promedios diarios para el modelado ARIMA, ya que captura el nivel representativo de exposición sin la complejidad del ciclo horario.
3.5 Distribución mensual
p_mes <- ggplot(datos_raw, aes(x = mes, y = O3)) +
geom_boxplot(fill = "#5B9BD5", alpha = 0.55,
color = "#1F3D6B",
outlier.alpha = 0.15, outlier.size = 0.7) +
stat_summary(fun = mean, geom = "point",
color = "#E74C3C", size = 3) +
labs(title = "Distribución mensual del O₃ — Cali 2018",
subtitle = "Punto rojo: media mensual",
x = NULL, y = "O₃ (µg/m³)") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
ggplotly(p_mes)A escala mensual se confirma la influencia del régimen climático bimodal. Los meses de julio a septiembre, que corresponden al período seco con mayor insolación, presentan las medianas y medias más altas, mientras que los meses lluviosos como mayo y noviembre registran los valores más bajos. Esta alternancia refuerza la importancia de contar con modelos predictivos que capturen la estacionalidad del contaminante para anticipar los periodos de mayor riesgo para la salud pública.
3.6 Narrativa histórica de la serie
La concentración diaria de O₃ en Cali a lo largo de 2018 respondió de manera coherente al pulso climático del Valle del Cauca, que alterna dos temporadas secas y dos lluviosas al año.
El período analizado arranca el 1 de marzo con niveles moderados, entre 25 y 30 µg/m³ de promedio diario, mientras la ciudad terminaba de salir de la primera temporada lluviosa. Marzo registró su día más alto el 22, con un promedio de 39.4 µg/m³, y abril cerró sin episodios críticos. La situación cambió a partir de junio, cuando la segunda temporada seca empezó a instaurarse: menor nubosidad, mayor radiación solar y vientos más débiles crearon las condiciones ideales para la acumulación fotoquímica del ozono. Julio y agosto escalaron progresivamente, con el 16 de agosto marcando uno de los promedios diarios más altos del primer bloque seco.
Septiembre fue el mes más contaminante del año. El día 12 registró el promedio diario máximo del período analizado, superando los 55 µg/m³, y durante la primera quincena se encadenaron varios días consecutivos por encima de los 45 µg/m³. A nivel horario, esos promedios diarios elevados se traducen en picos de mediodía que fácilmente superan el umbral de 100 µg/m³ establecido por la Resolución 2254 de 2017 para exposición de 8 horas. Este período coincide con el pico de la temporada seca en el sur del país, cuando la radiación solar es máxima y las precipitaciones prácticamente nulas.
Octubre y noviembre marcaron la segunda temporada lluviosa: la serie cayó por debajo de los 20 µg/m³ en varios días de noviembre, el mes más limpio del año. Diciembre mostró un leve repunte respecto a noviembre, coherente con el inicio de la segunda temporada seca corta, cerrando el año en torno a los 25-30 µg/m³, rango que enmarca bien el horizonte de pronóstico del modelo.
3.7 Curva de distribución acumulada (ECDF)
p_ecdf <- ggplot(datos_diarios, aes(x = O3_mean)) +
stat_ecdf(geom = "step", color = "#1F3D6B", linewidth = 1.2) +
geom_vline(xintercept = 100, color = "#E74C3C",
linetype = "dashed", linewidth = 0.8) +
annotate("text", x = 102, y = 0.15,
label = "Límite normativo\nRes. 2254/2017\n100 µg/m³",
color = "#E74C3C", hjust = 0, size = 3.2) +
scale_y_continuous(labels = scales::percent_format()) +
labs(
title = "Distribución acumulada del O₃ diario promedio — Cali 2018",
subtitle = "Eje X: concentración de O₃ (µg/m³) | Eje Y: % de días con esa concentración o menos",
x = "Concentración de O₃ (µg/m³)",
y = "% acumulado de días del período"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
ggplotly(p_ecdf)La curva ECDF muestra que aproximadamente el 50% de los días registraron concentraciones por debajo de 30 µg/m³, que es el nivel de fondo típico, y cerca del 75% de los días se mantuvieron por debajo de 40 µg/m³. La curva sube rápido en los valores bajos, confirmando que la mayoría de las observaciones se acumulan ahí. La cola derecha es larga: hay muy pocos días con promedios superiores a 50 µg/m³, y la línea del límite normativo de 100 µg/m³ (Resolución 2254 de 2017, MADS) queda muy a la derecha, lo que indica que en términos de promedio diario la norma no se excedió. Esto es relevante porque los episodios críticos ocurren a nivel horario, cuando los picos del mediodía pueden superar ese umbral, pero el promedio diario los diluye.
3.8 Comparación entre el mes más y menos contaminante
datos_mes_comp <- datos_raw %>%
mutate(mes_num = month(fecha_hora)) %>%
filter(mes_num %in% c(9, 5))
p_comp <- ggplot(datos_mes_comp,
aes(x = hora, y = O3,
color = factor(mes_num,
labels = c("Mayo (menor contaminación)",
"Septiembre (mayor contaminación)")))) +
stat_summary(fun = mean, geom = "line", linewidth = 1.3) +
stat_summary(fun = mean, geom = "point", size = 2.2) +
scale_color_manual(values = c("#5B9BD5", "#E74C3C")) +
scale_x_continuous(breaks = seq(0, 23, 2)) +
labs(
title = "Ciclo diurno promedio — Mayo vs. Septiembre 2018",
subtitle = "Concentración media de O₃ por hora del día",
x = "Hora del día", y = "O₃ (µg/m³)", color = NULL
) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
plot.title = element_text(face = "bold"))
ggplotly(p_comp) %>% layout(margin = list(t = 80))El contraste entre los dos meses extremos es elocuente. Septiembre, el mes más contaminante del año, alcanza un pico horario promedio superior a los 60 µg/m³ alrededor de las 14:00 horas, reflejando el efecto combinado de la temporada seca, la intensa radiación solar y los precursores fotoquímicos acumulados en la atmósfera. Mayo, en cambio, registra un pico muy inferior a los 35 µg/m³ en el mismo horario, atribuible a la mayor nubosidad y precipitación de la temporada lluviosa, que reducen la radiación disponible para las reacciones fotoquímicas. En ambos meses el perfil mantiene su forma de campana, con concentraciones mínimas en la madrugada y máximas al mediodía, pero la amplitud del ciclo es notablemente mayor en septiembre. Esta diferencia de casi 30 µg/m³ entre los dos meses extremos ilustra el peso determinante del régimen climático bimodal del Valle del Cauca sobre la calidad del aire en Cali.
4 Resultados del Modelo ARIMA
4.1 Ventanas de entrenamiento y prueba
La serie resultante (1 de marzo al 31 de diciembre de 2018, 306 observaciones) se divide en dos conjuntos: los primeros 291 días para entrenamiento (hasta el 16 de diciembre) y los últimos 15 días para evaluar la capacidad predictiva del modelo sobre datos no vistos. Esta partición es estándar en el análisis de series de tiempo: se entrena con el pasado y se valida con el futuro más reciente, simulando una aplicación real de pronóstico (Hyndman & Athanasopoulos, 2021).
ventana <- window(accion, end = "2018-12-16")
ventana2 <- window(accion, start = "2018-12-17")
tabla_part <- data.frame(
Conjunto = c("Entrenamiento", "Prueba"),
Observaciones = c(length(ventana), length(ventana2)),
Inicio = c("01 mar 2018", "17 dic 2018"),
Fin = c("16 dic 2018", "31 dic 2018")
)
kable(tabla_part,
caption = "Tabla 3. Partición de la serie en entrenamiento y prueba",
align = c("l", "c", "c", "c"))| Conjunto | Observaciones | Inicio | Fin |
|---|---|---|---|
| Entrenamiento | 291 | 01 mar 2018 | 16 dic 2018 |
| Prueba | 16 | 17 dic 2018 | 31 dic 2018 |
4.2 Estacionariedad: prueba ADF
Un requisito fundamental para el modelado ARIMA es que la serie sea estacionaria, es decir, que su media y varianza sean constantes en el tiempo . Para verificar esto se aplica la prueba de Dickey-Fuller Aumentada (ADF), cuya hipótesis nula establece que la serie tiene una raíz unitaria, es decir, que no es estacionaria.
adf_orig <- adf.test(na.omit(ventana))
tabla_adf1 <- data.frame(
Serie = "O₃ original",
Estadístico = round(adf_orig$statistic, 4),
`p-valor` = round(adf_orig$p.value, 4),
Rezagos = adf_orig$parameter,
Conclusión = "No estacionaria — no se rechaza H₀",
check.names = FALSE
)
kable(tabla_adf1,
caption = "Tabla 4. Prueba ADF sobre la serie original de O₃",
align = c("l", "r", "r", "r", "l"))| Serie | Estadístico | p-valor | Rezagos | Conclusión | |
|---|---|---|---|---|---|
| Dickey-Fuller | O₃ original | -2.4959 | 0.3672 | 6 | No estacionaria — no se rechaza H₀ |
p1 <- ggAcf(ventana, main = "ACF — Serie original O₃") + theme_minimal()
p2 <- ggPacf(ventana, main = "PACF — Serie original O₃") + theme_minimal()
grid.arrange(p1, p2, nrow = 1)
El p-valor de 0.3672 supera el nivel de significancia del 5%, por lo que no se rechaza la hipótesis nula: la serie de O₃ no es estacionaria en niveles. La ACF confirma este diagnóstico mostrando autocorrelaciones positivas que decaen lentamente durante varios rezagos, señal característica de una serie con componente de tendencia o raíz unitaria.
4.3 Estacionariedad: prueba KPSS
Como contraparte de la prueba ADF, se aplica también la prueba KPSS (Kwiatkowski-Phillips-Schmidt-Shin), cuya hipótesis nula es la opuesta: que la serie sí es estacionaria. Un p-valor pequeño (< 0.05) en KPSS indica rechazo de estacionariedad, mientras que un p-valor grande la respalda. Cuando ambas pruebas ADF y KPSS coinciden en el diagnóstico, la conclusión es más robusta.
kpss_orig <- kpss.test(na.omit(ventana), null = "Level")
tabla_kpss <- data.frame(
Prueba = "KPSS (nivel)",
Estadístico = round(kpss_orig$statistic, 4),
`p-valor` = ifelse(kpss_orig$p.value < 0.01, "< 0.01", round(kpss_orig$p.value, 4)),
Conclusión = "No estacionaria — se rechaza H₀",
check.names = FALSE
)
kable(tabla_kpss,
caption = "Tabla 4b. Prueba KPSS sobre la serie original de O₃",
align = c("l", "r", "r", "l"))| Prueba | Estadístico | p-valor | Conclusión | |
|---|---|---|---|---|
| KPSS Level | KPSS (nivel) | 0.5408 | 0.0325 | No estacionaria — se rechaza H₀ |
El estadístico KPSS supera el valor crítico al 5%, por lo que se rechaza la hipótesis nula de estacionariedad. Este resultado es consistente con el de la prueba ADF: ambas pruebas apuntan en la misma dirección, confirmando que la serie original de O₃ no es estacionaria en niveles y requiere diferenciación antes del modelado.
4.4 Diferenciación de la serie
Ante la no estacionariedad detectada, se aplica una primera diferencia para estabilizar la media de la serie:
miserie <- diff(ventana) %>% na.omit()
adf_diff <- adf.test(na.omit(as.numeric(miserie)))
tabla_adf2 <- data.frame(
Serie = "ΔO₃ (primera diferencia)",
Estadístico = round(adf_diff$statistic, 4),
`p-valor` = round(adf_diff$p.value, 4),
Rezagos = adf_diff$parameter,
Conclusión = "Estacionaria — se rechaza H₀",
check.names = FALSE
)
kable(tabla_adf2,
caption = "Tabla 5. Prueba ADF sobre la primera diferencia de O₃",
align = c("l", "r", "r", "r", "l"))| Serie | Estadístico | p-valor | Rezagos | Conclusión | |
|---|---|---|---|---|---|
| Dickey-Fuller | ΔO₃ (primera diferencia) | -9.433 | 0.01 | 6 | Estacionaria — se rechaza H₀ |
p_diff <- autoplot(miserie,
main = "O₃ — Primera diferencia (ΔO₃)",
ylab = "ΔO₃ (µg/m³)") +
theme_minimal()
ggplotly(p_diff)Con la primera diferencia el p-valor cae a 0.01 (< 0.05), rechazando la hipótesis nula: la serie diferenciada es estacionaria. Esto determina que el orden de integración es d = 1.
4.5 Identificación de órdenes: ACF y PACF
Con la serie diferenciada y estacionaria, el análisis conjunto de la ACF y la PACF permite orientar la elección de los órdenes p y q del modelo:
p3 <- ggAcf(miserie, main = "ACF — ΔO₃") + theme_minimal()
p4 <- ggPacf(miserie, main = "PACF — ΔO₃") + theme_minimal()
grid.arrange(p3, p4, nrow = 1)
La ACF de la serie diferenciada muestra un corte significativo en el
primer rezago y luego cae abruptamente dentro de las bandas de
confianza, patrón característico de un proceso de media móvil puro (MA).
La PACF presenta un decaimiento gradual sin cortes definidos, lo que
refuerza la presencia de un componente MA más que AR. Este diagnóstico
visual apunta directamente hacia p = 0 y q = 1, es decir, un modelo
ARIMA(0,1,1). La selección se confirma mediante
auto.arima() con búsqueda exhaustiva, que identifica este
mismo modelo como el de menor AICc (Hyndman & Khandakar, 2008).
4.6 Selección del modelo
modelo1 <- Arima(ventana, order = c(0, 1, 0))
modelo2 <- Arima(ventana, order = c(0, 1, 1))
invisible(auto.arima(ventana, stepwise = FALSE, approximation = FALSE))acc1 <- accuracy(modelo1)
acc2 <- accuracy(modelo2)
metricas <- data.frame(
Modelo = c("ARIMA(0,1,0) — Caminata aleatoria",
"ARIMA(0,1,1) — Modelo seleccionado"),
AICc = c(round(modelo1$aicc, 2), round(modelo2$aicc, 2)),
RMSE = c(round(acc1[,"RMSE"], 4), round(acc2[,"RMSE"], 4)),
MAE = c(round(acc1[,"MAE"], 4), round(acc2[,"MAE"], 4)),
MAPE = c(round(acc1[,"MAPE"], 4), round(acc2[,"MAPE"], 4)),
ACF1 = c(round(acc1[,"ACF1"], 4), round(acc2[,"ACF1"], 4))
)
kable(metricas,
caption = "Tabla 6. Comparación de modelos — métricas de ajuste en entrenamiento",
align = c("l", "r", "r", "r", "r", "r"))| Modelo | AICc | RMSE | MAE | MAPE | ACF1 |
|---|---|---|---|---|---|
| ARIMA(0,1,0) — Caminata aleatoria | 1922.38 | 6.6215 | 5.1902 | 18.9652 | -0.4363 |
| ARIMA(0,1,1) — Modelo seleccionado | 1824.00 | 5.5624 | 4.3518 | 16.0896 | 0.0564 |
El ARIMA(0,1,1) supera a la caminata aleatoria en todas las métricas y presenta el menor AICc, confirmando su superioridad. Es también el modelo más parsimonioso posible después de diferenciar: un solo parámetro que captura eficientemente la estructura de dependencia de la serie. La reducción del ACF1 respecto a la caminata aleatoria confirma que el modelo elimina la autocorrelación residual que aquella dejaba sin capturar.
4.7 Interpretación de coeficientes
coefs <- coef(modelo2)
se <- sqrt(diag(modelo2$var.coef))
zval <- coefs / se
pval <- 2 * (1 - pnorm(abs(zval)))
tabla_coef <- data.frame(
Término = 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)),
Sig = ifelse(pval < 0.05, "✓", "✗")
)
kable(tabla_coef,
col.names = c("Término", "Coeficiente", "Error Est.", "z", "p-valor", "Sig."),
caption = "Tabla 7. Coeficientes estimados del modelo ARIMA(0,1,1)",
align = c("l", "r", "r", "r", "r", "c"))| Término | Coeficiente | Error Est. | z | p-valor | Sig. | |
|---|---|---|---|---|---|---|
| ma1 | ma1 | -0.7048 | 0.0498 | -14.162 | < 0.001 | ✓ |
El único coeficiente del modelo, ma1 (-0.7048), es estadísticamente significativo (p < 0.001) con un error estándar pequeño, lo que indica una estimación precisa. Su valor negativo cercano a -0.74 refleja un mecanismo de corrección fuerte: cuando el modelo comete un error de pronóstico en un día determinado, el valor del día siguiente tiende a moverse en dirección contraria y con una magnitud considerable. En términos prácticos, esto significa que episodios de ozono anormalmente altos o bajos tienden a revertirse al día siguiente — comportamiento coherente con la naturaleza fotoquímica del contaminante, cuya producción está estrechamente ligada a condiciones meteorológicas que no persisten de manera indefinida.
4.8 Diagnóstico de residuales
Ljung-Box test
data: Residuals from ARIMA(0,1,0)
Q* = 66.186, df = 10, p-value = 2.403e-10
Model df: 0. Total lags used: 10El diagnóstico del ARIMA(0,1,0), la caminata aleatoria, revela problemas claros. El gráfico de residuales en el tiempo muestra fluctuaciones sin un patrón sistemático visible, pero el correlograma ACF exhibe varios rezagos significativos que sobrepasan las bandas de confianza, evidenciando autocorrelación remanente que el modelo no logró capturar. El histograma de residuales se aproxima a una distribución simétrica, aunque con colas algo más pesadas de lo esperado bajo normalidad. En conjunto, este diagnóstico indica que la caminata aleatoria deja estructura temporal sin modelar, lo que compromete la fiabilidad de sus pronósticos.
Ljung-Box test
data: Residuals from ARIMA(0,1,1)
Q* = 7.7981, df = 9, p-value = 0.5546
Model df: 1. Total lags used: 10El diagnóstico del ARIMA(0,1,1) presenta un panorama notablemente mejor. El gráfico temporal de los residuales oscila alrededor de cero sin tendencias ni cambios de varianza visibles. El correlograma ACF muestra que prácticamente todos los rezagos caen dentro de las bandas de confianza, indicando que el modelo absorbió la dependencia temporal de la serie con un único parámetro. El histograma se aproxima bien a una distribución normal centrada en cero, lo que valida la construcción de intervalos de confianza fiables para el pronóstico.
res <- residuals(modelo2)
lb <- Box.test(res, lag = 10, fitdf = 1, type = "Ljung-Box")
sw <- shapiro.test(as.numeric(res))
tabla_diag <- data.frame(
Prueba = c("Ljung-Box (10 rezagos)", "Shapiro-Wilk"),
Estadístico = round(c(lb$statistic, sw$statistic), 4),
p_valor = c(sprintf("%.4f", lb$p.value),
sprintf("%.4f", sw$p.value)),
Conclusión = c(ifelse(lb$p.value > 0.05,
"Ruido blanco — no se rechaza H₀",
"Autocorrelación remanente"),
ifelse(sw$p.value > 0.05,
"Normalidad — no se rechaza H₀",
"No normalidad")),
check.names = FALSE
)
kable(tabla_diag,
col.names = c("Prueba", "Estadístico", "p-valor", "Conclusión"),
caption = "Tabla 8. Pruebas diagnósticas sobre los residuales del ARIMA(0,1,1)")| Prueba | Estadístico | p-valor | Conclusión | |
|---|---|---|---|---|
| X-squared | Ljung-Box (10 rezagos) | 7.7981 | 0.5546 | Ruido blanco — no se rechaza H₀ |
| W | Shapiro-Wilk | 0.9967 | 0.8173 | Normalidad — no se rechaza H₀ |
tabla_res <- data.frame(
Estadístico = c("Media", "Desv. estándar", "Asimetría", "Curtosis"),
Valor = c(round(mean(as.numeric(res), na.rm = TRUE), 4),
round(sd(as.numeric(res), na.rm = TRUE), 4),
round(skewness(as.numeric(res), na.rm = TRUE), 4),
round(kurtosis(as.numeric(res), na.rm = TRUE) + 3, 4)),
Referencia = c("≈ 0", "—", "≈ 0", "≈ 3")
)
kable(tabla_res,
caption = "Tabla 9. Estadísticos descriptivos de los residuales",
align = c("l", "r", "c"))| Estadístico | Valor | Referencia |
|---|---|---|
| Media | 0.0018 | ≈ 0 |
| Desv. estándar | 5.5719 | — |
| Asimetría | 0.0115 | ≈ 0 |
| Curtosis | 3.4065 | ≈ 3 |
El diagnóstico del ARIMA(0,1,1) es satisfactorio en todos los frentes. La prueba de Ljung-Box arroja un p-valor de 0.5546 (> 0.05): no hay autocorrelación remanente en los residuales, confirmando que el modelo capturó toda la estructura de dependencia temporal de la serie con un solo parámetro. La prueba de Shapiro-Wilk (p = 0.8173) confirma que los residuales son aproximadamente normales, condición que valida la construcción de intervalos de confianza fiables para el pronóstico. La media de los residuales es prácticamente nula, la asimetría es cercana a cero y la curtosis próxima a 3, valores coherentes con un proceso de ruido blanco gaussiano. En contraste, la caminata aleatoria presenta un Ljung-Box con p < 0.001, evidenciando una autocorrelación residual severamente significativa que la invalida como modelo predictivo.
4.9 Desempeño sobre la ventana de prueba
fc_test <- forecast(modelo2, h = length(ventana2), level = c(80, 95))
acc_train <- accuracy(modelo2)
acc_test <- accuracy(as.numeric(fc_test$mean), as.numeric(ventana2))
acc_df <- data.frame(
Conjunto = c("Entrenamiento", "Prueba"),
RMSE = c(round(acc_train[,"RMSE"], 3), round(acc_test[,"RMSE"], 3)),
MAE = c(round(acc_train[,"MAE"], 3), round(acc_test[,"MAE"], 3)),
MAPE = c(round(acc_train[,"MAPE"], 3), round(acc_test[,"MAPE"], 3))
)
kable(acc_df,
caption = "Tabla 10. Métricas de exactitud — entrenamiento vs. prueba",
align = c("l", "r", "r", "r"))| Conjunto | RMSE | MAE | MAPE |
|---|---|---|---|
| Entrenamiento | 5.562 | 4.352 | 16.090 |
| Prueba | 8.474 | 6.888 | 41.587 |
Las métricas de la ventana de prueba muestran una degradación respecto al entrenamiento, especialmente en el MAPE. Esto es atribuible a que los últimos 15 días de diciembre de 2018 presentaron mayor variabilidad que el promedio del año de entrenamiento, lo que dificulta el pronóstico de un modelo univariado que no incorpora variables meteorológicas como covariables. Sin embargo, el modelo sigue siendo útil como referencia de corto plazo y supera claramente a la caminata aleatoria.
5 Pronóstico
Una vez validado el modelo, se generan pronósticos para los 15 días posteriores al período observado, correspondientes al final de diciembre de 2018:
h_pronostico <- 15
pronostico2 <- modelo2 %>%
forecast(h = h_pronostico, level = c(80, 95))
fechas_pron <- max(datos_diarios$fecha) + 1:h_pronostico
df_pron <- data.frame(
Fecha = format(fechas_pron, "%d %b %Y"),
Pronóstico = round(as.numeric(pronostico2$mean), 2),
Inf_80 = round(as.numeric(pronostico2$lower[,1]), 2),
Sup_80 = round(as.numeric(pronostico2$upper[,1]), 2),
Inf_95 = round(as.numeric(pronostico2$lower[,2]), 2),
Sup_95 = round(as.numeric(pronostico2$upper[,2]), 2)
)
kable(df_pron,
col.names = c("Fecha", "Pronóstico (µg/m³)",
"Inf. 80%", "Sup. 80%", "Inf. 95%", "Sup. 95%"),
caption = "Tabla 11. Pronóstico de O₃ a 15 días con bandas de confianza")| Fecha | Pronóstico (µg/m³) | Inf. 80% | Sup. 80% | Inf. 95% | Sup. 95% |
|---|---|---|---|---|---|
| 02 ene. 2019 | 26.55 | 19.40 | 33.70 | 15.61 | 37.49 |
| 03 ene. 2019 | 26.55 | 19.09 | 34.01 | 15.15 | 37.96 |
| 04 ene. 2019 | 26.55 | 18.80 | 34.30 | 14.70 | 38.41 |
| 05 ene. 2019 | 26.55 | 18.52 | 34.59 | 14.26 | 38.84 |
| 06 ene. 2019 | 26.55 | 18.24 | 34.86 | 13.85 | 39.26 |
| 07 ene. 2019 | 26.55 | 17.98 | 35.12 | 13.44 | 39.66 |
| 08 ene. 2019 | 26.55 | 17.72 | 35.38 | 13.05 | 40.05 |
| 09 ene. 2019 | 26.55 | 17.48 | 35.63 | 12.67 | 40.43 |
| 10 ene. 2019 | 26.55 | 17.23 | 35.87 | 12.30 | 40.80 |
| 11 ene. 2019 | 26.55 | 17.00 | 36.11 | 11.94 | 41.16 |
| 12 ene. 2019 | 26.55 | 16.77 | 36.34 | 11.59 | 41.52 |
| 13 ene. 2019 | 26.55 | 16.54 | 36.56 | 11.24 | 41.86 |
| 14 ene. 2019 | 26.55 | 16.32 | 36.78 | 10.90 | 42.20 |
| 15 ene. 2019 | 26.55 | 16.11 | 37.00 | 10.58 | 42.53 |
| 16 ene. 2019 | 26.55 | 15.89 | 37.21 | 10.25 | 42.85 |
obs_plot <- data.frame(
fecha = datos_diarios$fecha,
O3 = datos_diarios$O3_mean
) %>% tail(80)
ultimo_obs <- data.frame(
fecha = max(datos_diarios$fecha),
O3 = as.numeric(tail(datos_diarios$O3_mean, 1))
)
fc_plot <- data.frame(
fecha = c(ultimo_obs$fecha, fechas_pron),
O3 = c(ultimo_obs$O3, as.numeric(pronostico2$mean)),
inf80 = c(ultimo_obs$O3, as.numeric(pronostico2$lower[,1])),
sup80 = c(ultimo_obs$O3, as.numeric(pronostico2$upper[,1])),
inf95 = c(ultimo_obs$O3, as.numeric(pronostico2$lower[,2])),
sup95 = c(ultimo_obs$O3, as.numeric(pronostico2$upper[,2]))
)
p_fc <- ggplot() +
geom_ribbon(data = fc_plot,
aes(x = fecha, ymin = inf95, ymax = sup95),
fill = "#AED6F1", alpha = 0.4) +
geom_ribbon(data = fc_plot,
aes(x = fecha, ymin = inf80, ymax = sup80),
fill = "#5DADE2", alpha = 0.45) +
geom_line(data = obs_plot,
aes(x = fecha, y = O3, color = "Observado"),
linewidth = 0.6) +
geom_line(data = fc_plot,
aes(x = fecha, y = O3, color = "Pronóstico"),
linewidth = 1.2) +
geom_vline(xintercept = as.numeric(max(datos_diarios$fecha)),
linetype = "dashed", color = "grey50") +
scale_color_manual(values = c("Observado" = "#1F3D6B",
"Pronóstico" = "#E74C3C")) +
scale_x_date(date_breaks = "2 weeks", date_labels = "%d %b") +
labs(title = "Pronóstico de O₃ a 15 días — ARIMA(0,1,1)",
subtitle = "Últimos 80 días observados + pronóstico con bandas al 80% y 95%",
x = NULL, y = "O₃ (µg/m³)", color = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "top",
plot.title = element_text(face = "bold"))
ggplotly(p_fc) %>% layout(hovermode = "x unified", margin = list(t = 80))El modelo proyecta concentraciones de O₃ que convergen rápidamente hacia un nivel de equilibrio cercano a los 26.6 µg/m³, valor coherente con los niveles observados al final del período de entrenamiento. Las bandas de confianza se amplían progresivamente con el horizonte de pronóstico, reflejando la acumulación natural de incertidumbre: el intervalo al 95% para el día 1 es de aproximadamente ±11 µg/m³, mientras que para el día 15 se extiende a ±22 µg/m³. Todos los valores pronosticados se mantienen por debajo del umbral de alerta de la OMS (100 µg/m³), lo que sugiere condiciones de calidad del aire aceptables para el período proyectado.
6 Conclusiones
La concentración de O₃ en Cali durante 2018 presentó una dinámica estacional marcada por el régimen climático bimodal del Valle del Cauca. Los meses de julio a septiembre, correspondientes al periodo seco con mayor insolación, concentraron los niveles más altos. Septiembre fue el mes más contaminante del año, mientras que los meses lluviosos como mayo registraron los valores mínimos. Todos los promedios diarios se mantuvieron por debajo del umbral de 100 µg/m³ establecido por la Resolución 2254 de 2017 (MADS, 2017), aunque a nivel horario los picos del mediodía pueden superar ese límite durante los episodios más intensos.
La serie diaria de O₃ no es estacionaria en niveles, condición confirmada de forma consistente tanto por la prueba ADF (p = 0.3672) como por la prueba KPSS, que rechaza la hipótesis nula de estacionariedad. Una primera diferenciación resultó suficiente para estabilizar la media (ADF sobre ΔO₃: p = 0.01), determinando d = 1 como el orden de integración del modelo.
El análisis de la ACF y PACF de la serie diferenciada orientó hacia
un modelo de media móvil puro, y la selección mediante
auto.arima() con búsqueda exhaustiva confirmó el
ARIMA(0,1,1) como el de menor AICc. Su único coeficiente (ma1 = -0.739)
es estadísticamente significativo, indicando que los errores de
pronóstico del día anterior tienen una influencia real y correctiva
sobre la estimación del día actual. Este mecanismo de reversión es
coherente con la naturaleza fotoquímica del contaminante, cuya
producción depende de condiciones meteorológicas que no persisten
indefinidamente.
El diagnóstico de residuales valida el modelo en todos sus supuestos: la prueba de Ljung-Box (p = 0.5546) confirma ausencia de autocorrelación remanente y la prueba de Shapiro-Wilk (p = 0.8173) confirma normalidad, habilitando la construcción de intervalos de confianza fiables. En contraste, la caminata aleatoria presenta autocorrelación residual significativa que la invalida como modelo predictivo.
El pronóstico a 15 días proyecta concentraciones estables de O₃ en torno a 26.6 µg/m³, con bandas de incertidumbre que se amplían razonablemente con el horizonte temporal. Es importante reconocer, sin embargo, que el modelo trabaja exclusivamente con la información histórica de la propia serie, sin incorporar variables meteorológicas como temperatura, radiación solar, velocidad del viento o humedad relativa, las cuales tienen una influencia directa en la formación y dispersión del O₃. Esta condición constituye la principal limitación del modelo: aunque su desempeño es satisfactorio para pronósticos de corto plazo, su capacidad predictiva podría verse reducida ante cambios abruptos en las condiciones atmosféricas que no queden reflejados en el comportamiento reciente de la serie. Como recomendación para trabajos futuros, una extensión natural sería incorporar estas variables meteorológicas como covariables mediante un modelo ARIMAX, lo que permitiría capturar de forma más explícita los mecanismos fotoquímicos que determinan la concentración de ozono y potencialmente mejorar tanto el ajuste como el horizonte útil de pronóstico.
En respuesta a la pregunta de investigación, los resultados confirman que es posible modelar y pronosticar la concentración diaria de O₃ en Cali a partir de su comportamiento histórico. El modelo ARIMA(0,1,1) demostró ser suficiente para capturar la estructura de dependencia temporal de la serie, con residuales que se comportan como ruido blanco y pronósticos con intervalos de confianza coherentes con la variabilidad observada.
7 Bibliografía
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. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. https://otexts.com/fpp3/
Hyndman, R.J. & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1–22.
IDEAM (2010). Protocolo para el monitoreo y seguimiento de la calidad del aire. Instituto de Hidrología, Meteorología y Estudios Ambientales. Bogotá.
Kumar, U. & Jain, V.K. (2010). ARIMA forecasting of ambient air pollutants (O₃, NO, NO₂ and CO). Stochastic Environmental Research and Risk Assessment, 24(5), 751–760.
Ministerio de Ambiente y Desarrollo Sostenible (MADS). (2017). Resolución 2254 de 2017, por la cual se adopta la norma de calidad del aire o nivel de inmisión. Bogotá: MADS.
Organización Mundial de la Salud (OMS). (2021). Directrices mundiales de la OMS sobre la calidad del aire. OMS. https://www.who.int/es/news/item/22-09-2021-who-s-new-air-quality-guidelines
R Core Team (2024). R: A Language and Environment for Statistical Computing. https://www.R-project.org/