Análisis del Impacto de un Choque Monetario sobre el Desempleo usando VAR en R

 En este post, exploraremos cómo estimar un modelo VAR en R utilizando la base de datos Canada del paquete vars. Se intenta mostrar cómo usar impulso respuesta de un VAR y como generar escenarios a partir de la introducción de choques. Nos enfocaremos en analizar la respuesta del desempleo ante un choque monetario, graficar la función de impulso-respuesta (IRF) y proyectar la tasa de desempleo con y sin la presencia del choque.

Antes de comenzar, asegurémonos de tener instalados y cargados los paquetes necesarios:
 
library(vars)
library(ggplot2)
library(forecast)

Utilizaremos la base Canada, que contiene datos macroeconómicos de Canadá. Seleccionaremos la tasa de desempleo (U) y la tasa de interés (R) para modelar su relación mediante un VAR con dos rezagos.
 
# Cargar la base de datos
data(Canada)
Canada_ts <- as.data.frame(Canada)
 
# Estimar un modelo VAR con 2 rezagos
desempleo_VAR <- VAR(Canada_ts, p = 2, type = "const")
 
Funciones de Impulso-Respuesta (IRF)

Ahora obtenemos la función de impulso-respuesta para analizar el efecto de un choque monetario (en la tasa de interés R) sobre el desempleo (U).
 
irf_result <- irf(desempleo_VAR, impulse = "rw", response = "U", n.ahead = 10, boot = TRUE)
 
irf_df <- data.frame(
  Horizon = 1:11,
  Response = irf_result$irf$rw |> as.vector(),
  Lower = irf_result$Lower$rw |> as.vector(),
  Upper = irf_result$Upper$rw |> as.vector()
)
 
ggplot(irf_df, aes(x = Horizon, y = Response)) +
  geom_line(color = "blue", size = 1) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "blue", alpha = 0.2) +
  labs(title = "Respuesta del Desempleo a un Choque Monetario",
       x = "Horizonte",
       y = "Impacto en el Desempleo") +
  theme_minimal()
 

Proyección del Desempleo sin Choque
 
Ahora realizamos una proyección de la tasa de desempleo para los próximos 10 períodos sin ningún choque:
 
forecasts <- predict(desempleo_VAR, n.ahead = 10)
U_forecast <- data.frame(
  Horizon = 1:10,
  Forecast = forecasts$fcst$U[, 1]
)

Proyección del Desempleo con un Choque Monetario

En un modelo VAR (Vector Autoregressive), un choque se refiere a una perturbación o innovación inesperada que afecta a una o más variables del sistema. Estos choques representan eventos imprevistos que alteran la trayectoria normal de las variables económicas o financieras modeladas.

Un choque en un modelo VAR tiene las siguientes características:
·        Es una perturbación exógena que entra al sistema a través de los términos de error en las ecuaciones del modelo
·        Puede afectar simultáneamente a múltiples variables debido a la estructura multivariante del VAR
·        Permite analizar la propagación de efectos a lo largo del tiempo (mediante funciones de impulso-respuesta)
·        Ayuda a entender cómo un cambio inesperado en una variable afecta a las demás variables del sistema
 
Introducimos un choque en la tasa de interés (+0.5 puntos porcentuales en el primer período) y ajustamos la predicción:
 
forecasts <- predict(desempleo_VAR, n.ahead = 10)
 
U_forecast <- data.frame(
  Horizon = 1:10,
  Forecast = forecasts$fcst$U[, 1]
)
 
  shock <- rep(0, 10)
shock[1] <- 0.5  # Simulamos un aumento de 0.5 pp en la tasa de interés
response_with_shock <- irf_result$irf$rw[-11] + shock  # Aplicamos el choque a la IRF
 
U_forecast_shock <- data.frame(
  Horizon = 1:10,
  Forecast = U_forecast$Forecast + response_with_shock
)
 
ggplot() +
  geom_line(data = U_forecast, aes(x = Horizon, y = Forecast, color = "Sin Choque"), size = 1) +
  geom_line(data = U_forecast_shock, aes(x = Horizon, y = Forecast, color = "Con Choque"), size = 1) +
  labs(title = "Proyección del Desempleo con y sin Choque Monetario",
       x = "Horizonte",
       y = "Tasa de Desempleo") +
  scale_color_manual(name = "Escenario", values = c("Sin Choque" = "black", "Con Choque" = "red")) +
  theme_minimal()

