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. [video] Análisis descriptivo en R
2. Clase 2.2. [video] 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. Macroeconometría aplicada en R: Para consolidar todos los conocimientos en casos prácticos.
1. Clase 7.1. Introducción aplicada (R intermedio)
2. Análisis univariado
3. Modelos lineles
4. Vectores Autorregresivos (VAR)
5. Vector de corrección de errores (VEC)

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

Obtener el ciclo HP para múltiples series de forma simultánea en R

En la siguiente entrada se crea una función sencilla que combinada con la función vectorizada across, permite obtener el componente cíclico de las series numéricas contenidas en una base de datos. Para tales fines se usan datos de la economía de Canada, disponibles en R en la base de datos con el mismo nombre.

library(dplyr)
library(mFilter)
library(ggplot2)
library(mFilter)
 data("Canada")
Canada %>%  head()
               e     prod       rw    U
1980 Q1 929.6105 405.3665 386.1361 7.53
1980 Q2 929.8040 404.6398 388.1358 7.70
1980 Q3 930.3184 403.8149 390.5401 7.47
1980 Q4 931.4277 404.2158 393.9638 7.27
1981 Q1 932.6620 405.0467 396.7647 7.37
1981 Q2 933.5509 404.4167 400.0217 7.13

Partir de una representación grafica de las series siempre es una buena opción:

 Canada %>% autoplot(facet = T)

Ahora creamos la función cicloHP, que recupera el componente cíclico del logaritmo de un vector de datos, usando la función hpfilter de la librería mFilter.

 cicloHP <- function(x) {
  lnx <- log(x)
  clnx <- mFilter::hpfilter(x,freq=50,type="frequency",drift=T)["cycle"]
  unlist(clnx)
}
 
cdata <- Canada %>%
  data.frame() %>%
  mutate(across(where(is.numeric), ~cicloHP(.x)))

Finalmente usamos la referida función:

 ts(cdata) %>%
  autoplot(facet = T)

Representar múltiples gráficos (correlogramas) desde una lista en R

En la siguiente entrada se muestra como guardar múltiples gráficos en una lista, para posteriormente representarlos de manera conjunta. Puntualmente, deseamos graficar el correlograma asociado a cada una de las series disponibles en una base de datos, para posteriormente representarlos en un panel conjunto de gráficos. Para ello se utiliza la base EuStockMarkets, una data disponible en el ambiente base de R, siendo una serie de tiempo sobre cotizaciones de diversos índices financieros.
 
library(purrr)
library(tseries)
library(dplyr)
library(ggplot2)
 
EuStockMarkets %>% head()
 
Time Series:
Start = c(1991, 130)
End = c(1991, 135)
Frequency = 260
             DAX    SMI    CAC   FTSE
1991.496 1628.75 1678.1 1772.8 2443.6
1991.500 1613.63 1688.5 1750.5 2460.2
1991.504 1606.51 1678.6 1718.0 2448.2
1991.508 1621.04 1684.1 1708.1 2470.4
1991.512 1618.16 1686.6 1723.1 2484.7
1991.515 1610.61 1671.6 1714.3 2466.8
 
Puntualmente, queremos guardar un correlograma para cada una de las series. Como se tratan de cotizaciones o precios de diversos índices, obtenemos las tasas de crecimiento de cada columna para generar un correlograma de estas tasas de variaciones, lugar de las cotizaciones (recuerde que será muy seguro que las cotizaciones no se comporten como series estacionarias dado que tienen una clara tendencia. En caso de dudas comparar gráficos de la serie en nivel y en variaciones relativas).
 
my_data <- EuStockMarkets %>%
  as.tibble() %>%
  mutate(
    across(everything(),
           ~((./dplyr::lag(.)-1)*100),
           .names= "tc_{.col}"
    )) %>%
  dplyr::select(starts_with("tc_")) %>%
  na.omit()
 
my_data
# A tibble: 1,859 x 4
   tc_DAX tc_SMI  tc_CAC tc_FTSE
    <dbl>  <dbl>   <dbl>   <dbl>
 1 -0.928  0.620 -1.26    0.679
 2 -0.441 -0.586 -1.86   -0.488
 3  0.904  0.328 -0.576   0.907
 4 -0.178  0.148  0.878   0.579
 5 -0.467 -0.889 -0.511  -0.720
 6  1.25   0.676  1.18    0.855
 7  0.578  1.23   1.32    0.824
 8 -0.287 -0.358 -0.193   0.0837
 9  0.637  1.11   0.0171 -0.522
10  0.118  0.437  0.314   1.41 
# ... with 1,849 more rows
 
