Para llevar acabo el proyecto es necesario llamar las siguientes
librerias:
library(dygraphs)
library(readxl)
library(tseries) #Modelos Arma y prueba Dickey-Fuller
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(lmtest) #Análsis de los coeficientes ARIMA
## Warning: package 'lmtest' was built under R version 4.2.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(quantmod)
## Loading required package: xts
## Loading required package: TTR
library(rugarch)
## Warning: package 'rugarch' was built under R version 4.2.3
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
library(rmgarch)
## Warning: package 'rmgarch' was built under R version 4.2.3
##
## Attaching package: 'rmgarch'
## The following objects are masked from 'package:xts':
##
## first, last
library(forecast)
## Warning: package 'forecast' was built under R version 4.2.3
library(MTS)#prueba multiplicadores de Lagrange... Elementos ARCH-GARCH
## Warning: package 'MTS' was built under R version 4.2.3
##
## Attaching package: 'MTS'
## The following object is masked from 'package:TTR':
##
## VMA
ventas_mayoreo <- read_excel("Venta al mayores de vehículos automóviles.xlsx")
Ventas al mayoreo de automóviles en México. Fuente: Banxico
Los datos utilizados fueron obtenidos de Banxico y consisten en una
serie de tiempo de la venta de automóviles al mayoreo del año 2000 al
2012.
# Convertimos los datos a una serie de tiempo.
v_ma = ventas_mayoreo
v_ma_ts = ts(v_ma, start = c(2000,1), frequency = 12)
Patrones y tendencias en la serie de tiempo
Una serie de tiempo estacionaria es una serie de tiempo en la que
las propiedades estadísticas básicas son constantes en el tiempo y no
dependen del tiempo.En la Gráfica 1 no se aprecia un patrón claro, ni
una tendencia evidente. Sin embargo, se aprecian picos constantes todos
los diciembres, lo que podría indicar que la venta de automóviles
depende de la temporada del año y hace complicado que la serie sea
estacionaria.
# Realizamos distintas gráficas para visualizar el comportamiento de los datos y poder comprobar si son estacionarios.
### Gráfica de serie de tiempo
v_ma_tsg = dygraph(v_ma_ts, main = "Ventas al mayoreo de vehículos en unidades (2000 - 2012)", xlab = "Año", ylab = "Autos")%>% dyRangeSelector()
dyOptions(v_ma_tsg, colors = c("red", "red"), fillGraph = TRUE, strokeWidth = 2)
Histograma
El histograma tiene algunas caractrísticas de una distribución
normal, pero está muy apuntalado por el centro. Un histograma apuntalado
sugiere que la mayoría de los datos tienen un valor cercano a la media,
lo que implica que la distribución de los datos es estrecha o
concentrada en torno a un valor central. Esto es ideal para identificar
datos estacionarios, puesto se espera que tengan media y varianza
constantes.
Gráfico de densidad
Este tipo de distribución sugiere una distribución leptocúrtica, que
se caracteriza por una mayor concentración de datos alrededor de la
media y una cola más corta. Por lo tanto, si la gráfica de densidad
muestra una distribución leptocúrtica muy apuntalada en el centro, esto
indica que los datos son bastante homogéneos y no hay muchos valores
atípicos.
Una gráfica de densidad apuntalada en el centro no necesariamente
implica que los datos sean estacionarios, aunque puede ser un indicador
de que los datos están más estables y hay menos variabilidad en la
distribución de los mismos.
# Hacemos un gráfico de densidad.
density = plot(density(v_ma_ts), main = "Ventas al mayoreo de vehículos en unidades (2000 - 2012)")
Q-Qplot
En el caso de que algunos datos se alejen de la línea de referencia
en la parte derecha del QQ-plot, esto indica que hay algunos valores
extremos o atípicos en la distribución de los datos. Es decir, estos
valores son significativamente diferentes a los valores que se esperan
en una distribución normal, lo que sugiere que hay una cierta
heterogeneidad o variabilidad en los datos.
En cuanto a la estacionariedad de los datos, un QQ-plot con algunos
valores extremos puede indicar la presencia de cambios en la
distribución de los datos a lo largo del tiempo, lo que podría afectar
la estacionariedad de la serie.
# Hacemos un gráfico qqplot
qqnorm(v_ma_ts)
qqline(v_ma_ts)
Cáluculo de diferencias y diferencias estacionales
Antes de realizar las puebas de estacionareidad, se calculan las
diferencias que requiere el modelo, urilizando la función “ndiffs” de la
librería forecast. La función “ndiffs” utiliza una prueba de raíz
unitaria para determinar el número de diferencias necesarias. La prueba
de raíz unitaria se utiliza para evaluar si una serie de tiempo tiene
una raíz unitaria, lo que indica que la serie es no estacionaria. La
hipótesis nula de la prueba es que la serie tiene una raíz estacionaria,
mientras que la hipótesis alternativa es que la serie tiene una raíz
unitaria.
La estacionalidad se refiere a los patrones repetitivos que ocurren
en una serie de tiempo a intervalos regulares, como las estaciones del
año o los días de la semana. En muchas series de tiempo, especialmente
aquellas que exhiben un comportamiento estacional, es necesario tomar en
cuenta las diferencias estacionales para lograr la estacionariedad. Para
lograr la estacionareidad, es necesario comprobar cuántas diferencias
estacionales necesitamos. Se utiliza la función “nsdiffs” para
determinar las diferencias necesarias.
La función “nsdiffs” utiliza pruebas de raíz unitaria estacional
para determinar el número de diferencias estacionales requeridas. Estas
pruebas se utilizan para evaluar si una serie de tiempo tiene una raíz
unitaria estacional, lo que indica que la serie es no estacionaria en
términos de su componente estacional. La hipótesis nula de las pruebas
es que la serie tiene una raíz estacionaria en términos de su componente
estacional, mientras que la hipótesis alternativa es que la serie tiene
una raíz unitaria estacional.
# Cálculo de las diferencias normales y estacionales
ndiffs(v_ma_ts)
## [1] 1
nsdiffs(v_ma_ts)
## [1] 1
# Análisis de las diferencias
v_ma.1d<-(diff(v_ma_ts));v_ma.1d
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2000 5037 6366 -8044 8821 -8 -11015 5429 194 17498
## 2001 -9733 -527 7145 -15050 1203 940 8550 -4116 -333 9483
## 2002 -21726 -2364 736 8373 -2220 305 -2515 -13361 4650 20780
## 2003 -24304 -2520 -1629 4358 -5577 151 6410 -3472 7095 5628
## 2004 -45459 8207 4546 -5142 2589 -2985 -7414 2783 18782 4774
## 2005 -41073 3644 4827 -5208 -8096 5645 -9297 7141 7212 4955
## 2006 -33042 -2024 7814 -13353 971 2347 -13272 14187 5868 18123
## 2007 -35562 -1025 3550 -10768 2047 643 -3665 5144 6844 6656
## 2008 -6646 -10777 -1830 -1824 3041 -3596 459 4503 -2370 3818
## 2009 -21417 -6092 1065 -9477 1797 2276 3296 3193 1269 9671
## 2010 -17383 -1404 4584 -5683 2316 4879 -2426 2491 6553 7664
## 2011 -15708 -3391 13439 -18789 5327 6263 -6639 12694 -749 6771
## Nov Dec
## 2000 16063 -15319
## 2001 15784 2866
## 2002 850 7359
## 2003 9505 23712
## 2004 5141 18408
## 2005 8690 15277
## 2006 5695 6765
## 2007 4853 -1503
## 2008 -2502 5906
## 2009 6096 4425
## 2010 1066 3741
## 2011 2730 6356
plot(v_ma.1d,type="l",main= "Primeras diferencias")
par(mfrow=c(1,2))
acf(v_ma.1d [-1], main = "ACF Venta al mayoreo",sub="2001 - 2012")
pacf(v_ma.1d[-1],main = "PACF Venta al mayoreo", sub="2001 - 2012")
#Diferencia estacional / Graficos
v_ma.1ds<-(diff(v_ma_ts,lag = 12,diferences=1));v_ma.1ds
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2001 15289 9725 10504 3498 -4120 -3172 16393 6848 6321 -1694
## 2002 4219 2382 -4027 19396 15973 15338 4273 -4972 11 11308
## 2003 -1711 -1867 -4232 -8247 -11604 -11758 -2833 7056 9501 -5651
## 2004 -1798 8929 15104 5604 13770 10634 -3190 3065 14752 13898
## 2005 8616 4053 4334 4268 -6417 2213 330 4688 -6882 -6701
## 2006 1748 -3920 -933 -9078 -11 -3309 -7284 -238 -1582 11586
## 2007 -2441 -1442 -5706 -3121 -2045 -3749 5858 -3185 -2209 -13676
## 2008 6130 -3622 -9002 -58 936 -3303 821 180 -9034 -11872
## 2009 -26589 -21904 -19009 -26662 -27906 -22034 -19197 -20507 -16868 -11015
## 2010 136 4824 8343 12137 12656 15259 9537 8835 14119 12112
## 2011 8073 6086 14941 1835 4846 6230 2017 12220 4918 4025
## Nov Dec
## 2001 -1973 16212
## 2002 -3626 867
## 2003 3004 19357
## 2004 9534 4230
## 2005 -3152 -6283
## 2006 8591 79
## 2007 -14518 -22786
## 2008 -19227 -11818
## 2009 -2417 -3898
## 2010 7082 6398
## 2011 5689 8304
plot(v_ma.1ds,type="l")
#Diferencia regular estacional / Graficos
v_ma.1drs<-(diff(v_ma.1d,lag = 12,diferences=1));v_ma.1drs
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2001 -5564 779 -7006 -7618 948 19565 -9545 -527 -8015
## 2002 -11993 -1837 -6409 23423 -3423 -635 -11065 -9245 4983 11297
## 2003 -2578 -156 -2365 -4015 -3357 -154 8925 9889 2445 -15152
## 2004 -21155 10727 6175 -9500 8166 -3136 -13824 6255 11687 -854
## 2005 4386 -4563 281 -66 -10685 8630 -1883 4358 -11570 181
## 2006 8031 -5668 2987 -8145 9067 -3298 -3975 7046 -1344 13168
## 2007 -2520 999 -4264 2585 1076 -1704 9607 -9043 976 -11467
## 2008 28916 -9752 -5380 8944 994 -4239 4124 -641 -9214 -2838
## 2009 -14771 4685 2895 -7653 -1244 5872 2837 -1310 3639 5853
## 2010 4034 4688 3519 3794 519 2603 -5722 -702 5284 -2007
## 2011 1675 -1987 8855 -13106 3011 1384 -4213 10203 -7302 -893
## Nov Dec
## 2001 -279 18185
## 2002 -14934 4493
## 2003 8655 16353
## 2004 -4364 -5304
## 2005 3549 -3131
## 2006 -2995 -8512
## 2007 -842 -8268
## 2008 -7355 7409
## 2009 8598 -1481
## 2010 -5030 -684
## 2011 1664 2615
plot(v_ma.1drs,type="l",main="v_ma_ts diferencia estacional y regular")
Tanto la función “ndiffs” como la función “nsdiffs” arrojaron un
valor de 1 en una serie de tiempo, esto indica que se requiere una
diferencia de orden 1 y una diferencia estacional de orden 1 para lograr
la estacionariedad en la serie.
La diferencia de orden 1 implica restar cada observación de la serie
de tiempo a su observación anterior, lo que ayuda a eliminar la
tendencia o el patrón de crecimiento lineal presente en la serie. Por
otro lado, la diferencia estacional de orden 1 implica restar cada
observación de la serie de tiempo a su observación correspondiente en el
mismo período de tiempo en la temporada anterior, lo que ayuda a
eliminar los efectos de las fluctuaciones estacionales.
Prueba Dickey Fuller y eliminación de tendencias
La diferencia regular con la diferencia estacional de la serie de
tiempo pasaron el test Dickey-Fuller con un valor de p de 0.01, y al
observar los gráficos de las diferencias se aprecia que se ha eliminado
el patrón de crecimiento lineal y las fluctuaciones estacionales,
podemos concluir lo siguiente:
1. Estacionariedad lograda: El hecho de que ambas diferencias hayan
pasado el test Dickey-Fuller con un valor de p de 0.01 indica que se ha
logrado la estacionariedad en la serie de tiempo. Esto significa que las
propiedades estadísticas de la serie se mantienen constantes a lo largo
del tiempo, lo cual es esencial para el análisis de series de
tiempo.
2. Eliminación del patrón de crecimiento lineal: Al observar los
gráficos de las diferencias y notar que se ha eliminado el patrón de
crecimiento lineal, significa que la serie original presentaba una
tendencia o un patrón de crecimiento a lo largo del tiempo. La
diferencia regular aplicada ha permitido eliminar esta tendencia, lo
cual es importante para modelar adecuadamente la serie y realizar
pronósticos más precisos.
3. Eliminación de las fluctuaciones estacionales: El hecho de que se
haya observado la eliminación de las fluctuaciones estacionales en los
gráficos de las diferencias implica que la serie original presentaba
patrones de variación recurrentes en determinadas estaciones o períodos
de tiempo. La diferencia estacional aplicada ha permitido eliminar estas
fluctuaciones estacionales, lo cual es crucial para eliminar los efectos
de las estacionalidades y modelar adecuadamente la serie.
# Dickey Fuller
adf.test(v_ma.1d) #Regular
## Warning in adf.test(v_ma.1d): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: v_ma.1d
## Dickey-Fuller = -6.2431, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
adf.test(v_ma.1ds)#Estacional
##
## Augmented Dickey-Fuller Test
##
## data: v_ma.1ds
## Dickey-Fuller = -2.1923, Lag order = 5, p-value = 0.4965
## alternative hypothesis: stationary
adf.test(v_ma.1drs)#Regular estacional
## Warning in adf.test(v_ma.1drs): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: v_ma.1drs
## Dickey-Fuller = -4.5543, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
ACF y PACF
En ambas gráficas, se observa que el primer rezago sea el más
correlacionado debido a la diferencia regular, y los rezagos posteriores
tienen autocorrelaciones bajas. Sin embargo, algunos rezagos rebasan las
bandas de confianza, lo que sugiere la presencia de autocorrelación
residual.
# Gráfica de diferencia regular estacional
ggtsdisplay(v_ma.1drs)
Elección de modelo
El criterio de información de Akaike (AIC) es una medida utilizada
para comparar diferentes modelos estadísticos y determinar cuál de ellos
proporciona el mejor equilibrio entre ajuste y complejidad. El AIC se
basa en la idea de que un buen modelo debe ajustarse bien a los datos,
pero al mismo tiempo, debe ser lo más simple posible. Un modelo que
minimice el AIC, se puede obtener con la función “auto,arima”, por lo
que se procede a elegir el modelo.
Análisis de los residuos
En primer lugar, la prueba de normalidad de Shapiro-Wilk indica que
los residuales no siguen una distribución normal, ya que el valor p
asociado es significativamente menor que un nivel de significancia
común. Esto sugiere que los residuales no son consistentes con una
distribución normal.
A continuación, se realizaron dos pruebas de Box-Ljung para evaluar
la presencia de autocorrelación en los residuales. La primera prueba
muestra un valor p alto, lo que indica que no hay suficiente evidencia
para rechazar la hipótesis nula de no autocorrelación. Sin embargo, en
la segunda prueba, aunque el valor p es mayor a un nivel de
significancia común, está cerca del límite establecido. Esto sugiere que
podría existir cierta autocorrelación en los residuales, aunque no es
concluyente.
En resumen, los resultados indican que los residuales del modelo
SARIMA presentan desviaciones de la normalidad y podrían tener cierta
autocorrelación.
#Pruebas a los Residuales del mejor modelo
RE<-residuals(mejor)
par(mfrow=c(1,2))
acf(RE,main='ACF Residuales de MA(6)',xlab='Retardos')
pacf(RE,main='PACF Residuales de MA(6)',xlab='Retardos')
#Pruebas a los Residuales al cuadrado del mejor modelo
RE2<-RE^2
#par(mfrow=c(1,2))
acf(RE2,main='ACF Residuales de MA (6) al cuadrado',xlab='Retardos')
pacf(RE2,main='PACF Residuales de MA (6) al cuadrado',xlab='Retardos')
# Realizamos la prueba de normalidad Shapiro Wilk test
#H0:La distribución es normal
#H1:La distribución no es normal
shapiro.test(RE)
##
## Shapiro-Wilk normality test
##
## data: RE
## W = 0.97035, p-value = 0.003239
#pvalue<=0.05 rechazo Ho
#pvalue> 0.05 acepto Ho
#Buscamos valores más grandes a 0.05 en el pvalue para
#Afirmar que los residuos provienen de una distribución normal
##supuestos del error: ho= ruido blanco; h1: los errores cuadrados estan relacioandos
Box.test(RE,1,type="Ljung") # VERIFICA QUE LOS RESIDUOS SE COMPORTAN COMO RUIDO BLANCO
##
## Box-Ljung test
##
## data: RE
## X-squared = 0.0010439, df = 1, p-value = 0.9742
Box.test(RE2,1,type="Ljung")
##
## Box-Ljung test
##
## data: RE2
## X-squared = 3.5123, df = 1, p-value = 0.06091
# Pruebas de heterocedasticidad condicional en el modelo
resultado_lm <- lmtest::bptest(RE2 ~ 1 + lag(RE2, k = 1) + lag(RE2, k = 2))
# Eliminar los valores NA
v_ma.1drss <- na.omit(v_ma.1drs)
# Eliminar los espacios en blanco
v_ma.1drsss <- gsub("\\s+", "", v_ma.1drss)
print(resultado_lm)
##
## studentized Breusch-Pagan test
##
## data: RE2 ~ 1 + lag(RE2, k = 1) + lag(RE2, k = 2)
## BP = 0.2784, df = 1, p-value = 0.5977
# Imprimir los resultados
Heterocedasticidad condicional
Se realizaron dos pruebas para evaluar la presencia de
heterocedasticidad en los residuos de un modelo. En primer lugar, se
llevó a cabo la prueba de Breusch-Pagan, también conocida como prueba
BP. Los resultados de la prueba BP indicaron un valor de estadístico BP
de 0.2784, con un grado de libertad de 1 y un valor p asociado de
0.5977. Esta prueba busca detectar la presencia de heterocedasticidad en
los residuos del modelo. En este caso, debido a que el valor p es mayor
que el nivel de significancia comúnmente utilizado (0.05), no se
encontró evidencia suficiente para rechazar la hipótesis nula de
homocedasticidad.
En segundo lugar, se realizó una prueba ARCH multivariante para
evaluar la presencia de autocorrelación en los residuos al cuadrado del
modelo. Los resultados de la prueba mostraron un estadístico Q(m) de
15.21097, con un valor p asociado de 0.9884915. Esta prueba busca
detectar la presencia de autocorrelación en los residuos al cuadrado, lo
cual indicaría la presencia de heterocedasticidad condicional en el
modelo. Sin embargo, en este caso, el valor p es mayor que el nivel de
significancia establecido, lo que sugiere que no hay suficiente
evidencia para rechazar la hipótesis nula de ausencia de autocorrelación
en los residuos al cuadrado.Por esta razón, se concluye que es posible
modelar la varianza.
