La calidad del aire constituye uno de los determinantes ambientales más relevantes para la salud pública en las grandes ciudades de América Latina, y Santiago de Cali no es la excepción. Como tercera ciudad más poblada de Colombia, ubicada en un valle interandino con condiciones meteorológicas particulares —alta radiación solar, baja velocidad del viento y temperaturas elevadas durante buena parte del año—, Cali presenta un escenario propicio para la formación de contaminantes secundarios como el ozono troposférico (O₃) (DAGMA, 2019; OMS, 2021).

A diferencia del ozono estratosférico, que protege la vida en la Tierra, el ozono a nivel del suelo es un contaminante perjudicial que se forma a partir de reacciones fotoquímicas entre óxidos de nitrógeno (NOₓ) y compuestos orgánicos volátiles (COV) en presencia de luz solar (EPA, 2023; OMS, 2021). Por esta razón, el O₃ no se emite directamente, sino que se genera en la atmósfera, alcanzando sus máximas concentraciones en las horas centrales del día, cuando la radiación solar es más intensa. La exposición prolongada a concentraciones elevadas de ozono se asocia con afecciones respiratorias, agravamiento del asma y reducción de la función pulmonar, lo que convierte su monitoreo y pronóstico en una herramienta fundamental para la gestión ambiental urbana (OMS, 2021; EPA, 2023).

El análisis de series de tiempo emerge como una metodología analítica poderosa para transformar el registro histórico de mediciones ambientales en horizontes de predicción, permitiendo identificar patrones diurnos, ciclos estacionales y tendencias que orientan la toma de decisiones en materia de política pública ambiental (Box & Jenkins, 1970; Hyndman & Athanasopoulos, 2018). En particular, cuando se aplica a contaminantes atmosféricos, esta metodología proporciona insumos valiosos para anticipar episodios críticos de contaminación, emitir alertas tempranas a la población vulnerable y evaluar la efectividad de las medidas de mitigación (DAGMA, 2019).

Este trabajo analiza la serie temporal de la concentración de ozono troposférico (O₃, en µg/m³) registrada por la red de monitoreo de calidad del aire de Cali durante el año 2018, con resolución horaria. A partir de la agregación de los datos a promedios diarios, se construye y valida un modelo ARIMA siguiendo la metodología de Box-Jenkins, con el objetivo de caracterizar la dinámica temporal de la serie, evaluar la calidad del ajuste mediante pruebas diagnósticas sobre los residuales y generar pronósticos de corto plazo con sus respectivas bandas de confianza. De esta manera, el trabajo busca cerrar la brecha entre la teoría del análisis de series temporales y su aplicación práctica a un problema ambiental concreto y de alta relevancia para la ciudad.

2 Metodología

2.1 Serie de tiempo

Una serie de tiempo es una secuencia de observaciones de una variable específica registradas en intervalos de tiempo regulares y ordenadas cronológicamente (Box & Jenkins, 1970). Estas observaciones pueden estar espaciadas en forma horaria, diaria, mensual o anual, dependiendo del fenómeno bajo estudio. En el contexto ambiental, las series de tiempo permiten analizar el comportamiento histórico de variables como la concentración de contaminantes, la temperatura, la humedad o la radiación solar, con el objetivo de identificar patrones subyacentes y realizar pronósticos sobre su evolución futura (Hyndman & Athanasopoulos, 2018).

El análisis de series temporales se fundamenta en el estudio de cuatro componentes que caracterizan el comportamiento de los datos a través del tiempo. La tendencia representa el patrón de largo plazo que muestra la dirección general de la serie, ya sea creciente, decreciente o constante. La estacionalidad corresponde a fluctuaciones periódicas que se repiten en intervalos fijos de tiempo —en el caso del ozono, un marcado ciclo diurno y posibles variaciones semanales asociadas a la actividad humana—. El componente cíclico refleja oscilaciones de período no fijo, vinculadas en este caso al régimen climático bimodal del valle del Cauca (dos temporadas secas y dos lluviosas al año). Finalmente, el componente aleatorio o residual representa la variabilidad que permanece después de remover los componentes anteriores.

2.2 Modelo ARIMA

El modelo ARIMA (AutoRegressive Integrated Moving Average) es una de las herramientas estadísticas más utilizadas y versátiles para el análisis y pronóstico de series temporales (Hyndman & Athanasopoulos, 2018). Combina tres componentes fundamentales que trabajan de manera conjunta para capturar la estructura de dependencia temporal de los datos. El componente Autorregresivo (AR) establece que el valor actual de la serie puede explicarse como una función lineal de sus valores pasados. El componente Integrado (I) se refiere al proceso de diferenciación aplicado a la serie para lograr estacionariedad, removiendo tendencias. El componente Media Móvil (MA) modela la dependencia entre una observación y los errores de pronóstico de observaciones pasadas.

La notación del modelo se expresa como ARIMA(p,d,q). El componente AR(p) se describe en la ecuación 1:

\[ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \ldots + \phi_p y_{t-p} + \varepsilon_t \tag{1} \]

donde \(\varepsilon_t\) es ruido blanco. El parámetro \(d\) indica el número de diferenciaciones necesarias para alcanzar la estacionariedad. La primera diferenciación (\(d=1\)) se muestra en la ecuación 2:

\[ y'_t = y_t - y_{t-1} \tag{2} \]

El componente MA(q) se presenta en la ecuación 3:

\[ y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \ldots + \theta_q \varepsilon_{t-q} \tag{3} \]

y el modelo ARMA(\(p, q\)) general se expresa en la ecuación 4:

\[ y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q} \tag{4} \]

