Presentación

Esta publicación tiene dos objetivos, uno personal y otro general. El objetivo personal es poder tener un apunte en limpio sobre Tidy modeling y la librería seasonal. El segundo objetivo, el general, es simplemente compartirlo a todo aquel que guste.

No asumo que mi código sea el más eficiente o esté libre de errores, así que cualquier corrección, sugerencia o comentario házmelo saber por este medio o a través de mi twitter ecodiegoale.

En este blog haré uso de funciones personalizadas, si te interesa profundizar sobre la programación de funciones, te recomiendo leer R for Data Science de H. Wickham y G. Grolemund, disponible en este enlace. Asimismo, recomiendo revisar Tidy Modeling with R de Max Kuhn y Julia Silge (click en este enlace), pues la idea de este blog surgió de su lectura.

Además, nunca está de más revisar stackoverflow, github, rpubs o preguntarle a ChatGPT.

Aviso: Si algún link no abre con click izquierdo, trata con el derecho y abrir en otra ventana o pestaña.

Contenido

En este blog usaré la API de INEGI de una manera muy sencilla utilizando principalmente la librería de tidyverse. No me detendré a explicar mucho porque este tema ya lo traté en una entrada pasada, así que si no estás familiriazado te recomiendo revisar mi blog.

El principal objetivo de este blog es desestacionalizar múltiples series al mismo tiempo usando seasonal. El motivo de usar específicamente esta biblioteca es que es una interfaz fácil de usar para X-13-ARIMA-SEATS, el software de ajuste estacional de la Oficina del Censo de EE. UU., el cual ofrece acceso completo a casi todas las opciones y salidas del método X-13, incluyendo X-11 y SEATS, búsqueda automática de modelos ARIMA, detección de valores atípicos y soporte para variables de vacaciones definidas por el usuario.

Finalmente, como elemento adicional visualizo los datos procesados con la paleta de colores MetroBrewer (click en este enlace).

Preparación

Lo primero será definir tu directorio de trabajo y cargar las librerías necesarias para el análisis.

#setwd("C:/Users/tu carpeta")

rm(list = ls()) #para borrar todos los objetos en la memoria de la consola
options(scipen=999) #para desactivar la notación científica

#Si no cuentas con los paquetes, usa el comando install.packages("ejemplo")

library(inegiR)
library(tidyverse)
library(seasonal)
library(MetroBrewer)

Conforme avancemos en el blog haré mención de cuando estemos usando las librerías cargadas.

Construcción del tibble con la API de INEGI

El primer paso es descargar las series y construir el tibble que usaremos. Vamos a descargar las series originales del indicador mensual de actividad industrial por entidad federativa. ¿Qué significa que sean las series “originales”? Cuando se dice que una serie es “original” y no desestacionalizada, significa que los datos no han sido ajustados para eliminar las variaciones estacionales, y se presentan tal como se recopilaron inicialmente. Una variación estacional es un patrón de cambio en los datos que ocurre regularmente en un período específico de tiempo, como por ejemplo durante una temporada del año, un mes en particular, o incluso en días específicos de la semana. Estas variaciones están asociadas con factores estacionales como el clima, festividades, eventos culturales, o patrones de comportamiento humano que se repiten de manera predecible.

¿Qué mide y cómo se mide el indicador mensual de la actividad industrial de INEGI?

De acuerdo con el INEGI:

El Indicador Mensual de la Actividad Industrial por Entidad Federativa (IMAIEF) proporciona información estadística de corto plazo en el ámbito estatal que permite dar seguimiento al comportamiento de las actividades secundarias en los estados. Los datos se publican en forma de índices de volumen físico mensuales y acumulados y sus respectivas variaciones anuales, así como la contribución al crecimiento mensual tanto de cada estado al nacional para cada actividad, como de las actividades económicas en cada estado.

Así, queda claro que este es un indicador de corto plazo de cantidades reales de productos producidos en las actividades secundarias de cada entidad federativa. El término volumen físico contrasta con el valor monetario de la producción, ya que se enfoca específicamente en la cantidad física de bienes producidos, en lugar del valor monetario que representan. Ahora, siguiendo una lectura de Sistema de cuentas nacionales de México, año base 2018. Fuentes y metodologías de INEGI (click en este enlace), las actividades secundarias que se contemplan son Minería; generación y distribución de Energía eléctrica, suministro de agua y gas por ductos al consumidor final; Construcción e Industrias manufactureras. En este blog me enfocaré únicamente en las industrias manufactureras, las cuales incluyen actividades como industria, alimentaria, de bebidas, de madera, de papel, química, etc.