Curso completo de R: de 0 a avanzando

 El curso que he diseñado (claude.ia) para ti sigue una estructura progresiva que te permitirá dominar R desde los conceptos básicos hasta técnicas avanzadas. Lo he organizado en 7 módulos principales:

1. Introducción a R: Cubre instalación, tipos de datos, estructuras fundamentales e importación/exportación.
1. Clase 1.1. [video] Introducción a R. Tipo de datos, operadores, funciones, consola, editor, entorno, visualizador, entorno de trabajo, estructura de datos:  creación, indexación y operaciones. Importación de datos y funciones integradas.
2. Análisis Estadístico: Incluye estadística descriptiva, manejo de datos atípicos, imputación de valores perdidos y pruebas de independencia, tal como solicitaste.
1. Clase 2.1. Análisis descriptivo en R
2. Clase 2.2. Análisis de independencia
3. Clase 2.3. Valores atípicos
4. Clase 2.4. Impuetación de datos en R
3. Tidyverse: Abarca todo el ecosistema con especial énfasis en dplyr, ggplot2 y manipulación de datos relacionales (long, join, wide).
1. Clase 3.1. Dplyr. Análisis de base de hogares
2. Clase 3.2. Ggplot
3. Clase 3.3. Dplur II. Bases relacionales (join, long data, nest...)
4. Programación en R: Desarrolla bucles, funciones, estructuras y vectorización, con una clase completa dedicada al paquete purrr.
1. Clase 4.1. Introducción a la programación
2. Clase 4.2. Funciones propias
3. Clase 4.3. Operaciones vectorizadas + purr
5. Modelado Estadístico y Machine Learning: Módulo adicional que te permitirá aplicar R en contextos analíticos avanzados.
6. Desarrollo Avanzado en R: Para convertirte en un profesional completo, incluyendo optimización, desarrollo de paquetes y aplicaciones.
7. Proyectos Integradores: Para consolidar todos los conocimientos en casos prácticos.

Recodificación de variables usando dplyr en R

Una base de datos suele tener diversos tipos de variables del tipo cualitativo y cuantitativo. En función del tipo de variables aplicamos diversos tipos de herramientas estadísticas para el análisis de datos. Por ejemplo, el estudio de la independencia entre dos variables depende de la forma de la misma:

* Dos variables cuantitativas: análisis de regresión o correlación

* Una variable cuantitativa y una cualitativa: test de media o comparación de media.

-          * Dos variables cualitativas: test chi-cuadrado.

Por tanto, pueden existir casos donde estemos interesados en recodificar una de estas variables para aplicar algún tipo de análisis, como convertir una variable cuantitativa en cualitativa o transformar una variable cualitativa para recalificar grupos. Por ejemplo agregar edades por rango: 0-10, 11-20, 21-30, …; o reagrupar grupos dentro de una variable cualitativa, ejemplo, agrupar los estudiantes de economía y contabilidad como miembros de la facultad de económicas.

Para ilustrar estos ejemplos, simulamos la siguiente base de datos:

set.seed(1)  # Para reproducibilidad
library(tidyverse)
 
# Simulación de datos
data_simulada <- tibble(
  edad = sample(0:90, 100, replace = TRUE),
  sexo = sample(c("Masculino", "Femenino"), 100, replace = TRUE),
  ingresos = round(rnorm(100, mean = 30000, sd = 10000), 2),
  zona_residencia = sample(c("a", "b", "c"), 100, replace = TRUE))
 
# A tibble: 100 × 4
    edad sexo      ingresos zona_residencia
   <int> <chr>        <dbl> <chr>         
 1     7 Femenino    13113. a             
 2    24 Masculino   14276. a              
 3    19 Femenino    25950. a             
 4    83 Masculino   33193. a             
 5     4 Masculino   30404. c             
 6    73 Femenino    26100. a             
 7    19 Masculino   11808. c             
 8    54 Masculino   36592. b             
 9    55 Femenino    34596. c             