Cuando la serie presenta estacionalidad de período \(s\), el modelo se extiende a su versión estacional SARIMA(p,d,q)(P,D,Q)\(_s\), que incorpora términos autorregresivos y de media móvil sobre los rezagos estacionales.

2.3 Metodología de Box-Jenkins

La estimación de modelos ARIMA se realiza tradicionalmente mediante la metodología de Box-Jenkins, un procedimiento sistemático e iterativo que consta de cuatro etapas fundamentales (Box & Jenkins, 1970).

2.3.1 Primera Etapa: Identificación del Modelo

En esta fase se debe asegurar que la serie sea estacionaria, es decir, que su media, varianza y autocorrelación se mantengan constantes a través del tiempo. Para verificar la estacionariedad se aplica la prueba de Dickey-Fuller Aumentada (ADF), cuya hipótesis nula establece que la serie posee una raíz unitaria y por tanto no es estacionaria:

\(H_0\): La serie tiene raíz unitaria (no es estacionaria).
\(H_1\): La serie no tiene raíz unitaria (es estacionaria).

Si la prueba ADF arroja un p-valor inferior al nivel de significancia (típicamente 0.05), se rechaza \(H_0\) y se concluye que la serie es estacionaria. Cuando la serie no es estacionaria, se aplica diferenciación. Los órdenes \(p\) y \(q\) se identifican mediante el análisis de las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF).

2.3.2 Segunda Etapa: Estimación de Parámetros

Una vez identificados los órdenes candidatos, se estiman los coeficientes del modelo mediante máxima verosimilitud. Se comparan varias especificaciones utilizando criterios de información como el AIC, el AICc y el BIC, prefiriéndose el modelo con menor valor, en equilibrio con la parsimonia.

2.3.3 Tercera Etapa: Diagnóstico y Validación del Modelo

Se verifica que los residuales se comporten como ruido blanco: media cero, varianza constante y ausencia de autocorrelación. Para ello se aplican la prueba de Ljung-Box (autocorrelación de los residuales), la prueba de Shapiro-Wilk (normalidad) y el análisis gráfico de la ACF/PACF de los residuales.

2.3.4 Cuarta Etapa: Predicción o Pronóstico

Validado el modelo, se generan pronósticos para horizontes futuros junto con sus intervalos de confianza, y se evalúa su desempeño predictivo sobre una ventana de prueba mediante métricas como MAE, RMSE y MAPE.

2.4 Descripción de las variables y fuente de los datos

La base de datos (Compartir2.xlsx) contiene 8.760 registros horarios correspondientes al año 2018 para una estación de monitoreo de la ciudad de Cali. El cuadro siguiente describe las variables disponibles; la variable de interés para este taller es la concentración de ozono troposférico (O₃).

Variable	Descripción	Unidad
`Fecha & Hora`	Marca temporal de la medición (resolución horaria)	fecha-hora
`O3`	Concentración de ozono troposférico (variable objetivo)	µg/m³
`Vel Viento`	Velocidad del viento	m/s
`Dir Viento`	Dirección del viento	grados
`Temperatura`	Temperatura del aire	°C
`Humedad`	Humedad relativa	%
`Radiacion Solar`	Radiación solar global	Watt/m²
`Lluvia`	Precipitación acumulada	mm

3 Descripción de la serie temporal

3.1 Carga y preparación de los datos

A continuación se carga la base, se renombran las columnas a nombres manejables, se construye la marca temporal y se agrega la serie horaria a promedios diarios. La agregación diaria permite obtener una serie limpia de 365 observaciones para 2018, ideal para el modelado ARIMA, mientras que el ciclo diurno (intra-día) se analiza por separado en la sección descriptiva.

# --- Lectura del archivo Excel ---
ruta <- "Compartir2.xlsx"
crudo <- read_excel(ruta, sheet = 1)

# Renombramos las columnas a nombres seguros (en el mismo orden del archivo)
names(crudo) <- c("FechaHora", "O3", "VelViento", "DirViento",
                  "Temperatura", "Humedad", "Radiacion", "Lluvia")

# Aseguramos el tipo fecha-hora y extraemos la fecha (día) y la hora
crudo <- crudo %>%
  mutate(
    FechaHora = as.POSIXct(FechaHora),
    Fecha     = as.Date(FechaHora),
    Hora      = hour(FechaHora),
    Mes       = month(FechaHora, label = TRUE, abbr = TRUE),
    DiaSemana = wday(FechaHora, label = TRUE, abbr = TRUE, week_start = 1)
  ) %>%
  filter(Fecha >= as.Date("2018-01-01") & Fecha <= as.Date("2018-12-31"))

