Tarea 3 | Series de tiempo
Modelación ARIMA
Prieto Palacios Mauricio Alejandro
NOTA: La mayoría de los gráficos son interactivos.
Primero se carga los paquetes necesarios.
require(tidyverse)
require(lubridate)
require(scales)
require(janitor)
require(ggsci)
require(plotly)
require(timetk)
require(modeltime)
require(sweep)
require(tidymodels)
require(forecast)
require(tsibble)
require(fable)
require(feasts)
require(lmtest)
require(nortest)
require(tseries)
require(imputeTS)
require(tsdl)Se establece el tema y los colores.
theme_set(theme_bw() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)))color_prin <- "#008080" # Teal
color_sec <- "#C71585" # MediumVioletRed1 Análisis descriptivo
Se carga los datos:
datos_ts <- tsdl[[12]]Se opta por trabajar con tibble y tsibble:
datos <- datos_ts %>%
tk_tbl() %>%
rename(fecha = index, gasto = "value") %>%
mutate(fecha = lubridate::yq(fecha), año = lubridate::year(fecha), trimestre = lubridate::quarter(fecha)) %>%
relocate(gasto, .after = trimestre) %>%
mutate(año = as_factor(año), trimestre = as_factor(trimestre))
datos_tsb <- datos_ts %>%
as_tsibble() %>%
rename(fecha = index, gasto = value)
datosg_datos <- datos %>%
ggplot(aes(fecha, gasto)) + geom_area(fill = color_prin, alpha = 0.1) + geom_line(colour = color_prin) +
scale_x_date(breaks = scales::breaks_width("year"), labels = scales::label_date_short()) +
scale_y_continuous(labels = scales::dollar) + labs(x = "Fecha", y = "Gasto",
caption = "Fuente: Abraham & Ledolter (1983)", title = "Gastos en nuevo equipo de planta en E.U.",
subtitle = "Observaciones trimestrales")
g_datos %>%
ggplotly() %>%
config(displayModeBar = FALSE)Parece haber un patrón cada 4 observaciones, es decir, hay un ciclo estacionario de periodo \(s=4\) y esto es debido a que en todos los años los gastos en el primer trimestre son bajos, para el segundo hay un aumento considerable, en el tercero se mantienen estables y finalmente en el cuarto alcanzan su máximo.
Para confirmar esto se muestra el comportamiento de las observaciones a lo largo de los años de manera separada.
g_datos_anual <- datos %>%
ggplot(aes(trimestre, gasto, group = año, colour = año)) + geom_line(alpha = 0.7) +
labs(title = "Gastos por año", x = "Trimestre", y = "Gastos", colour = "Año") +
scale_y_continuous(labels = scales::dollar)
g_datos_anual %>%
ggplotly() %>%
config(displayModeBar = FALSE)Adicionalmente, se puede apreciar que conforme pasan los años los gastos aumentan, por lo que se afirma que hay una tendencia lineal a la alza.
También se muestra el ACF y el PACF de la serie.
datos_tsb %>%
feasts::gg_tsdisplay(gasto, plot_type = "partial") + ggtitle("Serie de tiempo original con correlogramas")
Parece haber una varianza creciente en los datos, para analizar esto esto se realiza la prueba de Breusch-Pagan.
bptest(datos$gasto ~ seq(1, length(datos$gasto), by = 1))##
## studentized Breusch-Pagan test
##
## data: datos$gasto ~ seq(1, length(datos$gasto), by = 1)
## BP = 5.3713, df = 1, p-value = 0.02047Se rechaza la hipótesis nula y por tanto la serie no es homocedástica. Sin embargo, al considerar el logaritmo de los datos,
datos <- datos %>%
mutate(log_gasto = log(gasto))
bptest(datos$log_gasto ~ seq(1, length(datos$gasto), by = 1))##
## studentized Breusch-Pagan test
##
## data: datos$log_gasto ~ seq(1, length(datos$gasto), by = 1)
## BP = 1.965, df = 1, p-value = 0.161se observa que la varianza se estabiliza, así que ya se tiene un proceso homocedástico.
Ahora se considera la descomposición de la serie homocedástica.
descomposicion_tbl <- log(datos_ts) %>%
decompose() %>%
sw_tidy_decomp() %>%
rename(fecha = index) %>%
mutate(fecha = lubridate::yq(fecha)) %>%
pivot_longer(cols = -fecha, names_to = "tipo", values_to = "valor") %>%
mutate(tipo = as_factor(tipo), tipo = fct_recode(tipo, observado = "observed",
tendencia = "trend", ciclo = "season", ciclo_ajuste = "seasadj", aleatorio = "random"),
tipo = fct_relevel(tipo, c("observado", "tendencia", "ciclo", "ciclo_ajuste",
"aleatorio")))g_descompose <- descomposicion_tbl %>%
ggplot(aes(fecha, valor)) + geom_line(color = color_prin) + facet_grid(vars(tipo),
scales = "free_y") + labs(x = NULL, y = NULL, title = "Descomposición clásica")
g_descompose %>%
ggplotly() %>%
config(displayModeBar = FALSE)2 Imputación
Primero se elimina los valores deseados:
datos_na = datos_ts
datos_na[c(29, 38, 39)] = NA
datos_na## Qtr1 Qtr2 Qtr3 Qtr4
## 1964 10.00 11.85 11.70 13.42
## 1965 11.20 13.63 13.65 15.93
## 1966 13.33 16.05 15.92 18.22
## 1967 14.46 16.69 16.20 18.12
## 1968 15.10 16.85 16.79 19.03
## 1969 16.04 18.81 19.25 21.46
## 1970 17.47 20.33 20.26 21.66
## 1971 NA 20.60 20.14 22.79
## 1972 19.38 22.01 21.86 25.20
## 1973 21.50 NA NA 28.48
## 1974 24.10 28.16 28.23 31.92
## 1975 25.82 28.43 27.79 30.74
## 1976 25.87 29.70 30.41 34.52Y se visualiza los valores faltantes
ggplot_na_distribution(datos_na)
Ahora se utiliza distintos métodos para estimar los valores eliminados. El análisis será gráfico.
Se inicia con el método Kalman.
imp_kal <- na_kalman(datos_na)
kal <- ggplot_na_imputations(datos_na, imp_kal, title = "Método Kalman", xlab = "Tiempo",
ylab = "Gasto")
kal %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Ahora se incorpora la componente estacional en el método.
imp_kal_est <- na_seadec(datos_na, algorithm = "kalman", find_frequency = TRUE)
kal_est <- ggplot_na_imputations(datos_na, imp_kal_est, title = "Método Kalman estacional",
xlab = "Tiempo", ylab = "Gasto")
kal_est %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Ahora se realiza una imputación lineal mediante promedios móviles.
imp_ma_lin <- na_ma(datos_ts, k = 4, weighting = "linear")
ma_lin <- ggplot_na_imputations(datos_na, imp_ma_lin, title = "Promedios móviles lineales",
xlab = "Tiempo", ylab = "Gasto")
ma_lin %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)También se usa interpolación lineal
imp_int_lin <- na_interpolation(datos_ts, option = "linear")
int_lin <- ggplot_na_imputations(datos_na, imp_int_lin, title = "Interpolación lineal",
xlab = "Tiempo", ylab = "Gasto")
int_lin %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Otra opción es imputar con la media
imp_media <- na_mean(datos_ts, option = "mean")
media <- ggplot_na_imputations(datos_na, imp_media, title = "Media", xlab = "Tiempo",
ylab = "Gasto")
media %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Y también la mediana
imp_mediana <- na_mean(datos_ts, option = "median")
mediana <- ggplot_na_imputations(datos_na, imp_mediana, title = "Mediana", xlab = "Tiempo",
ylab = "Gasto")
media %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Comparando el método Kalman estacional fue el mejor y logró captar mejor el ciclo estacional.
3 Modelo
Para proponer distintos modelos se comienza analizando el ACF y el PACF. Antes de esto se determina si los datos transformados son estacionarios,
adf.test(log(datos_ts))$p.value## [1] 0.4818714kpss.test(log(datos_ts))$p.value## [1] 0.01como ambas pruebas coinciden en la no estacionariedad, se toma diferencias,
adf.test(diff(log(datos_ts)))$p.value## [1] 0.7243431kpss.test(diff(log(datos_ts)))$p.value## [1] 0.1aunque la segunda prueba acepta que el proceso es estacionario, la primera no lo hace, de nabera que se modifica la diferenciación para considerar el ciclo estacional \(s = 4\), obteniéndose
adf.test(diff(log(datos_ts), lag = 4))$p.value## [1] 0.01kpss.test(diff(log(datos_ts), lag = 4))$p.value## [1] 0.1y entonces ya se tiene datos estacionarios. En seguida se muestra sus correlogramas.
datos_tsb %>%
feasts::gg_tsdisplay(difference(log(gasto), lag = 4), plot_type = "partial") +
ggtitle("Diferencia estacional de logaritmo de la serie con correlogramas")
Por inspección del \(ACF\) se estima que \(q = 8\) y por el PACF se toma \(p =2\), así que se comienza con un proceso \(ARIMA(2,0,8)\times(0,1,0)_4\), que claramente tiene muchos parámetros, pero servirá como referencia.
En segundo lugar, se utiliza auto.arima para proponer los parámetros:
auto.arima(log(datos_ts)) %>%
summary(modelo_auto)## Series: log(datos_ts)
## ARIMA(2,0,1)(0,1,2)[4] with drift
##
## Coefficients:
## ar1 ar2 ma1 sma1 sma2 drift
## 1.7416 -0.8694 -0.7121 -0.3879 -0.3412 0.0188
## s.e. 0.1087 0.0934 0.1917 0.1865 0.1897 0.0008
##
## sigma^2 estimated as 0.0004258: log likelihood=118.39
## AIC=-222.79 AICc=-219.99 BIC=-209.69
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001984431 0.01854458 0.01392654 0.07957364 0.4692953 0.1714503
## ACF1
## Training set -0.05781879Como el orden de \(p\) es 2 y el de \(Q\) es 2, se ajusta cuatro modelos más, donde se varían estos parámetros. Se trabajará con los 6 modelos de manera simultánea.
ajuste <- datos_tsb %>%
model(arima_auto = ARIMA(log(gasto), stepwise = FALSE), inspeccion_208_010 = ARIMA(log(gasto) ~
0 + pdq(2, 0, 8) + PDQ(0, 1, 0), stepwise = FALSE), alt1_101_012 = ARIMA(log(gasto) ~
0 + pdq(1, 0, 1) + PDQ(0, 1, 2), stepwise = FALSE), alt2_201_011 = ARIMA(log(gasto) ~
0 + pdq(2, 0, 1) + PDQ(0, 1, 1), stepwise = FALSE), alt3_100_011 = ARIMA(log(gasto) ~
0 + pdq(1, 0, 0) + PDQ(0, 1, 1), stepwise = FALSE), alt4_101_110 = ARIMA(log(gasto) ~
0 + pdq(1, 0, 1) + PDQ(1, 1, 0), stepwise = FALSE))
nom <- c("arima_auto", "inspeccion_208_010", "alt1_101_012", "alt2_201_011", "alt3_100_011",
"alt4_101_110")g_modelos <- ajuste %>%
augment() %>%
clean_names() %>%
rename(modelo = model, ajuste = fitted) %>%
mutate(modelo = as_factor(modelo), modelo = fct_relevel(modelo, nom)) %>%
ggplot(aes(fecha, ajuste, colour = modelo)) + geom_line(linetype = 3, alpha = 0.5) +
scale_y_continuous(labels = scales::dollar) + scale_color_lancet() + labs(x = NULL,
y = "Gasto estimado", colour = "Modelo", title = "Comparación de modelos")
g_modelos %>%
ggplotly() %>%
config(displayModeBar = FALSE)3.1 Significancia coeficientes
Claramente utilizar intervalos de confianza es equivalente a realizar pruebas de hipótesis mediante el \(p-value\), sin embargo se realizará los dos métodos.
Para los intervalos de confianza se tiene que
coeficientes <- ajuste %>%
coef()
int_conf <- coeficientes %>%
mutate(inf = estimate - 1.96 * std.error, sup = estimate + 1.96 * std.error,
no_sig = inf <= 0 & 0 <= sup, significante = !no_sig) %>%
select(.model:estimate, significante) %>%
rename(modelo = .model, parámetro = term, estimación = estimate)
int_confy el total de parámetros no significativos para cada modelo es:
total_no_sig <- int_conf %>%
group_by(modelo) %>%
summarise(total_no_significante = sum(!significante)) %>%
arrange(desc(total_no_significante))
total_no_sigLo mismo ocurre mediante el \(p-value\):
sig_p_value <- coeficientes %>%
mutate(significante = p.value < 0.05) %>%
relocate(significante, .after = estimate) %>%
select(.model:significante) %>%
rename(modelo = .model, parámetro = term, estimación = estimate)
sig_p_valuetotal_no_sig_p_value <- sig_p_value %>%
group_by(modelo) %>%
summarise(total_no_significante = sum(!significante)) %>%
arrange(desc(total_no_significante))
total_no_sig_p_valueSe puede concluir que todos los modelos a excepción de alt3_100_011, incluyen al menos un parámetro no significante, lo cual introduce ruido al modelo (incluido arima_auto). Sin embargo, habrá que determinar con el siguiente anáisis si alt3_100_011 tiene un buen desempeño.
3.2 Bondad de ajuste
Primero se analiza la verosimilitud,
info <- ajuste %>%
glance() %>%
clean_names() %>%
rename(modelo = model) %>%
select(modelo:bic, -sigma2)
info %>%
rename(log_verosimilitud = log_lik) %>%
select(modelo, log_verosimilitud) %>%
arrange(desc(log_verosimilitud))y claramente el modelo con auto.arima es el que maximiza la verosimilitud.
Ahora se compara los criterios de información
info %>%
select(-log_lik) %>%
arrange(aic)y se contrasta los datos de manera gráfica;
g_info <- info %>%
select(-c(log_lik, ai_cc)) %>%
pivot_longer(-modelo, names_to = "metrica", values_to = "valor") %>%
ggplot(aes(modelo, valor, colour = modelo)) + geom_segment(aes(xend = modelo,
yend = 0), size = 1, lineend = "round") + scale_color_lancet() + facet_wrap(vars(metrica),
scales = "free_y") + theme(legend.position = "None") + labs(title = "Criterios de información",
x = "", y = "", colour = "Modelo") + theme(legend.position = "None", axis.text.x = element_text(angle = -30))
g_info %>%
ggplotly() %>%
config(displayModeBar = FALSE)y resulta evidente que el modelo con auto.arima minimiza todos los criterios, por lo que hasta ahora es considerablemente mejor.
3.3 Análisis de residuos
Primero se compara todos los residuos,
g_residuos <- ajuste %>%
augment() %>%
ggplot(aes(fecha, .resid, colour = .model)) + geom_line() + labs(title = "Comparación de residuos",
x = NULL, y = NULL, colour = "Modelo")
g_residuos %>%
ggplotly() %>%
config(displayModeBar = FALSE)corr <- nom %>%
map(~select(ajuste, .x)) %>%
set_names(nom) %>%
map(augment) %>%
map_df(`[`, ".innov", .id = "modelo") %>%
pivot_wider(names_from = modelo, values_from = .innov) %>%
unnest() %>%
mutate(fecha = datos$fecha) %>%
as_tsibble() %>%
pivot_longer(-fecha, names_to = "modelo", values_to = "residuo") %>%
mutate(modelo = as_factor(modelo)) %>%
group_by(modelo)## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(.x)` instead of `.x` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.## Using `fecha` as index variable.g_residuos_facet <- corr %>%
ggplot(aes(fecha, residuo, colour = modelo)) + geom_line() + scale_y_continuous(labels = scales::percent) +
scale_color_lancet() + facet_wrap(vars(modelo), scales = "free_y") + theme(legend.position = "None") +
labs(title = "Residuos", x = "", y = "") + theme(legend.position = "None")
g_residuos_facet %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)corr %>%
ACF(residuo) %>%
autoplot()3.4 Verifcación de supuestos en los residuos
residuos <- ajuste %>%
augment() %>%
as_tibble() %>%
clean_names() %>%
rename(modelo = model, residuo = innov) %>%
mutate(modelo = as_factor(modelo))3.4.1 Media cero
Se utiliza las prueba \(t\) y Wilcoxon, cuyas hipótesis nulas son que la media es 0.
pruba_t <- residuos %>%
group_by(modelo) %>%
summarise(t = t.test(residuo, mu = 0)$p.value, wilocx = wilcox.test(residuo,
mu = 0)$p.value) %>%
pivot_longer(-modelo, names_to = "prueba", values_to = "p_value")
g_t <- pruba_t %>%
ggplot(aes(modelo, p_value, colour = modelo)) + geom_segment(aes(xend = modelo,
yend = 0), size = 1, lineend = "round") + scale_y_continuous(labels = scales::percent) +
geom_hline(yintercept = 0.05, colour = color_prin) + scale_y_continuous(labels = scales::percent) +
facet_wrap(vars(prueba), scales = "free_y") + scale_color_lancet() + labs(title = "Pruebas para media 0",
x = "Modelo", y = "p-value") + theme(legend.position = "None", axis.text = element_text(angle = -30))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.g_t %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Aquí todos los modelos pasan la prueba en favor de la hipótesis nula, a excepción del modelo inspeccion_208_010, el cual es el que tiene más parámetros.
3.4.2 Varianza constante
Se usa la prueba de Breusch-Pagan cuya hipótesis nula es la homocedasticidad.
aux <- 1:length(datos_tsb$fecha)
breusch_pagan <- residuos %>%
group_by(modelo) %>%
summarise(p_value = bptest(residuo ~ aux)$p.value)
g_bp <- breusch_pagan %>%
ggplot(aes(modelo, p_value, colour = modelo)) + geom_segment(aes(xend = modelo,
yend = 0), size = 1, lineend = "round") + geom_hline(yintercept = 0.05, colour = color_prin) +
scale_y_continuous(labels = scales::percent) + scale_color_lancet() + labs(title = "Prueba de homocedasticidad de Breusch-Pagan",
x = "Modelo", y = "p-value") + theme(legend.position = "None")
g_bp %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)en todos los casos el \(p-value\) es considerablemente mayor al \(5 \%\), así que se acepta la hipótesis de homocedasticidad.
3.4.3 No correlación
Se aplica la prueba ljung_box cuya hipótesis nula es la no correlación con un retraso igual al doble del ciclo estacioona, esto es, 8.
lbox <- ajuste %>%
augment() %>%
features(.innov, ljung_box, lag = 8) %>%
rename(modelo = .model, p_value = lb_pvalue)g_lbox <- lbox %>%
ggplot(aes(modelo, p_value, colour = modelo)) + geom_segment(aes(xend = modelo,
yend = 0), size = 1, lineend = "round") + geom_hline(yintercept = 0.05, colour = color_prin) +
scale_y_continuous(labels = scales::percent) + scale_color_lancet() + labs(title = "Prueba de no correlación",
x = "Modelo", y = "p-value") + theme(legend.position = "None")
g_lbox %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Se nota que todos los modelos salvo alt3_100_011 y alt4_101_110 pasaron la prueba y por tanto son no correlacionados.
total_parametros <- ajuste %>%
coef() %>%
rename(modelo = .model) %>%
group_by(modelo) %>%
summarise(num_parametros = n()) %>%
arrange(num_parametros)3.4.4 Normalidad
3.4.4.1 Gráfica
g_qq <- residuos %>%
ggplot(aes(sample = residuo, colour = modelo)) + geom_qq() + stat_qq_line() +
scale_color_lancet() + scale_y_continuous(labels = scales::percent) + facet_wrap(vars(modelo),
scales = "free_y") + labs(title = "Gráfico qq", x = NULL, y = NULL) + theme(legend.position = "None")
g_qq %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)3.4.4.2 Pruebas normalidad
pruebas_nor <- residuos %>%
group_by(modelo) %>%
summarise(jarque = jarque.bera.test(residuo)$p.value, shapiro = shapiro.test(residuo)$p.value,
anderson = ad.test(residuo)$p.value) %>%
pivot_longer(-modelo, names_to = "prueba", values_to = "p_value") %>%
mutate(prueba = as_factor(prueba), prueba = fct_relevel(prueba, c("jarque", "shapiro",
"anderson")))
g_pruebas_nor <- pruebas_nor %>%
ggplot(aes(prueba, p_value, colour = modelo)) + geom_segment(aes(xend = prueba,
yend = 0), size = 1, lineend = "round") + geom_hline(yintercept = 0.05, colour = color_prin) +
scale_y_continuous(labels = scales::percent) + scale_color_lancet() + facet_wrap(vars(modelo),
scales = "free_y") + theme(legend.position = "None") + labs(title = "Pruebas de normalidad",
x = "Prueba", y = "p-value", colour = "Modelo") + theme(legend.position = "None")
g_pruebas_nor %>%
ggplotly() %>%
layout(showlegend = FALSE) %>%
config(displayModeBar = FALSE)Resalta que el modelo arima_auto no pasa ninguna prueba de normalidad, mientras que todos los demás pasan al menos una (e incluso las 3), de modo que se acepta la hipótesis de normalidad para todos ellos.
3.5 Elección de modelo
Finalmente, se elige el modelo arima_auto ya que fue el mejor respecto a la verosimilitud u los criterios de información AIC y BIC, no tenía tantos parámetros a estimar y solo falló en la normalidad de los residuos, aunque lo cierto es que tiene un parámetro no significante que distorsiona la precisión de las predicciones.
4 Predicción
Primero se muestra la predicción obtenida por todos los modelos:
ajuste %>%
forecast(h = "2 years") %>%
autoplot(datos_tsb) + labs(title = "Predicciones a dos años", x = "Fecha", y = "Gasto")
Ahora solo se muestra las predicciones para el modelo elegido, \(ARIMA(2,0,8)\times(0,1,0)_4\):
pred <- ajuste %>%
select(arima_auto) %>%
forecast(h = "2 years")
pred %>%
rename(modelo = .model, prediccion = .mean, distribucion = gasto)pred %>%
autoplot(datos_tsb, color = color_sec) + labs(title = "Predicción a 2 años",
x = "Fecha", y = "Gasto", level = "Intervalo de confianza") + theme(legend.position = "bottom")
Finalmente se aprecia que el modelo logra replicar el comportamiento estacional de los dos ciclos de predicción. En la particular el último ciclo observado presentó un comportamiento atípico, el cual fue bien replicado en la predicción en el segundo ciclo.