Tidyverse y minería de datos en R: estimar regresiones por grupos

La siguiente entrada utiliza la potente librería tidyverse, de R, para estimar regresiones por grupo. Para el siguiente ejemplo se utiliza la base  bwght del libro de wooldridge.

library(tidyverse); library(purrr)

library(wooldridge)

attach(bwght)

En primer lugar, indicamos la base de datos con la que deseamos trabajar (bwght %>%); posteriormente, le indicamos para cuales grupos específicos requerimos estimar nuestra regresiones (split(.$white) %>%). En este caso, se desea repetir la estimación para cada grupo, según color de la piel; finalmente, especificamos la regresión que deseamos estimar (map(~ lm(lbwght ~ cigprice+faminc, data = .)) %>%) ayudados por la funciones map y lm; y llamamos la tradicional salida de la regresión de R (map(summary)). Note, que la regresión propuesta se estima para cada una de las categorías indicadas en la variable White.

bwght %>%

  split(.$white) %>%

  map(~ lm(lbwght ~ cigprice+faminc, data = .)) %>%

  map(summary)

$`0`

Call:

lm(formula = lbwght ~ cigprice + faminc, data = .)

Residuals:

     Min       1Q   Median       3Q      Max

-1.57471 -0.11323  0.02753  0.15088  0.48549

Coefficients:

             Estimate Std. Error t value Pr(>|t|)   

(Intercept) 4.522e+00  1.706e-01  26.503   <2e-16 ***

cigprice    1.483e-03  1.327e-03   1.117    0.265   

faminc      2.742e-05  8.369e-04   0.033    0.974   

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2212 on 296 degrees of freedom

Multiple R-squared:  0.004233, Adjusted R-squared:  -0.002495

F-statistic: 0.6292 on 2 and 296 DF,  p-value: 0.5337

$`1`

Call:

lm(formula = lbwght ~ cigprice + faminc, data = .)

Residuals:

     Min       1Q   Median       3Q      Max

-1.21048 -0.08434  0.01967  0.11253  0.81989

Coefficients:

             Estimate Std. Error t value Pr(>|t|)   

(Intercept) 4.7104582  0.0692942  67.978  < 2e-16 ***

cigprice    0.0002696  0.0005266   0.512  0.60876   

faminc      0.0008468  0.0002940   2.880  0.00406 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1789 on 1086 degrees of freedom

Multiple R-squared:  0.00803,  Adjusted R-squared:  0.006204

F-statistic: 4.396 on 2 and 1086 DF,  p-value: 0.01255

En el segundo ejemplo, únicamente se agrega la función map_dbl, para llamar un estadístico específico de los calculados por R de forma automática. En este ejemplo se llama el R cuadrado (map_dbl("r.squared")).  

# Ejemplo: obtiene el R2 para cada grupo

bwght %>%

  split(.$white) %>%

  map(~ lm(lbwght ~ cigprice+faminc, data = .)) %>%

  map(summary) %>%

  map_dbl("r.squared")

En caso de desear ver cómo saber cuáles estadísticos podemos llamar:

models1<-lm(lbwght ~ cigprice+faminc)

summary(models1)

ames(summary(models1))

> names( summary(models1))

 [1] "call"          "terms"         "residuals"     "coefficients"  "aliased"       "sigma"        

 [7] "df"            "r.squared"     "adj.r.squared" "fstatistic"    "cov.unscaled" 

Introducción a la econometría de Wooldridge: ejercicios resueltos en R. Capítulo 1

C1.1. Para este ejercicio emplee la base de datos WAGE1.RAW.

> install.packages(“wooldridge”)

> library(wooldridge)

> attach(wage1)

> ??wage1

i) Determine el nivel educativo promedio de la muestra. ¿Cuáles son los niveles de educación menor y mayor?

> mean(educ)

[1] 12.56274

> summary(educ)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

   0.00   12.00   12.00   12.56   14.00   18.00

ii) Determine el salario promedio por hora (wage) en la muestra. ¿Parece ser alto o bajo?

> mean(wage)

[1] 5.896103

iii) Los datos de los salarios están dados en dólares de 1976. Usando el Economic Report of the President (de 2004 o posterior o el Informe de Gobierno en países de habla hispana) obtenga y dé los índices de precios al consumidor (IPC) correspondientes a 1976 y 2003.

