4 Curso de Econometría Aplicada en R

Parte I. Introducción a la Economettría

Curso 1. Regresión y naturaleza del análisis econométrico

Clase 1. Naturaleza del modelo Econometríco

Clase 2. Repaso de inferencia estadística

Clase 3. Modelo de regresión

Clase 4. Modelo de regresión múltiple

Clase 5. Modelo de regresión en R

Curso 2. Análisis del modelo de regresión y econometría aplicada

Clase 6. Test de hipótesis 

Clase 7. Inferencia sobre el modelo de regresión en R

Clase 8. Formas funcionales y variables discretas

Clase 9. Formas funcionales en R

Clase 10. Selección de modelos en R

Parte II. Microeconometría y datos de Panel

Curso 3. Microeconometría y datos de panal

Clase 1. Causalidad estadística

Clase 2. Endogeneidad y variables instrumentales

Clase 3. Modelos de probabilidad

Clase 4. Panel estatico 

Parte III. Econometría Financiera y de Series Temporales

Curso 4. Econometría financiera en R

Clase 1. Modelos ARIMA en R

Clase 2. Modelos de volatilidad (GARCH en R)

Clase 3. Modelos factoriales y riesgo sistémico 

Clase 4. Valor en riesgo (VaR)

Clase 5. Vectores autoregresivos (VAR)

Modelos GARCH sobre múltiples series financieras

 En la siguiente entrada vamos a estimar un modelo de volatilidad condicional para cada una de las variables contenidas en una base de datos, independientemente al numero de variables. Para tal ejemplo usamos la data EuStockMarkets que continente algunos índices financieros. Primero revisamos hemos importado los datos correctamente y realizamos un análisis descriptivo. En este caso no lo hacemos, pero recuerde siempre verificar temas como cambios estructurales, datos perdidos, distribución de los datos, ….

library(dplyr)
library(tibble)
library(fGarch)
library(ggplot2)
library(stargazer)
library(tidyr)

head(EuStockMarkets)
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

Posteriormente usamos la función across para obtener la variación de cada serie. Es importante precisar que un análisis formal requeriría un test de estacionariedad sobre las series.  

data_returns <- EuStockMarkets %>%
  as.tibble() %>%
  mutate(
    across(everything(),
           ~((./dplyr::lag(.)-1)*100)
    )) %>%
  na.omit()
 
head(data_returns)
# A tibble: 6 × 4
     DAX    SMI    CAC   F+TSE
   <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

Ahora, guardamos en una lista el modelo GARCH(1,1) correspondiente a cada una de las columnas de mi base de datos. Cada modelo se guarda en una lista dado que permite posteriormente recuperar todos los elementos de la salida de la función garchFit.

models0 <- list()
 
for (i in 1:ncol(data_returns)){
  models0[[i]] <- garchFit( ~  garch(1, 1), trace = F, data = data_returns[,i])
}
 