10    78 Femenino    46166. C

 1.       Análisis de correlación para estudiar dependencia entre dos variables cuantitativas:

 ggplot(data_simulada, aes(x = edad, y = ingresos)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "blue", se = FALSE) +
  labs(x = "Edad", y = "Ingresos", title = "Relación entre Edad e Ingresos") +
  theme_minimal()

 2.       Tabla de comparación de promedios, para estudiar independencia entre una variable cualitativa y otra cuantitativa:

 data_simulada |>
  group_by(sexo) |>
  summarize(ingreso_promedio = mean(ingresos, na.rm = TRUE))
 
# A tibble: 2 × 2
  sexo      ingreso_promedio
  <chr>                <dbl>
1 Femenino            29794.
2 Masculino           29546.

 3.       Recodification de variables.

 Ahora, queremos estudiar la relación entre edad e ingreso vista en el apartado 1, pero con el fin estudiar la independencia a partir del ejemplo visto en el apartado 2. Es decir, necesitamos convertir una variable cuantitativa en cualitativa para realizar algún tipo de análisis. Un primer ejemplo de este tipo de análisis es transformar una variable cuantitativa en grupos específicos, con el objetivo de resaltar clúster de observaciones con algún objetivo.

 3.1. Recodificar una variable cuantitativa, creando grupos arbitrarios.

 Según edad, por ejemplo, podemos crear grupos segmentarlo según edad del trabajador.

 data_simulada |>
  mutate(grupo_edad = case_when(
    edad >= 0 & edad <= 14 ~ "0-14",
    edad >= 15 & edad <= 64 ~ "15-64",
    edad > 65 ~ ">65",
    TRUE ~ NA_character_  # Manejo de posibles NA
  )) |>
  group_by(grupo_edad) |>
  summarise(media = mean(ingresos))
 
# A tibble: 3 × 2
  grupo_edad  media
  <chr>       <dbl>
1 0-14       29160.
2 15-64      28094.
3 >65        32791

Otro ejemplo de este tipo de recodificación es crear una variable con distintos niveles de ingresos: bajo, medio, alto. Podemos ahora obtener edad promedio según nivel de ingresos.

 data_simulada %>%
  mutate(
    nivel_ingreso = case_when(
      ingresos < 20000 ~ "Bajo",
      ingresos >= 20000 & ingresos < 40000 ~ "Medio",
      ingresos >= 40000 ~ "Alto"
    )
  ) |> 
  group_by(nivel_ingreso) |>
  summarize(edad_m_ing = mean(edad, na.rm = TRUE))
 
# A tibble: 3 × 2
  nivel_ingreso edad_m_ing
  <chr>              <dbl>
1 Alto                51.6
2 Bajo                44.7
3 Medio               42.1

 3.2. Recodificar una variable cuantitativa, creando grupos estadísticos: quintiles.

 Ahora, si quisiéramos agrupar los ingresos, no tenemos de antemano una segmentación por grupos que sea representativa como la edad[1]. En estos casos, donde no tenemos una segmentación orientada a caracterizar una población, podemos optar por grupos estadísticos, como los percentiles.

 Por ejemplo: podemos ver la edad promedio pro percentil de ingresos.

 data_simulada |>
  mutate(quintil_ingresos = ntile(ingresos, 5)) |>  # Dividir ingresos en 5 quintiles
  group_by(quintil_ingresos) |>
  summarize(edad_promedio = mean(edad, na.rm = TRUE))
 
# A tibble: 5 × 2
  quintil_ingresos edad_promedio
             <int>         <dbl>
1                1          44.1
2                2          43.8
3                3          41.8
4                4          44.4
5                5          59.2

 ·         mutate(quintil_ingresos = ntile(ingresos, 5)): Crea una nueva columna quintil_ingresos que divide ingresos en 5 grupos de igual tamaño.

·         group_by(quintil_ingresos): Agrupa los datos por cada quintil de ingresos.

·         summarize(edad_promedio = mean(edad, na.rm = TRUE)): Calcula la edad promedio dentro de cada quintil.

 3.3. Recodificar una variable cualitativa para crear una variable binaria.

Este ejemplo aplica para variables cuantitativas, al recodificar una variable para crear un grupo que sea igual a 1. Por ejemplo, podemos querer que menores de edad sean 1, o que mujeres sean menores a 1.

