Funciones vectorizadas en R: tapply

Las funciones vectorizadas en R permiten resumir un conjunto de tareas repetitivas, con algún patrón o secuencia común, o aplicar una operación a un conjunto de elementos. En la presente entrada mostramos algunos ejemplos de cómo usar la función tapply en R, esta permite crear variables por alguna condición agrupada, o collapse de base de datos.

En primer lugar, generamos una base aleatoria:

 set.seed(1)
 
my_datos <- data.frame(
              provincia = sample(c("Ocoa","Azua"), size = 10, replace = TRUE),
              ingreso = round(rnorm(10, sd = 30, mean = 100)),
              edad = round(rnorm(10, sd = 6, mean = 20)))
my_datos
 
  provincia ingreso edad
1       Ocoa      75   20
2       Azua     115   20
3       Ocoa     122   26
4       Ocoa     117   25
5       Azua      91   24
6       Ocoa     145   26
7       Ocoa     112   25
8       Ocoa      81   20
9       Azua      34    8
10      Azua     134   24

 1.1.   Collapse de base de datos

 La función tapply permite realizar operaciones agregadas a partir de alguna clase o característica indicada por alguna variable. Un ejemplo de esta son los promedios condicionados o promedios por grupo.

 tapply(my_datos$ingreso, my_datos$provincia, mean)
Azua     Ocoa
 93.5000 108.6667

También, podemos utilizar funciones propias. Esto sirve para personalizar las operaciones que queremos realizar, sin que necesariamente dependamos de las funcione pre establecidas en R. Un ejemplo de esto es calcular el coeficiente de variación por provincia:

tapply(my_datos$ingreso, my_datos$provincia, function(x) 100*sd(x)/mean(x))
    Azua     Ocoa
46.41021 24.26844

Otro ejemplo es calcular el porcentaje de personas con ingresos por debajo del promedio en cada provincia, que en nuestro ejemplo es de 102.6. Siempre que se esté probando alguna función nueva, es recomendado validar los resultados obtenidos con algunos cálculos manuales, para evitar estar cometiendo algún error de semántica:

tapply(my_datos$ingreso, my_datos$provincia, function(x) mean(x<mean(x))*100)
    Azua     Ocoa
50.00000 33.33333

 1.2.   Variables agregadas

Note usted que la secuencia anterior se ha realizado un collapse de los datos, por lo que, tenemos tantas observaciones como categorías de variable cualitativas contemos (provincias), sin embargo, en muchas ocasiones ameritamos contar con un estadístico para cada una de las variables, es decir, crear una nueva variable con una observación para individuo. Esta observación corresponde al dato agregado. Para esto utilizamos la indexación, mediante el argumento as.character, en la función tapply:

tapply(my_datos$ingreso, my_datos$provincia, mean)[as.character(my_datos$provincia)]
 
   Ocoa     Azua     Ocoa     Ocoa     Azua     Ocoa     Ocoa     Ocoa     Azua     Azua
108.6667  93.5000 108.6667 108.6667  93.5000 108.6667 108.6667 108.6667  93.5000  93.5000

 Otro ejemplo es que podemos colocar el ingreso relativo de cada individuo respecto a su familia:

 Library(dplyr)
 
my_datos <- my_datos %>%
  mutate(ing_rel = ingreso/tapply(ingreso, provincia, mean)[as.character(provincia)])
 
my_datos

1.3.   Tablas de contingencias

Usar una lista para agregar múltiples factores. Para ilustrarlo agregamos la variable edad. y agregamos una tabla especial donde se muestre la edad promedio por rango de edad, según provincia. Primero creamos una variable dummy cuyo 1 corresponde aquellos individuos menores de edad, usando la función mutate del paquete dplyr, posteriormente ilustramos varios ejemplos de tablas de contingencia.

my_datos <- my_datos %>%
            mutate(menor = factor(edad<=17, label = c("PET","Menor")))
 
  provincia ingreso edad   ing_rel menor
1       Ocoa      75   20 0.6901840   PET
2       Azua     115   20 1.2299465   PET
3       Ocoa     122   26 1.1226994   PET
4       Ocoa     117   25 1.0766871   PET
5       Azua      91   24 0.9732620   PET
6       Ocoa     145   26 1.3343558   PET
7       Ocoa     112   25 1.0306748   PET
8       Ocoa      81   20 0.7453988   PET
9       Azua      34    8 0.3636364 Menor
10      Azua     134   24 1.4331551   PET
 
tapply(my_datos$ingreso, my_datos$provincia, mean)
Azua     Ocoa
 93.5000 108.6667

