Modelo ARIMA
Julieth Cerón
2023-03-27
ARIMA es un modelo estadístico que utiliza la información pasada y actual para predecir futuros valores en una serie de tiempo.
El modelo ARIMA tiene tres parámetros principales: p, d y q.
p es el orden del componente autoregresivo, que utiliza valores pasados de la serie de tiempo para predecir valores futuros.
d es el orden de diferenciación, que se utiliza para hacer que la serie de tiempo sea estacionaria (es decir, con una media constante y una varianza constante) mediante la eliminación de tendencias y estacionalidad.
q es el orden del componente de media móvil, que utiliza errores pasados para predecir errores futuros.
Para implementar un modelo ARIMA en R con series de tiempo, podemos seguir los siguientes pasos:
ANÁLISIS EXLORATORIO DE LOS DATOS
Cargar datos
Recomendación: Para cargar los datos desde un archivo de excel, se recomienda eliminar la columna fecha del archivo de excel, sólo dejar las columnas con la información de las dos variables a analizar.
#PASO 1: Instalar y cargar la librería "readxl" en R, para cargar los datos:
#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, en este caso "readxl"y dar click en install.
#Después, se carga la librería con el siguiente comando:
library(readxl)
#PASO 2: Carga el archivo de Excel utilizando la función read_excel() de la librería "readxl".
#Debes especificar la ruta donde se encuentra el archivo, y el nombre de la hoja que contiene los datos.
#Se debe evitar colocar tildes en el nombre del archivo y procurar guardarlo en una ubicación que no tenga muchas subcarpetas de acceso.
#En este caso, mi base de datos se llamará "datos":
datos <- read_excel("D:/datarima.xlsx")2. Declarar la base de datos en serie de tiempo
#Utilizamos la función ts
#Debes especificar tres argumentos, separados por comas
#En este ejemplo sería:
# argumento 1: datos$CEM : Hace referencia a la base de datos cargada previamente y el nombre de la variable o columna de mi excel relacionada con despachos de cemento.
#frequency: la frecuencia de los valores de la serie de tiempo (por ejemplo, si los valores son mensuales, la frecuencia es 12, si son trimestrales, la frecuencias sería 4)
#start: el año y el mes o trimestre en que empieza la serie de tiempo o en este caso CEM o IPIR
cem<-ts(datos$CEM[1:96], frequency=12, start=c(2014,1))
ipir<-ts(datos$IPIR, 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
#Recordar que previo al cargue de librerias estas se deben instalar primero.
#Luego si cargarlas:
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 estacionalidad4.Graficar las series de tiempo y análisis exploratorio general de la serie
Las variables de análisis son despachos de cemento y el Índice de producción Industrial para Colombia. Son datos en frecuencia mensual, disponibles desde enero 2014 hasta diciembre 2022.
De acuerdo a la figura 1 y 2, se puede apreciar que tanto el sector de la construcción como la industria nacional en los dos último años se vieron fuertemente afectados por dos choques sin precedentes: la pandemia del Covid-19, especialmente durante el mes de abril 2020 y elparo nacional en mayo 2021 (aunque se evidencia un menor impacto de este suceso).En el año 2022 se evidencia una tendencia de recuperación que durante los últimos meses se ha desacelerado.
#Utilizamos la función autoplot para graficar
autoplot(cem,frequency=12,xlab="Años",ylab="No.Despachos",main="Figura 1.Despachos de cemento") ## Warning in ggplot2::geom_line(na.rm = TRUE, ...): Ignoring unknown parameters:
## `frequency`
autoplot(ipir,frequency=12,xlab="Años",ylab="IPIR",main="Figura 2. Indice de Produccion Industrial") ## Warning in ggplot2::geom_line(na.rm = TRUE, ...): Ignoring unknown parameters:
## `frequency`
EXTRACCIÓN DE SEÑALES
1.Graficar los datos, analizar patrones y observaciones atìpicas. Apoye su anàlisis con una descomposiciòn o extración de señales.
Se debe obtener Componente estacional, tendencia-ciclo, componente irregular y serie ajustada por estacionalidad.
#Utilizamos la función decompose (del paquete cargado previamente "STATS")
cem_decomp<- decompose(cem)
ipir_decomp <- decompose(ipir)
# Graficar los componentes
#Despachos de cemento
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(cem_decomp$x, main = "Despachos de cemento-Original", col = "black", ylab = "Serie de tiempo")
plot(cem_decomp$trend, main = "Tendencia", col = "blue", ylab = "Valores")
plot(cem_decomp$seasonal, main = "Estacionalidad", col = "red", ylab = "Valores")
plot(cem_decomp$random, main = "Irregularidad", col = "green", ylab = "Valores")
#IPIR
par(mfrow = c(2, 2)) #Se utiliza para dividir la ventana gráfica en una matriz de 2 filas y 2 columnas
plot(ipir_decomp$x, main = "IPIR-Original", col = "black", ylab = "Serie de tiempo")
plot(ipir_decomp$trend, main = "Tendencia", col = "blue", ylab = "Valores")
plot(ipir_decomp$seasonal, main = "Estacionalidad", col = "red", ylab = "Valores")
plot(ipir_decomp$random, main = "Irregularidad", col = "green", ylab = "Valores")
APLICACIÓN DEL MODELO ARIMA
PASO 1:IDENTIFICACIÓN
1. Validamos estacionariedad de las series
Un modelo ARIMA requiere que la serie sea estacionaria. Se dice que una serie es estacionaria cuando su media, varianza y autocovarianza son invariantes en el tiempo.Esta suposición tiene un sentido intuitivo: dado que ARIMA usa retardos previos de series para modelar su comportamiento.
RECORDAR: En el paso anterior se observó graficamente que la serie original del IPIR tiene una tendencia y ademas no tiene media ni varianza constante.Con lo cual pudieramos afirmar que visualmente la serie pareciera ser NO ESTACIONARIA. Para reafirmarlo, hacemos el Test de Dickey Fuller:. Este test se basa en una regresión lineal que incluye la propia serie de tiempo y sus rezagos.Las hipótesis respectivas son:
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.09. Como el p-valor > 0.05, no rechazamos H0. Podemos concluir que el IPIR es una variable NO Estacionaria.
#Cargar el paquete tseries
library(tseries)
# Realizar prueba de raíz unitaria
adf.test(ipir)##
## Augmented Dickey-Fuller Test
##
## data: ipir
## Dickey-Fuller = -3.1537, Lag order = 4, p-value = 0.09922
## alternative hypothesis: stationaryCuando la serie es NO ESTACIONARIA, debemos diferenciarla hasta que logre serlo. Pero debemos recordar que si la serie tiene un componente estacional importante podemos trabajar con la serie ajustada por estacionalidad. Por ello, a continuación haremos la comparación de la serie original y desestacionalizada del IPIR, para validar si la esatcionaldiad es unfactor importante en esta variable:
#Para hacer este proceso se usa la libreria seasonal instalada y cargada previamente.
#Colocamos el nombre a la nueva variable: ipir_SA (serie ajustada por estacionalidad)
ipir_SA <- seasadj(ipir_decomp)
# Graficar la serie de tiempo original y desestacionalizada en un mismo gráfico
plot(ipir, main = "Figura 3. IPIR- original y desestacionalizada")
lines(ipir_SA, col = "red")
legend("topleft", legend = c("Serie original", "Serie desestacionalizada"), col = c("black", "red"), lty = 1)
Como se observa en la Figura 3, hay un componente estacional relativamente en el IPIR. Por ello, se decide trabajar con la serie ajustada por estacionalidad.
Ahora, diferenciamos la serie ajustada por estacionalidad, para volverla estacionaria:
RECORDAR QUE: Las diferencias se pueden hacer en niveles o aplicando logaritmo (avaces esta última opción es recomendada para tener modelos ARIMA más robustos). Por eso si graficamente vemos que la serie no se vuelve estacionaria con 1 o 2 dferencias (en niveles), procedemos a realizar diferencias con logaritmo.
#Normalmente, lo ideal es diferenciar 1 o 2 veces, intentamos primero con 1 diferencia (en niveles) y graficamos
#Para observar el comportamiento, debería mostrar una media y varianza constante
#(no tendencia ni volatilidad)
ipirSA_d1= diff(ipir, differences = 1)
ipirSA_d2= diff(ipir, differences = 2)
plot(ipirSA_d1, main = "Figura 4. IPIR- Diferenciación en (niveles)")
lines(ipirSA_d2, col = "red")
legend("topleft", legend = c(" Una diferencia", "Dos diferencias"), col = c("black", "red"), lty = 1)
En la figura 4, se puede observar que al trazar la serie diferenciada (una y dos veces), vemos un patron oscilante alrededor de 0, sin una tendencia fuerte visible. Esto sugiere que la diferenciacion de los terminos del orden 1 es suficiente y debe incluirse en el modelo.
Hacemos el Test Dickey-Fuller para corroborar:
H0: No estacionaria H1: Estacionaria
CONCLUSIÓN P-value= 0.01< 0.05. Rechazamos H0, *el IPIR ahora es estacionario**.
# Realizar prueba de raíz unitaria
adf.test(ipirSA_d1)## Warning in adf.test(ipirSA_d1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ipirSA_d1
## Dickey-Fuller = -5.1015, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryadf.test(ipirSA_d2)## Warning in adf.test(ipirSA_d2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ipirSA_d2
## Dickey-Fuller = -6.8396, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary2. Determinamos p y q
A continuación, los picos en rezagos particulares de la serie diferenciada (IPIRSA_d1) pueden ayudar a informar la elección de p o q para nuestro modelo.Para ellos miraremos dos gráficos (ACF Y PACF):
La función de autocorrelación (ACF): mide la correlación entre una serie de tiempo y sus valores retrasados. Para determinar el valor de q en un modelo ARIMA, se puede mirar el gráfico ACF y buscar el primer retraso que tiene una correlación significativa (es decir en donde se sobrepasen las lineas punteadas. Este valor corresponderá a q (MA=1 O 2 o 3 o….q valor)
RecomendaciónPara iniciar a plantear los modelos, se puede partir de un modelo con AR(0) Ó MA(0).
En la Figura 5, se observa que el q podria tomar el valor de 1 o 2, ya que en estos primeros rezagos se evidencia una correlación significativa no decreciente o convergente a cero. Es decir que podria ser MA (1) o tal vez un MA (2).
Acf(ipirSA_d1, main='Figura 5.Función de autocorrelación (ACF) -IPIR diferenciado')
Ahora, para determinar el p, observaremos la función de autocorrelación parcial (PACF) que mide la correlación entre una serie de tiempo y sus valores retrasados, pero teniendo en cuenta la influencia de los valores intermedios.
En la Figura 6, se observa que el q podria tomar el valor de 0, O 1 O 2 ya que en estos primeros rezagos se sobrepasan las líneas punteadas pero no de form atan significativa. Es decir que podria ser AR (0) O AR (1) O AR (2).
Pacf(ipirSA_d1, main='Figura 6.Función de autocorrelación parcial (PACF) -IPIR 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:
Model1= (0,1,1) Model2= (0,1,2) Model3= (1,1,1) Model4= (1,1,2) Model5= (2,1,1) Model6= (2,1,2)
CONCLUSIÓN: El modelo óptimo seria el modelo 2, ya que de todos los modelos posibles es el que menor AIC y BIC representa.
NOTA:Para detectar de forma automática el mejor modelo, se puede usar la función auto.arima:
#Corremos la función auto.arima:
auto.arima(ipir_SA)## Series: ipir_SA
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.3263 -0.2215
## s.e. 0.0939 0.0905
##
## sigma^2 = 24.99: log likelihood = -323.16
## AIC=652.32 AICc=652.55 BIC=660.34PASO 3:VALIDACIÓN
Analizamos que los residuos sean Ruido Blanco (los residuos se distribuyen normalmente y no hay autocorrelación entre ellos).
Con la prueba de Ljung-Box, se evalúa si hay o no autocorrelación en los residuos:
Hipótesis H0: No hay autocorrelación de los residuos H1: Existe autocorrelación de los residuos
CONCLUSIÓN: Como el P-value (0.46) 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(ipir_SA,order = c(0,1,2))
Box.test(Model2$residuals, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: Model2$residuals
## X-squared = 19.961, df = 20, p-value = 0.4604Para validar si lor residuos se distribuyen de forma normal, se aplica el test de normalidad Shapiro-Wilk, que es una prueba estadística utilizada para evaluar si una muestra de datos sigue una distribución normal. Las hipótesis nula y alternativa del test de Shapiro-Wilk son las siguientes:
H0: Los residuos siguen una distribución normal. H1: Los residuos NO siguen una distribución normal.
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.
shapiro.test(Model2$residuals)##
## Shapiro-Wilk normality test
##
## data: Model2$residuals
## W = 0.85778, p-value = 8.973e-09PASO 4:PRONÓSTICO
#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
## Jan 2023 117.2965 107.4980 127.0949
## Feb 2023 118.1304 106.3156 129.9452
## Mar 2023 118.1304 105.5121 130.7488Pronostico6=forecast(Model2,level= c(95), h=6)
plot(Pronostico6)
Pronostico6## Point Forecast Lo 95 Hi 95
## Jan 2023 117.2965 107.4980 127.0949
## Feb 2023 118.1304 106.3156 129.9452
## Mar 2023 118.1304 105.5121 130.7488
## Apr 2023 118.1304 104.7567 131.5041
## May 2023 118.1304 104.0418 132.2191
## Jun 2023 118.1304 103.3614 132.8994MÉTRICAS DE EVALUACIÓN MAE Y RMSE
Investigar como calcular estas métricas en el pronóstico ARIMA (puede ser en excel, R, Eviews,etc)
Pista: Se debe pronosticar dentro de la muestra y no fuera de ella para poder validar, es decir pronosticar por ejemplo en este ejemplo 2019 ,2020 o 2021 o 2022 comparar valores reales con pronóstico.**
By Julieth Cerón
Predición Económica y Big Data 2023