Pronostico de serie de tiempo
2023-06-27
1 Justificación
Para el desarrollo de la asignatura “Análisis de series temporales”, se trabajarán con datos anuales relacionados con la venta de llantas en el parque automotor, especificamente de la empresa “Llanta Andina S.A”. El objetivo de identificar tendencias en las ventas, cambios en los patrones. Para esto se pronosticará la venta de llantas en Colombia, teniendo en cuenta las ventas historicas y posible relación con otras variables que influyen de forma directa e indirecta.
El pronostico de esta variable permitirá que las empresas“Llanta Andina S.A” puedan identificar oportunidades de mercado y tomar decisiones estratégicas basadas en datos, como el lanzamiento de nuevos promociones y servicios1.1 Modelo ARIMA
1.1.1 Análisis exploratorio
Primero, se llevará a cabo un análisis exploratorio de nuestra serie de interés “venta de llantas” para comprender mejor los datos y lograr identificar patrones y características importantes de esta variable. Así, se podran visualizar los datos y detectar patrones de tendencia, estacionalidad, ciclos, y ruido en la serie de tiempo. Esto es fundamental para elegir el modelo adecuado y para tomar decisiones informadas basadas en los datos de la serie.
# Instalar y cargar librerias necesarias para el proceso
library(fpp2)
library(readxl)
library(forecast)# Contiene el modelo ARIMA y Holt-Winters
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
library(zoo) #Para calcular la serie ajustada de estacionalidad
library(prophet) #Para es una biblioteca de Facebook que implementa el algoritmo Prophet para el pronóstico de series de tiempo.
library(tidyverse)
library(lubridate)
library(tidyverse)
library(quantmod)1.1.1.1 Serie original
Las variables de análisis es la venta de llantas nuevos en “Llanta Andina S.A”. Son datos en frecuencia anual, disponibles desde enero 2016 hasta febrero 2020.
De acuerdo a la figura 1, se puede apreciar que el comercio de los dos último años se vió fuertemente afectado porla pandemia del Covid-19 y las restricciones de aislamiento y movilidad para contener el avance de la misma, especialmente en el año 2020.
base <- read_excel("C:/Users/Natalia/Desktop/Venta_Llantas.xlsx")
vel<-ts(base$Ventas[1:100], frequency=12, start=c(2016,1))
autoplot(vel,frequency=12,xlab="Años",ylab="No. de llantas",main="Figura 1. Venta de Llantas") ## Warning in ggplot2::geom_line(na.rm = TRUE, ...): Ignoring unknown parameters:
## `frequency`
1.1.1.2 Promedio Móvil
El cálculo del promedio móvil es una técnica común utilizada en el análisis de series de tiempo para suavizar los datos y reducir el ruido. El objetivo es reducir la variabilidad en los datos, lo que puede hacer que las tendencias subyacentes sean más visibles. En esta caso se usará un promedio móvil de orden 3 para no perder tanta información relevante, sobretdo en los últimos dos años de análisis (2016-2020). En la Figura 2. Se observa que el MA(3) suaviza la serie temporal original y elimina la mayoría de las fluctuaciones de corto plazo.
#Calcular promedio llantas vendidas
promovil<- ma(vel, order = 3)
# Graficar serie original y promedio
ggplot() +
geom_line(aes(x = index(vel), y = vel, color = "Serie original")) +
geom_line(aes(x = index(promovil), y = promovil, color = "Promedio ")) +
labs(title = "Figura 2. Venta de vehículos con promedio ",
x = "Mes-Año",
y = "Número de vehículos",
color = "") +
theme_minimal()+
scale_color_manual(values = c("Serie original" = "black", "Promedio móvil MA(3)" = "purple"))## Warning: Removed 2 rows containing missing values (`geom_line()`).
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
1.1.1.3 Análisis de rezagos
La PACF es una medida de la correlación entre una observación y una observación retrasada, controlando el efecto de las observaciones intermedias. Un retraso significativo en la PACF puede indicar que ese número de retrasos es importante para explicar la serie de tiempo.
Hay que analizar la naturaleza de los datos para saber cuántas veces debes rezagar una serie de tiempo,con el objetivo del análisis que se está llevando a cabo. Una forma de determinar la cantidad adecuada de retrasos es mediante la prueba de autocorrelación parcial (PACF), que permite identificar los retardos significativos en una serie de tiempo.
En la gráfica de la PAFC anterior, se observa que un rezago es importante para explicar la serie. Por ende, a continuación aplicamos 1 rezago a la serie original de la venta de llantas:
library(stats)
rez_vel <- stats::lag(vel, k = 1)
plot(vel, main = "Serie original y rezagada")
lines(rez_vel, col = "red")
legend("topleft", legend = c("Serie original", "Serie rezagada"), col = c("black", "red"), lty = c(1, 1))
#### Estacionalidad
La función acf que devuelve un gráfico que muestra los coeficientes de correlación para cada rezago. Si la serie de tiempo tiene estacionalidad, se esperaría ver picos en los coeficientes de correlación en los múltiplos de la frecuencia de la serie (por ejemplo, si la frecuencia es mensual, se esperaría ver picos en los coeficientes de correlación para los rezagos 12, 24, 36, etc.). Estos picos indicarían la presencia de patrones de repetición en la serie, lo que sugiere la presencia de estacionalidad.
### Extracción de señales
La extracción de señales en series de tiempo es el proceso de identificar patrones, tendencias y características importantes en los datos de la serie temporal. Es una técnica común utilizada en análisis de series de tiempo para modelar el comportamiento de la serie, predecir valores futuros y entender las relaciones entre las variables.
La extracción de señales incluye el análisis de tendencias, el análisis de ciclos, el análisis de estacionalidad, la identificación de puntos atípicos y la descomposición de series de tiempo.
El objetivo de la extracción de señales es resumir la información en la serie de tiempo de una manera significativa y comprensible, lo que permite a los analistas y tomadores de decisiones identificar patrones y tendencias a largo plazo, así como patrones a corto plazo en los datos. La extracción de señales también puede ser útil para identificar la relación entre las variables de una serie de tiempo y cómo están cambiando a lo largo del tiempo.
En resumen, la extracción de señales en series de tiempo es un proceso crítico para comprender el comportamiento de los datos a lo largo del tiempo y proporciona información valiosa para la toma de decisiones y la predicción de eventos futuros.
1.1.1.4 Descomposición de las series
Es importante descomponer las series de tiempo porque permite identificar los diferentes componentes que la conforman, es decir, la tendencia, la estacionalidad y la variabilidad aleatoria. Cada uno de estos componentes puede proporcionar información valiosa sobre el comportamiento de la serie a lo largo del tiempo y su relación con otros factores.
La tendencia indica que la venta de llantas refleja un punto de quiebre en marzo 2020 generado por el impacto de la pandemia, luego se puede observar un ligero cambio de tendencia hacia la desaceleración.
La estacionalidad, por otro lado, refleja patrones repetitivos en la serie a lo largo del tiempo.
Por último, el componente irregular, también conocido como ruido, es la parte de la serie que no se puede explicar por la tendencia y la estacionalidad, y puede ser causada por factores impredecibles y/o eventos aleatorios. En este caso,el covid19, un evento sin precedentes e inesperado.
#Utilizamos la función decompose (del paquete cargado previamente "STATS")
library(stats)
vel_decomp<- decompose(vel)
# 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(vel_decomp$x, main = "Venta de llantas-Original", col = "black", ylab = "Serie de tiempo")
plot(vel_decomp$trend, main = "Tendencia", col = "blue", ylab = "Valores")
plot(vel_decomp$seasonal, main = "Estacionalidad", col = "red", ylab = "Valores")
plot(vel_decomp$random, main = "Irregularidad", col = "green", ylab = "Valores")
### Aplicación del Modelo Arima
1.1.1.5 Validación de Estacionariedad
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 de la venta de llantas 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 validarlo, 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.01. Como el p-valor < 0.05, se rechaza H0. En conclusión, la venta de llantas es una variable Estacionaria.
##
## Augmented Dickey-Fuller Test
##
## data: vel
## Dickey-Fuller = -3.2104, Lag order = 4, p-value = 0.09004
## alternative hypothesis: stationary1.1.1.6 Diferenciación
Dado que la venta de llantas es una variable estacionaria, no se requiere diferenciar la serie.
1.1.1.7 Estimación del modelo Auto.arima
La función auto.arima() es una herramienta muy útil para ajustar modelos ARIMA automáticamente. La idea detrás de esta función es seleccionar automáticamente el mejor modelo ARIMA para una serie de tiempo dada, basándose en criterios estadísticos como el criterio de información de Akaike (AIC) o el criterio de información bayesiano (BIC).
La notación ARIMA(p,d,q)(P,D,Q)[m] se refiere a los parámetros del modelo ARIMA, donde:
p: orden de la parte autoregresiva (AR) d: orden de diferenciación (I) q: orden de la parte de media móvil (MA) P: orden de la parte estacional autoregresiva (SAR) D: orden de la diferenciación estacional (SI) Q: orden de la parte estacional de media móvil (SMA) m: número de períodos en una temporada
En el modelo ARIMA(1,0,0)(0,1,1)[12], el orden AR es 1, el orden MA es 0, el orden de diferenciación es 0, el orden de SAR es 0, el orden de diferenciación estacional es 1, el orden SMA es 1, y el número de períodos en una temporada es 12.
La interpretación de cada parámetro es la siguiente:
-El parámetro AR(1) indica que se está utilizando la observación más reciente y la observación anterior para predecir la siguiente observación en la serie.
-El parámetro MA(0) indica que no se está utilizando ningún término de media móvil para hacer la predicción.
El parámetro de diferenciación d=0 indica que no se aplicó ninguna diferenciación a la serie.
El parámetro de diferenciación estacional D=1 indica que se aplicó una diferenciación estacional de primer orden para corregir la estacionalidad en la serie.
El parámetro SMA(1) indica que se está utilizando la observación de hace 12 períodos y la observación de hace 1 período para predecir la siguiente observación en la serie.
En resumen, el modelo ARIMA(1,0,0)(0,1,1)[12] es un modelo que utiliza la observación más reciente y la observación anterior para predecir la siguiente observación en la serie, y también tiene en cuenta la estacionalidad con una diferencia estacional de primer orden y una media móvil estacional de orden 1.
## Series: vel
## ARIMA(0,1,0)
##
## sigma^2 = 1.941e+12: log likelihood = -1541.03
## AIC=3084.05 AICc=3084.09 BIC=3086.651.1.1.8 Validación de supuestos
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
##
## Box-Ljung test
##
## data: autoarimaveh$residuals
## X-squared = 14.937, df = 20, p-value = 0.781.1.1.9 Predicción corto plazo
# Generar pronósticos futuros para 4 años (48 meses)
pronostico <- forecast(autoarimaveh, h = 48)
# Graficar la serie original y el pronóstico
plot(vel, main = "Serie de datos original y pronóstico", xlab = "Fecha", ylab = "Número de vehículos", xlim=c(2016,2020), ylim = c(0, max(vel)))
lines(pronostico$mean, col = "red")
legend("bottomleft", legend = c("Serie original", "Pronóstico"), col = c("black", "red"), lty = 1)
#Calculo del AIC Y BIC
# Calcular el AIC y el BIC
aic <- AIC(autoarimaveh)
bic <- BIC(autoarimaveh)
# Mostrar los resultados
cat("AIC:", aic, "\n")## AIC: 3084.053## BIC: 3086.6491.2 Evaluación del Modelo - MAE y RMSE
# Calcular las medidas de precisión
accuracy_arima <- accuracy(autoarimaveh)
# Obtener el MSE y el MAE
mae_arima <- accuracy_arima[1]
rmse_arima <- accuracy_arima[2]
# Mostrar los resultados
cat("MAE:", mae_arima, "\n")## MAE: -968.5269## RMSE: 13860421.3 Modelo Holt-Winters
Aquí podemos realizar dos ajustes y graficarlos en comparación con los datos originales para ver la calidad de los ajustes variando los paremetros alpha, beta y gamma.
HW1 <- HoltWinters(vel)
# Custom HoltWinters fitting
HW2 <- HoltWinters(vel, alpha=0.2, beta=0.1, gamma=0.1)
#Visually evaluate the fits
plot(vel, ylab="Ventas de Vehiculos", xlim=c(2014,2023))
lines(HW1$fitted[,1], lty=2, col="blue")
lines(HW2$fitted[,1], lty=2, col="red")
legend("bottomleft", legend = c("Serie original", "HW1","HW2"), col = c("black", "blue", "red"), lty = 1)
# Predicciones
Ambos ajustes parecen seguir bastante bien nuestros datos, así que ahora es el momento de ver cómo se desempeñan al predecir las ventas futura Usando la función “predict”, necesitaremos especificar cuántos puntos de datos queremos predecir en el futuro. Aquí, usaremos un valor de 24 para proyectar 2 años hacia el futuro (recuerdemos, esta es una serie temporal mensual). También nos gustaría tener una noción de los “intervalos de error” asociados con la predicción para tener una idea de nuestra confianza en la predicción. Para hacer esto, establecemos “prediction.interval=TRUE” y el nivel del intervalo de confianza (seleccionado aquí como 0.95). Una vez más, mostraremos un gráfico que incluya la cola de nuestros datos existentes y las nuevas predicciones.
HW1.pred <- predict(HW1, 24, prediction.interval = TRUE, level=0.95)
#gráfica
plot(vel, ylab="Ventas de llantas", xlim=c(2016,2020))
lines(HW1$fitted[,1], lty=2, col="blue")
lines(HW1.pred[,1], col="red")
lines(HW1.pred[,2], lty=2, col="orange")
lines(HW1.pred[,3], lty=2, col="orange")
legend("bottomleft", legend = c("Serie original", "HW1","pronostico","límites"), col = c("black", "blue", "red","orange"), lty = 1)
## Ajuste de la estacionalidad
Cuando realizamos ajustes, también tenemos la opción de ajustar el comportamiento de la componente de estacionalidad. Los ajustes estándar de Holt-Winters utilizan una estacionalidad aditiva, lo que significa que asumen que la amplitud de cualquier componente de estacionalidad es relativamente constante a lo largo de la serie. Sin embargo, si utilizamos una estacionalidad multiplicativa, permitimos que las variaciones estacionales (amplitud) crezcan con el nivel general de los datos. Para ver cómo funciona esto en nuestro caso de ventas de llantas, crearemos un nuevo ajuste, realizaremos predicciones en el futuro y las compararemos con nuestro ajuste aditivo de HW1.
## Warning in HoltWinters(vel, seasonal = "multiplicative"): optimization
## difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCHHW3.pred <- predict(HW3, 24, prediction.interval = TRUE, level=0.95)
plot(vel, ylab="Ventas de llantas", xlim=c(2016,2020))
lines(HW3$fitted[,1], lty=2, col="blue")
lines(HW3.pred[,1], col="red")
lines(HW3.pred[,2], lty=2, col="orange")
lines(HW3.pred[,3], lty=2, col="orange")
legend("bottomleft", legend = c("Serie original", "HW-mult","pronostico","límites"), col = c("black", "blue", "red","orange"), lty = 1)
Como podemos ver, la predicción se parece bastante a nuestro resultado de HW1, los intervalos de confianza muestran una tendencia más conservadora. Para este conjunto de datos, parece que el ajuste multiplicativo podría ser una mejor opción.
1.4 Uso de Forecast:
Utilizando nuestro ajuste de Holt-Winters HW1 anterior, podemos utilizar “forecast” para hacer nuevas predicciones e incluir intervalos de confianza del 80% y 95%.
library(forecast)
HW1_for <- forecast(HW1, h=24, level=c(80,95))
#grafica
plot(HW1_for, xlim=c(2016, 2024))
lines(HW1_for$fitted, lty=2, col="purple")
# Evaluación del modelo
Ahora calculamos la calidad de nuestras predicciones al compilar los valores observados menos los valores predichos para cada punto de datos. Estos se agregan a nuestro modelo de predicción como “$residuals”. Para evaluar mejor las funciones de suavizado que utilizamos en nuestro modelo, queremos verificar que no haya correlaciones entre los errores de predicción. En pocas palabras, si los puntos vecinos en nuestros ajustes constantemente se desvían de los valores observados de manera similar, nuestra línea de ajuste principal no es lo suficientemente reactiva a los cambios en los datos. Para capturar esto, utilizamos la función “acf” para evaluar la correlación de los residuos del ajuste entre puntos con diferentes separaciones temporales en la serie de tiempo (lag). Idealmente, para un lag no nulo, las barras de ACF deben estar dentro del rango de barras azules que se muestran a continuación. Es importante utilizar “na.action=na.pass” porque el último valor de “’$residuals” siempre es NA y, de lo contrario, la función mostrará un error.
Una prueba de Ljung-Box también puede indicar la presencia de estas correlaciones. Si obtenemos un valor de p > 0.05, hay un 95% de probabilidad de que los residuos sean independientes. Por último, es útil verificar el histograma de los residuos para asegurarnos de que sigan una distribución normal. Si los residuos están muy sesgados, entonces nuestro modelo puede estar constantemente sobrestimando en una dirección.
##
## Box-Ljung test
##
## data: HW1_for$residuals
## X-squared = 12.759, df = 20, p-value = 0.8875
## MAE y RMSE
# Calcular las medidas de precisión
accuracy_hw <- accuracy(HW1_for)
# Obtener el MSE y el MAE
mae_hw <- accuracy_hw[1]
rmse_hw <- accuracy_hw[2]
# Mostrar los resultados
cat("MAE:", mae_hw, "\n")## MAE: -121616.3## RMSE: 16434311.5 Comparación de los modelos ARIMA y Holt-Winters:
Comparando los resultados de los modelos ARIMA y Holt-Winters en términos de las métricas de precisión, podemos hacer las siguientes conclusiones:
MAE (Mean Absolute Error):
El modelo ARIMA tiene un MAE de -121.616.3 , mientras que el modelo Holt-Winters tiene un MAE de 164.343.1 En términos del MAE, el modelo ARIMA muestra un menor error absoluto promedio en comparación con el modelo Holt-Winters. Un valor de MAE más bajo indica que el modelo ARIMA tiene una mejor precisión en los pronósticos en comparación con el modelo Holt-Winters.