La construcción del indicador es por el método de la producción, el cual se basa en la suma de la producción de bienes finales dentro de la actividad económica específica. En esencia, el indicador se puede expresar como una división del valor del bien producido en el año de base (precio promedio del año base por cantidad actual) entre el valor medio de dicho bien (precio promedio del año base por la cantidad media del año base).

Para este blog solo usaremos el índice correspondiente a las actividades manufactureras.

Procesamiento de los datos

Como dije arriba, ya tengo una entrada dónde explico detalladamente cómo usar la API de INEGI. La API no es la única forma existente para bajar las series, puedes hacerlo directamente del Banco de Información Económica (BIE) del INEGI, yo uso la API para hacer este código reproducible y porque considero es una herramienta que vale la pena aprender a usar.

Aquí puedes revisar la entrada de mi blog donde lo explico.

Lo primero es guardar tu token en un objeto:

token <- "El token que INEGI te da al registrarte"

Luego, almacenamos en una lista los códigos de las series que vamos a usar:

series <- list(
  baja_california = "739297",
  coahuila = "739327",
  edomex = "739427",
  nuevo_leon = "739467", 
  puebla = "739487",
  qro = "739497" 
)

Posteriormente, usamos la función que creé y expliqué en mi ya citado blog:

# Definir función para procesar cada serie
procesar_serie <- function(serie_id, token) {
  serie <- inegi_series(serie_id, token)
  serie <- as_tibble(serie) %>%
    select(-date_shortcut, -notes) %>%
    arrange(date)
  return(serie)
}

Ahora vamos a bajar todas las series de una vez:

# Procesar cada serie
for (nombre_serie in names(series)) {
  assign(nombre_serie, procesar_serie(series[[nombre_serie]], token))
}

Lo que haré ahora será unir todos los dataframes descargados, como todos tienen el mismo número de observaciones y columnas, además de que las columnas están en el mismo orden, puedo unirlas sin miedo a que no tengan la misma estructura:

ind_man <- bind_rows(baja_california, coahuila,
                     edomex, nuevo_leon, puebla, qro)

#Aquí creo un vector con los nombres de cada dataframe
nombres <- c("baja_california", "coahuila",
             "edomex", "nuevo_leon", "puebla", "qro")

Si revisamos el código anterior, he creado un vector con los nombres de cada dataframe en el mismo orden que los uní con bind_rows(). Si vemos el resultado hasta ahora tenemos una tabla como sigue:

date values
2003-01-01 85.87970
2003-02-01 81.50035
2003-03-01 80.49787
2003-04-01 87.16143
2003-05-01 81.35903

Solo hay dos columnas, una de fechas y otra del valor de la serie, pero no hay forma de distinguir qué valores son de cada entidad. Para eso hice el vector nombres. Los valores de cada entidad vienen ordenados de acuerdo al orden que usamos en bind_rows(), así, los primeros valores son los de Baja California y los últimos son los de Querétaro. Cada entidad tiene 251 observaciones, esto lo podemos comprobar rápidamente así para cada entidad:

glimpse(qro)
## Rows: 251
## Columns: 2
## $ date   <date> 2003-01-01, 2003-02-01, 2003-03-01, 2003-04-01, 2003-05-01, 20…
## $ values <dbl> 71.88772, 71.57856, 67.15565, 73.00931, 65.86555, 64.84441, 64.…

Vamos a crear una variable categórica que indique la entidad federativa a la que pertenece cada observación (fila):

ind_man <- ind_man %>% 
  mutate(
    nombres = rep(nombres, each = 251)
  )
date values nombres
2003-01-01 85.87970 baja_california
2003-02-01 81.50035 baja_california
2003-03-01 80.49787 baja_california
2003-04-01 87.16143 baja_california
2003-05-01 81.35903 baja_california