Recodificar Sexo a Valores Binarios. Convertimos sexo en variables binarias, codificando "Masculino" como 1 y "Femenino" como 0. Al grupo base se coloca el valor de 1.

 data_simulada %>%
  mutate(hombre_binario = if_else(sexo == "Masculino", 1, 0)) |> 
  group_by(hombre_binario) |>
  summarize(edad_sexo = mean(ingresos, na.rm = TRUE))
 
# A tibble: 2 × 2
  hombre_binario edad_sexo
           <dbl>     <dbl>
1              0    28807.
2              1    29658.

 3.4. Recodificar una variable cualitativa para recalificar grupos.

 Finalmente, recodificamos una variable cualitativa como zona, para verificar que la zona a y b sean consideradas como zona rural y c como urbana. Esto se hace para obtener un grupo cuya calificación sean mas familiar para el lector, o este mas acorde como los objetivos de la investigación. Primero creamos la variable en la base de datos, mediante una recodificación de la variable zona, usando case_when para crear la nueva columna zona en la que las categorías "a" y "b" se recodifican como "Rural" y "c" como "Urbana":

 data_simulada |>
  mutate(zona = case_when(
    zona_residencia %in% c("a", "b") ~ "Rural",
    zona_residencia == "c" ~ "Urbana"
  ))
 
# A tibble: 100 × 5
    edad sexo      ingresos zona_residencia zona 
   <int> <chr>        <dbl> <chr>           <chr>
 1     7 Femenino    13113. a               Rural
 2    24 Masculino   14276. a               Rural
 3    19 Femenino    25950. a               Rural
 4    83 Masculino   33193. a               Rural
 5     4 Masculino   30404. c               Urbana
 6    73 Femenino    26100. a               Rural
 7    19 Masculino   11808. c               Urbana
 8    54 Masculino   36592. b               Rural
 9    55 Femenino    34596. c               Urbana
10    78 Femenino    46166. c               Urbana
# ℹ 90 more rows
# ℹ Use `print(n = ...)` to see more rows

Ahora calculamos el ingreso promedio por la zona de residencia.

 data_simulada |>
  mutate(zona = case_when(
    zona_residencia %in% c("a", "b") ~ "Rural",
    zona_residencia == "c" ~ "Urbana"
  ))   |>
  group_by(zona) |>
  summarize(ingreso_promedio = mean(ingresos, na.rm = TRUE))
 
# A tibble: 2 × 2
  zona   ingreso_promedio
  <chr>             <dbl>
1 Rural            29684.
2 Urbana           28133.

---

[1] Salvo en el caso de la pobreza, que podemos agrupar según la línea de pobreza.

Evaluación de pronóstico en R: validación cruzada de modelos de series temporales

En el siguiente ejemplo se muestra una validación cruzada de series temporales en R. Necesita evaluar la capacidad de pronóstico de un modelo para elegir alguna alternativa en función de su rendimiento. Es decir, necesito verificar cual hubiese sido el pronóstico histórico de una serie temporal, por ejemplo con un modelo AR(1), para comparar con la precisión de otros modelos y decidir cual estrategia seguir en lo adelante. En este sentido, necesitaríamos realizar el ejercicio de proyección en distintos puntos temporales, para posteriormente recuperar los errores y comparar estas medidas de error entre distintos modelos.

 1.       Simulamos la serie temporal (en el caso de una base de datos, deberías importarlas en R).

 library(writexl)   # Exportar a Excel
library(readxl)
library(lubridate) # dates
library(tidyverse)
library(forecast)
library(timeSeries)
library(zoo)
 
# Importando datos ---------------
# ********************************************************************
# ********************************************************************
 
# Configuración inicial
set.seed(4)    # Asegura reproducibilidad
n <- 250       # Número de días (simulación de aproximadamente un año laboral)
mu <- 0        # Promedio del cambio diario (sin tendencia)
sigma <- 0.01  # Volatilidad diaria
 
# Generación de la caminata aleatoria
tipo_cambio <- cumsum(c(1, rnorm(n, mean = mu, sd = sigma)))  # Inicia en 100 y suma los cambios diarios
 