iv) Use los valores de los IPC del inciso iii) para determinar el salario promedio por hora en dólares de 2003. ¿Parece ahora más razonable el salario promedio por hora?

# Buscar el IPC correspondiente a cada uno de los años del problema y actualizar el salario.

v) ¿Cuántas mujeres (females) hay en la muestra? ¿Cuántos hombres?

> #mujeres

> sum(female)

[1] 252

>

> #hombres

> length(female)-sum(female)

[1] 274

C1.2. Para responder estas preguntas emplee la base de datos BWGHT.RAW.

i) ¿Cuántas mujeres hay en la muestra (male = 0) y cuántas de las informantes fumaron durante un embarazo?

> table(male)

male

  0   1

665 723

Pero, dado que todas las observaciones son mujeres, la pregunta sería cuantas mujeres tuvieron hijas. Para saber el total de observaciones en la muestra:

> length(male)

[1] 1388

> length(male[cigs>=1])

[1] 212

ii) ¿Cuál es la cantidad promedio de cigarros consumidos por día (cigs)? ¿Es el promedio,

en este caso, una medida representativa de la mujer “típica”? Explique.

> mean(cigs)

[1] 2.087176

Estudiar la asimetría de la serie, así como la presencia de valores atípicos influyentes, ayudan a identificar la validez o no de la media como medida representativa de las series.

> summary(cigs[male==0])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

   0.00    0.00    0.00    2.09    0.00   40.00

> hist(cigs)

Aunque, independientemente a las propiedades estadísticas de las series, se puede verificar que el 84% de las mujeres no consumieron cigarrillos durante el embarazo. Por lo que, el promedio no sería un indicador adecuado.

> mean(cigs==0)

[1] 0.8472622

> table(cigs)

cigs

   0    1    2    3    4    5    6    7   

1176    3    4    7    9   19    6    4    ...

iii) Entre las mujeres que fumaron durante el embarazo, ¿cuál es la cantidad promedio de

cigarros consumidos por día? ¿Cuál es la relación de esto con su respuesta al inciso ii) y

por qué?

> mean(cigs[cigs>=1])

[1] 13.66509

iv) Determine el promedio de fatheduc (años de educación del padre) en la muestra. ¿Por

qué se emplean sólo 1,192 observaciones para calcular este promedio?

Para calcular los años promedio de educación del padre, es necesario indicarle al programa que el vector tiene valores perdidos (na.rm = TRUE).

> mean(fatheduc, na.rm = TRUE)

[1] 13.18624

Tabulando los resultados, se verifica que se utiliza esta cantidad de observaciones porque son las únicas disponibles. Las demás corresponden a valores perdidos.

> table(is.na(fatheduc))

FALSE  TRUE

 1192   196

v) Dé el ingreso familiar promedio (famine) y su desviación estándar en dólares.

> mean(faminc)

[1] 29.02666

> sd(faminc)

[1] 18.73928

C1.3. Los datos de MEAP01.RAW pertenecen al estado de Michigan en el año 2001. Emplee estos datos para contestar las preguntas siguientes.

> #Borra todos los objetos en memoria

> rm(list=ls())

> attach(meap01)

i) Determine los valores mayor y menor de math4. ¿Es lógico este intervalo? Explique.

> ??meap01

> # math4: percent students satisfactory, 4th grade math

> summary(math4)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

   0.00   61.60   76.40   71.91   87.00  100.00

ii) ¿Cuántas escuelas tienen una tasa perfecta de aprobados en el examen de matemáticas?

¿A qué porcentaje del total de la muestra corresponde esta cantidad?

> # Escuela con porcentaje completos de aprobados

> sum(math4==100)

[1] 38

> # Porcentaje de escualas

> mean(math4==100)

[1] 0.02084476

iii) ¿En cuántas escuelas la tasa de aprobados en matemáticas es exactamente 50%?

> sum(math4==50)

[1] 17

iv) Compare el promedio de las tasas de aprobados en matemáticas y en lectura. ¿Cuál de estas pruebas es más difícil de aprobar?

> apply(cbind(math4, read4), 2, mean)

   math4    read4

71.90900 60.06188

v) Encuentre la correlación entre math4 y read4. ¿Qué concluye?

