El 2022 fue un año muy difícil para las importaciones. No se trato de un fenómeno que afecto a unos países específicos, por el contrario, la afectación fue casi de índole global. Los factores que tuvieron incidencia fueron muy variados, inesperados y de consecuencias insuperables en el corto plazo. En primera instancia, apenas se empezaba la recuperación de todo el caos económico causado por las medidas tomadas para contrarrestar la pandemia del Covid-19, los aislamientos masivos de la población que no solo causo impacto en la demanda de muchos productos, sino también en la oferta, ya que las empresas vieran afectadas sus cadenas de producción debido a la ausencia de un gran porcentaje de la mano de obra calificada. Esto causo represamientos tanto en los envíos como en los pedidos internacionales. Cuando las economías empezaron a reactivarse en 2021, el gran acumulamiento de pedidos, colapso muchas de las cadenas de interdependencias complejas que hacen parte del comercio entre países, una de ellas, por ejemplo, la de los sistemas de transporte en contenedores, y fue así como por varios meses vimos en casi todos los mayores puertos del mundo, un gran represamiento de barcos cargueros de contenedores. Apenas se estaba superando esta crisis, en febrero del 2022 estallo la guerra Ucrania Rusia, afectándose la importantísima cadena de suministros de granos de provenientes de esos dos países, pero también las cadenas energéticas, que tienen gran incidencia en los precios de producción, transporte e importación de muchos productos. Igualmente, a nivel nacional, el paro y las protestas sociales del segundo semestre del 2021, tuvo repercusiones que salieron a la luz en el primer semestre del 2022, pues por varias semanas, el puerto de Buenaventura, por donde ingresan la mayoría de las importaciones del país fue bloqueado. Todos estos factores, hicieron que las importaciones durante el 2022, tuvieran muchos altibajos, exacerbados aún más por la gran devaluación del peso colombiano.

library(readxl) 
Bd <- read_excel("basesBIG/BD_Taller2_BC.xlsx")

2. Declarar la base de datos en serie de tiempo

M <-ts(Bd$MCOL[1:110], frequency=12, start=c(2014,1))
X <-ts(Bd$XCOL[1:110], frequency=12, start=c(2014,1))

3. Instalar y cargar librerias necesarias para el proceso

Forma 1: Escribir en el código install.packages(“nombre de la libreria) Forma 2: Seleccionar de la ventana lateral derecha la opción Packages/Install, colocar el nombre de la libreria y click en install

library(forecast)# Contiene el modelo ARIMA
library(tseries) #Para series de tiempo
library(TSA)     #Para series de tiempo
library(urca)    #Para hacer el Test de Raiz Unitaria (detectar hay o no estacionariedad)
library(ggplot2) #Para hacer gráficos
library(dplyr)   #Para la manipulación de datos (filtrar, seleccionar, agregar, transformar)
library(stats)   #Se usa para diversas pruebas estadísticas (medias,varianza, arima,etc)
library(seasonal)#Para calcular la serie ajustada de estacionalidad

4.Graficar las series de tiempo y análisis exploratorio general de la serie

Según se observa en las Figuras 1 y 2, el sector de importaciones y exportaciones experimentó un impacto significativo debido a un evento sin precedentes: el paro nacional de mayo de 2021. En 2022, se ha identificado una tendencia de recuperación, aunque en los últimos meses ha mostrado signos de desaceleración.

Un paro nacional puede afectar las importaciones de diversas maneras debido a la interconexión de factores económicos y sociales. En primer lugar, un paro puede dar lugar a bloqueos y protestas que afectan la logística y la cadena de suministro. Esto puede dificultar o retrasar la llegada de productos importados, lo que a su vez impacta la disponibilidad de bienes en el mercado local. Además, las tensiones sociales y políticas generadas por un paro pueden llevar a la incertidumbre económica, lo que disuade a los inversores y a los consumidores de gastar o invertir, lo que a su vez reduce la demanda de productos importados.

autoplot(M,frequency=12,xlab="Años",ylab="No.Importaciones",main="Figura 1.Importacines Totales ")
## Warning in ggplot2::geom_line(na.rm = TRUE, ...): Ignoring unknown parameters:
## `frequency`