# Creación de un data.frame con fechas
fechas <- seq.Date(from = as.Date("2023-01-01"), by = "day", length.out = n + 1)
df <- data.frame(fecha = fechas, tipo_cambio = tipo_cambio)
 
# Gráfico de series temporales
library(ggplot2)
ggplot(df, aes(x = fecha, y = tipo_cambio)) +
  geom_line(color = "blue") +
  labs(title = "Simulación de Caminata Aleatoria del Tipo de Cambio",
       x = "Fecha",
       y = "Tipo de Cambio") +
  theme_minimal()
 
library(xts)
tipo_cambio_xts <- xts(tipo_cambio, order.by = fechas)

 
2.       Validación cruzada

En resumen, lo que se hace es definir una ventana de estimación y otra de validación (pronóstico). En este ejemplo son de 60 y 5 observaciones respectivamente. 60 observaciones serán usadas para estimar el modelo, mientras que 5 será la longitud de la ventana del pronóstico. De esta forma, en la primera iteración se va estimar los modelos con observaciones del 1:60 y otra ventana de pronóstico de la 61-65, luego 2:61 y 62:66, … y así respectivamente hasta recorrer toda la muestra. El documento de trabajo “Pronóstico del consumo privado en la República Dominicana” de Alejandro J. Balcácer y Nerys F. Ramírez, muestra en detalle esta estrategia.

 # Análisis de precisión. Error histórico ---------------

h.valid <- 5
horizonte <- paste0("h", 1:h.valid)
vent <- 60

Posteriormente, se genera un for para realizar este recorrido usando una serie. Primero se define el valor del ultimo dato observado en cada iteración del bucle (seq_for), este dato sirve de referencia para el resto de las ventanas. Además, se define un contador para el punto de inicio de la ventana de estimación (cont). Los puntos clave de esta secuencia son:

-          Indicar la secuencia con la posición de las observaciones que vamos a usar en cada punto de la serie temporal [obser in seq_for].

-          La serie temporal se segmenta en dos partes: la estimación (tipo_cambio_xts[cont:obser]) y la otra parte de pronósticos (tipo_cambio_xts[(obser+1):(obser+h.valid)]). Fíjese que la primera va desde el contador hasta obser, recuerde que obser es el valor final hasta donde irá la ventana de estimación y el valor que en el presente ejemplo sirve de referencia. Luego, la parte para evaluar las proyecciones se definen a partir de la última observación de la parte de estimación (obser+1):(obser+h.valid).

-          La tercera parte importante del bucle es incluir el modelo de pronóstico. Aquí debemos usar alguna estrategia para proyectar los periodos que nos resulten de interés. En este ejemplo son 5 observaciones, y se usa un modelo ar(1).

-          Finalmente, se crea un data.frame para guardar el pronóstico, el dato observado, los horizontes y el nombre del model. Además, se guarda cada data.frame en una lista creada anteriormente (data_Rbind). Primero se llaman todos los objetos que inician con “df_” (tener pendiente guardar cada nuevo modelo con ese nombre).  

 # Esto lo que hace, es que permite recorrer las ventanas
seq_for <- vent:(nrow(tipo_cambio_xts)-h.valid)
 
# OJO CORRER EL BULCE DESDE AQUI
cont <- 1
data_Rbind <- list() # Lista que guarda resultados
 
