Análisis y modelación ARIMA de la concentración de ozono (O3) en Santiago de Cali durante el año 2018

Introducción

La contaminación atmosférica es uno de los principales problemas ambientales de las ciudades contemporáneas por sus efectos sobre la salud pública, los ecosistemas y la calidad de vida. Dentro de los contaminantes criterio, el ozono troposférico (O₃) ocupa un lugar destacado por su naturaleza secundaria: no se emite directamente, sino que se forma a partir de reacciones fotoquímicas entre óxidos de nitrógeno (NOx) y compuestos orgánicos volátiles (COV) en presencia de radiación solar. Esta dependencia de las condiciones meteorológicas hace que sus concentraciones varíen de forma marcada a lo largo del día y del año.

Santiago de Cali cuenta con el Sistema de Vigilancia de la Calidad del Aire de Santiago de Cali (SVCASC), administrado por el Departamento Administrativo de Gestión del Medio Ambiente (DAGMA), que opera una red de estaciones de monitoreo continuo de contaminantes y variables meteorológicas. Entre ellas, la estación Compartir, ubicada en la zona oriental de la ciudad, registra de forma permanente las concentraciones de ozono.

El presente trabajo analiza la serie de concentración de ozono (O₃) registrada en la estación Compartir durante el año 2018. La información original tiene periodicidad horaria, por lo que se construyó una serie de frecuencia diaria a partir del promedio de las observaciones de cada día, con el fin de reducir la variabilidad intradiaria y resaltar los patrones temporales de mediano plazo. El objetivo es doble: (i) caracterizar de manera descriptiva el comportamiento del ozono diario en 2018 y (ii) ajustar un modelo ARIMA que permita modelar y pronosticar la serie. Para ello, el documento se organiza en una etapa descriptiva (estadísticos de resumen, distribución y variabilidad mensual), una etapa de modelación bajo la metodología de Box-Jenkins (identificación, estimación, validación y diagnóstico) y una etapa de pronóstico evaluada contra datos observados reservados como conjunto de prueba.

Metodología

El estudio adopta un enfoque cuantitativo de tipo descriptivo y de modelación de series de tiempo univariadas, a partir de datos secundarios de carácter oficial.

Fuente y descripción de los datos. Se utilizó la serie de concentración de ozono troposférico (O₃), expresada en microgramos por metro cúbico (µg/m³), registrada en la estación Compartir del SVCASC (DAGMA) durante el año 2018. La base original contiene la marca temporal (Fecha & Hora) y la concentración horaria de O₃ (O3 (ug/m3)).

Procesamiento y depuración. A partir de la información horaria se calculó el promedio diario de O₃. Durante la depuración se conservaron únicamente los días con registro válido, descartando los días sin observaciones (valores faltantes), lo que representó más del 94 % del periodo analizado y garantiza una representatividad temporal adecuada. La serie resultante se organizó como un objeto de serie de tiempo (ts) para su posterior modelación.

Análisis descriptivo. Se calcularon estadísticos de resumen (media, mediana, cuartiles, mínimo y máximo) y se construyeron visualizaciones interactivas: un diagrama de caja general para evaluar la dispersión y los valores atípicos, diagramas de caja mensuales para identificar patrones estacionales, y la serie temporal con una tendencia suavizada (LOESS) para describir su evolución a lo largo del año.

Modelación ARIMA (metodología Box-Jenkins). El procedimiento siguió cuatro fases: (i) identificación, evaluando la estacionariedad de la serie mediante la prueba Aumentada de Dickey-Fuller (ADF), determinando el orden de diferenciación necesario y analizando las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie diferenciada; (ii) estimación, ajustando y comparando varios modelos candidatos —ARIMA(0,1,0), (0,1,1), (1,1,0), (1,1,1) y (2,1,1)— mediante el criterio de información de Akaike corregido (AICc) y el principio de parsimonia; (iii) validación y diagnóstico, verificando la significancia de los coeficientes y el cumplimiento de los supuestos sobre los residuales (ausencia de autocorrelación y comportamiento aproximadamente normal de ruido blanco); y (iv) pronóstico, reservando los datos de enero a noviembre de 2018 como conjunto de entrenamiento y diciembre como conjunto de prueba (horizonte h = 24 días) para comparar las predicciones con los valores reales observados.

Herramientas. Todo el procesamiento, análisis y visualización se realizó en el lenguaje R dentro de un documento reproducible en RMarkdown, empleando, entre otros, los paquetes tidyverse y lubridate (manipulación de datos y fechas), readxl (lectura del archivo de Excel), tseries y urca (pruebas de estacionariedad), forecast (estimación y pronóstico ARIMA) y ggplot2, plotly y gridExtra (visualización interactiva y diagnóstico).

Estadísticas descriptivas y narración histórica de la serie

La serie analizada corresponde a la concentración promedio diaria de ozono (O₃) registrada en la estación Compartir del Sistema de Vigilancia de Calidad del Aire de Santiago de Cali durante el año 2018. La información original presentaba periodicidad horaria; sin embargo, para facilitar el análisis se calculó el promedio diario de las observaciones disponibles, obteniéndose un total de 346 registros válidos. Durante la etapa de depuración se identificaron algunos días sin información disponible debido a la ausencia de registros en la base de datos. No obstante, la cantidad de observaciones conservadas representa más del 94 % del periodo analizado, por lo que la serie resultante mantiene una adecuada representatividad temporal. Los principales estadísticos descriptivos obtenidos para la concentración promedio diaria de ozono se presentan en la Tabla 1.

library(tidyverse)
library(lubridate)
library(readxl)
library(tseries)
library(forecast)
library(ggplot2)
library(plotly)
library(ggfortify)
library(urca)
library(xts)
library(gridExtra)

datos <- read_excel("Compartir2.xlsx", sheet = "Hoja1")

Carga y limpieza de datos