# --- Agregación a promedios DIARIOS (serie objetivo) ---
diario <- crudo %>%
  group_by(Fecha) %>%
  summarise(
    O3          = mean(O3,          na.rm = TRUE),
    Temperatura = mean(Temperatura, na.rm = TRUE),
    Humedad     = mean(Humedad,     na.rm = TRUE),
    Radiacion   = mean(Radiacion,   na.rm = TRUE),
    VelViento   = mean(VelViento,   na.rm = TRUE),
    Lluvia      = sum(Lluvia,       na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(Fecha)

# Algunos días pueden quedar como NaN si no hubo registros; los marcamos como NA
diario <- diario %>% mutate(across(where(is.numeric),
                                   ~ ifelse(is.nan(.x), NA, .x)))

# Construimos el objeto ts con periodicidad semanal (frecuencia 7) e
# interpolamos los pocos huecos para garantizar una serie continua.
serie_ts <- ts(diario$O3, frequency = 7)
serie_ts <- na.interp(serie_ts)          # interpolación de valores faltantes
diario$O3 <- as.numeric(serie_ts)        # serie diaria de O3 ya completa

cat("Número de días (observaciones diarias):", nrow(diario), "\n")

Número de días (observaciones diarias): 365

cat("Rango de fechas:", format(min(diario$Fecha)), "a", format(max(diario$Fecha)), "\n")

Rango de fechas: 2018-01-01 a 2018-12-31

cat("Concentración media diaria de O3:", round(mean(diario$O3), 2), "ug/m3\n")

Concentración media diaria de O3: 29.59 ug/m3

3.2 Estadísticas descriptivas

# Función auxiliar para resumir cada variable diaria
resumen_var <- function(x, etiqueta, unidad) {
  data.frame(
    Variable = etiqueta,
    Unidad   = unidad,
    Media    = round(mean(x, na.rm = TRUE), 2),
    Mediana  = round(median(x, na.rm = TRUE), 2),
    Desv     = round(sd(x, na.rm = TRUE), 2),
    Minimo   = round(min(x, na.rm = TRUE), 2),
    Maximo   = round(max(x, na.rm = TRUE), 2),
    CV       = round(sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100, 1),
    stringsAsFactors = FALSE
  )
}

tabla_desc <- bind_rows(
  resumen_var(diario$O3,          "Ozono (O3)",        "ug/m3"),
  resumen_var(diario$Temperatura, "Temperatura",       "C"),
  resumen_var(diario$Humedad,     "Humedad relativa",  "%"),
  resumen_var(diario$Radiacion,   "Radiacion solar",   "W/m2"),
  resumen_var(diario$VelViento,   "Velocidad viento",  "m/s"),
  resumen_var(diario$Lluvia,      "Lluvia (diaria)",   "mm")
)

verde <- list(profundo = "#0B6E4F", medio = "#1B998B", claro = "#E8F5F1",
              texto = "#FFFFFF")

tabla_desc %>%
  gt() %>%
  tab_header(
    title    = md("**Tabla 1. Estadísticas descriptivas de las variables (promedios diarios, 2018)**"),
    subtitle = md("Estación de monitoreo de calidad del aire — Cali")
  ) %>%
  cols_label(
    Desv = "Desv. Est.", CV = "CV (%)"
  ) %>%
  fmt_number(columns = c(Media, Mediana, Desv, Minimo, Maximo, CV), decimals = 2) %>%
  tab_style(
    style = list(cell_fill(color = verde$profundo),
                 cell_text(color = verde$texto, weight = "bold")),
    locations = cells_column_labels(everything())
  ) %>%
  tab_style(
    style = list(cell_fill(color = verde$claro), cell_text(weight = "bold")),
    locations = cells_body(columns = Variable)
  ) %>%
  tab_options(table.font.size = px(13),
              heading.title.font.size = px(16),
              heading.subtitle.font.size = px(13))

Variable	Unidad	Media	Mediana	Desv. Est.	Minimo	Maximo	CV (%)
Tabla 1. Estadísticas descriptivas de las variables (promedios diarios, 2018)
Estación de monitoreo de calidad del aire — Cali
Ozono (O3)	ug/m3	29.59	29.40	7.55	12.70	58.68	25.50
Temperatura	C	28.28	28.43	1.21	24.34	30.93	4.30
Humedad relativa	%	57.66	57.06	5.53	42.50	73.20	9.60
Radiacion solar	W/m2	108.24	104.14	38.00	11.17	253.80	35.10
Velocidad viento	m/s	2.51	2.48	0.31	1.80	3.77	12.50
Lluvia (diaria)	mm	6.59	0.00	16.37	0.00	108.96	248.50

La concentración media diaria de ozono se ubica alrededor de los 30 µg/m³, con un coeficiente de variación elevado que evidencia una alta variabilidad día a día. Como es característico de este contaminante, la temperatura y la radiación solar —motores de la formación fotoquímica del O₃— presentan promedios elevados, coherentes con las condiciones de un valle tropical con alta insolación.

3.3 Evolución temporal de la serie (gráfico dinámico)

# Media móvil de 30 días para visualizar la tendencia de fondo
diario <- diario %>%
  mutate(MM30 = zoo::rollmean(O3, k = 30, fill = NA, align = "center"))

p_serie <- ggplot(diario, aes(x = Fecha)) +
  geom_line(aes(y = O3), color = "#1B998B", linewidth = 0.5, alpha = 0.8) +
  geom_line(aes(y = MM30), color = "#0B6E4F", linewidth = 1.1, na.rm = TRUE) +
  labs(
    title = "Concentración diaria de Ozono (O3) en Cali — 2018",
    subtitle = "Serie diaria (línea clara) y media móvil de 30 días (línea oscura)",
    x = NULL, y = expression(O[3]~"(µg/m"^3*")")
  ) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))

# Versión dinámica (interactiva) con plotly
ggplotly(p_serie) %>% layout(hovermode = "x unified")

La serie muestra una clara modulación climática: las concentraciones de ozono tienden a ser más altas durante las temporadas secas y de mayor radiación, y descienden en los periodos lluviosos, cuando la nubosidad reduce la actividad fotoquímica y la lluvia favorece la remoción del contaminante. Este comportamiento concuerda con el régimen bimodal de precipitación del valle del Cauca.

3.4 Estacionalidad: ciclo diurno y patrón mensual