El commando stargazer nos permite obtener una tabla con los coeficientes asociado a cada variable.
stargazer(models0, type="text")
===========================================================
                              Dependent variable:         
                    ---------------------------------------
                        [         [         [         [   
                       (1)       (2)       (3)       (4)  
-----------------------------------------------------------
mu                  0.070***  0.106***   0.049**  0.052***
                     (0.021)   (0.020)   (0.025)   (0.017)
                                                          
omega               0.044***  0.120***   0.081**   0.009* 
                     (0.012)   (0.024)   (0.037)   (0.005)
                                                           
alpha1              0.068***  0.127***  0.051***  0.046***
                     (0.015)   (0.024)   (0.015)   (0.012)
                                                          
beta1               0.891***  0.735***  0.882***  0.941***
                     (0.023)   (0.043)   (0.042)   (0.018)
                                                          
-----------------------------------------------------------
Observations          1,859     1,859     1,859     1,859 
Log Likelihood      2,587.920 2,412.241 2,788.477 2,136.479
Akaike Inf. Crit.     2.789     2.600     3.004     2.303 
Bayesian Inf. Crit.   2.800     2.611     3.016     2.315 
===========================================================
Note:                           *p<0.1; **p<0.05; ***p<0.01

Otra opción es recuperar la volatilidad histórica asociada a cada modelo. Para esto utilizamos la función vectorizada lapply que nos permite repetir una operación/función, para cada uno de los elementos de una lista, generando una lista de vectores de volatilidades que convertimos en una matriz a partir de la función do.call.

 vols_models <- lapply(models0, function(x) x@sigma.t)
dataVol <- data.frame(do.call(cbind, vols_models))
colnames(dataVol) <- colnames(EuStockMarkets)
 
dataVol %>%  head()
 
        DAX       SMI      CAC      FTSE
1 1.0281714 0.9246996 1.102485 0.7967414
2 1.0265176 0.8843931 1.114000 0.7901069
3 1.0004172 0.8695389 1.167243 0.7807891
4 0.9915619 0.8258928 1.141444 0.7849713
5 0.9614558 0.7884808 1.125005 0.7754941
6 0.9420005 0.8386133 1.101589 0.7760277

 A partir de estas series de volatilidades, podemos obtener una representación gráfica de las series de volatilidades:

dataVol %>%
  mutate(fecha = 1:nrow(dataVol)) %>%
  as.data.frame() %>%
  pivot_longer(!fecha, names_to = "var", values_to = "value") %>%
     ggplot(aes(x = fecha, y = value)) +
     geom_line(aes(color = var))+
     facet_wrap(~var, scales='free_y', ncol = 2) +
     theme_classic()+
  theme(legend.position = "none",
        axis.ticks = element_line(colour = "grey70", size = 0.2),
        strip.background = element_blank()) +
  scale_colour_grey(start = 0.1,end = 0.1)

Estudio de eventos: representar promedios de series temporales en ggplot2

En la siguiente entrada se utiliza la base de datos tipo wide de nombre EuStockMarkets, disponible en R, la misma cuenta con los índices asociados a los siguientes activos financieros (solo se presenta como obtener los datos sin profundizar sobre la forma de interpretación de los mismos):

library(tidyverse)
head(EuStockMarkets)
head(EuStockMarkets)
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

Transformamos la data a versión wide, porque al hacer un by_group.

wide_data <- data.frame(EuStockMarkets) %>%
  mutate(fecha = time(EuStockMarkets)) %>%
  gather(id, value, -fecha)
 
wide_data %>%  head()
 
> wide_data %>%  head()
     fecha  id   value
1 1991.496 DAX 1628.75
2 1991.500 DAX 1613.63
3 1991.504 DAX 1606.51
4 1991.508 DAX 1621.04
5 1991.512 DAX 1618.16
6 1991.515 DAX 1610.61

Gráfico resaltando un evento determinado:

wide_data %>%
  ggplot(aes(x = fecha, y = value)) +
  geom_line(aes(color = id), size = 1)+
  theme_minimal() +
  facet_wrap(~id, scales='free_y',ncol=1)+
  theme(legend.position = "none")+
  geom_vline(xintercept=1997,linetype=4)

Agregando promedio alrededor de un evento

wide_data %>%
  dplyr::mutate(dummyT = (as.numeric(fecha)>1997)*1) %>%
  group_by(id,dummyT) %>%
  mutate(medias = mean(value, na.rm = T)) %>%
  ggplot(aes(x=fecha,y=value))+
  geom_line(aes(color = id))+
  geom_line(aes(x=fecha,y=medias,color = id))+
  facet_wrap(~id, scales='free_y',ncol=1)+
  theme(legend.position = "none")+
  geom_vline(xintercept=1997,linetype=4)+
  theme_minimal()+
  theme(legend.position = "none")

Agregando promedio alrededor de varios eventos

wide_data %>%
  dplyr::mutate(dummyT = case_when(
    as.numeric(fecha)<1994 ~ 1,
    between(as.numeric(fecha),1994,1997) ~ 2,
    as.numeric(fecha)>1997 ~ 3
  )) %>%
  group_by(id,dummyT) %>%
  mutate(medias = mean(value, na.rm = T)) %>%
  ggplot(aes(x=fecha,y=value))+
  geom_line(aes(color = id))+
  geom_line(aes(x=fecha,y=medias,color = id))+
  facet_wrap(~id, scales='free_y',ncol=1)+
  theme(legend.position = "none")+
  geom_vline(xintercept=c(1994,1997),linetype=4)+
  theme_minimal()+
  theme(legend.position = "none")

Exportar gráficos al directorio de trabajo en R

En el siguiente ejemplo se muestra como exportar una lista de gráficos al directorio de trabajo, guardado cada uno con un nombre en específico. Primero se guardan los gráficos en una lista; posterior se crea una lista de nombres con los que se guardaran los gráficos en el directorio de trabajo; posteriormente se asigna a cada nombre un gráfico determinado; por último, se guardan estos gráficos con ggsave.

library(tidyverse)
 
cesar <- list()
cesar[[1]] <- ggplot2::ggplot(iris, aes(x=Sepal.Length))+geom_density()
cesar[[2]] <- ggplot2::ggplot(iris, aes(x=Sepal.Length,y=Petal.Length))+geom_point()
 
nombres <- paste("graph", 1:length(cesar), sep="")
 
for (i in 1:length(cesar)){
  assign(nombres[i], cesar[[i]])
 
  nombres2 <- paste(nombres[i], ".png", sep="")
   ggplot2::ggsave(nombres2)
}

Modelos ARIMA en R. Tips rápidos para su estimación y pronósticos

Para el siguiente ejemplo utilizamos el índice de alimentos disponible en la FRED, con este se intentan explicar algunos tips básicos para la estimación e interpretación de los modelos ARIMA en R. En las siguientes líneas mostramos la estructura de los datos, similar al formato tradicional en Excel, una columna con fechas y otras con la serie de datos asociado al índice con el que vamos a trabajar.

library(tidyverse)
library(readxl)
library(xts)
 
datos <- read_excel("~/food_index_fred.xls",
                              col_types = c("date", "numeric"))
 
head(datos)
 
# A tibble: 6 x 2
  observation_date    PFOODINDEXM
  <dttm>                    <dbl>
1 1992-01-01 00:00:00        57.7
2 1992-02-01 00:00:00        57.8
3 1992-03-01 00:00:00        57.5
4 1992-04-01 00:00:00        56.2
5 1992-05-01 00:00:00        57.7
6 1992-06-01 00:00:00        57.3

Vamos la evolución histórica de la serie. Declaramos que la serie es un objeto temporal a partir del paquete xts y usando el índice de precios. Esta definición es especialmente útil cuando trabajamos con datos diarios o donde la secuencia de fecha asociada no es x-distantes en el tiempo.

# asignando representation de una serie temporal
ts_food <- xts(x = datos$PFOODINDEXM, order.by = datos$observation_date)
 
autoplot(ts_food) +
  theme_minimal() +
  labs(y = "Global price of Food index", x = "Tiempo")

Estimamos el correlograma. Si los coeficientes están dentro del intervalo de confianza, no rechazamos la hipótesis nula de que estos sean iguales a cero. En caso contrario rechazamos. Claramente quedan fuera, por lo que, no rechazamos la hipótesis nula. Cuando rechazamos (están fuera del intervalo) decimos que la serie presenta autocorrelación.

Ahora obtenemos una representación de las series en niveles y sus diferencias (se transforma en logaritmo porque ayuna a cumplir supuestos como estabilidad de varianza, y permite que la diferenciación se puede interpretar como una tasa de variación que no depende de las unidades de escalas).

# Obtener cambio logaritmico
dlnts_food <- diff(log(ts_food))[-1]
 
par(mfrow=c(2,2))
plot(ts_food)
acf(ts_food)
 
plot(dlnts_food)
acf(dlnts_food)

Estimamos la estrategia ARIMA. En el siguiente ejemplo se estima un ARIMA para una serie integrada de primer orden (suponemos solo se necesita una diferenciación para obtener una serie estacionaria) de orden 1,1, para los coeficientes AR y MA, respectivamente. El primer número corresponde al AR, el segundo al orden de integración de la serie y el tercero al MA.  

 # estimacion ARIMA
model1 <- arima(ts_food,order=c(1,1,1)
model1
 
Call:
arima(x = ts_food, order = c(1, 1, 1))
 
Coefficients:
        ar1    ma1
      0.313  0.058
s.e.  0.174  0.188
 
sigma^2 estimated as 6.69:  log likelihood = -860,  aic = 1726

Evaluación del modelo. La evaluación se hace recuperando los residuos y verificando estos no tienen autocorrelación o no están asociado con ninguna otra variable. En caso contrario podría apuntar a que se omite información o que no se está modelando correctamente la estructura de la serie. Note en el siguiente correlograma que todos los valores están dentro, por tanto, no rechazamos la hipótesis nula de que estos sean iguales a cero. Aceptamos el modelo. En caso de que estos estén fuera apunta a que necesitamos otra alternativa para modelarlo.

acf(r)

Una alternativa de evaluación de la autocorrelación de las series se obtiene del test de Box. Cuya hipótesis nula es que estos coeficientes son iguales a cero. En ambos casos se obtienen p-valor elevados (mayor al valor critico seleccionado), esto hace que no se rechace la hipótesis nula asociada al test, por tanto, decimo que la serie no presenta autocorrelación.

Box.test(r,lag=1,type = "Box-Pierce")
       Box-Pierce test
 
data:  r
X-squared = 0.011, df = 1, p-value = 0.9
 
Box.test(r,lag=1,type = "Ljung-Box")
       Box-Ljung test
 
data:  r
X-squared = 0.011, df = 1, p-value = 0.9

Uso del modelo. El modelo se puede usar para predecir o descomponer la evolución de la serie en una parte sistemática y otra no (clave en la gestión del riesgo). Predict nos permite recuperar la predicción y el error estándar de dicho error.

predict(model1,n.ahead=12)
$pred
Time Series:
Start = 31449601
End = 32400001
Frequency = 0.0000115740740740741
 [1] 160.9 161.1 161.2 161.2 161.2 161.2 161.2 161.2 161.2 161.2 161.2 161.2
 
$se
Time Series:
Start = 31449601
End = 32400001
Frequency = 0.0000115740740740741
 [1]  2.586  4.389  5.836  7.042  8.084  9.010  9.851 10.625 11.347 12.026
[11] 12.669 13.280

Finalmente, recuerde que la función auto.arima de R, permite identificar este tipo de modelos.

# Auto.Arima
library(forecast)
model2 <-auto.arima(ts_food, seasonal=F, ic ="aic", trace=T)
 
 
 Fitting models using approximations to speed things up...
 
 ARIMA(2,1,2) with drift         : 1729
 ARIMA(0,1,0) with drift         : 1770
 ARIMA(1,1,0) with drift         : 1723
 ARIMA(0,1,1) with drift         : 1724
 ARIMA(0,1,0)                    : 1772
 ARIMA(2,1,0) with drift         : 1726
 ARIMA(1,1,1) with drift         : 1725
 ARIMA(2,1,1) with drift         : 1727
 ARIMA(1,1,0)                    : 1723
 ARIMA(2,1,0)                    : 1726
 ARIMA(1,1,1)                    : 1725
 ARIMA(0,1,1)                    : 1725
 ARIMA(2,1,1)                    : 1727
 
 Now re-fitting the best model(s) without approximations...
 
 ARIMA(1,1,0)                    : 1724
 
 Best model: ARIMA(1,1,0)

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