autoplot(X,frequency=12,xlab="Años",ylab="No.Exportaciones",main="Figura 2. Exportaciones Totales") 
## Warning in ggplot2::geom_line(na.rm = TRUE, ...): Ignoring unknown parameters:
## `frequency`

EXTRACCIÓN DE SEÑALES

M_decomp<- decompose(M)
X_decomp <- decompose(X)

# Graficar los componentes

par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(M_decomp$x, main = "Importaciones", col = "black", ylab = "Serie de tiempo")
plot(M_decomp$trend, main = "Tendencia", col = "blue", ylab = "Valores")
plot(M_decomp$seasonal, main = "Estacionalidad", col = "red", ylab = "Valores")
plot(M_decomp$random, main = "Irregularidad", col = "green", ylab = "Valores")

APLICACIÓN DEL MODELO ARIMA

Vamos a llevar a cabo una prueba para determinar si nuestros datos son estacionarios. La estacionariedad de los datos es un aspecto fundamental al trabajar con modelos ARIMA (Autoregressive Integrated Moving Average), ya que influye en la capacidad de predecir de manera efectiva. En un modelo ARIMA, se asume que la serie temporal es estacionaria, lo que significa que sus propiedades estadísticas, como la media y la varianza, permanecen constantes a lo largo del tiempo. Si los datos no son estacionarios, es más complicado modelar sus patrones y tendencias, lo que puede resultar en predicciones menos precisas. Por lo tanto, determinar la estacionariedad es un paso crucial para asegurar que nuestro modelo ARIMA sea válido y útil para la predicción de series temporales.

Contraste de hipótesis:

H0: Serie No estacionaria: Hay raiz unitaria H1: Serie Estacionaria: No hay raiz unitaria

Tras realizar la prueba aumentada de Dickey-Fuller (ADF), obtenemos un p-valor = 0.1774. Como el p-valor > 0.05, no rechazamos H0. Podemos concluir que el Importaciones (M) es una variable NO Estacionaria.

library(tseries)
# Realizar prueba de raíz unitaria
adf.test(M)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  M
## Dickey-Fuller = -2.9623, Lag order = 4, p-value = 0.1774
## alternative hypothesis: stationary

Cuando la serie es NO ESTACIONARIA, Es esencial aplicar diferenciación a la serie hasta que alcance estacionariedad. Sin embargo, es fundamental recordar que cuando la serie contiene un componente estacional significativo, podemos optar por trabajar con la serie que ha sido ajustada por estacionalidad.

Por lo tanto, en el siguiente paso, compararemos la serie original de Importaciones (M) con la versión desestacionalizada para determinar si la estacionalidad es un factor relevante en esta variable. Esto es crucial para comprender si las variaciones en las importaciones están influenciadas por patrones estacionales y, en caso afirmativo, cómo afectan la interpretación de los datos y las predicciones.:

M_SA <- seasadj(M_decomp)
plot(M, main = "Figura 3. Importaciones - original y desestacionalizada")
lines(M_SA, col = "red")
legend("topright", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)

Observando la Figura 3, podemos notar la presencia de un patrón estacional en la serie de Importaciones. Esto significa que a lo largo del tiempo, hay fluctuaciones regulares que se repiten, como estacionalidades anuales o mensuales, que afectan los datos. Por lo tanto, se toma la decisión de trabajar con la serie de Importaciones que ha sido ajustada para eliminar estas variaciones estacionales, lo que nos permite centrarnos en las tendencias y patrones subyacentes más que en las fluctuaciones regulares.

Luego, para lograr que la serie ajustada por estacionalidad sea más adecuada para análisis posteriores, se procede a realizar una operación de diferenciación. Esta operación consiste en tomar la diferencia entre valores sucesivos de la serie en un intento de eliminar cualquier tendencia o patrón no estacionario que pueda quedar. En otras palabras, se busca convertir la serie en estacionaria, donde las observaciones no dependen del tiempo y se comportan de manera más uniforme. Esto facilita el análisis estadístico y la modelización de la serie de Importaciones.

M_d1= diff(M, differences = 1)
M_d2= diff(M, differences = 2)