El ozono es un contaminante fuertemente gobernado por la luz solar; por ello su comportamiento intra-día es uno de los rasgos más característicos de la serie. El gráfico siguiente, construido a partir de los datos horarios originales, muestra el ciclo diurno promedio del O₃.

perfil_horario <- crudo %>%
  group_by(Hora) %>%
  summarise(O3_medio = mean(O3, na.rm = TRUE),
            O3_sd    = sd(O3,  na.rm = TRUE), .groups = "drop")

p_diurno <- ggplot(perfil_horario, aes(x = Hora, y = O3_medio)) +
  geom_ribbon(aes(ymin = O3_medio - O3_sd, ymax = O3_medio + O3_sd),
              fill = "#1B998B", alpha = 0.18) +
  geom_line(color = "#0B6E4F", linewidth = 1.2) +
  geom_point(color = "#0B6E4F", size = 1.8) +
  scale_x_continuous(breaks = seq(0, 23, 2)) +
  labs(title = "Ciclo diurno promedio del Ozono (O3) — Cali 2018",
       subtitle = "Promedio horario ± 1 desviación estándar",
       x = "Hora del día", y = expression(O[3]~"(µg/m"^3*")")) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))
p_diurno

p_mes <- ggplot(crudo, aes(x = Mes, y = O3)) +
  geom_boxplot(fill = "#1B998B", alpha = 0.55, outlier.alpha = 0.15,
               outlier.size = 0.6) +
  stat_summary(fun = mean, geom = "point", color = "#C0392B", size = 2) +
  labs(title = "Distribución mensual del Ozono (O3) — Cali 2018",
       subtitle = "Cajas: distribución horaria por mes · puntos rojos: media mensual",
       x = NULL, y = expression(O[3]~"(µg/m"^3*")")) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))
p_mes

El ciclo diurno confirma la naturaleza fotoquímica del ozono: las concentraciones son mínimas en la madrugada, comienzan a aumentar tras el amanecer y alcanzan su máximo hacia el mediodía y primeras horas de la tarde, cuando la radiación solar es más intensa, para luego decaer al anochecer. A escala mensual se aprecia la huella del régimen climático bimodal de Cali, con meses secos de mayor concentración y meses lluviosos con valores más bajos.

4 Resultados del modelo ARIMA

4.1 Ventanas de entrenamiento y prueba

Siguiendo la práctica estándar, dividimos la serie diaria en una ventana de entrenamiento (con la que se identifica y estima el modelo) y una ventana de prueba de las últimas 21 observaciones (3 semanas), reservada para evaluar la capacidad predictiva.

h <- 21                                  # horizonte de pronóstico (días)
n <- nrow(diario)

train_df <- diario[1:(n - h), ]
test_df  <- diario[(n - h + 1):n, ]

ventana  <- ts(train_df$O3, frequency = 7)   # entrenamiento
# La ventana de prueba continúa el índice temporal del entrenamiento, de modo
# que las fechas del pronóstico y de la prueba queden alineadas.
ventana2 <- ts(test_df$O3, frequency = 7, start = tsp(ventana)[2] + 1/7)  # prueba

cat("Observaciones de entrenamiento:", length(ventana), "\n")

Observaciones de entrenamiento: 344

cat("Observaciones de prueba:       ", length(ventana2), "\n")

Observaciones de prueba:        21

ggplot() +
  geom_line(data = train_df, aes(Fecha, O3, color = "Entrenamiento"), linewidth = 0.6) +
  geom_line(data = test_df,  aes(Fecha, O3, color = "Prueba"),        linewidth = 0.8) +
  scale_color_manual(values = c("Entrenamiento" = "#1B998B", "Prueba" = "#C0392B")) +
  labs(title = "Partición de la serie diaria de O3", x = NULL,
       y = expression(O[3]~"(µg/m"^3*")"), color = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top",
        plot.title = element_text(face = "bold", color = "#0B6E4F"))

4.2 Estacionariedad: prueba ADF y función de autocorrelación

# Prueba de Dickey-Fuller aumentada sobre la serie de entrenamiento
adf_original <- adf.test(ventana)
print(adf_original)


    Augmented Dickey-Fuller Test

data:  ventana
Dickey-Fuller = -3.1027, Lag order = 6, p-value = 0.1115
alternative hypothesis: stationary

grid.arrange(
  ggAcf(ventana)  + ggtitle("ACF — serie original")  + theme_minimal(base_size = 11),
  ggPacf(ventana) + ggtitle("PACF — serie original") + theme_minimal(base_size = 11),
  nrow = 1
)

La prueba ADF arroja un p-valor de 0.1115. Bajo el nivel de significancia de 0.05, la serie se considera no estacionaria. Por tanto, se aplica una diferenciación de primer orden para alcanzar la estacionariedad.

4.3 Diferenciación de la serie

miserie <- diff(ventana) %>% na.omit()     # primera diferencia

adf_diff <- adf.test(miserie)
print(adf_diff)


    Augmented Dickey-Fuller Test

data:  miserie
Dickey-Fuller = -9.9416, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary

grid.arrange(
  autoplot(miserie) + ggtitle("Serie diferenciada (d = 1)") +
    ylab("ΔO3") + theme_minimal(base_size = 11),
  ggAcf(miserie) + ggtitle("ACF — serie diferenciada") + theme_minimal(base_size = 11),
  nrow = 1
)

grid.arrange(
  ggAcf(miserie)  + ggtitle("ACF — ΔO3")  + theme_minimal(base_size = 11),
  ggPacf(miserie) + ggtitle("PACF — ΔO3") + theme_minimal(base_size = 11),
  nrow = 1
)

