El Banco Central de la República Dominicana opera bajo un esquema de metas de inflación, en el que establece una meta explícita durante un horizonte de política determinado y se compromete a tomar decisiones que lleven a la consecución del objetivo. Esta estrategia genera certidumbre en los agentes económicos y ha resultado ser unos de los pilares de la estabilidad macroeconómica del país en los últimos años.
La adopción de este esquema de política por parte de las autoridades monetarias hace de la inflación una variable interesante, alrededor de la cual existe todo un sistema de proyecciones. Por esto quise utilizarla para explicar algunos temas que desde hace un tiempo he querido compartir y de paso aportar enfoques novedosos para el manejo de múltiples modelos con R.
El objetivo del ejercicio será generar proyecciones de la inflación utilizando dos enfoques de pronóstico:
- Proyección general
- Agregación de pronósticos
En el primer enfoque la idea es aplicar varias metodologías autoregresivas a la inflación general, terminando con una serie y múltiples modelos.
El segundo enfoque cosiste en la agregación ponderada de las proyecciones de los grupos de artículos que componen la canasta del índice, de manera que el resultado contenga múltiples series y múltiples modelos.
Esta publicación se enfoca principalmente en la estrategia de programación para la realización de la tarea propuesta. Los aspectos metodológicos de los modelos utilizados no serán muy detallados, aunque habrá referencias metodológicas de publicaciones anteriores en este espacio o en fuentes alternativas. Por otro lado, los aspectos relacionados a la evaluación y selección de modelos quedará pendiente para otra publicación por cuestiones de volumen.
Para concluir con la introducción, es importante aclarar que el resultado de este ejercicio, aunque presenta proyecciones válidas de la inflación, no es necesariamente congruente con las proyecciones que evalúa el Banco Central para el diseño de la política monetaria. Esto es simplemente un ejercicio para compartir un enfoque para manejar múltiples modelos en R. De querer acceder a las proyecciones de inflación del Banco Centra lo ideal es consultar el Informe de Política Monetaria y el Programa Monetario y Financiero que publica la institución.
paquetes
library(tidyverse)
library(zoo)
library(forecast)
library(broom)
library(sweep)
library(lubridate)Los datos
Para seguir cada paso es importante acceder a un set de datos que he construido con las informaciones que publica el Banco Central. Pueden acceder al objetos RDS desde el repositorio en github de este post.
La base de datos tiene 9 variables, entre las que se encuentran el nombre del rubro de la canasta, los diferentes niveles de agregación de la canasta, las ponderaciones, la fecha, y las variaciones mensuales e interanuales.
# Archivos con las series a nivel de artículo
ipc_articulos <- readRDS("data/ipc_articulos_long.RDS")
# Glimpse de la data
glimpse(ipc_articulos, 70)## Rows: 56,264
## Columns: 9
## Groups: nombre [493]
## $ nombre <chr> "Indice General", "Alimentos y Bebidas ~
## $ ponderador <dbl> 100.00000, 25.10051, 22.50438, 5.11303,~
## $ division <chr> "Indice general", "Grupo", "Subgrupo", ~
## $ codigo <chr> "0", "1", "1.01", "01.01.01", "01.01.01~
## $ mes <date> 2011-01-01, 2011-01-01, 2011-01-01, 20~
## $ indice <dbl> 101.24, 101.71, 101.64, 101.94, 106.88,~
## $ fecha_ym <yearmon> Jan 2011, Jan 2011, Jan 2011, Jan 2~
## $ variacion_mes <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,~
## $ variacion_interanual <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,~La primera transformación a realizar será filtrar las observaciones correspondientes al IPC general y las series de los diferentes grupos de bienes y servicios. Aunque lo ideal sería probar la agregación de pronósticos desde el nivel más bajo (artículos), hacerlo a nivel de grupos puede que tenga mejores ventajas pedagógicas, porque al ser 12 grupos se pueden visualizar fácilmente los resultados. Aprovechando esta manipulación, resulta conveniente crear un factor con el nombre de los grupos, de esta manera se logra que el ipc general quede de primero en todas las visualizaciones.
Para crear el factor del nombre se pueden utilizar las funciones forcats::fct_inorder() que crea factores y organiza los niveles en orden de aparición, conveniente aquí porque el ipc general aparece de primero. Por otro lado, también está forcats::fct_reorder() que recibe un factor y una variable de organización como argumento. De usar esta segunda opción la variable ponderador serviría para organizar los factores y lograr el objetivo.
ipc_grupos <- ipc_articulos %>%
ungroup() %>%
filter(
division == "Grupo" | nombre == "Indice General") %>%
mutate(
nombre = str_wrap(nombre, width = 25),
nombre = fct_inorder(nombre)
)Visualizando la inflación
Es válido detallar en este punto que la estrategia de forecasting a utilizar consistirá en proyectar las variaciones mensuales del ipc general y de los diferentes grupos y con estas recuperar el IPC implícito para cada período. La inclinación por este procedimiento obedece a la evidencia empírica de que esta forma conlleva a menores errores de pronóstico que, por ejemplo, proyectar las diferencias logarítmicas de los índices.
El gráfico siguiente muestra los índices para el IPC general y los distintos grupos de la canasta. Con facilidad se puede apreciar las tendencias de cada uno, entre las que se destaca Alimentos y Bebidas no Alcohólicas, Bebidas Alcohólicas y Tabaco, Educación y Salud por acumular el mayor incremento desde el año base. Por otro lado, prendas de vestir traza una tendencia negativa muy peculiar.
ipc_grupos %>%
ggplot(aes(x = mes, y = indice, color = nombre)) +
geom_line(show.legend = FALSE) +
facet_wrap(~nombre, scales = "free_y") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45)) +
scale_y_continuous(breaks = ) +
labs(
x = "",
y = "índice",
title = "IPC e índice de los grupos de bienes y servicios"
)
ipc_grupos %>%
ggplot(aes(x = mes, y = variacion_mes * 100, color = nombre)) +
geom_line(show.legend = FALSE) +
facet_wrap(~nombre, scales = "free_y") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45)) +
labs(
x = "",
y = "Variación mensual",
title = "Variación mensual del IPC y de los grupos de bienes y servicios de la canasta"
)
Generando las proyecciones
Para la realización de este ejercicio se tomó en consideración el enfoque utilizado por Hadley Wickham y Garrett Grolemund en el capítulo 25 del libro R for data science. Aunque ellos realizaron algo muy distinto a esto, en general en ese capítulo los autores muestran como aprovechar las ventajas de los dataframe anidados (capacidad de tener listas de objetos como columnas de un data frame) en combinación con las herramientas del paquete purrr para la realización de operaciones iterativas.
Así, los pasos a dar para proyectar la inflación son los siguientes
- Anidar el data frame agrupando por el nombre del grupo, para tener series “separadas”.
- Agregar una columna con la transformación del las series a objetos de serie de tiempo.
- Agregar una columna con los modelos que le serán aplicados a cada serie y estimarlos
- Realizar un forecast de 6 meses de cada rubro
- Recuperar el índice implícito dada las proyecciones de las variaciones mensuales
Paso 1
La función tidyr::nest() permite anidar un dataframe en base a grupos. Este tipo de objeto facilita la tarea de aplicar acciones a subgrupos de datos, mejorando la manipulación de los resultados y evitando la duplicidad de códigos.
Para ilustrar un poco el concepto, una forma rápida de realizar el ejercicio sería hacer todos los pasos para uno de los grupos del IPC y copiar el código y pegarlo tantas veces como grupos hay, o bien hacer un loop usando los nombres de los grupos y aplicar la serie de pasos a todos en una sola ejecución. Pero hacer esto conllevaría, en el caso uno a la creación de muchos objetos y, en el caso dos, a la creación de varias listas de objetos (Complicados de manipular).
El procedimiento usando un data frame anidado ayuda al manejo de los diferentes objetos que deben ser creados, al mantenerlos en una estructura de filas y columnas. En esta estructura cada fila corresponde a uno de los grupos y las columnas pueden contener básicamente cualquier cosas, preferiblemente cosas relacionadas al grupo en cuestión o bien resultados de usar la data del grupo.
by_grupos <- ipc_grupos %>%
# agrupando
group_by(nombre, ponderador) %>%
# excluyendo la primera observación de cada grupo
# como es una variación mensual, la primera obs
# es un NA
slice(-1) %>%
# anidando las series
nest()
head(by_grupos)## # A tibble: 6 x 3
## # Groups: nombre, ponderador [6]
## nombre ponderador data
## <fct> <dbl> <list>
## 1 "Indice General" 100. <tibble [103 x 7]>
## 2 "Alimentos y Bebidas No\nAlcohólicas" 25.1 <tibble [103 x 7]>
## 3 "Bebidas Alcohólicas y\nTabaco" 2.32 <tibble [103 x 7]>
## 4 "Prendas de Vestir y\nCalzado" 4.56 <tibble [103 x 7]>
## 5 "Vivienda" 11.6 <tibble [103 x 7]>
## 6 "Muebles y Artículos para\nel Hogar" 6.46 <tibble [103 x 7]>En el objeto by_grupo, la fila 1 de la columna data tiene la serie del índice general.
# desagregando la data de la primera fila
# para mostrar el contenido
by_grupos %>%
slice(1) %>%
unnest(data) %>%
head()## # A tibble: 6 x 9
## # Groups: nombre, ponderador [1]
## nombre ponderador division codigo mes indice fecha_ym variacion_mes
## <fct> <dbl> <chr> <chr> <date> <dbl> <yearmo> <dbl>
## 1 Indice Ge~ 100. Indice ~ 0 2011-02-01 102. Feb 2011 0.0121
## 2 Indice Ge~ 100. Indice ~ 0 2011-03-01 104. Mar 2011 0.0116
## 3 Indice Ge~ 100. Indice ~ 0 2011-04-01 105. Apr 2011 0.00868
## 4 Indice Ge~ 100. Indice ~ 0 2011-05-01 105. May 2011 0.00220
## 5 Indice Ge~ 100. Indice ~ 0 2011-06-01 106. Jun 2011 0.0102
## 6 Indice Ge~ 100. Indice ~ 0 2011-07-01 107. Jul 2011 0.00850
## # ... with 1 more variable: variacion_interanual <dbl>Paso 2
El segundo paso es crear objetos ts con dada uno de los elementos de la columna data. Para hacer esto se utiliza alguna función que itere sobre los elementos de data y devuelva una lista con las series de tiempo. Las funciones de purrr cumplirán dicha tarea.
Con el siguiente código se añade la columna ts que almacena un objeto ts para cada grupo.
by_grupos <- by_grupos %>%
mutate(
ts = map(
data,
~.x %>%
select(variacion_mes) %>%
ts(
frequency = 12,
start = c(2011, 02)
)
)
)Paso 3
En este paso la idea es asignarle un dataframe a cada grupo en el que haya una columna con las funciones o modelos que se le van a aplicar a la serie y otra columna con los parámetros para estimar los modelos. En el caso del ejercicio actual se aplicarán 3 modelos a cada serie, un Arima automático, un suavizado exponencial y un BATS y no se especificarán parámetros adicionales a la serie, para simplificar el ejercicio.
by_grupos <- by_grupos %>%
# de vez en cu
mutate(
params = map(
ts,
~list(
auto.arima = list(y = .x),
ets = list(y = .x),
bats = list(y = .x))
),
params = map(params, enframe, name = "f", value = "params")
)En este punto la función purrr::invoke_map() será de utilidad, ya que permite acceder a los argumentos de múltiples funciones contenidos en una estructura tabular como la creada en el chunk de códigos anterior y aplicarlos en un proceso iterativo. Para leer un poco más sobre esta función consideren el apartado 21.7.1 del libro citado.
by_grupos <- by_grupos %>%
mutate(
models = map(
params,
~.x %>% mutate(
fit = invoke_map(f, params)
)
)
)Paso 4
A diferencia de algunas funciones del paquete forecast los modelos ets, auto.arima y bats no generan el forecast inmediatamente, haciendo necesario un paso adicional para obtenerlos. Por esto se agregará la proyección de cada uno de los grupos de la inflación utilizando la función forecast::forecast(). De igual forma, se utiliza la función sweep::sw_sweep() para llevar el summary de cada forecast a un formato tidy. Esta última función se creó para extender las funciones del paquete broom para manejar los resultados de los modelos de series de tiempo.
by_grupo_fcast <- by_grupos %>%
mutate(
# Agregando el forecast de los 3 modelos
# estimados para cada serie
models = map(
models,
~.x %>%
mutate(fcast = map(fit, forecast, h = 6))
),
#Summary de cada objeto forecast
models = map(
models,
~.x %>%
mutate(
sweep = map(
fcast,
sw_sweep,
fitted = FALSE,
timetk_index = TRUE,
rename_index = "date"
)
)
),
# Como eran 3 modelos para cada serie, se puede
# combinar el resultado del sw_sweep de cada uno
# en un solo tibble
sweep = map(
models, ~.x %>%
select(f, sweep) %>%
unnest(sweep)
)
)Paso 5
Para generar los índices implícitos y calcular las variaciones interanuales se realiza una extensa manipulación de datos. Es recomendable comprender el output de cada objeto previo para entender los pasos siguientes. También es válido ir corriendo parcialmente cada bloque del pipe line e ir viendo el proceso por etapas.
A lo largo de los chucks de códigos que vienen a continuación coloqué una serie de comentarios que podrían ayudar a comprender cada paso.
by_grupo_last <-
## creando un data frame con la data utilizada para
## hacer el forecast (variaciones mensuales) y los índices
## observados en cada periodo
left_join(
# Forecast con las variaciones e intervalos de confianza
by_grupo_fcast %>%
select(nombre, sweep, ponderador) %>%
unnest(sweep) %>%
ungroup() %>%
select(date, everything()),
# Agregando la columna con el indice de cada articulo
by_grupo_fcast %>%
select(nombre, data) %>%
unnest(data) %>%
ungroup() %>%
select(date = fecha_ym, nombre, division, indice),
by = c("date", "nombre")
) %>%
## Para recuperar el ipc en base al forecast dada la variacion
## mensual. Como lo que se poryectó fueron las variaciones mensuales,
## no tenemos índices los periodos del forecast. El primer paso
## para obtenerlos el ponerle al lado el último índice de las series observadas
group_by(nombre, f) %>%
# Yo sé que aquí podía usar la función dplyr::fill(), pero
# ya lo había hecho de la forma difícil cuando lo recordé
mutate(
indice = ifelse(
is.na(indice),
last(indice[key == "actual"]),
indice)
) %>%
group_by(nombre, f, key) %>%
## recuperando el ipc implicito data la proyección de la variación mensual.
## Esta es una formula sencilla
mutate(
cum_var = variacion_mes + 1,
cum_var = cumprod(cum_var),
indice = ifelse(key == "forecast", indice * (cum_var), indice)
) %>%
# Calculando la variación interanual implicita
group_by(nombre, f) %>%
mutate(
indice_vi = (indice - lag(indice, 12))/lag(indice, 12)
) %>%
ungroup()En vista de que la función sw_sweep retorna la serie original y las observaciones proyectadas para cada modelo, lo ideal es no tener repetida las series originales tantas veces. El código a continuación elimina una poco esa redundancia.
by_grupo_last <-
by_grupo_last %>%
filter(key == "actual", f == "auto.arima") %>%
select(-f) %>%
bind_rows(
by_grupo_last %>%
filter(key != "actual")
) %>%
mutate(
date = zoo::as.Date.yearmon(date),
f = ifelse(is.na(f), "observado", f),
key = ifelse(key == "actual", "observado", key)
) %>%
select(-cum_var)Visualizando los resultados
Primero el resultado del pronóstico de la inflación general con las distintas metodologías, para ver el resultado se usará el gráfico de la variación interanual implícita.
Estas visualizaciones podrían representar buena guía para la realización de gráficos en el futuro, recomiendo prestar atención a los detalles.
by_grupo_last %>%
filter(nombre == "Indice General") %>%
# para que los valores observados tengan línea continua
# y que ggplot no lo organice en orden alfabético
mutate(key = factor(key, c("observado", "forecast"))) %>%
# No es necesario mostrar toda la serie observada
filter(date > "2017-12-01") %>%
ggplot(aes(x = date, y = indice_vi, color = f, linetype = key)) +
geom_line(size = 1) +
theme_minimal() +
# No creo que sea necesario incluir una leyenda para
# el tipo de línea
guides(linetype = "none") +
labs(
x = "",
y = "Variación interanual",
color = "Modelo",
title = "Proyección de la variación interanual del IPC"
) +
theme(legend.position = "bottom")
Ahora las variaciones mensuales proyectadas para cada grupo de la canasta. Este es otro gráfico instructivo ya que combina dos geometrías. geom_ribbon y geom_line.
by_grupo_last %>%
# restringiendo el output para una sola metodología
filter(f %in% c("observado", "auto.arima")) %>%
mutate(key = factor(key, c("observado", "forecast"))) %>%
filter(date > "2017-12-01") %>%
# esteticas del gráfico
ggplot(
aes(
x = date, y = variacion_mes * 100,
color = nombre, linetype = key)
) +
# Geometría de área con el intervalos de confianza
# d ela proyección
geom_ribbon(
aes(
ymin = lo.80 * 100, ymax = hi.80 * 100),
fill = "gray", alpha = 0.8, color = NA
) +
# Líneas
geom_line() +
# Facets
facet_wrap(~nombre, scales = "free_y") +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text.x = element_text(angle = 45)
) +
# Modificando el formato de la etiquetas
# del eje x. Recordar que el argumentos labels
# recibe o un vector con las etiquetas o una
# función con la transformación que se le debe hacer
scale_x_date(
labels = function(x){format(x, "%b-%y")}
) +
guides(
color = "none",
linetype = "none"
) +
labs(
x = "",
y = "Variación mensual",
title = "Proyección de la variación mensual por grupo"
)## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
Por último el forecast con las dos estrategias de pronóstico principales. En este caso, salta a la vista que la metodología de agregación de pronósticos es un poco más volátil, quizás porque recoge mejor las fluctuaciones estacionales de la serie.
# Uniendo los resultados de los diferentes modelos
# para los 2 enfoques de pronóstico
bind_rows(
# resultado del forecast de la serie general
by_grupo_last %>%
filter(
nombre == "Indice General") %>%
select(date, key, f, variacion_mes),
# Resultado de la agregación del pronóstico
# de los grupos
by_grupo_last %>%
filter(
nombre != "Indice General") %>%
select(date, ponderador, variacion_mes, f, key) %>%
filter(key == "forecast") %>%
group_by(date, f, key) %>%
summarise(
variacion_mes = sum((ponderador/100) * variacion_mes)
) %>%
ungroup() %>%
mutate(f = paste0(f, " (agregación)"))
) %>%
ungroup() %>%
# Voy a tratar de que todas empiecen justo donde
# termina la serie observada
# (La verdad nunca había hecho esto,
# Es muy probable que haya una forma más
# elegante de hacerlo)
add_row(
date = ymd("2019-08-01"),
key = "forecast",
f = c("auto.arima", "ets", "bats",
"auto.arima (agregación)",
"bats (agregación)", "ets (agregación)"),
variacion_mes = 0.003415041
) %>%
# organizando el orden de las series
# para que las líneas punteadas correspondan
# al forecast
mutate(
key = factor(
key,
levels = c("observado", "forecast"))
) %>%
# restringiendo el espacio temporal
filter(date > "2017-12-01") %>%
ggplot(
aes(x = date, y = variacion_mes, linetype = key, color = f)
) +
geom_line(size = 1) +
# el argumento labels puede recibir una función
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
theme(
legend.position = "bottom"
) +
guides(linetype = "none") +
labs(
x = "",
y = "Variación mensual",
color = "",
title = "Forecast de la inflación mensual",
subtitle = "Resultados según metodología y enfoque de pronóstico"
)## `summarise()` has grouped output by 'date', 'f'. You can override using the
## `.groups` argument.
Comentarios finales
El ejercicio concluye sin dar muchos detalles sobre el resultados de los pronósticos. Advertí desde el principio que la selección y evaluación de los modelos se realizaría en otra publicación, por eso el set de datos utilizado solo incluye datos hasta agosto del 2019, el resto de las observaciones serán utilizadas para probar el desempeño de cada una de las metodologías.
comments powered by Disqus