O una lista de variables:

tapply(my_datos$ingreso, list(my_datos$provincia, my_datos$menor), mean)
         PET Menor
Azua 113.3333    34
Ocoa 108.6667    NA

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.

Ejemplos de funciones de impulso respuesta (ifr) usando VARS y ggplot2 en R

En esta primera parte usamos la función across para calcular las tasas de crecimiento de todas las variables en la base de datos. Usamos los datos de EuStockMarkets disponible en R (data()). Esto se hace porque recuerden que el VAR solo se estima con variables estacionarias, por lo que, un análisis formar requeriría hacer test sobre las variables que usaríamos.

Posteriormente, una vez seleccionada la serie que vamos a utilizar (var_data), usamos el paquete VARS para tres cuestiones puntuales:

 1.  VARselect identifica el orden del VAR, es decir, la cantidad a retardos utilizar en el VAR.

2.  VAR es la función que nos permite estimar el VAR.

3.  Irf nos permite obtener los datos asociados a la función impulso respuesta.

#ifr in R
data <- EuStockMarkets  %>%
  data.frame() %>%
  mutate(across(is.numeric, ~((./dplyr::lag(.)-1)*100))) %>%
  na.omit()
 
# var
var_data <- data.frame(x=EuStockMarkets[,"DAX"],y=EuStockMarkets[,"FTSE"])
 
plag <- VARselect(var_data, lag.max = 8, type = "both")$selection[1]
p1ct <- VAR(var_data, p = plag)
irf1 <- irf(p1ct, impulse = "x", response = "y",
            n.ahead = ifr_lag, boot = TRUE,
            cumulative = F, ortho = TRUE, ci = 0.95)

 Finalmente, solo es necesario utilizar la función plot para graficar la función impulso respuesta asociada.

#Opcion 1
plot(irf1)

Adicionalmente, podemos extraer los coeficientes y los límites de confianzas calculados para obtener las IFR estimadas en el paso anterior. Note que este se guarda en un marco de datos (data frame) que permitiría en lo adelante graficar nuestra función atendiendo los criterios de estéticas que prefiramos y usando ggplot2 del tidyverse.

ifr_lag <- 24
data_ifr <- data.frame(lagifr=1:(ifr_lag+1), irf1$irf, irf1$Lower, irf1$Upper)
names(data_ifr) <- c("retardos","ifr","li","ls")
head(data_ifr)
 
retardos      ifr       li       ls
1        1 20.63704 19.08058 22.93661
2        2 22.23282 19.40226 24.77678
3        3 22.47799 19.47766 25.11293
4        4 22.49421 19.39938 25.11469
5        5 22.47920 19.23570 25.06734
6        6 22.45973 19.07158 25.02734

Debajo se colocan tres ejemplos de ifr usando ggplot2. Note que en estos ejemplos se encuentra comentada la línea de intercepto en el 0. Esto para recortar el eje de las x que estamos usando, en sus ejemplos particulares pueden des comentar esta línea del comando borrando el símbolo de #.

#Opcion 2

data_ifr %>%
  ggplot(aes(x = retardos, y = ifr)) +
  geom_line(size=1)+
  geom_line(aes(x = retardos, y = li),linetype=4) +
  geom_line(aes(x = retardos, y = ls),linetype=4)+
  theme_minimal() +
  #geom_hline(yintercept=0) +
  theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank())

 

#Opcion 3
  ggplot(aes(x = retardos, y = ifr)) +
  geom_line(size=1, color="#21618C")+
  geom_ribbon(aes(x = retardos, ymin = li, ymax = ls), col = 'grey',
              fill = 'grey', alpha = 0.3, linetype=4) +
  theme_minimal() +
  #geom_hline(yintercept=0) +
  theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank()) +
  #scale_y_continuous(expand = c(0, 0), breaks = seq(min(data_ifr$li), max(data_ifr$ls, 2)))+
  scale_x_continuous(expand = c(0, 0), breaks = seq(0, ifr_lag, 2))

#Opcion 4
data_ifr %>%
ggplot(aes(x = retardos, y = ifr, ymin = li, ymax = ls)) +
  geom_line(size=1, color="#21618C")+
  geom_ribbon(aes(x = retardos), col = 'grey',
              fill = 'grey', alpha = 0.3, linetype=4) +
  # geom_hline(yintercept = 0, color="red") +
  geom_ribbon(fill="grey", alpha=0.2) +
  theme_light() +
  xlab("") +
  theme(plot.title = element_text(size = 11, hjust=0.5),
        axis.title.y = element_text(size=11))