> cor(math4, read4)

[1] 0.8427281

vi) La variable exppp es gasto por alumno. Determine el promedio y la desviación estándar

de exppp. ¿Parece haber una gran variación en el gasto por alumno?

> coef_var<-sd(exppp)/mean(exppp)

> coef_var

[1] 0.2101863

vii) Suponga que la escuela A gasta 6 000 dólares por alumno y la escuela B gasta 5 000 dólares por alumno. Dé el porcentaje en el que el gasto de la escuela A supera al gasto de la escuela B. Compare este porcentaje con 100 · [log(6 000) – log(5 500)], que es la diferencia porcentual aproximada basada en la diferencia de los logaritmos naturales. (Veáse la sección A.4 del apéndice A.)

> cat("Lo supera en \n",((6000-5000)/5000)*100, "\n")

Lo supera en

 20

> 100*(log(6000)-log(5000))

[1] 18.23216

C1.4. La base de datos de JTRAIN2.RAW proviene de un experimento de capacitación para el trabajo realizado para hombres con bajos ingresos durante 1976-1977; véase Lalonde (1986).

> rm(list=ls())

> attach(jtrain2)

i) Emplee la variable indicadora train para determinar la proporción de hombres a los que se les dio capacitación para el trabajo.

> mean(train)

[1] 0.4157303

ii) La variable re78 es ingresos desde 1978, dados en dólares de 1982. Determine el promedio de re78 para la muestra de hombres a los que se les dio capacitación laboral y para la muestra de hombres a los que no se les dio. ¿Es esta diferencia económicamente grande?

> tapply(re78, train, mean)

       0        1

4.554802 6.349145

Hacer un test t, de comparación de medias. Claramente, a partir del p-valor, se rechaza la hipótesis nula de que las medias sean iguales.

> t.test(re78[train==1],re78[train==0])

     Welch Two Sample t-test

data:  re78[train == 1] and re78[train == 0]

t = 2.6741, df = 307.13, p-value = 0.007893

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

 0.4740108 3.1146754

sample estimates:

mean of x mean of y

 6.349145  4.554802

iii) La variable unem78 indica si un hombre estuvo desempleado o no en 1978. ¿Qué proporción de los hombres a los que se les dio capacitación para el trabajo están desempleados? ¿Y de aquellos a los que no se les dio capacitación laboral? Comente la diferencia.

> tapply(unem78, train, mean)

        0         1

0.3538462 0.2432432

iv) Con base en los incisos ii) y iii), ¿parece haber sido efectivo el programa de capacitación laboral? ¿Qué haría que nuestra conclusión fuera más convincente?

Basados en que las diferencias entre promedios de ingresos y menor proporción de desempleo son estadísticamente significativas, es importante destacar que estas comparaciones no han considerado la evolución de otros factores externos que pudiesen influir en los resultados y que deberían ser considerados, además de poder garantizar la equivalencia entre los grupos, para que estas comparaciones, antes/después, puedan ser válidas.

Usar el bucle for para guardar resultados de la estimación de una regresión

En la siguiente entrada, se muestra cómo utilizar el buche for, para guardar resultados estimados en una regresión lineal simple (lm), sea en un vector o en una matriz. En el ejemplo siguiente, se utiliza la base de datos iris, que contiene las siguientes variables:

head(iris)

  Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1          5.1         3.5          1.4         0.2  setosa

2          4.9         3.0          1.4         0.2  setosa

3          4.7         3.2          1.3         0.2  setosa

4          4.6         3.1          1.5         0.2  setosa

5          5.0         3.6          1.4         0.2  setosa

6          5.4         3.9          1.7         0.4  setosa

En el ejemplo, se estima mediante lm, un modelo de regresión simple de la primera variable (iris[,1] == Sepal.Length) en función de la segunda, tercera y cuarta variable, respectivamente (iris[,i] == Sepal.Length). Véase de forma directa, como este ejemplo se puede extender a cualquier número de variables, o inclusive utilizar indexación para identificar de forma automática las variables de un data.frame que cumple con alguna condición.

# --------------------------------

head(iris)

for (i in 2:4){

  mod<-lm(iris[,1]~iris[,i])

}

summary(mod)

O de forma alternativa, se puede utilizar la opción get() para llamar variables directamente, a partir de un vector de nombres:

Sepal.Length = b0 + b1*Sepal.Width

Sepal.Length = b0 + b1*Petal.Length

Sepal.Length = b0 + b1*Petal.Width

# --------------------------------

attach(iris)

var<-names(iris)[2:4]

> var

[1] "Sepal.Width"  "Petal.Length" "Petal.Width" 

for (i in var){

  mod<-lm(Sepal.Length~get(i))

}

El problema con las secuencias anteriores, es que pese a realizar un número determinado de regresiones, los resultados estimados no son guardados. Por lo que, en los siguientes ejemplos se muestra cómo se pueden guardar los resultados. Note como se crea un vector vacío (sumResid) que posteriormente se va completando con la suma de los residuos (mod$residuals) al cuadrado que se obtiene en la estimación de cada regresión de manera sucesiva [c(sumResid, sum(mod$residuals))].

# Guardar la suma de los residuos de cada regresión

# --------------------------------

sumResid <-c()

for (i in var){

  mod<-lm(Sepal.Length~get(i))

  sumResid <-c(sumResid, sum(mod$residuals)) 

}

> sumResid

[1] -4.982126e-15  9.992007e-16 -5.516421e-16

En un ejemplo final, en vez de guardar un vector 1x1 en cada iteración del bucle, es posible guardar un vector de coeficientes estimados en cada una de las regresiones (mod$coefficients). Ahora, en vez de ir anexando un escalar a un vector determinado, como hicimos en el ejemplo anterior, se va agregando un vector de coeficientes (coeff), a la matriz de coeficientes estimados que se va completando, utilizando el comando rbind.

# Guardar los coeficientes

# --------------------------------

coeff <-c()

for (i in var){

  mod<-lm(Sepal.Length~get(i))

  coeff <-rbind(coeff, mod$coefficients) 

}

> coeff

     (Intercept)     get(i)

[1,]    6.526223 -0.2233611

[2,]    4.306603  0.4089223

[3,]    4.777629  0.8885803

Fundamentos de Análisis de Series Temporales en R

La siguiente entrada corresponde a una presentación  sobre series temporales en R, realizada en la comunidad Data Science | DR, en marzo 14 de 2019, sobre técnicas de análisis de series temporales, con aplicaciones en R [DESCARGAR DOCUMENTO].

Contenido del documento:

1. Series temporales: aspectos básicos: procesos estocásticos; raíces unitarias, estacionariedad, operador de retardos, componentes de  series temporales; proceso ruido blanco; AR; caminata aleatoria; MA; ARMA; ARIMA; SARIMA; RegARIMA.

2. Box-Jenkins: Box-cox; Dickey-Fuller; Phillip-Perron; Box-Pierce; Ljiung-Box, pronósticos; análisis residual; correlogramas; funciones de autocorrelación; criterios de información; estimación.

3. Modelos de volatilidad: Volatilidad condicional, efectos ARCH; EWMA; ventanas de volatilidad; GARCH gaussianos; ARMA-GARCH; GARCH asimetricos (gjrGARCH); Markov-GARCH; modelos de volatilidad estocástica.

4. Análisis multivariado: Vectores autoregresivos (VAR); Vectores de corrección de errores (VEC); Cointegración; Johansen; Causalidad; Funciones de respuesta al impulso; descomposición de varianza; modelos factoriales; VAR estructurales (SVAR).

5. Modelos multivariados de volatilidad: VEC;  CCC; DCC; BEKK; transferencia en volatilidad; ventanas de correlaciones.

6. Métodos bayesianos: modelos lineales dinámicos.

Referencias

Modelos GARCH multivariados en R

1.      DCC (igual especificación de la varianza)

En las siguientes instrucciones se simulant dos series (x,y), solo para mostrar el uso del paquete rmgarch, combinado con rugarch y parallel, para la estimación del modelo Dynamic Conditional Correlation (DCC).

# ------------------------------- *
# Activa paquetes
library(rmgarch)
library(rugarch)
library(parallel)

# Simula las series
x<-rnorm(1000,1,3)*c(1:1000)+c(1:1000)
y<-rnorm(1000,1,3)*c(1:1000)

data1<-as.matrix(cbind(x,y))
plot(x, type="l")