plot(M_d1, main = "Figura 4. M - Diferenciación en (niveles)")
lines(M_d2, col = "red")

En la Figura 4, se muestra que al graficar la serie después de aplicar una diferenciación, ya sea una vez o dos veces, se observa un patrón oscilante que se mantiene alrededor del valor 0, sin que se perciba una tendencia marcada. Esto indica que la diferenciación de los términos de primer orden es suficiente y debe ser incluida en el modelo. En otras palabras, al aplicar una sola diferencia a la serie de datos, hemos eliminado la mayoría de las tendencias no estacionarias, lo que sugiere que no es necesario aplicar una diferencia adicional para obtener una serie estacionaria para el análisis y modelización.

Hacemos el Test Dickey-Fuller para corroborar:

H0: No estacionaria H1: Estacionaria

CONCLUSIÓN P-value= 0.01< 0.05. Rechazamos H0, *Las Importaciones (M) ahora es estacionario**.

adf.test(M_d1)
## Warning in adf.test(M_d1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  M_d1
## Dickey-Fuller = -5.2616, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
adf.test(M_d2)
## Warning in adf.test(M_d2): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  M_d2
## Dickey-Fuller = -7.8857, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

2. Determinamos p y q

Ahora, vamos a analizar algunos picos en la serie que hemos ajustado para entender mejor qué valores necesitamos para construir nuestro modelo. Para hacerlo, utilizaremos dos gráficos importantes: la Función de Autocorrelación (ACF) y la Función de Autocorrelación Parcial (PACF).

En la Figura 5, notamos que en los primeros puntos de los gráficos, vemos que hay una conexión significativa que no se desvanece rápidamente. Esto sugiere que podríamos necesitar un valor de “q” que sea igual a 1 o quizás 2 en este caso “q” pordria ser 1, 2, 3 o 4. Esto significa que podríamos tener que incorporar uno o dos términos de promedios móviles (MA) en nuestro modelo para que sea más preciso.

Acf(M_d1, main='Figura 5.Función de autocorrelación (ACF) -M diferenciado')

Ahora, para decidir cuántos términos “p” necesitamos en nuestro modelo, vamos a mirar la Función de Autocorrelación Parcial (PACF). Esta función nos ayuda a medir la relación entre una serie de tiempo y sus valores pasados, teniendo en cuenta la influencia de los valores intermedios.

En la Figura 6, observamos que el valor de “p” podría ser igual a 0, 1 o 2. Esto se debe a que en los primeros puntos del gráfico, la línea cruza las líneas punteadas, pero no de manera muy significativa. En otras palabras, podríamos necesitar un modelo autorregresivo (AR) con 0, 1 o 2 términos para que el modelo sea adecuado y se ajuste a los datos.