Tras la primera diferenciación, la prueba ADF confirma la estacionariedad (p-valor muy pequeño) y las gráficas ACF/PACF presentan un decaimiento rápido. El patrón de cortes en la ACF y la PADF orienta la selección de los órdenes \(p\) y \(q\) de los modelos candidatos.

4.4 Identificación y estimación de modelos candidatos

Se combina la sugerencia automática de auto.arima con la estimación manual de varias especificaciones, comparándolas mediante los criterios de información AIC, AICc y BIC.

# Modelo sugerido automáticamente (búsqueda exhaustiva, incluye estacionalidad)
modelo_auto <- auto.arima(ventana, seasonal = TRUE,
                          stepwise = FALSE, approximation = FALSE)

# Modelos candidatos estimados manualmente (metodología Box-Jenkins)
modelo1 <- Arima(ventana, order = c(1, 1, 1))
modelo2 <- Arima(ventana, order = c(2, 1, 1))
modelo3 <- Arima(ventana, order = c(1, 1, 2))
modelo4 <- Arima(ventana, order = c(2, 1, 2))
modelo5 <- Arima(ventana, order = c(0, 1, 1))

modelos <- list(Auto = modelo_auto, M1 = modelo1, M2 = modelo2,
                M3 = modelo3, M4 = modelo4, M5 = modelo5)

# Etiqueta legible del orden ARIMA de cada modelo
etiqueta_arima <- function(m) {
  ar <- arimaorder(m)
  if (length(ar) == 7) {
    sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[%d]",
            ar[1], ar[2], ar[3], ar[4], ar[5], ar[6], ar[7])
  } else {
    sprintf("ARIMA(%d,%d,%d)", ar[1], ar[2], ar[3])
  }
}

tabla_modelos <- data.frame(
  Modelo = names(modelos),
  Especificacion = sapply(modelos, etiqueta_arima),
  AIC  = round(sapply(modelos, AIC), 2),
  AICc = round(sapply(modelos, function(m) m$aicc), 2),
  BIC  = round(sapply(modelos, BIC), 2),
  row.names = NULL
)

tabla_modelos %>%
  gt() %>%
  tab_header(title = md("**Tabla 2. Comparación de modelos ARIMA candidatos**"),
             subtitle = "Criterios de información (menor es mejor)") %>%
  fmt_number(columns = c(AIC, AICc, BIC), decimals = 2) %>%
  tab_style(style = list(cell_fill(color = verde$profundo),
                         cell_text(color = verde$texto, weight = "bold")),
            locations = cells_column_labels(everything())) %>%
  tab_style(
    style = list(cell_fill(color = "#E8F5F1"), cell_text(weight = "bold")),
    locations = cells_body(rows = which.min(tabla_modelos$AICc))
  ) %>%
  tab_options(table.font.size = px(13), heading.title.font.size = px(16))

Modelo	Especificacion	AIC	AICc	BIC
Tabla 2. Comparación de modelos ARIMA candidatos
Criterios de información (menor es mejor)
Auto	ARIMA(1,1,2)	2,154.30	2,154.42	2,169.65
M1	ARIMA(1,1,1)	2,156.43	2,156.50	2,167.94
M2	ARIMA(2,1,1)	2,155.10	2,155.21	2,170.45
M3	ARIMA(1,1,2)	2,154.30	2,154.42	2,169.65
M4	ARIMA(2,1,2)	2,157.36	2,157.54	2,176.55
M5	ARIMA(0,1,1)	2,156.33	2,156.37	2,164.01

# Seleccionamos el modelo con menor AICc
idx_mejor <- which.min(sapply(modelos, function(m) m$aicc))
modelo_sel <- modelos[[idx_mejor]]
nombre_sel <- etiqueta_arima(modelo_sel)
cat("Modelo seleccionado (menor AICc):", names(modelos)[idx_mejor],
    "->", nombre_sel, "\n")

Modelo seleccionado (menor AICc): Auto -> ARIMA(1,1,2)

4.5 Interpretación de los coeficientes

coefs <- coef(modelo_sel)
if (length(coefs) > 0) {
  se   <- sqrt(diag(modelo_sel$var.coef))
  zval <- coefs / se
  pval <- 2 * (1 - pnorm(abs(zval)))
  tabla_coef <- data.frame(
    Termino     = names(coefs),
    Coeficiente = round(coefs, 4),
    Error_Est   = round(se, 4),
    z           = round(zval, 3),
    p_valor     = ifelse(pval < 0.001, "< 0.001", sprintf("%.4f", pval)),
    row.names = NULL
  )

  tabla_coef %>%
    gt() %>%
    tab_header(title = md(paste0("**Tabla 3. Coeficientes del modelo seleccionado — ",
                                 nombre_sel, "**"))) %>%
    cols_label(Error_Est = "Error Est.", p_valor = "p-valor") %>%
    tab_style(style = list(cell_fill(color = verde$profundo),
                           cell_text(color = verde$texto, weight = "bold")),
              locations = cells_column_labels(everything())) %>%
    tab_options(table.font.size = px(13), heading.title.font.size = px(15))
} else {
  cat("El modelo seleccionado no contiene coeficientes estimados",
      "(p. ej. un paseo aleatorio ARIMA(0,1,0)).\n")
}

Termino	Coeficiente	Error Est.	z	p-valor
Tabla 3. Coeficientes del modelo seleccionado — ARIMA(1,1,2)
ar1	0.8027	0.1464	5.483	< 0.001
ma1	-1.4296	0.1685	-8.482	< 0.001
ma2	0.4658	0.1360	3.425	< 0.001