datos <- read_excel("Compartir2.xlsx", sheet = "Hoja1") %>%
  rename(Fecha_Hora = `Fecha & Hora`,
         O3 = `O3  (ug/m3)`) %>%
  mutate(Fecha_Hora = ymd_hms(Fecha_Hora)) %>%
  filter(year(Fecha_Hora) == 2018)

datos_diarios <- datos %>%
  mutate(Dia = as.Date(Fecha_Hora)) %>%
  group_by(Dia) %>%
  summarise(
    O3_diario = mean(O3, na.rm = TRUE),
    n_horas = n(),
    n_NA = sum(is.na(O3))
  ) %>%
  filter(!is.na(O3_diario))

datos_345 <- datos_diarios %>%
  filter(!is.na(O3_diario))

serie_345 <- datos_345$O3_diario
serie_ts <- ts(serie_345, start = c(2018, 1), frequency = 1)

cat("Observaciones:", length(serie_345), "\n")

## Observaciones: 345

Análisis descriptivo

La serie analizada corresponde a la concentración promedio diaria de ozono (O₃) registrada en la estación Compartir del Sistema de Vigilancia de Calidad del Aire de Santiago de Cali durante el año 2018. La información original presentaba periodicidad horaria; sin embargo, para facilitar el análisis se calculó el promedio diario de las observaciones disponibles, obteniéndose un total de 346 registros válidos. Durante la etapa de depuración se identificaron algunos días sin información disponible debido a la ausencia de registros en la base de datos. No obstante, la cantidad de observaciones conservadas representa más del 94 % del periodo analizado, por lo que la serie resultante mantiene una adecuada representatividad temporal. Los principales estadísticos descriptivos obtenidos para la concentración promedio diaria de ozono se presentan en la Figura 1.

Distibución del O3

media <- round(mean(datos_diarios$O3_diario, na.rm = TRUE), 1)
mediana <- round(median(datos_diarios$O3_diario, na.rm = TRUE), 1)
q1 <- round(quantile(datos_diarios$O3_diario, 0.25, na.rm = TRUE), 1)
q3 <- round(quantile(datos_diarios$O3_diario, 0.75, na.rm = TRUE), 1)
min_val <- round(min(datos_diarios$O3_diario, na.rm = TRUE), 1)
max_val <- round(max(datos_diarios$O3_diario, na.rm = TRUE), 1)

# Boxplot
g2 <- ggplot(datos_diarios,
             aes(x = "", 
                 y = O3_diario,
                 text = paste(
                   "O3:", round(O3_diario, 1), "µg/m³",
                   "<br>─────────────────",
                   "<br>Media:", media, "µg/m³",
                   "<br>Mediana:", mediana, "µg/m³",
                   "<br>Q1:", q1, "µg/m³",
                   "<br>Q3:", q3, "µg/m³",
                   "<br>Mínimo:", min_val, "µg/m³",
                   "<br>Máximo:", max_val, "µg/m³"
                 ))) +
  geom_boxplot(fill = "#9E4242", alpha = 0.6, color = "#6B2D2D",
               outlier.color = "#9E4242", outlier.alpha = 0.3, outlier.size = 1) +
  geom_jitter(width = 0.05, alpha = 0.15, color = "#9E4242", size = 0.8) +
  labs(title = "Distribución de O3 — Cali 2018",
       subtitle = paste("n =", nrow(datos_diarios), "días"),
       y = "O3 (µg/m³)",
       x = "") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 11, color = "gray40"))

# Convertir a interactivo
ggplotly(g2, tooltip = "text")

Se evidencia que la mayor parte de los registros se concentra alrededor de los 30 µg/m³. Asimismo, se identifican algunos valores atípicos tanto inferiores como superiores, destacándose un valor mínimo cercano a 5.50 µg/m³ y un valor máximo cercano a 58.68 µg/m³. Estos registros representan episodios poco frecuentes dentro del comportamiento anual de la serie.

Variabilidad mensual

datos_diarios <- datos_diarios %>%
  mutate(mes = month(Dia, label = TRUE, abbr = TRUE))