Pacf(M_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -M diferenciado')

PASO 2: ESTIMACIÓN DE LOS POSIBLES MODELOS

Recordar: Aquí la forma del modelo ARIMA es (p,d,q)

De acuerdo a lo observado en el paso 1, podemos probar los siguientes 6 modelos:

Model2= (2,1,1)

CONCLUSIÓN: El modelo óptimo seria el modelo 2, ya que de todos los modelos posibles es el que menor AIC y BIC representa.

auto.arima(M_SA)
## Series: M_SA 
## ARIMA(2,1,1) 
## 
## Coefficients:
##           ar1      ar2     ma1
##       -1.0375  -0.5383  0.3965
## s.e.   0.2028   0.1080  0.2377
## 
## sigma^2 = 8.773e+10:  log likelihood = -1526.74
## AIC=3061.48   AICc=3061.87   BIC=3072.25

PASO 3:VALIDACIÓN

Revisamos los residuos para asegurarnos de que se comporten como “Ruido Blanco”, que es un término utilizado para describir residuos que se comportan de manera aleatoria y sin patrones discernibles. En otras palabras, esperamos que los residuos se distribuyan normalmente y que no tengan una correlación significativa entre ellos. Si los residuos se parecen al ruido blanco, esto indica que el modelo se está ajustando bien a los datos y que no está dejando información importante sin explicar.

Para verificar si los residuos se asemejan al ruido blanco, utilizamos la prueba de Ljung-Box. Esta prueba examina si hay patrones de autocorrelación en los residuos, es decir, si los valores de los residuos en un momento dado están relacionados con los valores en momentos anteriores. La ausencia de autocorrelación significativa en los residuos es un indicativo de que se comportan como ruido blanco, lo cual es deseable en el análisis de series temporales. Esta evaluación es fundamental para garantizar la validez del modelo y para asegurarnos de que no se están dejando patrones no explicados en los datos.:

Hipótesis H0: No hay autocorrelación de los residuos H1: Existe autocorrelación de los residuos

CONCLUSIÓN: Como el P-value (0.98) es mayor a 0.05 no se rechaza H0. En ese caso si se cumple la condición de los residuos, son ruido blanco (no se correlacionan los errores).

Model2=Arima(M_SA,order = c(2,1,1))
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  Model2$residuals
## X-squared = 9.1708, df = 20, p-value = 0.9808

Utilizamos el test de normalidad Shapiro-Wilk para verificar si los residuos de nuestro análisis se asemejan a una distribución normal, es decir, si siguen una forma de campana típica. La prueba se basa en dos hipótesis: la hipótesis nula (H0), que asume que los datos siguen una distribución normal, y la hipótesis alternativa (H1), que sugiere que los datos no siguen una distribución normal y tienen un comportamiento diferente. En resumen, esta prueba nos permite evaluar si nuestros residuos se ajustan a una distribución normal, lo que es deseable en muchos análisis estadísticos, o si muestran un patrón de comportamiento diferente.

Si el valor p es menor que un nivel de significancia de 0.05, se puede rechazar la hipótesis nula y concluir que la muestra no sigue una distribución normal. Como el P-value (0.36) es mayor a 0.05 no se rechaza H0.

shapiro.test(Model2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  Model2$residuals
## W = 0.98701, p-value = 0.3692

PASO 4:PRONÓSTICO

La capacidad de anticipar y planificar el futuro es fundamental en una amplia variedad de contextos, desde el ámbito empresarial hasta el gobierno. Uno de los elementos clave para tomar decisiones informadas es la capacidad de prever eventos futuros, y una herramienta esencial para lograrlo es el análisis de series temporales. En este ensayo, exploraremos un ejemplo concreto de pronósticos de importaciones totales en Colombia y discutiremos por qué esta actividad es de importancia crítica.

La ciencia detrás de los pronósticos se basa en la idea de que los patrones observados en datos históricos pueden proporcionar información valiosa para predecir eventos futuros. Para lograr esto, se utilizan modelos de series temporales que examinan las relaciones entre las observaciones pasadas y utilizan esta información para proyectar el comportamiento futuro. En el ejemplo que estamos considerando, se ha aplicado un modelo de series temporales (posiblemente un modelo ARIMA) a los datos de importaciones totales en Colombia.

Una vez que se ha ajustado el modelo de series temporales, es posible utilizarlo para generar pronósticos. En este caso, se han realizado pronósticos para los próximos tres meses y también para los siguientes tres meses. Por ejemplo, se han proporcionado estimaciones de las importaciones totales esperadas para marzo, abril y mayo de 2023, así como para junio, julio y agosto del mismo año. Estos pronósticos se presentan en forma de punto de pronóstico, que representa la estimación más probable, junto con un intervalo de confianza del 95%, que indica el rango dentro del cual se espera que se encuentre el valor real.

#Se hace pronpostico a 3 y 06 meses. Se elije el parametro c(95) para que sea a un nivel de confianza del 95%

Pronostico3=forecast(Model2,level= c(95), h=3)
plot(Pronostico3)

Pronostico3
##          Point Forecast   Lo 95   Hi 95
## Mar 2023        3815892 3235375 4396410
## Apr 2023        3764518 3147722 4381314
## May 2023        3851138 3172948 4529328
Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)

Pronostico6
##          Point Forecast   Lo 95   Hi 95
## Mar 2023        3815892 3235375 4396410
## Apr 2023        3764518 3147722 4381314
## May 2023        3851138 3172948 4529328
## Jun 2023        3788927 2998518 4579336
## Jul 2023        3806841 2981456 4632226
## Aug 2023        3821744 2926958 4716531