Ya tenemos un dataframe de 3 columnas: la fecha, los valores de la serie y el Estado de la República al que pertenecen. El objetivo de este blog es desestacionalizar las 6 series al mismo tiempo, existen diferentes formas de hacerlo, yo optaré por hacerlo con la filosofía del tidy modeling. La filosofía del tidy modeling se resume en la idea de que los datos deben estar “diseñados para seres humanos”. Esto implica que los datos deben estar organizados de manera clara y comprensible, facilitando su manipulación, interpretación y comunicación. En lugar de enfrentarse a una maraña de datos confusos, el tidy modeling promueve la estructuración de datos de manera que se alinee con la forma en que los humanos piensan y trabajan, lo que mejora significativamente la eficiencia y la efectividad en el análisis de datos y modelado. Si te interesa este tema te recomiendo revisar los libros que cité al inicio, especialmente este blog se basa en los principios mostrados en Tidy Modeling with R de Max Kuhn y Julia Silge (click en este enlace).

Aclarado lo anterior, proseguimos:

ind_man_list <- ind_man %>%
  select(nombres, values) %>% 
  group_by(nombres) %>%
  group_nest()

En este fragmento de código estamos utilizando dos funciones del paquete tidyverse para organizar nuestros datos de manera efectiva. Primero, usamos la función group_by() para agrupar nuestros datos según la variable nombres. Esto significa que estamos creando grupos separados para cada valor único en la columna nombres. Luego, aplicamos la función group_nest(), la cual toma esos grupos y los encapsula en una estructura de datos especial llamada lista. Esta estructura denominada listas anidadas nos permite mantener los grupos de datos organizados de manera ordenada y fácilmente accesible para su posterior procesamiento. En resumen, al usar group_by() y group_nest(), estamos preparando nuestros datos de manera eficiente para aplicar técnicas de modelado de datos utilizando el enfoque tidy modeling.

glimpse(ind_man_list)
## Rows: 6
## Columns: 2
## $ nombres <chr> "baja_california", "coahuila", "edomex", "nuevo_leon", "puebla…
## $ data    <list<tibble[,1]>> [<tbl_df[251 x 1]>], [<tbl_df[251 x 1]>], [<tbl_d…

Como vemos, ahora tenemos un dataframe de 6 filas y 2 columnas (ya no es un dataframe de 3 columnas y 1506 filas), cada fila es uno de los estados de México; la primera columna tiene el nombre de cada estado y la segunda columna, data, es una columna donde cada observación es una lista que contiene las 251 observaciones del indicador de actividad industrial.

Con esto ya tenemos la estructura de datos necesaria para una modelación eficiente. Recordemos, buscamos desestacionalizar cada serie, es decir, modelar cada serie.

Ajuste estacional con seasonal

Como dije al inicio de este blog, la librería seasonal permite desestacionalizar las series de tiempo con la metodología desarrollado por el United States Census Bureau. No me detendré a explicar esta metodología pues no es mi propósito, así que si te interesa profundizar en cómo se desestacionaliza una serie, la matemática que hay detrás, cualquier manual universitario de Análisis de Series de tiempo es un buen inicio.

Ahora, volviendo al tema que nos ocupa, vamos a usar seasonal, así que recomiendo enfáticamente leer siempre la documentación de la librería que vamos a usar, sobretodo si es para modelación.

Aquí dejo dos links muy importantes que deberías leer si deseas saber más sobre esta librería:

  • Documentación (disponible aquí)
  • Artículo académico de los desarrolladores de seasonal donde con ejemplos prácticos muestran cómo funciona su paquete (disponible aquí)

Requisitos de seasonal

Como bien leíste en la documentación de seasonal ;) su principal función es seas(). Con solo poner el nombre de tu serie de tiempo la función buscará el ARIMA óptimo para desestacionalizar la serie. Ahí tenemos el primer reto, la función NECESITA un objeto de tipo serie de tiempo (ts). Un objeto ts se utiliza para almacenar datos secuenciales que están indexados por tiempo. Estos tipos de datos especifican la frecuencia y periodo de la serie. Bien sabemos que la frecuencia de nuestros datos es mensual y que inician en enero 2003 y terminan en noviembre 2023.

#como ejemplo mostraré solo los de Coahuila
coahuila %>% summarize(
  inicio = first(date),
  final = last(date)
) %>% knitr::kable()
inicio final
2003-01-01 2023-11-01

El primer reto es convertir a ts todas las series, las cuales almacenamos en una lista para cada entidad federativa.

Vamos a definir algunos parámetros de la función ts():

# Frecuencia
freq <- 12

# Fecha de inicio
fecha_inicio <- as.Date("2003-01-01")

Ahora, podemos convertir a serie de tiempo como sigue:

edomex_ts <- ts(data = edomex$values, start = c(2003, 1), frequency = freq)
edomex_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2003  82.06657  83.56723  85.51187  86.54361  89.98652  89.56908  88.90984
## 2004  81.86807  85.05927  90.98427  86.13722  87.75478  89.09542  86.25499
## 2005  84.00183  86.13434  86.66691  91.02925  90.30600  89.13341  85.27205
## 2006  87.46764  86.48714  94.38286  87.99085  96.18563  98.42538  90.16859
## 2007  88.64705  85.00947  93.81989  88.56174  92.70582  94.60240  90.35015
## 2008  83.84013  88.14217  86.02752  95.41724  90.37797  91.56486  90.20909
## 2009  76.57139  75.15352  81.32044  78.30549  76.36721  78.78806  82.80760
## 2010  76.66024  77.48581  87.74024  82.53508  84.46261  90.18658  89.41083
## 2011  81.71812  82.59004  94.99137  87.72493  92.53937  95.66055  91.30539
## 2012  89.99269  89.30402  98.52410  89.66291  96.10551  98.37682  94.75220
## 2013  90.77210  87.04797  87.99757  94.68793  94.44584  93.88732  94.91633
## 2014  89.30641  85.22599  95.40072  91.77638  97.03131  95.64972  91.68055
## 2015  89.21552  87.99743  95.44596  92.40799  94.10442  98.09590  98.20985
## 2016  90.35473  88.14896  89.75869  90.87920  92.87522 100.60682  93.34433
## 2017  96.06794  92.28546 101.56076  92.70588 100.13681 103.44637 100.49582
## 2018  97.41605  92.71646  99.99506  99.65229 105.07458 102.62981 102.15882
## 2019  94.19844  92.93225  97.16156  95.28422  96.84837  94.63141  98.52405
## 2020  90.31449  87.44889  90.85424  65.55554  64.52581  80.01410  88.45635
## 2021  87.61838  79.02689  90.99847  82.18677  83.47472  92.39254  91.19064
## 2022  84.63217  86.69966  99.78670  88.85563  92.82409  95.58056  96.95047
## 2023  87.34856  85.49378  93.56279  91.93383  99.88261  97.17991  93.72253
##            Aug       Sep       Oct       Nov       Dec
## 2003  86.06559  90.81075  91.35054  77.96056  80.47088
## 2004  86.00884  88.09036  90.12669  83.10003  83.87713
## 2005  90.76886  93.31606  90.85966  85.10485  81.86958
## 2006  95.11201  91.81196  94.07331  90.27026  84.52995
## 2007  92.92673  89.92507  94.62811  90.14051  78.92079
## 2008  91.42441  90.15705  93.86736  84.05587  76.67480
## 2009  81.76951  82.30821  88.37847  81.22423  78.12004
## 2010  90.24746  90.99329  92.32959  87.69883  83.60738
## 2011  95.41111  92.78656  96.33226  95.49337  88.62572
## 2012  97.17360  95.32026 100.21359  95.41444  86.77199
## 2013  96.12026  91.48687  98.30138  91.47415  86.04739
## 2014  89.58625  88.77279  94.38638  88.65575  86.74469
## 2015  95.36857  95.81981 101.14895  95.32455  90.07012
## 2016  98.23914  96.68794  95.70212  95.37733  91.33821
## 2017 107.40534 102.18246 110.22277 100.97312  90.82715
## 2018 107.65804  97.50920 105.64352 100.76134  88.78483
## 2019 100.86529  91.25960  93.50755  95.44079  83.80548
## 2020  77.81197  88.47291  95.92441  87.92218  89.25255
## 2021  96.50007  94.62744  96.61545  94.01701  97.99408
## 2022  99.26996  96.33652  97.35519  96.32795  94.58028
## 2023  99.27349  96.31190 101.47640  93.40722

Luego, podemos aplicar el ajuste estacional con seas() y después obtener la serie pronosticada (con ajuste estacional) mediante predict()

