Economía Aplicada
Notas casi aleatorias sobre economía, programas estadísticos y econometría
27 jul 2022
4 Curso de Econometría Aplicada en R
1 jul 2022
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(tibble)
library(fGarch)
library(ggplot2)
library(stargazer)
library(tidyr)
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.
as.tibble() %>%
mutate(
across(everything(),
~((./dplyr::lag(.)-1)*100)
)) %>%
na.omit()
# 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[[i]] <- garchFit( ~ garch(1, 1), trace = F, data = data_returns[,i])
}
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.
dataVol <- data.frame(do.call(cbind, vols_models))
colnames(dataVol) <- colnames(EuStockMarkets)
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
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)
22 jun 2022
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")
2 jun 2022
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)
}
27 may 2022
Modelos ARIMA en R. Tips rápidos para su estimación y pronósticos
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)
28 abr 2022
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:
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
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
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
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
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
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í...
-
Missing Values (NA) in R. En la siguiente entrada se muestran algunas operaciones básicas para la identificación y tratamiento de valor...
-
Las series económicas pueden estar influidas por una serie de procesos no determinista, ni conocidos para el analista y que pueden incidir ...
-
Tablas de frecuencia con condicionales: tabulate Las tablas de frecuencia absoluta nos permiten contabilizar casos que cumplen con d...