Examen Practico Medio Termino
Daniel Sandoval, A000000, 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 se verá dividido en tres 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”. De igual, forma se cargaron previamente todos las librerías necesarias como library(pastecs) , library(tseries), library(lmtest), library(forecast), library(FinTS) y library(TSstudio).
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 analis con el metodo de analisis grafico, estadistica descriptiva y pruebas de raiz unitaria.
1.1 Análisis Gráfico
ts.plot(examen)
Gráfica 1
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 metodologia usada en este pryecto que supone la regresión de la media, a partir del limite central se descarta que existan fenomenos de tendencia que distorsionen el comortamiento 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, asi como valores minimos 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 le componente autoregresivo se concentra principalmente hasta p=3. Es importante mencionar que uno de los supuestos de Media Moviles es que el componente es negativo al corregir el componente determinista del modelo. Por lo tanto, la grafica 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 de la estadística descriptiva donde evaluaremos principalmente valores como el sesgo y la curtosis con una demostración matemática, a pesar de ya contar con un analisis 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.44324165Tabla 1
Se 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 por lo tanto 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 platicurtica.
Continuaremos con las pruebas de raiz unitaria seguimos con las pruebas de raíz unitaria. En primer lugar, utilizaremos a 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:
Warning message:
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).
Con esta prueba queremos que el pvalue > al nivel de significancia (.05) para no rechazar H0.
Nuestro p-value: a .08151, por lo que se continua 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. Por lo tanto se continua una serie diferenciada de orden 1, mas no se descarta la validez de usar la serie de nivel.
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 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 significante en sus coeficientes y cumple la condicion de converegenecia. 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 condicion solo se suman los coeficientes asignados a los factores autoregresivos, debido a que se esta 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 coeificientes y es inferior a 1 en la suma de sus coeficientes autoregresivos.
Se obtiene un valor .839937. El modelo sigue teniendo significancia ya que la suma de los coeficientes del modelo son menores a 1.
La prueba de raices invertibles indica si los modelos tienen o no raíces invertibles a traves de una grafica.
plot(modelo1)
El criterio se basa en localizar los puntos, si éstos se encuentran al interior de la circunferencia, se dice que cumple con la condicion de raices invertibles, por el contrario, si se encuentran afuera o situados sobre los extremos de los ejes, se concluye no cumple la condicion. Por tanto, confirma el cumplimiento de la condicion de raices ivertibles e incluso se evidencia que los parametros se encuentran alejados de presentar raiz 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 incia 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 resagos, 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 heteroscedásticos, es decir la hipotesis 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.0364. #### Se observa que en el primer renglon del P-value se dan valores superiores a 1%, lo cual descarta la no linearidad en la serie. De manera que se reafirma la posibilidad de generar pronosticos 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 rudio 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 proposito de generar una alternativa que de peso ma’s a los aspectos estocasticos que a los deterministas para evaluar la eficiencia que podrian 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, seran 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 ' ' 1.
Se descarta el uso del coeificiente segun 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 ' ' 1.
Ahora todos los valores resultan significantes 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 condicion de convergencia
.
plot(modelo2)
#### De nuevo se observa que los puntos se contienen dentro del circulo unitario, por lo tanto se cumple la condicion de raices invertibles.
Verificacion de supuesto 2 para el modelo
Prueba de no autocorrelacion
.
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 hipotesis nula en ninguno de sus rezagos, por lo tanto se cumple la condicon de no autocorrelacion serial
.
Utilizando a prueba Ljung Box (box.test). Por definición la hiptesis ´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.5100281.
Debido a que solo existe un factor autorregresivo se evalua su que su valor sea menor a 1 en valor absoluto, lo cual se cumple. Por tanto, se cumple la condicion 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.
DESGLOSAR RESULTADOS DE CADA PRUEBA Y BREVE RESUMEN DE CADA UNA
.
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 Digir Span (BDS) se a examina la dependencia temporal de la serie. asi 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.
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
Iniciariamos 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.
print('Modelo 1')## [1] "Modelo 1"BIC(modelo1)## [1] 2026.18print('Modelo 2')## [1] "Modelo 2"BIC(modelo2)## [1] 2037.258Para el bic1 nos muestra un valor de 2021.8 y 2037.25. Por lo que el el modelo 1 es mucho mejor ajustando los datos. Es decir se sugiere un modelo ARIMA (3,0,1)
Se realizaran pronosticos in-sample para reforzar el criterio de seleccion entre los modelos segun la precision que sugiera el criterio MAPE con una relacion de datos 80/20.
A continuacion se grafica la serie de prueba
real <- window(examen, start = 561)
ts.plot(real)
#### En la siguiente grafica se evidencia la prediccion que genero el modelo arima (3,0,1) en la muestra y los valores que sugiere con diferente bandas de confianza # ESCRIBIR BANDAS EN DONDE CONVERGEN LOS PRONOSTICOS, SE VEN CON EL COMANDO prueba.modelo1
#### 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)
#lines(real, col="red").
MAPE significa error porcentual absoluto medio. La prueba de precision sugiere que el MAPE se encuentra en 237.22. Nosotros queremos el valor MAPE más pequeño posible ya que eso significara que la diferencia promedio entre el valor pronosticado y el valor real es el más bajo posible..
.
accuracy(prueba.modelo1,real)## ME RMSE MAE MPE MAPE MASE
## Training set -0.008760919 1.016773 0.8023825 110.68651 237.22342 0.8134851
## Test set -0.267405867 1.198649 0.9591855 99.70946 99.70946 0.9724578
## ACF1 Theil's U
## Training set 0.005466997 NA
## Test set 0.530939840 0.9930925Modelo 2
###La siguiente grafica presenta las bandas segun 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 precision segun 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 de criterio de infromacion BIC y la precision segun el crtierio de MAPE que el mejor modelo para generar el pronostico en 300 periodos hacia el futuro es arima (3,0,1). En cuanto a los valores de MAPE el modelo 1 tenia el MAPE más pequeño y se reafirma nuestro supuesto de que el modelo 1 es el más adecuado.
.
PRONOSTICO
La grafica presenta que (REPORTAR A QUE GRADO DE CONFIANZA SE CONVERGE A CUALES BANDAS)
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