En las siguientes líneas, se introduce la estructura del modelo a estimar, se verifica un modelo gjrGARCH con distribución t, para cada una de las series consideradas.

# ------------------------------- *
l <- 2 #número de variables
gjrtspec <- ugarchspec(
         mean.model=list(armaOrder=c(0,0)),
   variance.model =list(model = "gjrGARCH"),
          distribution="std")

dcc_spec = dccspec(
              uspec = multispec(replicate(l, gjrtspec)),
              distribution = "mvt")

# Fit DCC
garchdccfit = dccfit(dcc_spec, data1, fit.control=list(scale=TRUE))

Luego de estimar el modelo, podemos acceder a la lista de resultados y obtener los parámetros asociados a los modelos de volatilidad inivariados, información sobre la especificación estimada y los criterios de información.

*---------------------------------*

*          DCC GARCH Fit          *

*---------------------------------*

Optimal Parameters

-----------------------------------

                Estimate  Std. Error     t value Pr(>|t|)

[x].mu         88.018073  4.0934e-01  2.1503e+02 0.000000

[x].ar1         0.031399  4.3700e-04  7.1876e+01 0.000000

[x].omega       0.004366  2.1990e-03  1.9852e+00 0.047120

[x].alpha1      0.054871  1.9430e-03  2.8239e+01 0.000000

[x].beta1       0.945898  0.0000e+00  5.7637e+06 0.000000

[x].gamma1     -0.088033  1.3151e-01 -6.6939e-01 0.503249

[x].delta       0.091525  1.5678e-02  5.8378e+00 0.000000

[y].mu         81.424172  2.1262e+01  3.8295e+00 0.000128

[y].ar1         0.060300  3.7522e-02  1.6071e+00 0.108038

[y].ar2         0.021913  3.8664e-02  5.6675e-01 0.570888

[y].omega    1571.508757  1.0843e+03  1.4494e+00 0.147235

[y].alpha1      0.141311  2.5299e-02  5.5857e+00 0.000000

[y].beta1       0.870438  2.0259e-02  4.2965e+01 0.000000

[y].gamma1     -0.025498  3.5864e-02 -7.1094e-01 0.477120

[Joint]dcca1    0.002363  9.5070e-03  2.4860e-01 0.803672

[Joint]dccb1    0.956960  1.9858e-01  4.8191e+00 0.000001

A continuación se muestra como recuperar la serie condicional de matrices de varianzas covarianzas.

# ------------------------------- *
# Matriz Var-Cov
dcccov<-rcov(garchdccfit)

# Matrix de Corr
dcccor<-rcor(garchdccfit)

Ahora, en el siguiente bloque de instrucciones, se utiliza un blucle e indexación numérica de Arrays, para recuperar la serie de volatilidad condicional, guardada en Sigma.

rcov(garchdccfit)[,,1]
x y
x 484311.21 13967.32
y 13967.32 3262974.28

En este aspectos, como obtenemos T matrices, cada una de las cuales guarda el valor modelado de los momentos para el día correspondiente, [1,1,] corresponde a la varianza del primer a la varianza de la primera serie; por lo que, mediante indexación de arrays podemos acceder a los componentes de las matrices de varianzas covarinzas y de correlaciones que hemos guardado anteriormente.

plot(dcccov[1,1,], type="l")

Una alternativa a la anterior, es utilizar un bucle un contador i, que recorra las T observaciones, correspondiente a la dimensión de las variables en cuestión. Accediendo a la posición de la matriz que resulte de interés.

# ------------------------------- *
# Recuperar, serie temporal a partir de la estructura array
vol1<-c()
T<-1000

for (i in 1:T){
  vol1<-c(vol1, dcccov[1,1,i])
}

plot(vol1, type="l")

2.      DCC (diferentes especificaciones de la varianza)

En el ejemplo anterior se consideraba que ambas series seguían el mismo proceso de varianzas, cuestión que puede ser cuestionable en muchos procedimientos empíricos, por lo que, en el próximo ejemplo, se muestra cómo se pueden realizar especificaciones de varianzas distintas para cada una de las n-variables consideradas.

uspec1 <- ugarchspec(  
             mean.model = list(armaOrder = c(1,0)),
        variance.model = list(model = "apARCH"),
   distribution.model = "norm")

