SeriesTCode
Vásquez Guerra Carlos Fernando
12/2020
Datos
Así como se tienen objetos survival para un análisis de supervivencia, en este caso tendremos objetos referentes a series de tiempo en R. Hay varias formas de crear una serie de tiempo en R, pero la más sencilla y que basta para nuestros fines son los objetos ts. Para ver otros tipos véase el siguiente enlace. Antes de continuar se dejan aquí unas Cheat Sheets que pueden ser de utilidad cada vez que se realice un análisis de series de tiempo en R.
Este último enlace hace referencia al paquete lubridate, un paquete especializado en el manejo de fechas; para ver un poco de esta librería así como el uso y manipulación de fechas con objetos básicos de R véase el post Manejo de fechas con R. Véase lo sencillo que es crear una serie de tiempo con la función base::ts().
set.seed(123)
datos <- cumsum(rnorm(100))
(firstSerie <- ts(data = datos, start = 2020, frequency = 1))Time Series:
Start = 2020
End = 2119
Frequency = 1
[1] -0.5604756 -0.7906531 0.7680552 0.8385636 0.9678513 2.6829163
[7] 3.1438325 1.8787713 1.1919184 0.7462564 1.9703382 2.3301521
[13] 2.7309235 2.8416062 2.2857651 4.0726782 4.5705287 2.6039116
[19] 3.3052675 2.8324760 1.7646523 1.5466774 0.5206730 -0.2082182
[25] -0.8332575 -2.5199508 -1.6821638 -1.5287907 -2.6669276 -1.4131127
[31] -0.9866485 -1.2817199 -0.3865943 0.4915392 1.3131203 2.0017605
[37] 2.5556782 2.4937665 2.1878038 1.8073328 1.1126258 0.9047086
[43] -0.3606878 1.8082682 3.0162302 1.8931216 1.4902368 1.0235814
[49] 1.8035465 1.7201775 1.9734960 1.9449492 1.9020788 3.2706810
[55] 3.0449101 4.5613807 3.0126279 3.5972416 3.7210958 3.9370374
[61] 4.3166769 3.8143534 3.4811461 2.4625707 1.3907795 1.6943081
[67] 2.1425179 2.1955221 3.1177896 5.1678743 4.6768431 2.3676742
[73] 3.3734127 2.6642120 1.9762034 3.0017747 2.7170017 1.4962840
[79] 1.6775875 1.5386961 1.5444603 1.9297407 1.5590807 2.2034572
[85] 1.9829707 2.3147526 3.4115916 3.8467731 3.5208415 4.6696492
[91] 5.6631530 6.2115500 6.4502817 5.8223756 7.1830281 6.5827685
[97] 8.7701015 10.3027121 10.0670118 9.0405909¿Qué tipo de serie o proceso es descrito en el siguiente ejemplo? Tal vez ayude un poco al ver la gráfica correspondiente (La función plot es genérica y tenemos la correspondiente plot.ts() por lo que se pudo haber obtenido la gráfica anterior con plot(firstSerie)). A lo largo de este documento se mostraran diversas formas de graficación, pero siempre se hará mención de las gráficas que se pueden obtener con la función plot().
(tibble(time = time(firstSerie), serie = firstSerie) %>%
ggplot(aes(x = time, y = firstSerie)) +
geom_line(color = "dodgerblue3") +
labs(x = "Tiempo", y = NULL) +
general_theme) %>%
ggplotly()En el constructor de una serie de tiempo se tienen los siguientes argumentos
data: Los datos.start: Punto de inicio de la serie de tiempo, este puede ser un vector indicando el año y día correspondiente a dicho año.end: Análogo astartsalvo que este puede ser calculado por los datos. Nos sirve para acotar la serie si así lo deseamos.frecuency: Frecuencia con la que los datos son tomados en cuenta en el año
Hay que hacer algunas aclaraciones sobre este último parámetro. La frecuencia nos ayudará determinar cuantas observaciones se consideran por año, por lo que al indicar frecuency = 1, se está indicando que cada observación es un año. En la siguiente tabla se resumen algunas temporalidades comunes
| Temporalidad | Frecuencia |
|---|---|
| Anual | 1 |
| Semestral | 2 |
| Cuatrimestral | 3 |
| Trimestral | 4 |
| Bimestral | 6 |
| Mensual | 12 |
| Semanal | 52 |
| Diaria | 365 |
Además se pueden colocar otro valor que represente las horas, minutos, etc. En general este debe ser un número entero pero no hay problema en que no lo sea, así bien se podría colocar un valor más exacto para un periodo semanal: \(365.25/7=52.18\). Este parámetro es de suma importancia ya que tenga importantes repercusiones en la estacionalidad de la serie, por lo que siempre hay que pensar bien en que temporalidad se dan los datos.
Descomposición
En general podemos pensar que una serie de tiempo se puede descomponer en tres partes importantes: Tendencia, Estacionalidad y Remanente (o componente aleatoria) y de acuerdo al enfoque que se desee aplicar, se pueden utilizar dos tipos de modelos clásicos de descomposición: Aditivo y Multiplicativo. Para ver más acerca de este último consúltese el siguiente enlace.
Recordemos que el modelo aditivo supone que la serie de tiempo tiene la siguiente estructura: \(X_t = T_t+E_t+It\). En este caso, la forma clásica de descomposición consiste en obtener la tendencia mediante un modelo \(MA\), del cual se hablará más adelante, por lo que por el momento no se realizará esta descomposición con todos los pasos. En cambio, se utilizarán algunas funciones de R que nos dan una descomposición por distintos métodos. Antes de eso véamos una transformación que nos ayudará a obtener la tendencia de una serie de tiempo
Medias móviles
Recordemos que una media móvil de orden \(q\) es descrita por la siguiente ecuación
\[
\hat T_t=(2q+1)^{-1}\sum_{j=-q}^q X_{t+j} = \frac{1}{m}\sum_{j = -q}^qX_{t+j}
\] con \(\space q+1\leq t\leq n-q\) y \(m = 1/(2q+1)\). Es decir, sólo es un simple promedio de los valores de la serie en \(q\) periodos de tiempo. Al hacer esto, estamos pensando que las observaciones cercanas serán parecidas por lo que se elimina esa componente aleatoria dejando sólo el comportamiento de la serie referente a la tendencia. Podríamos hacer una función que calculará esto, pero se deja como ejercicio al lector, en este caso vamos a utilizar la función forecast::ma() con el siguiente ejemplo. Este contiene la información correspondiente al número mensual de pasajeros en distintas aerolíneas entre los años 1949-1960.
AirPassengers Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432Vamos a aplicar aplicar periodos de tres observaciones
(AirPassnger_MA3 <- forecast::ma(AirPassengers, 3)) Jan Feb Mar Apr May Jun Jul Aug
1949 NA 120.6667 126.3333 127.3333 128.3333 134.6667 143.6667 144.0000
1950 119.6667 127.3333 134.0000 133.6667 136.3333 148.0000 163.0000 166.0000
1951 145.0000 157.6667 163.6667 171.0000 171.0000 183.0000 192.0000 194.0000
1952 172.3333 181.3333 184.6667 185.6667 194.0000 210.3333 230.0000 227.0000
1953 195.3333 209.3333 222.3333 233.3333 235.6667 245.3333 259.6667 257.6667
1954 197.6667 209.0000 216.6667 232.0000 241.6667 266.6667 286.3333 284.6667
1955 234.6667 247.3333 256.3333 268.6667 284.6667 316.3333 342.0000 341.0000
1956 279.6667 292.6667 302.3333 316.0000 335.0000 368.3333 397.3333 391.0000
1957 307.3333 324.0000 335.0000 353.0000 375.0000 414.0000 451.3333 445.3333
1958 331.3333 340.0000 342.6667 357.6667 382.0000 429.6667 477.0000 466.6667
1959 346.3333 369.3333 381.3333 407.3333 429.3333 480.0000 526.3333 523.3333
1960 404.3333 409.0000 423.6667 450.6667 489.3333 543.0000 587.6667 578.6667
Sep Oct Nov Dec
1949 134.3333 119.6667 113.6667 112.3333
1950 153.6667 135.0000 129.0000 133.0000
1951 181.6667 164.0000 158.0000 161.0000
1952 214.0000 190.6667 185.6667 187.3333
1953 240.0000 209.3333 197.3333 195.0000
1954 260.3333 230.3333 220.3333 224.6667
1955 311.0000 274.3333 263.0000 266.3333
1956 355.3333 310.6667 294.3333 297.3333
1957 406.0000 352.0000 329.3333 327.0000
1958 422.6667 357.6667 335.3333 335.6667
1959 476.3333 410.6667 391.3333 394.6667
1960 525.0000 453.0000 427.6667 NABien, ahora véamos como se ve graficamente esto
forecast::autoplot(AirPassengers, color = "black", alpha = 0.5) +
forecast::autolayer(AirPassnger_MA3, color = "red") +
general_theme +
labs(x = "Años", y = "Número de pasajeros")
Y ahora para diferentes periodos
colors <- c(futurevisions::futurevisions("mars")[1:4], futurevisions::futurevisions("europa")[2])
AirPassengers_MAPlots <- autoplot(AirPassengers, series = "Datos", alpha = 0.5) +
autolayer(ma(AirPassengers, 1), series = "MA1") +
autolayer(ma(AirPassengers, 2), series = "MA2") +
autolayer(ma(AirPassengers, 3), series = "MA3") +
autolayer(ma(AirPassengers, 4), series = "MA4") +
autolayer(ma(AirPassengers, 5), series = "MA5") +
general_theme +
scale_color_manual(values = c("black",colors)) +
labs(x = "Años", y = "Número de pasajeros") +
guides(colour=guide_legend(title = "Modelo")) +
theme(legend.position = "right",
legend.direction = "vertical",
legend.box = "vertical")
ggplotly(AirPassengers_MAPlots)¿Que sucede al aumentar la cantidad de periodos?
STL & More
Existen muchas funciones en R que nos ayudan a descomponer una serie de tiempo, entre ellas dos de las más utilizadas son decompose() y stl(). La primera realiza una descomposición de la serie utilizando medias móviles y stl() es una acrónimo para Seasonal and Trend decomposition using Loess, el cual es un método propuesto en (???). Aquí se deja un enlace para más información sobre las ventajas de este último método. Estos métodos nos devolverán los valores para cada componente a través del tiempo, pero en este caso vamos a utilizar la información que nos proporcionan y realicemos directamente algunas gráficas (se puede usar directamente la función plot en cada uno de las descomposiciones).
AirPassengers_descomposition <- stats::decompose(AirPassengers, type = "additive")
(tibble(time = time(AirPassengers_descomposition$x),
Datos = AirPassengers_descomposition$x,
Tendencia = AirPassengers_descomposition$trend,
Estacionalidad = AirPassengers_descomposition$seasonal,
Remanente = AirPassengers_descomposition$random) %>%
gather("type", "value", 2:5) %>%
mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>%
ggplot(aes(x = time, y = value)) +
general_theme +
geom_line() +
facet_grid(type~., scales = "free_y") +
labs(x = "Años", y = "Número de pasajeros") +
ggtitle("Descomposición temporal de AirPassengers")) %>% ggplotly()Ahora veamos la gráfica que resultad de la descomposición por el método stl().
highcharter::hchart(stl(AirPassengers, s.window = "per"))Como se puede apreciar, los resultados son bastantes similares, por lo que nuestras conclusiones serán equivalentes. En general existe una homogeneidad entre los métodos de descomposición. Para ver otros véase Time series decomposition.
- ¿Para que nos sirve esto?
- ¿Que componente eliminarías para tener un mejor entendimiento del comportamiento de la serie?
Un caso curioso
El ejemplo tratado anteriormente es un ejemplo tratado y “bonito”, pero que pasa cuando tenemos que lidiar con datos no tratados. Por ejemplo, vamos a tomar los datos los datos del coronavirus de Our World in Data. Este estudio es válido hasta el día (2020-12-09).
data_coronavirus <- read_csv("https://covid.ourworldindata.org/data/owid-covid-data.csv", col_names = TRUE)
(data_coronavirus_EC <- data_coronavirus %>% filter(date <= as.Date("2020-12-09")))Y vamos a considerar específicamente los datos de nuevos contagios en México.
(coronavirus_Mexico <- data_coronavirus_EC %>%
dplyr::filter(location == "Mexico") %>%
dplyr::select(date, new_cases) %>%
dplyr::filter(!is.na(new_cases)) %>%
dplyr::filter(new_cases != 0))Por lo que ahora vamos a crear nuestra serie de tiempo con estos datos. Véase el como estamos instanciando el parámetro start y frequency, este último por que conocemos que esa es la frecuencia con la que se dan a conocer tales datos.
tseries <- ts(coronavirus_Mexico$new_cases,
start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%j"))),
frequency = 365)Pero, ¿Qué sucede al ejecutar el comando stl(tseries)? Bueno, se obtendría el siguiente error:
Error in stl(tseries) : series is not periodic or has less than two periods
Este error se debe a que no tenemos ni un ciclo (son necesarios al menos dos para obtener una estimación sobre la estacionalidad) ya que contamos con 279 días de un año. Para solucionar este inconveniente no utilizaremos la descomposión proporcionada por stl() o decompose(). Utilizaremos un estimación de la estacionalidad basada en series de Fourier, tal como se muestra en el post Seasonal decomposition of short time series. Donde se utiliza el hecho de que es posible utilizar un modelo de regresión lineal para descomponer una serie de tiempo, así como se explica en el siguiente capítulo Time series regression models del bookdown ya varias veces mencionado.
Para estos en particular, consideraremos que los datos tienen una estacionalidad semanal, por lo que al tener 52 semanas, el parámetro \(K\) mostrado en dicho modelo debe cumplir que se \(K<=52/2 = 26\), tal como lo establece la frecuencia Nyquist. Para seleccionar el mejor \(K\), es decir la cantidad de pares de funciones \(sen\left(\frac{2\pi tk}{m}\right)\) y \(cos\left(\frac{2\pi tk}{m}\right)\), utilizaremos el modelo que obtenga el mínimo valor de la función CV() (CV statistic).
#Función para obtener el valor número de componentes de Fourier de acuerdo al Coeficiente de Variación
k_min <- function(ts, kmax){
season_f <- map(1:kmax, ~ fourier(ts, .x)) #Obtenemos las apxomiaciones de estacionalidad por coef. de Fourier.
lineal_ajusts <- map(season_f, ~tslm(ts ~ trend + .x)) #Creamos los ajuste lineales
coef_fourier <- map(lineal_ajusts, ~CV(.x)) #Aplicamos la función CV a los anteriores modelos
coef_fourier <- map_dbl(coef_fourier, ~'[['(.x, 1)) #Obtenemos los resultados interesantes
which(coef_fourier == min(coef_fourier)) #Determinamos cual es el k que minimiza el CV
}
(min_casesT <- k_min(tseries, 26))[1] 4Entonces, vamos a realizar el ajuste lineal mediante la función forecast::tslm() y la función forecast::fourier() la cual nos devolverá los valores evaluados en las componentes de Fourier.
fitted_ts <- tslm(tseries ~ trend + fourier(tseries, min_casesT))
fitted_ts
Call:
tslm(formula = tseries ~ trend + fourier(tseries, min_casesT))
Coefficients:
(Intercept) trend
-11860.2 125.7
fourier(tseries, min_casesT)S1-365 fourier(tseries, min_casesT)C1-365
6532.1 8590.5
fourier(tseries, min_casesT)S2-365 fourier(tseries, min_casesT)C2-365
-2957.8 3635.5
fourier(tseries, min_casesT)S3-365 fourier(tseries, min_casesT)C3-365
-1865.4 394.3
fourier(tseries, min_casesT)S4-365 fourier(tseries, min_casesT)C4-365
-321.4 -740.5 Entonces, como la tendencia la estamos considerando como un modelo lineal, tomamos los valores correspondientes del anterior modelo para crear la tendencia de nuestra serie.
trend <- coef(fitted_ts)[1] + coef(fitted_ts)['trend']*seq_along(coronavirus_Mexico$new_cases)
trend %>% head()[1] -11734.52 -11608.83 -11483.13 -11357.44 -11231.74 -11106.05Ahora sí, vamos a crear nuevamente el gráfico para ver nuestros resultados. Véase como se calculo la estacionalidad.
(tibble(time = coronavirus_Mexico$date,
Datos = (coronavirus_Mexico$new_cases),
Tendencia = trend,
Estacionalidad = (coronavirus_Mexico$new_cases) - trend - residuals(fitted_ts),
Remanente = residuals(fitted_ts)) %>%
gather("type", "value", 2:5) %>%
mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>%
ggplot(aes(x = time, y = value)) +
geom_line() +
general_theme +
facet_grid(type~., scales = "free_y") +
labs(x = "Tiempo", y = "Población") +
ggtitle("Descomposición temporal. Casos nuevos diarios para México")) %>% ggplotly()¿Cuál fue nuestro error aquí?
Bueno, el problema aquí es el gran supuesto que estamos haciendo sobre el comportamiento lineal en la tendencia. Vamos a considerar entonces un comportamiento no lineal en esta. Para más información véase el siguiente capítulo. En este caso vamos a utilizar una función piecewise para la tendencia considerando como puntos importantes (nodos) las siguientes fechas: 2020-5-5, 2020-8-1, 2020-10-4 y 2020-11-17.
t <- time(tseries)
#Puntos de corte
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-5-5")
t.break1 <- time(tseries)[61]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-8-1")
t.break2 <- time(tseries)[149]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-10-4")
t.break3 <- time(tseries)[213]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-11-17")
t.break4 <- time(tseries)[257]Y creamos la función por partes como objetos ts().
tb1 <- ts(pmax(0, t - t.break1), start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%j"))), frequency = 365)
tb2 <- ts(pmax(0, t - t.break2), start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%j"))), frequency = 365)
tb3 <- ts(pmax(0, t - t.break3), start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%j"))), frequency = 365)
tb4 <- ts(pmax(0, t - t.break4), start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%j"))), frequency = 365)Ajustamos el modelo no lineal
fit.pw <- tslm(tseries ~ t + tb1 + tb2 + tb3 + tb4)
fit.pw <- fitted(fit.pw)
fit.pw %>% head(10)Time Series:
Start = c(2020, 59)
End = c(2020, 68)
Frequency = 365
[1] -587.4754 -548.4677 -509.4601 -470.4524 -431.4447 -392.4370 -353.4294
[8] -314.4217 -275.4140 -236.4064Y finalmente graficamos creando un nuevo modelo contemplando los coeficientes de Fourier para la estacionalidad considerando nuevamente que esta sea semanal.
fitted_ts <- tslm(tseries ~ fit.pw + fourier(tseries, k_min(tseries, 26)))
(tibble(time = coronavirus_Mexico$date,
Datos = (coronavirus_Mexico$new_cases),
Tendencia = as.vector(fit.pw),
Estacionalidad = (coronavirus_Mexico$new_cases) - as.vector(fit.pw) - residuals(fitted_ts),
Remanente = residuals(fitted_ts)) %>%
gather("type", "value", 2:5) %>%
mutate(type = factor(type, levels = c("Datos", "Tendencia", "Estacionalidad", "Remanente"))) %>%
ggplot(aes(x = time, y = value)) +
geom_line() +
general_theme +
facet_grid(type~., scales = "free_y") +
labs(x = "Tiempo", y = "Población") +
ggtitle("Descomposición temporal. Casos nuevos diarios para México")) %>% ggplotly()Véase que en este caso mejora un poco, ya que la tendencia ahora se ajusta mejor a los datos proporcionados y la estacionalidad marca un patrón, aunque es evidente la presencia de las funciones sinusoidales en la estacionalidad. Los remanentes siguen siendo parecidos al comportamiento de la gráfica aunque estos redujeron su varianza, lo cual es un punto a favor.
Algunas el comportamiento de los datos no estarán acorde a nuestras expectativas y tendremos que eliminar los supuestos o lo que esperamos ver. En este caso en particular sería interesante considerar otro tipo de estimación, ya sea sobre la tendencia, estacionalidad o ambos. Lo que se quiere dejar con esta sección son posibles problemas que el lector puede encontrar al hacer un análisis de series de tiempo, haciendo que este busque otro tipo de técnicas para dara solución a un problema dado.
Autocorrelación
Para obtener los valores y las gráficas (con hacer uso de la función plot() basta para ver el comportamiento) de autocorrelación y autocorrelación parcial se utilizan las funciones acf() y pacf() de manera correspondiente.
(autoplot(acf(AirPassengers,plot = FALSE), colour = "darkorchid4") +
geom_hline(yintercept = 0) +
ggtitle("Autocorrelación Simple") +
general_theme) +
(autoplot(pacf(AirPassengers,plot = FALSE), colour = "darkorchid4") +
geom_hline(yintercept = 0) +
ggtitle("Autocorrelación Parcial") +
labs(y = "Partial ACF")) +
plot_annotation(title = 'AirPassengers',theme = theme(plot.title = element_text(hjust = 0.5))) +
general_theme
(forecast::ggAcf(tseries) + ggtitle("Autocorrelación Simple") + general_theme) +
(ggPacf(tseries) + ggtitle("Autocorrelación Parcial")) +
plot_annotation(title = 'Casos diarios nuevos de Covid en México',
theme = theme(plot.title = element_text(hjust = 0.5))) +
general_theme
Hay que hacer unas cuantas aclaraciones:
- Las lineas azules indican si las correlaciones son significativamente diferentes de cero. En nuestros ejemplos se tiene una gran cantidad de correlaciones estadísticamente significativas, por lo que los valores de la serie de tiempo están correlacionados entre distintos desfaces.
- Si sólo se tuvieran valores entre las lineas azules, indicaría que las observaciones no tienen alguna autocorrelación, tal como sucede en el ruido blanco.
set.seed(11)
whitenoise <- tibble(x = 1:100, y = rnorm(100, mean = 0, sd = 1))
(whitenoise %>%
ggplot(aes(x = x, y = y)) +
geom_line(color = "darkslateblue") +
lims(y = c(-4, 4)) +
labs(x = "Tiempo") +
ggtitle(TeX("N(0,1)")) +
geom_hline(yintercept = 0) +
general_theme) +
(autoplot(acf(ts(whitenoise$y, start = 1, end = 100), plot = FALSE), colour = "darkslateblue") +
geom_hline(yintercept = 0) +
ggtitle("Autocorrelación") +
labs(x = "Rezago") +
general_theme) 
- En algunas ocasiones se pueden ver patrones en este tipo de gráficas indicando una fuerte estacionalidad o tendencia.
- Cuando existe una tendencia, las autocorrelaciones tenderán a ser grandes al inicio e irán decayendo a medida que aumentan los rezagos.
- Cuando existe estacionalidad, las autocorrelaciones serán más grandes para valores estacionarios (múltiplos) que en otros rezagos.
- Estos comportamientos se pueden combinar como en nuestro primer ejemplo: AirPassengers
APgraph <- autoplot(AirPassengers) + ggtitle("AirPassengers") +
labs(x = "Tiempo", y = "Número de pasajeros")+
general_theme
APACF <- ggAcf(AirPassengers) +
ggtitle("Autocorrelación Simple: AirPassengers")+
general_theme
subplot(ggplotly(APgraph), ggplotly(APACF), nrows = 2)La clara disminución de los coeficientes es debido a la tendencia que tiene la serie (los datos van aumentando) y aquellos valores como cada 12 rezagos es debido a la estacionalidad (cada 12 años hay se alcanzan valores más grandes).
Veamos otro ejemplo sencillo, una función sinusoidal.
tibble_sen <- tibble(x = 1:100, y = sin(x * pi/3))
sin_plot1 <- tibble_sen %>% ggplot(aes(x = x, y = y)) + geom_line() + ggtitle("Autocorrelación de una función sinusoidal") + general_theme
sin_plot2 <- ggAcf(tibble_sen$y, lag.max = 30) + ggtitle(NULL)+ general_theme
subplot(ggplotly(sin_plot1), ggplotly(sin_plot2), nrows = 2)Es notoria la falta de tendencia en la propia serie de tiempo y en la gráfica de autocorrelación se puede distinguir fácilmente por que los coeficientes significantes mantienen una cierta altura en todo el correlograma. Otro punto importante es que hay una clara estacionalidad (\(\sim3\) periodos) y esto fue por la construcción de la función. En el siguiente link se pueden modificar algunos parámetros de una función sinusoidal.
Otro punto interesante sobre la autocorrelación es que la suma de funciones periódicas es la suma de las autocorrelaciones individuales, para ver esto, vamos a tomar otra función armónica
tibble_sen_2 <- tibble(x = 1:100, y = sin(x * pi/10))
sin_plot3 <- tibble_sen_2 %>% ggplot(aes(x = x, y = y)) + geom_line() + ggtitle("Autocorrelación de una función sinusoidal") + general_theme
sin_plot4 <- ggAcf(tibble_sen_2$y, lag.max = 30) + ggtitle(NULL)+ general_theme
subplot(ggplotly(sin_plot3), ggplotly(sin_plot4), nrows = 2)Entonces, tendríamos el siguiente comportamiento para la función \(sen\left(x\frac{\pi}{3}\right)+sen\left(x\frac{\pi}{10}\right)\).
tibble_sen_sum <- tibble(x = 1:100, y = sin(x * pi/3) + sin(x * pi/10))
sin_sum_plot <- tibble_sen_sum %>% ggplot(aes(x = x, y = y)) + geom_line() + ggtitle("Autocorrelación de una función sinusoidal") + general_theme
sin_sumACF_plot <- ggAcf(tibble_sen_sum$y, lag.max = 30) + ggtitle(NULL)+ general_theme
subplot(ggplotly(sin_sum_plot), ggplotly(sin_sumACF_plot), nrows = 2)Pero, ¿Qué desventaja tiene el uso del ACF? ¿Qué ventaja tiene o para que utilizamos la PACF? Veamos la siguiente gráfica
(tibble(x = 1:100) %>% ggplot(aes(x = x, y = x)) +
geom_point() +
ggtitle("Datos no estacionarios") +
labs(x = "Tiempo", y = TeX("$f(x)$")) +
general_theme) /
((tibble(x = 1:100) %>% ggAcf()+
ggtitle("Autocorrelación simple") +
labs(x = "Lag", y = "ACF") +
general_theme) +
(tibble(x = 1:100) %>% ggPacf()+
ggtitle("Autocorrelación parcial") +
scale_y_continuous(name = "PACF", position = 'right') +
labs(x = "Lag") +
general_theme))
La ausencia de estacionariedad genera que la gráfica de autocorrelación simple sea engañosa y no nos de información útil, ya que es evidente la tendencia pero nos indica que todos los rezagos están correlacionados, lo cual es cierto, pero ¿esto que significa? Algo más interesante es lo que nos dice la gráfica de autocorrelación parcial, la cual nos esta indicando que, si no consideramos las correlaciones intermedias entre dos rezagos, sólo el valor inmediato anterior esta correlacionado con el estudiado.
Ahora, ya que se entendió un poco mejor el comportamiento de este tipo de gráficas, ¿Para que nos sirve determinar la correlación entre diferentes rezagos?
Identificación y modelación lineal
Recordando lo que indica el capítulo 7 del libro.
Los modelos de series de tiempo \(AR\), \(MA\) y \(ARMA\) se basan en el supuesto de estacionariedad del proceso […]. Sin embrago, muchas series de tiempo relacionadas con aplicaciones reales no son estacionarias, ya sea porque cambian de nivel en el tiempo (existe tendencia) o la varianza no es constante en el tiempo, a este tipo de procesos se les conoce como procesos integrados.
Para trabajar con estas series de tiempo lo que se hace es calcular las diferencias de la serie de tiempo \(d\) veces para hacerla estacionaria y posteriormente aplicar a la serie diferenciada un modelo \(ARMA(p,q)\) en este caso se diría que la serie original es un proceso \(ARIMA(p,d, q)\)
Y no esta demás recordar lo siguiente
Los modelos integrados se usan para reducir la NO estacionariedad de la serie al usar las diferentes.
Por lo que nos interesa como proponer un modelo \(ARMA\) o sus derivados.
La manera resumida se condensa en la siguiente tabla William and Wei (2006)
| Proceso | \(ACF\) | \(PACF\) |
|---|---|---|
| \(AR(p)\) | Decaimiento exponencial o amortiguamiento sinusoidal | Corte después del retraso \(p\) |
| \(MA(q)\) | Corte después del retraso \(q\) | Decaimiento exponencial o amortiguamiento sinusoidal |
| ARMA(p,q) | Sin cortes bruscos | Sin cortes bruscos |
\(AR(p)\)
Recordemos que en este proceso el pasado puede predecir el futuro; es decir, que dados ciertos \(p\) valores pasados (\(X_{t-1}, X_{t-2}, \dots, X_p\)), se conoce el valor \(X_t\).
Vamos a empezara utilizar algunas funciones de R y para ver como se comportan algunos procesos vamos a utilizar la función stats::arima.sim() la cual nos permitirá simular un proceso \(ARIMA\) y por lo tanto un \(AR\), \(MA\) y \(ARMA\) de manera adecuada. Cabe aclarar que se esta tomando los parámetros del ejemplo mostrado en el libro Holmes and Ward (2019).
set.seed(10)
#Diferentes coeficientes para los modelos AR
ARcoefficients <- c(0.7, 0.2, -0.1)
#Vamos a simular y guardar todo en una lista
AR_models <- list()
for(p in 1:3){
AR_models[[p]] <- stats::arima.sim(n = 10000, list(ar = ARcoefficients[1:p]))
}
#Utilizamos una función creada para este material
series_plot(AR_models[[1]], 100, TRUE) /
series_plot(AR_models[[2]], 100) /
series_plot(AR_models[[3]], 100, FALSE, TRUE)
Las gráficas corresponden a modelos \(AR(1)\), \(AR(2)\) y \(AR(3)\). Véase que no tenemos una información relevante por parte de los ACF pero sí por parte de los segundos correlogramas dando a notar que el número de correlaciones parciales significantes indican el parámetro \(p\) del modelo \(AR\). Podemos tomar la función de R stats::ar() la cual selecciona, mediante el criterio Akaike, el mejor modelo \(AR\) para la serie de tiempo dada.
ar(AR_models[[1]], method = "mle")
Call:
ar(x = AR_models[[1]], method = "mle")
Coefficients:
1
0.7055
Order selected 1 sigma^2 estimated as 1.011- Por último, se recomienda ver unos puntos interesantes sobre el modelo \(AR\) en este enlace.
\(MA(q)\)
Para el modelo \(MA(q)\) recordemos que estamos suponiendo que los valores pasados son independientes y causados por una fuente externa. Para este caso, vamos a hacer algo similar a modelo anterior para justificar de manera empírica la tabla con la que inicio esta sección.
set.seed(11)
#Diferentes coeficientes para los modelos AR
AMcoefficients <- c(0.7, 0.2, -0.1)
#Vamos a simular y guardar todo en una lista
AM_models <- list()
for(p in 1:3){
AM_models[[p]] <- stats::arima.sim(n = 10000, list(ma = ARcoefficients[1:p]))
}
series_plot(AM_models[[1]], 100, TRUE) /
series_plot(AM_models[[2]], 100) /
series_plot(AM_models[[3]], 100, FALSE, TRUE)
De manera similar, las gráficas corresponden a modelos \(MA(1)\), \(MA(2)\) y \(MA(3)\). En este caso, los correlogramas de los coeficientes parciales parecen no ofrecernos mucha información, ya que en este caso es el ACF quien nos proporciona una primera idea para modelar nuestros datos. ¿Por qué sucede esto?
Si lo pensamos, en un modelo \(AR\), los valores pasados influyen directamente con los valores más actuales, por lo que habrá una gran correlación en los retrasos que irán disminuyendo con el paso del número de retrasos y gracias a esa correlación que se va acumulando, el \(PACF\) es el indicado. En un modelo \(MA\), no se tiene dicha correlación y los eventos pasados son independientes entre sí, por lo que ver la correlación “real” entre distintos rezagos no es de gran ayuda, en este caso el \(ACF\) nos ayuda a ver cuantos rezagos significativos tienen influencia con el valor actual.
\(ARMA(p,q)\)
Para este punto es claro que el modelaje dependerá de su ACF como de su PACF. Vamos a crear dos modelos para ver dicho comportamiento.
set.seed(15)
ARIMAmodel1 <- arima.sim(model = list(order = c(1, 0, 1), ar = .7,ma =.2), n = 1000)
ARIMAmodel2 <- arima.sim(model = list(order = c(2, 0, 2), ar = c(.7, 0.2),ma = c(-0.1, -0.3)), n = 1000)
series_plot(ARIMAmodel1, 100, TRUE) /
series_plot(ARIMAmodel2, 100, FALSE, TRUE)
Como se puede notar, no hay algún rezago en el cual cambien drásticamente los valores del ACF o el PACF de manera secuencial. En general, si el comportamiento del ACF decrece lentamente a medida que los rezagos transcurren, se tiene un proceso no estacionario, y como ya se sabe, se utiliza un método \(ARIMA\) para solucionar este problema.
\(ARIMA(p,d,q)\)
Sólo para recordar, muchas de las series no cumplirán ser estacionarias, por lo que en la búsqueda de ajustar un modelo anterior optamos por aplicar diferencias para volver a la serie más estable y después aplicar un modelo \(ARMA(p,q)\) de algún tipo.
Recordemos también que los propios correlogramas (y en especial el ACF) nos darán indicios de tendencia y estacionalidad, por lo que desde este punto podemos partir para pensar en un modelo \(ARIMA\). Tomemos el caso de AirPassengers
Vamos a aplicar múltiples diferencias para ver el comportamiento de la serie resultante junto con sus gráficas correspondientes ACF y PACF. Para esto, se utilizará la función diff() la cual tiene tres argumentos: diff(x, lag = 1, differences = 1). El primero son los datos (la serie de tiempo), el segundo indica el número de rezagos para aplicar la diferencia, esto ayudará a eliminar la estacionalidad y el último indica el número de diferencias a aplicar, en el modelo \(ARIMA(p,d,q)\), este último argumento corresponde al parámetro \(d\), lo cual ayudará a eliminar la tendencia (con \(d = 1\) se eliminará la tendencia lineal y con \(d = 2\) se eliminará tendencia cuadrática).
AirPassengers_d1 <- diff(AirPassengers, differences = 1)
AirPassengers_d2 <- diff(AirPassengers, differences = 2)
diff_plots(AirPassengers, AirPassengers_d1, AirPassengers_d2)
Es rápido ver desde las nuevas series de tiempo que la varianza aumenta conforme el tiempo lo hace, más adelante veremos algunas transformaciones para solucionar esto, pero por el momento podemos ver que con la aplicación de 1 diferencia es suficiente. Ahora, desde la serie inicial se había notado una estacionalidad que el ACF indicaba que se daba cada 12 periodos, la cual se sigue manteniendo en estos nuevos modelos. Para tratar esto, véase como se modifican las funciones diff() indicando que se debe realizar una diferenciación cada \(m\) periodos, es decir: \(X_t^{'} = X_t-X_{t-m}\). Para aclarar, sólo se esta realizando una diferencia estacional, no una diferencia no estacional y estacional a la vez.
AirPassengers_d1_12 <- diff(AirPassengers, differences = 1, lag = 12)
AirPassengers_d2_12 <- diff(AirPassengers, differences = 2, lag = 12)
diff_plots(AirPassengers, AirPassengers_d1_12, AirPassengers_d2_12)
Ahora tenemos una serie temporal más estable, sin olvidar que tenemos algunos inconvenientes previos. Lo que es importante hacer notar, es el comportamiento que tenemos en los gráficos de correlación. Por una parte, en el ACF tenemos un decaimiento que, junto al PACF, nos indica un modelo \(AR(1)\), o bien un modelo \(AR(2)\), lo cual es más fácil de ver con la primera diferencia. ¿Cuántas diferencias necesitamos para volver a la serie estacionaria? eso lo podremos determinar creando las diferencias y utilizar alguna prueba de hipótesis para ver si la serie, con la diferencia aplicada, es estacionaria o no.
Por el momento, podemos utilizar la función forecast::ndiffs(), la cual realiza las etapas antes mencionadas y alguna prueba estadística de la que se hablará más adelante; y para determinar si era necesaria una diferenciación pero para eliminar la estacionalidad se puede utilizar la función forecast::nsdiffs().
cat(paste("Número de diferencias: ", ndiffs(AirPassengers),
"\nNúmero de diferencias estacionales: ", nsdiffs(AirPassengers)))Número de diferencias: 1
Número de diferencias estacionales: 1Véase que al aplicar alguna diferencia de algún tipo no se elimina el comportamiento estacional o no-estacional. Por ejemplo, véase los resultados obtenidos por las funciones ndiffs y nsdiffs.
diff(AirPassengers) %>% nsdiffs(): 1diff(AirPassengers, lag = 12) %>% ndiffs(): 1
Es decir, las funciones, mediante el uso de pruebas de hipótesis, dicen que se necesita una diferencia estacional de 12 periodos y otra de no estacional para volver a la serie estacionaria.
¿Qué podemos concluir de los resultados de las siguientes funciones?
AirPassengers %>% diff(lag = 12) %>% ndiffs()AirPassengers %>% diff(lag = 12) %>% nsdiffs()
Vamos a aplicar una diferencia a los datos de manera estacional para aplicar un modelo \(AR\). Obtener los coeficientes puede ser complicado, por fortuna podemos utilizar herramientas computacionales para esto, en este caso vamos a utilizar la función forecast::Arima(), la cual obtiene los coeficientes mediante métodos numéricos (?Arima). Vamos a obtener el ajuste para los modelos \(AR(1)\) y \(AR(2)\) (existe otra función del paquete stats pero puede causar problemas al tratar con modelos \(ARIMA\) con \(d\neq0\)).
forecast::Arima(diff(AirPassengers, differences = 1, lag = 12), order = c(1,0,0))Series: diff(AirPassengers, differences = 1, lag = 12)
ARIMA(1,0,0) with non-zero mean
Coefficients:
ar1 mean
0.7566 31.0184
s.e. 0.0569 4.0403
sigma^2 estimated as 135.3: log likelihood=-510.62
AIC=1027.25 AICc=1027.44 BIC=1035.9Entonces, bajo este método y los nuevos datos, nuestro modelo sería \(X_t = 31.0184 + 0.7566X_{t-1}+\epsilon_t\) el cual sería un \(ARIMA(1,0,0)\). No hay que olvidar que estos nuevos datos en realidad tienen una diferencia estaciona aplicada, por lo que la representación algebráica no es la correcta, más adelante se hablará de la expresión real.
forecast::Arima(diff(AirPassengers, differences = 1, lag = 12), order = c(2,0,0))Series: diff(AirPassengers, differences = 1, lag = 12)
ARIMA(2,0,0) with non-zero mean
Coefficients:
ar1 ar2 mean
0.5967 0.2130 30.5274
s.e. 0.0847 0.0852 4.9826
sigma^2 estimated as 130.1: log likelihood=-507.58
AIC=1023.16 AICc=1023.47 BIC=1034.69Ahora, el modelo sería \(X_t = 30.05724 + 0.5967X_{t-1}+0.2130X_{t-2}+\epsilon_t\) el cual sería un \(ARIMA(2,0,0)\); de igual manera, esta expresión sólo correspondería a los nuevos datos, más no a los originales. De acuerdo al estadístico \(AIC\), este modelo minimiza la pérdida de información así que se elegiría este.
Otra función muy útil para estos modelos es la función forecast::auto.arima(). Veamos un ejemplo para posteriormente dar una explicación
forecast::auto.arima(AirPassengers)Series: AirPassengers
ARIMA(2,1,1)(0,1,0)[12]
Coefficients:
ar1 ar2 ma1
0.5960 0.2143 -0.9819
s.e. 0.0888 0.0880 0.0292
sigma^2 estimated as 132.3: log likelihood=-504.92
AIC=1017.85 AICc=1018.17 BIC=1029.35Lo que se indica con ARIMA(2,1,1)(0,1,0)[12] es un modelo conocido como \(SARIMA(p,d,q)(P,D,Q)_m\) donde básicamente se trata un comportamiento de estacional en los datos. Este tipo de modelos es el que en realidad se ha hecho en los dos últimos modelos antes de la función auto.arima y de hecho, como mencionaban las funciones ndiffs y nsdiff, se necesitaban dos diferencias para volver a la serie estacionaria tal como el algoritmo lo esta determinando.
En el siguiente enlace se puede encontrar más información, por el momento es bueno saber que este tipo de modelos se pueden identificar con el ACF y PACF ya que se tendrán comportamientos similares para identificar un modelo \(ARMA\) pero ahora pensando que ahora habrá valores más grandes cada cierta cantidad de periodos (la estacionalidad). De hecho, se estuvo mencionando en los ejemplos anteriores que los modelos eran sobre los datos ya transformados, ya que al aplicar diferencias para reducir la estacionalidad estamos hablando de un modelo \(SARIMA\).
Véase que si aplicamos la misma función sobre la serie diferenciada considerando la estacionalidad, ya no aparece un modelo \(SARIMA\) en el output aunque, de hecho, se obtienen los mismos resultados. Otro punto que no esta de más, es que la diferencia tiene efecto sobre el parámetro \(D\) del modelo.
forecast::auto.arima(diff(AirPassengers, differences = 1, lag = 12))Series: diff(AirPassengers, differences = 1, lag = 12)
ARIMA(2,1,1)
Coefficients:
ar1 ar2 ma1
0.5960 0.2143 -0.9819
s.e. 0.0888 0.0880 0.0292
sigma^2 estimated as 132.3: log likelihood=-504.92
AIC=1017.85 AICc=1018.17 BIC=1029.35Otro punto que hay que notar es la aparición de un modelo \(ARIMA(2, 1, 1)\), contrario a nuestro modelo propuesto \(ARIMA(2,1,0)\) con índices \(AIC\) de 1017.85 y 1023.16, lo que indicaría que este modelo tiene una perdida menor de información con respecto al modelo \(ARIMA(2,1,0)\). Como bien se menciono, el ACF y PACF nos dan un primera idea sobre que modelo debemos debemos ajustar y es común colocar variaciones de esta propuesta y elegir aquel que tenga un mejor estadístico.
La función auto.arima() trabaja muy bien; trata de encontrar el mejor modelo utilizando pruebas de hipótesis para determinar la estacionariedad y minimizaciones al AIC y a la verosimilitud. Además de que se ajusta a los datos, los cuales podrían requerir algún tipo de tratamiento previo, puede que la función no obtenga el modelo con el AIC mínimo por el mismo hecho de las aproximaciones que realiza. Para más información puede consultarse el siguiente enlace.
Finalmente, sabemos que la expresión para un modelo \(ARIMA(p, d, q)\) se puede plasmar como
\[ X_t^d=c+\phi_1X_{t-1}^d+...+\phi_pX_{t-p}^d-\theta_1\epsilon_{t-1}^d-\theta_2\epsilon_{t-2}^d-...-\theta_q\epsilon_{t-q}^d+\epsilon_t \] Y en terminos de operadores de retardo de la siguiente manera
\[ \begin{array}{c} (1-\phi_1B-\phi_2B^2-...-\phi_pB^p)X_t^d=c+(1-\theta_1B-\theta_2B^2-...-\theta_qB^q)\epsilon_t\\ \phi_p(B)(1-B)^dX_t=c+\theta_q(B)\epsilon_t \end{array} \]
Y ya que el modelo \(SARIMA\) contiene retrocesos por la estacionalidad, este quedará expresado en términos de operados de retardo, expresión que bien podría llevarse a los términos comunes. Por ejemplo, para un modelo \(ARIMA(1, 1, 1, )(1, 1, 1)_[12]\) se tiene lo siguiente
\[ (1-\phi_1 B)(1-\Phi_1B^{12})(1-B)(1-B^{12})X_t = (1+\theta_1B)(1+\Theta_1B^{12})\epsilon_t \]
Donde los coeficientes \(\Phi\) y \(\Theta\) corresponden a la parte estacional de los modelos \(AR\) y \(MA\) de manera correspondiente.
Transformaciones
Nunca esta demás hacer notar lo siguiente: los métodos \(ARIMA(p,d,q)\) (y esto involucra todos los métodos lineales antes vistos) son modelos de ajuste lineal sobre los datos. Puede suceder que los datos (la serie de tiempo) sea muy caótica para que sea ajustada a alguno de los anteriores modelos o algún otro, por lo que muchas veces es mejor tratar los datos con algún tipo de transformación y a esos nuevos datos aplicar el modelo \(ARIMA\) y no olvidar que si se desea predecir, se esta realizando una predicción sobre esa transformación, por lo que se podría aplicar una transformación inversa para tener una correcta predicción.
Estas transformaciones son utilizadas para reducir algún comportamiento en la serie que impida que esta sea estacionaria. Antes, se vio que las diferencias ayudan a la serie con la media pero en una de las gráficas anteriores se notó que la varianza aumentaba a medida que transcurría el tiempo, o el general era grande. Para evitar este tipo de comportamientos se sugiere utilizar alguna de las siguientes transformaciones:
- Logaritmo: \(X^{'}_t = \log(X_t)\)
- Potencias: \(X^{'}_t = (X_t)^p\)
- Transformación Box-Cox: \(X_t^{'} = \begin{cases}\log(X_t) & \mbox{si } \lambda = 0\\(X_t^\lambda-1)/\lambda & \mbox{e.o.c }\end{cases}\)
Véase que esta última utiliza las primeras dos transformaciones; ya depende de lo que requiera el problema, pero se puede utilizar cualquier de estas para obtener un mejor modelo. No se adentrará más a este tema pero se sugiere leer el siguiente capitulo del libro Hyndman and Athanasopoulos (2018) y las funciones útiles en este caso son forecast::BoxCox() y forecast::BoxCox.lambda(). Finalmente, no hay que olvidar condiciones implícitas antes de aplicar una transformación.
Vamos a realizar la primera transformación sobre los datos para reducir la variación en los datos y veamos que efecto tiene sobre el modelo.
subplot(ggplotly(autoplot(AirPassengers) + general_theme + labs(x = "Tiempo", y = NULL)),
ggplotly(autoplot(log(AirPassengers)) + general_theme + ggtitle("AirPassengers") + labs(x = "Tiempo", y = NULL)), nrows = 2, shareX = T)La tendencia aún persiste, así que vamos a determinar las diferencias que debemos aplicar para mejorar la estabilidad de la serie
cat(paste("Número de diferencias: ", ndiffs(log(AirPassengers)),
"\nNúmero de diferencias estacionales: ", nsdiffs(log(AirPassengers))))Número de diferencias: 1
Número de diferencias estacionales: 1Así que vamos a aplicar dicho modelo con las diferencias sugeridas
AirPassengers_d1_12_log <- diff(log(AirPassengers), differences = 1, lag = 12) %>% diff()
diff_plots(log(AirPassengers), AirPassengers_d1_12_log)
Y de manera similar, viendo las gráficas para la serie diferenciada de manera estacional podemos proponer un modelo \(AR(p)\), aunque ya sabiendo que el comportamiento en el decaimiento del ACF con indicios de estacionalidad sugiere un modelo \(SARIMA(p,1,0)(0,1,1)_{12}\). Vamos a sugerir varios valores para \(p\) y a obtener el modelo con el \(AIC\) más pequeño, esto sin modificar alguno de los otros parámetros del modelo.
AP_Sarima <- function(p){
forecast::Arima(diff(log(AirPassengers), lag = 12) %>% diff(), order = c(p, 0, 0), seasonal = c(0,0,1))
}
#Diferentes modelos
AP_ModelsSarima <- map(1:4, ~AP_Sarima(.x))
which.min(map_dbl(AP_ModelsSarima, ~`[[`(.x, "aic")))[1] 4Este resultado es curioso, ya que podía intuirse desde la serie sin diferencias, que el ACF y PACF indican un modelo \(AR\), después se nota que cada 12 periodos hay picos en ambas gráficas por lo que indica que se necesita aplicar diferencias estacionales de 12 periodos y el ACF indica tendencia por lo que es necesaria una diferencia no estacional. Entonces veamos el modelo con 4
AP_ModelsSarima[[4]]Series: diff(log(AirPassengers), lag = 12) %>% diff()
ARIMA(4,0,0)(0,0,1)[12] with non-zero mean
Coefficients:
ar1 ar2 ar3 ar4 sma1 mean
-0.3865 -0.1195 -0.1704 -0.1808 -0.5697 -2e-04
s.e. 0.0865 0.0911 0.0912 0.0866 0.0753 9e-04
sigma^2 estimated as 0.001362: log likelihood=246.91
AIC=-479.82 AICc=-478.91 BIC=-459.69Por cuestiones de interpretabilidad, veámos también para \(p = 2\)
AP_ModelsSarima[[2]]Series: diff(log(AirPassengers), lag = 12) %>% diff()
ARIMA(2,0,0)(0,0,1)[12] with non-zero mean
Coefficients:
ar1 ar2 sma1 mean
-0.3616 -0.0637 -0.5617 -0.0002
s.e. 0.0875 0.0870 0.0738 0.0011
sigma^2 estimated as 0.001404: log likelihood=244.02
AIC=-478.04 AICc=-477.56 BIC=-463.66La diferencia entre los modelos no es grande, así que vamos a seguir con el modelo con \(p = 2\), por lo que tendríamos un modelo \(SARIMA(2, 1, 0)(0,1,1)_{12}\). ¿Qué hubiera sucedido si utilizamos auto.arima()?
(good_model_AP <- auto.arima(log(AirPassengers)))Series: log(AirPassengers)
ARIMA(0,1,1)(0,1,1)[12]
Coefficients:
ma1 sma1
-0.4018 -0.5569
s.e. 0.0896 0.0731
sigma^2 estimated as 0.001371: log likelihood=244.7
AIC=-483.4 AICc=-483.21 BIC=-474.77Obtenemos un modelo \(SARIMA(0, 1, 1)(0,1,1)_{12}\), el cual tiene un AIC más bajo que nuestro modelo anterior y mucho menor al modelo sin la aplicación del logaritmo sobre los datos. Otro punto importante es que tenemos un modelo \(MA\) presente en lugar de un \(AR\) como nos había sugerido la intuición. Ay que recordar que podemos escribir un modelo \(AR\) como una suma infinita de proceso \(MA\), aunque no parece ser el caso. Hay que recordar que la función auto.arima() aplica métodos numéricos y bajo diferentes modelos, por lo que podríamos haber obtenido este modelo sí hubiéramos propuesto algunas variantes a nuestro modelo propuesto.
- ¿Puede el modelo dado por
auto.arima()no ser el mejor? Revisen el post Why doesn’t auto.arima() return the model with the lowest AICc value?. - ¿Qué hubiera pasado si hubiéramos aplicado una transformación Box-Cox?
Evaluación y pruebas de hipótesis
Cómo cualquier modelo supervisado, no queremos tener problemas con el overfitting y el underfitting por lo que podemos aplicar aquí las técnicas comunes para validar nuestro modelo, ya sea tomando un porcentaje de los datos, aprender de ellos, crear el modelo y después evaluar con nuestras predicciones el remanente de la información; esto último se tratará en la sección de Predicción y por el momento se hablará sólo de los modelos lineales utilizados.
Como bien se menciona en el libro de referencia, en el capítulo 10 se menciona lo siguiente
Verificación: Si el modelo es el adecuado, es decir los valores de \(p\) y \(q\) han sido correctamente especificados, entonces el modelo deberá ajustar bien a los datos y los residuales (la diferencia entre lo observado y lo estimado con el modelo) deberán comportarse como ruido blanco, esto se puede comprobar con la prueba de Ljung-Box. Si el modelo es adecuado, la función de autocorrelación de los residuales no debe de tener ninguna estructura[…].
Es evidente que los residuales están expresados por \(\epsilon = X_t-\hat{X}_t\) donde \(\hat{X}_t\) es el modelo. La siguiente gráfica muestra como el modelo se ajusta a los datos.
autoplot(log(AirPassengers), color = "black") +
autolayer(fitted(good_model_AP), color = "royalblue1") +
labs(y = "log(Número de pasajeros)", x = "Tiempo") +
general_theme
Dado esto, vamos a obtener los residuales de la serie AirPassengers y el mejor modelo que hemos conseguido: \(ARIMA(0, 1, 1)(0,1,1)_{12}\), esto considerando la serie log(AirPassengers). La propia función auto.arima() devuelve una lista y entre sus elementos se encuentran los residuales entre los datos y el modelo, aunque igual se puede utilizar la función stats::residuals(model) para obtenerlos.
(autoplot(good_model_AP$residuals) + labs(x = NULL, y = NULL) + ggtitle(TeX("Residuales del modelo $ARIMA(0,1,1)(0,1,1)_{12}$")) + general_theme) /
((ggAcf(good_model_AP$residuals) + ggtitle(NULL) + general_theme) + (gghistogram(good_model_AP$residuals) + labs(x = "Residuales", y = "Conteo") + general_theme))
forecast::checkresiduals(good_model_AP, plot = F)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
Q* = 26.446, df = 22, p-value = 0.233
Model df: 2. Total lags used: 24La anterior gráfica y la prueba se pudieron haber obtenido con la anterior función forecast::checkresiduals(good_model_AP), pero sin modificar el argumento plot o test.
La función anterior nos da un análisis gráfico de los residuales junto con la prueba de hipótesis Ljung-Box, la cual evalúa por grupo los datos (en este caso los residuales) y determina si estos son independientes entre sí; como bien se menciona en las notas de referencia, bajo esta hipótesis, las autocorrelaciones deberán ser simultáneamente 0. En este caso el p-value es lo suficientemente grande, por lo que no tenemos suficiente información para rechazar que los residuales sean simultáneamente cero, además podemos decir lo siguiente:
- Los residuales no son distinguibles de un ruido blanco.
- No es posible rechazar que los residuales no están correlacionados.
Dado esto, podemos decir que nuestro modelo se ajusta bien a nuestros datos. Es importante mencionar que la prueba Ljung-Box es una mejora a la prueba Box-Pierce, pero aún así se puede utilizar la anterior prueba si se desea. Ambas pruebas se pueden obtener del modelo con la función Box.test(). A este tipo de pruebas se les conoce como Portmanteau test.
Antes se había mencionado que las funciones auto.arima() y ndiffs() utilizan una prueba de hipótesis sobre la estacionariedad de la serie. A este tipo de pruebas se les llama de pruebas de raíz unitaria donde básicamente se busca que la raíces de la ecuación característica sean distintas de 1. Puede consultarse más sobre esto en los enlaces dejados en este párrafo pero la idea central radica en las raíces que se obtienen del polinomio derivado de utilizar los operadores de retardo
\[ \phi_p(B)(1-B)^dX_t=c+\theta_q(B)\epsilon_t^d \]
Cómo ya se estudio, las condiciones para que los modelos sean estacionarios es que los coeficientes cumplan ciertas restricciones en su dominio y básicamente tienen que ser menor a 1, lo cual hace sentido para que este tipo de pruebas determinen estadísticamente si la serie es estacionaria o no. Una de estas pruebas es la prueba Augmented Dickey-Fuller, la cual, como menciona su nombre, es una versión aumentada de la prueba Dickey-Fuller, ya que esta última esta definida para un modelo autorregresivo.
De acuerdo a la versión, la prueba se puede establecer de la diferentes maneras, pero en este caso queda establecida de la siguiente manera:
\[
\begin{array}{cc}
H_0:\mbox{Hay una raíz unitaria}\\
H_1: \mbox{La serie de tiempo es estacionaria}
\end{array}
\] Vamos a utilizar la función urca::ur.df() para realizar esta prueba, aunque bien podría utilizarse la función tseries::adf.test(). ¿En que datos vamos a realizar la prueba? El modelo que fue otorgado por el modelo auto.arima() ya es estacionario, así que vamos a realizar esta prueba sobre los residuales ya que nos interesa que estos se comporten como ruido blanco, el cual es estacionario.
library(urca)
urca::ur.df(good_model_AP$residuals) %>% summary()
###############################################
# 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.11724 -0.01909 -0.00037 0.02217 0.10772
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.96169 0.11893 -8.086 2.63e-13 ***
z.diff.lag -0.02387 0.08472 -0.282 0.779
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03553 on 140 degrees of freedom
Multiple R-squared: 0.4926, Adjusted R-squared: 0.4854
F-statistic: 67.96 on 2 and 140 DF, p-value: < 2.2e-16
Value of test-statistic is: -8.0858
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62Entonces, el ruido es estacionario. ¿Que resultados esperaríamos de urca::ur.df(log(AirPassengers))?
Otra prueba importante es la prueba KPSS (Kwiatkowski-Phillips-Schmidt-Shin) en la que también se busca dar evidencia estadística de la estacionariedad de la serie. En este caso, las hipótesis no están sujetas a las raíces de la ecuación característica, si no que la hipótesis nula establece que la serie es estacionaria y la hipótesis alternativa lo contrario. Puede leerse más sobre esto aquí o en la publicación original de la prueba. Esta prueba la podemos aplicar con la función urca::ur.kpss() (existe también tseries::kpss.test()).
urca::ur.kpss(good_model_AP$residuals) %>% summary()
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 4 lags.
Value of test-statistic is: 0.1491
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739De acuerdo a la prueba, con los residuales se tiene un estadístico menor al marcado por un nivel de significancia del 10%, por lo que el estadística no está dentro de la región de rechazo y por lo tanto no tenemos evidencia suficiente para rechazar que la serie es estacionara, y en este caso los residuales.
- ¿Qué sucede al ejecutar
autoplot(good_model_AP)? - En el siguiente enlace se pueden encontrar referencias a más pruebas de raíz unitaria junto una pequeña explicación de este tipo de pruebas.
Predicción
El objetivo de un modelo lineal será predecir, esto aprendiendo de la información obtenida que servirá para proponer valores futuros. En este caso, ya que se tiene el modelo que se ajusta a la serie de tiempo, ahora se aplicarán las funciones correspondientes para realizar predicciones, lo cual será sencillo y más en R.
good_model_AP %>% forecast(h = 30) %>% autoplot() + general_theme + labs(x = "Tiempo") + ggtitle(TeX("Pronóstico desde $ARIMA(0,1,1)(0,1,1)_{12}$"))
Obviamente hay que aplicar el efecto inverso al logaritmo ya que el modelo esta en dichos términos
prediction <- good_model_AP %>% forecast(h = 30)
prediction$x <- exp(prediction$x)
prediction$mean <- exp(prediction$mean)
prediction$lower <- exp(prediction$lower)
prediction$upper <- exp(prediction$upper)
autoplot(prediction) + general_theme + labs(x = "Tiempo", y = "AirPassengers") + ggtitle(TeX("Pronóstico desde $ARIMA(0,1,1)(0,1,1)_{12}$"))
- Aún faltan muchos métodos por ver para otorgar una predicción.
- ¿Cómo funciona el algoritmo para dar esta predicción?
- ¿Qué tan buena es esta predicción? Obviamente los intervalos de predicción aumentan conforme los periodos avanzan.
- ¿Los intervalos de predicción llevan algún supuesto?
- ¿No estaremos realizando overfitting? Hay que recordar que el mejor modelo no se basa en nuestra intuición, tal vez ajustando un modelo a una parte de los datos tengamos una mayor certeza de la predicción y de un modelo más adecuado y sencillo de interpretar.
- Revisen las funciones
window(),head(),tail()ysubset()para realizar filtros sobre objetostsy la funciónforecast::accuracy()evaluar la predicción.
Como comentario final, este es sólo un tipo de modelos, actualmente existen muchos más que sólo modelos lineales, aún así este tipo de modelos tienen gran popularidad debido a su sencilla interpretación y ejecución, aunque nunca estará de más obtener el conocimiento y adaptar mejores modelos a nuestros datos.
Bibliografía
Cryer, Jonathan D, and Kung-Sik Chan. 2008. Time Series Analysis: With Applications in R. Springer Science & Business Media.
Holmes, E, and EJ Ward. 2019. Applied Time Series Analysis for Fisheries and Environmental Sciences. University of Washington, Lecture Material.
Hyndman, Rob J, and George Athanasopoulos. 2018. Forecasting: Principles and Practice. OTexts.
Nielsen, Aileen. 2019. Practical Time Series Analysis: Prediction with Statistics and Machine Learning. " O’Reilly Media, Inc.".
William, WS, and S Wei. 2006. Time Series Analysis: Univariate and Multivariate Methods. Second. Pearson/Addison-Wesley Reading.
A work by Carlos Vásquez
carlosfvasquez@ciencias.unam.mx