El estudio de las emisiones de dióxido de carbono (CO₂) se ha convertido en un tema central dentro de las discusiones ambientales a nivel mundial. Entre las distintas fuentes de emisiones, la aviación focaliza interés y preocupación, especialmente en naciones que dependen mucho del transporte aéreo, ya que este sector representa anualmente cerca del 2% de las emisiones globales de CO₂ y aporta aproximadamente un 4,9% al calentamiento global (Perlman, 2018). Además, se proyecta que el crecimiento del tráfico aéreo continúe incrementando estas emisiones, a pesar de las mejoras tecnológicas en eficiencia, lo que ha llevado a organismos internacionales como la Organización de Aviación Civil Internacional (OACI) y la Unión Europea a establecer programas específicos de mitigación.
En este trabajo se analiza la serie temporal mensual “Monthly CO₂ emissions from international and domestic flights” en Canadá, un país cuya geografía, economía y desarrollo histórico han hecho de la aviación un factor importante para su transporte y conectividad. La nación ha ampliado rápidamente su red de aeropuertos desde mediados del siglo XX, debido a su amplio territorio, el crecimiento del comercio internacional y la necesidad de incorporar comunidades remotas en el norte del país. Esto ha generado un aumento constante en el tráfico aéreo, haciendo que los vuelos sean una de las actividades con mayor crecimiento relativo en términos de emisiones como CORSIA y el EU ETS.
Reducir las emisiones de CO₂ es especialmente importante para Canadá. Según Environment and Climate Change Canada (2019), el país se está calentando más rápidamente que el promedio mundial, fenómeno que genera efectos significativos como la acelerada pérdida de hielo en el Ártico, el deshielo del permafrost y un aumento en la magnitud de los incendios forestales, lo que repercute en las comunidades vulnerables, ecosistemas y a la infraestructura. A este desafío se suma el hecho de que la aviación intensifica su impacto climático, ya que los vuelos liberan no sólo CO₂, sino también vapor de agua y óxidos de nitrógeno en altitudes elevadas, amplificando el efecto de calentamiento.
Al considerar los compromisos climáticos que Canadá ha asumido a nivel mundial, esta serie temporal adquiere relevancia. El país al firmar el Acuerdo de París, se comprometió a disminuir en gran medida sus emisiones de gases de efecto invernadero y a alcanzar la neutralidad de carbono para el año 2050 (Minister of Transport, 2022). Adicionalmente, toma parte en CORSIA (Carbon Offsetting and Reduction Scheme for International Aviation), un programa internacional cuyo objetivo es reducir y compensar las emisiones de vuelos internacionales. Al mismo tiempo, a nivel nacional el gobierno federal ha implementado medidas como el precio al carbono, regulaciones sobre la eficiencia de las aerolíneas y estrategias para promover el uso de combustibles sostenibles de aviación (SAF).
En este contexto, el análisis de series temporales es fundamental para entender el comportamiento de las emisiones, su dinámica, así como para anticipar posibles escenarios futuros. Este estudio permite detectar patrones estacionales, como los crecimientos en verano por mayor actividad turística o en invierno caracterizado por periodos festivos y viajes a lugares más cálidos. Además, facilita identificar tendencias a largo plazo junto con cambios asociados al crecimiento del tráfico aéreo o a modificaciones en las políticas de transporte. Reconocer y modelar estas conductas es esencial para distinguir entre variaciones temporales y patrones propios de cada estación del año, especialmente en un país que ha asumido compromisos climáticos.
A través de modelos ARIMA es posible capturar componentes como tendencia, estacionalidad y ruido para generar pronósticos confiables que permitan prever el comportamiento de las emisiones en los próximos meses. Contar con esta información es fundamental para evaluar si el desempeño del sector es coherente con los objetivos medioambientales de Canadá, así como para respaldar decisiones vinculadas a la adopción de tecnologías emergentes como los combustibles sostenibles de aviación, ajustes normativos, mejoras en eficiencia o acciones correctivas dentro del sector.
Una serie de tiempo es un conjunto de observaciones registradas de manera secuencial en intervalos regulares, su estructura permite identificar patrones como tendencias, estacionalidad, ciclos y variaciones aleatorias. El estudio de este tipo de datos facilita comprender el comportamiento histórico de una variable y anticipar cómo podría evolucionar. En este trabajo se utiliza la serie mensual de emisiones de CO₂ provenientes de vuelos internacionales que salen de Canadá, lo que permite observar cómo varía el comportamiento del sector aeronáutico a lo largo de distintos periodos y en diversas circunstancias.
Para modelar este comportamiento se empleó la metodología ARIMA (Autoregressive Integrated Moving Average), la cual se usa mucho en el análisis de series temporales. Este tipo de modelos combina componentes autorregresivos (AR), una parte de diferenciación (I) para lograr estacionariedad y un término de medias móviles (MA), en conjunto, estos elementos permiten capturar la dependencia temporal de la serie y modelar sus variaciones. Su utilidad radica en que posibilitan generar predicciones basadas en la relación entre valores pasados y presentes, lo que resulta importante cuando se busca anticipar cambios o evaluar si una trayectoria es coherente con metas o escenarios planteados.
Dado que la serie analizada presenta un patrón estacional marcado, asociado al incremento de vuelos en meses como los de verano, se optó por utilizar un modelo SARIMA (ARIMA con componente estacional). Este modelo incorpora parámetros estacionales (P,D,Q)s que permiten capturar la repetición de patrones en intervalos regulares, en este caso con una periodicidad mensual. La elección del orden del modelo, tanto en los componentes estacionales como no estacionales (p,d,q), se realizó mediante la función auto.arima() y la revisión de los gráficos ACF y PACF. Para seleccionar el modelo más adecuado se emplearon diferentes criterios como AIC, BIC y AICc.
CComo se mencionó anteriormente, para el análisis se utilizó la serie mensual “Monthly CO₂ emissions from international and domestic flights”, publicada por Our World in Data y elaborada a partir del conjunto de datos experimentales de emisiones del transporte aéreo desarrollado por la OCDE (2025). Este indicador, el internacional, estima las toneladas mensuales de CO₂ emitidas por los vuelos comerciales de pasajeros que salen desde Canadá, siguiendo el criterio internacional que atribuye las emisiones al país de origen del vuelo. Cada punto del tiempo representa el volumen total de CO₂ asociado a la actividad aérea del país durante ese mes, lo que permite interpretar estas observaciones como un reflejo del movimiento turístico, comercial y social que caracteriza a la población en Canadá.
library(DT)
library(dplyr)
tabla_serie <- co2_aviones %>%
mutate(
CO2 = round(CO2, 3)
) %>%
select(Fecha, CO2)
datatable(
tabla_serie,
rownames = FALSE,
colnames = c("Fecha", "CO₂"),
options = list(
pageLength = 8,
scrollX = TRUE,
autoWidth = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0),
list(className = 'dt-right', targets = 1),
list(width = "150px", targets = 0),
list(width = "120px", targets = 1)
)
)
) %>%
formatRound(columns = 2, digits = 3)El periodo registrado y tomado para analizar va desde enero de 2019 hasta diciembre de 2024, una ventana interesante porque permite observar la dinámica aérea antes de la pandemia, su colapso durante 2020 y parte de 2021, y luego la reactivación progresiva del sector. Cada dato mensual proviene del conjunto “Air transport CO₂ emissions (experimental)” de la OCDE, que OWID procesa para estandarizarlo, ajustar nombres, verificar unidades y presentarlo de manera accesible. La elaboración del indicador se basa en estimaciones construidas a partir del consumo de combustible, el número de vuelos y los factores de emisión característicos de la aviación civil.
library(plotly)
Serie <- plot_ly(
data = co2_aviones,
x = ~Fecha,
y = ~CO2,
type = 'scatter',
mode = 'lines+markers',
line = list(color = "#2F4F6A", width = 3),
marker = list(color = "#2F4F6A", size = 6)
) %>%
layout(
title = list(
text = "Serie Temporal de Emisiones de CO₂ (Aviación - Canadá)",
x = 0.5,
font = list(size = 20, color = "#2F4F6A")
),
xaxis = list(
title = "Fecha",
tickfont = list(size = 12),
titlefont = list(size = 14, color = "#2F4F6A")
),
yaxis = list(
title = "CO₂",
tickfont = list(size = 12),
titlefont = list(size = 14, color = "#2F4F6A")
),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
SerieLos cambios observados en la serie temporal de las emisiones de CO₂ entre 2019 y 2024 pueden explicarse por sucesos históricos, transformaciones en la movilidad aérea y elementos estacionales. Durante 2019 e inicios de 2020 las emisiones se mantienen elevadas debido a que Canadá es un país donde una parte importante de la población viaja a menudo al exterior, especialmente hacia Estados Unidos, México, el Caribe y Europa. En estos meses se aprecian picos típicos de épocas con alta movilidad; en verano, cuando el clima cálido impulsa los viajes turísticos; y en invierno cuando muchas personas tienden a viajar a lugares más templados o visitan familiares en el extranjero.
A partir de marzo de 2020 la serie muestra un descenso abrupto y significativo. Esta caída coincide con la llegada de la pandemia de COVID-19, que provocó el cierre de fronteras, restricciones en los desplazamientos, cancelación de vuelos comerciales y confinamientos tanto en Canadá como en gran parte del mundo. El sector aeronáutico fue uno de los más perjudicados a nivel mundial, y esto se refleja claramente en la serie, con niveles mínimos que se mantuvieron durante casi todo este año.
En 2021 comienza una recuperación lenta pero sostenida. Con el levantamiento progresivo de restricciones, la reapertura de fronteras, la reactivación del turismo y las actividades económicas, la demanda de vuelos vuelve a aumentar. La serie muestra un crecimiento gradual de las emisiones, acompañado de variaciones estacionales anuales como los incrementos durante los meses de verano y los periodos festivos, y caídas en épocas con menor movilidad. Para 2022 y 2023 el crecimiento es más marcado, impulsado por el turismo reprimido, el retorno de los vuelos internacionales y la normalización del transporte aéreo.
Finalmente, en 2024 las emisiones se estabilizan en valores altos, aunque con variaciones intermensuales habituales del sector. Las cifras de este año indican que la aviación canadiense ha recuperado, casi por completo, su actividad previa a la pandemia, impulsada por la reapertura mundial, el aumento del turismo, los viajes de negocios y la necesidad de mantener conexiones internacionales. Las oscilaciones observadas son el resultado de la estacionalidad del tráfico aéreo, a las dinámicas propias del mercado y a factores climáticos, como el aumento del turismo local en verano o los desplazamientos hacia climas cálidos durante la época invernal. Esta dinámica evidencia cómo la actividad aérea incide directamente en los niveles de emisiones y cómo eventos extraordinarios, como una pandemia, pueden modificar drásticamente la tendencia histórica del transporte internacional.
Para profundizar en el comportamiento de la serie temporal, se calcularon las principales medidas de tendencia central y dispersión de las emisiones de CO₂. La siguiente tabla resume los valores más relevantes de la variable:
library(DT)
library(htmltools)
x <- co2_aviones$CO2
n <- sum(!is.na(x))
tabla_vertical <- tibble(
Estadística = c(
"Mínimo",
"Q1 (25%)",
"Mediana",
"Media",
"Q3 (75%)",
"Máximo",
"Desviación Estándar",
"N° Observaciones"
),
Valor = c(
min(x, na.rm = TRUE),
quantile(x, 0.25, na.rm = TRUE),
median(x, na.rm = TRUE),
mean(x, na.rm = TRUE),
quantile(x, 0.75, na.rm = TRUE),
max(x, na.rm = TRUE),
sd(x, na.rm = TRUE),
n
)
)
tabla_vertical$Valor[1:7] <- round(tabla_vertical$Valor[1:7], 2)
datatable(
tabla_vertical,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '250px', targets = 0),
list(width = '150px', targets = 1)
)
),
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)En la tabla se observa que la media de la serie se ubica cerca de 769540 toneladas mensuales, mientras que la mediana alcanza un valor mayor de 872688 toneladas. Esta diferencia indica que la distribución está influenciada por los valores extremadamente bajos registrados durante 2020, lo que hace que la media caiga sin reflejar el comportamiento habitual del resto de los meses. El valor mínimo registrado, 139269 toneladas, corresponde precisamente a uno de los momentos más críticos de la pandemia, cuando el tráfico aéreo prácticamente se detuvo. Por otro lado, el valor máximo, 1327710 toneladas, coincide con un periodo de alta movilidad, seguramente asociado a los meses de verano o a fechas festivas, en un contexto previo al impacto de la pandemia.
La desviación estándar, cercana a 324280 toneladas, confirma que las emisiones presentan una variabilidad considerable, lo que es coherente con un sector que depende mucho de factores como las temporadas turísticas, condiciones climáticas, cambios regulatorios y, como se vio en 2020, eventos extraordinarios capaces de alterar completamente su comportamiento. Estas medidas permiten comprender que la serie no se mueve dentro de un rango estrecho, sino que responde a ciclos amplios y a situaciones externas que generan cambios bruscos.
library(ggplot2)
library(plotly)
library(scales)
p_co2 <- ggplot(co2_aviones, aes(y = CO2)) +
geom_boxplot(
fill = "#89ACC6",
color = "#2c3e50",
outlier.color = "#284C64",
outlier.size = 3,
alpha = 0.85
) +
labs(
title = "",
y = "Emisiones de CO₂",
x = ""
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold", color = "#2F4F6A", hjust = 0.5),
axis.title.y = element_text(size = 13, color = "#2F4F6A", face = "bold"),
axis.text = element_text(size = 11, color = "#2F4F6A")
)
ggplotly(p_co2, tooltip = c("y"))El boxplot permite visualizar rápidamente la dispersión de las emisiones de CO₂. La amplitud de la caja muestra que los valores varían bastante entre los meses, lo que coincide con lo ya señalado en la tabla. La mediana se ubica en una zona relativamente alta dentro del rango intercuartílico, lo cual refleja que, en la mayoría de los meses, las emisiones se mantienen en niveles elevados de actividad aérea.
El valor atípico que aparece separado del resto corresponde al periodo en el que la pandemia detuvo gran parte de las operaciones aéreas y como resultado las emisiones cayeron a niveles extremadamente bajos. Su presencia destaca porque rompe por completo con el comportamiento normal de la serie y marca un punto de quiebre dentro de su historia reciente.
En conjunto, la tabla descriptiva y el boxplot muestran que las emisiones asociadas a la aviación en Canadá presentan una variabilidad notable, con meses de alta actividad, otros de actividad moderada y un conjunto de valores extremadamente bajos que corresponden a un evento extraordinario. Todo lo anterior ayuda a confirmar que el comportamiento de la serie está fuertemente marcado por patrones estacionales y, al mismo tiempo, por shocks externos que interrumpen la dinámica tradicional del sector.
La descomposición STL permite separar la serie en tres componentes: la tendencia, la estacionalidad y el componente aleatorio. Esta herramienta facilita entender cómo se comportan las emisiones mensuales de CO₂ a lo largo del tiempo y qué parte de las variaciones responde a patrones regulares y cuáles a cambios inesperados.
stl_descomp <- stl(Serie_ts, s.window = "periodic")
df_descomp <- data.frame(
Fecha = as.Date(time(Serie_ts)),
Observado = as.numeric(Serie_ts),
Tendencia = as.numeric(stl_descomp$time.series[, "trend"]),
Estacional = as.numeric(stl_descomp$time.series[, "seasonal"]),
Aleatorio = as.numeric(stl_descomp$time.series[, "remainder"])
)
p1 <- plot_ly(df_descomp, x = ~Fecha, y = ~Observado,
type = 'scatter', mode = 'lines',
line = list(color = "#2F4F6A")) %>%
layout(title = "Observado")
p2 <- plot_ly(df_descomp, x = ~Fecha, y = ~Tendencia,
type = 'scatter', mode = 'lines',
line = list(color = "#2F4F6A")) %>%
layout(title = "Tendencia")
p3 <- plot_ly(df_descomp, x = ~Fecha, y = ~Estacional,
type = 'scatter', mode = 'lines',
line = list(color = "#2F4F6A")) %>%
layout(title = "Estacional")
p4 <- plot_ly(df_descomp, x = ~Fecha, y = ~Aleatorio,
type = 'scatter', mode = 'lines',
line = list(color = "#2F4F6A")) %>%
layout(title = "Aleatorio")
subplot(p1, p2, p3, p4, nrows = 4, shareX = TRUE, titleY = TRUE) %>%
layout(
title = list(
text = "Descomposición de la Serie Temporal",
x = 0.5,
font = list(size = 22, color = "#2F4F6A")
),
showlegend = FALSE,
margin = list(t = 90)
)El componente observado reproduce la forma general de la serie original, niveles altos antes de 2020, una caída brusca durante la pandemia y una recuperación progresiva que se estabiliza entre 2023 y 2024. Esta primera lectura confirma lo ya visto en el análisis descriptivo, pero la descomposición permite profundizar en cómo se construye ese comportamiento.
La tendencia muestra con claridad la ruta de largo plazo de las emisiones. Antes de 2020 se mantiene en valores elevados y relativamente estables; sin embargo, con la llegada de la pandemia se observa un descenso significativo que refleja la paralización casi total del tráfico aéreo internacional. La tendencia alcanza su punto más bajo hacia finales de 2020 y, a partir de 2021, comienza una recuperación gradual que se extiende durante 2022 y 2023. Para 2024 la tendencia parece haberse estabilizado nuevamente en niveles altos, lo que indica que el sector ha retomado gran parte de su actividad anterior a la pandemia, la que se consideraba normal.
En cuanto al componente estacional, se aprecian patrones que se repiten cada año de forma marcada. Los picos de mayor actividad se concentran en verano, cuando el turismo internacional aumenta, y también en diciembre, coincidiendo con las vacaciones de invierno y los desplazamientos por festividades. Por el contrario, se observan valores negativos en meses de primavera temprana y otoño, periodos en los que la demanda aérea suele disminuir. Este comportamiento confirma la presencia de una estacionalidad fuerte y constante a lo largo de los años analizados, evidenciando que las subidas y bajadas mensuales no son aleatorias, sino que obedecen a un patrón estrechamente vinculado con el comportamiento de los pasajeros.
Por último, el componente aleatorio abarca las variaciones que no pueden explicarse ni por la tendencia ni por la estacionalidad. En este caso, destacan los valores inusualmente altos en 2020, lo cual es coherente con la irrupción repentina de la pandemia, un evento que no forma parte de ningún patrón histórico de la serie. Después de ese año, el componente aleatorio vuelve a mostrar oscilaciones más moderadas, asociadas a fluctuaciones habituales del sistema, como variaciones en la demanda internacional, ajustes regulatorios o cambios en las condiciones de viaje.
El gráfico estacional permite visualizar cómo se comportan las emisiones de CO₂ a lo largo del año y qué meses concentran los mayores y menores niveles de actividad aérea internacional. Esta representación facilita comparar directamente las curvas de cada año y observar si los patrones se repiten o si existen desviaciones marcadas asociadas a eventos extraordinarios.
library(dplyr)
library(lubridate)
library(plotly)
df_season <- co2_aviones %>%
mutate(
Year = year(Fecha),
Month_num = month(Fecha)
)
colores_CO2 <- c(
"#A9B9C6",
"#8299AB",
"#668299",
"#546B7D",
"#415362",
"#2F3C46"
)
años <- sort(unique(df_season$Year))
fig_season <- plot_ly()
for (i in seq_along(años)) {
df_temp <- df_season %>% filter(Year == años[i])
fig_season <- fig_season %>%
add_trace(
data = df_temp,
x = ~Month_num,
y = ~CO2,
type = "scatter",
mode = "lines+markers",
name = paste0(años[i]),
line = list(color = colores_CO2[i], width = 3),
marker = list(color = colores_CO2[i], size = 7)
)
}
fig_season <- fig_season %>%
layout(
title = list(
text = "Patrón Estacional de las Emisiones de CO₂ en Canadá",
x = 0.5,
font = list(size = 20, color = "#2F4F6A")
),
xaxis = list(
title = "Mes",
tickvals = 1:12,
ticktext = month.abb
),
yaxis = list(
title = "Emisiones CO₂"
),
legend = list(
orientation = "v",
x = 1.05, y = 1,
font = list(color = "#2F4F6A")
),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
fig_seasonComo se puede observar, las emisiones siguen un ciclo anual consistente. Se identifica un periodo de mayor actividad hacia mediados de año, asociado al incremento de los desplazamientos internacionales durante la temporada de vacaciones, cuando muchas personas viajan aprovechando el clima cálido y el calendario escolar. Hacia finales del año aparece un segundo subidón, relacionado con los viajes por celebraciones de diciembre, visitas a familiares y desplazamientos hacia destinos de clima más templado. Entre estos dos momentos de mayor actividad también hay meses más tranquilos, especialmente durante la primavera temprana y parte del otoño. Estos períodos representan transiciones esperables entre temporadas, cuando el flujo de viajeros disminuye debido a que no coinciden con vacaciones extendidas ni con condiciones climáticas que incentiven un aumento notable de la movilidad aérea. Esta dinámica cíclica se repite casi todos los años y confirma la existencia de una estacionalidad sólida.
La excepción evidente es el año 2020. La curva correspondiente presenta una caída abrupta que rompe totalmente con la estructura estacional habitual. Desde marzo de ese año, las emisiones se desploman y permanecen en niveles mínimos debido al cierre de fronteras, la cancelación masiva de vuelos y la paralización del turismo internacional provocada por la pandemia. A diferencia del resto de los años, 2020 no muestra ni los picos asociados al verano ni los incrementos de diciembre, lo que lo convierte en un caso atípico dentro de toda la serie.
A partir de 2021, las curvas comienzan a recuperar su forma característica. Aunque los niveles siguen siendo más bajos que los de 2019, se observa nuevamente una estructura cíclica reconocible. En 2022 y, cn mayor claridad, en 2023 y 2024, el patrón estacional retoma su forma habitual con mayor claridad, lo que sugiere que la demanda internacional de vuelos se ha estabilizado y que la movilidad aérea regresó a comportamientos similares a los previos a la pandemia. Además, se aprecia una tendencia progresiva hacia niveles cada vez más altos, acercándose nuevamente a la estabilidad que se tenía antes del impacto sanitario.
En conjunto, este análisis confirma que la serie presenta un ciclo anual muy marcado, donde los cambios responden tanto a factores turísticos y climáticos como a las dinámicas sociales propias del país. La estacionalidad es un componente fundamental de la serie y solo se ve alterada por eventos externos extraordinarios, como ocurrió en 2020.
El análisis del comportamiento de la serie temporal es un paso indispensable antes de aplicar modelos como ARIMA o SARIMA. La forma en que los datos evolucionan en el tiempo determina si un modelo puede ajustarse adecuadamente y si los pronósticos generados serán confiables. Para comprender esa estructura es necesario identificar tendencias, estacionalidades, ciclos repetitivos o quiebres estructurales que modifican el comportamiento habitual de la serie. También se debe verificar si los datos cumplen con la condición de estacionariedad o si necesitan transformaciones para lograrla. Por todo esto, estudiar la dinámica temporal es una etapa fundamental en el proceso de modelado.
Para este proceso, la serie se dividió en dos ventanas. La ventana de entrenamiento abarca desde enero de 2019 hasta febrero de 2024, cubriendo 62 períodos, este tramo captura distintas fases del comportamiento de las emisiones, como la estabilidad previa a 2020, el colapso provocado por la pandemia y la recuperación posterior, donde la estacionalidad vuelve a hacerse evidente. La ventana de prueba corresponde a los diez meses restantes de 2024 y se utiliza para comparar los valores reales con los pronósticos y verificar si el modelo logra reproducir la forma general de la serie.
library(plotly)
Serie <- xts(co2_aviones$CO2,
order.by = co2_aviones$Fecha)
ventana1 <- window(Serie, end = "2024-03-15")
ventana2 <- window(Serie, start = "2024-03-15", end = "2024-12-15")
df1 <- data.frame(
Fecha = index(ventana1),
CO2 = as.numeric(ventana1)
)
df2 <- data.frame(
Fecha = index(ventana2),
CO2 = as.numeric(ventana2)
)
ventanas <- plot_ly() %>%
add_trace(
data = df1,
x = ~Fecha,
y = ~CO2,
type = "scatter",
mode = "lines+markers",
line = list(color = "#2F3C46", width = 3),
marker = list(color = "#2F3C46", size = 6),
name = "Enero 2019–Febrero 2024"
) %>%
add_trace(
data = df2,
x = ~Fecha,
y = ~CO2,
type = "scatter",
mode = "lines+markers",
line = list(color = "#89ACC6", width = 3),
marker = list(color = "#89ACC6", size = 6),
name = "Marzo 2024–Diciembre 2024"
) %>%
layout(
title = list(
text = "",
x = 0.5,
font = list(size = 22, color = "#2F3C46")
),
xaxis = list(
title = "Fecha",
titlefont = list(color = "#2F3C46", size = 16),
tickfont = list(color = "#2F3C46", size = 12)
),
yaxis = list(
title = "CO₂",
titlefont = list(color = "#2F3C46", size = 16),
tickfont = list(color = "#2F3C46", size = 12)
),
plot_bgcolor = "white",
paper_bgcolor = "white",
legend = list(
x = 0.98,
y = 0.98,
xanchor = "right",
yanchor = "top",
bgcolor = "rgba(255,255,255,0.8)",
font = list(color = "#2F3C46", size = 12)
),
margin = list(t = 90)
)
ventanasUna vez segmentada la serie, se analizaron los gráficos de autocorrelación (ACF) y autocorrelación parcial (PACF) para identificar la dependencia temporal presente en los datos. El ACF muestra cuánto influyen los valores pasados en el comportamiento actual y si esa influencia se mantiene durante muchos rezagos. El PACF revela la relación directa entre un rezago y el valor presente dejando fuera el efecto de los rezagos intermedios.
Estos gráficos ayudan a identificar si la serie presenta tendencia o estacionalidad. Cuando la ACF disminuye de forma lenta, suele ser señal de que la serie no es estacionaria; y cuando se repite un patrón alrededor de rezagos específicos, como el 12, indica estacionalidad anual en series mensuales. Los rezagos destacados en la PACF orientan la selección de los parámetros p y q del modelo ARIMA.
library(plotly)
v <- as.numeric(ventana1)
acf_obj <- acf(v, plot = FALSE, lag.max = 60)
pacf_obj <- pacf(v, plot = FALSE, lag.max = 60)
ci <- qnorm((1 + 0.90)/2) / sqrt(length(v))
y_acf_min <- min(acf_obj$acf) - 0.1
y_pacf_min <- min(pacf_obj$acf) - 0.1
acf_plot <- plot_ly(showlegend = FALSE) %>%
add_bars(
x = acf_obj$lag[,1,1],
y = acf_obj$acf[,1,1],
marker = list(color = "#2F3C46")
) %>%
add_lines(x = c(0,60), y = c(ci,ci),
line = list(color = "#89ACC6", dash="dash")) %>%
add_lines(x = c(0,60), y = c(-ci,-ci),
line = list(color = "#89ACC6", dash="dash")) %>%
add_annotations(
text = "ACF",
x = 30,
y = y_acf_min,
xref = "x",
yref = "y",
showarrow = FALSE,
font = list(size = 14, color = "#2F3C46")
) %>%
layout(
xaxis = list(title = "Lag"),
yaxis = list(showticklabels = TRUE),
plot_bgcolor = "white",
paper_bgcolor = "white"
)
pacf_plot <- plot_ly(showlegend = FALSE) %>%
add_bars(
x = pacf_obj$lag[,1,1],
y = pacf_obj$acf[,1,1],
marker = list(color = "#2F3C46")
) %>%
add_lines(x = c(0,60), y = c(ci,ci),
line = list(color = "#89ACC6", dash="dash")) %>%
add_lines(x = c(0,60), y = c(-ci,-ci),
line = list(color = "#89ACC6", dash="dash")) %>%
add_annotations(
text = "PACF",
x = 30,
y = y_pacf_min,
xref = "x",
yref = "y",
showarrow = FALSE,
font = list(size = 14, color = "#2F3C46")
) %>%
layout(
xaxis = list(title = "Lag"),
yaxis = list(showticklabels = TRUE),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
subplot(acf_plot, pacf_plot, nrows = 1, shareY = FALSE)En este caso la ACF muestra una caída progresiva que confirma la presencia de una tendencia y la falta de estacionariedad; también se observa un agrupamiento alrededor del rezago doce que coincide con una estacionalidad anual presente en la serie original. En la PACF destaca únicamente el rezago uno como el de mayor influencia mientras los demás se mantienen cercanos a cero; esto señala que la relación directa más relevante es con el mes anterior y que la serie original no muestra una estructura autorregresiva simple. Ambos gráficos confirman que la serie no es estacionaria y necesita ser diferenciada para estabilizar su comportamiento. Además, el patrón que reaparece cada doce meses señala la necesidad de considerar también una diferenciación estacional.
Después se aplicó la prueba Augmented Dickey-Fuller (ADF) para evaluar de manera numérica si la serie era o no estacionaria. Una serie estacionaria mantiene su nivel y su variabilidad a lo largo del tiempo mientras que una serie no estacionaria suele mostrar tendencias claras, cambios de nivel o variaciones que se amplían con el tiempo. La prueba ADF permite formalizar esta idea al verificar la presencia de raíz unitaria. Cuando la prueba no es significativa se concluye que la serie no es estacionaria y requiere de una transformación mediante diferenciación. En los modelos ARIMA este proceso se resume en el parámetro d, el cual indica cuántas veces se diferencia la serie para estabilizarla, los otros parámetros p y q provienen de la PACF y el ACF respectivamente.
library(DT)
library(tibble)
library(htmltools)
tabla_adf <- tibble(
Estadística = c(
"Estadístico Dickey–Fuller",
"Lag order (rezagos utilizados)",
"p-value",
"Resultado"
),
Valor = c(
-1.7088,
3,
0.6927,
"No estacionaria"
)
)
tabla_adf$Valor[1:3] <- round(as.numeric(tabla_adf$Valor[1:3]), 4)
datatable(
tabla_adf,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '270px', targets = 0),
list(width = '160px', targets = 1)
)
),
caption = "TABLA ADF SIN DIFERENCIAR",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)El resultado de la prueba ADF aplicada a la serie original arrojó un valor no significativo, lo cual confirma que la serie no es estacionaria. El resultado coincide con lo observado tanto en la gráfica como en los correlogramas donde aparecen patrones repetitivos que reflejan tanto tendencia como estacionalidad, esta primera conclusión permitió avanzar hacia la diferenciación como forma de estabilizar la serie.
El primer camino consistió en diferenciar la serie con un rezago uno para retirar la tendencia que se apreciaba en la gráfica original. El comportamiento creciente y la caída abrupta de 2020 indicaban que la serie tenía variaciones sostenidas en el nivel, por lo que este paso era necesario para intentar estabilizar su dinámica.
diff_lag1 <- diff(ventana1) %>% na.omit()
library(plotly)
df_diff <- data.frame(
Fecha = as.Date(time(diff_lag1)),
Valor = as.numeric(diff_lag1)
)
Serie_diff <- plot_ly(
data = df_diff,
x = ~Fecha,
y = ~Valor,
type = "scatter",
mode = "lines+markers",
line = list(color = "#2F3C46", width = 3),
marker = list(color = "#2F3C46", size = 6)
) %>%
layout(
title = list(
text = "Serie diferenciada - primera diferencia",
x = 0.5,
font = list(size = 20, color = "#2F3C46")
),
xaxis = list(
title = "Fecha",
titlefont = list(color = "#2F3C46"),
tickfont = list(color = "#2F3C46")
),
yaxis = list(
title = "Valor",
titlefont = list(color = "#2F3C46"),
tickfont = list(color = "#2F3C46")
),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
Serie_diffUna vez aplicada la diferenciación, se volvió a calcular la prueba ADF. En esta ocasión el p-value fue de 0.01, lo que confirmó que la serie diferenciada con rezago uno sí cumplía con la condición de estacionariedad.
library(DT)
library(tibble)
library(htmltools)
tabla_adf <- tibble(
Estadística = c(
"Estadístico Dickey–Fuller",
"Lag order (rezagos utilizados)",
"p-value",
"Resultado"
),
Valor = c(
-4.3544,
3,
0.01,
"Estacionaria"
)
)
tabla_adf$Valor[1:3] <- round(as.numeric(tabla_adf$Valor[1:3]), 4)
datatable(
tabla_adf,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '270px', targets = 0),
list(width = '160px', targets = 1)
)
),
caption = "Resumen prueba ADF (Augmented Dickey–Fuller)",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)A pesar de haber obtenido una serie estacionaria por medio de la diferenciación con rezago uno no se dio por terminado el análisis. La estacionalidad presente en la serie original seguía siendo un hallazgo demasiado importante como para no tomarlo en cuenta. En los datos sin transformar los picos y valles se repetían con bastante claridad año tras año, especialmente fuera del impacto de la pandemia, mientras que al aplicar la primera diferenciación ese comportamiento desaparecía. La serie se estabiliza, sí, pero sin tomar en cuenta una característica fundamental de su comportamiento.
Por ese motivo se tomó un segundo camino de trabajo partiendo nuevamente de la serie original. Se aplicó una diferenciación con rezago doce para retirar directamente el componente anual y luego una segunda diferenciación con rezago uno para ajustar la tendencia residual. Esta alternativa no surgió porque la primera diferenciación hubiera fallado sino porque la estructura estacional era tan marcada en la serie original que necesitaba ser tratada de manera específica para construir después un modelo capaz de reproducir con mayor precisión los ciclos anuales de las emisiones.
diff_lag12 <- diff(ventana1, lag=12) %>% na.omit()
library(plotly)
df_diff12 <- data.frame(
Fecha = as.Date(time(diff_lag12)),
Valor = as.numeric(diff_lag12)
)
fig_diff12 <- plot_ly(
data = df_diff12,
x = ~Fecha,
y = ~Valor,
type = "scatter",
mode = "lines+markers",
line = list(color = "#2F3C46", width = 3),
marker = list(color = "#2F3C46", size = 6)
) %>%
layout(
title = list(
text = "Serie diferenciada - diferencia estacionaria con lag 12",
x = 0.5,
font = list(size = 20, color = "#2F3C46")
),
xaxis = list(title = "Fecha"),
yaxis = list(title = "Valor"),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
fig_diff12diff2_lag12 <- diff(diff_lag12) %>% na.omit()
df_diff12_2 <- data.frame(
Fecha = as.Date(time(diff2_lag12)),
Valor = as.numeric(diff2_lag12)
)
fig_diff12_2 <- plot_ly(
data = df_diff12_2,
x = ~Fecha,
y = ~Valor,
type = "scatter",
mode = "lines+markers",
line = list(color = "#2F3C46", width = 3),
marker = list(color = "#2F3C46", size = 6)
) %>%
layout(
title = list(
text = "Serie diferenciada 2 - diferencia estacionaria con lag 12",
x = 0.5,
font = list(size = 20, color = "#2F3C46")
),
xaxis = list(title = "Fecha"),
yaxis = list(title = "Valor"),
plot_bgcolor = "white",
paper_bgcolor = "white",
margin = list(t = 90)
)
fig_diff12_2Después de esta segunda ruta se evaluó nuevamente la estacionariedad por medio de la prueba ADF. La serie diferenciada solo con rezago doce no resultó estacionaria y presentó un p-value de 0.8001. Luego de aplicar la segunda diferenciación con rezago uno la prueba ADF arrojó un p-value de 0.0217 confirmando que la serie ya era estacionaria.
library(DT)
library(tibble)
library(htmltools)
tabla_adf <- tibble(
Diferenciación = c(
"Diferenciación lag 12",
"Diferenciación lag 1"
),
Estadístico_DF = c(
-1.4385,
-3.8825
),
Lag_order = c(
3,
3
),
`p-value` = c(
0.8001,
0.02175
),
Resultado = c(
"No estacionaria",
"Estacionaria"
)
)
tabla_adf$Estadístico_DF <- round(tabla_adf$Estadístico_DF, 4)
tabla_adf$`p-value` <- round(tabla_adf$`p-value`, 4)
datatable(
tabla_adf,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:4),
list(width = '200px', targets = 0),
list(width = '140px', targets = 1),
list(width = '120px', targets = 2),
list(width = '120px', targets = 3),
list(width = '140px', targets = 4)
)
),
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:4,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)Posterior a obtener las dos series estacionarias mediante las rutas de diferenciación propuestas, el siguiente paso consistió en examinar estructuras de modelado adecuadas para capturar la dinámica temporal de las emisiones. Para lograrlo fue necesario retomar los hallazgos previos del análisis de la serie original, con sus oscilaciones anuales y su tendencia pronunciada, y combinarlos con el comportamiento resultante después de la diferenciación.
La serie original mostraba una estacionalidad marcada, con picos y caídas que se repetían cada año y que se vuelven más claras después de 2022. En cambio, las dos series diferenciadas revelaron comportamientos distintos. La serie diferenciada con rezago uno eliminó por completo la tendencia y dejó una dinámica mucho más estable, pero al mismo tiempo borró la estacionalidad que era tan evidente en los datos originales. Por el otro camino, la diferenciación estacional con rezago doce quitó directamente el ciclo anual, aunque todavía quedaban variaciones amplias propias del periodo de pandemia, después de aplicar la segunda diferenciación con rezago uno la serie terminó de estabilizarse alrededor de un nivel constante. En ambas diferenciaciones salieron picos que representan el abrupto valor de 2020 que no se pudo apaciguar.
Tener estas dos formas estacionarias permitió avanzar en dos direcciones para el modelado. La primera, usando la serie diferenciada con rezago uno, condujo al análisis de modelos ARIMA sin componente estacional. La segunda se construyó a partir de la serie diferenciada estacionalmente, que era la base adecuada para desarrollar los modelos SARIMA.
Inicialmente, antes de seleccionar cualquier modelo, se recurrió a la función auto.arima(), cuyo objetivo era ofrecer un punto de partida. Si bien esta herramienta debe usarse con precaución porque selecciona modelos basándose únicamente en criterios automáticos como el AIC y no siempre reconoce rupturas estructurales, valores atípicos o estacionalidad compleja, resulta útil como referencia preliminar. Su objetivo es generar un modelo base que permita comparar los resultados con los modelos construidos manualmente y verificar si la estructura sugerida coincide con lo observado en la serie, el ACF y el PACF.
Al aplicar la función, el procedimiento sugirió como primera opción un ARIMA(0,1,0). Este resultado indica que el algoritmo reconoció la tendencia como el componente dominante de la serie y consideró suficiente una sola diferenciación para estabilizar su comportamiento, sin necesidad de incluir términos autorregresivos ni de medias móviles. En otras palabras, el modelo interpretó que la serie no mantenía una dependencia fuerte con sus valores pasados ni con los errores previos, y que su estructura principal podía capturarse únicamente a través de la diferenciación.
library(DT)
library(tibble)
library(htmltools)
tabla_arima010 <- tibble(
Elemento = c(
"Estructura",
"sigma² (Varianza de residuos)",
"Log-Likelihood",
"AIC",
"AICc",
"BIC"
),
Valor = c(
"ARIMA(0,1,0)",
1.273e+10,
-796.21,
1594.42,
1594.49,
1596.53
)
)
tabla_arima010$Valor[2:6] <- round(as.numeric(tabla_arima010$Valor[2:6]), 2)
datatable(
tabla_arima010,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '280px', targets = 0),
list(width = '150px', targets = 1)
)
),
caption = "Resumen del modelo auto.arima",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)El comportamiento de los indicadores confirma esta lectura. El valor de sigma cuadrado, que corresponde a la varianza de los residuos del modelo, es decir, la diferencia entre los valores reales y las predicciones, resultó ser alto. En modelos ARIMA, valores más bajos indican que el ajuste explica bien la dinámica de la serie porque deja residuos pequeños; sin embargo, en este caso la varianza residual es considerable. Esto indica que, aunque la diferenciación eliminó la tendencia principal, queda aún una porción importante de variabilidad que el modelo no logra capturar. Esta variación restante sugiere la presencia de otras dependencias temporales que no fueron recogidas por el modelo y que, por tanto, indican la necesidad de explorar estructuras más complejas o incorporar factores adicionales fuera del modelo.
El Log-Likelihood ofrece una medida sobre qué tan probable es que los datos provengan de la estructura estimada. En este caso, el valor obtenido, -796.21, es negativo y se encuentra lejos de cero, lo cual sugiere que el modelo no representa de manera especialmente sólida el proceso que genera la serie. Aunque este criterio por sí solo no determina la calidad del modelo, sí señala que logra explicar solo parte del comportamiento observado. En general, valores más cercanos a cero indican una mayor plausibilidad del modelo, por lo que este resultado refuerza la idea de que existen alternativas capaces de capturar la dinámica temporal con mayor precisión.
Los criterios AIC, AICc y BIC también aportaron señales en la misma dirección.
El AIC es una medida del ajuste que reduce la complejidad del modelo, obtuvo un valor de 1594.42, buen ajuste general. Este criterio establece que, a medida que el AIC disminuye, el modelo mejora en términos de equilibrio entre ajuste y simplicidad. En este caso, el valor indica que el modelo logra capturar la tendencia general, pero no es suficientemente eficiente frente a alternativas más complejas que podrían explicar mejor la estructura interna de la serie.
El AICc es una variante del AIC diseñada para muestras pequeñas y que restringe aún más la complejidad cuando el número de observaciones es limitado, su resultado es prácticamente igual al AIC. Esto confirma que el modelo no está sobreajustado a la muestra empleada; simplemente refleja una capacidad explicativa reducida, coherente con su estructura simple.
El BIC restringe la complejidad de forma más estricta que el AIC y para el cual valores bajos son preferibles, también presentó un valor elevado. Esto señala que, a pesar de la parsimonia del ARIMA(0,1,0), su simplicidad no compensa la falta de componentes que permitan modelar patrones presentes en los datos, como la dependencia a corto plazo o la estacionalidad.
En conjunto, estos criterios muestran que el modelo dado por el auto.arima() funciona únicamente como un punto de referencia inicial. No incorpora ninguna de las dinámicas que el análisis previo había revelado como la estacionalidad anual. Por ese motivo, avanzar hacia modelos más elaborados no sólo era pertinente, sino necesario para representar de manera adecuada el comportamiento real de las emisiones.
Después de revisar los resultados del auto.arima y confirmar que el modelo ARIMA(0,1,0) funcionaba únicamente como un punto de referencia inicial, se examinó nuevamente la estructura de autocorrelación, esta vez aplicada a la serie diferenciada con rezago uno. Este paso es esencial para identificar qué rezagos seguían influyendo directamente en el comportamiento actual de las emisiones y, por tanto, cuáles debían considerarse en los modelos ARIMA construidos manualmente.
diff_lag1 <- diff(ventana1) %>% na.omit()
library(plotly)
v_diff <- as.numeric(diff_lag1)
acf_raw <- acf(v_diff, plot = FALSE, lag.max = 61, type="correlation")
pacf_raw <- pacf(v_diff, plot = FALSE, lag.max = 61)
ci <- qnorm((1 + 0.90)/2) / sqrt(length(v_diff))
lags_acf <- acf_raw$lag[,1,1][-1]
vals_acf <- acf_raw$acf[,1,1][-1]
y_acf_min <- min(vals_acf) - 0.05
acf_plot <- plot_ly(showlegend=FALSE) %>%
add_bars(
x = lags_acf,
y = vals_acf,
marker=list(color="#2F3C46")
) %>%
add_lines(x=c(1,61), y=c(ci,ci), line=list(color="#89ACC6", dash="dash")) %>%
add_lines(x=c(1,61), y=c(-ci,-ci), line=list(color="#89ACC6", dash="dash")) %>%
add_annotations(
text="ACF", x=31, y=y_acf_min, xref="x", yref="y",
showarrow=FALSE, font=list(size=14, color="#2F3C46")
) %>%
layout(
xaxis=list(title="Lag"),
yaxis=list(showticklabels=TRUE),
plot_bgcolor="white",
paper_bgcolor="white"
)
lags_pacf <- pacf_raw$lag[,1,1]
vals_pacf <- pacf_raw$acf[,1,1]
y_pacf_min <- min(vals_pacf) - 0.05
pacf_plot <- plot_ly(showlegend=FALSE) %>%
add_bars(
x = lags_pacf,
y = vals_pacf,
marker=list(color="#2F3C46")
) %>%
add_lines(x=c(1,61), y=c(ci,ci), line=list(color="#89ACC6", dash="dash")) %>%
add_lines(x=c(1,61), y=c(-ci,-ci), line=list(color="#89ACC6", dash="dash")) %>%
add_annotations(
text="PACF", x=31, y=y_pacf_min, xref="x", yref="y",
showarrow=FALSE, font=list(size=14, color="#2F3C46")
) %>%
layout(
xaxis=list(title="Lag"),
yaxis=list(showticklabels=TRUE),
plot_bgcolor="white",
paper_bgcolor="white"
)
subplot(acf_plot, pacf_plot, nrows=1, shareY=FALSE)En los correlogramas del primer camino de diferenciación, donde se aplicó únicamente el rezago uno, la serie perdió por completo la tendencia que dominaba los datos originales. En la ACF la autocorrelación cae con rapidez y desaparece la prolongación que caracterizaba a la serie sin transformar, lo que confirma que la tendencia fue eliminada de manera adecuada. Sin embargo, destaca un rezago dominante alrededor del quinto periodo, que sobresale tanto en la ACF como en la PACF. Este comportamiento indica que la serie mantiene una dependencia estructural con respecto a los valores observados cinco meses atrás, una dinámica que no fue capturada por el modelo sugerido automáticamente y que sugiere una estructura interna más compleja. Además, este rezago puede relacionarse con patrones operativos propios del sector, como ciclos en la oferta y demanda de vuelos o ajustes mensuales recurrentes, lo que le da coherencia tanto estadística como contextual dentro del comportamiento de las emisiones.
Siguiendo la posibilidad de que haya efectos a corto plazo, se examinó además el modelo ARIMA(1,1,1), que incorpora un término autorregresivo y uno de medias móviles de primer orden para captar relaciones a corto plazo. Sin embargo, la presencia significativa del rezago 5 llevó a tener en cuenta modelos más extensos, en especial los modelos ARIMA(5,1,4) y ARIMA(5,1,5). Incluir órdenes cinco en los componentes autorregresivos y de medias móviles se debe directamente al descubrimiento empírico de un rezago dominante, lo cual coincide con el enfoque teórico que sugiere incluir aquellos rezagos que muestran mayor significancia. No obstante, el riesgo de sobreajuste se incrementa al ampliarse la cantidad de parámetros, por lo que estos modelos fueron analizados comparativamente con alternativas más parsimoniosas para asegurar un equilibrio apropiado entre ajuste interno, estabilidad y capacidad de predicción fuera de muestra.
library(DT)
library(tibble)
library(htmltools)
tabla_arima_lag1 <- tibble(
Modelo = c(
"Modelo 1",
"Modelo 2",
"Modelo 3",
"Modelo 4"
),
`Orden (p,d,q)` = c(
"(0,1,0)",
"(1,1,1)",
"(5,1,5)",
"(5,1,4)"
)
)
datatable(
tabla_arima_lag1,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '200px', targets = 0),
list(width = '150px', targets = 1)
)
),
caption = "Modelos ARIMA lag 1",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)En el segundo camino, la estructura estacional se retiró mediante una diferenciación con rezago doce, y luego se aplicó una segunda diferenciación con rezago uno para eliminar la variación restante de tendencia. El resultado en los correlogramas es distinto tanto al primer camino como a la serie original. En la ACF desaparece el patrón repetitivo que se observaba cada doce meses, lo que confirma que el componente estacional fue eliminado de manera adecuada. En la PACF aparecen algunos rezagos aislados que sobresalen ligeramente, pero no siguen una secuencia reconocible ni mantienen la regularidad propia de un patrón estacional. Esto sugiere que las fluctuaciones restantes corresponden al comportamiento propio del ruido después de haber aplicado las dos diferenciaciones, y no a una estacionalidad residual.
library(plotly)
v2 <- as.numeric(diff2_lag12)
acf_raw <- acf(v2, plot = FALSE, lag.max = 49, type = "correlation")
pacf_raw <- pacf(v2, plot = FALSE, lag.max = 49)
ci <- qnorm((1 + 0.90)/2) / sqrt(length(v2))
lags_acf <- acf_raw$lag[,1,1][-1]
vals_acf <- acf_raw$acf[,1,1][-1]
lags_pacf <- pacf_raw$lag[,1,1]
vals_pacf <- pacf_raw$acf[,1,1]
y_acf_min <- min(vals_acf) - 0.05
y_pacf_min <- min(vals_pacf) - 0.05
acf_diff12 <- plot_ly(showlegend=FALSE) %>%
add_bars(
x = lags_acf,
y = vals_acf,
marker=list(color="#2F3C46")
) %>%
add_lines(x=c(1,49), y=c(ci,ci), line=list(color="#89ACC6", dash="dash")) %>%
add_lines(x=c(1,49), y=c(-ci,-ci), line=list(color="#89ACC6", dash="dash")) %>%
add_annotations(
text="ACF",
x=25, y=y_acf_min,
xref="x", yref="y",
showarrow=FALSE,
font=list(size=14, color="#2F3C46")
) %>%
layout(
xaxis=list(title="Lag"),
yaxis=list(showticklabels=TRUE),
plot_bgcolor="white",
paper_bgcolor="white"
)
pacf_diff12 <- plot_ly(showlegend=FALSE) %>%
add_bars(
x = lags_pacf,
y = vals_pacf,
marker=list(color="#2F3C46")
) %>%
add_lines(x=c(1,49), y=c(ci,ci), line=list(color="#89ACC6", dash="dash")) %>%
add_lines(x=c(1,49), y=c(-ci,-ci), line=list(color="#89ACC6", dash="dash")) %>%
add_annotations(
text="PACF",
x=25, y=y_pacf_min,
xref="x", yref="y",
showarrow=FALSE,
font=list(size=14, color="#2F3C46")
) %>%
layout(
xaxis=list(title="Lag"),
yaxis=list(showticklabels=TRUE),
plot_bgcolor="white",
paper_bgcolor="white"
)
subplot(acf_diff12, pacf_diff12, nrows=1, shareY=FALSE)El comportamiento de estos correlogramas indica que la diferenciación estacional cumplió su objetivo y que la estacionalidad ya no se encuentra en la parte irregular de la serie. Por lo tanto, el modelo que se construya desde este segundo camino podrá capturar la estacionalidad de forma explícita a través de los parámetros estacionales del SARIMA.
De acuerdo con la teoría de Box–Jenkins, los picos fuertes en la ACF permiten identificar muchos MA estacionales, mientras que los picos en la PACF ayudan a distinguir los términos AR. Al mismo tiempo, los rezagos notables de baja frecuencia indicaron que se incluyeran términos adicionales AR y MA no estacionales. En conjunto, estos factores permitieron definir la estructura de ciertos modelos SARIMA adecuados para captar el comportamiento anual y la dinámica mensual de la serie.
En tal sentido, se calcularon tres modelos SARIMA propuestos. El primer modelo, SARIMA(0,1,0)(1,1,0)[12], fue sugerido de manera preliminar por auto.arima(). Este modelo es una propuesta compacta que incluye únicamente un componente autorregresivo estacional (P = 1) junto con diferenciación ordinaria y diferenciación estacional. Seguidamente, se desarrolló un segundo modelo, SARIMA(1,1,1)(1,1,0)[12], que incluye dinámicas a corto plazo mediante términos no estacionales AR y MA (p = 1, q = 1), pero mantiene el componente autorregresivo estacional. Por último, se calculó un tercer modelo, SARIMA(1,1,1)(1,1,1)[12], que agrega un componente MA estacional extra (Q = 1), de acuerdo con la existencia de rezagos estacionales relevantes en la PACF y la ACF. Estos tres modelos, de dificultad progresiva, permiten evaluar de forma sistemática cómo mejora el ajuste al incluir componentes estacionales adicionales.
library(DT)
library(tibble)
library(htmltools)
tabla_sarima_lag12 <- tibble(
Modelo = c(
"Modelo 5 (Sarima)",
"Modelo 6 (Sarima)",
"Modelo 7 (Sarima)"
),
`Orden (p,d,q)(P,D,Q)[12]` = c(
"(0,1,0)(1,1,0)[12]",
"(1,1,1)(1,1,0)[12]",
"(1,1,1)(1,1,1)[12]"
)
)
datatable(
tabla_sarima_lag12,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '220px', targets = 0),
list(width = '220px', targets = 1)
)
),
caption = "Modelos SARIMA lag 12",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)Para culminar esta sección, la selección del modelo final se realizó mediante la comparación de los criterios de información (AIC), el análisis de residuos con la prueba Ljung-Box y las medidas de rendimiento. Los hallazgos indicaron que el modelo SARIMA(1,1,1)(1,1,1)[12], con residuos no autocorrelacionados (p-valor > 0,05), cumplió con los supuestos del modelo y tuvo el AIC más bajo de todos los candidatos. Este modelo también obtuvo los mejores indicadores de exactitud, con los valores más bajos de MAE, RMSE y MAPE. En resumen, estos componentes permiten concluir que el modelo SARIMA(1,1,1)(1,1,1)[12] refleja correctamente la estructura temporal de la serie, reproduciendo tanto la dinámica mensual como el patrón estacional anual. Por lo tanto, se escoge como el modelo más sólido y confiable para llevar a cabo proyecciones futuras de emisiones. A continuación se desglosa este análisis.
El criterio de Información Akaike (AIC) se utiliza como una medida estadística para la comparación de modelos, en particular para los pronósticos SARIMA Y ARIMA. Su principal propósito es equilibrar el ajuste de un modelo con la complejidad del mismo. Un modelo que tiene un AIC más reducido utiliza menos parámetros para representar de manera efectiva la estructura de la serie. Esto es esencial en la modelación de series temporales, ya que los modelos demasiado complejos pueden sobre ajustar los datos históricos y, por lo tanto, disminuir su capacidad predictiva. Fundamental, el AIC restringe la complejidad que no es necesaria y prefiere las soluciones eficientes; por ende, valores de AIC bajos señalan un mejor desempeño del modelo.
tabla_aic <- tibble(
Modelo = c(
"Modelo 1",
"Modelo 2",
"Modelo 3",
"Modelo 4",
"Modelo 5 (Sarima)",
"Modelo 6 (Sarima)",
"Modelo 7 (Sarima)"
),
AIC = c(
1594.419,
1597.191,
1606.375,
1603.019,
1302.027,
1300.79,
1294.228
)
)
tabla_aic$AIC <- round(tabla_aic$AIC, 3)
datatable(
tabla_aic,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:1),
list(width = '250px', targets = 0),
list(width = '250px', targets = 1)
)
),
caption = "Criterio AIC por modelo",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:1,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)Al analizar los valores de AIC que aparecen en las tablas de ARIMA Y SARIMA, se nota que algunos modelos tienen valores mucho más altos que otros. Esto señala que, a pesar de que tienen la capacidad de adecuarse a los datos desde un punto de vista matemático, necesitan una cantidad de parámetros que no se justifica por la calidad del ajuste alcanzado. Modelos como ARIMA (5,1,4) y ARIMA (5,1,5), por ejemplo, tienen un AIC alto; esto indica que están sobre parametrizados: añaden una complejidad innecesaria para explicar patrones que podrían describirse con un modelo más sencillo. En el modelo SARIMA (0,1,0)(1,1,0)[12], se aprecia una conducta similar, su AIC más alto señala que su capacidad de captar la estacionalidad es menos efectiva en comparación con otras opciones más completas.
Con base en el criterio de información de Akaike el modelo SARIMA (1,1,1)(1,1,1)[12] es el que mejor se ajusta porque tiene el AIC más bajo entre todos los modelos analizados. Este hallazgo señala que el modelo mencionado logra captar la estructura estacional y la dinámica temporal no estacional de la serie, con una relación adecuada entre ajuste y cantidad de parámetros. Este modelo, en el marco de la serie estudiada que presenta patrones estacionales fuertes y cambios sistemáticos durante todo el año, ofrece la capacidad predictiva más robusta. Esto permite que se produzcan estimaciones más acorde con la conducta real del fenómeno y, por lo tanto, más eficaces para el análisis y la toma de decisiones.
La función checkresiduals analiza los residuos, que son los faltantes que el modelo no es capaz de explicar dentro de la serie de tiempo, para valorar la idoneidad de un modelo tipo ARIMA o SARIMA. Para que un modelo se considere idóneo, los residuos deben ser ruido blanco. Ruido blanco es una serie sin autocorrelación, con media cero y varianza constante. La prueba de Ljung-Box, parte de checkresiduals, evalúa la existencia de autocorrelación para los residuos. En este punto, se tendrá un p-value alto (p > 0,05) por lo tanto, no se podrá rechazar la hipótesis nula que dice que los residuos son ruido blanco. Esta condición implica que el modelo de ARIMA o SARIMA logra interpretar correctamente la estructura temporal de la serie por lo que no hay patrones que se encuentren por debajo de la normalidad. En el caso contrario, un valor p bajo indica que se trata de un modelo no satisfactorio.
library(DT)
library(tibble)
library(htmltools)
tabla_check_lag1 <- tibble(
Modelo = c(
"Modelo 1",
"Modelo 2",
"Modelo 3",
"Modelo 4",
"Modelo 5 (Sarima)",
"Modelo 6 (Sarima)",
"Modelo 7 (Sarima)"
),
`ARIMA (p,d,q)` = c(
"(0,1,0)",
"(1,1,1)",
"(5,1,5)",
"(5,1,4)",
"(0,1,0)(1,1,0)[12]",
"(1,1,1)(1,1,0)[12]",
"(1,1,1)(1,1,1)[12]"
),
`Q*` = c(
7.1993,
6.7534,
1.8619,
3.2836,
7.8952,
4.2202,
4.9221
),
`df` = c(
12,
10,
3,
3,
11,
9,
8
),
`p-value` = c(
0.8442,
0.7485,
0.6015,
0.3499,
0.7227,
0.8963,
0.7659
),
`Model df` = c(
0,
2,
10,
9,
1,
3,
4
),
`Lags usados` = c(
12,
12,
13,
12,
12,
12,
12
)
)
tabla_check_lag1$`Q*` <- round(tabla_check_lag1$`Q*`, 7)
tabla_check_lag1$`p-value` <- round(tabla_check_lag1$`p-value`, 7)
datatable(
tabla_check_lag1,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:6),
list(width = '180px', targets = 0),
list(width = '150px', targets = 1)
)
),
caption = "Checkresiduals – Prueba de Ljung–Box",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:6,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
textAlign = 'center',
fontWeight = 'bold'
)library(forecast)
library(plotly)
library(dplyr)
modelos <- list(
"Modelo 1" = modelo1_1,
"Modelo 2" = modelo2_1,
"Modelo 3" = modelo3_1,
"Modelo 4" = modelo4_1,
"Modelo 5" = modelo5_sarima12,
"Modelo 6" = modelo6_sarima12,
"Modelo 7" = modelo7_sarima12
)
n_mod <- length(modelos)
trazas_por_mod <- 6
fig_res <- plot_ly()
for (i in seq_along(modelos)) {
md <- modelos[[i]]
nm <- names(modelos)[i]
res <- residuals(md)
fechas <- as.Date(time(res))
v <- as.numeric(res)
visible_i <- ifelse(i == 1, TRUE, FALSE)
fig_res <- fig_res %>%
add_lines(
x = fechas,
y = v,
name = paste(nm, "- residuales"),
xaxis = "x",
yaxis = "y",
visible = visible_i,
showlegend = FALSE,
line = list(color = "#2F3C46", width = 3)
)
acf_obj <- acf(v, plot = FALSE, lag.max = 24, type = "correlation")
ci <- qnorm((1 + 0.95)/2) / sqrt(length(v))
lag <- acf_obj$lag[,1,1][-1]
val <- acf_obj$acf[,1,1][-1]
fig_res <- fig_res %>%
add_bars(
x = lag,
y = val,
name = paste(nm, "- ACF"),
xaxis = "x2",
yaxis = "y2",
visible = visible_i,
showlegend = FALSE,
marker = list(color = "#2F3C46"),
width = 0.4
) %>%
add_lines(
x = c(min(lag), max(lag)),
y = c(ci, ci),
xaxis = "x2",
yaxis = "y2",
visible = visible_i,
showlegend = FALSE,
line = list(color = "#89ACC6", dash = "dash")
) %>%
add_lines(
x = c(min(lag), max(lag)),
y = c(-ci, -ci),
xaxis = "x2",
yaxis = "y2",
visible = visible_i,
showlegend = FALSE,
line = list(color = "#89ACC6", dash = "dash")
)
h <- hist(v, plot = FALSE, breaks = 20)
den <- density(v)
fig_res <- fig_res %>%
add_bars(
x = h$mids,
y = h$density,
name = paste(nm, "- hist"),
xaxis = "x3",
yaxis = "y3",
visible = visible_i,
showlegend = FALSE,
marker = list(color = "#2F3C46"),
width = diff(h$breaks)[1] * 0.9
) %>%
add_lines(
x = den$x,
y = den$y,
name = paste(nm, "- densidad"),
xaxis = "x3",
yaxis = "y3",
visible = visible_i,
showlegend = FALSE,
line = list(color = "#89ACC6", width = 3)
)
}
fig_res <- fig_res %>%
layout(
title = "Diagnóstico de Residuales ARIMA / SARIMA",
xaxis = list(domain = c(0, 1), anchor = "y",
title = "Fecha"),
yaxis = list(domain = c(0.55, 1), title = ""),
xaxis2 = list(domain = c(0, 0.5), anchor = "y2",
title = "Lag"),
yaxis2 = list(domain = c(0, 0.45), title = ""),
xaxis3 = list(domain = c(0.5, 1), anchor = "y3",
title = "Residuales"),
yaxis3 = list(domain = c(0, 0.45), title = ""),
margin = list(t = 110)
)
buttons <- lapply(seq_len(n_mod), function(i){
vis <- rep(FALSE, n_mod * trazas_por_mod)
idx_ini <- (i-1) * trazas_por_mod + 1
idx_fin <- i * trazas_por_mod
vis[idx_ini:idx_fin] <- TRUE
list(
method = "restyle",
args = list("visible", vis),
label = names(modelos)[i]
)
})
fig_res <- fig_res %>%
layout(
updatemenus = list(
list(
type = "dropdown",
x = 0.02,
y = 1.10,
xanchor = "left",
yanchor = "top",
showactive = TRUE,
buttons = buttons,
bgcolor = "rgba(255,255,255,0.8)"
)
)
)
fig_resLos modelos ARIMA (0,1,0) Y ARIMA (1,1,1) demuestran un comportamiento aceptable (valores p 0,84 y 0,74 respectivamente). Sin embargo, estos modelos muestran un Q relativamente alto en comparación con los modelos SARIMA, lo que sugiere que el comportamiento estacional no está siendo capturado. Los modelos de ARIMA más complicados (5,1,5) y (5,1,4) si muestran un valor p aceptable, pero tienen un df de modelo muy alto (9 y 10 parámetros, respectivamente) y requieren más rezagos para la prueba (hasta 13), lo que sugiere un riesgo elevado de sobreajustar el modelo a los patrones específicos del conjunto de entrenamiento y, por consiguiente, perder capacidad de generalización y utilidad debido a mayor complejidad. Con los modelos SARIMA estacionales, (0,1,0)(1,1,0)[12] fue también omitido, en el modo en que aunque su p-value (0.72) está acorde, su Q* Se encuentra en una posición superior a los otros modelos estacionales y su estructura tan simple resulta en patrones estacionales sin modelar. Además, estos resultados reflejan que aunque en varios modelos se cumple que p > 0.05, no todos lo hacen con el mismo nivel de estabilidad, eficiencia y parsimonia, motivo por el cual fueron descartados en favor de mejores alternativas que se adecuarán de manera robusta y consistente con la estacionalidad de la serie.
La SARIMA (1,1,1)(1,1,0)[12] fue la que mostró el mejor desempeño al contar con un balance ideal en el ajuste y la complejidad. Se puebla que no hay residuos autocorrelacionados, a lo que también apunta que su p-valor (0.76) no muestra autocorrelación significativa, algo que también versa sobre el estadístico Q, el cual fue de los más bajos de la serie. El número de parámetros (model df = 4) es el adecuado para poder capturar la dinámica regular y estacional. Además, usa 12 rezagos, lo que es ideal ya que corresponde a la periodicidad mensual de la serie. Este comportamiento es un indicativo de que el modelo captura perfectamente los picos, los valles y los ciclos estacionales de la serie históricamente, y que no es sobreajustado o que no está excesivamente sobrecargado con la estructura matemática. Por todo esto, el modelo presenta residuos que se comportan de forma isotrópica. Además, el modelo logra capturar la dinámica temporal, y la estacionalidad que es evidente dentro de la serie. Desde un punto de vista contextual, esto quiere decir que el modelo tiene la capacidad de reproducir de manera apropiada la tendencia y estacionalidad recurrente en las emisiones, lo que produce pronósticos más fiables y realistas para los períodos futuros. En conclusión, este SARIMA logra describir de manera completa la estructura temporal del fenómeno sin perder generalidad, lo que lo convierte en la opción más sólida para estimar al futuro.
El análisis de precisión es un indicador clave al momento de medir el desempeño predictivo de los modelos ARIMA y SARIMA. Las métricas representan el error que existe al comparar los valores reales y los valores pronosticados, y permiten determinar qué tan bien el modelo es capaz de captar la dinámica de la serie. El ME, RMSE, MAE, MPE, MAPE, MASE y ACF1 miden distintas cualidades de un error. Por ejemplo, de sesgo, de magnitud de error, de autocorrelación, entre otros. En primer lugar, y en la generalidad, el error que arrojan los pronósticos es menor y de mayor distorsión en el sistema y, por lo tanto, posee una mayor capacidad predictiva del sistema. En segundo lugar, un ACF1 que se acerque a cero es una condición deseable, pues se obtiene un comportamiento aleatorio de los residuos, lo que indica que el modelo no ha dejado estructura sin captar.
library(DT)
library(tibble)
library(htmltools)
tabla_accuracy_lag1 <- tibble(
Modelo = c(
"Modelo 1",
"Modelo 2",
"Modelo 3",
"Modelo 4",
"Modelo 5 (Sarima)",
"Modelo 6 (Sarima)",
"Modelo 7 (Sarima)"
),
ME = c(
-2571.158,
-2312.896,
-2379.114,
-2159.185,
6369.583,
5661.034,
7932.631
),
RMSE = c(
111920.5,
110781.4,
102112.2,
99073.52,
115549.5,
107396.6,
83354.29
),
MAE = c(
70759.44,
71051.38,
68951.17,
67474.24,
66049.69,
65573.04,
50599.10
),
MPE = c(
-5.573973,
-4.923518,
-4.205787,
-3.957538,
-2.551123,
-1.849983,
-1.112302
),
MAPE = c(
16.26096,
16.53539,
16.42924,
16.55925,
18.81022,
18.12702,
14.74639
),
MASE = c(
0.1926912,
0.1934863,
0.187767,
0.183745,
0.1798657,
0.1785677,
0.1377908
),
ACF1 = c(
0.1398579,
-0.001598831,
0.001107127,
-0.0183251,
0.3003756,
0.008998942,
-0.00458025
)
)
tabla_accuracy_lag1[,2:8] <- lapply(tabla_accuracy_lag1[,2:8], function(x) round(x, 7))
datatable(
tabla_accuracy_lag1,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:7),
list(width = '200px', targets = 0)
)
),
caption = "Accuracy – Métricas de desempeño de modelos ARIMA",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:7,
fontSize = '13px',
fontFamily = 'Arial, sans-serif',
textAlign = 'center',
fontWeight = 'bold'
)En comparación con los modelos ARIMA, se ha notado que los modelos que poseen estructuras ya sea más simples o más complejas, tienden a tener resultados menos favorables. En caso de los modelos ARIMA (0,1,0), este tiene el RMSE y MAE más alto, por lo que se puede decir que tiene un bajo desempeño en cuanto a la estimación de la variabilidad de la serie de tiempo en cuestión. En el caso de los modelos ARIMA(5,1,5) y ARIMA(5,1,4), aunque si tienen resultados algo mejores, se corren el riesgo de tener un sobreajuste de los datos debido a la cantidad de parámetros de los que se dispone, y por lo tanto, se tiene un MASE más alto y un grado de complejidad que no se necesita en el modelo. En el caso de los modelos SARIMA, el modelo (0,1,0)(1,1,0)[12] también tiene un mayor MAU y MAP y esto sugiere que el no tener un término autorregresivo en la parte no estacional, hace que su modelación sea deficiente.
De todos los modelos analizados, el SARIMA (1,1,1)(1,1,0)[12] se destaca con el mejor rendimiento, debido a que presenta los menores valores de RMSE, MAE, MASE, MAPE y un ACF1 de cero que refleja que no existe autocorrelación residual y de un buen ajuste a la estacionalidad de la serie. Esto quiere decir que el modelo es ágil para captar con efectividad el comportamiento temporal de la serie del CO₂ y las oscilaciones anuales que se presentan de manera cíclica. A partir de este estudio, se traduce en una mayor estabilidad, precisión y confianza para pronosticar tendencias en las emisiones y así, optimizar la planificación y la toma de decisiones en el ámbito energético o ambiental, dentro de un contexto de alta volatilidad y variabilidad estacional.
Se compararon los valores reales de las emisiones de CO₂ para los meses de 2024, con las predicciones para cada uno de los modelos, ARIMA Y SARIMA, con diversas configuraciones de parámetros. Se observa el comportamiento de cada uno de estos modelos y se tiene en cuenta la capacidad de cada modelo de capturar la estacionalidad, la tendencia y el comportamiento de las emisiones en un periodo de tiempo.
library(DT)
library(tibble)
library(dplyr)
library(htmltools)
tabla_comp <- tibble(
Mes = c(
"Mar 2024", "Abr 2024", "May 2024", "Jun 2024",
"Jul 2024", "Ago 2024", "Sep 2024", "Oct 2024",
"Nov 2024", "Dic 2024"
),
Real = c(
1003407.7, 991918.9, 1041493.2, 1072367.6, 1161195.5,
1159170.9, 1077230.6, 1001764.3, 917366.4, 1031531.2
),
Modelo_1 = c(
898700.9, 898700.9, 898700.9, 898700.9,
898700.9, 898700.9, 898700.9, 898700.9,
898700.9, 898700.9
),
Modelo_2 = c(
890186.5, 888588.9, 888289.1, 888232.9,
888222.3, 888220.3, 888220.0, 888219.9,
888219.9, 888219.9
),
Modelo_3 = c(
869807.7, 829279.9, 854575.7, 874333.9,
880808.5, 891637.1, 886616.0, 870697.7,
858944.1, 854849.7
),
Modelo_4 = c(
878651.0, 826500.5, 871652.1, 853472.4,
899426.7, 904126.4, 882502.9, 896075.6,
838114.4, 869212.4
),
Modelo_5 = c(
1024323, 1031430, 1068165, 1101235,
1258462, 1271478, 1208694, 1141659,
1069056, 1173788
),
Modelo_6 = c(
1034589, 1043634, 1081297, 1108427,
1273849, 1288219, 1226130, 1159291,
1084734, 1188088
),
Modelo_7 = c(
927054.4, 787021.2, 826448.6, 854707.6,
960275.8, 989156.3, 954734.7, 934600.2,
886784.1, 971137.0
)
)
datatable(
tabla_comp,
rownames = FALSE,
class = 'cell-border stripe hover compact',
options = list(
dom = 't',
scrollX = TRUE,
paging = FALSE,
ordering = FALSE,
columnDefs = list(
list(className = 'dt-center', targets = 0:8),
list(width = '120px', targets = 0),
list(width = '100px', targets = 1:8),
list(targets = 0:8,
render = JS(
"function(data, type, row, meta) {
return '<div style=\"white-space:nowrap;\">' + data + '</div>';
}"
))
)
),
caption = "Comparación de valores reales y pronósticos",
callback = JS("
$(table.table().header()).css({
'background-color':'#2F4F6A',
'color':'white',
'font-weight':'bold',
'text-align':'center'
});
")
) %>%
formatStyle(
columns = 0:8,
fontSize = '12px',
fontFamily = 'Arial, sans-serif',
fontWeight = 'bold',
textAlign = 'center'
)Los modelos ARIMA (Modelo 1, Modelo 2, Modelo 3 y Modelo 4) tienen limitaciones considerables al tratar de predecir las emisiones de CO₂ de esta serie temporal. La principal razón por la que los modelos no logran predecir con exactitud es que ARIMA no tiene en cuenta la estacionalidad, un componente clave en esta serie.
El Modelo 1 y el Modelo 2 subestiman significativamente las emisiones, particularmente en los meses de máxima actividad como junio, julio y diciembre. Al no capturar la estacionalidad, los modelos muestran una línea de tendencia más plana y no tienen en cuenta adecuadamente los picos de emisiones durante los periodos de alta demanda. Las predicciones en estos meses se desvían significativamente de los valores reales. Los modelos 3 y 4, a pesar de tener un mejor ajuste que los primeros modelos ARIMA, siguen sin captar correctamente la estacionalidad pronunciada de esta serie. Los errores de predicción siguen siendo altos, especialmente en los meses mas estacionalmente variables, como verano e invierno. Por lo tanto, a pesar de su mayor complejidad, los modelos siguen siendo insuficientes para evidenciar la dinámica real de las emisiones del sector.
En contraste de los modelos ARIMA, dentro del grupo SARIMA si se incluye la estacionalidad de manera explícita. El modelo 7, a pesar de que muestra un ajuste adecuado en las pruebas estadísticas, tiene demasiados parámetros estacionales y no estacionales. La gran complejidad del modelo hace que comprenda no sólo la conducta estacional natural, sino también las anomalías, como el pico de 2020 por el COVID. Por lo tanto su predicción se ve afectada por esa distorsión y tiende a generar estimaciones excesivamente adaptadas al pasado o erráticas, lo que hace que pierda su capacidad de predecir. El modelo 5 tiene un mejor rendimiento que el modelo 7, aunque todavía no es el más exacto. Este modelo tiene la capacidad de captar la estacionalidad y el patrón de la serie, lo que permite realizar pronósticos más ajustados a los valores reales, en particular durante los meses con gran actividad aérea. Sin embargo no es el modelo que mejor predice porque tiene ciertas imprecisiones, sobre todo en la segunda mitad del año 2024. Por el contrario, el Modelo 6 hace predicciones más precisas entre todos los modelos analizados. Este modelo puede captar la tendencia y estacionalidad de la serie, lo que permite hacer predicciones mucho más precisas sobre los niveles reales de emisiones de CO₂ en 2024 y a diferencia de otros modelos no se ve afectado por los picos atípicos del 2020 y logra adaptarse a la recuperación post-pandemia.
library(plotly)
library(forecast)
library(dplyr)
pronostico1 <- forecast(modelo1_1, h = 9)
pronostico2 <- forecast(modelo2_1, h = 9)
pronostico3 <- forecast(modelo3_1, h = 9)
pronostico4 <- forecast(modelo4_1, h = 9)
pronostico5 <- forecast(modelo5_sarima12, h = 9)
pronostico6 <- forecast(modelo6_sarima12, h = 9)
pronostico7 <- forecast(modelo7_sarima12, h = 9)
ventana2_ts <- ts(as.numeric(ventana2),
start = c(2024, 3),
frequency = 12)
df_real <- data.frame(
Fecha = as.Date(time(ventana2_ts)),
Valor = as.numeric(ventana2_ts)
)
convert_fcast <- function(pr){
data.frame(
Fecha = as.Date(time(pr$mean)),
Media = as.numeric(pr$mean),
Lower = as.numeric(pr$lower[,2]),
Upper = as.numeric(pr$upper[,2])
)
}
fc_list <- list(
"Modelo 1" = convert_fcast(pronostico1),
"Modelo 2" = convert_fcast(pronostico2),
"Modelo 3" = convert_fcast(pronostico3),
"Modelo 4" = convert_fcast(pronostico4),
"Modelo 5" = convert_fcast(pronostico5),
"Modelo 6" = convert_fcast(pronostico6),
"Modelo 7" = convert_fcast(pronostico7)
)
base <- fc_list[[1]]
fig_pronostico <- plot_ly()
fig_pronostico <- fig_pronostico %>% add_lines(
data = df_real,
x = ~Fecha, y = ~Valor,
name = "Datos reales",
line = list(color = "#2F4F6A", width = 3),
legendgroup = "real",
showlegend = TRUE
)
fig_pronostico <- fig_pronostico %>% add_lines(
x = base$Fecha, y = base$Lower,
name = "IC 95% (lower)",
line = list(color = "rgba(0,0,0,0)"),
showlegend = FALSE,
legendgroup = "ic"
)
fig_pronostico <- fig_pronostico %>% add_lines(
x = base$Fecha, y = base$Upper,
name = "Intervalo 95%",
fill = "tonexty",
fillcolor = "rgba(137,172,198,0.25)",
line = list(color = "rgba(0,0,0,0)"),
showlegend = TRUE,
legendgroup = "ic"
)
fig_pronostico <- fig_pronostico %>% add_lines(
x = base$Fecha, y = base$Media,
name = "Pronóstico",
line = list(color = "#89ACC6", width = 3),
showlegend = TRUE,
legendgroup = "forecast"
)
buttons <- lapply(seq_along(fc_list), function(i){
df <- fc_list[[i]]
list(
method = "restyle",
args = list(
list(
x = list(
df_real$Fecha,
df$Fecha,
df$Fecha,
df$Fecha
),
y = list(
df_real$Valor,
df$Lower,
df$Upper,
df$Media
)
)
),
label = names(fc_list)[i]
)
})
fig_final <- fig_pronostico %>%
layout(
title = list(
text = "Pronóstico vs Datos reales",
x = 0.5,
font = list(size = 22, color = "#2F4F6A")
),
updatemenus = list(
list(
type = "dropdown",
x = 0.04, y = 1.12,
showactive = TRUE,
buttons = buttons
)
),
xaxis = list(title = "Tiempo"),
yaxis = list(title = "Valor"),
plot_bgcolor = "white",
paper_bgcolor = "white",
legend = list(
x = 0.02, y = 0.98,
bgcolor = "rgba(255,255,255,0.8)"
)
)
fig_finalCon base en el análisis anterior se escogió el modelo SARIMA(1,1,1)(1,1,0)[12] como la opción más adecuada para describir el comportamiento de la serie. Este modelo no solo muestra un mayor acercamiento a los valores reales en la perspectiva proyectada, sino que también enfatiza en los criterios técnicos utilizados para su validación tiene uno de los valores AIC más bajos, lo cual demuestra un equilibrio eficaz entre ajuste y complejidad, presenta métricas precisas con errores más pequeños y estables, lo que refleja una mejor capacidad de predicción y supera de manera adecuada la prueba Checkresiduals, mostrando residuos que se comportan como ruido blanco. Esto significa que logra capturar correctamente la dinámica temporal sin dejar patrones estructurados sin modelar.
Aunque este modelo muestra un gran rendimiento en comparación con las otras especificaciones analizadas, sus predicciones no coinciden a la perfección con todos los valores observados en la ventana de prueba de 2024. Esto es esperable en modelos de series temporales, puesto que su propósito no consiste en reproducir cada observación individualmente, sino en captar de manera coherente la estructura subyacente del proceso que genera los datos. El modelo, a pesar de que es capaz de replicar el nivel general de emisiones, la dirección de la tendencia y el patrón estacional anual en lo que concierne a las emisiones de CO₂ en la aviación internacional canadiense, presenta diferencias en algunos meses en los que las emisiones reales superan el pronóstico basado en el promedio histórico.
Al comparar los valores reales y pronosticados de un mes a otro, se observa que el modelo sigue adecuadamente los picos y valles de la serie, es decir, las predicciones tienden a estar en el orden de magnitud adecuado y a mantener la forma general del ciclo anual, sobre todo durante los meses en que la estacionalidad habitual determina el comportamiento, como invierno o verano. Sin embargo, las diferencias entre los datos estimados y observados incrementan en momentos concretos, especialmente cuando ocurren aumentos o descensos más significativos que los registrados en años anteriores.
library(plotly)
library(forecast)
df_full <- data.frame(
Fecha = co2_aviones$Fecha,
Valor = co2_aviones$CO2
)
pronostico_full6 <- forecast(modelo6_sarima12, h = 9)
df_fcast6 <- data.frame(
Fecha = as.Date(time(pronostico_full6$mean)),
Media = as.numeric(pronostico_full6$mean),
Lower = as.numeric(pronostico_full6$lower[,2]),
Upper = as.numeric(pronostico_full6$upper[,2])
)
fig_full6 <- plot_ly()
fig_full6 <- fig_full6 %>% add_lines(
data = df_full,
x = ~Fecha, y = ~Valor,
name = "Serie completa",
line = list(color="#2F4F6A", width=3)
)
fig_full6 <- fig_full6 %>% add_ribbons(
data = df_fcast6,
x = ~Fecha,
ymin = ~Lower,
ymax = ~Upper,
name = "IC 95%",
fillcolor = "rgba(137,172,198,0.25)",
line = list(color="transparent")
)
fig_full6 <- fig_full6 %>% add_lines(
data = df_fcast6,
x = ~Fecha, y = ~Media,
name = "Pronóstico Modelo 6",
line = list(color="#89ACC6", width=3)
)
fig_full6 <- fig_full6 %>% layout(
title = list(
text = "Serie Completa + Pronóstico (Modelo 6)",
x = 0.5,
font = list(size=22, color="#2F4F6A")
),
xaxis = list(title = "Fecha"),
yaxis = list(title = "CO₂"),
plot_bgcolor = "white",
paper_bgcolor = "white",
legend = list(x=0.02, y=0.98),
margin = list(t = 90)
)
fig_full6Teniendo en cuenta el contexto canadiense, estas diferencias tienen aún más sentido. La aviación en Canadá se ve fuertemente afectada por elementos climáticos y geográficos, la dependencia con conexiones hacia Europa, Estados Unidos y lugares cálidos, así como una notable sensibilidad estacional invernal. A parte de la evidente estacionalidad anual, con veranos de gran actividad y finales de año con picos relacionados con celebraciones, este ámbito continúa vulnerable a impactos externos como olas de calor, fuertes tormentas de nieve, incendios forestales a gran escala o alteraciones en la percepción del riesgo de viajar. Muchos de estos acontecimientos no siguen un patrón estrictamente repetitivo y, en consecuencia, un modelo que se fundamenta en la repetición histórica de la serie no puede captarlos en su totalidad. Que el modelo 6 no se adecue de manera óptima en algunos meses no implica que esté mal especificado, sino que la realidad del sector económico y climático incluye factores de variabilidad que superan lo que una estructura estadística puede captar.
Un modelo que replicara con precisión cada cambio en la serie probablemente estaría sobreajustado, reproduciendo no solamente la estacionalidad y la tendencia, sino también los efectos particulares de años concretos como el 2020, que no reflejan un comportamiento normal en el sistema. Este tipo de sobreajuste disminuye la capacidad del modelo para generalizar y generar predicciones útiles en el futuro. Este modelo no se centra en analizar cada pequeña variación mensual, sino que se ajusta para mantener una estructura flexible que represente adecuadamente la dinámica interna de las emisiones aéreas en Canadá. Esto es especialmente importante en un país donde la demanda de transporte aéreo es vulnerable ante modificaciones base, como el desarrollo de combustibles sostenibles para la aviación, el endurecimiento de políticas climáticas o la recuperación post-pandemia.
El hecho que el modelo 6 no replique con exactitud todos los datos de la ventana de prueba no demerita su calidad, por el contrario, confirma que es un modelo realista que capta la estructura fundamental. Esto implica que el modelo es lo suficientemente sensible para mostrar el patrón de emisiones de CO₂ y la tendencia de recuperación del sector en un año, pero también lo suficientemente cuidadoso como para no replicar eventos inusuales como la pandemia del 2020. Desde el punto de vista de la planificación climática y del monitoreo de los objetivos de reducción de dióxido de carbono, contar con un modelo que represente con precisión la dinámica general, adapte un margen de error mensual y mantenga una buena capacidad para generalizar es más valioso que uno que se ajuste perfectamente a los acontecimientos recientes pero no aplica en el momento en que las condiciones del sistema cambien.
El análisis de la serie temporal de emisiones de CO₂ del sector aéreo canadiense permitió identificar una dinámica compleja, marcada por la estacionalidad anual, tendencias estructurales y choques excepcionales como la caída abrupta de 2020 asociada a la pandemia. A lo largo de los distintos modelos evaluados, se observó que la estacionalidad constituye uno de los elementos más determinantes del comportamiento de la serie, lo que hace que los modelos ARIMA tradicionales resulten insuficientes para reproducirla adecuadamente. En cambio, los modelos SARIMA lograron capturar con mayor precisión tanto la variabilidad mensual como los ciclos estacionales característicos del sector, destacándose el SARIMA(1,1,1)(1,1,0)[12] como la alternativa más robusta, estable y equilibrada entre complejidad y capacidad predictiva.
Los pronósticos obtenidos con este modelo ofrecen una representación coherente de la evolución esperada de las emisiones en el corto plazo, reproduciendo la recuperación post-pandemia y los picos estacionales propios de la aviación canadiense. Esto convierte al modelo en una herramienta valiosa para anticipar escenarios futuros y apoyar la formulación de políticas públicas orientadas a la mitigación climática, la planeación energética, y el monitoreo de compromisos internacionales como CORSIA. La capacidad del modelo para capturar patrones persistentes sin replicar anomalías excepcionales, como el colapso de la actividad aérea en 2020, refuerza su utilidad en el diseño de estrategias de mitigación flexibles, que respondan a la estructura real del sistema y no a eventos irrepetibles.
Asimismo, el análisis evidenció la importancia de contar con modelos capaces de predecir de manera confiable las emisiones futuras, dado que la efectividad de las políticas climáticas depende en gran medida de la precisión de estas proyecciones. En un sector influido por factores tecnológicos, económicos y regulatorios a escala global, los pronósticos brindados por un modelo como el SARIMA permiten anticipar cambios, ajustar políticas y planificar intervenciones de forma oportuna y sostenible. Su flexibilidad favorece que las políticas públicas se adapten a las variaciones del sector aeronáutico, en lugar de reaccionar únicamente a los cambios una vez ocurridos.
No obstante, el modelo presenta ciertas limitaciones, ya que sus predicciones pueden verse afectadas por factores externos no incluidos formalmente, tales como innovaciones tecnológicas, fluctuaciones macroeconómicas, eventos climáticos extremos o cambios inesperados en la demanda turística. Aun así, al centrarse en los patrones estructurales y estacionales recurrentes, el modelo proporciona una base sólida para la vigilancia y proyección de emisiones, permitiendo una comprensión más profunda del fenómeno y apoyando la toma de decisiones en un contexto global dinámico y comprometido con la descarbonización.
En conjunto, los resultados muestran que el uso de modelos SARIMA no solo es adecuado para describir y proyectar las emisiones del sector aéreo canadiense, sino que también constituye un instrumento fundamental para la formulación de políticas climáticas efectivas, adaptativas y alineadas con los objetivos de sostenibilidad a corto, mediano y largo plazo.
Barandica, J. (s.f.). Joaqui Barandica - Recursos de programación y análisis de datos. Recuperado de https://www.joaquibarandica.com/
Canada, E. A. C. C. (2019, 2 de abril). Canada’s climate is warming twice as fast as global average. Recuperado de https://www.canada.ca/en/environment-climate-change/news/2019/04/canadas-climate-is-warming-twice-as-fast-as-global-average.html
Datosmacro. (s. f.). Canadá – Emisiones de CO₂ 2023. Recuperado de https://datosmacro.expansion.com/energia-y-medio-ambiente/emisiones-co2/canada
Transport Canada. (2022). Canada’s Aviation Climate Action Plan 2022–2030. Recuperado de https://tc.canada.ca/sites/default/files/2022-11/canada-aviation-climate-action-plan-2022-2030.pdf
Ministerio para la Transición Ecológica y Solidaria de Francia. (2019). Acuerdo entre Canadá y Francia sobre cooperación en cambio climático y energía limpia. Recuperado de https://www.ecologie.gouv.fr/accords-climatiques-canada-france
Joaqui-Barandica, O. (2024). Gestión de Datos: Series de tiempo y pronóstico. Universidad del Valle.
Xie, Y. (2024). R Markdown: The Definitive Guide. Bookdown. Recuperado de https://bookdown.org/yihui/rmarkdown/html-document.html
Ukpanah, I. (2024). Does Air Travel Significantly Impact the Environment? A Comprehensive Analysis. GreenMatch Blog. Recuperado de https://www.greenmatch.co.uk/blog/aviation-environmental-impact
Kiracı, K. (2025). Determinants of air transport CO₂ emissions in OECD countries. Sustainable Futures, 10, 100877. https://doi.org/10.1016/j.sftr.2025.100877
Ritchie, H. (2020). Where in the world do people have the highest CO₂ emissions from flying? Our World in Data. Recuperado de https://archive.ourworldindata.org/20251125-173858/carbon-footprint-flying.html
Zheng, S. (2024). Demand response to aviation carbon pricing in Canada (Working Paper). International Council on Clean Transportation (ICCT).
Transport Canada. (2022). Canada’s aviation climate action plan 2022–2030. Recuperado de https://www.tc.gc.ca/aviation-emissions/
Perlman, K. (2018). Contribución del sector de la aviación mundial al logro de los objetivos climáticos del Acuerdo de París.
Naciones Unidas. (2020). Las emisiones de CO₂ rompen otro récord: un calentamiento global catastrófico amenaza el planeta. Noticias ONU.