PROCESO DE BOX - JENKINS
Jennifer Dayana Calderón Fonseca
Mayo 2024
1. INTRODUCCIÓN
ANOMALÍAS EN LA TEMPERATURA
Las anomalías en la temperatura son una herramienta fundamental en el análisis climático y meteorológico, permitiendo una comprensión más precisa de los cambios y de las tendencias a lo largo del tiempo. Desde una perspectiva estadística, las anomalías en la temperatura se refieren a las desviaciones de las temperaturas observadas respecto a un valor de referencia, generalmente media histórica o climática.
Para calcular una anomalía de temperatura, primero se determina un período base o de referencia. Este período base es típicamente un intervalo de tiempo largo, que se considera representativo dentro de las condiciones climáticas normales. La temperatura media de este período se utiliza como referencia.
1.1. Contexto
Colombia, es un país situado en la franja ecuatorial, por ello expermienta una variedad de climas debido a su geografía diversa que incluye costas, montañas y selvas tropicales. El análisis de anomalías climáticas en este contexto es crucial para entender los efectos del cambio climático, identificar patrones de variabilidad climática y desarrollar estrategias de mitigación y adaptación. Las anomalías de temperatura en particular, indican desviaciones respecto a una media climática de referencia y son indicadoras clave del calentamiento global y sus impactos locales.
El presente estudio, busca analizar un conjunto de datos que registra las anomalías de temperatura mensuales en Colombia desde el mes de Julio en 2015 hasta el mes de septiembre de 2023. Este análisis es esencial para identificar tendencias a largo plazo, variaciones estacionales y posibles eventos climáticos extremos que afectan a diversas regiones del país.
1.2. Objetivo
El objetivo principal de este proyecto es, realizar el proceso completo de Box-Jenkins para la identificación, estimación, verificación y uso de un modelo SARIMA, a través de un conjunto de datos reales sobre anomalías de temperaturas en Colombia. Este enfoque permitirá entender mejor las dinámicas de las anomalías de temperatura y generar pronósticos precisos basados en el modelo SARIMA ajustado.
Se pretende determinar la tendencia general de las anomalías de tempratura identificando si existe un aumento o una disminución significativa a largo plazo. Asi mismo, identificar variaciones estacioanales y mensuales donde se aprecie a simple vista cómo varían las anomalías entre los diferentes meses y estaciones del año.
1.3. Metodología
Para alcanzar el objetivo propuesto se recopilará y se preparará adecuadamente los datos seleccionados, que pueden ser encontrados en https://www.copernicus.eu/es/node/43003
* Se utilizarán datos de anomalías de temperatura mensuales registrados desde Julio de 2015 hasta Septiembre de 2023.
* Se generarán gráficos de la serie temporal original de las anomalías de temperatura.
* Se realizará una descripción de los datos, identificando la presencia de tendencias, estacionalidad y heterocedasticidad.
* Se aplicarán transformaciones necesarias, como diferenciación para lograr la estacionariedad de la serie.
* Se evaluará la estacionariedad mediante gráficos de autocorrelación y pruebas estadisticas.
* Se identificarán por lo menos dos posibles modelos SARIMA utilizando gráficos de autocorrelación FAC & FACP.
* Se estimarán los parámetros de los modelos identificados utilizando métodos de máxima verosimilitud.
* Se compararán los modelos propuestos utilizando criterios de información como AIC, BIC, RMSE, & MAPE.
* Se analizarán los residuos de los modelos seleccionados para verificar los supuestos de dependencia y de normalidad.
* Se estimarán los parámetros del modelo final seleccionado se presentará utilizando la notación adecuada.
* Se realizarán pronósticos de la serie temporal original.
* Se elaborará un resumen de los hallazgos, los problemas encontrados durante el ajuste del modelo y las razones por las cuales el modelo final es el más adecuado.
2. IDENTIFICACIÓN DEL MODELO.
2.1. Carga de información y variables a considerar
Las siguientes líneas de código permiten cargar las librerías necesarias para el análisis que se va a realizar y así poder visualizar los datos.
library(readr)## Warning: package 'readr' was built under R version 4.2.3cli_colombia <- read_csv("CLIMA.CSV")## Rows: 105 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): Year, Anomaly
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.ts_datos<-ts(cli_colombia$Anomaly, start=c(2015,1),
frequency=12)
ts_datos## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2015 1.20 0.99 1.17 1.46 1.35 1.64 1.76 2.03 1.90 1.57 1.41 1.46
## 2016 1.72 1.85 1.41 1.53 0.88 0.48 1.42 1.64 0.94 1.23 0.98 1.09
## 2017 0.93 0.93 1.08 0.96 1.45 1.42 1.39 1.88 1.69 1.37 1.02 1.16
## 2018 0.81 0.63 1.07 1.36 1.04 0.41 0.90 0.57 1.66 1.38 1.20 0.92
## 2019 1.27 1.14 0.97 1.50 1.39 1.94 1.16 1.04 1.70 1.25 1.72 1.23
## 2020 1.40 1.16 1.86 1.35 1.01 1.62 1.27 1.49 1.83 1.79 1.43 1.16
## 2021 0.90 0.85 0.87 1.16 0.85 0.92 1.06 1.57 2.09 1.48 1.00 1.15
## 2022 1.39 1.05 0.73 0.90 0.25 0.44 1.87 0.93 0.88 1.02 0.78 0.99
## 2023 0.84 1.01 1.37 1.10 1.55 1.51 2.12 2.35 2.482.2. Gráfico de la serie original
plot(ts_datos, main="ANOMALÍAS CLIMÁTICAS EN COLOMBIA (2015 - 2023)",
xlab="Mes", ylab="Temperatura")
El anterior gráfico, muestra las anomalías climáticas en Colombia desde el año 2015 hasta 2023, modeladas mediante un modelo SARIMA.
2.3. Análisis de autocorrelación
La función de autocorrelación acf proporciona la autocorrelación en todos los retrasos posibles. La autocorrelación en el retraso 0 se incluye por defecto, que siempre toma el valor 1, ya que representa la correlación entre los datos y ellos mismos.
acf(ts_datos, main="Serie Original : Función de Autocorrelación Simple")
La función de autocorrelación pacf, en el rezago k es una función que describe la correlación entre todos los puntos de datos que son exactamente k pasos separados, después de tener en cuenta su correlación con los datos entre esos k pasos. Ayuda a identificar el número de coeficientes de autorregresión (valor p) en un modelo SARIMA.
pacf(ts_datos, main="Serie Original : Función de Autocorrelación Parcial")
3. DESCOMOPOSICIÓN DE LA SERIE.
3.1. Gráficos
ts_datos_dec<-decompose(ts_datos)
par(mfrow=c(2,1))
plot(ts_datos_dec$x,main="Serie original",col="black",ylab="Serie
de Tiempo")
plot(ts_datos_dec$trend, main = "Tendencia de la serie", col = "blue",
ylab = "Valores")
plot(ts_datos_dec$seasonal, main = "Estacionalidad de la serie", col =
"red", ylab = "Valores")
plot(ts_datos_dec$random, main = "Heterocedasticidad de la serie", col =
"green", ylab = "Valores")
plot(decompose(ts_datos))
Observando el gráfico, no parece haber una tendencia clara y persistente en el tiempo. Las anomalías climáticas fluctúan entorno a un rango constante sin mostrar un aumento o una disminución sostenida. Lo que sugiere, que a lo largo del periodo analizado, las variaciones en las temperaturas anómalas no siguen una dirección única y continua.
Del mismo modo, los datos presentan cierta estacionalidad, que se evidencia por los picos y valles que se repiten en algunos intervalos regulares. Estos patrones ciclícos indican que hay factores estacionales que afectan las temperaturas de manera recurrente cada cierto tiempo. Hay que aclarar que este tipo de comportamientos son bastante común en datos climáticos, donde las variaciones estacionales influyen en las temperaturas.
Finalmente, observamos que en cuanto a la heterocedasticidad, que es la que se refiere a la variabilidad cambiante de los datos a lo largo del tiempo, el gráfico nos muestra fluctuaciones en la amplitud de las anomalías. Durante ciertos periodos, las desviaciones de la temperatura parecen ser un poco más amplias o un poco más estrechas que en otros periodos. Esto indica que hay evidencia de variabilidad en la magnitud de las anomalias climáticas en Colombia, para este periodo de tiempo. Por lo tanto, esto sugiere la presencua de heterocedasticidad en los datos observados.
4. DESESTACIONALIZACIÓN DE LA SERIE.
library(forecast)
ts_datos_dec_SA <- seasadj(ts_datos_dec)
ts_datos_dec_SA## Jan Feb Mar Apr May Jun Jul
## 2015 1.2752176 1.1487593 1.2409468 1.4239974 1.5769736 1.8089974 1.6396968
## 2016 1.7952176 2.0087593 1.4809468 1.4939974 1.1069736 0.6489974 1.2996968
## 2017 1.0052176 1.0887593 1.1509468 0.9239974 1.6769736 1.5889974 1.2696968
## 2018 0.8852176 0.7887593 1.1409468 1.3239974 1.2669736 0.5789974 0.7796968
## 2019 1.3452176 1.2987593 1.0409468 1.4639974 1.6169736 2.1089974 1.0396968
## 2020 1.4752176 1.3187593 1.9309468 1.3139974 1.2369736 1.7889974 1.1496968
## 2021 0.9752176 1.0087593 0.9409468 1.1239974 1.0769736 1.0889974 0.9396968
## 2022 1.4652176 1.2087593 0.8009468 0.8639974 0.4769736 0.6089974 1.7496968
## 2023 0.9152176 1.1687593 1.4409468 1.0639974 1.7769736 1.6789974 1.9996968
## Aug Sep Oct Nov Dec
## 2015 1.8679260 1.5465718 1.4157385 1.4486551 1.5465197
## 2016 1.4779260 0.5865718 1.0757385 1.0186551 1.1765197
## 2017 1.7179260 1.3365718 1.2157385 1.0586551 1.2465197
## 2018 0.4079260 1.3065718 1.2257385 1.2386551 1.0065197
## 2019 0.8779260 1.3465718 1.0957385 1.7586551 1.3165197
## 2020 1.3279260 1.4765718 1.6357385 1.4686551 1.2465197
## 2021 1.4079260 1.7365718 1.3257385 1.0386551 1.2365197
## 2022 0.7679260 0.5265718 0.8657385 0.8186551 1.0765197
## 2023 2.1879260 2.1265718Al desestacionalizar la serie de datos, se eliminó el componente estacional de la serie, lo que permite analizar las tendencias y las fluctuaciones subyacentes sin la influencia de patrones estacionales recurrentes.
Por lo tanto, la serie ajustada ya no presenta los patrones regulares repetitivos que se observan tipicamente en los datos estacionales.
En 2023, se observa un notable aumento en la anomalías desestacionalizadas, especialmente en los meses de Agosto y Septiembre, con valores de 2.1879260 y 2.1265718 respectivamente, lo que indica eventos climáticos inusuales o extremos durante estos meses.
De esta manera, en el año 2018 se muestran valores más bajos en general, como en el mes de Junio 0.5789974 y Agosto 0.4079260, lo que podría sugerir un periodo de menor variabilidad climática.
4.1 Gráfico de la serie original y desestacionalizada
plot(ts_datos, main="Serie Original VS. Serie Desestacionalizada", col="blue")
lines(ts_datos_dec_SA, col="red")
Se puede observar que, la desestacionalización ha permitido obtener una visión un poco más clara de las tendencias y de las variaciones, libre de los efectos de la estacionalidad.
Esto facilita la identificación de patrones anómalos y la evaluación de tendencias a largo plazo, lo que es crucial para el desarrollo de modelos predictivos.
5. TEST DE DIKEY FULLER AUMENTADO (ADF)
Hipótesis:
\(H_0:\) La serie no es estacionaria, tiene raíz unitaria.
\(H_1:\) La serie es estacionaria, no tiene raíz unitaria.
5.1. Estadístico de prueba
library(tseries)
adf.test(ts_datos)##
## Augmented Dickey-Fuller Test
##
## data: ts_datos
## Dickey-Fuller = -2.426, Lag order = 4, p-value = 0.4
## alternative hypothesis: stationaryEl valor del estadistico es -2.426, con un p.valor de 0.4. Como este p valor es mayor que un nivel de significancia de 0.05, no hay suficiente evidencia estadística para rechazar la hipótesis nula. Por lo tanto, no se puede concluir que la serie de tiempo original sea estacionaria en este nivel de significancia.
En definitiva, la serie de anomalías climáticas en Colombia , tal como está, no es estacionaria.
5.2. Implicaciones para el MODELO SARIMA
La no estacionariedad de la serie indica que se debe realizar transformaciones adicionales para estabilizar la media y/o la varianza antes de ajustar un modelo SARIMA.
Entonces se realizará la diferenciación, que será aplicar una o más diferencias a la serie para eliminar tendencias y convertir la serie en estacionaria. Y la desestacionalización, que ya fue aplicada en pasos anteriores, que es la que elimina el componente estacional de la serie.
6. DIFERENCIAS
library(forecast)
ndiffs(ts_datos) # Número de diferencias necesarias para estacionalizar la serie ## [1] 0nsdiffs(ts_datos) # Número de diferencias estacionales necesarias. ## [1] 0plot(ts_datos,type="o",
main="Serie de tiempo Original",
xlab="Año",
ylab="Anomalías")
plot(diff(ts_datos),type="o",
main="Serie de tiempo Diferenciada",
xlab="Año",
ylab="Primera diferencia")
plot(diff(ts_datos,12),type="o",
main="Serie de tiempo con diferencia Estacional",
xlab="Año",
ylab="Diferencia Estacional")
acf(ts_datos, lag.max=40,
main="Función de Autocorrelación (ACF), Serie Original",
xlab="",
ylab="ACF")
pacf(ts_datos, lag.max=40,
main="Función de Autocorrelación Parcial (ACF), Serie Original",
xlab="",
ylab="PACF")
acf(diff(ts_datos), lag.max=40,
main="Función de Autocorrelación (ACF), Serie Diferenciada",
xlab="",
ylab="ACF")
pacf(diff(ts_datos), lag.max=40,
main="Función de Autocorrelación Parcial (PACF), Serie Diferenciada",
xlab="",
ylab="ACF")
acf(diff(ts_datos,12), lag.max=100,
main="Función de Autocorrelación (ACF), Serie con diferencia Estacional
",
xlab="",
ylab="ACF")
pacf(diff(ts_datos,12), lag.max=40,
main="Función de Autocorrelación (ACF), Serie con diferencia Estacional",
xlab="",
ylab="ACF")
Así mismo, el correlograma confirma un declinamiento rápido en el tiempo de la FAC Y FACP.
6.1. Prueba de Dickey - Fuller a la serie en primera diferencia
Hipótesis:
\(H_0:\) La serie no es estacionaria, tiene raíz unitaria.
\(H_1:\) La serie es estacionaria, no tiene ráiz unitaria.
library(tseries)
adf.test(diff(ts_datos)) # Prueba de Dickey-Fuller en la Serie de Tiempo Diferenciada.## Warning in adf.test(diff(ts_datos)): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff(ts_datos)
## Dickey-Fuller = -5.7588, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationarydiff(ts_datos,12)# Serie de Tiempo con Diferencia Estacional.## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2016 0.52 0.86 0.24 0.07 -0.47 -1.16 -0.34 -0.39 -0.96 -0.34 -0.43 -0.37
## 2017 -0.79 -0.92 -0.33 -0.57 0.57 0.94 -0.03 0.24 0.75 0.14 0.04 0.07
## 2018 -0.12 -0.30 -0.01 0.40 -0.41 -1.01 -0.49 -1.31 -0.03 0.01 0.18 -0.24
## 2019 0.46 0.51 -0.10 0.14 0.35 1.53 0.26 0.47 0.04 -0.13 0.52 0.31
## 2020 0.13 0.02 0.89 -0.15 -0.38 -0.32 0.11 0.45 0.13 0.54 -0.29 -0.07
## 2021 -0.50 -0.31 -0.99 -0.19 -0.16 -0.70 -0.21 0.08 0.26 -0.31 -0.43 -0.01
## 2022 0.49 0.20 -0.14 -0.26 -0.60 -0.48 0.81 -0.64 -1.21 -0.46 -0.22 -0.16
## 2023 -0.55 -0.04 0.64 0.20 1.30 1.07 0.25 1.42 1.60ts_cor(diff(ts_datos), lag.max=40) # Gráfico de Correlación (ACF y PACF) de la Serie de Tiempo Diferenciada. ts_cor(ts_datos, lag.max=40) # Gráfico de Correlación (ACF y PACF) de la Serie de Tiempo Original.El valor del estadístico es de -5.7588, con un valor p de 0.01, Como este valor (p) es menor que el nivel de significancia (0.05), indica que hay suficiente evidencia estadística para rechazar la hipótesis nula.
Por lo tanto, se puede concluir que la serie diferenciada es estacionaria en este nivel de significancia. Es decir, la aplicación de la diferenciación ha logrado hacer que la serie de anomalías climáticas en Colombia sea estacionaria.
En los correlogramas podemos observar que la gráfica FAC muestra la autocorrelación en la serie temporal diferenciada y se observa una clara autocorrelación significativa en los retardos no estacionales 1, 2, 3 y 12. Esto indica que los valores actuales de la serie temporal están correlacionados con los valores de la serie en los periodos 1, 2, 3 y 12 periodos atrás. La autocorrelación se va desvaneciendo con el tiempo, por lo tanto, se afirma que la serie temporal es estacionaria.
La gráfica FACP muestra la autocorrelación parcial en la serie temporal diferenciada. Se observa una clara autocorrelación parcial significativa en los retardos no estacionales 1, 2 y 12. Esto indica que los valores actuales de la serie temporal están correlacionados con los valores de la serie en los periodos 1, 2 y 12 periodos atrás, incluso después de controlar por la autocorrelación en periodos anteriores.Después se va desvaneciendo con el tiempo, y por lo tanto, la serie temporal es estacionaria.
Esto nos lleva a pensar de un posible modelo SARIMA(3, 1, 1)(0, 1, 1)12 Explicación del modelo:
(3, 1, 1): Este componente ARIMA(3, 1, 1) modela la autocorrelación no estacional en la serie temporal diferenciada.
7. AJUSTE DEL MODELO
7.1 Modelo propuesto en R Studio
Se implementa un Modelo Autorregresivo Integrado de Promedio Movil Estacional (SARIMA), para la serie temporal de anamolías climáticas en Colombia. Este modelo tiene en cuenta tanto la estructura autoregresiva (AR) como la estructura de promedio móvil (MA), así como los componentes de estacionalidad.
modelo<-auto.arima(ts_datos,
stepwise=FALSE,approximation=FALSE)
modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07Tenemos un modelo ARIMA(1,0,0) con un omponente autoregresivo de orden 1 (AR) y sin componentes de diferenciación ni promedio móvil. El modelo incluye un término de media no nula con un valor estimado de 1.2830.
Coeficiente AR1: 0.5263 Término de Media: 1.2830 Errores Estándar (s.e.): 0.0868 para AR1 y 0.0726 para la media.
summary(modelo)## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0005361043 0.3554915 0.2871918 -11.0671 27.86595 0.647644
## ACF1
## Training set -0.02067759modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07coeftest(modelo)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.526304 0.086811 6.0627 1.339e-09 ***
## intercept 1.283022 0.072586 17.6758 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07confint(modelo)## 2.5 % 97.5 %
## ar1 0.3561582 0.6964503
## intercept 1.1407554 1.4252888modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07El modelo SARIMA ajustado es de la forma ARIMA(1,0,0) con un término de media no nula. Esto indica que el modelo tiene en cuenta la dependencia temporal a través de un componente autoregresivo de orden 1 y un término de media constante. Los coeficientes estimados sugieren una asociación positiva entre los valores pasados y futuros de la serie de anomalías climáticas, así como un nivel medio de 1.2830.
7.2. Prueba de Ljung-Box (Evaluación del Modelo)
error<-residuals(modelo)
Box.test(error,type="Ljung-Box")##
## Box-Ljung test
##
## data: error
## X-squared = 0.046189, df = 1, p-value = 0.8298plot(error)
La prueba de Ljung-Box se ha aplicado a los residuos del modelo SARIMA para evaluar si hay autocorrelación significativa.
El valor del estadístico X-squared es 0.046189 con 1 grado de libertad.
El valor p de 0.8298 indica que no hay suficiente evidencia para rechazar la hipótesis nula de que no hay autocorrelación en los residuos del modelo SARIMA. Esto sugiere que los residuos no muestran autocorrelación significativa después de haber sido ajustados por el modelo.
8. PRÓNOSTICOS
pronostico<-forecast::forecast(modelo,h=5)
pronostico ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Oct 2023 1.912997 1.4530141 2.372979 1.2095141 2.616479
## Nov 2023 1.614580 1.0947805 2.134380 0.8196151 2.409546
## Dec 2023 1.457523 0.9223348 1.992710 0.6390234 2.276022
## Jan 2024 1.374862 0.8354898 1.914235 0.5499632 2.199762
## Feb 2024 1.331358 0.7908320 1.871884 0.5046948 2.158021* Octubre 2023: Intervalo de confianza (1.453, 2.373), Pronóstico: 1.913 * Noviembre 2023: Intervalo de confianza (1.095, 2.134), Pronóstico: 1.615 * Diciembre 2023: Intervalo de confianza (0.922, 1.993), Pronóstico: 1.458 * Enero 2024: Intervalo de confianza (0.835, 1.914), Pronóstico: 1.375 * Febrero 2024: Intervalo de confianza (0.791, 1.872), Pronóstico: 1.331
Estos resultados proporcionan estimaciones de la anomalía climática esperada para los próximos meses, junto con intervalos de confianza que indican la incertidumbre asociada con cada pronóstico. Es importante tener en cuenta estos intervalos para comprender la variabilidad esperada en las predicciones.
plot(pronostico)
9. MODELOS PROPUESTOS
9.1 Modelo 1
ts_cor(diff(ts_datos), lag.max=40)ts_cor(ts_datos, lag.max=40)modelo1 <- stats::arima(ts_datos, order = c(0, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 12))Se calcula la autocorrelación parcial (PACF) para la serie temporal diferenciada (diff(ts_datos)) y la serie temporal original (ts_datos) hasta un retraso máximo de 40. La autocorrelación parcial es útil para identificar el orden de los términos autorregresivos (AR) en un modelo ARIMA.
summary(modelo1)##
## Call:
## stats::arima(x = ts_datos, order = c(0, 1, 1), seasonal = list(order = c(0,
## 0, 1), period = 12))
##
## Coefficients:
## ma1 sma1
## -0.4732 0.0610
## s.e. 0.1069 0.0927
##
## sigma^2 estimated as 0.1396: log likelihood = -45.33, aic = 96.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02046692 0.3718321 0.2981068 -7.86989 27.5897 0.9205199
## ACF1
## Training set 0.06275459modelo1##
## Call:
## stats::arima(x = ts_datos, order = c(0, 1, 1), seasonal = list(order = c(0,
## 0, 1), period = 12))
##
## Coefficients:
## ma1 sma1
## -0.4732 0.0610
## s.e. 0.1069 0.0927
##
## sigma^2 estimated as 0.1396: log likelihood = -45.33, aic = 96.66print(modelo1)##
## Call:
## stats::arima(x = ts_datos, order = c(0, 1, 1), seasonal = list(order = c(0,
## 0, 1), period = 12))
##
## Coefficients:
## ma1 sma1
## -0.4732 0.0610
## s.e. 0.1069 0.0927
##
## sigma^2 estimated as 0.1396: log likelihood = -45.33, aic = 96.66coeftest(modelo1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.473162 0.106936 -4.4247 9.657e-06 ***
## sma1 0.060991 0.092670 0.6581 0.5104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07confint(modelo1)## 2.5 % 97.5 %
## ma1 -0.6827521 -0.2635713
## sma1 -0.1206393 0.2426210modelo## Series: ts_datos
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.5263 1.2830
## s.e. 0.0868 0.0726
##
## sigma^2 = 0.1288: log likelihood = -40.55
## AIC=87.11 AICc=87.35 BIC=95.07Los intervalos de confianza proporcionan una medida de la incertidumbre asociada con las estimaciones de los coeficientes. Para el coeficiente MA(1), por ejemplo, el intervalo de confianza del 95% va desde -0.6827521 hasta -0.2635713. Esto significa que tenemos un 95% de confianza de que el verdadero valor del coeficiente MA(1) está dentro de este rango. Un intervalo de confianza que no incluye cero indica que el coeficiente es estadísticamente significativo. En este caso, el intervalo de confianza para MA(1) no incluye cero, lo que sugiere que es significativo en el modelo.
9.2. Prueba Ljung-Box (Evaluación del Modelo)
error<-residuals(modelo1)
Box.test(error,type="Ljung-Box")##
## Box-Ljung test
##
## data: error
## X-squared = 0.42543, df = 1, p-value = 0.5142plot(error)
El valor del estadístico X-squared es 0.42543 con 1 grado de libertad.
El valor p indica la probabilidad de observar la estadística de prueba bajo la suposición de que la hipótesis nula es verdadera. En este caso, el valor de p es 0.5142, lo que sugiere que no hay suficiente evidencia para rechazar la hipótesis nula de no autocorrelación. Por lo tanto, los residuos parecen ser ruido blanco o aleatorio.
9.3. Prónostico
pronostico2<-forecast::forecast(modelo1,h=5)
pronostico2 ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Oct 2023 2.298549 1.819741 2.777356 1.566276 3.030821
## Nov 2023 2.285608 1.744416 2.826800 1.457926 3.113290
## Dec 2023 2.297803 1.700710 2.894897 1.384627 3.210979
## Jan 2024 2.287701 1.639509 2.935893 1.296377 3.279025
## Feb 2024 2.299326 1.603779 2.994873 1.235579 3.363073plot(pronostico2)
La anterior tabla, nos arroja los límites inferiores y superiores del intervalo de confianza del 80%. Esto significa que hay un 80% de probabilidad de que el valor real caiga dentro de este rango. De igual manera, apreciamos los límites inferiores y superiores del intervalo de confianza del 95%. El intervalo de confianza del 95% es más amplio que el del 80%, ya que proporciona una cobertura más amplia para el valor real, con un 95% de probabilidad.
Para octubre de 2023, el pronóstico puntual es 2.298549, con un intervalo de confianza del 80% entre 1.819741 y 2.777356, y un intervalo del 95% entre 1.566276 y 3.030821. Se proporcionan resultados similares para los meses subsiguientes (noviembre de 2023 a febrero de 2024).
9.4 Modelo 2, segunda diferencia en la serie general.
ts_datos2 = diff(ts_datos,differences = 2)acf(ts_datos2, lag.max=40)
pacf(ts_datos2, lag.max=40)
acf(ts_datos, lag.max=40)
pacf(ts_datos, lag.max=40)
modelo2 <- stats::arima(ts_datos, order = c(2, 2, 1),
seasonal = list(order = c(1, 0, 1), period = 12))
modelo2##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.98summary(modelo2)##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.98
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.005736462 0.3715231 0.2922434 -7.9628 26.82681 0.9024142
## ACF1
## Training set -0.02741017modelo2##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.98print(modelo2)##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.98coeftest(modelo2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.360438 0.099012 -3.6403 0.0002723 ***
## ar2 -0.231431 0.096163 -2.4067 0.0160987 *
## ma1 -0.999999 0.034086 -29.3375 < 2.2e-16 ***
## sar1 0.310210 0.984412 0.3151 0.7526685
## sma1 -0.217390 0.999455 -0.2175 0.8278121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo2##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.98confint(modelo2)## 2.5 % 97.5 %
## ar1 -0.5544987 -0.16637712
## ar2 -0.4199066 -0.04295634
## ma1 -1.0668069 -0.93319206
## sar1 -1.6192011 2.23962205
## sma1 -2.1762853 1.74150560modelo2##
## Call:
## stats::arima(x = ts_datos, order = c(2, 2, 1), seasonal = list(order = c(1,
## 0, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## -0.3604 -0.2314 -1.0000 0.3102 -0.2174
## s.e. 0.0990 0.0962 0.0341 0.9844 0.9995
##
## sigma^2 estimated as 0.1407: log likelihood = -47.99, aic = 107.989.5. Prueba Ljung-Box (Evaluación del Modelo)
error<-residuals(modelo2)
Box.test(error,type="Ljung-Box")##
## Box-Ljung test
##
## data: error
## X-squared = 0.081164, df = 1, p-value = 0.7757plot(error)
9.6. Prónostico
pronostico<-forecast::forecast(modelo2,h=5)
pronostico ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Oct 2023 2.384434 1.901469 2.867400 1.645803 3.123066
## Nov 2023 2.376726 1.801030 2.952423 1.496274 3.257178
## Dec 2023 2.429677 1.795043 3.064311 1.459089 3.400265
## Jan 2024 2.421388 1.707279 3.135497 1.329252 3.513524
## Feb 2024 2.440954 1.658141 3.223766 1.243745 3.638162plot(pronostico)
9.7. Modelo 3, una diferencia en el componente estacional.
acf(diff(ts_datos), lag.max=40)
pacf(diff(ts_datos), lag.max=40)
acf(diff(ts_datos,12), lag.max=40)
pacf(diff(ts_datos,12), lag.max=40)
ts_cor(diff(ts_datos), lag.max=40)ts_cor(diff(ts_datos,12), lag.max=40)modelo3 <- stats::arima(ts_datos, order = c(2, 1, 1),
seasonal = list(order = c(1, 1, 1), period = 12))
modelo3##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.16summary(modelo3)##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.16
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.007545159 0.3280768 0.2436925 -7.747747 23.82477 0.7524946
## ACF1
## Training set -0.004050947modelo3##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.16print(modelo3)##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.16coeftest(modelo3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.085857 0.242298 0.3543 0.723081
## ar2 -0.088296 0.154969 -0.5698 0.568837
## ma1 -0.598833 0.231817 -2.5832 0.009789 **
## sar1 -0.031083 0.115998 -0.2680 0.788732
## sma1 -0.999999 0.167652 -5.9647 2.45e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo3##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.16confint(modelo3)## 2.5 % 97.5 %
## ar1 -0.3890388 0.5607524
## ar2 -0.3920288 0.2154373
## ma1 -1.0531868 -0.1444796
## sar1 -0.2584350 0.1962698
## sma1 -1.3285913 -0.6714075modelo3##
## Call:
## stats::arima(x = ts_datos, order = c(2, 1, 1), seasonal = list(order = c(1,
## 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sma1
## 0.0859 -0.0883 -0.5988 -0.0311 -1.0000
## s.e. 0.2423 0.1550 0.2318 0.1160 0.1677
##
## sigma^2 estimated as 0.1228: log likelihood = -47.58, aic = 107.169.8. Prueba Ljung-Box (Evaluación del Modelo)
error<-residuals(modelo3)
Box.test(error,type="Ljung-Box")##
## Box-Ljung test
##
## data: error
## X-squared = 0.0017728, df = 1, p-value = 0.9664plot(error)
9.9. Prónostico
pronostico3 <- forecast::forecast(modelo3,h=5)
pronostico3## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Oct 2023 2.112401 1.635783 2.589020 1.3834767 2.841326
## Nov 2023 1.913561 1.383433 2.443690 1.1027992 2.724324
## Dec 2023 1.870680 1.314264 2.427095 1.0197155 2.721644
## Jan 2024 1.885285 1.299680 2.470889 0.9896795 2.780890
## Feb 2024 1.781102 1.165144 2.397060 0.8390756 2.723128plot(pronostico3)
10. CRITERIOS DE INFORMACIÓN
10.1. BIC
Es un criterio de selección de modelo que penaliza la complejidad del modelo. Cuanto menor sea el valor de BIC, mejor será el modelo.
BIC(modelo1)## [1] 104.5887BIC(modelo2)## [1] 123.7909BIC(modelo3)## [1] 122.293310.2. MSE
print(class(pronostico3))## [1] "forecast"length(ts_datos)## [1] 105length(pronostico$mean)## [1] 5mse_modelo1 = mean((ts_datos[1:5] - pronostico$mean)^2) ; mse_modelo1## [1] 1.405426mse_modelo2 = mean((ts_datos[1:5] - pronostico2$mean)^2) ; mse_modelo2## [1] 1.148732mse_modelo3 = mean((ts_datos[1:5] - pronostico3$mean)^2) ; mse_modelo3## [1] 0.50862210.3. MAE
El MSE y el MAE son medidas adicionales de precisión y error, respectivamente. Al igual que con las métricas anteriores, valores más bajos indican un mejor rendimiento del modelo.
mae_modelo1 = mean(abs(ts_datos[1:5] - pronostico$mean)) ; mae_modelo1## [1] 1.176636mae_modelo2 = mean(abs(ts_datos[1:5] - pronostico2$mean)) ; mae_modelo2## [1] 1.059797mae_modelo3 = mean(abs(ts_datos[1:5] - pronostico3$mean)) ; mae_modelo3## [1] 0.678605810.4. MAPE
Es una medida de precisión que expresa el error como un porcentaje del valor real. Cuanto menor sea el MAPE, mejor será la precisión del modelo.
mape_modelo1 <- mean(abs((ts_datos[1:5] - pronostico$mean) /ts_datos[1:5])) * 100 ;
mape_modelo1## [1] 98.62016mape_modelo2 <- mean(abs((ts_datos[1:5] - pronostico2$mean) /ts_datos[1:5])) * 100 ; mape_modelo2## [1] 89.16419mape_modelo3 <- mean(abs((ts_datos[1:5] - pronostico3$mean) /ts_datos[1:5])) * 100 ; mape_modelo3## [1] 58.0544410.5. CRITERIO DE AIC
Al igual que el BIC, el AIC también penaliza la complejidad del modelo. Al igual que el BIC, cuanto menor sea el valor de AIC, mejor será el modelo.
modelos <- list(modelo1, modelo2, modelo3)
obtener_AIC <- function(modelos) {
return(AIC(modelos))
}
AICs <- sapply(modelos, obtener_AIC)
mejor_modelo <- modelos[[which.min(AICs)]]
print("Valores de AIC para los modelos:")## [1] "Valores de AIC para los modelos:"print(AICs)## [1] 96.65551 107.98250 107.16262print("El mejor modelo según el AIC:")## [1] "El mejor modelo según el AIC:"print(mejor_modelo)##
## Call:
## stats::arima(x = ts_datos, order = c(0, 1, 1), seasonal = list(order = c(0,
## 0, 1), period = 12))
##
## Coefficients:
## ma1 sma1
## -0.4732 0.0610
## s.e. 0.1069 0.0927
##
## sigma^2 estimated as 0.1396: log likelihood = -45.33, aic = 96.66| MODELO | BIC | AIC | RMSE | MAPE | MSE | MAE | | :——: | :——–: | :———:| :——-: |:———-:|:———-:|:———:| | 1 | 104.5887 | 96.65551 | 0.3718321 | 27.5897 | 1.405426 |1.176636 | | 2 | 123.7909 | 107.98250 | 0.3715231 | 26.826811 | 1.148732 |1.059797 | | 3 | 122.2933 | 107.16262 | 0.3280768 | 23.82477 | 0.508622 |0.6786058 |
Basado en estos criterios, el modelo 1 parece ser el más adecuado para el conjunto de datos de anomalías de temperatura en Colombia, ya que tiene los valores más bajos de BIC, AIC, RMSE, MAPE, MSE y MAE, lo que sugiere un mejor ajuste y mayor precisión en comparación con los otros modelos.
10.6. Justificación:
Modelo 1 (ARIMA(1,0,0) tiene el valor más bajo de AIC (96.65551), lo que lo hace el mejor modelo según este criterio. Modelo 3 tiene valores más bajos en RMSE, MAPE, MSE y MAE, lo que sugiere un mejor rendimiento en términos de ajuste y precisión. En el contexto de la serie de tiempo de anomalías de temperatura en Colombia, la elección del mejor modelo depende de qué métrica se prioriza. Si el objetivo es minimizar el AIC, el Modelo 1 es el adecuado. Si se busca un modelo con mejores métricas de error (precisión y ajuste), el Modelo 3 es el más adecuado.
11. RESUMEN:
El análisis y modelado de las anomalías de temperatura en Colombia se llevaron a cabo utilizando series de tiempo y métodos ARIMA. El objetivo principal fue entender y predecir estas anomalías, que son variaciones en las temperaturas que se desvían de los valores promedio históricos.
El modelo ARIMA(1,0,0) con media no nula fue seleccionado como el más adecuado, basado en el AIC más bajo, indicando un buen equilibrio entre ajuste y complejidad. Aunque otros modelos mostraron métricas de error menores, el AIC y la prueba de Ljung-Box apoyan la selección del Modelo 1 como el más parsimonioso y adecuado para pronosticar las anomalías de temperatura en Colombia. Los pronósticos son confiables y ajustados a la realidad dentro de los límites de incertidumbre proporcionados por los intervalos de confianza. Este modelo captura la dinámica de las anomalías de temperatura con una buena precisión y balance entre ajuste y parsimony.Del mismo modo, la prueba de Dickey-Fuller Aumentada confirmó que la serie diferenciada es estacionaria.Y la prueba de Ljung-Box aplicada a los residuos del modelo no mostró autocorrelaciones significativas, lo que valida la adecuación del modelo.
Los pronósticos generados para los próximos cinco meses muestran una continuidad en las anomalías de temperatura, con intervalos de confianza que reflejan la incertidumbre inherente a las predicciones. Lo que implica que las anomalías de temperatura en Colombia pueden tener múltiples implicaciones, desde impactos en la agricultura y la biodiversidad hasta efectos en la salud pública y la economía. Un modelo estadístico robusto como el ARIMA(1,0,0) proporciona una herramienta valiosa para anticipar estos cambios y preparar respuestas adecuadas. La capacidad de prever anomalías de temperatura permite a los tomadores de decisiones implementar medidas preventivas y adaptativas, mejorando la resiliencia del país ante los cambios climáticos.
Para cerrar, el modelo desarrollado no solo ayuda a entender las variaciones históricas de temperatura en Colombia, sino que también ofrece una base sólida para realizar pronósticos futuros. Estos pronósticos son cruciales para la planificación y gestión de recursos en sectores vulnerables al cambio climático, contribuyendo así a la mitigación de sus efectos adversos y a la adaptación estratégica en el país.