Cada coeficiente autorregresivo (\(\phi\)) cuantifica la persistencia de la serie: indica en qué medida el nivel de ozono de un día depende de los días previos. Los coeficientes de media móvil (\(\theta\)) capturan el efecto de los choques (errores) recientes sobre el valor actual. Un p-valor inferior a 0.05 indica que el término es estadísticamente significativo y aporta a la explicación de la dinámica de la serie.

4.6 Desempeño predictivo sobre la ventana de prueba

# Pronóstico del modelo de entrenamiento sobre el horizonte de prueba
fc_test <- forecast(modelo_sel, h = h, level = c(80, 95))

# Métricas de exactitud comparando contra la ventana de prueba
acc <- accuracy(fc_test, ventana2)
acc_df <- as.data.frame(round(acc, 3))
acc_df <- cbind(Conjunto = rownames(acc_df), acc_df)
rownames(acc_df) <- NULL

acc_df %>%
  gt() %>%
  tab_header(title = md("**Tabla 4. Métricas de exactitud del modelo seleccionado**"),
             subtitle = "Entrenamiento vs. prueba (últimos 21 días)") %>%
  tab_style(style = list(cell_fill(color = verde$profundo),
                         cell_text(color = verde$texto, weight = "bold")),
            locations = cells_column_labels(everything())) %>%
  tab_options(table.font.size = px(12), heading.title.font.size = px(15))

Conjunto	ME	RMSE	MAE	MPE	MAPE	MASE	ACF1	Theil's U
Tabla 4. Métricas de exactitud del modelo seleccionado
Entrenamiento vs. prueba (últimos 21 días)
Training set	0.05	5.513	4.258	-3.618	15.995	0.656	-0.001	NA
Test set	2.95	6.598	5.488	5.519	19.307	0.846	0.466	0.825

El MAPE (error porcentual absoluto medio) sobre la ventana de prueba resume el error típico del pronóstico en términos relativos; el RMSE y el MAE lo expresan en las unidades originales (µg/m³). Que las métricas de entrenamiento y prueba sean de magnitud comparable indica que el modelo no está sobreajustado y generaliza razonablemente a datos no vistos.

5 Pronóstico

Una vez validado el modelo, se reestima sobre la serie completa (los 365 días de 2018) y se proyectan las 21 observaciones futuras (primeras tres semanas de 2019), con bandas de confianza al 80 % y 95 %.

# Reestimamos el modelo seleccionado sobre TODA la serie
serie_full <- ts(diario$O3, frequency = 7)
orden_sel  <- arimaorder(modelo_sel)

if (length(orden_sel) == 7) {
  modelo_final <- Arima(serie_full,
                        order    = orden_sel[1:3],
                        seasonal = list(order = orden_sel[4:6],
                                        period = orden_sel[7]))
} else {
  modelo_final <- Arima(serie_full, order = orden_sel[1:3])
}

fc_fut <- forecast(modelo_final, h = h, level = c(80, 95))

# Fechas futuras asociadas al pronóstico
fechas_fut <- max(diario$Fecha) + 1:h

tabla_fc <- data.frame(
  Fecha       = fechas_fut,
  Pronostico  = round(as.numeric(fc_fut$mean), 2),
  Lo80        = round(as.numeric(fc_fut$lower[, 1]), 2),
  Hi80        = round(as.numeric(fc_fut$upper[, 1]), 2),
  Lo95        = round(as.numeric(fc_fut$lower[, 2]), 2),
  Hi95        = round(as.numeric(fc_fut$upper[, 2]), 2)
)

datatable(tabla_fc,
          rownames = FALSE,
          colnames = c("Fecha", "Pronóstico (µg/m³)", "Inf. 80%", "Sup. 80%",
                       "Inf. 95%", "Sup. 95%"),
          caption = "Tabla 5. Pronóstico de O3 a 21 días con bandas de confianza",
          options = list(pageLength = 10, dom = "tip")) %>%
  formatStyle(columns = "Pronostico", fontWeight = "bold")

autoplot(fc_fut, include = 90) +
  labs(title = paste0("Pronóstico de O3 a 21 días — ", nombre_sel),
       subtitle = "Últimos 90 días observados + pronóstico con bandas al 80% y 95%",
       x = "Tiempo (semanas)", y = expression(O[3]~"(µg/m"^3*")")) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#0B6E4F"))

6 Pruebas sobre los residuales

La validación del modelo (tercera etapa de Box-Jenkins) exige comprobar que los residuales se comporten como ruido blanco.

# Diagnóstico gráfico integral + prueba de Ljung-Box
checkresiduals(modelo_final)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,2)
Q* = 8.0076, df = 11, p-value = 0.7126

Model df: 3.   Total lags used: 14

res <- residuals(modelo_final)

# Prueba de Ljung-Box (autocorrelación de los residuales)
lb <- Box.test(res, lag = 14, fitdf = length(coef(modelo_final)),
               type = "Ljung-Box")

# Prueba de normalidad de Shapiro-Wilk
sw <- shapiro.test(as.numeric(res))

