Santiago Dominguez - 2416770 Marianne Rojas Palacios - 2417218 Daniela Alexandra Carlosama - 2515804
library(readxl)
library(dplyr)
library(lubridate)
library(xts)
library(ggplot2)
library(fpp2)
library(kableExtra)
library(moments)
library(forecast)
library(zoo)
library(gridExtra)
library(plotly)
library(tseries)La contaminación atmosférica constituye uno de los principales problemas ambientales y de salud pública en las ciudades, debido a la exposición continua de la población a contaminantes capaces de afectar tanto la salud humana como el funcionamiento de los ecosistemas urbanos.
Entre los contaminantes de mayor interés científico y regulatorio se encuentra el ozono troposférico (O₃), un contaminante secundario que no es emitido directamente a la atmósfera, sino que se forma mediante reacciones fotoquímicas entre óxidos de nitrógeno (NOₓ) y compuestos orgánicos volátiles (COV) en presencia de radiación solar.
La exposición a concentraciones elevadas de O₃ se ha asociado en la literatura científica con efectos adversos sobre el sistema respiratorio y cardiovascular, lo que hace de su monitoreo continuo una fuente de información relevante para apoyar la gestión de la calidad del aire.
En la ciudad de Cali, este seguimiento adquiere especial pertinencia dado su régimen climático bimodal, caracterizado por dos temporadas secas y dos lluviosas a lo largo del año, con variaciones en radiación solar, temperatura y precipitación que pueden favorecer cambios temporales en la concentración de O₃. Estas condiciones hacen pertinente estudiar su comportamiento mediante herramientas estadísticas que permitan describir su evolución e identificar posibles patrones.
En este contexto, el análisis de series de tiempo ofrece una metodología adecuada para modelar la dinámica temporal de variables ambientales y generar pronósticos cuantitativos. La pregunta central que orienta este trabajo es determinar si la concentración diaria de O₃ registrada en Cali durante 2018 puede describirse adecuadamente mediante un modelo ARIMA y si dicho modelo ofrece capacidad predictiva útil.
El objetivo de este trabajo es modelar y pronosticar la concentración diaria de O₃ en Cali mediante un modelo ARIMA ajustado sobre los registros de una estación de monitoreo durante el año 2018. Para ello, el estudio comienza con una caracterización descriptiva de la serie temporal de concentraciones diarias de O₃. Posteriormente se ajusta un modelo ARIMA siguiendo la metodología Box-Jenkins, evaluando tanto el cumplimiento de sus supuestos como la calidad del ajuste mediante el análisis de los residuales. Finalmente, el modelo seleccionado se emplea para generar pronósticos sobre los últimos doce días de la serie y evaluar su capacidad predictiva.
Debido a que el estudio utiliza únicamente los registros correspondientes al año 2018 de una sola estación de monitoreo, los resultados deben interpretarse dentro de ese alcance temporal y espacial, sin pretender generalizar el comportamiento del ozono para otros años o para toda la ciudad.
Los datos utilizados en este análisis provienen de una estación de monitoreo de calidad del aire ubicada en la ciudad de Cali, Colombia, y comprenden registros horarios correspondientes a las 24 horas de cada uno de los 365 días de 2018, para un total de 8.760 observaciones. La base de datos contiene mediciones de ocho variables ambientales y meteorológicas: concentración de ozono (O₃), velocidad del viento, dirección del viento, temperatura, humedad relativa, radiación solar y precipitación. Para este análisis se trabajó exclusivamente con la variable concentración de ozono troposférico (O₃), expresada en microgramos por metro cúbico (µg/m³).
Con el fin de obtener una serie adecuada para el modelado, los datos horarios fueron sometidos a las siguientes transformaciones:
Agregación diaria: las 24 lecturas horarias de cada día fueron resumidas mediante su promedio aritmético, obteniendo una serie de 365 observaciones diarias. Esta decisión permite capturar el nivel general de exposición a O₃ durante cada jornada, evitando el sesgo que introduciría seleccionar una hora específica del día — como el mediodía, donde la concentración tiende a ser máxima por efecto de la radiación solar.
Tratamiento de valores faltantes: tras construir la serie diaria, se identificaron 20 observaciones con valores faltantes en la concentración promedio de O₃. Estos fueron tratados mediante interpolación lineal, que estima cada valor ausente trazando una línea recta entre el último dato válido previo y el primero posterior al hueco. La proporción de valores imputados representa aproximadamente el 5.5% de la serie diaria (20 de 365 observaciones), por lo que la interpolación lineal constituye una alternativa razonable para preservar la continuidad temporal sin alterar sustancialmente la dinámica general.
Construcción de los objetos de serie de tiempo:
la serie diaria fue organizada como una serie de tiempo para el análisis
de autocorrelación y el ajuste del modelo ARIMA. Para la descomposición
STL se construyó adicionalmente un objeto ts con frecuencia
semanal (frequency = 7), permitiendo identificar un posible
patrón dentro de cada semana.
El modelo empleado para el análisis y pronóstico de la serie es el ARIMA(p, d, q) — siglas en inglés de AutoRegressive Integrated Moving Average —, un enfoque estadístico ampliamente utilizado para modelar series de tiempo univariadas. El modelo combina tres componentes:
El componente autorregresivo AR(p) establece que el valor actual de la serie puede modelarse como una combinación lineal de sus p valores pasados. Captura la inercia o memoria de la serie: el hecho de que lo que ocurrió ayer influye sobre lo que ocurre hoy.
El componente integrado I(d) indica el número de diferenciaciones necesarias para transformar la serie en estacionaria, es decir, para que su nivel promedio y su variabilidad se mantengan constantes en el tiempo. En el enfoque ARIMA clásico, una serie no estacionaria requiere ser transformada mediante diferenciación antes del ajuste del modelo.
El componente de medias móviles MA(q) expresa el valor actual como función lineal de los q errores de predicción pasados. Captura las correcciones que el modelo hace sobre sus propios errores recientes.
La identificación del modelo siguió la metodología Box-Jenkins, que consta de las siguientes etapas: verificación de estacionariedad mediante la prueba Dickey-Fuller Aumentada (ADF); diferenciación cuando fue necesaria; identificación preliminar de los órdenes p y q mediante los gráficos ACF y PACF, cuyos resultados sirvieron para proponer varios modelos candidatos; selección mediante el criterio de información AICc, que penaliza la complejidad para evitar el sobreajuste; y finalmente diagnóstico residual mediante las pruebas de Ljung-Box (independencia) y Shapiro-Wilk (normalidad).
Se reservó el último tramo de la serie como conjunto de prueba para simular un escenario de pronóstico real, en el que únicamente se dispone de observaciones históricas para estimar el modelo y predecir valores futuros:
Esta división permite calcular métricas de evaluación fuera de la muestra — RMSE, MAE, MAPE y MASE — que ofrecen una medida más objetiva del desempeño predictivo que las métricas calculadas sobre los datos de entrenamiento.
loess_fit <- loess(o3 ~ as.numeric(fecha), data = datos_o3, span = 0.3)
datos_o3$loess_pred <- predict(loess_fit)
plot_ly(datos_o3, x = ~fecha) %>%
add_lines(y = ~o3, name = "O₃ diario",
line = list(color = "#1A6FA8", width = 1.5)) %>%
add_lines(y = ~loess_pred, name = "Tendencia LOESS",
line = list(color = "#E07B39", width = 2.5)) %>%
layout(
title = list(text = "Concentración diaria de O₃ — Cali (2018)",
font = list(color = "#1A6FA8", size = 15)),
xaxis = list(title = "Mes",
tickformat = "%b",
dtick = "M1"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "x unified",
legend = list(orientation = "h", x = 0, y = -0.2)
)El gráfico presenta la concentración diaria de O₃ en Cali durante el año 2018 mediante dos capas complementarias: la línea azul, que refleja los valores diarios con toda su variabilidad, y la línea naranja (suavizado loess), que captura la tendencia de fondo suavizando la variabilidad diaria. La serie inicia el año en torno a 30 µg/m³, desciende levemente entre marzo y mayo hasta alcanzar su punto mínimo aproximado de 26–27 µg/m³, periodo asociado a mayores precipitaciones y menor disponibilidad de radiación solar para la formación fotoquímica del ozono. A partir de junio se observa un ascenso sostenido en la curva de tendencia que se mantiene hasta finales de agosto, alcanzando su nivel más alto del año en torno a 37–38 µg/m³, lo cual es coherente con la temporada seca principal del Valle del Cauca (junio–agosto), periodo en que las precipitaciones disminuyen y la radiación solar se intensifica. En ese mismo periodo se registra el valor diario máximo absoluto de la serie: 58.7 µg/m³ el 27 de agosto, un pico puntual que podría estar asociado a condiciones meteorológicas particulares de ese día, acumulación temporal de contaminantes o episodios locales de emisión — factores que no es posible determinar con certeza a partir de un único año de datos. Desde septiembre la concentración desciende de forma progresiva hasta cerrar el año en su nivel más bajo, en torno a 23 µg/m³. Este comportamiento sugiere un patrón asociado al régimen climático bimodal de la región, caracterizado por menores concentraciones durante los periodos lluviosos y mayores durante la temporada seca. Sin embargo, la presencia de variaciones intraanuales no implica necesariamente una estacionalidad estadísticamente demostrada, ya que el análisis se limita a un único año de observaciones, por lo que esta interpretación debe considerarse exploratoria.
plot_ly(datos_o3, x = ~mes, y = ~o3, type = "box",
fillcolor = "rgba(26,111,168,0.7)",
line = list(color = "#1A6FA8"),
marker = list(color = "#E07B39", size = 6)) %>%
layout(
title = list(text = "Distribución mensual de O₃ — Cali (2018)",
font = list(color = "#1A6FA8", size = 15)),
xaxis = list(title = "Mes"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "closest"
)El análisis de la distribución mensual del O₃ mediante diagramas de caja sugiere un comportamiento estacional con diferencias sistemáticas entre meses del año. Los meses de julio, agosto y septiembre presentan las medianas más elevadas —33.2, 36.7 y 38.7 µg/m³ respectivamente—, mientras que marzo también presenta una mediana relativamente elevada (32.1 µg/m³). En el extremo opuesto, mayo, noviembre y diciembre registran las medianas más bajas, por debajo de los 24 µg/m³, coherentes con periodos de mayor nubosidad y precipitación en la región.
En cuanto a los valores atípicos, se identificaron cinco outliers en la serie, que conviene distinguir por su magnitud. El único outlier verdaderamente extremo corresponde al 27 de agosto con 58.7 µg/m³, superando en más de 7 µg/m³ el límite superior de su mes (51.7 µg/m³) y constituyendo el valor máximo absoluto de toda la serie. En contraste, los outliers de abril (46.2 µg/m³), junio (43.3 µg/m³) y julio (17.6 µg/m³) son casos marginales, ya que superan sus límites por menos de 1 µg/m³, lo que sugiere que se trata de observaciones limítrofes más que de eventos verdaderamente anómalos. Finalmente, el 13 de mayo registró 12.7 µg/m³, por debajo del límite inferior de ese mes (14.5 µg/m³), siendo el mínimo absoluto de la serie. Estos eventos puntuales podrían estar asociados a condiciones meteorológicas particulares o episodios locales de emisión, aunque con un único año de datos no es posible establecer una causa definitiva.
resumen <- data.frame(
Estadístico = c("Media", "Mediana", "Desv. Estándar",
"Mínimo", "Máximo", "Asimetría", "Curtosis"),
Valor = c(29.58, 29.40, 7.49, 12.70, 58.68, 0.20, 2.95)
)
kable(resumen,
caption = "Estadísticas descriptivas — O₃ diario Cali 2018",
align = c("l", "r")) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
position = "center") %>%
row_spec(0, bold = TRUE, color = "white", background = "#1A6FA8") %>%
row_spec(seq(2, nrow(resumen), 2), background = "#EAF4FB") %>%
column_spec(1, bold = TRUE, color = "#2D2D2D") %>%
column_spec(2, color = "#E07B39", bold = TRUE)| Estadístico | Valor |
|---|---|
| Media | 29.58 |
| Mediana | 29.40 |
| Desv. Estándar | 7.49 |
| Mínimo | 12.70 |
| Máximo | 58.68 |
| Asimetría | 0.20 |
| Curtosis | 2.95 |
Los patrones observados tanto en la evolución temporal como en la distribución mensual quedan respaldados por los estadísticos globales de la serie. La media de 29.58 µg/m³ y la mediana de 29.40 µg/m³ son prácticamente idénticas, lo que sugiere que el comportamiento central de la serie no está significativamente afectado por los valores extremos identificados previamente. La desviación estándar de 7.49 µg/m³ refleja una variabilidad moderada a lo largo del año, mientras que el rango entre el mínimo de 12.70 µg/m³ —registrado en mayo— y el máximo de 58.68 µg/m³ —registrado en agosto— recoge la amplitud total de las fluctuaciones observadas. Finalmente, la asimetría de 0.20 y la curtosis de 2.95 indican una distribución prácticamente simétrica y con un comportamiento de colas muy cercano al de una distribución normal, proporcionando una caracterización general consistente con lo observado en los análisis gráficos previos.
# Calcular densidad manualmente
dens <- density(datos_o3$o3)
dens_df <- data.frame(x = dens$x, y = dens$y)
# Escalar la densidad a la frecuencia del histograma
n <- nrow(datos_o3)
ancho <- (max(datos_o3$o3) - min(datos_o3$o3)) / 30
escala <- n * ancho
plot_ly() %>%
add_histogram(data = datos_o3, x = ~o3,
nbinsx = 30,
marker = list(color = "rgba(26,111,168,0.7)",
line = list(color = "#1A6FA8", width = 1)),
name = "Frecuencia") %>%
add_lines(data = dens_df, x = ~x, y = ~(y * escala),
line = list(color = "#8B0000", width = 2.5),
name = "Densidad") %>%
layout(
title = list(text = "Distribución de O₃ — Cali (2018)",
font = list(color = "#1A6FA8", size = 15)),
xaxis = list(title = "O₃ (µg/m³)"),
yaxis = list(title = "Frecuencia"),
hovermode = "x unified",
legend = list(orientation = "h", x = 0, y = -0.2)
)La distribución de la concentración diaria de O₃ en Cali durante 2018 presenta una forma aproximadamente simétrica y unimodal, con la mayor concentración de observaciones ubicada en el rango de 25 a 30 µg/m³, coherente con una media de 29.58 µg/m³ y una mediana de 29.40 µg/m³. La cercanía entre estos dos estadísticos confirma que la distribución no está fuertemente desplazada hacia ningún extremo, lo cual es una señal de estabilidad en el comportamiento central de la serie.
La desviación estándar de 7.49 µg/m³ indica una dispersión moderada alrededor de la media, equivalente a aproximadamente el 25% de su valor, lo que refleja una variabilidad moderada en las concentraciones observadas a lo largo del año. El rango total de la serie va desde 12.7 µg/m³ —registrado el 13 de mayo— hasta 58.68 µg/m³ —el 27 de agosto— evidenciando una diferencia de casi 46 µg/m³ entre los extremos, consistente con lo observado en los análisis previos.
El coeficiente de asimetría de 0.20 es consistente con lo observado en el histograma: una leve cola hacia la derecha explicada por los valores altos observados entre julio y septiembre. Al ser positivo pero cercano a cero, indica que la distribución es prácticamente simétrica, sin que los valores extremos superiores distorsionen significativamente su forma general. Finalmente, la curtosis de 2.95 es muy cercana al valor teórico de una distribución normal (3), lo que indica que la distribución observada no presenta colas excesivamente pesadas ni concentraciones extremas de probabilidad respecto a dicho referente.
acf_vals <- acf(serie, lag.max = 30, plot = FALSE)
acf_df <- data.frame(
lag = as.numeric(acf_vals$lag[-1]),
acf = as.numeric(acf_vals$acf[-1])
)
ic <- qnorm(0.975) / sqrt(length(serie))
plot_ly(acf_df, x = ~lag, y = ~acf, type = "bar",
marker = list(color = "rgba(26,111,168,0.8)",
line = list(color = "#1A6FA8", width = 1)),
name = "ACF") %>%
add_lines(x = c(0, 31), y = c(ic, ic),
line = list(color = "#E07B39", dash = "dash", width = 1.5),
name = "Banda superior", showlegend = FALSE) %>%
add_lines(x = c(0, 31), y = c(-ic, -ic),
line = list(color = "#E07B39", dash = "dash", width = 1.5),
name = "Banda inferior", showlegend = FALSE) %>%
layout(
title = list(text = "Función de Autocorrelación — O₃ diario Cali (2018)",
font = list(color = "#1A6FA8", size = 15)),
xaxis = list(title = "Lag", dtick = 5),
yaxis = list(title = "ACF", range = c(-0.15, 0.7)),
hovermode = "x unified",
bargap = 0.3
)La función de autocorrelación de la serie original de O₃ muestra un patrón de decaimiento lento y persistente, donde todas las barras permanecen muy por encima de las bandas de confianza durante los 30 rezagos analizados. Este comportamiento es consistente con una posible no estacionariedad de la serie, entendida como la presencia de un nivel promedio que no se mantiene constante a lo largo del tiempo. Esta es una evidencia visual que será contrastada formalmente más adelante mediante la prueba ADF.
Los primeros rezagos presentan las autocorrelaciones más altas: lag 1 = 0.621, lag 2 = 0.567, lag 3 = 0.506, lag 4 = 0.493, descendiendo de forma gradual hasta lag 7 = 0.436. Esto indica que cada observación diaria de O₃ está fuertemente relacionada con sus valores recientes, lo cual tiene sentido físicamente: las condiciones atmosféricas que determinan la concentración de ozono en un día dado —temperatura, radiación solar, viento— no cambian abruptamente de un día para otro, sino que evolucionan de forma continua, haciendo que el O₃ conserve memoria de sus valores pasados.
El decaimiento gradual de la ACF es característico de series con fuerte dependencia temporal y posible no estacionariedad, lo que sugiere la necesidad de aplicar diferenciación antes de identificar adecuadamente la estructura ARIMA.
plot_ly(comp, x = ~tiempo, y = ~tendencia,
type = "scatter", mode = "lines",
line = list(color = "#2ECC71", width = 2.5),
name = "Tendencia STL") %>%
layout(
title = list(text = "Componente de tendencia — O₃ diario Cali (2018)",
font = list(color = "#2ECC71", size = 15)),
xaxis = list(title = "Mes", tickformat = "%b", dtick = "M1"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "x unified"
)La tendencia extraída por la descomposición STL revela un comportamiento claramente no lineal a lo largo del año, con tres fases bien diferenciadas.
La serie arranca el 1 de enero en 22.10 µg/m³ y desciende hasta su mínimo absoluto el 10 de febrero con 18.25 µg/m³. Esta caída inicial coincide temporalmente con la primera temporada de lluvias del Valle del Cauca, aunque la tendencia STL únicamente refleja el movimiento suave de la serie sin incorporar información meteorológica directa.
Desde ese mínimo, la tendencia inicia una recuperación sostenida y pronunciada: pasa por 30.68 µg/m³ el 1 de julio y alcanza su pico máximo el 25 de agosto con 42.97 µg/m³, un incremento de 24.72 µg/m³ respecto al punto más bajo del año. Este ascenso es temporalmente consistente con la temporada seca principal de la región, periodo en que una mayor radiación solar favorece los procesos fotoquímicos asociados a la generación de ozono, si bien esta observación es de carácter exploratorio.
A partir del pico de agosto la tendencia corrige parcialmente, sin llegar a los niveles del inicio del año: al 31 de diciembre se sitúa en 35.95 µg/m³, es decir, 13.85 µg/m³ por encima del valor de enero. El patrón real del año describe entonces una trayectoria de descenso inicial → ascenso pronunciado → corrección parcial, lo que implica que el nivel tendencial del O₃ en Cali cerró 2018 considerablemente más alto que al inicio del período analizado.
En conjunto, la amplitud total de la tendencia —de 18.25 a 42.97 µg/m³— indica que esta componente captura una variación de fondo relevante dentro de la serie, cuyo comportamiento es consistente con el régimen climático estacional característico de Cali, aunque la información disponible no permite cuantificar su influencia específica.
plot_ly(comp, x = ~tiempo, y = ~estacional,
type = "scatter", mode = "lines",
line = list(color = "#1ABC9C", width = 1.5),
name = "Estacionalidad STL") %>%
layout(
title = list(text = "Componente de estacionalidad — O₃ diario Cali (2018)",
font = list(color = "#1ABC9C", size = 15)),
xaxis = list(title = "Mes", tickformat = "%b", dtick = "M1"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "x unified"
)La componente estacional presenta un patrón semanal perfectamente repetitivo, resultado directo del procedimiento de descomposición STL al especificar una frecuencia de 7 observaciones. Esto explica por qué la componente estacional aparece como un patrón perfectamente repetitivo a lo largo de todo el año: el mismo ciclo de 7 días se repite de forma matemáticamente idéntica durante las 52 semanas.
La amplitud total del ciclo es de apenas 1.08 µg/m³ — de -0.63 a +0.46 µg/m³ — lo que representa menos del 4% de la media global de la serie (29.58 µg/m³) y es considerablemente menor que la desviación estándar de 7.49 µg/m³. Estos valores indican que la estacionalidad semanal tiene un peso marginal en la explicación del comportamiento del O₃, siendo la tendencia y los residuos los componentes que concentran la mayor parte de la variabilidad de la serie.
Dentro del ciclo, el sábado registra el efecto positivo más alto (+0.4552 µg/m³), seguido del domingo (+0.3220) y el lunes (+0.2228), mientras que el martes presenta el efecto negativo más pronunciado (-0.6271 µg/m³) y el viernes el segundo más bajo (-0.4252). Si bien este contraste entre fin de semana y días entre semana podría ser consistente con diferencias en los patrones de actividad urbana de la ciudad, la magnitud del efecto es tan reducida que no permite atribuirle una importancia práctica significativa.
En conjunto, esta componente refleja principalmente una característica del método de descomposición más que una estacionalidad semanal con sustancia física demostrable en los datos.
plot_ly(comp, x = ~tiempo, y = ~residuo,
type = "scatter", mode = "lines",
line = list(color = "#8E44AD", width = 1.5),
name = "Residuo STL") %>%
add_lines(x = range(comp$tiempo), y = c(0, 0),
line = list(color = "#F39C12", dash = "dash", width = 1.2),
name = "Cero", showlegend = FALSE) %>%
layout(
title = list(text = "Componente de residuos — O₃ diario Cali (2018)",
font = list(color = "#8E44AD", size = 15)),
xaxis = list(title = "Fecha", tickformat = "%b", dtick = "M1"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "x unified"
)La componente residual representa lo que queda de la serie una vez extraídas la tendencia y la estacionalidad: la variación que el modelo de descomposición no logra explicar estructuralmente.
La media de los residuos es prácticamente cero (0.006 µg/m³), lo que indica que la descomposición no dejó ningún sesgo sistemático sin capturar. La desviación estándar es de 4.19 µg/m³, y el 79.5% de las observaciones se mantiene dentro del rango ±5 µg/m³, lo que refleja un comportamiento mayoritariamente contenido a lo largo del año.
Sin embargo, el gráfico revela tres picos que se salen claramente de ese rango: el más pronunciado ocurre el 27 de agosto con +15.50 µg/m³, seguido del 29 de septiembre con +11.07 µg/m³ y el 24 de julio con -12.88 µg/m³. Estos eventos representan apenas el 2.5% de la serie, correspondientes a observaciones que no fueron adecuadamente explicadas por las componentes de tendencia y estacionalidad identificadas por el modelo.
Comparando la dispersión entre períodos, la desviación estándar de los residuos en febrero–marzo es de 3.96 µg/m³, mientras que en julio–septiembre asciende a 4.54 µg/m³. Esta diferencia sugiere una mayor dispersión residual durante el segundo semestre, aunque no es lo suficientemente marcada como para concluir la existencia de heterocedasticidad.
Visto en perspectiva con las otras dos componentes, los números cuentan una historia clara sobre el peso relativo de cada una:
| Componente | Medida | Valor |
|---|---|---|
| Tendencia | Amplitud | 24.72 µg/m³ |
| Residuo | Desv. estándar | 4.19 µg/m³ |
| Estacionalidad | Amplitud | 1.08 µg/m³ |
La tendencia domina ampliamente la variación de la serie, el residuo conserva una magnitud relevante, y la estacionalidad semanal resulta casi anecdótica en comparación. En conjunto, los residuos presentan una media cercana a cero y una dispersión moderada, lo que sugiere que la descomposición STL logró capturar una parte importante de la estructura sistemática de la serie. Sin embargo, la presencia de fluctuaciones persistentes y algunos eventos extremos indica que aún podría existir dependencia temporal remanente, aspecto que será evaluado posteriormente mediante las herramientas de identificación y diagnóstico del modelo ARIMA.
acf_vals <- acf(serie, lag.max = 25, plot = FALSE)
acf_df <- data.frame(
lag = as.numeric(acf_vals$lag[-1]),
acf = as.numeric(acf_vals$acf[-1])
)
ic <- qnorm(0.975) / sqrt(length(serie))
plot_ly(acf_df, x = ~lag, y = ~acf, type = "bar",
marker = list(color = "rgba(46,204,113,0.8)",
line = list(color = "#2ECC71", width = 1)),
name = "ACF") %>%
add_lines(x = c(0, 26), y = c(ic, ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
add_lines(x = c(0, 26), y = c(-ic, -ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
layout(
title = list(text = "ACF — O₃ original Cali (2018)",
font = list(color = "#2ECC71", size = 15)),
xaxis = list(title = "Lag", dtick = 5),
yaxis = list(title = "ACF", range = c(-0.15, 0.7)),
hovermode = "x unified",
bargap = 0.3
)Antes de identificar un modelo ARIMA es necesario verificar si la serie es estacionaria, es decir, si su nivel promedio y su variabilidad se mantienen constantes a lo largo del tiempo. Esta condición es un requisito fundamental para que el modelo funcione correctamente.
La prueba Dickey-Fuller Aumentada (ADF) evalúa la presencia de una raíz unitaria en la serie. Su hipótesis nula establece que la serie posee raíz unitaria y, por tanto, es no estacionaria. Con un estadístico de -2.89 y un p-valor de 0.20, no se rechaza la hipótesis nula. Para que este resultado fuera concluyente en favor de la estacionariedad, el p-valor debería ser menor a 0.05 — al ser 0.20, la evidencia no es suficiente para descartar la presencia de raíz unitaria.
Este resultado se refuerza con lo observado en la ACF de la serie original. Las autocorrelaciones arrancan en 0.621 en el lag 1 y descienden de forma lenta y persistente: en el lag 4 aún se mantienen en 0.493 y en el lag 25 siguen siendo positivas y significativas en torno a 0.22, lo que indica que observaciones separadas por casi un mes conservan todavía una relación estadística entre sí. Este patrón de decaimiento lento es consistente con una posible no estacionariedad, aunque por sí solo no la demuestra.
Tomadas en conjunto, ambas evidencias — el p-valor del ADF y el patrón de decaimiento prolongado de la ACF — apuntan en la misma dirección y constituyen una base sólida para concluir que la serie no es estacionaria en su forma original, lo que hace necesario aplicar una diferenciación antes de identificar adecuadamente la estructura ARIMA.
miserie_df <- data.frame(
fecha = o3_diario$fecha[-1],
delta = as.numeric(miserie)
)
plot_ly(miserie_df, x = ~fecha, y = ~delta,
type = "scatter", mode = "lines",
line = list(color = "#2ECC71", width = 1.2),
name = "ΔO₃") %>%
add_lines(x = range(miserie_df$fecha), y = c(0, 0),
line = list(color = "#F39C12", dash = "dash", width = 1.2),
showlegend = FALSE) %>%
layout(
title = list(text = "Serie O₃ diferenciada — Cali (2018)",
font = list(color = "#2ECC71", size = 15)),
xaxis = list(title = "Fecha", tickformat = "%b", dtick = "M1"),
yaxis = list(title = "ΔO₃ (µg/m³)"),
hovermode = "x unified"
)Ante la evidencia de no estacionariedad, se aplica una diferenciación de orden 1, que consiste en reemplazar cada observación por el cambio respecto al día anterior. El resultado es una nueva serie que ya no mide la concentración de O₃ sino su variación diaria.
El gráfico de la serie diferenciada muestra un comportamiento visualmente distinto al de la serie original: las fluctuaciones oscilan de forma irregular alrededor de cero sin mostrar una tendencia sostenida en ninguna dirección, lo que es una primera señal favorable hacia la estacionariedad. Se identifican algunos picos puntuales que superan los ±15 µg/m³, siendo el más extremo una caída cercana a -30 µg/m³ alrededor de abril. Estos valores reflejan cambios diarios abruptos en la concentración de O₃ y sugieren la presencia de episodios excepcionales dentro de la serie.
La prueba ADF sobre la serie diferenciada proporciona evidencia estadística consistente con este resultado: el estadístico cae a -9.64 y el p-valor es menor a 0.01, por lo que se rechaza con contundencia la hipótesis nula de raíz unitaria. El aviso de R indica que el p-valor real es incluso menor al 0.01 reportado, lo que refuerza aún más la conclusión. En comparación con la serie original, cuyo p-valor fue 0.20, la reducción a un valor inferior a 0.01 muestra que la diferenciación eliminó gran parte de la dependencia asociada a la raíz unitaria. Con una sola diferenciación se obtiene evidencia estadística consistente con estacionariedad, por lo que resulta razonable establecer d = 1 dentro del modelo ARIMA.
acf_diff <- acf(miserie, lag.max = 30, plot = FALSE)
pacf_diff <- pacf(miserie, lag.max = 30, plot = FALSE)
ic <- qnorm(0.975) / sqrt(length(miserie))
acf_df <- data.frame(
lag = as.numeric(acf_diff$lag[-1]),
acf = as.numeric(acf_diff$acf[-1])
)
pacf_df <- data.frame(
lag = as.numeric(pacf_diff$lag),
pacf = as.numeric(pacf_diff$acf)
)
p_acf <- plot_ly(acf_df, x = ~lag, y = ~acf, type = "bar",
marker = list(color = "rgba(46,204,113,0.8)",
line = list(color = "#2ECC71", width = 1)),
name = "ACF") %>%
add_lines(x = c(0, 31), y = c(ic, ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
add_lines(x = c(0, 31), y = c(-ic, -ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
layout(
title = list(text = "ACF — O₃ diferenciado",
font = list(color = "#2ECC71", size = 13)),
xaxis = list(title = "Lag", dtick = 5),
yaxis = list(title = "ACF", range = c(-0.5, 0.3)),
bargap = 0.3
)
p_pacf <- plot_ly(pacf_df, x = ~lag, y = ~pacf, type = "bar",
marker = list(color = "rgba(26,188,156,0.8)",
line = list(color = "#1ABC9C", width = 1)),
name = "PACF") %>%
add_lines(x = c(0, 31), y = c(ic, ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
add_lines(x = c(0, 31), y = c(-ic, -ic),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
layout(
title = list(text = "PACF — O₃ diferenciado",
font = list(color = "#1ABC9C", size = 13)),
xaxis = list(title = "Lag", dtick = 5),
yaxis = list(title = "PACF", range = c(-0.5, 0.3)),
bargap = 0.3
)
subplot(p_acf, p_pacf, nrows = 1, shareY = FALSE, titleX = TRUE) %>%
layout(hovermode = "x unified")Una vez confirmada la estacionariedad de la serie diferenciada, el análisis de la ACF y la PACF permite identificar la estructura de dependencia temporal que guiará la selección de los órdenes p (componente autorregresiva) y q (componente de medias móviles) del modelo ARIMA.
La ACF muestra un único lag significativo fuera de la banda de confianza (±0.103) en los rezagos iniciales: el lag 1 con -0.429. A partir del lag 2 los valores caen dentro de la banda y permanecen allí hasta los rezagos 28, 29 y 30, donde aparecen tres valores marginalmente significativos (-0.121, +0.133, -0.134). Estos últimos, al estar aislados en rezagos muy alejados, se interpretan con cautela ya que en series de 364 observaciones es esperable que algunos lags crucen la banda por azar. El corte abrupto de la ACF después del primer rezago es compatible con la presencia de una componente de medias móviles de orden 1, por lo que resulta razonable considerar q = 1 como candidato inicial.
La PACF presenta varios rezagos significativos en los primeros lags, con valores de -0.429 en el lag 1, -0.209 en el lag 2, -0.191 en el lag 3 y -0.121 en el lag 5. Aunque el patrón no muestra un corte perfectamente definido — dado que el lag 4 no resulta significativo pero el lag 5 sí —, la persistencia de autocorrelaciones parciales significativas en los primeros rezagos sugiere explorar componentes autorregresivas de bajo orden, especialmente p = 1, 2 y 3.
En conjunto, la ACF sugiere q = 1 como candidato inicial y la PACF sugiere explorar p entre 1 y 3, lo que define el espacio de modelos candidatos que serán comparados formalmente en la siguiente sección mediante criterios de información.
tabla_modelos <- data.frame(
Modelo = c("ARIMA(1,1,2)", "ARIMA(2,1,1)", "ARIMA(1,1,1)",
"ARIMA(0,1,1)", "ARIMA(2,1,2)"),
AICc = c(2208.089, 2208.312, 2208.780, 2209.395, 2210.359),
RMSE = c(5.493, 5.495, 5.515, 5.537, 5.495),
MAE = c(4.265, 4.265, 4.307, 4.344, 4.264)
)
kable(tabla_modelos,
caption = "Comparación de modelos candidatos — O₃ Cali 2018",
align = c("l", "r", "r", "r")) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed"),
position = "center") %>%
row_spec(0, bold = TRUE, color = "white", background = "#2ECC71") %>%
row_spec(1, bold = TRUE, background = "#D5F5E3") %>%
row_spec(seq(3, nrow(tabla_modelos), 2), background = "#F0FBF4") %>%
column_spec(1, bold = TRUE, color = "#2D2D2D") %>%
column_spec(2, color = "#2ECC71", bold = TRUE)| Modelo | AICc | RMSE | MAE |
|---|---|---|---|
| ARIMA(1,1,2) | 2208.089 | 5.493 | 4.265 |
| ARIMA(2,1,1) | 2208.312 | 5.495 | 4.265 |
| ARIMA(1,1,1) | 2208.780 | 5.515 | 4.307 |
| ARIMA(0,1,1) | 2209.395 | 5.537 | 4.344 |
| ARIMA(2,1,2) | 2210.359 | 5.495 | 4.264 |
Con base en las señales de la ACF y la PACF, se estimaron cinco
modelos candidatos, todos con d = 1 ya establecido,
variando los órdenes p y q dentro del rango sugerido. Adicionalmente,
auto.arima() fue consultado como referencia, sugiriendo un
ARIMA(1,1,2), que coincide con uno de los candidatos evaluados.
La comparación se realizó mediante el criterio de información
AICc, que penaliza la complejidad del modelo para evitar el
sobreajuste. Un valor menor de AICc indica un mejor balance entre ajuste
y parsimonia. Aunque las diferencias entre los modelos ARIMA(1,1,2),
ARIMA(2,1,1) y ARIMA(1,1,1) son pequeñas, el
ARIMA(1,1,2) obtuvo simultáneamente el menor
AICc (2208.089), el menor RMSE (5.493) y fue
además el modelo seleccionado por auto.arima(), por lo que
se tomó como candidato principal para la fase de diagnóstico
residual.
res <- as.numeric(residuals(modelo3))
idx <- seq_along(res)
acf_r <- acf(res, lag.max = 25, plot = FALSE)
acf_res_df <- data.frame(
lag = as.numeric(acf_r$lag[-1]),
acf = as.numeric(acf_r$acf[-1])
)
ic_r <- qnorm(0.975) / sqrt(length(res))
p1 <- plot_ly(x = idx, y = res,
type = "scatter", mode = "lines",
line = list(color = "#2ECC71", width = 1),
name = "Residuales") %>%
add_lines(x = range(idx), y = c(0, 0),
line = list(color = "#F39C12", dash = "dash", width = 1.2),
showlegend = FALSE) %>%
layout(
title = list(text = "Residuales en el tiempo",
font = list(color = "#2ECC71", size = 13)),
xaxis = list(title = "Observación"),
yaxis = list(title = "Residual")
)
p2 <- plot_ly(acf_res_df, x = ~lag, y = ~acf, type = "bar",
marker = list(color = "rgba(26,188,156,0.8)",
line = list(color = "#1ABC9C", width = 1)),
name = "ACF") %>%
add_lines(x = c(0, 26), y = c(ic_r, ic_r),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
add_lines(x = c(0, 26), y = c(-ic_r, -ic_r),
line = list(color = "#E74C3C", dash = "dash", width = 1.5),
showlegend = FALSE) %>%
layout(
title = list(text = "ACF residuales",
font = list(color = "#1ABC9C", size = 13)),
xaxis = list(title = "Lag"),
yaxis = list(title = "ACF"),
bargap = 0.3
)
p3 <- plot_ly(x = res, type = "histogram",
nbinsx = 30,
marker = list(color = "rgba(142,68,173,0.8)",
line = list(color = "#8E44AD", width = 1)),
name = "Residuales") %>%
layout(
title = list(text = "Distribución de residuales",
font = list(color = "#8E44AD", size = 13)),
xaxis = list(title = "Residual"),
yaxis = list(title = "Frecuencia")
)
subplot(
p1,
subplot(p2, p3, nrows = 1, shareY = FALSE, titleX = TRUE),
nrows = 2, titleX = TRUE, titleY = TRUE
) %>%
layout(hovermode = "x unified",
showlegend = FALSE)El diagnóstico de residuales permite verificar si el modelo capturó adecuadamente la estructura de dependencia de la serie. Un modelo bien especificado se espera que produzca residuales que se comporten como ruido blanco: sin autocorrelación, con media cercana a cero y distribución aproximadamente normal.
Gráfico de residuales en el tiempo
Los residuales oscilan de forma irregular alrededor de cero a lo largo de las 353 observaciones de entrenamiento, sin mostrar tendencias, ciclos ni cambios sistemáticos en su nivel. Se identifican algunos valores extremos puntuales — el más pronunciado ocurre alrededor de la observación 110, con una caída cercana a -25 — correspondientes a los episodios atípicos ya identificados en etapas previas del análisis. Fuera de estos eventos puntuales, la dispersión se mantiene relativamente estable a lo largo del tiempo.
ACF de los residuales
La ACF de los residuales muestra que todos los coeficientes de autocorrelación se mantienen dentro de las bandas de confianza, sin patrones sistemáticos visibles. Esto indica que la dependencia temporal observada en la serie original fue capturada adecuadamente por el modelo y que los residuales se comportan de manera compatible con ruido blanco.
Histograma de los residuales
La distribución de los residuales presenta una forma aproximadamente simétrica y campaniforme, coherente con una distribución normal. La cola derecha muestra una ligera extensión asociada a los valores extremos ya mencionados, pero el comportamiento central es regular.
La prueba de Ljung-Box evalúa formalmente si existe autocorrelación en los residuales considerando los primeros 10 rezagos en conjunto, descontando los 3 parámetros estimados del modelo ARIMA(1,1,2), por lo que la prueba evalúa autocorrelación remanente que no fue explicada por el ajuste. Su hipótesis nula establece que los residuales son independientes. Con un estadístico Q* = 5.262 y un p-valor de 0.628, no se rechaza la hipótesis nula, lo que constituye evidencia estadística de que los residuales no presentan autocorrelación significativa. Un p-valor de 0.628 está muy lejos del umbral de 0.05, lo que refuerza con solidez esta conclusión.
En conjunto, los resultados del gráfico temporal de los residuales, la ACF residual, el histograma y la prueba de Ljung-Box son consistentes con la ausencia de estructura remanente, lo que permite concluir que el ARIMA(1,1,2) está adecuadamente especificado para esta serie y puede utilizarse para realizar pronósticos.
sw <- shapiro.test(residuals(modelo3))
cat("Estadístico W:", round(sw$statistic, 5), "\n")## Estadístico W: 0.99434cat("p-valor:", round(sw$p.value, 3), "\n")## p-valor: 0.216La prueba de Shapiro-Wilk evalúa formalmente si los residuales del modelo siguen una distribución normal, complementando el diagnóstico visual del histograma. Su hipótesis nula establece que los datos provienen de una distribución normal, por lo que valores de p por encima de 0.05 indican que no hay evidencia suficiente para rechazarla.
El estadístico W = 0.994 es notablemente cercano a 1 — el valor máximo posible, que correspondería a una distribución perfectamente normal. Esta proximidad es consistente con una distribución aproximadamente normal de los residuales, sin desviaciones importantes respecto a la normalidad. El p-valor de 0.216, muy por encima del umbral de 0.05, confirma que no existe evidencia estadística para rechazar la normalidad de los residuales.
Tomados en conjunto con el resultado de la prueba de Ljung-Box (p = 0.628), los diagnósticos indican que los residuales del ARIMA(1,1,2) son compatibles con un proceso de ruido blanco, al no presentar evidencia de autocorrelación significativa y mostrar una distribución aproximadamente normal. Esto respalda la validez de los supuestos del modelo y la confiabilidad de los pronósticos e intervalos de confianza que se presentan a continuación.
pron_df <- data.frame(
dia = 354:365,
pred = as.numeric(pronostico$mean),
lo95 = as.numeric(pronostico$lower),
hi95 = as.numeric(pronostico$upper)
)
hist_df <- data.frame(
dia = 293:353,
o3 = as.numeric(tail(ventana, 61))
)
plot_ly() %>%
add_lines(data = hist_df, x = ~dia, y = ~o3,
line = list(color = "#1A6FA8", width = 1.5),
name = "Serie observada") %>%
add_ribbons(data = pron_df, x = ~dia,
ymin = ~lo95, ymax = ~hi95,
fillcolor = "rgba(46,158,143,0.25)",
line = list(color = "transparent"),
name = "IC 95%") %>%
add_lines(data = pron_df, x = ~dia, y = ~pred,
line = list(color = "#E07B39", width = 2.5, dash = "dash"),
name = "Pronóstico") %>%
layout(
title = list(text = "Pronóstico O₃ — ARIMA(1,1,2) — Cali (2018)",
font = list(color = "#1A6FA8", size = 15)),
xaxis = list(title = "Día del año"),
yaxis = list(title = "O₃ (µg/m³)"),
hovermode = "x unified",
legend = list(orientation = "h", x = 0, y = -0.2)
)pron_tabla <- data.frame(
Día = 354:365,
Pronóstico = round(as.numeric(pronostico$mean), 2),
`Límite inf. 95%` = round(as.numeric(pronostico$lower), 2),
`Límite sup. 95%` = round(as.numeric(pronostico$upper), 2)
)
kable(pron_tabla,
caption = "Valores pronosticados — O₃ Cali, días 354 a 365",
align = c("c","r","r","r")) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped","hover","condensed"),
position = "center") %>%
row_spec(0, bold = TRUE, color = "white", background = "#1A6FA8") %>%
row_spec(seq(2, 12, 2), background = "#EAF4FB") %>%
column_spec(2, color = "#E07B39", bold = TRUE)| Día | Pronóstico | Límite.inf..95. | Límite.sup..95. |
|---|---|---|---|
| 354 | 32.60 | 21.78 | 43.43 |
| 355 | 31.89 | 20.34 | 43.43 |
| 356 | 31.35 | 19.27 | 43.42 |
| 357 | 30.94 | 18.44 | 43.44 |
| 358 | 30.63 | 17.79 | 43.47 |
| 359 | 30.40 | 17.26 | 43.54 |
| 360 | 30.22 | 16.83 | 43.62 |
| 361 | 30.09 | 16.46 | 43.72 |
| 362 | 29.99 | 16.15 | 43.83 |
| 363 | 29.92 | 15.87 | 43.96 |
| 364 | 29.86 | 15.63 | 44.09 |
| 365 | 29.82 | 15.40 | 44.23 |
El modelo ARIMA(1,1,2) genera un pronóstico de 12 días hacia adelante (días 354 al 365 de la serie, correspondientes al 20 al 31 de diciembre de 2018). Los valores pronosticados descienden de forma gradual desde 32.60 µg/m³ en el día 354 hasta 29.82 µg/m³ en el día 365, lo cual resulta consistente con el comportamiento observado en las últimas semanas de la serie, donde las concentraciones mostraban una ligera disminución.
El intervalo de confianza al 95% se amplía progresivamente con el horizonte de pronóstico: el límite inferior desciende de 21.78 a 15.40 µg/m³, mientras que el límite superior aumenta ligeramente de 43.43 a 44.23 µg/m³, dando lugar a un ensanchamiento progresivo del intervalo de confianza. Esto refleja que la incertidumbre crece con la distancia temporal al último dato observado, comportamiento esperable en cualquier modelo de pronóstico y que no constituye una señal de debilidad del modelo sino una característica inherente a la naturaleza estocástica de las series de tiempo.
acc <- accuracy(pronostico, ventana2)
acc_tabla <- data.frame(
Conjunto = c("Entrenamiento", "Prueba"),
ME = round(acc[,"ME"], 3),
RMSE = round(acc[,"RMSE"], 3),
MAE = round(acc[,"MAE"], 3),
MAPE = round(acc[,"MAPE"], 3),
MASE = round(acc[,"MASE"], 3)
)
kable(acc_tabla,
caption = "Métricas de evaluación del pronóstico — ARIMA(1,1,2)",
align = c("l","r","r","r","r","r")) %>%
kable_styling(full_width = FALSE,
bootstrap_options = c("striped","hover","condensed"),
position = "center") %>%
row_spec(0, bold = TRUE, color = "white", background = "#1A6FA8") %>%
row_spec(2, background = "#D5F5E3", bold = TRUE) %>%
column_spec(3, color = "#E07B39", bold = TRUE)| Conjunto | ME | RMSE | MAE | MAPE | MASE | |
|---|---|---|---|---|---|---|
| Training set | Entrenamiento | 0.100 | 5.493 | 4.265 | 15.979 | 0.850 |
| Test set | Prueba | -0.544 | 5.451 | 4.167 | 15.873 | 0.831 |
Las métricas obtenidas en el conjunto de prueba son muy similares a las del entrenamiento (RMSE: 5.451 frente a 5.493; MAE: 4.167 frente a 4.265), lo que indica que el desempeño predictivo del modelo se mantuvo estable al aplicarse a observaciones no utilizadas durante la estimación. No obstante, dado que la evaluación se realizó sobre únicamente 12 observaciones, esta evidencia debe interpretarse con cautela.
El MAPE de 15.87% en prueba indica que, en promedio, el pronóstico se aleja un 15.87% del valor real — un nivel de error moderado, consistente con la variabilidad intrínseca de una serie ambiental diaria influenciada por múltiples factores meteorológicos no incorporados en el modelo. El MASE de 0.831, al ser menor que 1, indica que el modelo supera en precisión al pronóstico ingenuo de referencia, lo que respalda su utilidad para realizar pronósticos de corto plazo.
En conjunto, el modelo ARIMA(1,1,2) logró capturar adecuadamente la estructura temporal de la serie diaria de O₃ en Cali durante 2018. Los diagnósticos residuales respaldan la validez de sus supuestos principales y las métricas de evaluación muestran un desempeño predictivo estable en el horizonte analizado, aunque los resultados corresponden únicamente a un año de observaciones y deben interpretarse dentro de ese alcance temporal.
El análisis descriptivo de la serie diaria de O₃ en Cali durante 2018 permitió identificar una dinámica temporal caracterizada por un descenso al inicio del año, un ascenso sostenido hasta agosto y una disminución progresiva hacia finales del período. La tendencia resultó ser el componente dominante de la serie, capturando este ciclo de mediano plazo de forma clara. En contraste, el componente estacional semanal obtenido mediante la descomposición STL mostró una amplitud muy reducida (inferior a 1.1 µg/m³), lo que sugiere que el comportamiento del O₃ durante el período estudiado estuvo determinado principalmente por cambios graduales de mediano plazo más que por patrones asociados al día de la semana.
La distribución global de la serie presentó características aproximadamente simétricas, con una media de 29.58 µg/m³ y una mediana de 29.40 µg/m³, una asimetría positiva leve (0.20) y una curtosis cercana a la de una distribución normal (2.95). Se identificaron cinco valores atípicos distribuidos entre abril y agosto, siendo el más extremo el registrado el 27 de agosto (58.7 µg/m³).
La prueba ADF confirmó que la serie original no era
estacionaria (p = 0.20), condición corregida mediante una
diferenciación de primer orden que eliminó la evidencia
de raíz unitaria, aportando evidencia estadística consistente
con la estacionariedad (p < 0.01). Con base en los gráficos
ACF y PACF de la serie diferenciada y la referencia de
auto.arima(), se evaluaron cinco modelos candidatos,
seleccionando el ARIMA(1,1,2) por obtener
simultáneamente el menor AICc (2208.089) y los errores
de ajuste más bajos entre los modelos comparados.
El diagnóstico de residuales respaldó la validez del modelo: la prueba de Ljung-Box (p = 0.628) descartó la presencia de autocorrelación significativa y la prueba de Shapiro-Wilk (W = 0.994, p = 0.216) indicó que no existe evidencia estadística para rechazar la normalidad de los residuales. En conjunto, ambas pruebas son compatibles con el supuesto de ruido blanco, lo que valida los supuestos principales del modelo.
Los pronósticos para los últimos doce días de la serie mostraron una ligera disminución en las concentraciones proyectadas, de 32.60 a 29.82 µg/m³, coherente con el comportamiento observado al final del período de entrenamiento. El ensanchamiento progresivo de los intervalos de confianza al 95% — cuyo límite inferior desciende de 21.78 a 15.40 µg/m³ mientras el superior aumenta ligeramente de 43.43 a 44.23 µg/m³ — refleja el aumento natural de la incertidumbre conforme se amplía el horizonte de pronóstico. Las métricas de evaluación sobre la ventana de prueba (RMSE = 5.451, MAE = 4.167, MASE = 0.831) fueron muy similares a las del entrenamiento, lo que sugiere que el desempeño predictivo se mantuvo estable entre el conjunto de entrenamiento y el de prueba, aunque esta evidencia debe interpretarse con cautela dado que la evaluación se realizó sobre únicamente 12 observaciones.
Es importante reconocer que el modelo ARIMA se fundamenta exclusivamente en el comportamiento histórico de la serie y no incorpora variables externas que también pueden influir en la concentración de O₃, como las condiciones meteorológicas, las emisiones contaminantes o eventos extraordinarios. Para estudios futuros, se recomienda considerar series históricas más extensas y explorar modelos con variables exógenas — como temperatura, radiación solar, humedad y viento — que podrían mejorar tanto el poder predictivo como la comprensión de los factores que determinan el comportamiento del ozono troposférico en Cali.
En síntesis, los resultados obtenidos muestran que un modelo ARIMA(1,1,2) proporciona una representación adecuada de la dinámica temporal de la concentración diaria de O₃ en Cali durante 2018 y ofrece un desempeño predictivo razonable para pronósticos de corto plazo dentro del alcance de los datos analizados, respondiendo afirmativamente a la pregunta de investigación planteada.
Series de tiempo y modelado ARIMA
Box, G. E. P., & Jenkins, G. M. (1976). Time Series Analysis: Forecasting and Control (2nd ed.). Holden-Day.
Hyndman, R. J., & Athanasopoulos, G. (2018). Forecasting: Principles and Practice (2nd ed.). OTexts. Recuperado de https://otexts.com/fpp2/
Ozono troposférico y salud
Observatorio Europeo del Clima y la Salud — Climate-ADAPT. (2022). Efectos del ozono troposférico sobre la salud humana en un clima cambiante. Agencia Europea de Medio Ambiente. Recuperado de https://climate-adapt.eea.europa.eu/es/observatory/evidence/health-effects/ground-level-ozone
Madrid Salud. (2025). Ozono troposférico y salud. Ayuntamiento de Madrid. Recuperado de https://madridsalud.es/ozono-y-salud/
Clima del Valle del Cauca y régimen de lluvias
Ideam. (s.f.). Climatología trimestral de Colombia. Instituto de Hidrología, Meteorología y Estudios Ambientales. Recuperado de http://www.ideam.gov.co/documents/21021/21789/Climatolog%C3%ADa+Trimestral+para+Colombia
Vélez-Upegui, J. J., et al. (2022). Análisis de la influencia de la variabilidad climática en la precipitación de la cuenca del río Cali, Colombia. Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, 46(179), 168–182. Recuperado de http://www.scielo.org.co/scielo.php?script=sci_arttext&pid=S0012-73532022000200168
Toda Colombia. (2019). Clima del departamento del Valle del Cauca. Recuperado de https://www.todacolombia.com/departamentos-de-colombia/valle-del-cauca/clima.html