g3 <- ggplot(datos_diarios,
             aes(x = mes, 
                 y = O3_diario,
                 fill = mes,
                 text = paste(
                   "Mes:", mes,
                   "<br>O3:", round(O3_diario, 1), "µg/m³"
                 ))) +
  geom_boxplot(alpha = 0.7, color = "gray30", outlier.alpha = 0.2) +
  scale_fill_manual(values = rep("#9E4242", 12)) +
  labs(title = "Distribución mensual de O3 — Cali 2018",
       x = "Mes",
       y = "O3 (µg/m³)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

ggplotly(g3, tooltip = "text")

Se permite observar que agosto y septiembre presentan no solo las mayores concentraciones promedio, sino también una mayor dispersión de los datos y la presencia de observaciones extremas. En contraste, mayo y noviembre muestran concentraciones más bajas y una menor variabilidad. Con el fin de validar los resultados obtenidos, se comparó el comportamiento mensual de la serie analizada con la información publicada por el Sistema de Vigilancia de Calidad del Aire de Santiago de Cali (SVCASC) en su Informe Anual de Calidad del Aire 2018.

Aunque los indicadores no son exactamente equivalentes, ya que el DAGMA reporta el ozono máximo mensual calculado mediante promedio móvil de 8 horas, mientras que en este trabajo se empleó el promedio diario de la concentración de O₃, ambos análisis presentan una evolución temporal similar. En particular, se observa un incremento progresivo entre junio y septiembre, seguido de una disminución durante octubre y noviembre y una recuperación en diciembre. Esta coincidencia constituye una evidencia de consistencia entre la serie utilizada y la información oficial reportada para la estación Compartir. Adicionalmente, la coincidencia entre el nombre de la base de datos, las variables registradas, el periodo analizado y la evolución temporal observada respecto a los reportes oficiales del SVCASC sugiere que la información utilizada corresponde a registros provenientes de la estación Compartir o de la misma red de monitoreo utilizada por el DAGMA.

En términos generales, los resultados descriptivos indican que la concentración de ozono en Cali durante 2018 presentó una variabilidad moderada y un comportamiento temporal claramente diferenciado a lo largo del año. Los niveles más elevados se concentraron entre julio y septiembre, mientras que los valores más bajos se registraron principalmente durante mayo y noviembre. Estos hallazgos son consistentes con la información reportada por los boletines mensuales y el informe anual de calidad del aire del Sistema de Vigilancia de Calidad del Aire de Santiago de Cali, lo que proporciona confianza sobre la representatividad de la serie utilizada para el posterior desarrollo del modelo ARIMA.

Comportamiento temporal de la serie

media_O3 <- mean(datos_345$O3_diario, na.rm = TRUE)

plot_ly() %>%
  add_trace(
    data = datos_345,
    x = ~Dia,
    y = ~O3_diario,
    type = 'scatter',
    mode = 'lines',
    line = list(color = '#9E4242', width = 2),
    name = 'O3 diario',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>O3:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  add_trace(
    data = datos_345,
    x = ~Dia,
    y = ~media_O3,
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'red', width = 1.5, dash = 'dash'),
    name = paste('Media =', round(media_O3, 1), 'µg/m³'),
    hovertemplate = paste(
      '<b>Media:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  add_trace(
    data = datos_345,
    x = ~Dia,
    y = ~loess(O3_diario ~ as.numeric(Dia), span = 0.1)$fitted,
    type = 'scatter',
    mode = 'lines',
    line = list(color = '#CB7DFF', width = 2),
    name = 'Tendencia suavizada',
    hovertemplate = paste(
      '<b>Tendencia:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  layout(
    title = list(
      text = "Ozono (O3) - Promedio diario - Cali 2018",
      font = list(size = 16)
    ),
    xaxis = list(
      title = "Fecha",
      gridcolor = "lightgray",
      rangeslider = list(visible = TRUE)
    ),
    yaxis = list(
      title = "O3 (µg/m³)",
      gridcolor = "lightgray"
    ),
    hoverlabel = list(
      bgcolor = "white",
      font = list(size = 12)
    ),
    legend = list(
      orientation = "h",
      yanchor = "bottom",
      y = 1.02,
      xanchor = "center",
      x = 0.5
    ),
    margin = list(t = 80, r = 30, b = 50, l = 60)
  )

Se observa que la serie presenta fluctuaciones continuas durante todo el año, sin mantener niveles constantes de concentración. Durante los primeros meses del año (enero a abril), los valores oscilaron principalmente entre 20 y 35 µg/m³, con algunos incrementos puntuales cercanos a 45 µg/m³. Posteriormente, entre mayo y junio se evidencia una disminución de las concentraciones promedio diarias, registrándose algunos de los valores más bajos del periodo de estudio. Este comportamiento coincide con los boletines mensuales de calidad del aire emitidos por el DAGMA para dichos meses, en los cuales se realiza seguimiento al comportamiento temporal de los contaminantes atmosféricos monitoreados por el SVCASC. A partir de julio comienza un incremento progresivo de las concentraciones de ozono que se mantiene durante agosto y septiembre. Este periodo corresponde a los niveles más altos observados en la serie, alcanzando valores máximos cercanos a 59 µg/m³. El comportamiento observado es consistente con los boletines mensuales de julio, agosto y septiembre, donde se reporta el seguimiento de este contaminante y sus excedencias dentro del sistema de vigilancia de la calidad del aire de Cali. Después de septiembre se observa una disminución progresiva de las concentraciones durante octubre y noviembre, seguida de una ligera recuperación hacia diciembre. Este comportamiento también es coherente con la información reportada en los boletines mensuales correspondientes al último trimestre del año.

Modelación ARIMA

Ventana de entrenamiento del modelo ARIMA

fecha_limite <- ymd("2018-11-30")

datos_entrenamiento <- datos_345 %>%
  filter(Dia <= fecha_limite)

media_O3_entrenamiento <- mean(datos_entrenamiento$O3_diario, na.rm = TRUE)

loess_fit <- loess(O3_diario ~ as.numeric(Dia), 
                   data = datos_entrenamiento, 
                   span = 0.1)
datos_entrenamiento$tendencia <- predict(loess_fit)

plot_ly() %>%
  add_trace(
    data = datos_entrenamiento,
    x = ~Dia,
    y = ~O3_diario,
    type = 'scatter',
    mode = 'lines',
    line = list(color = '#9E4242', width = 2),
    name = 'O3 diario (entrenamiento)',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>O3:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  layout(
    title = list(
      text = "Ozono (O3) - Entrenamiento: Enero a Noviembre 2018",
      font = list(size = 16)
    ),
    xaxis = list(
      title = "Fecha",
      gridcolor = "lightgray",
      rangeslider = list(visible = TRUE),
      range = c(as.Date("2018-01-01"), as.Date("2018-11-30"))
    ),
    yaxis = list(
      title = "O3 (µg/m³)",
      gridcolor = "lightgray"
    ),
    hoverlabel = list(
      bgcolor = "white",
      font = list(size = 12)
    ),
    legend = list(
      orientation = "h",
      yanchor = "bottom",
      y = 1.02,
      xanchor = "center",
      x = 0.5
    ),
    margin = list(t = 80, r = 30, b = 50, l = 60)
  )

La serie de ozono diario en Cali durante 2018 oscila entre aproximadamente 13 y 59 µg/m³, sin una tendencia lineal clara a lo largo del año, lo que es coherente con la necesidad de una diferenciación de orden 1 para estacionarizar la serie. Sin embargo, sí se distinguen fluctuaciones asociadas al régimen climático bimodal de la ciudad: los niveles tienden a ser más altos en los períodos secos (enero–febrero y junio–agosto), cuando mayor radiación solar favorece la formación fotoquímica del ozono, y más bajos durante las temporadas de lluvias (marzo–mayo y octubre–noviembre), cuando la nubosidad y la precipitación inhiben este proceso. El pico máximo (~59 µg/m³) se registra en septiembre, posiblemente bajo condiciones puntuales de alta temperatura y baja humedad.

La serie presenta además una alta volatilidad día a día, con variaciones de hasta 20 µg/m³ entre días consecutivos, lo que se refleja en una desviación estándar de ~14 µg/m³ (cerca del 44% de la media). Este comportamiento sugiere que el ozono en Cali responde principalmente a condiciones meteorológicas de corto plazo, más que a dinámicas de mediano o largo plazo, lo cual justifica un modelo de memoria corta como el ARIMA(0,1,1). No obstante, la estacionalidad intra-anual visible en el gráfico no es capturada por este modelo, lo que representa una limitación a considerar.

Prueba ADF

library(plotly)
library(dplyr)

df_comparacion <- data.frame(
  Dia = datos_345$Dia[2:nrow(datos_345)],
  Original = serie_ts[2:length(serie_ts)],
  Diferenciada = as.numeric(diff(serie_ts))
)

plot_ly() %>%
  add_trace(
    data = df_comparacion,
    x = ~Dia,
    y = ~Original,
    type = 'scatter',
    mode = 'lines',
    line = list(color = '#9E4242', width = 2),
    name = 'Serie original',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>O3:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  add_trace(
    data = df_comparacion,
    x = ~Dia,
    y = ~Diferenciada,
    type = 'scatter',
    mode = 'lines',
    line = list(color = '#9B62A6', width = 2),
    name = 'Serie diferenciada',
    yaxis = 'y2',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>ΔO3:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  layout(
    title = list(
      text = "Comparación: Serie original vs Diferenciada",
      font = list(size = 16)
    ),
    xaxis = list(
      title = "Fecha",
      gridcolor = "lightgray",
      rangeslider = list(visible = TRUE)
    ),
    yaxis = list(
      title = "O3 (µg/m³)",
      gridcolor = "lightgray",
      color = '#9E4242'
    ),
    yaxis2 = list(
      title = "ΔO3 (µg/m³)",
      overlaying = 'y',
      side = 'right',
      gridcolor = "lightgray",
      color = '#1E90FF'
    ),
    legend = list(
      orientation = "h",
      yanchor = "bottom",
      y = 1.02,
      xanchor = "center",
      x = 0.5
    ),
    hoverlabel = list(
      bgcolor = "white",
      font = list(size = 12)
    ),
    margin = list(t = 80, r = 60, b = 50, l = 60)
  )

El gráfico superpone la serie original de O3 (rojo) y su primera diferencia ΔO3 (morado), permitiendo evaluar visualmente el efecto de la diferenciación. La serie original mantiene los patrones ya descritos: niveles cambiantes a lo largo del año, alta volatilidad día a día y ausencia de tendencia lineal sostenida, pero con una media que oscila según el régimen climático, lo que confirma su no estacionariedad en niveles.

La serie diferenciada, en cambio, oscila de forma estable alrededor de cero durante todo el período, sin desplazamientos sostenidos hacia arriba o hacia abajo. Esto indica que una sola diferenciación es suficiente para remover la no estacionariedad de la serie, justificando el parámetro d=1 del modelo ARIMA(0,1,1). La amplitud de las fluctuaciones de ΔO3 se mantiene relativamente constante a lo largo del año, con valores que raramente superan ±20 µg/m³, lo cual es consistente con una varianza estable, condición necesaria para la validez del modelo.

Vale notar que los picos extremos de la serie diferenciada (visibles alrededor de marzo y septiembre) corresponden precisamente a las transiciones abruptas observadas en la serie original: la caída brusca hacia valores mínimos en febrero–marzo y el pico máximo de septiembre. Estos choques puntuales son captados por el término MA(1) del modelo, que precisamente modela la influencia del error del día anterior sobre el valor actual, corrigiendo en un 74% la magnitud de dichos saltos al día siguiente.

tabla_adf <- data.frame(
  Serie = c("Original (O3 diario)", "Diferenciada (ΔO3)"),
  `Estadístico ADF` = c(-2.5417, -11.157),
  `Valor p` = c(0.349, 0.01),
  `Conclusión` = c("No estacionaria", "Estacionaria")
)

knitr::kable(tabla_adf, 
             caption = "Resultados de la prueba ADF para estacionariedad",
             digits = 4)

Resultados de la prueba ADF para estacionariedad

Serie	Estadístico.ADF	Valor.p	Conclusión
Original (O3 diario)	-2.5417	0.349	No estacionaria
Diferenciada (ΔO3)	-11.1570	0.010	Estacionaria

La prueba Aumentada de Dickey-Fuller (ADF) se aplicó para evaluar la estacionariedad de la serie de ozono diario en Cali durante 2018. Esta prueba tiene como hipótesis nula que la serie presenta una raíz unitaria, es decir, que no es estacionaria. Para la serie original, se obtuvo un estadístico ADF de -2.5417 con un valor p de 0.349, lo que indica que no se rechaza la hipótesis nula al nivel de significancia del 5%. Este resultado confirma que la serie original no es estacionaria, comportamiento esperado dadas las fluctuaciones estacionales asociadas al régimen climático bimodal de Cali y la alta variabilidad diaria observada en los datos.

Posteriormente, se aplicó una diferenciación de orden 1 a la serie y se repitió la prueba ADF. En este caso, el estadístico ADF fue de -11.157 con un valor p de 0.01, lo que permite rechazar la hipótesis nula y concluir que la serie diferenciada es estacionaria. La diferenciación logró eliminar la tendencia y estabilizar la media, generando una serie que oscila alrededor de cero con varianza aproximadamente constante. Este resultado justifica la inclusión del componente integrador I(1) en el modelo ARIMA(0,1,1), que fue seleccionado automáticamente por el criterio AICc. La evidencia estadística confirma que el ozono en Cali responde principalmente a dinámicas de corto plazo, lo cual valida la estructura del modelo elegido y su capacidad para generar pronósticos confiables.

ACF y PACF de la serie diferenciada

serie_diff <- diff(serie_ts)

acf_diff <- acf(serie_diff, lag.max = 30, plot = FALSE)
pacf_diff <- pacf(serie_diff, lag.max = 30, plot = FALSE)

df_acf <- data.frame(
  Lag = acf_diff$lag[, , 1],
  ACF = acf_diff$acf[, , 1]
)

df_pacf <- data.frame(
  Lag = pacf_diff$lag[, , 1],
  PACF = pacf_diff$acf[, , 1]
)

limite <- 1.96 / sqrt(length(serie_diff))

acf_plot <- plot_ly() %>%
  add_trace(
    data = df_acf,
    x = ~Lag,
    y = ~ACF,
    type = 'bar',
    marker = list(color = '#9E4242'),
    name = 'ACF',
    hovertemplate = paste(
      '<b>Lag:</b> %{x}<br>',
      '<b>ACF:</b> %{y:.3f}<br>',
      '<extra></extra>'
    )
  ) %>%
  add_trace(
    x = c(0, max(df_acf$Lag)),
    y = c(limite, limite),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'blue', dash = 'dash'),
    name = 'Límite superior',
    showlegend = FALSE
  ) %>%
  add_trace(
    x = c(0, max(df_acf$Lag)),
    y = c(-limite, -limite),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'blue', dash = 'dash'),
    name = 'Límite inferior',
    showlegend = FALSE
  ) %>%
  layout(
    title = "ACF - Serie diferenciada",
    xaxis = list(title = "Lag"),
    yaxis = list(title = "ACF", range = c(-1, 1))
  )

pacf_plot <- plot_ly() %>%
  add_trace(
    data = df_pacf,
    x = ~Lag,
    y = ~PACF,
    type = 'bar',
    marker = list(color = '#1E90FF'),
    name = 'PACF',
    hovertemplate = paste(
      '<b>Lag:</b> %{x}<br>',
      '<b>PACF:</b> %{y:.3f}<br>',
      '<extra></extra>'
    )
  ) %>%
  add_trace(
    x = c(0, max(df_pacf$Lag)),
    y = c(limite, limite),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'blue', dash = 'dash'),
    showlegend = FALSE
  ) %>%
  add_trace(
    x = c(0, max(df_pacf$Lag)),
    y = c(-limite, -limite),
    type = 'scatter',
    mode = 'lines',
    line = list(color = 'blue', dash = 'dash'),
    showlegend = FALSE
  ) %>%
  layout(
    title = "ACF y PACF de la serie diferenciada",
    xaxis = list(title = "Lag"),
    yaxis = list(title = "PACF", range = c(-1, 1))
  )

subplot(acf_plot, pacf_plot, nrows = 1, titleX = TRUE, titleY = TRUE)

El ACF de la serie diferenciada muestra un corte abrupto tras el lag 1, con una autocorrelación negativa considerable (~−0.5) que cae inmediatamente dentro de las bandas de confianza a partir del lag 2. Los lags posteriores oscilan de forma irregular alrededor de cero, con algún spike aislado cerca del lag 27–28 que apenas toca el límite, pero sin un patrón sistemático que sugiera estructura adicional. Este comportamiento es la firma clásica de un proceso MA(1), donde solo el error del período inmediatamente anterior tiene efecto significativo sobre el valor actual.

El PACF, por su parte, muestra un decaimiento gradual desde el lag 1, con varios coeficientes negativos significativos en los primeros lags (especialmente el lag 1, que supera −0.5) antes de estabilizarse dentro de las bandas. Este patrón de decaimiento paulatino —en contraste con el corte brusco del ACF— es precisamente lo que distingue un proceso MA puro de uno AR, confirmando que no se requiere ningún componente autorregresivo en el modelo.

En conjunto, la estructura ACF con corte en lag 1 y PACF con decaimiento gradual constituye la evidencia gráfica que justifica la selección del modelo ARIMA(0,1,1): una diferenciación para lograr estacionariedad y un término de media móvil de orden 1 para capturar la única dependencia significativa que persiste en la serie diferenciada.

Comparación de modelos ARIMA

modelo_010 <- Arima(serie_ts, order = c(0,1,0)) 
modelo_011 <- Arima(serie_ts, order = c(0,1,1))  
modelo_110 <- Arima(serie_ts, order = c(1,1,0)) 
modelo_111 <- Arima(serie_ts, order = c(1,1,1)) 
modelo_211 <- Arima(serie_ts, order = c(2,1,1))  



tabla_modelos <- data.frame(
  Modelo = c("ARIMA(0,1,0)", "ARIMA(0,1,1)", "ARIMA(1,1,0)", 
             "ARIMA(1,1,1)", "ARIMA(2,1,1)"),
  AICc = c(modelo_010$aicc, modelo_011$aicc, modelo_110$aicc,
           modelo_111$aicc, modelo_211$aicc),
  LogLik = c(modelo_010$loglik, modelo_011$loglik, modelo_110$loglik,
             modelo_111$loglik, modelo_211$loglik),
  Parámetros = c(1, 2, 2, 3, 4)
)

tabla_modelos <- tabla_modelos %>%
  arrange(AICc)

knitr::kable(tabla_modelos, 
             caption = "Comparación de modelos ARIMA según AICc",
             digits = 2)

Comparación de modelos ARIMA según AICc

Modelo	AICc	LogLik	Parámetros
ARIMA(0,1,1)	2188.40	-1092.18	2
ARIMA(2,1,1)	2189.10	-1090.49	4
ARIMA(1,1,1)	2189.11	-1091.52	3
ARIMA(1,1,0)	2237.14	-1116.55	2
ARIMA(0,1,0)	2320.93	-1159.46	1

La tabla presenta la comparación de cinco modelos ARIMA ajustados para la serie de ozono diario en Cali durante 2018, evaluados según el criterio de información de Akaike corregido (AICc), el cual equilibra la bondad de ajuste con la complejidad del modelo, penalizando la inclusión de parámetros adicionales. El modelo ARIMA(0,1,1) presenta el menor AICc (2188.40), lo que lo posiciona como la mejor opción entre los modelos evaluados. Le siguen muy de cerca los modelos ARIMA(2,1,1) con 2189.10 y ARIMA(1,1,1) con 2189.11, cuyas diferencias con respecto al mejor modelo son inferiores a 1 unidad, lo que indica que estadísticamente son equivalentes en términos de ajuste.

Sin embargo, al aplicar el principio de parsimonia, el ARIMA(0,1,1) resulta ser el más adecuado porque logra un ajuste prácticamente igual al de los modelos más complejos pero con un número significativamente menor de parámetros (2 parámetros versus 3 o 4). Esta elección es coherente con el diagnóstico realizado mediante las funciones ACF y PACF sobre la serie diferenciada, que sugerían un proceso de media móvil de orden 1 como estructura subyacente. Los modelos ARIMA(1,1,0), ARIMA(0,1,0) presentan valores de AICc considerablemente más altos (2237.14 y 2320.93 respectivamente), confirmando que la inclusión del componente MA(1) es necesaria para capturar adecuadamente la dinámica de corto plazo del ozono en Cali, y que un modelo de caminata aleatoria o puramente autorregresivo no sería suficiente para explicar el comportamiento de la serie.

Resultados ARIMA (0,1,1)

coefs <- coef(modelo_011)
se <- sqrt(diag(vcov(modelo_011)))

acc <- accuracy(modelo_011)


tabla_completa <- data.frame(
  Modelo = "ARIMA(0,1,1)",
  
  ma1 = paste0(round(coefs[1], 4), " (", round(se[1], 4), ")"),
  
  AICc = round(modelo_011$aicc, 2),
  BIC = round(modelo_011$bic, 2),
  LogLik = round(modelo_011$loglik, 2),
  sigma2 = round(modelo_011$sigma2, 3),
  
  RMSE = round(acc[1, 2], 4),
  MAE = round(acc[1, 3], 4),
  MAPE = round(acc[1, 5], 4),
  MASE = round(acc[1, 6], 4),
  ACF1 = round(acc[1, 7], 4)
)


knitr::kable(tabla_completa, 
             caption = "Resultados del modelo ARIMA(0,1,1) seleccionado",
             align = "c",
             digits = 4)

Resultados del modelo ARIMA(0,1,1) seleccionado

Modelo	ma1	AICc	BIC	LogLik	sigma2	RMSE	MAE	MAPE	MASE	ACF1
ARIMA(0,1,1)	-0.7398 (0.0426)	2188.4	2196.04	-1092.18	33.536	5.7742	4.505	16.3567	0.8242	0.0448

El modelo ARIMA(0,1,1) ajustado para la serie de ozono diario en Cali presenta un coeficiente de media móvil ma1 = -0.7398 (error estándar: 0.0426), el cual es estadísticamente significativo. Este coeficiente indica que un error de pronóstico positivo del día anterior se asocia con una corrección a la baja de aproximadamente el 74% en la concentración de ozono del día siguiente, reflejando un mecanismo de autorregulación en la dinámica del contaminante. Desde una perspectiva ambiental, este comportamiento sugiere que los episodios de alta contaminación tienden a ser seguidos por días de menores concentraciones, posiblemente debido a cambios en las condiciones meteorológicas o a la propia dinámica fotoquímica del ozono, que tiene una vida media relativamente corta en la atmósfera.

Las métricas de ajuste y error confirman la idoneidad del modelo. El AICc de 2188.40 es el más bajo entre los modelos evaluados, lo que indica un equilibrio óptimo entre bondad de ajuste y complejidad. El RMSE de 5.77 µg/m³ y el MAE de 4.51 µg/m³ reflejan que, en promedio, el modelo se equivoca por aproximadamente 4.5 a 5.8 µg/m³ al pronosticar el ozono diario, mientras que el MAPE de 16.36% sugiere un error porcentual moderado, aceptable para este tipo de variables ambientales. El MASE de 0.824, al ser menor que 1, confirma que el modelo es superior a un modelo ingenuo (Naïve) en aproximadamente un 17.6%, demostrando que el componente MA(1) sí aporta valor predictivo. Finalmente, la ACF1 de 0.0448, cercana a cero, indica que los residuales no presentan autocorrelación significativa en el primer rezago, lo que valida la correcta especificación del modelo y respalda la confianza en los pronósticos generados. Estos resultados posicionan al ARIMA(0,1,1) como una herramienta útil para el monitoreo y la gestión de la calidad del aire en Cali.

Diagnóstico de residuales

residuos <- residuals(modelo_011)


df_res <- data.frame(
  tiempo   = 1:length(residuos),
  residuos = as.numeric(residuos)
)


p1 <- ggplot(df_res, aes(x = tiempo, y = residuos)) +
  geom_line(color = "#9E4242", linewidth = 0.4) +
  geom_hline(yintercept = 0, color = "#5C4033", linewidth = 0.5) +
  labs(title = "Residuos del modelo ARIMA(0,1,1)",
       y = "Residuos (µg/m³)", x = NULL) +
  theme_minimal() +
  theme(
    plot.title       = element_text(face = "bold", color = "#5C4033"),
    axis.line        = element_line(color = "#9E4242"),
    panel.grid.major = element_line(color = "grey90")
  )


p2 <- ggAcf(residuos, colour = "steelblue4", lag.max = 30) +
  labs(title = "ACF de residuos") +
  theme_minimal() +
  theme(
    plot.title       = element_text(face = "bold", color = "#5C4033"),
    axis.line        = element_line(color = "#9E4242"),
    panel.grid.major = element_line(color = "grey90")
  )


p3 <- ggplot(df_res, aes(x = residuos)) +
  geom_histogram(fill = "#9E4242", color = "white", bins = 40, alpha = 0.8) +
  geom_vline(xintercept = 0, color = "#5C4033", linewidth = 0.7) +
  labs(title = "Distribución de residuos",
       y = "Frecuencia", x = "Residuos (µg/m³)") +
  theme_minimal() +
  theme(
    plot.title       = element_text(face = "bold", color = "#5C4033"),
    axis.line        = element_line(color = "#9E4242"),
    panel.grid.major = element_line(color = "grey90")
  )


grid.arrange(p1, p2, p3, 
             layout_matrix = rbind(c(1, 1), c(2, 3)))

La gráfica superior de residuos en el tiempo muestra una serie que fluctúa de forma irregular alrededor de cero durante todo el período, sin tendencia visible, sin cambios en la varianza y sin patrones cíclicos reconocibles. Los valores se mantienen en su mayoría dentro del rango ±10 µg/m³, con algunos picos aislados que alcanzan ±15 µg/m³, correspondientes a los choques extremos ya identificados en la serie original (transición febrero–marzo y pico de septiembre). Este comportamiento es consistente con un proceso de ruido blanco.

El ACF de los residuos confirma esta conclusión: todos los coeficientes se encuentran dentro de las bandas de confianza (±0.10), sin ningún spike significativo en ningún lag. No hay estructura de correlación remanente, lo que indica que el modelo ARIMA(0,1,1) ha capturado adecuadamente toda la dependencia temporal de la serie y que no se justifica agregar términos AR o MA adicionales.

El histograma de la distribución de residuos muestra una forma aproximadamente simétrica y centrada en cero, con una ligera asimetría hacia la derecha por los valores extremos positivos. La distribución se asemeja razonablemente a una normal, aunque con colas algo más pesadas de lo esperado, lo que sugiere la presencia de valores atípicos puntuales más que una desviación sistemática. En conjunto, los tres paneles validan el modelo: los residuos son independientes, centrados en cero y aproximadamente normales, cumpliendo los supuestos necesarios para que los intervalos de confianza del pronóstico sean confiables.

Pronóstico

datos_entrenamiento <- datos_diarios %>%
  filter(Dia <= ymd("2018-11-30"))

datos_prueba <- datos_diarios %>%
  filter(Dia >= ymd("2018-12-01"))

serie_entrenamiento <- ts(datos_entrenamiento$O3_diario, 
                          start = c(2018, 1), 
                          frequency = 1)



modelo_entrenamiento <- Arima(serie_entrenamiento, order = c(0,1,1))

h <- 24
pronostico <- forecast(modelo_entrenamiento, h = h)


fechas_pronostico <- seq(ymd("2018-12-01"), by = "day", length.out = h)

df_pronostico <- data.frame(
  Dia = fechas_pronostico,
  Pronostico = as.numeric(pronostico$mean),
  Lo_80 = as.numeric(pronostico$lower[, 1]),
  Hi_80 = as.numeric(pronostico$upper[, 1]),
  Lo_95 = as.numeric(pronostico$lower[, 2]),
  Hi_95 = as.numeric(pronostico$upper[, 2])
)


datos_prueba_24 <- datos_prueba %>%
  filter(Dia <= ymd("2018-12-24"))

plot_ly() %>%
  
  add_trace(
    data = datos_entrenamiento,
    x = ~Dia, y = ~O3_diario,
    type = 'scatter', mode = 'lines',
    line = list(color = '#9E4242', width = 2),
    name = 'Histórico (Ene-Nov)',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>O3:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  
  add_trace(
    data = datos_prueba_24,
    x = ~Dia, y = ~O3_diario,
    type = 'scatter', mode = 'lines+markers',
    line = list(color = '#2E8B57', width = 2),
    marker = list(color = '#2E8B57', size = 5),
    name = 'Real (Diciembre)',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>O3 real:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  
  add_trace(
    data = df_pronostico,
    x = ~Dia, y = ~Pronostico,
    type = 'scatter', mode = 'lines',
    line = list(color = '#9E4242', width = 2, dash = 'dash'),
    name = 'Pronóstico',
    hovertemplate = paste(
      '<b>Fecha:</b> %{x|%d-%b-%Y}<br>',
      '<b>Pronóstico:</b> %{y:.1f} µg/m³<br>',
      '<extra></extra>'
    )
  ) %>%
  
  add_ribbons(
    data = df_pronostico,
    x = ~Dia, ymin = ~Lo_95, ymax = ~Hi_95,
    color = I('rgba(158, 66, 66, 0.15)'),
    name = 'IC 95%',
    hoverinfo = 'skip'
  ) %>%
  
  add_ribbons(
    data = df_pronostico,
    x = ~Dia, ymin = ~Lo_80, ymax = ~Hi_80,
    color = I('rgba(158, 66, 66, 0.25)'),
    name = 'IC 80%',
    hoverinfo = 'skip'
  ) %>%
  
  add_segments(
    x = as.Date("2018-11-30"), xend = as.Date("2018-11-30"),
    y = 0, yend = 60,
    line = list(color = '#5C4033', dash = 'dot', width = 1),
    showlegend = FALSE
  ) %>%
  
  layout(
    title = list(
      text = "Ozono diario en Cali - Pronóstico vs Real (Diciembre 2018)",
      font = list(size = 16)
    ),
    xaxis = list(
      title = "Fecha",
      gridcolor = "lightgray",
      rangeslider = list(visible = TRUE)
    ),
    yaxis = list(
      title = "O3 (µg/m³)",
      gridcolor = "lightgray",
      range = c(0, 60)
    ),
    hoverlabel = list(
      bgcolor = "white",
      font = list(size = 12)
    ),
    legend = list(
      orientation = "h",
      yanchor = "bottom",
      y = 1.02,
      xanchor = "center",
      x = 0.5
    ),
    margin = list(t = 80, r = 30, b = 50, l = 60)
  )

El gráfico compara el pronóstico del modelo ARIMA(0,1,1) contra los valores reales observados durante diciembre de 2018, período que fue reservado como conjunto de prueba. El pronóstico puntual (línea roja discontinua) es completamente plano en aproximadamente 21 µg/m³, comportamiento esperado para un modelo IMA(1,1) sin drift, cuya predicción a más de un paso converge al último valor observado. Los intervalos de confianza al 80% y 95% se expanden progresivamente con el horizonte de pronóstico, reflejando la acumulación de incertidumbre a medida que aumentan los días proyectados.

Los valores reales de diciembre (línea verde) muestran una dinámica claramente distinta: oscilan entre ~20 y ~40 µg/m³ con una tendencia ascendente hacia mediados de mes, antes de caer nuevamente hacia finales de diciembre. Si bien la mayoría de los valores reales se mantienen dentro del intervalo de confianza al 95%, el modelo no logra anticipar ni la dirección ni la magnitud de estas fluctuaciones, lo cual es una limitación inherente al pronóstico puntual plano del ARIMA(0,1,1).

Esta comparación ilustra con claridad la principal restricción del modelo: es adecuado para capturar la estructura de dependencia de corto plazo de la serie, pero no incorpora información sobre las condiciones meteorológicas (radiación, temperatura, humedad) que modulan el ozono en diciembre, un mes que en Cali corresponde al inicio de la segunda época seca y tiende a registrar concentraciones en recuperación.

Conclusiones

Hallazgos

El análisis descriptivo mostró que la concentración de ozono en Cali durante 2018 presentó una variabilidad moderada, con la mayoría de los registros concentrados alrededor de 30 µg/m³ y valores extremos puntuales entre ~5.5 y ~59 µg/m³. Se identificó un patrón estacional consistente con el régimen climático bimodal de la ciudad: las mayores concentraciones promedio se registraron entre julio y septiembre —con la mayor dispersión y presencia de valores extremos en agosto y septiembre—, mientras que los niveles más bajos ocurrieron en mayo y noviembre. Esta evolución temporal coincide con la reportada en los boletines mensuales y en el Informe Anual de Calidad del Aire 2018 del SVCASC, lo que respalda la representatividad de la serie utilizada.

Desde la modelación, la prueba ADF confirmó que la serie original no es estacionaria (p = 0.349) y que una sola diferenciación es suficiente para estacionarizarla (p = 0.01), justificando el componente I(1). El análisis de la ACF (corte en el lag 1) y la PACF (decaimiento gradual) apuntó a una estructura de media móvil de orden 1. Entre los modelos evaluados, el ARIMA(0,1,1) obtuvo el menor AICc (2188.40) y, por parsimonia, resultó el más adecuado: con prácticamente el mismo ajuste que modelos más complejos pero un número menor de parámetros. Su coeficiente ma1 = −0.7398 fue estadísticamente significativo e indica un mecanismo de autorregulación, en el que un episodio de alta concentración tiende a ser corregido a la baja al día siguiente.

El modelo presentó un desempeño aceptable para una variable ambiental (RMSE = 5.77 µg/m³, MAE = 4.51 µg/m³, MAPE = 16.36 % y MASE = 0.824, es decir, cerca de un 17.6 % mejor que un modelo ingenuo). El diagnóstico de residuales confirmó que se comportan como ruido blanco —independientes, centrados en cero y aproximadamente normales—, validando la especificación. En la evaluación contra diciembre de 2018, los valores reales se mantuvieron en su mayoría dentro del intervalo de confianza al 95 %, aunque el pronóstico puntual resultó plano (~21 µg/m³), como es propio de un modelo IMA(1,1) sin deriva.

Limitaciones

El modelo ARIMA(0,1,1), pese a su buen ajuste de corto plazo, presenta limitaciones relevantes. En primer lugar, captura únicamente la dependencia de corto plazo y no incorpora la estacionalidad intra-anual asociada al régimen climático bimodal, visible en los datos pero no modelada. En segundo lugar, el pronóstico puntual es plano y converge al último valor observado, por lo que no anticipa la dirección ni la magnitud de las fluctuaciones reales, como se evidenció en diciembre. En tercer lugar, se trata de un modelo univariado que no incorpora variables exógenas meteorológicas (radiación solar, temperatura, humedad, precipitación), determinantes en la formación fotoquímica del ozono; su inclusión mediante modelos ARIMAX o SARIMA podría mejorar sustancialmente la capacidad predictiva.

Adicionalmente, existen limitaciones en los datos: cerca del 6 % de los días no contaban con registros válidos y fueron excluidos, lo que afecta la continuidad de la serie. La comparación con la información oficial del DAGMA es solo indicativa, ya que el SVCASC reporta el ozono máximo mensual mediante promedio móvil de 8 horas, mientras que aquí se empleó el promedio diario, indicadores no directamente equivalentes. Finalmente, el horizonte de evaluación fue corto (24 días) y los residuales mostraron colas algo más pesadas que las de una normal, atribuibles a episodios atípicos puntuales. Estas limitaciones señalan líneas de mejora para trabajos futuros orientados a la gestión de la calidad del aire en la ciudad.

Bibliografía

• Departamento Administrativo de Gestión del Medio Ambiente (DAGMA). (2018). Boletines mensuales de calidad del aire de Santiago de Cali. Sistema de Vigilancia de la Calidad del Aire de Santiago de Cali (SVCASC). Alcaldía de Santiago de Cali. Recuperado de https://www.cali.gov.co/documentos/271/boletines-calidad-del-aire-dagma/

• Departamento Administrativo de Gestión del Medio Ambiente (DAGMA). (2019). Informe Anual de Calidad del Aire de Santiago de Cali 2018. Sistema de Vigilancia de la Calidad del Aire de Santiago de Cali (SVCASC). Alcaldía de Santiago de Cali. Recuperado de https://www.cali.gov.co/dagma/

• Box, G. E. P., Jenkins, G. M., Reinsel, G. C. & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.

• Peña, D. (2005). Análisis de series temporales. Alianza Editorial, Madrid.

• Hyndman, R. J. & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. https://otexts.com/fpp3/

• Dickey, D. A. & Fuller, W. A. (1979). Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association, 74(366), 427–431.

• R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Viena, Austria. https://www.R-project.org/