Time series analysis
Santillán Martínex Axl
2023-04-04
Las librerías y paquetes necesarios para la ejecuación del análisis son las siguientes:
library(dygraphs)
library(readxl)
library(lmtest)## 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.numericlibrary(quantmod)## Loading required package: xts## Loading required package: TTR## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(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':
##
## sigmalibrary(rmgarch) ## Warning: package 'rmgarch' was built under R version 4.2.3##
## Attaching package: 'rmgarch'## The following objects are masked from 'package:xts':
##
## first, lastlibrary(forecast)## Warning: package 'forecast' was built under R version 4.2.3library(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Índice_de_residencia <- read_excel("Índice de residencia.xlsx")Los datos utilizados son el índice de inversión fija bruta en construcción residencial. El periodo de tiempo que abarca es de 2005 a 2020.
# Convertimos los datos a una serie de tiempo.
rsdn = Índice_de_residencia
rsdn_ts = ts(rsdn, start = c(2005,1), frequency = 12)Se espera que una serie de tiempo estcionaria se mantenga constante en el tiempo. Es decir, que no se obsrven patrones ni tendencias, por lo que las pruebas visuales son fundamentales para evaluar su comportamiento. En la Gráfica 1 no se aprecia un patrón claro, ni una tendencia evidente.Sin embargo, a partir de 2014 se obserba un crecimiento contínuo de la serie.
# Realizamos distintas gráficas para visualizar el comportamiento de los datos y poder comprobar si son estacionarios.
### Gráfica de serie de tiempo
rsdn_tsg = dygraph(rsdn_ts, main = "Índice de inversión en construcción residencial (2005 - 2020)", xlab = "Año", ylab = "Índice")%>% dyRangeSelector()
dyOptions(rsdn_tsg, colors = c("purple", "red"), fillGraph = TRUE, drawPoints = TRUE, pointSize = 2, strokeWidth = 2)En el histograma no se puede apreciar una distribución normal, además de que la media se localiza hacia la derecha. No se puede percibir un sesgo positivo ni negativo, aunque hay una caída brusca hacia la izquierda. En el centro hay una especie de planicie.
# Graficamos un histograma
hist(rsdn_ts, main = "Índice de inversión en construcción residencial (2005 - 2020)")
El diagrama de densidad nos indica que hay un cierto grado de heterogeneidad en los datos, puesto que hay dos prominencias. Debido a esto, la distribución no se aprecia como una normal.
# Hacemos un gráfico de densidad.
density = plot(density(rsdn_ts), main = "Índice de inversión en construcción residencial (2005 - 2020)")
Podemos observar que una gran candidad de datos se concentra cerca de la línea de normalidad, sin embargo los extremos no parecen apegarse a la línea. El gréfico qqplot tampoco denota normalidad a simple vista, lo que podría indicar que la serie no es estacionaria.
# Hacemos un gráfico qqplot
qqnorm(rsdn_ts)
qqline(rsdn_ts)
En las gráficas de puede apreciar claramente que la primera y segunda diferencia presentan menos tendencias y patrones que la serie original. Esto podría indicar que el modelo por diferencias podría ser más adecuado pues parece tener una media y varianza consantes.
# Ahora graficamos la serie de tiempo original y la primera y segunda diferencia.
par(mfrow=c(2,2))
plot(diff(rsdn_ts), type = "l", main = "Primera diferencia")
plot(diff(diff(rsdn_ts)), type = "l", main = "Segunda diferencia")
plot(rsdn_ts, type = "l", main = "Serie de tiempo original")
Determinar si la serie de tiempo requiere diferencias
La librería “nsdiffs” proporciona una función que utiliza pruebas de raíz unitaria estacionales para determinar la cantidad de diferencias estacionales necesarias para convertir una serie de tiempo “x” en estacionaria. Esta función también tiene en cuenta la posibilidad de requerir diferencias regulares con retraso uno. En resumen, nos ayuda a identificar el número de diferencias estacionales y regulares que debemos aplicar a una serie de tiempo para lograr la estacionariedad.
Por otro lado, la librería también ofrece la función “ndiffs”, la cual estima el número de diferencias regulares necesarias para convertir una serie de tiempo en estacionaria. Esta función se enfoca específicamente en las diferencias regulares y no considera las diferencias estacionales. Al proporcionar un análisis más enfocado en las diferencias regulares, nos permite determinar cuántas diferencias simples debemos aplicar a una serie de tiempo para lograr la estacionariedad.
En conjunto, estas dos librerías nos brindan herramientas útiles para identificar y calcular las diferencias necesarias en una serie de tiempo con el objetivo de lograr la estacionariedad. De acuerdo con la librería “ndiffs”, al utilizar la función correspondiente, se observa que solo se requiere una diferencia regular para lograr la estacionariedad en la serie de tiempo. Esto implica que aplicar una simple diferencia a los datos será suficiente para eliminar cualquier tendencia o patrón no estacionario y convertir la serie en estacionaria.
# Cálculo de las diferencias normales y estacionales
ndiffs(rsdn_ts)## [1] 1nsdiffs(rsdn_ts)## [1] 0# Análisis de las diferencias
rsdn.1d<-(diff(rsdn_ts));rsdn.1d## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2005 -1.87 1.08 2.26 -6.82 4.16 1.03 -2.26 4.70 -4.08
## 2006 -2.47 -0.67 8.99 -1.55 -2.56 0.09 8.08 1.11 4.79 -5.32
## 2007 -10.15 2.19 6.11 -3.94 -3.98 -0.78 9.42 1.21 4.19 1.89
## 2008 -16.76 4.52 -1.06 7.84 -4.37 -0.63 6.81 -9.66 -0.37 8.02
## 2009 -16.06 -3.56 5.29 -0.68 3.38 -1.94 11.74 -8.98 3.75 -3.50
## 2010 5.58 -13.76 0.43 8.48 -16.94 3.11 8.45 -11.51 6.10 1.51
## 2011 -0.44 -3.09 -4.09 10.21 -2.01 6.49 2.56 -6.29 -1.74 6.52
## 2012 -9.32 1.43 -4.01 2.44 1.69 -6.52 18.87 -13.43 1.39 0.32
## 2013 -0.74 10.95 -9.21 2.92 6.65 -13.35 3.66 -1.52 -0.41 3.02
## 2014 -11.53 15.22 -0.92 2.63 0.91 -5.44 -1.59 -2.77 0.30 10.61
## 2015 -8.76 5.93 -1.62 11.00 -6.16 -1.61 11.50 -3.90 6.83 -13.70
## 2016 -4.90 7.18 -4.41 6.95 1.93 -3.75 2.84 -7.11 4.28 -1.58
## 2017 -3.33 4.84 3.58 -8.99 10.54 1.47 -1.22 -2.41 6.16 -4.34
## 2018 -8.55 7.37 -9.82 4.88 5.89 -4.95 1.95 -15.48 4.19 4.03
## 2019 9.48 3.65 -12.88 -0.14 -0.84 0.68 -0.49 -2.80 -4.60 -0.76
## 2020 -1.08
## Nov Dec
## 2005 -1.07 6.54
## 2006 -5.20 6.75
## 2007 -4.07 10.09
## 2008 -12.40 5.34
## 2009 -4.53 11.06
## 2010 -4.52 4.93
## 2011 5.28 0.03
## 2012 6.27 -14.89
## 2013 -0.91 3.82
## 2014 -8.42 3.50
## 2015 1.34 3.57
## 2016 6.48 -9.06
## 2017 0.24 9.13
## 2018 -8.42 8.88
## 2019 6.00 4.55
## 2020# Gráfica de diferencia regular estacional
ggtsdisplay(rsdn.1d)
ACF y PACF con una diferencia
Después de aplicar la diferencia regular a la serie de tiempo, se examinó el correlograma ACF y se observó que el primer rezago presentaba una correlación significativamente mayor que los rezagos posteriores. Sin embargo, varios de los rezagos posteriores aún rebasaban las bandas de confianza, lo cual indica la presencia de autocorrelación en esos rezagos. Esto sugiere que aún puede haber información relevante en los rezagos más allá del primero.
Por otro lado, al analizar el PACF, se encontró que el primer rezago tenía una correlación considerablemente alta, mientras que solo uno de los rezagos posteriores rebasaba las bandas de confianza. Esto sugiere que una vez que se ha aplicado la diferencia regular, la mayoría de la autocorrelación se explica mediante el primer rezago y que la influencia de los rezagos posteriores es mínima.
En resumen, los resultados indican que la diferencia regular ha sido efectiva para eliminar la tendencia y convertir la serie en estacionaria.
Elección de modelo
En el análisis de la serie de tiempo, se utilizó una herramienta llamada auto.arima para elegir el modelo que mejor se ajusta a los datos. Esta herramienta utiliza un criterio llamado AIC para comparar diferentes modelos y seleccionar aquel que logra un buen equilibrio entre el ajuste a los datos y la simplicidad del modelo.
El AIC es una medida que tiene en cuenta tanto la capacidad del modelo para ajustarse a los datos como su complejidad. Al minimizar el AIC, se elige un modelo que logre un buen ajuste con la menor cantidad de parámetros posible. Esto implica que se prefieren modelos más simples en lugar de modelos complicados con muchos parámetros.
b_m = auto.arima(rsdn_ts)
summary(b_m)## Series: rsdn_ts
## ARIMA(0,1,1)(0,0,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.5838 0.2275
## s.e. 0.0686 0.0687
##
## sigma^2 = 30.75: log likelihood = -563.25
## AIC=1132.51 AICc=1132.64 BIC=1142.09
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1713791 5.498926 4.348359 -0.04414486 4.068242 0.6091855
## ACF1
## Training set 0.01223802Resultados del modelo
El modelo ARIMA(0,1,1)(0,0,1) no tiene componentes de estacionalidad. En este modelo, los términos (0,0,1) indican que no hay componentes estacionales presentes. El primer conjunto de números (0,1,1) indica que hay una diferencia regular de orden 1 (d=1) y un componente de promedio móvil de orden 1 (q=1).
Por lo tanto, este modelo ARIMA se utiliza para modelar datos que no tienen patrones estacionales y solo considera la tendencia y la dependencia serial en los datos a través de la diferencia regular y el componente de promedio móvil.
Análisis de residuos
#Pruebas a los Residuales del mejor modelo
RE<-residuals(b_m)
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 bp',xlab='Retardos')
pacf(RE2,main='PACF Residuales de bp',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.9916, p-value = 0.3749#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.02756, df = 1, p-value = 0.8681Box.test(RE2,1,type="Ljung")##
## Box-Ljung test
##
## data: RE2
## X-squared = 1.125, df = 1, p-value = 0.2889# Pruebas de heterocedasticidad condicional en el modelo
resultado_bp <- lmtest::bptest(RE2 ~ 1 + lag(RE2, k = 1) + lag(RE2, k = 2))
print(resultado_bp)##
## studentized Breusch-Pagan test
##
## data: RE2 ~ 1 + lag(RE2, k = 1) + lag(RE2, k = 2)
## BP = 0.40977, df = 1, p-value = 0.5221Resultados del análisis de los residuos
1. Prueba de normalidad de Shapiro-Wilk: Con un valor p de 0.3749, no hay evidencia suficiente para rechazar la hipótesis nula de que los residuos siguen una distribución normal.
2. Prueba de Box-Ljung (residuos): Con un valor p de 0.8681, no hay evidencia suficiente para rechazar la hipótesis nula de que no hay autocorrelación en los residuos.
3. Prueba de Box-Ljung (residuos al cuadrado): Con un valor p de 0.2889, no hay evidencia suficiente para rechazar la hipótesis nula de que no hay autocorrelación en los residuos al cuadrado.
4.Prueba de Breusch-Pagan: Con un valor p de 0.5221, no hay evidencia suficiente para rechazar la hipótesis nula de que no hay heterocedasticidad condicional en los residuos.
De acuerdo con la prueba BP (studentized Breusch-Pagan test), el p-valor obtenido es 0.5221, lo cual indica que no hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad condicional. Esto significa que no se encontró una relación significativa entre los residuos y sus rezagos al cuadrado, lo que sugiere que no es necesario aplicar un modelo adicional para modelar la volatilidad o varianza de los datos.
Modelo GARCH
dis_garch10<-ugarchspec(mean.model = list(armaOrder = c(1,1), include.mean=F),
variance.model = list(garchOrder = c(1,1)))
ugfit = ugarchfit (spec = dis_garch10, data = rsdn.1d);ugfit #Estimación del modelo##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.04344 0.139356 0.31172 0.75525
## ma1 -0.66471 0.114614 -5.79951 0.00000
## omega 1.68200 2.091819 0.80409 0.42135
## alpha1 0.00000 0.016855 0.00000 1.00000
## beta1 0.95004 0.100954 9.41061 0.00000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.04344 0.136173 0.31901 0.74972
## ma1 -0.66471 0.129461 -5.13440 0.00000
## omega 1.68200 3.062962 0.54914 0.58291
## alpha1 0.00000 0.008313 0.00000 1.00000
## beta1 0.95004 0.099016 9.59478 0.00000
##
## LogLikelihood : -568.0928
##
## Information Criteria
## ------------------------------------
##
## Akaike 6.3677
## Bayes 6.4564
## Shibata 6.3662
## Hannan-Quinn 6.4037
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.001901 0.9652
## Lag[2*(p+q)+(p+q)-1][5] 1.661477 0.9934
## Lag[4*(p+q)+(p+q)-1][9] 4.388490 0.5982
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 1.735 0.18775
## Lag[2*(p+q)+(p+q)-1][5] 5.036 0.15022
## Lag[4*(p+q)+(p+q)-1][9] 10.091 0.04796
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.4693 0.500 2.000 0.49330
## ARCH Lag[5] 2.4676 1.440 1.667 0.37686
## ARCH Lag[7] 6.8802 2.315 1.543 0.09215
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.9862
## Individual Statistics:
## ar1 0.09396
## ma1 0.19010
## omega 0.11645
## alpha1 0.12398
## beta1 0.11707
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.4656 0.6421
## Negative Sign Bias 0.3098 0.7571
## Positive Sign Bias 0.4299 0.6678
## Joint Effect 1.2838 0.7330
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 22.44 0.2627
## 2 30 34.67 0.2157
## 3 40 43.11 0.2998
## 4 50 53.33 0.3112
##
##
## Elapsed time : 0.1543021```