Tidy Modeling y la API de la Reserva Federal

Con aplicación de {purrr} y {seasonal}

Presentación

Esta publicación tiene dos objetivos, uno personal y otro general. El objetivo personal es darle continuación a mi entrada anterior, donde por medio del uso de la API de INEGI con {inegiR} escribí diferentes funciones para poder desestacionalizar múltiples series de tiempo a la vez con la librería {seasonal}, la cual es una interfaz para usar X-13-ARIMA-SEATS, el software de ajuste estacional de la Oficina del Censo de EE. UU., todo bajo el paradigma del Tidy Modeling. Mi estudio y práctica de este enfoque de análisis de datos con R continúa y he decidido mejorar dicho blog, sin embargo, considero que mejor que actualizarlo es crear esta continuación, la cual usará la API de la Reserva Federal de Estados Unidos, esto con la librería {fredr}. Esta entrada además ofrece una alternativa a la graficación automatizada; en esta entrada lo hice por medio de una función personalizada, ahora ofrezco cómo hacerlo con {purrr} por medio de la manipulación de datos contenidos en listas en lugar de dataframes.

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 la Reserva Federal de Estados Unidos de una manera muy sencilla utilizando principalmente la librería de tidyverse. No me detendré a explicar mucho sobre la manipulación de datos con esta librería, de esto tengo varias entradas, además que puedes hacer uso de herramientas de Inteligencia Artificial para que te ayude con eso. Lo importante en este blog será la manipulación de los datos por medio de listas y no las operaciones sobre los datos en sí mismas.