uspec2 = ugarchspec(
            mean.model = list(armaOrder = c(2,0)),
       variance.model = list(model = "gjrGARCH"),
  distribution.model = "norm")

uspec = c(uspec1, uspec2)
   

 spec = dccspec(uspec  = multispec( uspec ),
                     dccOrder = c(1,1),
                 distribution = "mvlaplace")

garchdccfit = dccfit(spec, d ata1,  fit.control=list(scale=TRUE))

garchdccfit

3.      MGARCH-BEKK

El paquete mgarchBEKK permite la estimación del modelo mGARCH-BEKK, para las matrices C, A y B respectivamente. Por tanto, como B recoge el efecto hoy de la matriz sigma en el periodo anterior, los coeficientes fuera de su diagonal principal brindan una aproximación de la interacción entre las volatilidades del mercado. Los estadísticos de prueba de significación se pueden obtener utilizando los errores asintóticos, guardados en:

# ------------------------------- *
library(mgarchBEKK)
estimated <- BEKK(data1)

estimated$est.params
$`1`
[,1] [,2]
[1,] 3.663613 -8.176183
[2,] 0.000000 24.878929

$`2`
[,1] [,2]
[1,] 0.1539364 0.09550136
[2,] -0.1657275 0.09612336

$`3`
[,1] [,2]
[1,] -0.1222077 0.9859749
[2,] 0.9703889 0.1276075

#Errores de los coeficientes
estimated$asy.se.coef

Los momentos condicionales de la relación entre las variables (estimated$cor), así como las series de volatilidad condicional (estimated$sd), se pueden recuperar:

# ------------------------------- *
estimated$cor[[1]][[2]]
estimated$sd[[1]]

estimated$uncond.cov.matrix

estimated$aic
estimated$order

4.      gjrBEKK GARCH

El paquete mgarchBEKK también permite la estimación del modelo BEKK asimétrico, que incluye una variables binarias (1 cuando el residuo es negativo) que interactúan con el residuo correspondiente de forma tal que permiten considerar un efecto diferenciado del retorno asociado, dependiendo del signo del mismo. En el caso de los coeficientes (estimated1$est.params), verifique que aparece una nueva matriz, que corresponde a la matriz cuadradada asociada a los diferenciales de pendiente.  

# ------------------------------- *
estimated1 <- mGJR(x,y)

# Parámetros y errores asintóticos
estimated1$est.params
$`1`
[,1] [,2]
[1,] 2.335441 -2.395821
[2,] 0.000000 30.064757

$`2`
[,1] [,2]
[1,] 0.1503465 -0.01745010
[2,] -0.1305680 -0.03214134

$`3`
[,1] [,2]
[1,] 1.0032073 0.1680627
[2,] -0.2073904 0.9359224

$`4`
[,1] [,2]
[1,] 0.1016553 0.2527250
[2,] 0.1800089 -0.1394926

$`5`
[1] -0.3211898

estimated1$asy.se.coef

Las correlaciones y las varianzas, tantos condicionales como incondicionales, se pueden obtener a partir de la siguiente representación:

# ------------------------------- *
estimated1$cor
estimated1$sd1
estimated1$sd2
estimated1$uncond.cov.matrix
estimated1$resid1

Referencias

·        Ghalanos, Alexios (2019). Multivariate GARCH Models. Package ‘rmgarch’. Repository CRAN. Date/Publication 2019-01-15 05:40:03 UTC. https://cran.r-project.org/web/packages/rmgarch/rmgarch.pdf.  

·        Ghalanos, Alexios (2019). The rmgarch models: background and properties. (version 1.3.0). https://cran.r-project.org/web/packages/rmgarch/vignettes/The_rmgarch_models.pdf.

·        Orskaug, Elisabeth (2009). Multivariate DCC-GARCH Model -With Various Error Distributions. Norwegian University Of Science and Technology. Department of Mathematical Sciences.

·        Schmidbauer, Harald; Roesch, Angi; and Tunalioglu, Vehbi (2016). Simulating, Estimating and Diagnosing MGARCH (BEKK and mGJR) Processes. Package ‘mgarchBEKK’. Date/Publication 2016-04-10 00:49:46. https://cran.r-roject.org/web/packages/mgarchBEKK/mgarchBEKK.pdf.