Una vez tenemos las tasas de variaciones diarias, ahora convertimos cada fila en una serie temporal (map_df); agregamos una columna para representar los rezagos (bacdf$lag); y reorganizaos la base de datos (dplyr::select(lag, everything())). Note la estructura en que queda representada la data. En caso de tener variables no numérica, como fechas o nombres, o dummy representando eventos o temporadas, estas deben omitirse de esta data que servirá de base para generar los gráficos.
 
bacdf <- map_df(my_data, function(ts) acf(ts, plot = FALSE)$acf)
bacdf$lag <- 0:(nrow(bacdf) - 1)
bacdf <- bacdf %>% dplyr::select(lag, everything())
 
bacdf
# A tibble: 33 x 5
     lag        tc_DAX       tc_SMI       tc_CAC       tc_FTSE
   <int>   <dbl[,1,1]>  <dbl[,1,1]>  <dbl[,1,1]>   <dbl[,1,1]>
 1     0  1        ...  1       ...  1       ...  1        ...
 2     1 -0.000771 ...  0.0475  ...  0.0294  ...  0.0924   ...
 3     2 -0.0265   ... -0.0203  ...  0.00330 ... -0.00811  ...
 4     3 -0.0113   ... -0.0178  ... -0.0461  ...  0.00100  ...
 5     4  0.000469 ...  0.00677 ...  0.00565 ... -0.0240   ...
 6     5 -0.0321   ... -0.0459  ... -0.0313  ... -0.0302   ...
 7     6  0.00196  ... -0.0288  ...  0.00812 ... -0.0524   ...
 8     7 -0.0305   ... -0.00663 ... -0.0486  ... -0.0477   ...
 9     8 -0.00909  ...  0.0260  ... -0.0315  ... -0.000642 ...
10     9  0.0229   ...  0.00554 ...  0.0240  ...  0.0278   ...
# ... with 23 more rows
 
Ahora generamos tres argumentos: el nombre de la variable donde tenemos el rezago (lag); el nombre de las variables a representar en mi base de datos y una lista vacía donde guardaremos correlograma asociado a cada serie.
 
myaxis <- colnames(bacdf[1])
mynames <- colnames(bacdf[-1])
graficos <- list()
 
Posteriormente usamos una secuencia for usando la función seq_along(), note que en cada vuelta del bucle va a llamar el nombre de la variable que estamos usando. Ya dentro del bucle se usa ggplot para generar el correlograma, siendo la única diferencia que la variable y se llama como mynames[[i]], de forma tal que este indexado al valor de i, y por tanto se vaya modificando en cada serie. Estos se van guardando en la posición (graficos[[i]]) del mismo nombre de la lista que hemos creado.
 
for (i in seq_along(mynames)) {
  print(i)
 
  g  <-  ggplot(bacdf, aes_string(x = myaxis, y = mynames[[i]])) +
    geom_segment(mapping = aes(xend = lag, yend = 0)) +
    geom_point() +
    geom_hline(yintercept = c(significance_level, -significance_level), lty = 3, color = "blue") +
    ggtitle(mynames[[i]]) +
    theme_minimal() +
    coord_cartesian(xlim = c(2, 35), ylim=c(-.1,0.1))
 
  graficos[[i]] <- ggplotGrob(g)
 
}
 
Finalmente, generamos el gráfico.
 
grid.arrange(grobs=graficos, ncol=2)

Agradecimientos y créditos: el código base usado dentro del bucle para generar el grafico final se obtuvo de la entrada Plotting multiple ACF with ggplot disponible en stackoverflow.com.

Aplicaciones de la función Rollaplay (de R) en la econometría financiera

 1.       Lógica detrás de la función

La función crea una ventana móvil sobre la cual va calculando algo (mediante una función) y posteriormente podemos obtener información histórica sobre ese cálculo. Imaginemos que todos los días calculáramos la varianza del precio del arroz, si lo hacemos todos los días y vamos guardando este precio, terminaremos teniendo una serie promedio de precios del arroz. Ahora, si tenemos información de 10 años necesitamos repetir esta operación para cada día dentro de la muestra de datos, es aquí donde surge la importancia de esta función.  

Supongamos, tenemos el vector x=(x1,x2,x3,x4,x5), podemos obtener un promedio móvil usando una ventada de tres días. Note que la ventana de datos móvil que va utilizando el programa va dejando una observación detrás y tomando la observación siguiente (en caso de solo ir tomando una observación adicional, partiendo siempre desde la primera observación, se llaman datos recursivos). Los datos vacios pueden omitirse o quedar como NA.

x	Datos móvil
X1
X2
X3	x1,x2,x3
X4	x2,x3,x4
X5	x3,x4,x5
X6	x4,x5,x6