Otro objetivo de este blog es desestacionalizar múltiples series al mismo tiempo usando {seasonal} y {fredr}. El motivo de emplear específicamente la primera librería es que es una interfaz amigable para usar el 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 un gráfico de waffle (click en este para un tutorial completo 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(tidyverse)
library(seasonal)
library(scales)
library(ggh4x)

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

Cómo comenzar a usar la API de la Reserva Federal

Paso 1: instalación y token

Lo primero será instalar y llamar a la librería {fredr}.

#install.packages("fredr")
library(fredr)

Recomiendo dar una revisada a la documentación de la librería. Puedes revisarla en este enlace.

Después entraremos al siguiente enlace para poder registrarnos y así generar un token con el que podremos tener acceso a las series que nos interesan.

Accediendo al enlace, veremos la página de la FRED de la siguiente manera y hacemos click donde está señalado.

Una vez que te registres con tu correo, te vas a las opciones de tu cuenta (My Account) y das click en el botón de API Keys, ulteriormente aparecerá Request API Key:

Finalmente, después de llenar unos datos podremos solicitar el token.

Después, copiamos y almacenamos dicho token en un objeto como sigue y usamos la función de fredr_set_key() :

#key <- "TU TOKEN"
fredr_set_key(key)

Paso 2: Selección de variables

Hasta el momento no tengo una manera más eficiente de seleccionar los identificadores de las series de interés, pero comparto cómo lo hago yo.

Accediendo al sitio web de la FED ponemos en el buscador las series que nos interesan. Para este blog serán los establecimientos para todas las industrias, específicamente en el condado de Austin en Texas.

Después del menú que se despliega tras la búsqueda elegimos la serie de interés. Como nos interesan todos los establecimientos de Austin, vamos a seleccionar las cinco entidades que la conforman (Travis, Caldwell, Hays, Bastrop y Williamson). En este ejemplo solo muestro la serie para Travis.

Lo que tenemos que hacer es copiar esa cadena números y letras, pues es el identificador de esta serie.

Una vez hecho lo anterior para cada serie, vamos a guardarlas en una lista:

# Number of Private Establishments for All Industries, Austin Texas

Austin <- list(Travis = "ENU4845320510",
               Caldwell = "ENU4805520510",
               Hays = "ENU4820920510",
               Bastrop = "ENU4802120510",
               Williamson = "ENU4849120510")

Para poder acceder a dicha serie, {fredr} cuenta con la función fredr(). Esta se utiliza así;

Travis <- fredr(
    series_id = "ENU4845320510",
    observation_start = as.Date("2008-01-01"), #Elijo este inicio
    observation_end = Sys.Date() #Para tener la más reciente
    )

glimpse(Travis)

## Rows: 63
## Columns: 5
## $ date           <date> 2008-01-01, 2008-04-01, 2008-07-01, 2008-10-01, 2009-0…
## $ series_id      <chr> "ENU4845320510", "ENU4845320510", "ENU4845320510", "ENU…
## $ value          <dbl> 28299, 28527, 28932, 29119, 28979, 29068, 29182, 29256,…
## $ realtime_start <date> 2024-06-05, 2024-06-05, 2024-06-05, 2024-06-05, 2024-0…
## $ realtime_end   <date> 2024-06-05, 2024-06-05, 2024-06-05, 2024-06-05, 2024-0…

Lo anterior es rápido y práctico si trabajamos pocas series, pero imagina el caso donde queremos usar muchas series, crear un objeto por cada serie no es muy eficiente. Para solucionar esto propongo la siguiente función:

# Funcion 

multiple_fred <- function(x) {
  fredr(
    series_id = x,
    observation_start = as.Date("2008-01-01"),
    observation_end = Sys.Date() #Para tener la más reciente
  )
}

La función solo necesita un argumento, el identificador de la serie.

Ahora, recordemos que tenemos una lista donde cada elemento es el identificador de cada serie. La lista tiene un nombre por cada elemento y luce así:

# Lista 

Austin

## $Travis
## [1] "ENU4845320510"
## 
## $Caldwell
## [1] "ENU4805520510"
## 
## $Hays
## [1] "ENU4820920510"
## 
## $Bastrop
## [1] "ENU4802120510"
## 
## $Williamson
## [1] "ENU4849120510"

Paso 3: Dataframe anidado

Lo que haré ahora será crear un tibble o dataframe anidado, lo anterior significa que una o más columnas de nuestro dataframe es una lista de dataframes.

# Obtener los datos para cada código FRED
Austin <- tibble(
  County = names(Austin),
  Data = map(Austin, multiple_fred)
  )

Vamos a desglosar el código anterior:

• tibble(…): crea un dataframe optimizado para el uso de {tidyverse}. Mi tibble será un dataframe con dos columnas, County y Data, y cinco filas. County tendrá los territorios que conforman Austin (Travis, Caldwell, Hays, Bastrop y Williamson). Data tendrá las series completas correspondientes a cada territorio.
• County = names(Austin): devuelve los nombres de los elementos de la lista Austin original, por lo que names() devuelve un vector con los nombres de los condados.
• Data = map(Austin, multiple_fred):

  • map(): esta es una función de {purrr}, la cual aplica una función a cada elemento de una lista y devuelve una lista del mismo tamaño. El uso de {purrr} requiere paciencia, pues existen muchas versiones de map() para cada necesidad, aquí dejo una cheat sheet.

  • multiple_fred(): es la función que creamos.

Austin 

## # A tibble: 5 × 2
##   County     Data             
##   <chr>      <named list>     
## 1 Travis     <tibble [63 × 5]>
## 2 Caldwell   <tibble [63 × 5]>
## 3 Hays       <tibble [63 × 5]>
## 4 Bastrop    <tibble [63 × 5]>
## 5 Williamson <tibble [63 × 5]>

Ahora, Austin es un tibble de dos columnas por cinco filas, la columna Data tendrá un dataframe en cada celda, cada dataframe será de 5 columnas y 63 filas, justo como el dataframe que hicimos de ejemplo para el territorio de Travis. Las columnas que nos interesan dentro de cada dataframe anidado en el tibble Austin son date y value.

Como expliqué en el blog Cómo desestacionalizar múltiples series a la vez, el uso de {seasonal} requiere que los valores de la serie estén en formato de serie de tiempo, ts(), y con los datos directos de la API, estos están en formato numérico. Debemos convertir cada serie a ts.

Aplicación de {seasonal}

Aquí propongo la creación de una función que se pueda usar directamente a toda la columna Data, así podremos convertir todas las series a ts sin importar el número de filas o condados. Por ejemplo, aquí solo usamos 5, pero imagina fueran 20, 40 o 100 series.

Paso 4: Convertir a ts

# Definir la función para convertir a series de tiempo
my_ts <-  function(df) {
  ts(df$value, start = c(2008, 1), frequency = 4)
}

Vamos a desglosar la función anónima que hice:

• ts(…): es una función de R para convertir vectores numéricos a series de tiempo. el primer argumento es el vector, luego está el inicio y la frecuencia de la serie (mensual, trimestral, anual…).

  • Esta función solo requiere un argumento, df, que es cualquier dataframe.

  • df$value: extrae la columna value del dataframe anónimo, es decir, para cualquier dataframe se extrae la columna llamada value. Recordemos que la columna que tiene las series bajadas de la API se llama así.

  • start = c(2008, 1): Nuestra serie inicia en el primer trimestre de 2008.

  • frequency = 4: Son series trimestrales.

Ahora aplicamos nuestra función al dataframe anidado, para esto volveremos a usar map(). map() tomará la lista-columna Data y le aplicará la función my_ts().

# Convertir a series de tiempo
Austin <- Austin %>% 
  mutate(
      value_ts = map(Data, my_ts)
  )
head(Austin)

## # A tibble: 5 × 3
##   County     Data              value_ts    
##   <chr>      <named list>      <named list>
## 1 Travis     <tibble [63 × 5]> <ts [63]>   
## 2 Caldwell   <tibble [63 × 5]> <ts [63]>   
## 3 Hays       <tibble [63 × 5]> <ts [63]>   
## 4 Bastrop    <tibble [63 × 5]> <ts [63]>   
## 5 Williamson <tibble [63 × 5]> <ts [63]>

Ahora nuestro dataframe anidado, Austin, tiene una nueva columna, donde nuevamente cada elemento es una lista, pero cada lista es un vector de 63 elementos en formato ts.

Paso 5: Ajuste estacional

Lo que sigue es desestacionalizar cada serie. Para esto usaremos {purrr} y {seasonal}. El ajuste estacional va a requerir las siguientes funciones:

• {seasonal}:

  • seas(): estimar el modelo que mejor captura el efecto estacional.

  • predict(): calcular la serie pronosticada por el modelo computado en seas().

• {purrr}:

  • compose(): crea una función compuesta. Vamos a combinar las dos funciones de seas() que describí anteriormente.

  • safely(): crea una versión segura de la función, así si hay algún error en los cómputos, la máquina se sigue de largo con el siguiente sin interrumpir el proceso. Donde haya habido algún error marcará NULL.

my_seas <- compose(predict, seas) #Se compone de derecha a izquierda
safe_seas <- safely(my_seas)

Ya con esto aplicamos la función a nuestro dataframe anidado. Vamos a crear una nueva columna, value_sa, que también será una columna-lista. Aplicaremos la función safe_seas() a la columna-lista value_ts que ya está en formato ts

# aplico seas para generar el arima
Austin <- Austin %>%
  mutate(
    value_sa = map(value_ts, ~ safe_seas(.x)))
Austin

## # A tibble: 5 × 4
##   County     Data              value_ts     value_sa        
##   <chr>      <named list>      <named list> <named list>    
## 1 Travis     <tibble [63 × 5]> <ts [63]>    <named list [2]>
## 2 Caldwell   <tibble [63 × 5]> <ts [63]>    <named list [2]>
## 3 Hays       <tibble [63 × 5]> <ts [63]>    <named list [2]>
## 4 Bastrop    <tibble [63 × 5]> <ts [63]>    <named list [2]>
## 5 Williamson <tibble [63 × 5]> <ts [63]>    <named list [2]>

Como vemos, value_sa es una columna-lista de 2 elementos. Esto porque al emplear una función segura crea 2 resultados, el resultado exitoso y el error, en el caso de haber.

Los elementos de la columna-lista value_sa para el primer elemento de nuestro dataframe anidado se puede ver así:

Austin$value_sa[[1]]

## $result
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2008 28284.07 28589.65 28890.05 29112.85
## 2009 28963.70 29131.84 29139.69 29249.83
## 2010 29416.45 29553.76 29705.87 29907.70
## 2011 30274.98 30449.73 30869.18 31180.44
## 2012 31630.26 31945.00 32163.30 32529.16
## 2013 32899.59 33149.64 33500.36 33879.88
## 2014 34355.83 34721.08 35136.98 35504.54
## 2015 35945.99 36546.07 37209.96 37540.11
## 2016 38003.91 38347.02 38892.53 39279.74
## 2017 39798.97 40123.90 40408.34 40911.39
## 2018 41219.23 41414.71 41596.61 42291.10
## 2019 42828.40 43249.72 43771.44 44102.72
## 2020 44622.48 44873.25 45634.72 46635.19
## 2021 47605.95 48953.16 49715.77 50422.40
## 2022 50971.22 51505.72 52071.33 51984.07
## 2023 51695.86 52073.96 51758.78         
## 
## $error
## NULL

Los 2 elementos de esta columna lista se llaman result y error. Nosotros queremos result, que es donde viene el valor pronosticado sin efecto estacional para cada serie.

Nuevamente usaremos map() para aplicar una función a una columna-lista. La función que usaré será pluck() de {purrr} que se utiliza para extraer un elemento de una lista o un subelemento de un dataframe en una lista. El elemento a extraer se llama result, como ya vimos.

Austin <- Austin %>% 
  mutate(value_sa = map(value_sa, pluck, "result"))
Austin

## # A tibble: 5 × 4
##   County     Data              value_ts     value_sa    
##   <chr>      <named list>      <named list> <named list>
## 1 Travis     <tibble [63 × 5]> <ts [63]>    <ts [63]>   
## 2 Caldwell   <tibble [63 × 5]> <ts [63]>    <ts [63]>   
## 3 Hays       <tibble [63 × 5]> <ts [63]>    <ts [63]>   
## 4 Bastrop    <tibble [63 × 5]> <ts [63]>    <ts [63]>   
## 5 Williamson <tibble [63 × 5]> <ts [63]>    <ts [63]>

Este código es útil para limpiar y transformar datos, simplificando la estructura de las listas dentro de la columna value_sa al extraer y mantener solo los resultados deseados. value_sa ya solo es una lista de un elemento que contiene la serie desestacionalizada.

Las siguientes líneas van a juntar todas las columnas-listas en una sola columna-lista. Para una explicación detallada copia y pega en ChatGPT y pídele que desglose.

Austin <- Austin %>%
  unnest(-County) %>% 
  group_by(County) %>% 
  nest()

Austin

## # A tibble: 5 × 2
## # Groups:   County [5]
##   County     data             
##   <chr>      <list>           
## 1 Travis     <tibble [63 × 7]>
## 2 Caldwell   <tibble [63 × 7]>
## 3 Hays       <tibble [63 × 7]>
## 4 Bastrop    <tibble [63 × 7]>
## 5 Williamson <tibble [63 × 7]>

Graficar con un dataframe-list

En esta parte vamos a graficar las series desestacionalizadas para cada condado, todo al mismo tiempo, la lógica será como la anterior, a partir de una columna-lista vamos a aplicar una función, en este caso ggplot(), y obtendremos una nueva columna-lista donde cada observación será un gráfico. Así el proceso de graficar será automático.

Este fragmento no lo explicaré. Solo es la creación de variables.

Austin <- Austin %>%
  mutate(
    data = map(data, ~mutate(.x, 
                             year = year(date),
                             month = month(date),
                             month = case_when(
                               month == 1 ~ "I",
                               month == 4 ~ "II",
                               month == 7 ~ "III",
                               month == 10 ~ "IV",
                             ))
               )
  ) 

Lo que sigue será crear una nueva columna-lista llamada plots. Un gráfico de ggplot() es en sí mismo una lista, por lo que podemos usar algún miembro de la familia map().

Ahora usaremos map2(). Esta función permite aplicar una función a dos listas o vectores en paralelo. Es decir, toma dos listas (o vectores) y aplica una función que acepta dos argumentos, devolviendo una nueva lista con los resultados. La diferencia con map(), es que map() solo acepta 1 lista.

Bien podríamos usar simplemente map(), donde la lista a utilizar sería data e igualmente tendríamos nuestros gráficos. Yo uso map2() porque también deseo automatizar el título, si usamos map() podemos graficar pero el texto del título sería estático. Con map2() utilizo la columna-lista data y la columna County, la primera tiene los datos, la segunda es solo el nombre del condado.

Austin <- Austin %>%
  mutate(
    plots = map2(data, County, ~{
      ggplot(.x %>% 
               filter(year>2018), 
             aes(x = interaction(month, year), group = 1)) +
        geom_bar(aes(y = value, color = "Not Seasonally Adjusted"),
                 stat="identity", fill="#FF9800", alpha = 0.5) +
        labs(title = "Number of Private Establishments for All Industries",
             subtitle = paste(.y, "County,", "Austin Texas"),
             y = "Seasonally Adjusted",
             caption = "API, FRED | @ecodiegoale")+
        scale_x_discrete(NULL, guide = "axis_nested") +
        scale_y_continuous(label = comma)+
        theme_da # mi propio tema
    })
  )
Austin

## # A tibble: 5 × 3
## # Groups:   County [5]
##   County     data              plots 
##   <chr>      <list>            <list>
## 1 Travis     <tibble [63 × 9]> <gg>  
## 2 Caldwell   <tibble [63 × 9]> <gg>  
## 3 Hays       <tibble [63 × 9]> <gg>  
## 4 Bastrop    <tibble [63 × 9]> <gg>  
## 5 Williamson <tibble [63 × 9]> <gg>

Vamos a desglosar el código anterior.

1. Hemos creado la columna-lista plots.

2. Para su creación usamos map2(x, y, ~{fn}). La cual necesita 2 elementos y una función a aplicar. En este caso la función es ggplot(), pero como esta función a su vez requiere los elementos de data y el nombre del condado para que el título sea automática, fungirá como una función anónima.

   1. El primer elemento es data (x), el cual es el nombre de la columna-lista que contiene los dataframes para cada condado. El segundo elemento es County (y), el cual solo es una columna con los nombres de cada condado.

   2. ~{ggplot(data = .x) + labs(subtitle = .y)}: la función ggplot() tomará los 5 dataframes de la columna lista data y sus respectivos nombres County y graficará, cada gráfico se almacenará como lista en la nueva columna plots.

Como se ve en el cuadro anterior, la columna-lista plots es tiene elementos gg, los cuales son las listas que almacenan todo lo que se ve en un gráfico. Sin embargo, una lista no nos sirve como visualización. Vamos a imprimir todos los gráficos.

Usaremos walk() de {purrr}, la cual es similar a map(), pero se utiliza cuando quieres aplicar una función a cada elemento de una lista o vector por sus efectos secundarios y no necesitas que se devuelva un resultado, es decir, no devuelve una lista de resultados. En cambio, devuelve el input original (invisiblemente), pero es utilizado principalmente por sus efectos secundarios (por ejemplo, imprimir, guardar en archivos, etc.).

walk(Austin$plots, print)

A continuación, una función para guardar todas las imágenes en una carpeta. No explicaré el código.

output_dir <- "values_sa" #nombre de la carpeta donde almacenaré los gráficos
if (!dir.exists(output_dir)) {
  dir.create(output_dir)
}
# Iterar a través de los gráficos y guardarlos en el directorio
for (i in seq_along(Austin$plots)) {
  county_name <- Austin$County[i]
  plot <- Austin$plots[[i]]
  filename <- file.path(output_dir, paste0(county_name, ".png"))
  ggsave(filename, plot = plot, width = 14, height = 10, units = "cm")
}

Un waffle

Para concluir, mostraré cómo desanidar el dataframe-list y con esto poder hacer un gráfico de waffle. No me detendré mucho con lo último, si deseas un tutorial completo visita The R Graph Gallery, al inicio de este blog puse el link.

unnest_Austin <- Austin %>% 
  ungroup() %>% 
  select(-plots) %>% 
  unnest(-County) %>% 
  select(-starts_with("realtime_"))

unnest_Austin <- unnest_Austin %>% 
  mutate(
    year = as.character(year),
    value_sa = round(value_sa, 0)
  )

Ahora el código del waffle:

library(waffle)

p <- unnest_Austin %>% 
  filter(year == "2022") %>% 
  ggplot(aes(fill = County, values = value_sa/100))+
  geom_waffle(color = "white", size = .25, n_rows = 15, flip = TRUE) +
  facet_wrap(~ month, nrow = 1,
             strip.position = "bottom")+
  scale_y_continuous(labels = function(x) x * 15,
                     expand = c(0,0))+
  MetBrewer::scale_fill_met_d("Hiroshige", direction = -1)+
  theme_da+
  theme(
    axis.text.x = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    legend.position = "right",
    strip.background = element_blank(),
    strip.placement = "outside",
    panel.spacing.x = unit(0, "cm")
  )+
  labs(title = "Number of Private Establishments for All Industries*",
       subtitle = "Austin, Texas - Quarters of 2022",
       caption = "*Each square of the waffle chart represents 100 establishments. Seasonally adjusted figures. \nAPI, FRED | @ecodiegoale")+
  coord_equal()

p