m <- seas(edomex_ts)
predict(m)
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2003  84.76141  87.14721  84.66384  86.69921  87.52454  87.05927  87.02160
## 2004  85.92221  87.02904  87.22456  86.33619  86.82344  85.20491  85.84341
## 2005  89.41756  90.01059  87.35216  88.24798  88.11472  85.20316  86.28551
## 2006  91.55559  90.73499  90.57331  91.23665  92.77542  94.35247  91.06474
## 2007  91.50959  89.67008  91.43810  90.53970  89.59102  91.95223  89.66932
## 2008  86.91025  90.05573  89.58137  91.68268  88.96844  89.04052  87.90931
## 2009  81.11967  80.54929  78.99873  79.38193  76.63865  74.98126  80.37454
## 2010  82.56378  83.17207  83.90933  83.87572  84.77866  86.28338  88.28789
## 2011  87.52519  88.54692  91.07537  90.72614  91.42245  91.61884  91.60670
## 2012  94.18821  92.30864  96.02477  92.87479  93.55831  95.66570  93.63384
## 2013  93.35351  93.27649  91.47425  92.13640  91.85632  92.57039  92.34914
## 2014  91.69766  91.49828  94.45442  93.80569  95.82773  92.91674  89.03362
## 2015  92.84825  94.35540  93.26395  94.69039  94.51343  94.03029  95.38918
## 2016  95.28095  91.50772  90.86533  90.41268  92.15488  96.74860  93.22910
## 2017  99.55586  98.83688  98.25904  98.31231  98.20724  99.92415 100.22860
## 2018  99.61685  99.35443 102.70059  99.41835 103.34218 100.94117 100.31157
## 2019  96.45766  99.62171  96.95636  98.14914  95.40346  94.80381  95.11697
## 2020  92.57837  92.59670  89.41360  68.55009  66.25476  77.71153  85.01267
## 2021  92.83700  86.07356  88.21443  85.24002  85.28448  90.20165  89.15370
## 2022  90.18831  93.89922  97.06166  93.17374  93.00579  93.36662  96.32280
## 2023  91.73676  92.82349  91.00919  97.50718  98.40123  94.93650  93.12341
##            Aug       Sep       Oct       Nov       Dec
## 2003  85.90250  87.79455  85.66760  83.17271  85.60270
## 2004  84.36999  85.14344  87.25921  85.21791  88.97645
## 2005  87.52707  90.43312  87.95550  86.83141  88.44724
## 2006  91.61508  90.44673  89.63556  91.45045  92.57787
## 2007  89.15247  90.06758  88.53345  90.76259  87.01438
## 2008  90.22993  87.51715  87.46700  87.03192  81.81240
## 2009  80.48254  79.83516  83.25121  82.45697  83.06193
## 2010  87.51280  88.70976  88.64438  87.29341  88.39503
## 2011  91.26707  90.70979  92.73141  94.85811  94.67866
## 2012  92.97426  96.27218  93.84561  94.62428  94.16123
## 2013  93.24968  91.24250  92.02183  92.02619  91.96959
## 2014  88.01031  87.27110  88.13304  90.55047  91.19903
## 2015  93.47609  94.23073  96.18245  95.64867  94.43746
## 2016  93.11432  95.00200  92.04677  94.19975  97.09184
## 2017 101.84555 101.82696 104.98393  99.64295  97.96661
## 2018 101.75383  98.55062  98.95194  99.20300  95.69057
## 2019  96.18842  90.77669  86.74012  95.12794  88.82435
## 2020  74.53267  86.35764  90.32039  87.55444  92.29448
## 2021  91.87569  92.37224  92.31891  92.35228 100.70612
## 2022  93.32874  93.99858  93.01563  94.80083  98.51678
## 2023  93.39881  95.32431  95.63763  92.00022

El objeto automáticamente reconoce a qué año y mes corresponde cada valor, sin embargo, hacer esto cada vez para cada serie y luego desestacionalizarlos no suena muy eficiente. Vamos a convertir todos los elementos en ts.

Hemos definido una función para dicha tarea:

# Función para convertir un tibble en un objeto ts
convertir_a_ts <- function(tibble_data) {
  ts_data <- ts(tibble_data$values, frequency = freq, start = c(year(fecha_inicio), month(fecha_inicio)))
  return(ts_data)
}

seas() y map(), una aplicación del Tidy modeling

Si leemos la documentación de seasonal, podremos ver que la librería permite desestacionalizar mútiples series a la vez, la forma en que lo hacen es con lapply(), esta es una función que se utiliza para aplicar una función a cada elemento de una lista o un vector, y luego devuelve una lista con los resultados de aplicar la función a cada elemento individual. Lo que tendríamos que hacer es tener todas nuestras series de tiempo (en formato ts), almacenarlos en una lista y luego aplicar seas() mediante lapply(), esto puede lucir más directo y sencillo, y puede serlo, ¿pero qué pasaría si fueran más de 100 series, las cuales en origen no están en formato ts?

