Examen Practico Medio Termino
Daniel Sandoval, A01733438, Karen Silva A01023307
06/05/2022
Introducción
El propósito del proyecto es formular e identificar modelos capaces de describir series de acuerdo al procedimiento de Box Jenkins. Se identificarán características de las serie tales como estacionariedad, estacionalidad, tendencias, etc. El modelo formulado debe cumplir con las condiciones de un proceso estacionario y podemos hacer estimaciones y sugerencias futuras sobre la información proporcionada.
El proyecto será dividido en 4 bloques: identificación, estimación, verificación y justificación. Se utilizó una base de datos con 700 observaciones bajo el nombre “SdT_MT_2”.
Identificación
Para el bloque de identificación es importante tener una idea aproximada sobre la serie y cuáles serán los pasos adecuados para su estimación. Por ende, iniciaremos el análisis con el método de análisis gráfico, estadística descriptiva y pruebas de raíz unitaria.
1.1 Análisis Gráfico
ts.plot(examen)
Podemos concluir de forma visual que la gráfica obtenida por el comando ts.plot(examen) no tiene patrones de tendencias, estacionalidad o ciclos, por lo que no tiene similitud con la caminata aleatoria. A partir de la metodología usada en este proyecto que supone la regresión de la media, a partir del teorema límite central se descarta que existan fenómenos de tendencia que distorsionen el comportamiento de la serie hacia su media.
boxplot(examen)
Por parte de las gráfica obtenida por el boxplot(examen), observamos que tenemos algunos datos fuera de lo normal en la serie pero en todo caso nuestra media se sigue centrando en el cero.
.
hist(examen)
Continuando con el análisis gráfico, en el histograma anterior se sugiere que la mayoría de los valores se centran en la parte media de la gráfica en media = 0, tenemos índices de normalidad, y esto reafirma el supuesto de estacionariedad, así como valores mínimos en su sesgo.
.
acf.examen<-acf(examen, main="ACF de la Serie")
pacf.examen<-pacf(examen, main="PACF de la Serie generada")
La función de autocorrelación parcial sugiere que el componente autorregresivo se concentra principalmente hasta p=3. Es importante mencionar que uno de los supuestos de Medias Móviles es que el componente es negativo al corregir el componente determinista del modelo. Por lo tanto, la gráfica sugiere que tenemos un q=4, o sea a el número de retardos de la serie.
1.2 Estadística Descriptiva
Continuando con la metodología, el siguiente bloque se conforma con la evaluación principalmente de valores como el sesgo y la curtosis con una demostración matemática, a pesar de ya contar con un análisis gráfico.
En línea con lo anterior, utilizaremos la siguiente tabla:
stat.desc(examen,norm = T)## x
## nbr.val 700.00000000
## nbr.null 0.00000000
## nbr.na 0.00000000
## min -4.19803815
## max 3.63624802
## range 7.83428617
## sum -54.33256127
## median -0.06708635
## mean -0.07761794
## SE.mean 0.05170340
## CI.mean.0.95 0.10151257
## var 1.87126905
## std.dev 1.36794337
## coef.var -17.62406067
## skewness -0.08378120
## skew.2SE -0.45343778
## kurtosis -0.07756026
## kurt.2SE -0.21018114
## normtest.W 0.99767494
## normtest.p 0.44324165Se encuentran valores que en su mayoría son positivos y hay una media cercana a 0, que nos indica que el sesgo es hacia la izquierda (sesgo negativo). Para evaluar la significancia y la normalidad, tenemos un p-value (normtest.p) de .443 con un significancia de 95%. Concluimos que cómo .443 > .05, no se rechaza la hipótesis nula y la serie se distribuye de manera normal. Por la parte de la kurtosis obtuvimos un valor de -0.776. A pesar de que es un valor negativo, no contamos con un valor significativo que ayude a concluir que la distribución de la serie es platicúrtica.
Continuaremos con las pruebas de raíz unitaria. En primer lugar, utilizaremos la prueba de Augmented Dickey-Fuller Test (adf.test). El ADF indica que la hipótesis nula es phi = 1, es decir, la serie no es estacionaria, y la hipótesis alternativa nos indica que phi ≠ 1, es decir la serie si es estacionaria. Consideramos, entonces, un nivel de confianza del 95%, si el p-value es mayor a alpha (=0.05), no rechazamos la hipótesis nula y la serie es estacionaria.
Al correr este comando, el programa nos proyecta el siguiente mensaje:
In adf.test(examen) : p-value smaller than printed p-value
Seguido de los valores dentro de la de Augmented Dickey-Fuller Test que es -7.409,. la orden de rezago = 8, y el p-value = 0.01.
Como dice el mismo programa, p-value es menor al nivel de significancia y por ende rechazamos la hipótesis nula y sugerimos que la serie es estacionaria.
.
adf.test(examen)## Warning in adf.test(examen): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: examen
## Dickey-Fuller = -7.409, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary.
En el caso de la prueba Phillips-Perron (*pp.test¨), las hipótesis y la regla de decisión son iguales. Se mantiene el mismo criterio de decisión. Los valores son Dickey-Fuller = -7.409, orden de rezago = 8 y el p-value = 0.01. El p-value es menor al nivel de significancia y por ende rechazamos la hipótesis nula y sugerimos que la serie es estacionaria.
.
pp.test(examen)## Warning in pp.test(examen): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: examen
## Dickey-Fuller Z(alpha) = -282.71, Truncation lag parameter = 6, p-value
## = 0.01
## alternative hypothesis: stationary.
Para la prueba KPSS (del nombre de los autores, Kwiatkowski, Phillips, Schmidt y Shin) (kpss.test) es otra prueba para verificar la estacionariedad de una serie de tiempo. Las hipótesis nula y alternativa para la prueba KPSS son opuestas a las de la prueba ADF. Es decir, H0: El proceso es estacionario en tendencia y H1: La serie tiene una raíz unitaria (la serie no es estacionaria).
Nuestro p-value: a .08151, por lo que se continúa el mismo criterio de decisión.
kpss.test(examen)##
## KPSS Test for Level Stationarity
##
## data: examen
## KPSS Level = 0.3899, Truncation lag parameter = 6, p-value = 0.08151.
Las dos primeras pruebas de estacionariedad sugieren que la serie es estacionaria mientras que la prueba KPSS sugiere que la serie no es estacionaria a un nivel de 0.1, no obstante, se mantiene como estacionaria a un nivel de 0.5 y 0.01.
2.0 Estimación
Gracias al análisis gráfico de la función de autocorrelación, se tiene un primer indicio para esta segunda parte de la metodología. Sobre todo gracias a la función parcial de autocorrelación que nos indica que nuestro modelo debe contener, por lo menos, un MA de orden 4 y minimo un RA de 3 para corregir la alteración de los rezagos. Para el primer modelo se sugiere un modelo ARIMA (3,0,2) con el propósito de generar una alternativa que de mayor peso a los aspectos deterministas que a los estocásticos. Asimismo, se elige este criterio porque las funciones de PACF dan mayor peso estadístico a los componentes autoregresivos. A este modelo se le denominó *”modelo1””.
.
modelo1<-arima(examen,order=c(3,0,2),include.mean=T)
coeftest(modelo1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.492356 0.105715 4.6574 3.203e-06 ***
## ar2 -0.322405 0.089698 -3.5943 0.0003252 ***
## ar3 0.385068 0.060596 6.3546 2.089e-10 ***
## ma1 0.275494 0.111939 2.4611 0.0138511 *
## ma2 0.124480 0.108720 1.1450 0.2522243
## intercept -0.080976 0.118835 -0.6814 0.4956062
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1.
Se descarta incluir el coeficiente y el ma(2) debido a su baja significancia. Por lo tanto, se requiere comprobar si el modelo sigue siendo significativo en sus coeficientes y cumple la condición de convergencia. Por lo tanto, se repite la prueba de significancia individual y se suman los coeficientes asignados a los factores autorregresivos. Cabe destacar que en esta condición solo se suman los coeficientes asignados a los factores autorregresivos, debido a que se está desarrollando un modelo tipo ARIMA, el cual solo requiere dichos coeficientes para comprobar que sus dos primeros momentos sean finitos y constantes.
Prueba de significancia
.
modelo1<-arima(examen,order=c(3,0,1),include.mean=F)
coeftest(modelo1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.554846 0.087845 6.3162 2.68e-10 ***
## ar2 -0.259068 0.071535 -3.6216 0.0002928 ***
## ar3 0.330200 0.039678 8.3219 < 2.2e-16 ***
## ma1 0.214300 0.091654 2.3382 0.0193792 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1.
Prueba de convergencia
modelo1$coef[1]+modelo1$coef[2]+modelo1$coef[3]## ar1
## 0.6259779Se confirma que el modelo arima con p=3 y q=1 cumple con la significancia individual en sus coeficientes y es inferior a 1 en la suma de sus coeficientes autorregresivos. Se obtiene un valor .6259779. El modelo sigue teniendo significancia ya que la suma de los coeficientes del modelo son menores a 1.
plot(modelo1)
El criterio se basa en localizar los puntos que nos muestra la gráfica, si éstos se encuentran al interior de la circunferencia, se cumple la condición de raíces invertibles, por el contrario, si se encuentran afuera o situados sobre los extremos de los ejes, se concluye que no se cumple la condición. Por tanto, confirma el cumplimiento de la condición de raíces invertibles e incluso se evidencía que los parámetros se encuentran alejados de presentar una raíz unitaria.
3.0 Verificación de supuestos
El siguiente paso en la metodología es verificar los supuestos, esto nos permitirá conocer si los modelos cumplen con las restricciones teóricas para poder hacer predicciones de la variable en cuestión, midiendo el comportamiento del componente aleatorio o “ruido blanco”.
3.1 Prueba de Autocorrelación Serial
Se inicia haciendo la prueba de autocorrelación serial conjunta, utilizando a prueba Ljung Box (box.test). Por definición la hipótesis nula señala que los datos se distribuyen de forma independiente (es decir, las correlaciones son 0, por lo que cualquier correlación observada en los datos resulta de la aleatoriedad del proceso de muestreo), y la hipótesis alternativa señala que los datos no se distribuyen de forma independiente y exhiben correlación serial.
.
h<-c(1:length(examen)-1)
for (i in 1:length(examen)-1) {
Q=Box.test(modelo1$residuals, lag=i, type='Ljung-Box')
h[i]<-Q$p.value
}
head(h)## [1] 0.9272350 0.9283267 0.9787708 0.9952582 0.9914426 0.9185487sum(h<.05)## [1] 0.
Se encuentra que el p-value es mayor a alpha en cada uno de sus rezagos, por lo tanto, el modelo cumple el supuesto de no autocorrelación.
La prueba de ARCH Engle (ArchTest) se construye con base a rechazar que los residuos sean heterocedásticos, es decir la hipótesis alternativa, y si los residuos al cuadrado están autocorrelacionados. La prueba consiste en examinar si los cuadrados de los residuos son una secuencia de ruido blanco, que se denomina prueba Q de Portmanteau y es similar a la prueba de Ljung-Box en los cuadrados de los residuos.
.
ArchTest(modelo1$residuals,lags=1)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: modelo1$residuals
## Chi-squared = 0.5592, df = 1, p-value = 0.4546.
En caso de que el p-value sea menor al nivel de significancia, tendríamos que rechazar la hipótesis nula, lo cual indica la existencia de heterocedasticidad. Para este caso, el p-value es mayor al alpha de .05, por lo tanto no rechazamos H0 y decimos que nuestros modelos cumplen, estadísticamente con el supuesto de homocedasticidad.
La prueba de BDS (bds.test) es una de las herramientas más importantes y más utilizadas para la detección de no linealidad en series de tiempo.
.
bds.test(modelo1$residuals)##
## BDS Test
##
## data: modelo1$residuals
##
## Embedding dimension = 2 3
##
## Epsilon for close points = 0.5024 1.0047 1.5071 2.0094
##
## Standard Normal =
## [ 0.5024 ] [ 1.0047 ] [ 1.5071 ] [ 2.0094 ]
## [ 2 ] -1.7861 -1.8569 -1.6694 -1.5402
## [ 3 ] -2.2659 -2.1744 -2.0855 -2.0920
##
## p-value =
## [ 0.5024 ] [ 1.0047 ] [ 1.5071 ] [ 2.0094 ]
## [ 2 ] 0.0741 0.0633 0.095 0.1235
## [ 3 ] 0.0235 0.0297 0.037 0.0364Se observa que en el primer renglón del P-value se dan valores superiores a 1%, lo cual descarta la no linealidad en la serie. De manera que se reafirma la posibilidad de generar pronósticos mediante el uso del modelo.
.
jarque.bera.test(modelo1$residuals)##
## Jarque Bera Test
##
## data: modelo1$residuals
## X-squared = 1.2199, df = 2, p-value = 0.5434.
En la prueba de Jarque Bera (jarque.bera.test) se observa la distribución normal sobre el ruido blanco de los datos. La hipótesis nula señala que el ruido blanco se distribuye de manera normal y la alternativa señala lo opuesto. La implicación de este supuesto es garantizar que los intervalos para los pronósticos sean confiables. Los datos arrojados señalan que el p-value es mayor al nivel de significancia del 5% (.5434 > .05) por lo que no rechazamos la hipótesis nula y concluimos que la distribución del ruido blanco es de manera normal.
Evaluación para el segundo modelo
Para el segundo modelo se sugiere un modelo arima(1,0,4) con el propósito de generar una alternativa que de peso ma’s a los aspectos estocásticos que a los deterministas para evaluar la eficiencia que podrían tener estas aspectos en el modelo como lo sugiere la funcion de autocorrelacion. A este modelo se le denominó “modelo2”.
Los procesos realizados para la estimación, serán los mismos.
.
modelo2<-arima(examen,order=c(1,0,4),include.mean=T)
coeftest(modelo2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.506029 0.125660 4.0270 5.650e-05 ***
## ma1 0.261832 0.123036 2.1281 0.033329 *
## ma2 -0.199883 0.100302 -1.9928 0.046282 *
## ma3 0.150562 0.050781 2.9649 0.003028 **
## ma4 0.193123 0.042614 4.5319 5.846e-06 ***
## intercept -0.080187 0.108049 -0.7421 0.458006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se descarta el uso del coeficiente según sugiere la prueba de significancia
.
modelo2<-arima(examen,order=c(1,0,4),include.mean=F)
coeftest(modelo2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.510028 0.125684 4.0580 4.949e-05 ***
## ma1 0.258727 0.123102 2.1017 0.035578 *
## ma2 -0.202042 0.100509 -2.0102 0.044410 *
## ma3 0.150272 0.050982 2.9476 0.003203 **
## ma4 0.192897 0.042777 4.5094 6.501e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ahora todos los valores resultan significativos por ser menores al p-value.
Para la prueba de convergencia
.
modelo2$coef[1]## ar1
## 0.5100281.
Bajo los criterios mencionados anteriormente se concluye que se cumple la condición de convergencia
.
plot(modelo2)
#### De nuevo se observa que los puntos se contienen dentro del círculo unitario, por lo tanto se cumple la condición de raíces invertibles.
Verificación de supuesto 2 para el modelo
Prueba de no autocorrelación
.
h<-c(1:length(examen)-1)
for (i in 1:length(examen)-1) {
Box.test(modelo2$residuals, lag=1, type = "Ljung-Box")
h[i]<-Q$p.value
}
sum(h<.05)## [1] 0head(h)## [1] 0.8362662 0.8362662 0.8362662 0.8362662 0.8362662 0.8362662La prueba no rechaza la hipótesis nula en ninguno de sus rezagos, por lo tanto se cumple la condición de no autocorrelación serial
.
Utilizando a prueba Ljung Box (box.test). Por definición la hipótesis ´nula señala que los datos se distribuyen de forma independiente (es decir, las correlaciones son 0, por lo que cualquier correlación observada en los datos resulta de la aleatoriedad del proceso de muestreo), y la hipótesis alternativa señala que los datos no se distribuyen de forma independiente y exhiben correlación serial.
.
modelo2$coef[1]## ar1
## 0.5100281Debido a que solo existe un factor autorregresivo se evalúa su que su valor sea menor a 1 en valor absoluto, lo cual se cumple. Por tanto, se cumple la condición de convergencia.
ArchTest
.
ArchTest(modelo1$residuals,lags=1)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: modelo1$residuals
## Chi-squared = 0.5592, df = 1, p-value = 0.4546.
El p-value vuelve a ser mayor a alpha (.4546>-05) y por lo tanto tenemos homocedasticidad.
.
bds.test(modelo2$residuals)##
## BDS Test
##
## data: modelo2$residuals
##
## Embedding dimension = 2 3
##
## Epsilon for close points = 0.5040 1.0079 1.5119 2.0158
##
## Standard Normal =
## [ 0.504 ] [ 1.0079 ] [ 1.5119 ] [ 2.0158 ]
## [ 2 ] -1.6713 -1.7917 -1.588 -1.4090
## [ 3 ] -1.3798 -2.1213 -2.056 -1.9821
##
## p-value =
## [ 0.504 ] [ 1.0079 ] [ 1.5119 ] [ 2.0158 ]
## [ 2 ] 0.0947 0.0732 0.1123 0.1588
## [ 3 ] 0.1676 0.0339 0.0398 0.0475stat.desc(modelo2$residuals, norm=T)## x
## nbr.val 700.00000000
## nbr.null 0.00000000
## nbr.na 0.00000000
## min -3.43368777
## max 3.26502216
## range 6.69870992
## sum -19.21363395
## median -0.01469015
## mean -0.02744805
## SE.mean 0.03809545
## CI.mean.0.95 0.07479521
## var 1.01588408
## std.dev 1.00791075
## coef.var -36.72067074
## skewness -0.01377546
## skew.2SE -0.07455510
## kurtosis 0.25549389
## kurt.2SE 0.69236488
## normtest.W 0.99725088
## normtest.p 0.29065117jarque.bera.test(modelo2$residuals)##
## Jarque Bera Test
##
## data: modelo2$residuals
## X-squared = 2.0676, df = 2, p-value = 0.3557Para la prueba del Backward Digit Span (BDS) se examina la dependencia temporal de la serie. así como los puntos "cercanos ’’ .Mientras más pequeños, más cerca de la dependen. El estadístico de prueba BDS es normal asintóticamente estándar. La primera línea de la matriz de p-value está por encima del 1%, por lo que no rechazamos la H0. Es decir, no existe dependencia temporal.
Los datos arrojados señalan que el p-value en la prueba de Jarque Bera es mayor al nivel de significancia del 5% (.3547 > .05) por lo que no rechazamos la hipótesis nula y concluimos que la distribución del ruido blanco es de manera normal.
4.- Justificación
Iniciamos calculando los valores de BIC. El comando de BIC calcula el criterio de información bayesiano, también conocido como criterio bayesiano de Schwarz (SBC), este nos va a ayudar a comparar la bondad de ajuste de diferentes modelos de regresión.Se pretende elegir el modelo con el valor BIC más bajo como el modelo que mejor se ajusta a los datos.
BIC(modelo1)## [1] 2026.18BIC(modelo2)## [1] 2037.258Para el BIC nos muestra un valor de 2021.8 y 2037.25 que correpsonderia l BIC2. Por lo que el el modelo 1 es mucho mejor ajustando los datos. Es decir se sugiere un modelo ARIMA (3,0,1).
Se realizarán pronosticos in-sample para reforzar el criterio de selección entre los modelos según la precisión que sugiera el criterio MAPE con una relación de datos 80/20.
real <- window(examen, start = 560)
ts.plot(real)
#### En la siguiente gráfica se evidencia la predicción que genero el modelo ARIMA (3,0,1) en la muestra y los valores que sugiere con diferente bandas de confianza
#### MODELO 1
modelo1.train<-arima(examen[1:560], order=c(3,0,1),include.mean=F)
prueba.modelo1 <- forecast(modelo1.train, h = 140, level = c(90,95,99))
plot(prueba.modelo1)
.
MAPE significa error porcentual absoluto medio. La prueba de precisión sugiere que el MAPE se encuentra en 237.22. Nosotros queremos el valor MAPE más pequeño posible ya que eso significaba que la diferencia promedio entre el valor pronosticado y el valor real es el más bajo posible.
Modelo 2
.
La siguiente gráfica presenta las bandas según los intervalos de confianza que estima el modelo ARIMA.
.
#### MODELO 2
modelo2.train<-arima(examen[1:560], order=c(1,0,4),include.mean=F)
prueba.modelo2 <- forecast(modelo2.train, h = 140, level = c(90,95,99))
plot(prueba.modelo2)
La prueba de precisión según MAPE sugiere un valor de 239.57.
accuracy(prueba.modelo2,real)## ME RMSE MAE MPE MAPE MASE
## Training set -0.009497745 1.022149 0.805749 116.81306 239.54713 0.8168982
## Test set -0.273507187 1.199356 0.960393 99.83772 99.83772 0.9736820
## ACF1 Theil's U
## Training set 0.004382661 NA
## Test set 0.530228157 0.9932497Por lo tanto se concluye a partir del criterio de información BIC y la precisión según el criterio de MAPE que el mejor modelo para generar el pronóstico en 300 periodos hacia el futuro es ARIMA (3,0,1). En cuanto a los valores de MAPE el modelo 1 tenía el MAPE más pequeño y se reafirma nuestro supuesto de que el modelo 1 es el más adecuado.
.
Debido a que el propósito de las regresiones anteriores era determinar la precisión de cada dentro de la muestra, no se va a reportar el comportamiento de las bandas en los intervalos de confianza. No obstante, dicho análisis se desarrollará en el apartado de Pronósticos a fuera de la muestra.
PRONÓSTICO
modelo1.train<-arima(examen, order=c(3,0,1),include.mean=F)
pronostico.modelo1 <- forecast(modelo1.train, h = 300, level = c(90,95,99))
plot(pronostico.modelo1)
accuracy(pronostico.modelo1)## ME RMSE MAE MPE MAPE MASE
## Training set -0.02418051 1.004283 0.7951675 92.47113 223.1619 0.8180625
## ACF1
## Training set 0.003444345