tabla_tests <- data.frame(
  Prueba = c("Ljung-Box (lag = 14)", "Shapiro-Wilk"),
  Estadistico = round(c(lb$statistic, sw$statistic), 4),
  p_valor = c(ifelse(lb$p.value < 0.001, "< 0.001", sprintf("%.4f", lb$p.value)),
              ifelse(sw$p.value < 0.001, "< 0.001", sprintf("%.4f", sw$p.value))),
  Hipotesis = c("H0: residuales no autocorrelacionados (ruido blanco)",
                "H0: residuales con distribución normal"),
  Conclusion = c(ifelse(lb$p.value > 0.05,
                        "No se rechaza H0: sin autocorrelación",
                        "Se rechaza H0: autocorrelación remanente"),
                 ifelse(sw$p.value > 0.05,
                        "No se rechaza H0: normalidad",
                        "Se rechaza H0: no normalidad")),
  row.names = NULL
)

tabla_tests %>%
  gt() %>%
  tab_header(title = md("**Tabla 6. Pruebas diagnósticas sobre los residuales**")) %>%
  cols_label(p_valor = "p-valor") %>%
  tab_style(style = list(cell_fill(color = verde$profundo),
                         cell_text(color = verde$texto, weight = "bold")),
            locations = cells_column_labels(everything())) %>%
  tab_options(table.font.size = px(12), heading.title.font.size = px(15))

Prueba	Estadistico	p-valor	Hipotesis	Conclusion
Tabla 6. Pruebas diagnósticas sobre los residuales
Ljung-Box (lag = 14)	8.0076	0.7126	H0: residuales no autocorrelacionados (ruido blanco)	No se rechaza H0: sin autocorrelación
Shapiro-Wilk	0.9945	0.2111	H0: residuales con distribución normal	No se rechaza H0: normalidad

cat("Media de los residuales      :", round(mean(res, na.rm = TRUE), 4), "\n")

Media de los residuales      : 0.1194

cat("Desviación estándar          :", round(sd(res, na.rm = TRUE), 4), "\n")

Desviación estándar          : 5.5236

cat("Asimetría (ideal = 0)        :", round(skewness(res, na.rm = TRUE), 4), "\n")

Asimetría (ideal = 0)        : -0.1469

cat("Curtosis  (ideal = 3)        :", round(kurtosis(res, na.rm = TRUE), 4), "\n")

Curtosis  (ideal = 3)        : 3.5279

grid.arrange(
  ggAcf(res)  + ggtitle("ACF de los residuales")  + theme_minimal(base_size = 11),
  ggPacf(res) + ggtitle("PACF de los residuales") + theme_minimal(base_size = 11),
  nrow = 1
)

La prueba de Ljung-Box arroja un p-valor de 0.7126: no se rechaza la hipótesis de independencia, lo que indica que el modelo ha capturado adecuadamente la estructura de autocorrelación de la serie. La prueba de Shapiro-Wilk (p-valor 0.2111) no permite rechazar la hipótesis de normalidad; como es habitual en series ambientales, pueden persistir colas algo más pesadas que las de una distribución normal, sin que ello invalide la utilidad del modelo para el pronóstico puntual.

7 Gráfico profesional integrado

El siguiente gráfico, listo para incluir en el informe, resume en una sola figura los elementos esenciales del análisis: la serie original, su tendencia de fondo, la componente estacional (extraída por descomposición STL) y el pronóstico con sus bandas de confianza. La figura se guarda automáticamente como grafico_profesional_O3_Cali.png (300 dpi).

# --- Descomposición STL para extraer tendencia y estacionalidad ---
stl_fit <- stl(serie_full, s.window = "periodic")
comp <- as.data.frame(stl_fit$time.series)
diario$Tendencia <- as.numeric(comp$trend)
diario$Estacional <- as.numeric(comp$seasonal)

# --- Data frame del pronóstico con fechas reales ---
df_fc <- data.frame(
  Fecha = fechas_fut,
  mean  = as.numeric(fc_fut$mean),
  lo80  = as.numeric(fc_fut$lower[, 1]), hi80 = as.numeric(fc_fut$upper[, 1]),
  lo95  = as.numeric(fc_fut$lower[, 2]), hi95 = as.numeric(fc_fut$upper[, 2])
)

col_serie <- "#1B998B"; col_trend <- "#0B6E4F"
col_fc <- "#C0392B";    col_band95 <- "#E59866"; col_band80 <- "#D98880"

# Panel superior: serie + tendencia + pronóstico con bandas
p_main <- ggplot() +
  geom_line(data = diario, aes(Fecha, O3, color = "Serie observada"),
            linewidth = 0.45, alpha = 0.85) +
  geom_line(data = diario, aes(Fecha, Tendencia, color = "Tendencia (STL)"),
            linewidth = 1.2) +
  geom_ribbon(data = df_fc, aes(Fecha, ymin = lo95, ymax = hi95),
              fill = col_band95, alpha = 0.35) +
  geom_ribbon(data = df_fc, aes(Fecha, ymin = lo80, ymax = hi80),
              fill = col_band80, alpha = 0.45) +
  geom_line(data = df_fc, aes(Fecha, mean, color = "Pronóstico"),
            linewidth = 1.1) +
  geom_vline(xintercept = as.numeric(max(diario$Fecha)),
             linetype = "dashed", color = "grey40") +
  scale_color_manual(values = c("Serie observada" = col_serie,
                                "Tendencia (STL)" = col_trend,
                                "Pronóstico" = col_fc)) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b %Y") +
  labs(title = "Ozono troposférico (O3) en Cali: serie, tendencia y pronóstico",
       subtitle = paste0("Modelo ", nombre_sel,
                         " · bandas de confianza al 80% y 95% · pronóstico a 21 días"),
       x = NULL, y = expression(O[3]~"(µg/m"^3*")"), color = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top",
        plot.title = element_text(face = "bold", color = col_trend, size = 14),
        panel.grid.minor = element_blank())