Recordemos que convertimos nuestro dataframe en uno agrupado por entidad federativa donde sus respectivas series fueron anidadas en una lista, esto tuvo un motivo, el cual es usar la función map(). map() es una función clave del paquete purrr en R que se utiliza para aplicar una función a cada elemento de una lista, vector u otro tipo de objeto iterativo, y luego devuelve una lista con los resultados de aplicar la función a cada elemento individual. Ahora, en el contexto del tidy modeling y el uso de map(), hay algunas consideraciones importantes:

  • Iteración sobre modelos: Una aplicación común de map() en tidy modeling es iterar sobre diferentes modelos. Por ejemplo, puedes tener una lista de modelos de regresión lineal, regresión logística, y así sucesivamente, y usar map() para ajustar cada modelo a tus datos. Esto permite una forma eficiente y flexible de ajustar múltiples modelos y comparar sus resultados.

  • Preprocesamiento de datos: Antes de ajustar modelos, a menudo es necesario realizar operaciones de preprocesamiento en los datos, como imputación de valores faltantes, estandarización de variables, o creación de nuevas características. map() puede ser útil para aplicar estas transformaciones a diferentes conjuntos de datos de manera consistente.

Estas dos consideraciones que menciono son el objetivo del blog; primero preprocesamos nuestros datos para convertirlos en series de tiempo, luego vamos a iterar el modelo del ajuste estacional a nuestras series.

Enseguida me centraré en cómo usar la función que definimos.

# Convertir cada elemento de la lista de tibbles en un objeto ts
ind_man_list <- ind_man_list %>%
  mutate(data = map(data, convertir_a_ts))
  • mutate(data = map(data, convertir_a_ts)): Aquí es donde se utiliza map() en conjunto con mutate(), una función de dplyr que permite crear o modificar columnas en un dataframe.
    • data: Se está aplicando map() a la columna data de ind_man_list. Esta columna es la que tiene la lista de datos.
    • convertir_a_ts: Esta es una función que se está aplicando a cada elemento de la columna data. convertir_a_ts es la función definida por nosotros que convierte los datos a objetos de series de tiempo (ts).
    • mutate(data = …): La función mutate() está reemplazando la columna data original con los resultados de aplicar map(convertir_a_ts) a cada elemento de esa columna.

En resumen, cada conjunto de datos dentro de la lista data se convertirá en un objeto de series de tiempo (ts), lo que será de utilidad para usar seasonal.

Ahora, ya con cada elemento de la lista convertida a ts podemos continuar, del ejemplo de la serie del Estado de México aprendimos dos cosas: necesitamos aplicar la función seas() al vector de serie de tiempo y luego usar predict() para obtener la serie ya ajustada, es decir, necesitamos usar 2 funciones para lograr nuestro objetivo. Vamos de nuevo a utilizar map():

ind_man_list <- ind_man_list %>%
  mutate(data_sa = map(data, ~ seas(.x))) %>% 
  mutate(data_sa_final = map(data_sa, ~ predict(.x)))

En el primer mutate() creamos una columna nueva llamada data_sa, la cual contendrá toda la información generada por la función seas(), esta columna tendrá en cada una de sus observaciones (los estados de la República) una lista con todos lo generado por la modelación, es decir, cuadros de regresión, coeficientes, errores del modelo, serie ajustada, etc. Luego a partir de esta columna data_sa vamos a obtener la serie desestacionalizada con predict(), tal como lo hicimos con el ejemplo de Estado México, pero ahora lo estamos haciendo simultáneamente con todos los estados almacenados en el tibble ind_man_list. El resultado final es un dataframe o tibble de columnas:

glimpse(ind_man_list)
## Rows: 6
## Columns: 4
## $ nombres       <chr> "baja_california", "coahuila", "edomex", "nuevo_leon", "…
## $ data          <list> <85.87970, 81.50035, 80.49787, 87.16143, 81.35903, 86.2…
## $ data_sa       <list> [0.797547606, 5.505786863, 1.021021444, -0.601887069, -…
## $ data_sa_final <list> <87.47524, 86.79851, 84.79225, 85.87499, 85.45261, 85.8…

