11 may 2023

Descomposición de una serie temporal en R: componente transitorio vs permanente

La siguiente entrada utiliza la inflación mensual de la República Dominicana para obtener una descomposición histórica de la evolución de la inflación trimestral en un componente transitorio vs. componente subyacente. Primero se activan las librearas requeridas y se importa la data (dataset), que incluye la fecha mensual y el IPC con datos usando como base la canasta 2019-2020.

# Load the data
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
 
library(readxl)
library(forecast)
library(lubridate)
library(tidyverse)
 
dataset <- read_excel("data.xlsx")
 
# A tibble: 472 x 2
   fecha                 ipc
   <dttm>              <dbl>
 1 1984-01-01 00:00:00  1.38
 2 1984-02-01 00:00:00  1.42
 3 1984-03-01 00:00:00  1.44
 4 1984-04-01 00:00:00  1.46
 5 1984-05-01 00:00:00  1.48
 6 1984-06-01 00:00:00  1.54
 7 1984-07-01 00:00:00  1.56
 8 1984-08-01 00:00:00  1.57
 9 1984-09-01 00:00:00  1.64
10 1984-10-01 00:00:00  1.68
# ... with 462 more rows

Posteriormente, se agrega la serie trimestralmente usando el promedio del IPC. Se identifica el trimestre de cada mes (quarter), posteriormente se crea un agregado del promedio trimestral del IPC, y sobre este IPC trimestral, calculamos la inflación trimestral (q_inf). En una entrada anterior se explicó en detalle la agregación temporal y las transformaciones de series temporales

 q_dataset <- dataset |>
  mutate(quarter = zoo::as.yearqtr(fecha, "%Q")) |> # create a new variable for the quarter
  group_by(quarter)  |>  # group by the quarter variable
  summarise(q_ipc = mean(ipc)) |>
  mutate(q_inf = ((q_ipc/dplyr::lag(q_ipc))-1)*100) |>
  slice(-1)
 
# A tibble: 157 x 3
   quarter   q_ipc q_inf
   <yearqtr> <dbl> <dbl>
 1 1984 Q2    1.49  5.66
 2 1984 Q3    1.59  6.53
 3 1984 Q4    1.76 11.0
 4 1985 Q1    2.11 19.6
 5 1985 Q2    2.24  6.43
 6 1985 Q3    2.32  3.53
 7 1985 Q4    2.41  3.67
 8 1986 Q1    2.44  1.23
 9 1986 Q2    2.39 -2.01
10 1986 Q3    2.42  1.40
# ... with 147 more rows

Ahora, definimos una inflación trimestral como una serie temporal usando la función ts.

 ts_data <- ts(q_dataset$q_inf, start = c(1984,2), frequency = 4)
 
            Qtr1        Qtr2        Qtr3        Qtr4
1984              5.66392672  6.52527285 10.99767476
1985 19.64583811  6.43449238  3.53425010  3.66683217
1986  1.23093134 -2.00846916  1.39780793  4.35942001
1987  1.45783902  4.75590834  5.30846177  6.77347683