En R, este procedimiento se puede realizar utilizando la función rollapply, la misma segmenta las observaciones de determinados grupos de tamaño (width = 3), y posteriormente realiza el cálculo indicado en la función que se le brinda de argumento a la función.

library(dplyr)
library(lubridate)
library(zoo)
 
x_ej <- zoo(c(3,4,2,8,9,3))
 
rollapply(data = x_ej, width = 3, FUN = mean)
       2        3        4        5
3.000000 4.666667 6.333333 6.666667

También podemos crear nuestras propias funciones para alimentar los procedimientos dentro de la función rollapply:

mean2 = function(x){
  sum(x)/length(x)
}
 
rollapply(data = x_ej, width = 3, FUN = mean2)
   2        3        4        5
3.000000 4.666667 6.333333 6.666667

Alinear los datos arriba o debajo del vector de datos sobre el cual estamos aplicando la función, además de determinar cuál valor queremos sustituir dentro de las observaciones que quedan sin datos:

media<-rollapply(data = x_ej, width = 3, FUN = mean2, align = "right", fill=NA)
media2<-rollapply(data = x_ej, width = 3, FUN = mean2, align = "left", fill=NA)
 
tabla1 <- cbind(x_ej,media,media2)
tabla1
 
  x_ej    media   media2
1    3       NA 3.000000
2    4       NA 4.666667
3    2 3.000000 6.333333
4    8 4.666667 6.666667
5    9 6.333333       NA
6    3 6.666667       NA

También podemos realizar el cálculo a partir de x grupos no superpuesto. En otras palabras, calcular nuestro promedio cada x cantidad de observaciones utilizando el argumento by dentro de la función.

x	Datos móvil
X1	.
X2	.
X3	x1,x2,x3
X4	.
X5	.
X6	x4,x5,x6

media_by <- rollapply(data = x_ej, width = 3, FUN = mean2, align = "right", fill=NA, by=3)
 
tabla1 <- cbind(tabla1,media_by)
tabla1
x_ej    media   media2 media_by
1    3       NA 3.000000       NA
2    4       NA 4.666667       NA
3    2 3.000000 6.333333 3.000000
4    8 4.666667 6.666667       NA
5    9 6.333333       NA       NA
6    3 6.666667       NA 6.666667

Los valores NA también se pueden completar con información parcial, en este caso se usarán las ventadas disponibles para la estimación. Esto solo es recomendable en algunos procedimientos estadísticos, porque otros están influenciados en la cantidad de observaciones usadas. 

x	Datos móvil
X1	x1
X2	x1,x2
X3	x1,x2,x3
X4	x2,x3,x4
X5	x3,x4,x5
X6	x4,x5,x6

media_par <- rollapply(data = x_ej, width = 3, FUN = mean, fill=NA, align = "right", partial = T)
 
tabla1 <- cbind(tabla1,media_par)
tabla1
x_ej    media   media2 media_by media_par
1    3       NA 3.000000       NA  3.000000
2    4       NA 4.666667       NA  3.500000
3    2 3.000000 6.333333 3.000000  3.000000
4    8 4.666667 6.666667       NA  4.666667
5    9 6.333333       NA       NA  6.333333
6    3 6.666667       NA 6.666667  6.666667

También se puede modificar el argumento width para modificar el tamaño de la ventana a voluntad y generar ventanas recursivas, es decir, donde el tamaño de la ventana vaya creciendo con cada estimación.

x	Datos móvil
X1	x1
X2	x1,x2
X3	x1,x2,x3
X4	x1,x2,x3, x4
X5	x1,x2,x3,x4,x5
X6	x1,x2,x3,x4,x5,x6

obs<-length(x_ej)
mean_re <- rollapply(data = x_ej, width = 1:obs, FUN = mean, align = "right")
 
tabla1 <- cbind(tabla1,mean_re)
tabla1
 
x_ej    media   media2 media_by media_par  mean_re
1    3       NA 3.000000       NA  3.000000 3.000000
2    4       NA 4.666667       NA  3.500000 3.500000
3    2 3.000000 6.333333 3.000000  3.000000 3.000000
4    8 4.666667 6.666667       NA  4.666667 4.250000
5    9 6.333333       NA       NA  6.333333 5.200000
6    3 6.666667       NA 6.666667  6.666667 4.833333

2.       Aplicaciones de la función rollapply a la econometría financiera

2.1. Ventana de volatilidad (ventana móvil)

Usamos la data EuStockMarkets para estos ejemplos. El primero calcula las tasas de variación para estas series (recuerde siempre trabajar con series estacionarias) para posteriormente estimar la volatilidad de la misma en una ventana de datos de 100 días. Es un proceso muy utilizado en la gestión del riesgo para analizar la volatilidad condicional.

plot(EuStockMarkets)
 