ind_man_list es un dataframe de 6 filas (por las 6 entidades federativas) y 4 columnas; la primera columna contiene los nombres de cada una, la segunda columna, data, contiene los valores de las variables en formato ts, data_sa contiene la información generada por la modelación de seas() y data_sa_final contiene los valores finales con el ajuste mediante predict. Cabe señalar que estas últimas 3 columnas son cada una listas que contienen los valores.

Lo que sigue es extraer los datos finales de dichas listas.

ind_man_list_final <- ind_man_list %>%
  select(-data_sa) %>%
  unnest_longer(-nombres) %>% 
  bind_cols(ind_man %>% select(date))

Este código tiene como objetivo finalizar la preparación de los datos en el contexto del tidy modeling, especialmente en la manipulación y organización de datos para su posterior análisis. Ahora, desglosemos cada paso:

  • select(-data_sa): La función select() selecciona columnas específicas del dataframe ind_man_list. Aquí, -data_sa indica que se deben seleccionar todas las columnas excepto data_sa, lo que significa que la columna data_sa será eliminada del dataframe. Esto porque por el momento no nos interesa la demás información de los modelos realizados.
  • unnest_longer(-nombres): La función unnest_longer() de la librería tidyr se utiliza para desanidar listas en columnas, expandiéndolas a través de filas. El argumento -nombres indica que se desaniden todas las columnas excepto la columna nombres.
  • bind_cols(ind_man %>% select(date)): La función bind_cols() se utiliza para combinar columnas de diferentes dataframes por columnas. En este caso, se están combinando las columnas del dataframe original ind_man con las columnas seleccionadas del dataframe resultante de unnest_longer(). Lo que hacemos es pegar la columna de las fechas del dataframe original con el dataframe resultante, que solo tenía la columna de los nombres de cada entidad y sus datos originales y desestacionalizados.
nombres data data_sa_final date
baja_california 85.87970 87.47524 2003-01-01
baja_california 81.50035 86.79851 2003-02-01
baja_california 80.49787 84.79225 2003-03-01
baja_california 87.16143 85.87499 2003-04-01
baja_california 81.35903 85.45261 2003-05-01
baja_california 86.28048 85.80168 2003-06-01
baja_california 90.11158 87.26540 2003-07-01
baja_california 88.27101 87.94751 2003-08-01
baja_california 96.37173 87.84376 2003-09-01
baja_california 94.44599 87.66900 2003-10-01
baja_california 82.71260 86.43384 2003-11-01
baja_california 88.23116 88.71133 2003-12-01
baja_california 82.83335 86.55605 2004-01-01
baja_california 80.13604 85.44866 2004-02-01
baja_california 85.73944 85.99582 2004-03-01

Visualización de los resultados

De manera rápida podemos ver la diferencia de nuestras series originales con las ajustadas:

Sin embargo, la visualización dependerá mucho de lo que queremos comunicar o de nuestra pregunta u objetivo de investigación. En la visualización anterior mostré la diferencia entre las series originales y las ajustadas mediante seasonal, ahora vamos a relajarnos y visualizar con más esmero. Mostraré la evolución a partir de 2019 del indicador ya ajustado por estacionalidad para los 6 estados seleccionados.

ind_man_list_final <- ind_man_list_final %>% 
  filter(date > as.Date("2018-12-01")) %>% 
  mutate(
    entidades = recode(nombres,
      "baja_california" = "Baja California",
      "coahuila" = "Coahuila",
      "edomex" = "Estado de Méx.",
      "nuevo_leon" = "Nuevo León",
      "puebla" = "Puebla",
      "qro" = "Querétaro"
    )
  )

Ya con las transformaciones pertinenes al dataframe, toca visualizar, para esto usaré la paleta de colores del metro de la Ciudad de México disponible en MetroBrewer (repositorio aquí).

plot <- ind_man_list_final %>% 
  ggplot(aes(x = date, y = data_sa_final,
             color = entidades)) +
  theme_da + #este es mi propio tema o diseño
  geom_line(linewidth = 1.3) +
  scale_color_manual(values=metro.brewer("Mexico", 6)) +
  scale_x_date(date_labels = "%Y",
               date_breaks = "1 year") + 
  labs(title = "Indicador mensual de actividad industrial",
       subtitle = "Total de manufacturas, 2018=100",
       caption = "Series ajustadas mediante X-13ARIMA-SEATS \n API, INEGI | @ecodiegoale")