SeriesTCode
Vásquez Guerra Carlos Fernando
08/2024
Datos
R es un lenguaje de programación en el que uno de sus paradigmas es
el orientado a objetos y es común encontrar objetos para distintos tipos
de temas, 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, es relevante mencionar que para la creación de este
contenido se tomó como referencia el libro virtual Forecasting: Principles and Practice
(2nd ed). Actualmente existe una tercera edición en la cual se tienen
cambios importantes sobre el uso de paqueterías de R y se deja al lector
explorar este contenido. El propio autor del anterior creó el paquete forecast de R
y por esta razón se utiliza la segunda edición. En diversas ocasiones se
colocarán hipervínculos con contenido adicional en la lectura por lo que
se recomienda explorar el contenido referenciado, podrías distinguir
estos enlaces cuando un conjunto de palabras tenga un formato similar a
este formato.
También, 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().
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. 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 tiene 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.
Recordemos que el modelo aditivo supone que la serie de tiempo tiene la siguiente estructura: \(X_t = T_t+E_t+I_t\). En este caso, la forma clásica de descomposición consiste en obtener la tendencia mediante un modelo de suavizamiento de medias móviles, 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 veamos 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.
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
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 & más
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 Cleveland et al. (n.d.). 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().
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 semanales del coronavirus de Our World in Data y vamos a limitarnos 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 como estamos instanciando el parámetro start y
frequency. Este último lo conocemos ya que esa es la
frecuencia con la que se dan a conocer tales datos. %U es
el código para dar el formato de
la semana del año.
tseries <- ts(coronavirus_Mexico$new_cases,
start = c(2020, as.numeric(format(coronavirus_Mexico$date[1], "%U"))),
frequency = 52)Pero, ¿Qué sucede al ejecutar el comando stl(tseries)?
(sucede algo similar al realizar decompose(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 42 semanas de un año. Para solucionar este inconveniente no
utilizaremos la descomposició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 se ha mencionado.
Para estos en particular, supongamos que los datos tienen una
estacionalidad mensual, por lo que al tener 12 meses, el parámetro \(K\) mostrado en dicho modelo debe cumplir
que se \(K<=12/2 = 6\), 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, 6))[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.
Call:
tslm(formula = tseries ~ trend + fourier(tseries, min_casesT))
Coefficients:
(Intercept) trend
-47965.6 3691.8
fourier(tseries, min_casesT)S1-52 fourier(tseries, min_casesT)C1-52
33961.0 18017.6
fourier(tseries, min_casesT)S2-52 fourier(tseries, min_casesT)C2-52
-4793.9 16976.2
fourier(tseries, min_casesT)S3-52 fourier(tseries, min_casesT)C3-52
-937.3 6623.1
fourier(tseries, min_casesT)S4-52 fourier(tseries, min_casesT)C4-52
-3042.4 1709.4 Entonces, como la tendencia la estamos considerando como un modelo lineal, tomamos los valores correspondientes del modelo anterior 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] -44273.84 -40582.04 -36890.23 -33198.43 -29506.62 -25814.82Ahora sí, vamos a crear nuevamente el gráfico para ver nuestros resultados. Véase como se construyo 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()¿Qué está sucediendo? La tendencia no es del todo precisa para algunos periodos y en la estacionalidad no se aprecia algún ciclo o comportamiento relevante.
Hay que reconocer que el gran supuesto que estamos haciendo es que la tendencia se comporta de manera lineal. 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-4-12, 2020-7-26, 2020-9-20 y 2020-11-22 (semanas en que las existen cambios de dirección relevantes u algún otro comportamiento a considerar).
t <- time(tseries)
#Puntos de corte
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-4-12")
t.break1 <- time(tseries)[8]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-7-26")
t.break2 <- time(tseries)[23]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-9-20")
t.break3 <- time(tseries)[31]
# coronavirus_Mexico %>% filter(!is.na(new_cases)) %>% mutate(index = row_number()) %>% filter(date == "2020-11-22")
t.break4 <- time(tseries)[40]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
Time Series:
Start = c(2020, 8)
End = c(2020, 17)
Frequency = 52
[1] -1010.5523 -370.2226 270.1071 910.4369 1550.7666 2191.0963
[7] 2831.4260 3471.7558 6644.7048 9817.6538Y 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, 6)))
(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 un poco su varianza, lo cual es un punto a favor.
Algunas veces 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 realizar un análisis de series de tiempo, haciendo que este busque otro tipo de técnicas para dará solución a un problema dado.
¿Qué sucede si volvemos a tomar hacer la descomposición de la serie el día de hoy?
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 hasta 2020-12-09',
theme = theme(plot.title = element_text(hjust = 0.5))) +
general_theme
¿Qué puedes concluir de las anteriores gráficas?
Algunas 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) 
Al mencionar que las lineas azules indican si las correlaciones son significativamente diferentes de cero es equivalente a buscar que las correlaciones tengan un valor relevante y que no sea generado por pura aleatoriedad, como pasaría en el ruido blanco. Recordemos que el ruido blanco \(\epsilon_t\sim N(0, \sigma^2)\). Sin perdida de generalidad podemos considerar el caso \(\epsilon_t\sim N(0, 1)\).
\[ ci \leftarrow \frac{z_{1-\alpha}}{\sqrt{n}} \]
Si se accede al código
de la función autoplot.acf(), la cual es invocada cuando se
utiliza la función ggAcf() o ggPacf(), y con
el comentario anterior podemos entender como se calculan los intervalos
de confianza:
autoplot.acf <- function(object, ci=0.95, ...) {
...
# Add ci lines (assuming white noise input)
ci <- qnorm((1 + ci) / 2) / sqrt(object$n.used)
p <- p + ggplot2::geom_hline(yintercept = c(-ci, ci), colour = "blue", linetype = "dashed")Más adelante se verán las pruebas de hipótesis pertinentes a este tema, aunque es buen momento para mostrar el uso de la prueba de Ljung-Box. Recordemos que la prueba de hipótesis establece lo siguiente:
\[ H_0: \mbox{Los coeficientes de autocorrelación son simultáneamente iguales a cero} \\ H_1: \mbox{Alguno de los coeficientes de autocorrelación es distinto de cero} \]
Primero, veamos lo correspondiente al caso de COVID
Box-Ljung test
data: tseries
X-squared = 36.609, df = 1, p-value = 1.444e-09Y ahora al ruido blanco
Box-Ljung test
data: whitenoise$y
X-squared = 3.5137, df = 1, p-value = 0.06086Lo cual es consistente con las gráficas.
- 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
La clara disminución en los valores de los coeficientes es debido a la tendencia que tiene la serie y aquellos valores que parecen picos, 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 + theme(axis.ticks = element_blank())) /
((tibble(x = 1:100) %>% ggAcf()+
ggtitle("Autocorrelación simple") +
labs(x = "Lag", y = "ACF") +
general_theme + theme(axis.ticks = element_blank())) +
(tibble(x = 1:100) %>% ggPacf()+
ggtitle("Autocorrelación parcial") +
scale_y_continuous(name = "PACF", position = 'right') +
labs(x = "Lag") +
general_theme + theme(axis.ticks = element_blank())))
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 (PACF), 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 9 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.seeds <- c(10,11,1)
#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){
set.seed(set.seeds[p])
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 simluados \(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.
Call:
ar(x = AR_models[[1]], method = "mle")
Coefficients:
1
0.7055
Order selected 1 sigma^2 estimated as 1.011Otro punto interesante es lo siguiente: En los modelos simulados \(AR(1)\) y \(AR(2)\), en el PACF, los coeficientes significativos son adyacentes pero en el caso del modelo \(AR(3)\) tenemos 2 coeficientes significativos: en lag = 1 y lag = 3, por eso tenemos un \(AR(3)\), aunque la correlación parcial con lag = 2 no sea fuerte, la tenemos que considerar para darle peso al tercer rezago. Veamos que si, mediante métodos computacionales, ajustamos un modelo \(AR(3)\) obtenemos los coeficientes esperados:
Call:
arima(x = AR_models[[3]], order = c(3, 0, 0), method = "ML")
Coefficients:
ar1 ar2 ar3 intercept
0.7128 0.2011 -0.1161 -0.038
s.e. 0.0099 0.0121 0.0099 0.050
sigma^2 estimated as 1.024: log likelihood = -14307.82, aic = 28625.63Si pensamos que el hiperparámetro del modelo queda determinado por el número de rezagos significativos más que por el último rezago signficativo, propondríamos un modelo \(AR(2)\)
Call:
arima(x = AR_models[[3]], order = c(2, 0, 0), method = "ML")
Coefficients:
ar1 ar2 intercept
0.6989 0.1200 -0.0380
s.e. 0.0099 0.0099 0.0562
sigma^2 estimated as 1.038: log likelihood = -14375.7, aic = 28759.39Las métricas de ajuste nos indican que el modelo \(AR(3)\) tiene un mejor ajuste a la información. Es importante mencionar que esta heurística de ajuste solo es una opción que no garantiza el mejor ajuste, ya sea de manera interna o al momento de realizar predicciones.
- 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.seeds <- c(10,14,9)
#Diferentes coeficientes para los modelos AM
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){
set.seed(set.seeds[p])
AM_models[[p]] <- stats::arima.sim(n = 10000, list(ma = AMcoefficients[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\)).
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 = 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.
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 = 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
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 = 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.
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 = 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.
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
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 = 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\)
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 = 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()?
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 = 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
Como bien se menciona en el libro de referencia, en el capítulo 12 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[…].
Además de lo anterior debemos considerar poner a prueba la estacionariedad de la serie. A manera de resumen puedes consultar la siguiente liga en donde encontrarás una tabla con las principales pruebas de hipótesis en series de tiempo univariadas. En este capítulo se verán los puntos más relevantes de estas pruebas
Autocorrelación
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 + theme(axis.ticks = element_blank())) /
((ggAcf(good_model_AP$residuals) + ggtitle(NULL) + general_theme + theme(axis.ticks = element_blank())) + (gghistogram(good_model_AP$residuals) + labs(x = "Residuales", y = "Conteo") + general_theme))
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.
Podría utilizarse otro tipo de pruebas como la prueba Durbin-Watson o la prueba Breush-Godfrey, aunque hay que tener en cuenta lo siguiente:
- La prueba Durbin-Watson es solo válida para el modelo autorregresivo de primer orden (\(AR(1)\))
- La prueba Breush-Godfrey está diseñada específicamente para modelos
de regresión. De hecho, cuando se utiliza la función
forecast::checkresiduals()en un objeto de regresión, automáticamente aplicará la prueba Breush-Godfrey en lugar de la prueba Ljung-Box sobre los residuales.
Estacionariedad e invertibilidad
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 manera rápida de ver esto es con el uso de la
función forecast::autoplot(model) en donde se grafican las
raíces de la ecuación característica inversas, de esta
manera podemos apreciar visualmente la estacionariedad de la serie y la
invertibilidad de la misma ya que, como se mencionó en el capítulo 8 del
libro de referencia:
Un proceso \(ARMA(p, q)\) es estacionario si y sólo si el modulo de las raíces del polinomio \(\phi_p(B)\) está fuera del círculo unitario. Un proceso \(ARMA(p, q)\) es invertible si y sólo si el modulo de las raíces del polinomio \(\theta_q(B)\) está fuera del círculo unitario.
Veamos como se visualiza lo anterior con nuestro mejor modelo
Recordemos que nuestro mejor modelo no tiene un factor autorregresivo:
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 = 0.001371: log likelihood = 244.7
AIC=-483.4 AICc=-483.21 BIC=-474.77Solo para complementar, veamos como se ve con un modelo anterior:
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 = 132.3: log likelihood = -504.92
AIC=1017.85 AICc=1018.17 BIC=1029.35
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 y no para modelos más complejos.
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.
###############################################
# 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.62Como puedes ver, la función nos arrojó un resumen similar a una regresión, eso sucede por que la prueba ADF re interpreta al modelo de series de tiempo como un modelo lineal y realiza un ajuste lineal. Considerando un modelo AR, este se puede escribir como \(\Delta y_t=(p-1)y_{t-1}+u_t=\delta y_{t-1}+u_t\), lo cual haría equivalente determinar la raíz unitaria a poner a prueba si \(\delta=0\). Ya sea que se vean los valores críticos o el p-value, el modelo es estacionario.
¿Que resultados esperaríamos de
urca::ur.df(log(AirPassengers))?
La prueba Dickey-Fuller esta diseñada para un modelo autorregresivo y ajusta una regresión por mínimos cuadrados ordinarios, por lo que la correlación serial presentará un problema. La regresión de la prueba Dickey-Fuller aumentada incluye rezagos de las primeras diferencias de \(y_t\)
Otra prueba de este estilo es la prueba Phillips-Perron. Esta también ajusta una regresión pero realiza una corrección no paramétrica del estadístico de la prueba \(t\). Esta es robusta con respecto a la autocorrelación no especificada y la heterocedasticidad en el proceso de perturbación de la ecuación de prueba.
Phillips-Perron Unit Root Test
alternative: stationary
Type 1: no drift no trend
lag Z_rho p.value
4 0.217 0.738
-----
Type 2: with drift no trend
lag Z_rho p.value
4 -4.7 0.478
-----
Type 3: with drift and trend
lag Z_rho p.value
4 -46.8 0.01
---------------
Note: p-value = 0.01 means p.value <= 0.01 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, la hipótesis nula no está sujeta a las raíces de la ecuación
característica, si no que la hipótesis nula establece que la serie es
estacionaria. 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()).
#######################
# 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.739En este caso, véase que la prueba devuelve un estadístico y no un p-value. De 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.
- 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.
- Breusch–Godfrey test.
- How to Perform a Durbin-Watson Test in R.
- Using the Durbin-Watson (DW) test for testing time-series autocorrelation.
- Phillips–Perron unit-root test. Stata documentation.
- Phillips–Perron test.
- pp.test: Phillips-Perron Test
- KPSS test.
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::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::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.
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; puedes encontrar información sobre esto en el siguiente enlace en donde podrás entender como aplicar una separación de datos (train-test) en series de tiempo y hacer validación cruzada en este tipo de informació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
A work by Carlos Vásquez
carlosfvasquez@ciencias.unam.mx