for (obser in seq_for){
  
  tc_train <- tipo_cambio_xts[(obser-vent+1):obser]
  tc_test <- tipo_cambio_xts[(obser+1):(obser+h.valid)]
 
  fecha_ciclo <- fechas[obser]
  fecha_hat <-  fechas[(obser+1):(obser+h.valid)]
 
ordena_data <- function(objetos,...){
  # objetos = lista de nombre de vectore
 
  # Iterar sobre los nombres de los objetos y cambiar los nombres de las columnas
  for (i in seq_along(objetos_df)) {
    xdf_name <- get(objetos_df[i])
    new_column_names <- c("fechaCiclo","fechaHat","horizonte","obs","pro","model")
    names(xdf_name) <- new_column_names
    assign(objetos_df[i], xdf_name)
  }
 
  lista_objetos <- lapply(objetos_df, get)
  df1_final <- dplyr::bind_rows(lista_objetos)
  return(df1_final)
}
    ### 1. ESTRATEGIA AR(1) . Ventana móvil ************************
 
    model_ar1 <- Arima(tipo_cambio_xts, order=c(1,0,0))
    ht_f_ar1 <- forecast::forecast(model_ar1, h=h.valid)$mean                 
   
    df_ar1 <- data.frame(fecha_ciclo, fecha_hat, horizonte,  #fechas
                         obsT=tc_test, # dato observado
                         proT=ht_f_ar1,
                         model="1.Ar(1)_móvil")
 
    ### 1. ESTRATEGIA AR(2,1) . Ventana móvil ************************
   
    model_ar21 <- Arima(tipo_cambio_xts, order=c(2,0,1))
    ht_f_ar21 <- forecast::forecast(model_ar21, h=h.valid)$mean                 
   
    df_ar21 <- data.frame(fecha_ciclo, fecha_hat, horizonte,  #fechas
                         obsT=tc_test, # dato observado
                         proT=ht_f_ar21,
                         model="2.Ar(21)_móvil")
   
### ******************************************************************
### ******************************************************************
objetos <- ls()
objetos_df <- objetos[grep("^df_", objetos)]
 
data_Rbind[[(obser-vent+1)]] <- ordena_data(objetos_df)
   
# Esto al final siempre       
cont <- cont+1  # cont y esto es equivalente (obser-vent+1):
}

Luego de culminar, se utiliza una función propia para ordenar la data (ordena_data) y convertirla en un único data.frame con toda la data:

data_forecast_final
# A tibble: 1,870 × 6
   fechaCiclo fechaHat   horizonte   obs  proy model        
   <date>     <date>     <chr>     <dbl> <dbl> <chr>        
 1 2023-03-01 2023-03-02 h1        0.969 1.01  1.Ar(1)_móvil
 2 2023-03-01 2023-03-03 h2        0.981 1.01  1.Ar(1)_móvil
 3 2023-03-01 2023-03-04 h3        0.972 1.01  1.Ar(1)_móvil
 4 2023-03-01 2023-03-05 h4        0.978 1.00  1.Ar(1)_móvil
 5 2023-03-01 2023-03-06 h5        0.987 1.00  1.Ar(1)_móvil
 6 2023-03-01 2023-03-02 h1        0.969 1.01  1.Ar(21)_móvil
 7 2023-03-01 2023-03-03 h2        0.981 1.01  1.Ar(21)_móvil
 8 2023-03-01 2023-03-04 h3        0.972 1.01  1.Ar(21)_móvil
 9 2023-03-01 2023-03-05 h4        0.978 1.00  1.Ar(21)_móvil
10 2023-03-01 2023-03-06 h5        0.987 0.999 1.Ar(21)_móvil

 Con esta data, ahora se puede realizar cualquier análisis sobre los errores promedio, su evolución histórica o el comportamiento de estos errores condicionados a eventos económicos, por ejemplo, los errores durante el COVID o periodos de alta incertidumbre.

 # Cuadrando la base de los errores
# Crea la data final con todos los los dtos necesarios
data_forecast_final <- dplyr::bind_rows(data_Rbind) |> as_tibble()
names(data_forecast_final) <- c("fechaCiclo","fechaHat","horizonte","obs","proy","model")
 
# Errores total
# Total
data_forecast_final |>
  #dplyr:: filter(fechaCiclo > as.Date("2018-01-01"))  |>
  mutate(error=obs-proy) |>
  group_by(horizonte, model) |> # (horizonte, model, COVID)
  summarize(
    mae_e_abs = mean(abs(error), na.rm = TRUE),
  ) |>
  pivot_wider(names_from = horizonte, values_from = mae_e_abs) |>
  arrange(h1) |>
  dplyr::select(model,horizonte, everything())
 
`summarise()` has grouped output by 'horizonte'. You can override using the `.groups` argument.
# A tibble: 2 × 6
  model              h1     h2     h3     h4     h5
  <chr>           <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 1.Ar(21)_móvil 0.0479 0.0442 0.0405 0.0377 0.0350
2 1.Ar(1)_móvil  0.0483 0.0443 0.0409 0.0378 0.0353