# Panel inferior: componente estacional (un ciclo semanal representativo)
p_season <- ggplot(diario[1:35, ], aes(Fecha, Estacional)) +
  geom_hline(yintercept = 0, color = "grey60") +
  geom_line(color = col_trend, linewidth = 0.9) +
  geom_area(fill = col_serie, alpha = 0.20) +
  labs(title = "Componente estacional (semanal) extraída por STL",
       x = NULL, y = "Efecto estacional") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", color = col_trend, size = 11))

# Combinamos ambos paneles en una sola figura
g_final <- arrangeGrob(p_main, p_season, ncol = 1, heights = c(2.2, 1))

# Guardamos la figura para incluirla en el informe
ggsave("grafico_profesional_O3_Cali.png", g_final,
       width = 11, height = 8, dpi = 300)

grid::grid.newpage()
grid::grid.draw(g_final)

Esta figura integra de manera profesional los cuatro elementos solicitados: la serie observada de ozono (verde claro), su tendencia de largo plazo estimada por descomposición STL (verde oscuro), la componente estacional semanal (panel inferior) y el pronóstico de las próximas tres semanas (rojo) con sus bandas de confianza al 80 % y 95 %. La línea vertical punteada marca el límite entre los datos observados y el horizonte pronosticado.

8 Conclusiones

El análisis de la serie temporal de ozono troposférico (O₃) en Cali durante el año 2018, desarrollado mediante la metodología de Box-Jenkins y modelos ARIMA, permite extraer varias conclusiones relevantes tanto en el plano ambiental como en el metodológico.

En primer lugar, la serie de ozono exhibe una estructura temporal rica y multiescalar. A escala intra-día se observa un marcado ciclo diurno de naturaleza fotoquímica, con concentraciones mínimas durante la madrugada y máximas hacia el mediodía, coherente con el papel central de la radiación solar en la formación del contaminante. A escala anual, la serie de promedios diarios refleja la modulación climática bimodal del valle del Cauca, con concentraciones más altas en los periodos secos y de mayor insolación, y descensos en las temporadas lluviosas. Esta caracterización descriptiva, por sí sola, ya constituye un insumo valioso para la gestión de la calidad del aire en la ciudad.

En segundo lugar, desde la perspectiva del modelado, la serie diaria requirió una diferenciación de primer orden para alcanzar la estacionariedad, confirmada por la prueba de Dickey-Fuller aumentada. La comparación sistemática de modelos candidatos mediante los criterios AIC, AICc y BIC, complementada con la sugerencia de auto.arima, permitió seleccionar una especificación parsimoniosa que equilibra ajuste y simplicidad. El análisis de residuales —mediante las pruebas de Ljung-Box y Shapiro-Wilk y la inspección de la ACF/PACF— validó que el modelo captura adecuadamente la estructura de dependencia temporal de la serie, comportándose los residuales aproximadamente como ruido blanco en términos de autocorrelación.

En tercer lugar, el ejercicio de pronóstico demostró una capacidad predictiva razonable sobre la ventana de prueba, con métricas de error en entrenamiento y prueba de magnitud comparable, lo que sugiere ausencia de sobreajuste. Los pronósticos a 21 días, acompañados de sus bandas de confianza, proporcionan un horizonte de anticipación útil para la planeación de medidas preventivas, aunque deben interpretarse con cautela: como ocurre con la mayoría de las variables ambientales, la presencia de colas algo más pesadas que las de una distribución normal advierte sobre la posibilidad de episodios extremos que un modelo lineal univariado no logra anticipar plenamente.

Finalmente, este trabajo reconoce limitaciones que abren líneas de mejora. El modelo ARIMA es univariado e ignora variables exógenas que son determinantes físicos del ozono —la radiación solar, la temperatura, la humedad y el viento—, por lo que una extensión natural sería un modelo ARIMAX/Regresión dinámica que incorpore estos predictores meteorológicos. Asimismo, dado que la serie original es horaria con una fuerte estacionalidad diaria, futuros estudios podrían explorar modelos con estacionalidad múltiple (por ejemplo, tbats o modelos de estacionalidad horaria y semanal simultánea) para aprovechar plenamente la resolución temporal de los datos. En conjunto, el trabajo ilustra cómo las técnicas de series de tiempo transforman un registro ambiental en conocimiento accionable para la gestión de la calidad del aire en Cali.

9 Bibliografía

Box, G. E. P., & Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. Holden-Day.

DAGMA. (2019). Informe de calidad del aire del municipio de Santiago de Cali. Departamento Administrativo de Gestión del Medio Ambiente, Alcaldía de Santiago de Cali.

EPA. (2023). Ground-level Ozone Basics. United States Environmental Protection Agency. https://www.epa.gov/ground-level-ozone-pollution

Hyndman, R. J., & Athanasopoulos, G. (2018). Forecasting: Principles and Practice (2.ª ed.). OTexts. https://otexts.com/fpp2/

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1-22.

Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303.

Organización Mundial de la Salud (OMS). (2021). Directrices mundiales de la OMS sobre la calidad del aire: partículas (PM2.5 y PM10), ozono, dióxido de nitrógeno, dióxido de azufre y monóxido de carbono. OMS.

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

Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3/4), 591-611.

Análisis y pronóstico de la concentración de Ozono troposférico (O₃) en Cali mediante modelos de series de tiempo: metodología Box-Jenkins y ARIMA

Taller de Series de Tiempo — Gestión de Datos

Universidad del Valle (Cali) · Profesor: Orlando Joaqui Barandica

12 junio 2026

1 Introducción