Ahora, usamos la función stl (ver referencia: https://otexts.com/fpp2/stl.html) para obtener una descomposición de la serie, recuperando los componentes estacionales, la tendencia y el componente aleatorio de las series.

# Decompose the time series using the STL function
decomp <- stl(ts_data, s.window = "periodic")
 
# Extract the seasonal, trend, and remainder components
seasonal <- decomp$time.series[, "seasonal"]
trend <- decomp$time.series[, "trend"]
remainder <- decomp$time.series[, "remainder"]

Ahora, agregamos los componentes (remainder) transitorios vs. permanentes (trend + seasonal).

q_dataset$permanent <- trend + seasonal
q_dataset$transient <- remainder
 
# A tibble: 157 x 5
   quarter   q_ipc q_inf permanent transient
   <yearqtr> <dbl> <dbl>     <dbl>     <dbl>
 1 1984 Q2    1.49  5.66    5.67    -0.00924
 2 1984 Q3    1.59  6.53    8.72    -2.19  
 3 1984 Q4    1.76 11.0    10.7      0.325 
 4 1985 Q1    2.11 19.6    11.4      8.23  
 5 1985 Q2    2.24  6.43    8.27    -1.83  
 6 1985 Q3    2.32  3.53    5.81    -2.28  
 7 1985 Q4    2.41  3.67    3.29     0.376 
 8 1986 Q1    2.44  1.23    1.57    -0.340 
 9 1986 Q2    2.39 -2.01   -0.0908  -1.92  
10 1986 Q3    2.42  1.40    1.39     0.00589
# ... with 147 more rows

Finalmente, se crea un gráfico apilado asumiendo la incidencia de cada componente.  

q_dataset |>
  dplyr::filter(quarter >= "2000 Q1") |>
  gather(id, value, -c(quarter,q_ipc,q_inf)) |>
  ggplot(aes(x = quarter, y = value, fill = id)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("#56B4E9", "#E69F00"), name = element_blank(), labels = c("Permantente", "Transitorio")) +
  geom_line(aes(x = quarter, y = q_inf), size=0.8) +
  theme_classic() +
  theme(legend.position = "bottom")

Sin embargo, vemos que la metodología anterior asume los componentes estacionales como transitorios. Por lo que, la política monetaria estaría reaccionando a elementos estacionales al considerarlos como elementos permanentes, sin embargo, entendemos que este elemento debe agregarse como elemento transitorio.

20 abr 2023

Trabajar encuestas de hogares con R (usando dplyr)

Para ilustrar los siguientes ejemplos se utiliza las siguientes bases de miembros simuladas. Los códigos de simulación se anexan al final de la presente entrada. Son dos bases ejemplos: una con información sobre los miembros y otra de las viviendas. Como siempre, el primer paso es estudiar nuestra base de datos.

data_miembros
# A tibble: 14 x 6
   hogar_id  edad sexo  zona   empleo ing_lab_h
      <int> <int> <chr> <chr>   <int>     <dbl>
 1        1    33 M     Rural       1     286.
 2        1    34 F     Rural       1     199.
 3        1    34 M     Rural       0       0 
 4        1    17 M     Rural      NA       0 
 5        1    35 M     Rural       0       0 
 6        1    35 M     Rural       1      87.8
 7        2    13 F     Urbana     NA       0 
 8        2     7 F     Urbana     NA       0 
 9        2    33 F     Urbana      1      94.2
10        3    29 F     Urbana      1     257.
11        3    17 M     Urbana     NA       0 
12        4    18 M     Rural       1     477.
13        5    12 M     Urbana     NA       0 
14        5    22 F     Urbana      0       0 
 
data_vivienda
  hogar_id  pared estufa habitaciones
1        1  Block      0            3
2        2 Madera      0            1
3        3  Block      1            3
4        4  Block      1            2
5        5  Block      1            1

1.       Agregando información a nivel de miembros

Podemos contabilizar la cantidad de miembros por cada hogar combinando las funciones group_by+summarise. Summarise hace un collapse de la base de datos a nivel de los miembros identificados.

data_miembros |>
  group_by(hogar_id) |>
  summarise(n()) |>
  dplyr::select(`n()`)
# A tibble: 5 x 1
  `n()`
  <int>
1     6
2     3
3     2
4     1
5     2

Note que al final se usa la función n() para contabilizar la cantidad de miembros y solo aparece esta variable en la data final. En caso que quisiéramos poner esta variable delante de la base de miembros, usaríamos mutate.

data_miembros |>
  group_by(hogar_id) |>
  mutate(miembro = 1:n())

# A tibble: 14 x 7
# Groups:   hogar_id [5]
   hogar_id  edad sexo  zona   empleo ing_lab_h miembro
      <int> <int> <chr> <chr>   <int>     <dbl>   <int>
 1        1    33 M     Rural       1     286.        1
 2        1    34 F     Rural       1     199.        2
 3        1    34 M     Rural       0       0         3
 4        1    17 M     Rural      NA       0         4
 5        1    35 M     Rural       0       0         5
 6        1    35 M     Rural       1      87.8       6
 7        2    13 F     Urbana     NA       0         1
 8        2     7 F     Urbana     NA       0         2
 9        2    33 F     Urbana      1      94.2       3
10        3    29 F     Urbana      1     257.        1
11        3    17 M     Urbana     NA       0         2
12        4    18 M     Rural       1     477.        1
13        5    12 M     Urbana     NA       0         1
14        5    22 F     Urbana      0       0         2

Note que podemos contabilizar otras características de los miembros. Podemos contar el número de miembros y trabajadores de cada hogar, sumando las variables de empleo. Al ser binaria (0,1), la suma de esta indica la cantidad de observaciones que cumplen dichas características. Esto se cumple con cualquier condición lógica, siempre que sumemos estos vectores, obtendremos la cantidad de elementos que cumplen dichas condiciones. Finalmente, podemos verificar que podemos crear nuevas variables usando las variables previamente creadas. Puntualmente, se crea la variable percent_trab que indica que porcentaje de los miembros de cada hogar esta trabajando.

data_miembros |>
  group_by(hogar_id) |>
  summarise(miembro = max(1:n()),
            trabajo_hogar = sum(as.numeric(empleo), na.rm=TRUE),
            porcent_trab= trabajo_hogar/miembro)
# A tibble: 5 x 4
  hogar_id miembro trabajo_hogar porcent_trab
     <int>   <int>         <dbl>        <dbl>
1        1       6             3        0.5 
2        2       3             1        0.333
3        3       2             1        0.5 
4        4       1             1        1   
5        5       2             0        0 

Adicionalmente, podemos agregar la edad promedio y la edad mínima de cada hogar.

data_miembros |>
  group_by(hogar_id) |>
  mutate(edad_m = mean(edad),
         edad_min = min(edad))
# A tibble: 14 x 8
# Groups:   hogar_id [5]
   hogar_id  edad sexo  zona   empleo ing_lab_h edad_m edad_min
      <int> <int> <chr> <chr>   <int>     <dbl>  <dbl>    <int>
 1        1    33 M     Rural       1     286.    31.3       17
 2        1    34 F     Rural       1     199.    31.3       17
 3        1    34 M     Rural       0       0     31.3       17
 4        1    17 M     Rural      NA       0     31.3       17
 5        1    35 M     Rural       0       0     31.3       17
 6        1    35 M     Rural       1      87.8   31.3       17
 7        2    13 F     Urbana     NA       0     17.7        7
 8        2     7 F     Urbana     NA       0     17.7        7
 9        2    33 F     Urbana      1      94.2   17.7        7
10        3    29 F     Urbana      1     257.    23         17
11        3    17 M     Urbana     NA       0     23         17
12        4    18 M     Rural       1     477.    18         18
13        5    12 M     Urbana     NA       0     17         12
14        5    22 F     Urbana      0       0     17         12

2.       Estimar la incidencia de la pobreza monetaria

Ahora se muestra un ejemplo de como estimar la pobreza monetaria. Asumiendo una línea de pobreza de 5,000RD$, como tenemos el ingreso por hora, necesitamos estimar los salarios el ingreso mensual ing_t_h y luego dividiendo ese ingreso entre el total de miembros, obtenemos el ingreso per cápita del hogar para finalmente obtener una variable binaria que verifica si este ingreso per cápita está por encima o por debajo de la línea de pobreza.

data_pobre <- data_miembros |>
  group_by(hogar_id) |>
  mutate(miembro = max(1:n()),
   ing_t_h = sum(ing_lab_h)*8*23.8,
  ing_pc_h = ing_t_h/max(1:n()),
  pobres = (ing_pc_h < 5000)*1)
 
data_pobre
# A tibble: 12 x 9
# Groups:   hogar_id [4]
   hogar_id  edad sexo  empleo ing_lab_h miembro ing_t_h ing_pc_h pobres
      <int> <int> <chr>  <int>     <dbl>   <int>   <dbl>    <dbl>  <dbl>
 1        1    30 M          1      428.       2  81554.   40777.      0
 2        1    18 M          0        0        2  81554.   40777.      0
 3        2    16 F         NA        0        6 206198.   34366.      0
 4        2    27 F          1      409.       6 206198.   34366.      0
 5        2    26 M          1      452.       6 206198.   34366.      0
 6        2    10 F         NA        0        6 206198.   34366.      0
 7        2     8 M         NA        0        6 206198.   34366.      0
 8        2    27 M          1      222.       6 206198.   34366.      0
 9        3    10 F         NA        0        3      0        0       1
10        3     6 F         NA        0        3      0        0       1
11        3    11 M         NA        0        3      0        0       1
12        4    11 F         NA        0        1      0        0       1

Ahora, como tenemos la variable pobre, podemos obtener la incidencia de la pobreza y la cantidad de personas pobres.

data_pobre |>
  ungroup() |>
  summarise(num_pob = sum(pobres),
            por_pobr = mean(pobres))
# A tibble: 1 x 2
  num_pob por_pobr
    <dbl>    <dbl>
1       4    0.333

Ahora, podemos segmentar la según una categoría determinada. Esto lo hacemos en los casos donde necesitamos obtener incidencia de la pobreza por grupos específicos de mi población. Según los datos el 50% de las mujeres está en condición de pobreza, frente a un 17% en caso de los hombres.

data_pobre |>
  group_by(sexo) |>
  summarise(num_pob = sum(pobres),
            por_pobr = mean(pobres))
# A tibble: 2 x 3
  sexo  num_pob por_pobr
  <chr>   <dbl>    <dbl>
1 F           3    0.5 
2 M           1    0.167

3.       Combinando data

Finalmente, dado que tememos una basa de vivienda donde obtenemos información sobre las características de esta, podemos estar interesados en estimar que porcentaje de mi población vive en casa con paredes de zinc. En tal sentido necesitamos llevar esta información a la base de miembros, o por el contrario llevar el numero de miembros de cada vivienda a la base de viviendas. Ambas opciones pueden hacerse mediante un join usando el id del hogar como llave. Esto es especialmente útil cuando necesitamos combinar información contenida en distintas bases de datos.

Alt. 1. Llevo la información de la data de vivienda a la data de miembros y luego contabilizo mis datos.

left_join(data_1, data_2, by="id")

Noten como opera la función, vamos a pesar la información referida en la data_2, sobre la data_1 (que esta a la izquierda), usando el id como identificador. Por ejemplo, para el id 1 le vamos a colocar el material de la vivienda.

left_join(data_miembros, data_vivienda, by="hogar_id")
# A tibble: 14 x 7
   hogar_id  edad sexo  zona   empleo ing_lab_h pared
      <int> <int> <chr> <chr>   <int>     <dbl> <chr>
 1        1    33 M     Rural       1     286.  Block
 2        1    34 F     Rural       1     199.  Block
 3        1    34 M     Rural       0       0   Block
 4        1    17 M     Rural      NA       0   Block
 5        1    35 M     Rural       0       0   Block
 6        1    35 M     Rural       1      87.8 Block
 7        2    13 F     Urbana     NA       0   Madera
 8        2     7 F     Urbana     NA       0   Madera
 9        2    33 F     Urbana      1      94.2 Madera
10        3    29 F     Urbana      1     257.  Block
11        3    17 M     Urbana     NA       0   Block
12        4    18 M     Rural       1     477.  Block
13        5    12 M     Urbana     NA       0   Block
14        5    22 F     Urbana      0       0   Block

Ahora, como tememos la información de la pared en la base de miembros, podemos obtener el porcentaje de individuos que viven en casas con un determinado tipo de material (si quisiéramos saber el porcentaje de vivienda, repetimos este proceso en la base de viviendas). [Nota: tenga pendiente que ahora con esta base podemos obtener cualquier información adicional sobre el perfil de la población según el material de la vivienda].

left_join(data_miembros, data_vivienda, by="hogar_id") |>
  group_by(pared) |>
  summarise(num_v = n()) |>
  mutate(percent = num_v/sum(num_v))
 
# A tibble: 2 x 3
  pared  num_v percent
  <chr>  <int>   <dbl>
1 Block     11   0.786
2 Madera     3   0.214

Alt. 2. Llevo información de la data de miembros a la data de vivienda. Por tanto, debo antes de, agregar la data de miembros para poder pasar la data. Primero hacemos esta agregación y obtenemos el numero de miembros por hogar.

data_miembros |>
    group_by(hogar_id) |>
    summarise(n())
# A tibble: 5 x 2
  hogar_id `n()`
     <int> <int>
1        1     6
2        2     3
3        3     2
4        4     1
5        5     2

Luego anexamos estas cantidades en al base de viviendas usando el mismo left join.

left_join(data_vivienda,
          data_miembros |> group_by(hogar_id) |> summarise(n()),
          by="hogar_id")
 
  hogar_id  pared n()
1        1  Block   6
2        2 Madera   3
3        3  Block   2
4        4  Block   1
5        5  Block   2

Ahora con esta base podemos calcular los porcentajes anteriores, solo debemos sumar las frecuencias de los miembros del hogar, contenida en la data anterior.

left_join(data_vivienda,
          data_miembros |> group_by(hogar_id) |> summarise(n()),
          by="hogar_id") |>
  group_by(pared) |>
  summarise(num_v = sum(`n()`))
# A tibble: 2 x 2
  pared  num_v
  <chr>  <int>
1 Block     11
2 Madera     3

4.       Porcentaje de personas por sexo según rango de edad 

De las herramientas vistas hasta aquí, primero se recodifica la variable edad por rangos usando la función cut. Posteriormente es cuestión de contabilizar categorías y obtener las frecuencias relativas.

 left_join(data_miembros, data_vivienda, by="hogar_id") |>
    mutate(r_edad = cut(edad, seq(0,100,20))) |>
    group_by(r_edad, sexo) |>
    summarise(num_v = n()) |>
    mutate(percent = num_v/sum(num_v))

# A tibble: 4 x 4
# Groups:   r_edad [2]
  r_edad  sexo  num_v percent
  <fct>   <chr> <int>   <dbl>
1 (0,20]  F         2   0.333
2 (0,20]  M         4   0.667
3 (20,40] F         4   0.5 
4 (20,40] M         4   0.5

Ahora se agrega la opción pivot_inder que se utiliza para convertir una data de formato long a wide. Ahora se lleva cada sexo a una columna para obtener mejor representación de la información.

left_join(data_miembros, data_vivienda, by="hogar_id") |>
    mutate(r_edad = cut(edad, seq(0,100,10))) |>
    group_by(r_edad, sexo) |>
    summarise(num_v = n()) |>
    mutate(percent = num_v/sum(num_v)) |>
    dplyr::select(-num_v) |>
    pivot_wider(names_from = sexo, values_from = percent)
 
# A tibble: 4 x 3
# Groups:   r_edad [4]
  r_edad      F      M
  <fct>   <dbl>  <dbl>
1 (0,10]  1     NA   
2 (10,20] 0.2    0.8 
3 (20,30] 1     NA   
4 (30,40] 0.333  0.667

5.       Data simulada

# simulando datos de miembros
 
set.seed(20)
library(tidyverse)
#simula data 
 
# Define number of households and members per household
num_hogares <- 5
max_miembros <- 6
 
num_miembro_hogar <- sample(1:max_miembros, num_hogares, replace = TRUE)
hogar_id <- rep(1:num_hogares, num_miembro_hogar)
 
data_miembros <- data.frame(hogar_id=hogar_id)
 
n<-nrow(data_miembros)
data_miembros$edad <- sample(5:35, n, replace = TRUE)
data_miembros$sexo <- sample(c("M", "F"), n, replace = TRUE)
data_miembros$zona <- ifelse(hogar_id %in% sample(1:num_hogares, round(num_hogares*0.3)), "Rural", "Urbana")
 
# Simulate employment status for each adult
 
data_miembros <- data_miembros |>
  mutate(
    ramdon_var = runif(n, min = 0, max = 1),
    empleo = case_when(
      edad >= 18 & ramdon_var>0.4~ "1",
      edad >= 18 & ramdon_var<0.4~ "0",
      TRUE ~ NA_character_
    ),
    empleo = as.integer(empleo),
    ing_lab_h = case_when(
      empleo == 1 ~ runif(n, min = 85, max = 500),
      empleo == 0 ~  0,
      TRUE ~ 0
    )
  ) |>
  as.tibble() |>
  dplyr::select(-ramdon_var)
 
# simulando datos de viviendas
 
data_vivienda <- data.frame(hogar_id=1:num_hogares)
data_vivienda$pared <- sample(c("Madera", "Concreto", "Block"), num_hogares, replace = TRUE)
data_vivienda$estufa <-  sample(c(0, 1), num_hogares, replace = TRUE)
data_vivienda$habitaciones <-  sample(c(1:3), num_hogares, replace = TRUE)

Descomposición de una serie temporal en R: componente transitorio vs permanente

La siguiente entrada utiliza la inflación mensual de la República Dominicana para obtener una descomposición histórica de la evolución de la...