Tópicos de Econometría II - Sesión 3
Tópicos de Econometría II - Sesión 3
1. Introducción
Esta sesión corresponde al procedimiento para realizar un modelo de pronóstico ARIMA con el lenguaje de programación R. Para esto se aplica la metodología de Box-Jenkins con la transformación Box-Cox de Victor Guerrero (2003). Aunque el modelo ARIMA pertenece al campo de la estadística paramétrica, se aplican las pruebas Anderson-Darling y el criterio de información bayesiano (BIC).
El objetivo del trabajo es hacer un pronóstico doce meses adelante de la muestra con un modelo ARIMA(p,d,q). Todas las pruebas tienen una confianza del 95%.
a. Descripción de los datos
La serie de tiempo (\(ts\)) es el precio del barril del petróleo WTI. Los datos son numéricos discretos de razón, tienen una periodicidad mensual con un tamaño muestral de 324 observaciones, inicia en enero del año 1993 y finaliza en diciembre del 2019. La serie se obtiene de la página web de la Administración Energética de los Estados Unidos (EIA). Para este ejercicio los datos de la \(ts\) están en formato PDF.
b. Extracción de tablas en un archivo con formato PDF.
Con la librería \(tabulizer\) aplicamos la función \(extract\_tables()\) para poder importar las tablas contenidas dentro de las ocho páginas que hay en el PDF.
base_1 <- extract_tables("precio_wti.pdf",
output = "matrix")La función arroja una lista (base_1) que contiene ocho elementos de clase matriz.
Estableciendo un vector de nombres para los elementos de la lista
Con la función \(lapply()\) podemos aplicar un vector de caracteres, que contiene los nombres de las columnas, en la primera fila de cada una de las matrices contenidas en la lista.
base_1 <- lapply(base_1, setNames,
c("fecha", "precio_wti"))Uniendo las matrices mediante las filas
El paquete \(dplyr\) nos permite llamar la función \(bind\_rows()\) que es una implementación más eficiente respecto a la función base \(rbind()\).
base_1 <- bind_rows(base_1)Estudiamos la clase de las columnas del objeto base_1
sapply(base_1, class)## ...1 ...2
## "character" "character"La función \(sapply()\) es una simplificación de \(lapply()\). Esta permite una mejor implementación de \(apply\) para objetos que no son listas. Ambas columnas (fecha y precio) contienen caracteres y tienen como título …1 y …2 respectivamente.
Cambiamos los nombres de las columnas de base_1
names(base_1) <- matrix(base_1[1, ])Con la función \(names()\) podemos cambiar los nombres del objeto. Le indicamos que los nuevos nombres de las columnas sean los caracteres que están en la primera de la matriz.
Elimino la primera fila de base_1
base_1 <- base_1[-1, ]Como los títulos del objeto ya corresponden a primera fila procedemos a removerla.
c. Implementando un data-frame
Es necesario convertir completamente \(base\_1\) a un objeto de clase data-frame (df). Para esto se usa la función base \(as.data.frame()\).
base_1 <- as.data.frame(base_1)Primeras seis filas del df
La función \(head()\) permite observar en la consola las seis primeras filas del df estudiado
head(base_1)## fecha precio_wti
## 1 ene-1993 19,03
## 2 feb-1993 20,09
## 3 mar-1993 20,32
## 4 abr-1993 20,25
## 5 may-1993 19,95
## 6 jun-1993 19,09Podemos ver que en la segunda columna el separador de decimales es una coma por lo que debemos cambiarlo a un punto.
Columna fecha en formato Date
En la primera columna es necesario crear una secuencia de fechas mensuales que inicie el primero de enero de 1993.
base_1$fecha <- seq(as.Date("1993-01-01"),
length.out = nrow(base_1),
by = "month")Convertir la segunda columna del df a clase numérica
Cambiamos el caracter coma (,) con un punto (.) a través de la función \(gsub\), transformamos la columna a formato numérico y por último convertimos toda la columna a una serie de tiempo mensual.
base_1$precio_wti <- ts(as.numeric(gsub(",", ".", base_1$precio_wti)),
start = c(1993, 01), frequency = 12)Ahora ya tenemos un df con fechas (Date) y el precio del barril como una \(ts\) mensual.
## fecha precio_wti
## 1 1993-01-01 19.03
## 2 1993-02-01 20.09
## 3 1993-03-01 20.32
## 4 1993-04-01 20.25
## 5 1993-05-01 19.95
## 6 1993-06-01 19.092. Metodología Box-Jenkins
a. Identificación
Es necesario identificar si la \(ts\) del precio mensual del barril WTI es estacionaria. Lo anterior implica que la media y la varianza son constantes a través del tiempo, además que la covarianza entre dos periodos de tiempo dependa únicamente de la distancia entre estas (Villavicencio, 2014).
Primero vamos a realizar la gráfica de la serie de tiempo estudiada, por medio de los paquetes ggplot2 y scales. El siguiente es el código de la gráfica:
graf_5 <- ggplot(data = base_1, aes(x = fecha,
y = precio_wti)) +
geom_line(col="firebrick2", size=1) +
ylab("Dólares americanos") + xlab("") +
scale_y_continuous(breaks = seq(5, 140, 10))+
theme(axis.text.x = element_text(angle = 90,
size = 14),
plot.title = element_text(hjust = 0.5,
size = 19),
axis.title.y = element_text(size = 14),
axis.text.y = element_text(size = 14)) +
scale_x_date(breaks = date_breaks("1 year"),
expand = c(0,0), date_labels = "%Y") +
xlab("Años") +
ggtitle("Precio mensual del barril de petróleo WTI \n (1993-2019)")
graf_5
La gráfica de la ts muestra volatilidad en el precio mensual del barril WTI. Aparentemente la ts no tiene presencia de estacionalidad, sin embargo, vamos a comprobar con métodos gráficos y pruebas de hipótesis si esta ts tiene un comportamiento estacional.
Descomposición de la ts
La librería \(ggseas\) es un complemento de series de tiempo para \(ggplot2\). La función \(ggsdc()\) arroja las gráficas de tendencia, estacionalidad e irregularidad de la ts estudiada.
ggsdc(base_1, mapping = aes(x = fecha, y = precio_wti),
method = "decompose", start = c(1993, 01),
frequency = frequency(base_1$precio_wti)) +
geom_line()
La ts del WTI no tiene tendencia y en la irregularidad podemos ver la alta volatilidad del precio mensual del barril de petróleo.
Gráfica de estacionalidad con coordenadas polares
En el paquete \(forecast\) está la función \(ggseasonplot\) que tiene la opción de graficar una ts, en coordenadas polares, versus su respectiva periodicidad y de esta manera examinar si la serie presenta estacionalidad.
ggseasonplot(base_1$precio_wti, polar = TRUE) +
ggtitle("Coordenadas polares") +
theme(plot.title = element_text(hjust = 0.5,
size = 14)) +
xlab("Meses")
La gráfica de coordenadas polares no muestra comportamiento estacional en ninguno de los doce meses.
Método de las variables dummy para detectar estacionalidad
Aplicamos una regresión con variables dummy que tiene el mismo proceso que está en la guía número dos. La salida de la regresión lineal es la siguiente:
##
## Call:
## lm(formula = base_1$precio_wti ~ ., data = base_dummies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.389 -15.111 -3.913 13.791 76.422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.33730 4.65077 2.868 0.00442 **
## tendencia 0.21566 0.01292 16.686 < 2e-16 ***
## dummy_2 0.41693 5.91813 0.070 0.94388
## dummy_3 2.28904 5.91817 0.387 0.69918
## dummy_4 3.74782 5.91824 0.633 0.52703
## dummy_5 3.97290 5.91834 0.671 0.50254
## dummy_6 4.00686 5.91847 0.677 0.49890
## dummy_7 4.47564 5.91862 0.756 0.45010
## dummy_8 3.56331 5.91881 0.602 0.54759
## dummy_9 3.11616 5.91902 0.526 0.59894
## dummy_10 1.92198 5.91926 0.325 0.74563
## dummy_11 0.21520 5.91953 0.036 0.97102
## dummy_12 -1.35157 5.91982 -0.228 0.81955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.74 on 311 degrees of freedom
## Multiple R-squared: 0.4746, Adjusted R-squared: 0.4543
## F-statistic: 23.41 on 12 and 311 DF, p-value: < 2.2e-16Tomando como fijo la dummy del mes de enero (\(dummy\_1\)), ninguno de los otros once meses es significativo en la regresión con una confianza del 95%. A partir de estos resultados es posible concluir que la \(ts\) no presenta estacionalidad en ningún mes.
Prueba Dickey-Fuller aumentada (urca) para la serie original
Es necesario identificar si la \(ts\) del precio mensual del barril WTI es estacionaria. Una serie estacionaria implica que la media y la varianza son constantes a través del tiempo, además que la covarianza entre dos periodos de tiempo dependa únicamente de la distancia entre estas (Villavicencio, 2014).
El test de Dickey-Fuller aumentado permite examinar si una \(ts\) presenta raíz unitaria. La hipótesis nula (\(H_{0}\)) de esta prueba es raíz unitaria en la \(ts\). La función \(ur.df()\), de la librería \(urca\), implementa la prueba en la \(ts\) que se desea examinar.
summary(ur.df(base_1$precio_wti,
type = "none"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.5176 -1.6161 0.2519 2.5379 14.8033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.003456 0.004183 -0.826 0.409
## z.diff.lag 0.368041 0.052089 7.066 1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.39 on 320 degrees of freedom
## Multiple R-squared: 0.1354, Adjusted R-squared: 0.13
## F-statistic: 25.05 on 2 and 320 DF, p-value: 7.792e-11
##
##
## Value of test-statistic is: -0.8262
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62Como el valor crítico de la prueba está ubicado a la izquierda del estadístico no hay suficiente evidencia estadística para rechazar \(H_{0}\), con una significancia del 5%. Se concluye que la ts presenta raíz unitaria y no es estacionaria.
Prueba Dickey-Fuller aumentada (tseries) para la serie original
El paquete \(tseries\) tiene la función \(adf.test()\) que también corre la prueba Dickey-Fuller aumentada.
adf.test(base_1$precio_wti)##
## Augmented Dickey-Fuller Test
##
## data: base_1$precio_wti
## Dickey-Fuller = -2.1384, Lag order = 6, p-value = 0.5181
## alternative hypothesis: stationaryComo el \(p-value\) de la prueba es mayor que la significancia elegida (0,05), no se rechaza \(H_{0}\).
Aplicando la transformación Box-Cox (método de Víctor Guerrero)
Para poder estabilizar la \(ts\) en varianza vamos a aplicar la transformación boxcox con el método de Guerrero en donde se selecciona el lambda que minimiza el coeficiente de variación para la \(ts\) estudiada (Guerrero, 2003).
Lambda para la ts estudiada
Vamos a guardar el valor del lambda en un objeto para poder usarlo posteriormente.
lambda_1 <- BoxCox.lambda(base_1$precio_wti)
lambda_1## [1] -0.1655381Transformación Box-Cox
Con la función \(BoxCox()\) transformamos la \(ts\) original e indicamos el lambda seleccionado anteriormente.
base_1$bcx <- BoxCox(base_1$precio_wti,
lambda = lambda_1)Gráfica de la serie transformada con Box-Cox
graf_6 <- ggplot(base_1, aes(x = fecha,
y = bcx))+
geom_line(col="forestgreen", size=1) +
scale_y_continuous(breaks = seq(0, 10, 0.1)) +
ylab("Escala Box-Cox") + xlab("Años") +
theme(axis.text.x = element_text(angle = 90),
plot.title = element_text(hjust = 0.5,
size = 21)) +
scale_x_date(breaks = "1 year", expand = c(0,0),
date_labels = "%Y") +
ggtitle("Transformación Box-Cox para el precio \n del barril WTI")
graf_6
Esta nueva \(ts\) está ajustada en varianza. A continuación se va a aplicar la prueba Dickey-Fuller aumentada para poder determinar si la transformación Box-Cox es estacionaria.
Prueba Dickey-Fuller aumentada (urca) para la transformación Box-Cox
summary(ur.df(base_1$bcx, type = "none"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.141274 -0.027457 0.004137 0.029262 0.136658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0003362 0.0008667 0.388 0.698
## z.diff.lag 0.2367040 0.0543290 4.357 1.78e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04335 on 320 degrees of freedom
## Multiple R-squared: 0.057, Adjusted R-squared: 0.05111
## F-statistic: 9.671 on 2 and 320 DF, p-value: 8.352e-05
##
##
## Value of test-statistic is: 0.3879
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62Como el valor crítico es menor que el estadístico de prueba no hay suficiente evidencia estadística para rechazar \(H_{0}\), con una confianza del 95%. Se concluye que la serie transformada por Box-Cox no es estacionaria.
Prueba Dickey-Fuller aumentada (tseries) para la transformación Box-Cox
adf.test(base_1$bcx)##
## Augmented Dickey-Fuller Test
##
## data: base_1$bcx
## Dickey-Fuller = -2.0403, Lag order = 6, p-value = 0.5595
## alternative hypothesis: stationaryComo el p-valor es mayor que la significancia elegida (0,05) no se puede rechazar \(H_{0}\). Por lo tanto la \(ts\) transformada con Box-Cox presenta raíz unitaria.
Aplicando primera diferencia
Para establizar la serie Box-Cox en media es necesario la aplicación de la primera diferencia con la pérdida del primer dato de la \(ts\). Con la función \(diff()\) se implementa la primera diferencia a la serie estudiada.
diff_1 <- diff(base_1$bcx)Será necesario construir un nuevo df que tenga la nueva serie diferenciada con sus respectivo vector de fechas.
Segundo data-frame
El nuevo df tiene el nombre base_dif con la pérdida de la primera fecha de la \(ts\) original.
base_dif <- data.frame(fecha_1 = base_1$fecha[-1], diff_1)El objeto base_dif ya permite graficar la primera diferencia con el paquete ggplot2.
graf_7 <- ggplot(base_dif, aes(x = fecha_1,
y = diff_1)) +
geom_line(col="blueviolet", size = 0.7) +
scale_y_continuous(breaks = seq(-0.3, 0.3, 0.01)) +
xlab("Años") + ylab("Variación porcentual") +
theme(axis.text.x = element_text(angle = 90),
plot.title = element_text(hjust = 0.5,
size = 21)) +
scale_x_date(breaks = date_breaks("1 year"),
expand = c(0,0),
date_labels = "%Y") +
ggtitle("Primera diferencia de la serie transformada")
graf_7
La tercera gráfica tiene estabilidad alrededor de cero. A continuación se aplica, de nuevo, la prueba Dickey-Fuller aumentado.
Prueba Dickey-Fuller aumentada (urca) para la primera diferencia
summary(ur.df(base_dif$diff_1, type = "none"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.136927 -0.027191 0.005359 0.030179 0.135317
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.74032 0.06915 -10.706 <2e-16 ***
## z.diff.lag -0.02870 0.05597 -0.513 0.608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04341 on 319 degrees of freedom
## Multiple R-squared: 0.3814, Adjusted R-squared: 0.3775
## F-statistic: 98.34 on 2 and 319 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -10.7064
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62En este caso el valor del estadístico es menor que el del valor crítico se rechaza la \(H_{0}\) y se concluye que la \(ts\) de la primera diferencia es estacionaria con una significancia del 5%.
Prueba Dickey-Fuller aumentada (tseries) para la primera diferencia
adf.test(base_dif$diff_1)## Warning in adf.test(base_dif$diff_1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: base_dif$diff_1
## Dickey-Fuller = -7.247, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationaryEl p-valor es menor que la significancia elegida (0,05), entonces se rechaza \(H_{0}\) y se concluye que la serie de la primera diferencia no tiene raíz unitaria.
b. Estimación de los parámetros
Para elegir el modelo ARIMA(p, d, q) más adecuado para el pronóstico es necesario usar los criterios de información de Akaike y Schwarz para elegir el valor de los rezagos \(p\) y \(q\).
Loop con los criterios de información AIC y SBC
Usamos un bucle con los criterios de información con valores de cero a dos para los rezagos:
resultados <- c("p", "d", "q", "AIC", "SBC")
for (i in 0:2) {
for (j in 0:2) {
prueba <- arima(base_1$bcx, order = c(i, 1, j))
resultados <- rbind(resultados, as.numeric(c(i, 1, j,
AIC(prueba),
BIC(prueba))))
}
}
resultados## [,1] [,2] [,3] [,4] [,5]
## resultados "p" "d" "q" "AIC" "SBC"
## "0" "1" "0" "-1092.27530313372" "-1088.4976508105"
## "0" "1" "1" "-1106.90060897063" "-1099.34530432418"
## "0" "1" "2" "-1108.0774484474" "-1096.74449147773"
## "1" "1" "0" "-1109.02217915318" "-1101.46687450673"
## "1" "1" "1" "-1107.19338046354" "-1095.86042349387"
## "1" "1" "2" "-1106.07746486857" "-1090.96685557568"
## "2" "1" "0" "-1107.28670070654" "-1095.95374373687"
## "2" "1" "1" "-1105.5380878322" "-1090.42747853931"
## "2" "1" "2" "-1105.08470886242" "-1086.1964472463"Los resultados del AIC y SBC muestran que la combinación 1,1,0 es la más adecuada para un modelo ARIMA(p,d,q).
Alternativa con la función lapply para los criterios de información AIC y SBC
Las funciones apply son una alternativa a los loops en R. Implementar estas funciones en lugar de un bucle puede ahorrar memoria RAM del procesador.
Creamos tres vectores: p, d, q
p <- 0:2
d <- 1
q <- 0:2Con la función \(expand.grid()\) generamos un nuevo df (df_apply) con todas las combinaciones posibles de p,d,q.
df_apply <- expand.grid(p, d, q)Ahora voy a cambiar el nombre de las tres columnas de df_apply
colnames(df_apply) <- c("p", "d", "q")Ahora vamos a generar una lista que contiene todas las combinaciones de p,d,q. Con la función \(lapply()\) separamos las filas.
tester <- lapply(split(df_apply, seq(NROW(df_apply))),
unlist)Ya podemos crear un df con el AIC y SBC correspondiente a las posibles combinaciones de p y q (con la primera diferencia) para un modelo ARIMA.
tester <- lapply(tester,
function(x) arima(base_1$bcx,
order = x) %>%
glance) %>% do.call(rbind, .)Pegamos el df tester con el anterior df_apply. Después ordenamos las filas por los criterios de información:
df_apply <- cbind(df_apply, tester[, 3:4])
df_apply[order(df_apply$AIC), ]## p d q AIC BIC
## 2 1 1 0 -1109.022 -1101.467
## 7 0 1 2 -1108.077 -1096.744
## 3 2 1 0 -1107.287 -1095.954
## 5 1 1 1 -1107.193 -1095.860
## 4 0 1 1 -1106.901 -1099.345
## 8 1 1 2 -1106.077 -1090.967
## 6 2 1 1 -1105.538 -1090.427
## 9 2 1 2 -1105.085 -1086.196
## 1 0 1 0 -1092.275 -1088.498El df con la función lapply y la matriz que arroja el loop nos muestran que el modelo más adecuado de predicción, según los criterios de información, es un ARIMA(1,1,0).
Implementando el modelo ARIMA(1,1,0)
Con la función \(arima()\) guardamos el modelo de predicción con el orden de los rezagos seleccionados.
modelo_1 <- arima(base_1$bcx,
order = c(1, 1, 0))
coeftest(modelo_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.237402 0.054011 4.3954 1.106e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El paquete \(lmtest\) tiene disponible la función \(coeftest()\) que nos permite comprobar que el rezago 1 de p es significativo con una confianza del 95%.
Podemos graficar la raíz inversa del modelo ARIMA(1,1,0)
autoplot(modelo_1)
La raíz inversa del modelo está dentro del círculo unitario.
c. Verificación de los supuestos
Normalidad de los residuos del ARIMA(1,1,0)
Vamos a graficar el histograma de los residuos del modelo. Además, esta gráfica contiene una línea recta para la media y mediana de los residuos.
hist(modelo_1$residuals, col = "royalblue",
xlab = "Residuos", ylab = "Frecuencia",
main = "Histograma de los residuos")
abline(v = mean(modelo_1$residuals),
col = "firebrick2", lwd = 3)
abline(v = median(modelo_1$residuals),
col = "darkgoldenrod2", lwd = 3)
legend("topright", legend = c("Media", "Mediana"),
col = c("firebrick2", "darkgoldenrod2"), lwd = c(3, 3))
El histograma de los residuos tiene aparentemente una distribución normal y las dos medidas de tendencia central tienen valores cercanos entre ellas.
Con el paquete \(car\) podemos obtener la gráfica Q-Q con un intervalo de confianza (IC) del 95%.
car::qqPlot(modelo_1$residuals, xlab = "Cuantiles",
ylab = "Residuos",
main = "QQ-plot de los residuos")
## [1] 192 75La gráfica QQ muestra que la mayoría de los residuos están dentro del IC seleccionado y solo dos residuos (75 y 192) pueden considerarse como datos atípicos.
Aplicamos cuatro pruebas de hipótesis para comprobar si los residuos tienen una distribución normal: Anderson-Darling (AD), Lilliefors (Lill), Shapiro-Wilk (SW) y Jarque-Bera (JB). Creamos una lista que contiene las funciones de estas cuatro pruebas y posteriormente extraemos únicamente el p-valor de cada una de estas:
normalidad <- function(x){
list(AD = ad.test(x), Lill = lillie.test(x),
SW = shapiro.test(x), JB = jarque.bera.test(x))
}
sapply(normalidad(modelo_1$residuals), "[[", "p.value")## AD Lill SW JB
## 0.08513986 0.14881278 0.10067253 0.05067999Como en las cuatro pruebas el p-valor es mayor que el nivel de significancia escogido (0,05) no se rechaza \(H_{0}\) y concluimos que los residuos siguen una distribución normal.
No autocorrelación de los residuos
Graficamos la Función de autocorrelación (ACF) y la Función de autocorrelación parcial (PACF) de los residuos del modelo.
par(mfrow = c(1, 2))
acf(modelo_1$residuals, xlab = "Rezagos",
main = "Función de autocorrelación")
pacf(modelo_1$residuals, xlab = "Rezagos", main = "Función de autocorrelación parcial")
Como en ambas gráficas las líneas verticales no transpasan los límites horizontales podemos concluir que los residuos del modelo no presentan autocorrelación con una significancia del 5%.
Ahora aplicaremos la prueba de Ljung-Box para los primeros 25 residuos del modelo tal y como recomienda Tsay (2010). La hipótesis nula (\(H_{0}\)) es no hay presencia de autocorrelación. En primer lugar vamos a aplicar este test con un loop:
columnas <- c("Rezago", "p_value")
for (k in 1:25) {
test <- Box.test(modelo_1$residuals,
type = "Ljung-Box", lag = k)
columnas <- rbind(columnas,
as.numeric(c(k, test$p.value)))
}
columnas ## [,1] [,2]
## columnas "Rezago" "p_value"
## "1" "0.89105877241453"
## "2" "0.758067236090059"
## "3" "0.785829350084746"
## "4" "0.677413011061437"
## "5" "0.802132536187996"
## "6" "0.865757478082586"
## "7" "0.875215113300186"
## "8" "0.847781601855228"
## "9" "0.86708196965185"
## "10" "0.857779595931882"
## "11" "0.719661364639362"
## "12" "0.786986754602244"
## "13" "0.608938855608978"
## "14" "0.668425615639136"
## "15" "0.57710820735862"
## "16" "0.585779023040026"
## "17" "0.654277334389172"
## "18" "0.713532971611072"
## "19" "0.769223192429207"
## "20" "0.758937404251339"
## "21" "0.795170909487892"
## "22" "0.788439199961225"
## "23" "0.752501906983909"
## "24" "0.658029225772182"
## "25" "0.592632680894079"Ninguno de los p-value es menor que la significancia del modelo, por lo que no se puede rechazar la \(H_{0}\) y se concluye que no hay autocorrelación en los primeros 25 residuos del modelo ARIMA(1,1,0).
Prueba Ljung-Box con la familia de funciones apply
Con la función \(lapply()\) podemos pasar el test Ljung-Box a los primeros 25 residuos estudiados
lb_alt <- lapply(1:25,
function(x) Box.test(modelo_1$residuals,
type = "Ljung-Box",
lag = x))Ahora creamos un df extrayendo los p-valores de la lista anterior con la función \(sapply()\)
lb_alt <- data.frame(Rezago = 1:25,
P_value = sapply(lb_alt,
"[[", "p.value"))
lb_alt## Rezago P_value
## 1 1 0.8910588
## 2 2 0.7580672
## 3 3 0.7858294
## 4 4 0.6774130
## 5 5 0.8021325
## 6 6 0.8657575
## 7 7 0.8752151
## 8 8 0.8477816
## 9 9 0.8670820
## 10 10 0.8577796
## 11 11 0.7196614
## 12 12 0.7869868
## 13 13 0.6089389
## 14 14 0.6684256
## 15 15 0.5771082
## 16 16 0.5857790
## 17 17 0.6542773
## 18 18 0.7135330
## 19 19 0.7692232
## 20 20 0.7589374
## 21 21 0.7951709
## 22 22 0.7884392
## 23 23 0.7525019
## 24 24 0.6580292
## 25 25 0.5926327Para el df de las funciones apply la conclusión es exactamente la misma que la del loop.
Gráficos con la librería forecast
La función \(checkresiduals()\) arroja algunos gráficos que realizamos arriba:
checkresiduals(modelo_1$residuals)
d. Uso del modelo para pronóstico
Para pronosticar usamos la función \(forecast\), en esta indicamos que el pronóstico sea doce meses adelante y que haga dos IC (80% y 95% de significancia).
pronostico <- forecast(modelo_1, h = 12,
level = c(0.80, 0.95))Para poder tener un pronóstico con las unidades de la \(ts\) inicial (dólares americanos) vamos a usar la función \(InvBoxCox()\) en: los datos pronosticados, la serie histórica y los límites superior e inferior del pronóstico.
pronostico$x <- InvBoxCox(pronostico$x,
lambda = lambda_1)
pronostico$mean <- InvBoxCox(pronostico$mean,
lambda = lambda_1)
pronostico$upper <- InvBoxCox(pronostico$upper,
lambda = lambda_1)
pronostico$lower <- InvBoxCox(pronostico$lower,
lambda = lambda_1)Con la escala original ya podemos hacer la gráfica del pronóstico
autoplot(base_1$precio_wti, series = "WTI") +
autolayer(pronostico,
series = "Pronóstico", showgap = FALSE) +
scale_colour_manual("", values = c("royalblue",
"firebrick2")) +
xlab("Años") + ylab("Dólares americanos") +
scale_y_continuous(breaks = seq(5, 140, 5)) +
theme(plot.title = element_text(hjust = 0.5, size = 21)) +
ggtitle("Pronóstico del precio mensual \n del barril WTI")
3. Referencias
Guerrero, V. (2003). Análisis estadístico de series de tiempo económicas. México D.F. México: editorial Thompson.
Tsay, R. (2010). Analysis of Financial Time Series (Third edition). New Jersey. Estados Unidos de América: John Willie & Sons.
Quintana, L. y Mendoza, M. (2017). Econometría aplicada utilizando R. Ciudad de México. México: Universidad Nacional Autónoma de México (UNAM). Recuperado de: http://www.saree.com.mx/unam/EbR.
Villavicencio, J. (2014). Introducción a series de tiempo. Recuperado de: http://www.estadisticas.gobierno.pr/iepr/LinkClick.aspx?fileticket=4_BxecUaZmg%3D.
4. Comentarios
Las base de datos y el código están disponibles en el respositorio de GitLab: https://gitlab.com/fevidals/topicos-de-econometria-ii