data <- EuStockMarkets  %>%
  data.frame() %>%
  mutate(across(is.numeric, ~((./dplyr::lag(.)-1)*100))) %>%
  na.omit()
 
r_data0 <- data[,"DAX"]
 
vol_t <- rollapply(r_data0, width=100, FUN=sd)
plot(vol_t, type="l")

2.2. Ventana de volatilidad (ventana recursiva)

 ***

2.3. Calcular el VaR histórico de un activo

Aquí también podemos obtener la evolución del VaR de determinada serie financiera. Se muestra además la flexibilidad de la función al permitirnos colocar dos argumentos a la vez.

rollapply(data, width, FUN, …, )

var_dax <- rollapply(r_data0, 100, quantile, probs = 0.01, align = "right")
 
plot(r_data0, type="l")
lines(var_dax)

#backtesting
mean(r_data0<var_dax)
[1] 0.01667563

2.4. Correlación histórica entre activos

La función no solo permite incluir más de un argumento de una función, sino que facilita el trabajo con más de una variable. Aquí se muestra como trabajar con varias columnas. Este es necesario cuando quisiéramos ver cuestiones como la evolución histórica de la correlación entre dos series y como vimos anteriormente, cuando necesitamos incluir funciones propias para reservar los datos que vamos obteniendo. Esta función la podemos crear fuera o colocar dentro de la misma función.

r_data <- data[,c(1,2)]
 
corr_t <- rollapply(r_data, width=100, function(x) cor(x[,1],x[,2]), by.column=FALSE)
plot(corr_t, type="l")

 

2.5. Obtener el beta histórico de una cartera

 Sumamente útil para descomponer el riesgo de un activo e identificar el componente sistémico de su volatilidad, además de identificar fuentes de riesgo y la sensibilidad de la serie a factores puntuales de riesgo. El procedimiento consiste en guardar el coeficiente de regresión de un modelo factorial sencillo, pero note que modificando la función podemos obtener modelos factoriales con mayor cantidad de factores de riesgo.  

beta_sis = function(x){
  lm(x[,2] ~ x[,1])[[1]][[2]]
}
 
Beta <- rollapply(r_data, width=100, by.column = F,  beta_sis)
plot(Beta, type="l")

Verifique que se usa el argumento by.column, cuyo valor por default es TRUE, esto para evitar que la función se aplique a cada una de las columnas con las que estamos trabajando. En los casos anteriores, como solo teníamos un vector como argumento de datos, no era necesario preocuparnos por este argumento de la función.

Note que esta función la podemos modificar para obtener más de un resultado de la función con la que estamos interesados. Por ejemplo, podemos estar enterados en recuperar el intercepto y la pendiente del modelo de regresión, o adicionalmente recuperar el R2.

beta_sis2 = function(x){
  lmr <-lm(x[,2] ~ x[,1])
  c(beta=coef(lmr)[[2]], r2=summary(lmr)$r.squared)
}
 
Beta_r2 <- rollapply(r_data, width=100, by.column = F,  beta_sis2)
 
par(mfrow=c(2,1))
plot(Beta_r2[,1], type="l")
plot(Beta_r2[,2], type="l")

Una alternativa mostrada en data.camp es:

FUN = function(z) coef(lm(y ~ y1 + y12, data = as.data.frame(z))

3.       Condiciones de la estimación

3.1. Tamaño de la ventana

 Note que no es trivial elegir el tamaño de la ventana. Pues la sensibilidad al contexto, o inclusive la dinámica de la serie puede cambiar de forma significativa en función del tamaño de la ventana seleccionada.

 #diferentes ventanas de estimation
vol_t <- rollapply(r_data0, width=100, FUN=sd,align = "right", fill=NA)
vol_t500 <- rollapply(r_data0, width=500, FUN=sd,align = "right", fill=NA)
vol_t1000 <- rollapply(r_data0, width=1000, FUN=sd, align = "right", fill=NA)
 
plot(vol_t, type="l")
lines(vol_t500)
lines(vol_t1000)

3.2. Ventana recursiva vs. ventana móvil

 ***

4.       Referencias

DataCamp. 2022. rollapply: Apply Rolling Functions.
 
Dumitrescu, E.; Hurlin, C.; Pham., V. Backtesting Value-at-Risk: From Dynamic Quantile to Dynamic Binary Tests.
 
Grothendieck, K. (2020). Rollapply for backtesting the value at risk. stackoverflow.
 
Grothendieck, K. (2021). Using Rollapply to return both the Coefficient and RSquare. stackoverflow.
 
Kleiber, C. (2017). Applied Econometrics whith R.
 
Simaan, M. (2018). Tip of the Month: rollapply.
 
Wickham., H. (2021). Advanced R. Functionals.