library(readxl)
library(xts)
library(zoo)
library(dplyr)
library(lubridate)
library(fpp2)
library(tseries)
library(urca)
library(ggplot2)
library(plotly)
library(DT)
library(kableExtra)
library(urca)El ozono troposférico (O₃) es uno de los contaminantes atmosféricos criterio más importantes en entornos urbanos. A diferencia del ozono estratosférico que protege a la Tierra de la radiación ultravioleta, el O₃ a nivel del suelo es un contaminante secundario: no se emite directamente, sino que se forma mediante reacciones fotoquímicas entre el dióxido de nitrógeno (NO₂) y los compuestos orgánicos volátiles (COV) en presencia de radiación solar. Sus principales fuentes precursoras son los gases de escape vehicular, la quema de combustibles fósiles y los procesos industriales.
Desde el punto de vista normativo, Colombia regula este contaminante mediante la Resolución 2254 de 2017 del Ministerio de Ambiente y Desarrollo Sostenible (MADS), vigente a partir del 1 de enero de 2018, que establece un nivel máximo permisible de 100 µg/m³ para una exposición de 8 horas. Superar este umbral de forma sostenida representa un riesgo significativo para la salud pública, especialmente para la población con enfermedades respiratorias o cardiovasculares previas.
La ciudad de Santiago de Cali ubicada a 995 m.s.n.m., enclavada en el corazón del Valle del Cauca y rodeada por las cordilleras Central y Occidental, presenta condiciones geográficas y climáticas únicas que la convierten en un escenario particularmente sensible a la acumulación de contaminantes fotoquímicos como el ozono troposférico. Su régimen climático bimodal, su alta radiación solar durante todo el año y la limitada dispersión atmosférica provocada por las barreras montañosas que enmarcan el valle, crean condiciones propicias para que el O₃ se acumule en determinadas épocas del año, representando un riesgo latente para la salud de sus más de 2.2 millones de habitantes.
En este contexto, el presente trabajo analiza la serie de concentración diaria promedio de O₃ (µg/m³) registrada por una estación del SVCASC durante el año 2018, con el objetivo de ajustar un modelo ARIMA (AutoRegressive Integrated Moving Average) que permita describir el comportamiento temporal del contaminante y generar pronósticos confiables, aportando así una herramienta cuantitativa de utilidad para la gestión ambiental urbana.
El desarrollo de este trabajo partió de una pregunta concreta: ¿es posible anticipar el comportamiento del ozono troposférico en Cali a partir de su propio historial? Para responderla, se siguió un proceso de análisis por etapas que va desde la preparación de los datos hasta la generación y evaluación de pronósticos, todo desarrollado en el entorno estadístico R.
La base de datos con la que se trabajó contiene registros horarios de diversas variables ambientales y meteorológicas recolectadas en la ciudad de Cali a lo largo del año 2018: ozono (O₃), temperatura, humedad relativa, radiación solar, velocidad y dirección del viento, y precipitación. De todas ellas, se seleccionó el ozono troposférico (O₃ en µg/m³) como variable de análisis, dada su relevancia para la salud pública y su regulación directa por la normativa ambiental colombiana.
Dado que trabajar con datos horarios generaría una serie de más de 8.700 observaciones, los registros se agregaron calculando el promedio diario, obteniendo una serie que resume de forma más limpia el comportamiento del contaminante a lo largo del año. Los valores faltantes presentes en la serie fueron tratados mediante interpolación lineal, técnica que estima los datos ausentes trazando una línea recta entre el último valor conocido antes del hueco y el primero disponible después de él. Sin embargo, cuando el número de días consecutivos sin datos es muy elevado, esta técnica pierde confiabilidad al no poder capturar la variabilidad real del contaminante en ese período. Por esta razón, se definió cuidadosamente la ventana de entrenamiento excluyendo los períodos con ausencias masivas de datos, garantizando que el modelo se ajuste únicamente sobre observaciones reales y representativas.
La serie se dividió en dos partes: una ventana de entrenamiento (marzo - noviembre 2018) para construir el modelo, y una ventana de prueba (diciembre 2018 - enero 2019) para evaluar su capacidad predictiva comparando los pronósticos con los valores reales observados.
Antes de ajustar cualquier modelo ARIMA, es necesario verificar que la serie sea estacionaria es decir, que su comportamiento promedio y su variabilidad no cambien de forma sistemática a lo largo del tiempo. Esta condición es fundamental porque los modelos ARIMA están diseñados para series cuya estructura es estable en el tiempo.
Para verificar esta condición se recurrió a tres herramientas de forma simultánea: la inspección visual de la serie y su función de autocorrelación, la prueba ADF (Augmented Dickey-Fuller) que evalúa si la serie tiene raíz unitaria, y la prueba KPSS (Kwiatkowski-Phillips-Schmidt-Shin) que parte del supuesto contrario. Usar ambas pruebas juntas es una práctica recomendada porque, al tener hipótesis opuestas, cuando coinciden en su conclusión la evidencia es mucho más sólida. En caso de que la serie no sea estacionaria, se aplicará el número de diferenciaciones necesarias para lograr esta condición.
Con la serie transformada para cumplir la condición de estacionariedad, el siguiente paso es identificar qué tipo de modelo ARIMA se ajusta mejor a los datos. Para ello se analizarán las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie diferenciada, que permiten identificar visualmente los órdenes candidatos para los componentes autorregresivo (p) y de medias móviles (q).
A partir de este análisis se propondrán varios modelos candidatos, que se compararán entre sí utilizando el criterio AICc una medida que equilibra el ajuste del modelo con su complejidad, penalizando los modelos con demasiados parámetros. Adicionalmente se realizará una búsqueda automática exhaustiva como referencia para validar la selección manual.
Una vez seleccionado el modelo, se verificará que los residuales se comporten como ruido blanco errores aleatorios sin memoria ni estructura temporal remanente. Esta verificación se realizará de forma visual, inspeccionando el gráfico de residuales y su función de autocorrelación, y de forma formal mediante la prueba de Ljung-Box, que evalúa si existe autocorrelación remanente, y la prueba de Shapiro-Wilk, que evalúa si los residuales siguen una distribución normal.
Con el modelo validado, se generarán pronósticos para la ventana de prueba con sus respectivos intervalos de confianza al 95%, que se compararán con los valores reales para evaluar la precisión del modelo. La evaluación se realizará con tres métricas complementarias: el MAE (error promedio en las mismas unidades de la serie), el RMSE (más sensible a errores grandes) y el MAPE (error expresado como porcentaje del valor real).
# Cargar datos desde Excel
Data <- read_excel("Compartir2.xlsx")
# Renombrar columnas
names(Data) <- c("fecha", "O3", "vel_viento", "dir_viento",
"temperatura", "humedad", "radiacion", "lluvia")
# Agregar a nivel diario
Data_diaria <- Data %>%
mutate(dia = as.Date(fecha)) %>%
group_by(dia) %>%
summarise(O3_diario = mean(O3, na.rm = TRUE)) %>%
ungroup()
# Crear objeto xts
Serie <- xts(Data_diaria$O3_diario, order.by = Data_diaria$dia)
# Interpolar datos faltantes
Serie <- na.approx(Serie, na.rm = FALSE)
# Tabla de verificación
data.frame(
Descripción = c("Fecha inicial", "Fecha final",
"Total observaciones", "Valores faltantes"),
Valor = c(as.character(index(Serie)[1]),
as.character(index(Serie)[nrow(Serie)]),
as.character(nrow(Serie)),
as.character(sum(is.na(Serie))))
) %>%
kable(caption = "Verificación de la serie",
align = c("l", "r")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Descripción | Valor |
|---|---|
| Fecha inicial | 2018-01-01 |
| Fecha final | 2019-01-01 |
| Total observaciones | 366 |
| Valores faltantes | 0 |
data.frame(
Estadístico = c("Mínimo", "Percentil 25", "Mediana", "Media",
"Percentil 75", "Máximo", "Desv. Estándar", "CV (%)"),
`Valor (µg/m³)` = c(
round(min(Serie, na.rm = TRUE), 2),
round(quantile(Serie, 0.25, na.rm = TRUE), 2),
round(median(as.numeric(Serie), na.rm = TRUE), 2),
round(mean(Serie, na.rm = TRUE), 2),
round(quantile(Serie, 0.75, na.rm = TRUE), 2),
round(max(Serie, na.rm = TRUE), 2),
round(sd(as.numeric(Serie), na.rm = TRUE), 2),
round((sd(as.numeric(Serie), na.rm = TRUE) /
mean(Serie, na.rm = TRUE)) * 100, 2)
)
) %>%
datatable(
caption = "Estadísticas descriptivas - O₃ diario (µg/m³)",
rownames = FALSE,
options = list(
pageLength = 10,
dom = "t",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
) %>%
formatStyle(
"Estadístico",
backgroundColor = styleEqual(
c("Media", "Máximo", "Mínimo"),
c("#d4edff", "#ffd4d4", "#d4ffd4")
),
fontWeight = styleEqual("Media", "bold")
)La concentración diaria promedio de O₃ en Cali durante 2018 presentó una media de 29.52 µg/m³ y una mediana de 29.36 µg/m³, valores muy cercanos entre sí, lo que sugiere una distribución aproximadamente simétrica sin sesgos extremos pronunciados. La desviación estándar de 7.58 µg/m³ refleja una variabilidad moderada en torno a la media, explicable en gran medida por el régimen climático bimodal característico del Valle del Cauca: el año presenta dos temporadas secas la principal entre junio y agosto, y una secundaria entre diciembre y febrero y dos temporadas lluviosas la primera entre marzo y mayo, y la segunda entre septiembre y noviembre. Durante las temporadas secas, la menor nubosidad favorece una mayor radiación solar incidente, acelerando las reacciones fotoquímicas que producen O₃ y elevando sus concentraciones. Por el contrario, en las temporadas lluviosas la mayor cobertura nubosa reduce la radiación disponible, inhibiendo parcialmente la formación del contaminante y diluyéndolo mediante el lavado atmosférico de sus precursores.
El valor mínimo registrado fue de 5.50 µg/m³, correspondiente a días con alta nubosidad o precipitación intensa que inhiben la formación fotoquímica del ozono. Por su parte, el valor máximo fue de 58.68 µg/m³, el cual, si bien es elevado, se mantiene por debajo del límite máximo permisible de 100 µg/m³ para exposición de 8 horas establecido por la Resolución 2254 de 2017, lo que indica que durante 2018 la concentración de O₃ en la estación monitoreada no superó los umbrales de alerta normativa. El 50% central de los datos se concentra entre 24.52 y 34.75 µg/m³,rango que refleja las condiciones típicas de ozono en la ciudad bajo condiciones meteorológicas normales.
El coeficiente de variación de 25.68% refleja una variabilidad moderada-alta de la serie respecto a su media, lo cual contrasta con variables climáticas como la temperatura que en Cali presentan CV inferiores al 5%. Esta mayor dispersión relativa es característica del ozono troposférico, cuya formación depende de la combinación simultánea de múltiples factores fotoquímicos y meteorológicos que varían considerablemente día a día.
serie_df <- data.frame(
fecha = as.Date(index(Serie)),
O3 = as.numeric(Serie)
)
p <- ggplot(serie_df, aes(x = fecha, y = O3)) +
geom_line(color = "black", linewidth = 0.5) +
geom_smooth(method = "loess", color = "#e74c3c", se = FALSE,
linewidth = 0.8, linetype = "dashed") +
scale_x_date(date_breaks = "1 month", date_labels = "%B") +
ggtitle("Concentración diaria promedio de O₃ en Cali - 2018") +
xlab("Mes") +
ylab("O₃ (µg/m³)") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)La serie de concentración diaria de O₃ en Cali durante 2018 revela un comportamiento dinámico que la curva de tendencia suavizada (LOESS), representada en rojo punteado, permite interpretar con mayor claridad al eliminar el ruido de las fluctuaciones diarias y mostrar el patrón estructural subyacente.
Durante el período enero - mayo se observan valores moderados y relativamente estables, oscilando en torno a la media de 29.52 µg/m³, con fluctuaciones frecuentes pero sin una tendencia definida. Este comportamiento es coherente con la alternancia entre la temporada seca de inicio de año y la primera temporada lluviosa del régimen bimodal del Valle del Cauca, que limitan la acumulación sostenida del contaminante. Cabe señalar que el segmento de línea con pendiente positiva uniforme visible en febrero no representa un comportamiento real de la serie, sino el resultado de la interpolación lineal para estimar los 20 días sin registro en ese período, atribuibles a interrupciones técnicas en el equipo sensor de la estación de monitoreo. Esta técnica estima los valores faltantes trazando una línea recta entre el último dato conocido antes del hueco y el primero disponible después de él, garantizando la continuidad temporal necesaria para el modelamiento ARIMA.
Entre junio y septiembre la curva LOESS muestra un incremento gradual y sostenido, alcanzando los valores más altos del año. El máximo absoluto de 58.68 µg/m³ se registró el 27 de agosto, precisamente en la temporada seca principal del Valle del Cauca (junio - agosto), cuando la menor nubosidad y la mayor disponibilidad de radiación solar aceleran las reacciones fotoquímicas entre el NO₂ vehicular y los compuestos orgánicos volátiles (COV), favoreciendo la formación de ozono troposférico. Este patrón es consistente con lo reportado por el DAGMA en sus boletines de calidad del aire, donde los episodios de mayor concentración de O₃ en Cali se asocian históricamente con los meses de menor precipitación y mayor irradiación solar.
A partir de octubre y hasta diciembre la curva LOESS evidencia un descenso sostenido de las concentraciones, asociado al inicio de la segunda temporada lluviosa, que reduce la radiación solar disponible y favorece el lavado atmosférico de los precursores del ozono. Este descenso es especialmente marcado en noviembre y diciembre, donde los valores se acercan nuevamente a los niveles observados a inicio de año.
En cuanto a la variabilidad, la serie presenta oscilaciones frecuentes y pronunciadas entre días consecutivos a lo largo de todo el año, con diferencias de hasta 20 µg/m³ entre días adyacentes, lo que evidencia la alta sensibilidad del O₃ a las condiciones meteorológicas locales de corto plazo especialmente la radiación solar, la velocidad del viento y la precipitación que un modelo univariado como el ARIMA no captura directamente como variables explicativas.
serie_df_mes <- data.frame(
fecha = as.Date(index(Serie)),
O3 = as.numeric(Serie)
) %>%
mutate(mes = month(fecha, label = TRUE, abbr = FALSE))
p_box <- ggplot(serie_df_mes, aes(x = mes, y = O3, fill = mes)) +
geom_boxplot(show.legend = FALSE) +
scale_fill_manual(values = colorRampPalette(
c("#2c7bb6", "#abd9e9", "#fdae61", "#d7191c"))(12)) +
labs(
title = "Distribución mensual de O₃ - Cali 2018",
x = "Mes",
y = "O₃ (µg/m³)"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(angle = 45, hjust = 1)
)
ggplotly(p_box)El boxplot mensual permite profundizar en el comportamiento del O₃ a lo largo del año, revelando patrones que la serie temporal diaria no muestra con tanta claridad.
En términos de concentración central, se observa una progresión clara a lo largo del año: los meses de marzo y abril presentan las medianas más altas del primer semestre, superando los 30 µg/m³, mientras que mayo registra la mediana más baja de todo el año, cercana a los 24 µg/m³, coincidiendo con el pico de la primera temporada lluviosa. A partir de julio las medianas escalan progresivamente, alcanzando su punto máximo en agosto y septiembre los meses más cálidos y secos del año con medianas que superan los 35 µg/m³, coherente con las condiciones fotoquímicas más favorables para la formación de O₃.
En cuanto a la variabilidad interna, los meses de septiembre y octubre presentan las cajas más amplias de toda la serie, lo que refleja una alta dispersión en las concentraciones diarias durante este período. Esto es consistente con la transición entre la temporada seca y la segunda temporada lluviosa, donde la alternancia entre días despejados y días con precipitación genera diferencias notables entre días consecutivos. Por el contrario, noviembre presenta la caja más compacta del segundo semestre, sugiriendo un comportamiento más homogéneo asociado a condiciones lluviosas más estables.
Respecto a los valores atípicos, el punto más extremo corresponde al mes de agosto con el valor máximo de 58.68 µg/m³ ya identificado anteriormente. También se observan valores atípicos en abril, mayo y junio, asociados a episodios puntuales de alta radiación solar en medio de períodos predominantemente lluviosos. En el extremo inferior, enero presenta un valor atípico bajo que corresponde al dato interpolado del inicio de la serie.
En conjunto, el boxplot confirma visualmente el patrón estacional del O₃ en Cali durante 2018 y refuerza la necesidad de un modelo que capture las dependencias temporales entre observaciones consecutivas, como el modelo ARIMA seleccionado.
# Identificar valores atípicos
datos_df <- data.frame(fecha = index(Serie), O3 = as.numeric(Serie))
datos_df %>% filter(O3 > 45 | O3 < 10) %>% arrange(desc(O3))Respecto a los valores atípicos, la serie presenta varios días con concentraciones notablemente superiores a la media. Los picos más pronunciados se registraron el 27 de agosto de 2018 (58.68 µg/m³), el 29 de septiembre (50.06 µg/m³) y el 20 de septiembre (47.95 µg/m³). Estos episodios coinciden con el período de transición entre la temporada seca principal (junio - agosto) y el inicio de la segunda temporada lluviosa, cuando la combinación de alta radiación solar acumulada, baja velocidad del viento y mayor concentración de precursores fotoquímicos como el NO₂ proveniente del tráfico vehicular, crea condiciones ideales para la acumulación del ozono troposférico. Adicionalmente, el pico del 19 de abril (46.22 µg/m³) corresponde al final de la primera temporada lluviosa, donde episodios de cielo despejado intercalados con días lluviosos pueden generar picos puntuales de O₃ por la alta irradiación solar disponible.
En el extremo opuesto, el valor mínimo de 5.50 µg/m³ registrado el 1 de enero de 2019 corresponde al último punto de la serie, en plena temporada seca de inicio de año, lo cual podría estar asociado a condiciones locales de alta humedad o precipitación puntual que inhibieron temporalmente la formación del contaminante.
# Ventana de entrenamiento: marzo - noviembre 2018 (sin datos interpolados)
ventana <- window(Serie, start = "2018-03-01", end = "2018-11-30")
# Ventana de prueba: diciembre 2018
ventana2 <- window(Serie, start = "2018-12-01")
# Tabla resumen de ventanas
data.frame(
Ventana = c("Entrenamiento", "Prueba"),
`Fecha inicio` = c("2018-03-01", "2018-12-01"),
`Fecha fin` = c("2018-11-30", "2019-01-01"),
Observaciones = c(length(ventana), length(ventana2))
) %>%
kable(caption = "División de la serie en ventanas de entrenamiento y prueba",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Ventana | Fecha.inicio | Fecha.fin | Observaciones |
|---|---|---|---|
| Entrenamiento | 2018-03-01 | 2018-11-30 | 275 |
| Prueba | 2018-12-01 | 2019-01-01 | 32 |
La serie se divide en dos ventanas: una de entrenamiento con 275 observaciones (marzo - noviembre 2018), que excluye deliberadamente los meses de enero y febrero debido a la presencia de 20 días consecutivos completamente sin registro en febrero, cuya interpolación lineal no es confiable para períodos tan extensos. De esta forma se garantiza que el modelo se ajuste únicamente sobre datos reales y representativos del comportamiento del O₃ en Cali. La ventana de prueba cuenta con 32 observaciones (diciembre 2018 - enero 2019), sobre las cuales se evaluará la capacidad predictiva del modelo comparando los valores pronosticados con los valores reales observados.
p2 <- autoplot(ventana) +
ggtitle("Ventana de entrenamiento - O₃ diario (mar - nov 2018)") +
xlab("Fecha") +
ylab("O₃ (µg/m³)") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p2)La ventana de entrenamiento comprende el período marzo - noviembre de 2018, excluyendo los meses de enero y febrero para garantizar que el modelo se ajuste únicamente sobre datos reales y confiables. La serie muestra claramente las tres fases identificadas en el análisis descriptivo: valores moderados con alta variabilidad durante marzo y abril, un descenso hacia los valores mínimos del período en mayo coincidiendo con el pico de la primera temporada lluviosa, un incremento gradual y sostenido entre junio y septiembre alcanzando los valores máximos del año, y un descenso progresivo a partir de octubre asociado al inicio de la segunda temporada lluviosa.
La media no es constante a lo largo del período oscila claramente entre fases de mayor y menor concentración y la variabilidad de las fluctuaciones diarias parece aumentar en el período agosto - septiembre respecto al resto de la ventana. Estas características constituyen evidencia visual de no estacionariedad, lo que justifica la aplicación de las pruebas formales a continuación.
p3 <- ggAcf(ventana) +
ggtitle("ACF - Serie O₃ original") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p3)La ACF de la serie original presenta autocorrelaciones significativas y decrecientes de forma lenta hasta rezagos muy lejanos, lo cual es un patrón clásico de una serie no estacionaria. En una serie estacionaria, la ACF caería rápidamente a cero después de los primeros rezagos. Este comportamiento confirma la inspección visual y justifica la aplicación de las pruebas formales de raíz unitaria.
# Prueba ADF - H0: la serie NO es estacionaria (tiene raíz unitaria)
adf <- adf.test(ventana)
# Prueba KPSS - H0: la serie ES estacionaria
kpss <- ur.kpss(ventana)
# Tabla resumen
data.frame(
Prueba = c("ADF (Augmented Dickey-Fuller)", "KPSS"),
Hipótesis_nula = c("La serie NO es estacionaria", "La serie ES estacionaria"),
Estadístico = c(round(adf$statistic, 4), round(kpss@teststat, 4)),
`P-valor` = c(round(adf$p.value, 4), "Ver valor crítico")
) %>%
kable(caption = "Pruebas formales de estacionariedad",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Prueba | Hipótesis_nula | Estadístico | P.valor | |
|---|---|---|---|---|
| Dickey-Fuller | ADF (Augmented Dickey-Fuller) | La serie NO es estacionaria | -2.1635 | 0.5071 |
| KPSS | La serie ES estacionaria | 0.5650 | Ver valor crítico |
Los resultados de las pruebas formales confirman lo observado visualmente:
Prueba ADF (Augmented Dickey-Fuller):
El estadístico obtenido fue -2.1635 con un p-valor de 0.5071. Para entender este resultado es importante recordar que la prueba ADF funciona bajo la lógica de que mientras más negativo sea el estadístico, mayor es la evidencia en contra de la raíz unitaria. Los umbrales críticos habituales son -2.57 al 10%, -2.86 al 5% y -3.43 al 1%, por lo que el estadístico obtenido no supera ninguno de estos umbrales. Adicionalmente, el p-valor de 0.5071 es superior al nivel de significancia del 5%, lo que significa que no existe evidencia estadística suficiente para rechazar la hipótesis nula. En términos prácticos, la serie presenta raíz unitaria: su media no es constante en el tiempo y los valores pasados tienen una influencia persistente sobre los valores futuros, características que impiden ajustar directamente un modelo ARIMA sin antes transformar la serie.
Prueba KPSS (Kwiatkowski-Phillips-Schmidt-Shin):
El estadístico obtenido fue 0.5650, valor que supera el umbral crítico al 5% (0.463), indicando que a ese nivel de significancia se rechaza la hipótesis nula de estacionariedad. Es importante señalar que el estadístico se ubica muy cerca del valor crítico al 2.5% (0.574), lo que refuerza la solidez de esta conclusión no se trata de un resultado marginal sino de una evidencia clara de no estacionariedad.
En este caso ambas pruebas apuntan en la misma dirección: la serie de O₃ diario no es estacionaria y requiere al menos una diferenciación antes de ajustar el modelo ARIMA.
nd <- ndiffs(ventana)
data.frame(
Resultado = "Número de diferencias requeridas (ndiffs)",
Valor = nd
) %>%
kable(align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Resultado | Valor |
|---|---|
| Número de diferencias requeridas (ndiffs) | 1 |
Mediante la función ndiffs() confirma que se requiere
una diferenciación (d = 1) para lograr la
estacionariedad de la serie, resultado consistente con las pruebas ADF y
KPSS realizadas anteriormente.
# Aplicar primera diferencia
miserie <- diff(ventana) %>% na.omit()
p4 <- autoplot(miserie) +
ggtitle("Serie O₃ diferenciada (d = 1)") +
xlab("Fecha") +
ylab("Δ O₃ (µg/m³)") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p4)Dado que las pruebas anteriores de ADF y KPSS indicaron que la serie original de \(O_3\) no era estacionaria, se determinó que se requería una primera diferenciacion (\(d = 1\)). Al aplicar esta transformación, se observa gráficamente que la serie corregida fluctúa alrededor de una media constante de cero y muestra una varianza estable a lo largo del tiempo, logrando así la estacionariedad necesaria para la fase de modelado.
adf2 <- adf.test(miserie)
data.frame(
Prueba = "ADF - Serie diferenciada",
Hipótesis_nula = "La serie NO es estacionaria",
Estadístico = round(adf2$statistic, 4),
`P-valor` = round(adf2$p.value, 4)
) %>%
kable(caption = "Verificación de estacionariedad tras diferenciar",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Prueba | Hipótesis_nula | Estadístico | P.valor | |
|---|---|---|---|---|
| Dickey-Fuller | ADF - Serie diferenciada | La serie NO es estacionaria | -9.1857 | 0.01 |
El resultado de la prueba ADF sobre la serie diferenciada arroja un estadístico de -9.186 con un p-valor de 0.01, el cual es menor al nivel de significancia del 5%, por lo que se rechaza la hipótesis nula de no estacionariedad, confirmando que la serie diferenciada sí es estacionaria.
Un estadístico de -9.186 es un valor considerablemente negativo que supera con amplio margen los umbrales críticos habituales de la prueba (-2.57 al 10%, -2.86 al 5% y -3.43 al 1%), lo que indica que la estacionariedad lograda tras la diferenciación es sólida y no un resultado marginal. El p-valor reportado de 0.01 corresponde al límite mínimo, sugiriendo que el p-valor real podría ser incluso menor, lo que refuerza aún más el rechazo de la hipótesis nula.
Este resultado es completamente consistente con lo observado en el gráfico de la serie diferenciada, donde la media se mantiene estable alrededor de cero y no se aprecian cambios sistemáticos en la varianza. Con esto queda definitivamente confirmado que d = 1 es suficiente para lograr la estacionariedad, y se procede a identificar los órdenes p y q mediante el análisis de la ACF y PACF de la serie diferenciada.
p5 <- ggAcf(miserie) +
ggtitle("ACF - Serie O₃ diferenciada") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
p6 <- ggPacf(miserie) +
ggtitle("PACF - Serie O₃ diferenciada") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
subplot(ggplotly(p5), ggplotly(p6), nrows = 1, shareY = FALSE)El análisis conjunto de la ACF y la PACF de la serie diferenciada permite identificar la estructura de dependencia temporal que guiará la selección del modelo ARIMA.
En la ACF, lo primero que llama la atención es el spike negativo pronunciado en el rezago 1 (aproximadamente -0.45), que cae de forma abrupta dentro de las bandas de confianza en los rezagos siguientes. Este comportamiento un corte brusco tras el primer rezago es la firma típica de un proceso de medias móviles de orden 1, MA(1). En otras palabras, los cambios diarios en la concentración de O₃ parecen estar influenciados principalmente por el error del día inmediatamente anterior, sin que esa influencia se prolongue más allá.
En la PACF, el panorama es algo diferente. Se observan spikes significativos y negativos en los primeros 3 rezagos, con una caída hacia las bandas de confianza a partir del rezago 4. Este patrón es indicativo de un componente autorregresivo de orden 3, donde el valor actual de la serie depende de sus tres valores pasados más recientes.
Combinando ambas lecturas, los modelos candidatos más razonables son:
modelo_auto <- auto.arima(ventana, stepwise = FALSE, approximation = FALSE)
data.frame(
Modelo = paste0("ARIMA(", paste(modelo_auto$arma[c(1,6,2)], collapse = ","), ")"),
AIC = round(modelo_auto$aic, 4),
AICc = round(modelo_auto$aicc, 4),
BIC = round(modelo_auto$bic, 4),
`Log-Likelihood` = round(modelo_auto$loglik, 4)
) %>%
kable(caption = "Resultado de auto.arima() - búsqueda exhaustiva",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Modelo | AIC | AICc | BIC | Log.Likelihood |
|---|---|---|---|---|
| ARIMA(0,1,1) | 1727.123 | 1727.167 | 1734.349 | -861.5613 |
# Modelo sugerido por auto.arima
modelo1 <- Arima(ventana, order = c(0, 1, 1))
# Modelo parsimonioso identificado visualmente
modelo2 <- Arima(ventana, order = c(1, 1, 1))
# Modelo con mayor estructura AR identificado en PACF
modelo3 <- Arima(ventana, order = c(3, 1, 1))
# Tabla comparativa
data.frame(
Modelo = c("ARIMA(0,1,1)", "ARIMA(1,1,1)", "ARIMA(3,1,1)"),
AIC = round(c(modelo1$aic, modelo2$aic, modelo3$aic), 4),
AICc = round(c(modelo1$aicc, modelo2$aicc, modelo3$aicc), 4),
BIC = round(c(modelo1$bic, modelo2$bic, modelo3$bic), 4),
`Log-Likelihood` = round(c(modelo1$loglik, modelo2$loglik, modelo3$loglik), 4)
) %>%
kable(caption = "Comparación de modelos ARIMA candidatos",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center") %>%
row_spec(which.min(c(modelo1$aicc, modelo2$aicc, modelo3$aicc)),
bold = TRUE, background = "#d4edff")| Modelo | AIC | AICc | BIC | Log.Likelihood |
|---|---|---|---|---|
| ARIMA(0,1,1) | 1727.123 | 1727.167 | 1734.349 | -861.5613 |
| ARIMA(1,1,1) | 1727.442 | 1727.531 | 1738.281 | -860.7208 |
| ARIMA(3,1,1) | 1730.592 | 1730.816 | 1748.658 | -860.2962 |
La comparación de los tres modelos candidatos muestra resultados muy interesantes. El ARIMA(0,1,1) obtiene el menor AICc (1727.167), siendo el modelo más adecuado según este criterio y confirmando que la estructura de la serie, una vez eliminados los datos interpolados de febrero, es capturada de forma óptima por un modelo puramente de medias móviles de orden 1.
Vale la pena analizar las diferencias entre los modelos con algo más de detalle. El ARIMA(1,1,1) presenta un AICc de 1727.531, apenas 0.364 puntos por encima del ganador. En la práctica, diferencias de AICc menores a 2 puntos se consideran marginales, lo que significa que este modelo también sería una opción estadísticamente válida. Sin embargo, el principio de parsimonia favorece al ARIMA(0,1,1) al lograr prácticamente el mismo ajuste con un parámetro menos.
Por su parte, el ARIMA(3,1,1) presenta el AICc más alto (1730.816) a pesar de tener la mayor log-verosimilitud (-860.296), lo que indica que los parámetros adicionales del componente autorregresivo no justifican su complejidad. El AICc penaliza este exceso de parámetros, descartándolo como modelo óptimo.
Es notable que al trabajar con datos más limpios y confiables excluyendo el período con interpolaciones masivas de febrero el modelo seleccionado resulta más parsimonioso que en análisis previos, lo cual es una señal positiva de que la calidad de los datos de entrenamiento influye directamente en la simplicidad y robustez del modelo ajustado.
En conclusión, el modelo seleccionado es el ARIMA(0,1,1), que establece que el cambio diario en la concentración de O₃ depende únicamente del error del día inmediatamente anterior, sin componente autorregresivo. Este modelo será sometido al diagnóstico de residuales en la siguiente sección.
# Extraer residuales
resid <- residuals(modelo1)
resid_df <- data.frame(
fecha = index(resid),
residual = as.numeric(resid)
)
# Gráfico 1: Residuales en el tiempo
p_resid <- ggplot(resid_df, aes(x = fecha, y = residual)) +
geom_line(color = "#2c7bb6", linewidth = 0.5) +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
ggtitle("Residuales del modelo ARIMA(0,1,1)") +
xlab("Fecha") + ylab("Residual") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
# Gráfico 2: Histograma de residuales
p_hist <- ggplot(resid_df, aes(x = residual)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "#2c7bb6", color = "white", alpha = 0.8) +
geom_density(color = "red", linewidth = 1) +
ggtitle("Distribución de residuales") +
xlab("Residual") + ylab("Densidad") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
# Gráfico 3: ACF de residuales
p_acf <- ggAcf(resid) +
ggtitle("ACF de residuales") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p_resid)Gráfico de residuales en el tiempo:
El gráfico de residuales muestra que estos oscilan aleatoriamente alrededor de cero a lo largo de todo el período de entrenamiento, sin mostrar tendencias, ciclos ni patrones sistemáticos visibles. La línea roja punteada en cero sirve como referencia y confirma visualmente que los errores del modelo no presentan sesgo sistemático hacia valores positivos o negativos. Este comportamiento es precisamente lo que se espera de un modelo bien ajustado: que los errores sean impredecibles y no contengan información estructural que el modelo haya omitido.
Se observan algunos valores atípicos puntuales especialmente un pico positivo pronunciado alrededor de la observación 200, correspondiente al período de agosto - septiembre donde la serie presentaba sus concentraciones máximas del año. Estos residuales de mayor magnitud son esperables dado que los episodios extremos de O₃ son inherentemente difíciles de predecir con un modelo univariado que no incorpora variables meteorológicas como la radiación solar o la velocidad del viento como variables explicativas.
Histograma de residuales:
El histograma muestra que los residuales se distribuyen de forma aproximadamente simétrica y unimodal alrededor de cero, con una forma acampanada que se ajusta razonablemente bien a la curva de densidad normal superpuesta en rojo. La mayor concentración de residuales se ubica entre -5 y 5 µg/m³, con colas que se extienden simétricamente hacia ambos lados hasta aproximadamente ±20 µg/m³. No se observan asimetrías marcadas ni colas excesivamente pesadas, lo que anticipa un resultado favorable en la prueba formal de normalidad de Shapiro-Wilk.
ACF de residuales:
La ACF de los residuales muestra que prácticamente todas las autocorrelaciones se mantienen dentro de las bandas de confianza, con valores muy cercanos a cero en la mayoría de los rezagos. El rezago 15 presenta una barra ligeramente más alta que el resto, aunque sin superar las bandas de confianza. Este comportamiento es el esperado para residuales que se comportan como ruido blanco: no existe estructura temporal remanente que el modelo haya dejado de capturar, lo que valida el ajuste del ARIMA(0,1,1).
# Prueba Ljung-Box
ljung <- Box.test(resid, lag = 10, type = "Ljung-Box", fitdf = 3)
# Prueba Shapiro-Wilk
shapiro <- shapiro.test(as.numeric(resid))
# Tabla resumen
data.frame(
Prueba = c("Ljung-Box", "Shapiro-Wilk"),
Hipótesis_nula = c("Los residuales son ruido blanco",
"Los residuales siguen una distribución normal"),
Estadístico = c(round(ljung$statistic, 4), round(shapiro$statistic, 4)),
`P-valor` = c(round(ljung$p.value, 4), round(shapiro$p.value, 4))
) %>%
kable(caption = "Pruebas sobre los residuales del modelo ARIMA(1,1,2)",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center") %>%
row_spec(1, background = ifelse(ljung$p.value > 0.05, "#d4ffd4", "#ffd4d4")) %>%
row_spec(2, background = ifelse(shapiro$p.value > 0.05, "#d4ffd4", "#ffd4d4"))| Prueba | Hipótesis_nula | Estadístico | P.valor | |
|---|---|---|---|---|
| X-squared | Ljung-Box | Los residuales son ruido blanco | 7.2291 | 0.4054 |
| W | Shapiro-Wilk | Los residuales siguen una distribución normal | 0.9964 | 0.7918 |
Los resultados de ambas pruebas sobre los residuales del modelo ARIMA(0,1,1) son altamente satisfactorios y confirman la validez del modelo ajustado.
Prueba de Ljung-Box:
Esta prueba evalúa si los residuales presentan autocorrelación significativa en conjunto hasta el rezago 10. El estadístico obtenido fue 7.2291 con un p-valor de 0.4054, ampliamente superior al nivel de significancia del 5%. Esto significa que no se rechaza la hipótesis nula, confirmando que los residuales se comportan como ruido blanco — no existe información adicional que el modelo haya dejado de capturar en la serie. Este resultado es consistente con lo observado en la ACF de residuales, donde ningún rezago supera las bandas de confianza de forma significativa.
Prueba de Shapiro-Wilk:
Esta prueba evalúa si los residuales siguen una distribución normal. El estadístico obtenido fue 0.9964 valor muy cercano a 1, lo que indica una distribución prácticamente normal con un p-valor de 0.7918, muy superior al 5%. Por lo tanto, no se rechaza la hipótesis nula de normalidad. Este es un resultado excepcionalmente bueno: un p-valor de 0.79 indica que los residuales se ajustan muy bien a la distribución normal, mucho más de lo estrictamente necesario para validar el modelo. Este resultado es coherente con el histograma, que mostró una distribución simétrica y unimodal centrada en cero con una curva de densidad bien ajustada.
En conjunto, los tres elementos del diagnóstico, gráfico de residuales, ACF y pruebas formales confirman de forma contundente que el modelo ARIMA(0,1,1) es adecuado para describir el comportamiento de la serie de O₃ diario en Cali. Al trabajar con datos de mayor calidad excluyendo el período con interpolaciones masivas de febrero el modelo no solo resulta más parsimonioso sino que también obtiene mejores resultados en el diagnóstico de residuales, validando la decisión de ajustar la ventana de entrenamiento.
# Pronosticar 24 observaciones
pronostico <- forecast(modelo1, h = 24, level = 0.95)
# Convertir a data frame para graficar con ggplot
pronostico_df <- data.frame(
fecha = seq(as.Date("2018-12-01"), by = "day", length.out = 24),
prediccion = as.numeric(pronostico$mean),
lower = as.numeric(pronostico$lower),
upper = as.numeric(pronostico$upper)
)
# Datos históricos recientes (últimos 60 días)
historico_df <- data.frame(
fecha = index(ventana),
valor = as.numeric(ventana)
) %>% tail(60)
# Gráfico
p_forecast <- ggplot() +
geom_line(data = historico_df, aes(x = fecha, y = valor),
color = "black", linewidth = 0.6) +
geom_ribbon(data = pronostico_df, aes(x = fecha, ymin = lower, ymax = upper),
fill = "#2c7bb6", alpha = 0.2) +
geom_line(data = pronostico_df, aes(x = fecha, y = prediccion),
color = "#2c7bb6", linewidth = 0.8, linetype = "dashed") +
ggtitle("Pronóstico de O₃ diario - 24 días (ARIMA(0,1,1))") +
xlab("Fecha") + ylab("O₃ (µg/m³)") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p_forecast)El modelo ARIMA(0,1,1) genera pronósticos para los 24 días siguientes al período de entrenamiento, correspondientes a las primeras tres semanas de diciembre de 2018.
El pronóstico muestra una estabilización inmediata de las concentraciones alrededor de 21 µg/m³, valor que se encuentra por debajo de la media del período de entrenamiento. Este comportamiento es característico del modelo ARIMA(0,1,1): al no tener componente autorregresivo, el pronóstico converge rápidamente hacia un valor constante, la última observación ajustada de la serie sin oscilaciones ni tendencias hacia el futuro. Esto es coherente con el patrón histórico observado, donde los meses de noviembre y diciembre presentan una reducción progresiva de las concentraciones asociada al inicio de la segunda temporada lluviosa del Valle del Cauca.
Los intervalos de confianza al 95% se amplían progresivamente a medida que el horizonte de pronóstico se aleja, comportamiento natural e inevitable en cualquier modelo de pronóstico. Para los primeros días el intervalo es relativamente estrecho, mientras que hacia el día 24 el rango se amplía considerablemente, reflejando la incertidumbre creciente inherente a proyecciones de más largo plazo.
Desde el punto de vista ambiental, los valores pronosticados se mantienen en todo momento muy por debajo del límite normativo de 100 µg/m³ establecido por la Resolución 2254 de 2017, lo que sugiere que no se anticipan episodios de alerta por ozono en el período pronosticado bajo condiciones meteorológicas normales.
pronostico_df %>%
mutate(
fecha = as.character(fecha),
prediccion = round(prediccion, 2),
lower = round(lower, 2),
upper = round(upper, 2)
) %>%
rename(
Fecha = fecha,
`Pronóstico (µg/m³)` = prediccion,
`Límite inferior (95%)` = lower,
`Límite superior (95%)` = upper
) %>%
datatable(
caption = "Valores pronosticados de O₃ - 24 días",
rownames = FALSE,
options = list(
pageLength = 10,
dom = "tip",
columnDefs = list(list(className = "dt-center", targets = "_all"))
)
)La tabla de valores pronosticados revela una característica fundamental del modelo ARIMA(0,1,1): el pronóstico puntual es constante e igual a 21.09 µg/m³ para los 24 días del horizonte de predicción. Este comportamiento es completamente esperado y matemáticamente consistente con la estructura del modelo al no tener componente autorregresivo (p=0), el modelo no puede proyectar oscilaciones futuras basadas en valores pasados, sino que converge inmediatamente hacia el último nivel ajustado de la serie como mejor estimación para todos los períodos futuros.
Lo que sí evoluciona a lo largo del horizonte de pronóstico es la incertidumbre asociada a las predicciones. El intervalo de confianza al 95% se amplía progresivamente día a día: mientras que el 1 de diciembre el rango es de 10.07 a 32.10 µg/m³ (amplitud de 22.03 µg/m³), hacia el 24 de diciembre el intervalo se extiende desde 2.21 hasta 39.96 µg/m³ (amplitud de 37.75 µg/m³). Esta expansión refleja la incertidumbre creciente inherente a cualquier pronóstico entre más lejos en el tiempo se proyecta, mayor es el rango de valores posibles que podría tomar la serie.
Es importante destacar que en ningún momento del horizonte de pronóstico el límite superior del intervalo de confianza supera el umbral normativo de 100 µg/m³ establecido por la Resolución 2254 de 2017, lo que indica que incluso en el escenario más desfavorable contemplado por el modelo, las concentraciones de O₃ en Cali durante diciembre de 2018 no representarían un riesgo de alerta ambiental.
# Tomar los primeros 24 días de la ventana de prueba
ventana2_24 <- head(ventana2, 24)
# Data frame comparativo
comparacion_df <- data.frame(
fecha = index(ventana2_24),
real = as.numeric(ventana2_24),
prediccion = as.numeric(pronostico$mean),
lower = as.numeric(pronostico$lower),
upper = as.numeric(pronostico$upper)
)
# Gráfico comparativo
p_comp <- ggplot(comparacion_df, aes(x = fecha)) +
geom_ribbon(aes(ymin = lower, ymax = upper),
fill = "#2c7bb6", alpha = 0.2) +
geom_line(aes(y = real, color = "Real"), linewidth = 0.8) +
geom_line(aes(y = prediccion, color = "Pronóstico"),
linewidth = 0.8, linetype = "dashed") +
scale_color_manual(values = c("Real" = "black", "Pronóstico" = "#2c7bb6")) +
ggtitle("Pronóstico vs Valores reales - Diciembre 2018") +
xlab("Fecha") + ylab("O₃ (µg/m³)") +
labs(color = "") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"))
ggplotly(p_comp)El gráfico de comparación entre el pronóstico y los valores reales de diciembre de 2018 permite evaluar visualmente la capacidad predictiva del modelo ARIMA(0,1,1).
Lo primero que llama la atención es que la línea negra de valores reales presenta una variabilidad considerable a lo largo del mes, oscilando entre aproximadamente 16 y 37 µg/m³, mientras que el pronóstico puntual representado por la línea azul punteada permanece constante en 21.09 µg/m³. Esta diferencia es completamente esperada dado que el modelo ARIMA(0,1,1), al no tener componente autorregresivo, no puede anticipar las fluctuaciones diarias de la serie y converge inmediatamente hacia un valor constante como mejor estimación del futuro.
Sin embargo, el aspecto más relevante del gráfico es que los 24 valores reales observados se mantienen dentro del intervalo de confianza al 95% durante todo el período pronosticado, sin que ningún valor real supere las bandas establecidas por el modelo. Esto indica que aunque el pronóstico puntual no replica la dinámica diaria de la serie, el modelo calibra correctamente la incertidumbre asociada al pronóstico que es precisamente el objetivo fundamental de un intervalo de confianza bien construido.
Se observa que durante la segunda y tercera semana de diciembre los valores reales tienden a alejarse más del pronóstico puntual, con un pico notable alrededor del 17 y 18 de diciembre donde la concentración real superó los 35 µg/m³. Este episodio podría estar asociado a condiciones meteorológicas particulares de esos días mayor radiación solar o menor velocidad del viento que favorecieron puntualmente la acumulación de ozono, factores que un modelo univariado como el ARIMA no puede anticipar al no incorporarlos como variables explicativas.
# Métricas de error
accuracy_df <- data.frame(
Métrica = c("MAE", "RMSE", "MAPE (%)"),
Valor = c(
round(mean(abs(comparacion_df$real - comparacion_df$prediccion)), 4),
round(sqrt(mean((comparacion_df$real - comparacion_df$prediccion)^2)), 4),
round(mean(abs((comparacion_df$real - comparacion_df$prediccion) /
comparacion_df$real)) * 100, 4)
)
)
accuracy_df %>%
kable(caption = "Métricas de error - Pronóstico vs Valores reales",
align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "center")| Métrica | Valor |
|---|---|
| MAE | 6.8491 |
| RMSE | 8.1810 |
| MAPE (%) | 23.4816 |
Sobre el comportamiento de la serie:
El análisis de la concentración diaria de O₃ durante 2018 reveló un comportamiento marcadamente estacional, coherente con la dinámica climática del Valle del Cauca. La serie presentó una media de 29.52 µg/m³ y un coeficiente de variación de 25.68%, reflejando una variabilidad moderada-alta característica de los contaminantes fotoquímicos. Se identificaron tres fases claramente diferenciadas: valores moderados entre marzo y mayo, un incremento gradual y sostenido entre junio y septiembre coincidiendo con la temporada seca principal cuando la mayor disponibilidad de radiación solar acelera las reacciones fotoquímicas entre el NO₂ vehicular y los COV alcanzando el máximo absoluto de 58.68 µg/m³ el 27 de agosto, y un descenso progresivo en octubre y noviembre asociado al inicio de la segunda temporada lluviosa. En ningún momento del año las concentraciones superaron el límite normativo de 100 µg/m³ establecido por la Resolución 2254 de 2017, lo que indica condiciones de calidad del aire aceptables para este contaminante durante el período analizado.
Sobre la preparación de los datos:
Una decisión metodológica clave en este análisis fue la exclusión de los meses de enero y febrero de la ventana de entrenamiento, debido a la presencia de 20 días consecutivos completamente sin registro en febrero. Esta decisión, aunque implica una reducción en el número de observaciones de entrenamiento, resultó determinante para garantizar la confiabilidad del modelo: al trabajar únicamente con datos reales sin segmentos artificialmente suavizados por interpolación lineal, el modelo ajustado resultó más parsimonioso, con mejores métricas de diagnóstico y residuales que se ajustan de forma casi perfecta a una distribución normal. Este resultado ilustra un principio fundamental del modelamiento estadístico: la calidad de los datos prima sobre la cantidad.
Sobre la estacionariedad y la identificación del modelo:
La serie de O₃ diario demostró ser no estacionaria en su forma original, resultado confirmado de forma consistente por la inspección visual, la ACF con caída lenta, la prueba ADF (estadístico = -2.1635, p-valor = 0.5071) y la prueba KPSS (estadístico = 0.5650, superior al valor crítico al 5%). Una única diferenciación de orden 1 fue suficiente para lograr la estacionariedad, la prueba ADF sobre la serie diferenciada (estadístico = -9.186, p-valor = 0.01). El análisis de la ACF y PACF de la serie diferenciada orientó la identificación hacia modelos con componente MA(1) y posible estructura AR de orden 1 a 3.
Sobre el modelo seleccionado:
La comparación de tres modelos candidatos: ARIMA(0,1,1), ARIMA(1,1,1) y ARIMA(3,1,1), mediante el criterio AICc determinó que el ARIMA(0,1,1) es el modelo más adecuado (AICc = 1727.167), resultado que coincide con la recomendación de “Auto-Arima” en modo exhaustivo. Este modelo establece que el cambio diario en la concentración de O₃ depende únicamente del error del día inmediatamente anterior, sin necesidad de incorporar valores pasados de la propia serie. La interpretación ambiental de este resultado es coherente: la concentración de O₃ en Cali responde principalmente a perturbaciones aleatorias de corto plazo, como episodios puntuales de alta radiación, cambios en el viento o eventos de lluvia más que a una estructura autorregresiva determinista de largo alcance.
El diagnóstico de residuales validó plenamente el modelo: los residuales se comportan como ruido blanco (Ljung-Box, p-valor = 0.4054), siguen una distribución normal de forma casi perfecta (Shapiro-Wilk, estadístico = 0.9964, p-valor = 0.7918) y no presentan estructura temporal remanente en su ACF. Estos resultados confirman que el modelo captura adecuadamente la dinámica de la serie sin dejar información sin explotar.
Sobre el pronóstico y su evaluación:
El modelo ARIMA(0,1,1) proyectó una concentración constante de 21.09 µg/m³ para los 24 días de diciembre de 2018, valor coherente con el descenso histórico observado en este período. La evaluación frente a los valores reales arrojó un MAE de 6.85 µg/m³, un RMSE de 8.18 µg/m³ y un MAPE de 23.48%. Si bien el pronóstico puntual no replica la variabilidad diaria de la serie limitación inherente a los modelos ARIMA univariados en horizontes amplios los 24 valores reales observados se mantuvieron dentro del intervalo de confianza al 95% durante todo el período, demostrando que el modelo calibra correctamente la incertidumbre asociada al pronóstico. Desde el punto de vista de la pregunta de investigación planteada, el modelo ARIMA sí es posible y estadísticamente válido para pronosticar el O₃ diario en Cali, con una precisión moderada que puede mejorarse incorporando variables meteorológicas exógenas.
Limitaciones y recomendaciones:
El principal limitante del modelo es su naturaleza univariada: al considerar únicamente los valores pasados del O₃, no incorpora variables meteorológicas como la radiación solar, la velocidad del viento o la precipitación, que son determinantes en la formación fotoquímica del contaminante y explican gran parte de la variabilidad diaria no capturada por el modelo. Para investigaciones futuras se recomienda explorar modelos que incorporen estas variables como regresores externos, o modelos de aprendizaje automático como redes neuronales recurrentes (LSTM) capaces de capturar relaciones no lineales entre el O₃ y sus factores precursores. Adicionalmente, contar con series de datos más largas y continuas idealmente varios años de registros sin interrupciones que permitiría identificar patrones interanuales y construir modelos con mayor capacidad predictiva a largo plazo.
Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
Departamento Administrativo de Gestión del Medio Ambiente – DAGMA. (2018). Informe de calidad del aire Santiago de Cali. Alcaldía de Santiago de Cali.
Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. https://otexts.com/fpp3/
Instituto de Hidrología, Meteorología y Estudios Ambientales – IDEAM. (2018). Boletín climatológico mensual. Ministerio de Ambiente y Desarrollo Sostenible de Colombia.
Ministerio de Ambiente y Desarrollo Sostenible. (2017). Resolución 2254 de 2017: Por la cual se adopta la norma de calidad del aire ambiente y se dictan otras disposiciones. Diario Oficial de Colombia.
R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
Hyndman, R. J. (2024). forecast: Forecasting functions for time series and linear models. R package version 8.21. https://pkg.robjhyndman.com/forecast/