Para trabajar este taller usaremos la base de datos conejomalo que habíamos revisado en la monitoría. En primer lugar construimos las matrices x e y

# leemos la base de excel
library(readxl)
conejomalo <- read_excel("C:/Users/teoli/Documents/ECONOMETRIA_AVANZADA_2020/Monitorias R/conejomalo.xlsx")
View(conejomalo)
# construimos y
y<-conejomalo$Rep
# constuimos x
x<-cbind(conejomalo$danceability,conejomalo$energy,conejomalo$loudness,conejomalo$speechiness,conejomalo$valence)
# como se pierden los nombres se los volvemos a colocar
colnames(x)<-c("danceability","energy","loudness","speechiness","valence")

Ejercicio 1

Realice una rutina, function, con nombre MiRegresion en R tal que

  1. Entre
  1. Salga
MiRegresion<-function(Y,X,int=T,anova=T) {
  # se puede especificar en los argumentos los valores por defecto de verdadero (T) para el intercepto y la tabla de análisis de varianza, en caso de no querer incluirlos se especificará FALSE al momento de llamar la función.
  # Como X es una matriz, el numero de observaciones será insuficiente si el numero de parametros a estimar es mayor, esto es, si no hay grados de libertad. El numero de parametros será igual al número de columnas de X, por tanto
  if(nrow(X)<ncol(X)){stop("El número de observaciones es insuficiente")
  }
  # tambien es valido usar print para enviar el mensaje, sin embargo, stop detiene la ejecución controlando la generación de errores
  # la matriz X'X no es invertible si su determinante es cero
  if(det(t(X)%*%X)==0) {stop("la matriz X'X no es invertible")
  }
  # En caso de haber incluido constante se agrega una columna de unos a la matriz X, esto se puede hacer empleando cbind, para combinar un vector de unos con la matriz de regresores X, el vector de 1 se realiza con rep, función que replica elementos de un vector, en este caso se le dice que replique un 1 el número de veces que tiene de observaciones la matriz X. Finalmente se le pone un nombre con colnames:
  if(int==T){X <- cbind(rep(1, nrow(X)), X)
  colnames(X)[1] <- "Const"
  }
  # El vector de betas estimados se obtiene facilmente con:
  b<-(solve(t(X)%*%X)%*%t(X)%*%Y)
  print("betas estimados")
  print(b)
  # hallamos los errores del modelo
  e<-(Y-(X%*%b))
  # Calculamos la varianza del error
  s2<-(t(e)%*%e)/(nrow(X)-ncol(X))
  print("Varianza del error")
  print(s2)
  # Con ella el vector de errores estándar. se debe añadir la operacion diag porque cuando se trabaja con más de una variable regresora, el resultado de X'X no es un número sino una matriz.
  SE<-sqrt(diag(solve(t(X)%*%X))%*%s2)
  # Usar la diagonal hace que se pierdan los nombres de las variables, se los volvemos a colocar con rownames
  row.names(SE)<-colnames(X)
  print("errores estandar")
  print(SE)
  # El calculo del estadístico t, en el que se espera que los parametros sean distintos de cero
  ti<-(b-0)/SE
  print("estadístico t")
  print(ti)
  # Escalar de los grados de libertad de todo el modelo
  if(int==T) {g<-nrow(X)-ncol(X)+1
  print("grados de libertad")
  print(g)
  } else {g<-nrow(X)-ncol(X)
  print("grados de libertad")
  print(g)
  }
  # Vector de p-values
  pv<-2*pt(-abs(ti),df=g)
  print("p-values")
  print(pv)
  # Ahora tenemos los analisis de anova, que se pueden hacer 
  if(anova==T) {
    # sumas de los cuadrados
    SRR<-t(e)%*%e
    SRE<-s2
    SST<-SRR+SRE
    squares<-rbind(SRR,SRE,SST)
    row.names(squares)<-c("SRR","SRE","SST")
    print("sumas de los cuadrados")
    print(squares)
    # Prueba F
    Fp<- SRE/SRR
    print("estadístico F")
    print(Fp)
    Ft <- pf(Fp,nrow(X)-ncol(X),nrow(X)-1)
    print("valor p de la prueba F")
    print(Ft)
  }
  # finalmente el R2
    R2<-1-((t(e)%*%e)/(t(Y)%*%Y))
    print("R^2")
    print(R2)
  # como los valores no se guardan sino que únicamente se imprimen en la consola, se crea una lista con todos los resultados
    resultados<-list()
    resultados$b<-b
    resultados$s2<-s2
    resultados$SE<-SE
    resultados$ti<-ti
    resultados$g<-g
    resultados$pv<-pv
    resultados$squares<-squares
    resultados$Fp<-Fp
    resultados$Ft<-Ft
    resultados$R2<-R2
    #esta lista se puede guardar permanentemente con
    resultados<<-resultados
}

Evaluamos el resultado con la base de Bad Bunny

MiRegresion(y,x)
## [1] "betas estimados"
##                    [,1]
## Const         547.00048
## danceability  177.18502
## energy       -485.29214
## loudness       37.48488
## speechiness   619.15327
## valence      -164.31249
## [1] "Varianza del error"
##          [,1]
## [1,] 16929.53
## [1] "errores estandar"
##                   [,1]
## Const        333.42565
## danceability 206.77004
## energy       259.81529
## loudness      19.47388
## speechiness  244.64685
## valence       97.62221
## [1] "estadístico t"
##                    [,1]
## Const         1.6405471
## danceability  0.8569182
## energy       -1.8678352
## loudness      1.9248792
## speechiness   2.5308042
## valence      -1.6831465
## [1] "grados de libertad"
## [1] 38
## [1] "p-values"
##                    [,1]
## Const        0.10914443
## danceability 0.39686211
## energy       0.06951133
## loudness     0.06175469
## speechiness  0.01564158
## valence      0.10054645
## [1] "sumas de los cuadrados"
##          [,1]
## SRR 626392.54
## SRE  16929.53
## SST 643322.07
## [1] "estadístico F"
##            [,1]
## [1,] 0.02702703
## [1] "valor p de la prueba F"
##              [,1]
## [1,] 1.865678e-20
## [1] "R^2"
##           [,1]
## [1,] 0.5706075

Ejercicio 2

Realice un par de rutinas MiResumen y MiAnalDeVar que realicen, respectivamente:

  1. Imprima la tabla de los resultados de la regresión, es decir: (i) los nombres de los regresores, (ii) \(\hat{\beta}\), (iii) los \(\hat{\sigma}_\beta\), (iv) las estadísticas \(t_i\) y (v) los p-valores
MiResumen <- function(a){
  # el objeto que entra será la lista con los resultados que se guardó permanente en el environment como "resultados". Estos se pueden convertir en una tabla usando data.frame, colocando nombres a las columnas
  tabla1 <- data.frame("beta estimado" = a$b,
    "standard error" = a$SE, 
    "t statistic" = a$ti,
    "p-values" = a$pv)
  print(tabla1)
}

Ensayamos en con los datos de Bad Bunny

MiResumen(resultados)
##              beta.estimado standard.error t.statistic   p.values
## Const            547.00048      333.42565   1.6405471 0.10914443
## danceability     177.18502      206.77004   0.8569182 0.39686211
## energy          -485.29214      259.81529  -1.8678352 0.06951133
## loudness          37.48488       19.47388   1.9248792 0.06175469
## speechiness      619.15327      244.64685   2.5308042 0.01564158
## valence         -164.31249       97.62221  -1.6831465 0.10054645
  1. Imprima la tabla de análisis de varianza, es decir (i) rótulo de la fuente, (ii) las sumas de cuadrados, (iii) los grados de libertad, (iv) los cuadrados medios, (v) La estadística F y (vi) que imprima debajo de la tabla el \(R^2\)
MiAnalDeVar <- function(a){
  #En primer lugar creamos un vector con los grados de libertad. Los grados de libertad de la ecuación son k-1, los grados de libertad de los residuos son n-k, que es el que está en "g" dento de la lista "resultados", y el de los cuadrados totales es n-1
  gr<-a$g
  ge<-nrow(a$b)-1
  gt<-a$g+nrow(a$b)-1
  gl<-rbind(gr,ge,gt)
  #Calculamos un vector de cuadrados medios:
  cm<-a$squares/gl
  #Finalmente creamos otra tabla con data.frame. Este se llena columna a columna, la primera será la fuente, es decir, residuos, ecuación, y total. se conserva el mismo orden en el que se creo en el punto 1 con rbind
  tabla2 <- data.frame("Fuente" = c("Residuos", "Ecuación", "Total"),"Suma de cuadrados" = a$squares,"Grados de libertad" = gl,"Cuadrados medios" = cm, "Estadística F" = c(a$Fp, "", ""), "R^2"=c(a$R2, "", ""))
  #
  row.names(tabla2) <- NULL
  print(tabla2)
}

Probamos con el conejo malo

MiAnalDeVar(resultados)
##     Fuente Suma.de.cuadrados Grados.de.libertad Cuadrados.medios
## 1 Residuos         626392.54                 38        16484.014
## 2 Ecuación          16929.53                  5         3385.906
## 3    Total         643322.07                 43        14960.978
##       Estadística.F               R.2
## 1 0.027027027027027 0.570607513126973
## 2                                    
## 3

Funciona!!!

Ejercicio 3

En primer lugar creamos los elementos que el profesor nos da en el punto

beta <- c(1.5, -0.3, 0.6)
s_e <- 0.64
s_xx <- matrix(c(0.4, 0, 0, 0.3), ncol = 2)
m_x <- c(1,2)

Posteriormente cargamos el paquete MASS y creamos la función, aunque este se puede hacer que se lea dentro de la función y ello podría evitar errores

SimReg<-function(beta,s_e,s_xx,m_x,t){
  library(MASS)
  #simulamos la matriz de regresores
  x2<-mvrnorm(t,m_x,s_xx)
  #añadimos el intercepto
  x3<-cbind(rep(1,t),x2)
  #simulamos los errores con el supuesto de normalidad
  er<-rnorm(t,0,s_e)
  #hallamos el valor del regresando
  ys<-x3%*%beta+er
  #empleamos la rutina MiRegresion, notese que marcamos intercepto = Falso porque se lo pusimos aquí
  MiRegresion(ys,x3,int=F)
}

Ensayamos los resultados:

SimReg(beta,s_e,s_xx,m_x,1000)
## Warning: package 'MASS' was built under R version 3.6.3
## [1] "betas estimados"
##            [,1]
## [1,]  1.4777378
## [2,] -0.3625564
## [3,]  0.6247583
## [1] "Varianza del error"
##           [,1]
## [1,] 0.4155851
## [1] "errores estandar"
##            [,1]
## [1,] 0.08434297
## [2,] 0.03186910
## [3,] 0.03831584
## [1] "estadístico t"
##           [,1]
## [1,]  17.52058
## [2,] -11.37642
## [3,]  16.30548
## [1] "grados de libertad"
## [1] 997
## [1] "p-values"
##              [,1]
## [1,] 4.015009e-60
## [2,] 2.788482e-28
## [3,] 3.644731e-53
## [1] "sumas de los cuadrados"
##            [,1]
## SRR 414.3383360
## SRE   0.4155851
## SST 414.7539211
## [1] "estadístico F"
##             [,1]
## [1,] 0.001003009
## [1] "valor p de la prueba F"
##      [,1]
## [1,]    0
## [1] "R^2"
##           [,1]
## [1,] 0.9327389

Ahora hagamos la rutina de la gran simulación. Para generar multiples simulaciones debemos crear un loop, esto lo hacemos con un for. para acumular los resultados de las simulaciones recordemos que nuestros resultados de MiRegresion se guardan en la lista “resultados” por tanto creamos una lista que guarde desde “resultados” los parametros que nos pide el punto con cbind y rbind

GranSimul = function(beta, s_e, s_xx, m_x, t, n){
  # Se crean los objetos en los que se almacenarán los resultados
  betas<-NULL
  Std_Err<-NULL
  R_2<-NULL
  F_value<-NULL
  # En cada una de las n iteraciones se realiza una simulación y se almacenan los valores de interés
  for (i in 1:n) {SimReg(beta, s_e, s_xx, m_x, t)
    betas<-cbind(betas , resultados$b)
    Std_Err<-cbind(Std_Err , resultados$SE)
    R_2<-rbind(R_2 , resultados$R2)
    F_value<-rbind(F_value , resultados$Fp)
  }  
  #convertimos los objetos en matrices para poder trasponerlos
  betas<-as.matrix(betas)
  Std_Err<-as.matrix(Std_Err)
  #Los trasponemos dado que los betas estan en filas mas no en columnas, y les asignamos nombres
  betas<-t(betas)
  colnames(betas)<-c("beta1","beta2","beta3")
  Std_Err<-t(Std_Err)
  colnames(Std_Err)<-c("beta1","beta2","beta3")
  #este proceso no se hizo para R_2 ni F_value dado que estos se llenaron con rbind, porque solo eran un dato, mientas que si se llenaba con r bind los betas o los errores estándar, se obtendría solamente un listado de 1500. Ahora los convertimos en data.frames
  betas<-data.frame(betas)
  Std_Err<-data.frame(Std_Err)
  # El siguiente paso es crear los histogramas con la linea de densidad, para ello usamos ggplot2. Dentro de geom_histogram se debe especificar que el eje y sea la densidad, de lo contrario será la frecuencia y la escala no será comparable con la linea de densidad. Los histogramas se podrían plotear desde adentro de la función, pero los hacemos permanentes para poderlos plotear desde afuera, así como las pruebas de normalidad
  library(ggplot2)
  h1<<-ggplot(betas,aes(x=beta1))+geom_histogram(aes(y =..density..),col="blue",fill="cyan")+geom_density(col="orange",linetype="dashed")+ ggtitle("Histograma de Beta 1")
  h2<<-ggplot(betas,aes(x=beta2))+geom_histogram(aes(y =..density..),col="blue",fill="cyan")+geom_density(col="orange",linetype="dashed")+ ggtitle("Histograma de Beta 2")
  h3<<-ggplot(betas,aes(x=beta3))+geom_histogram(aes(y =..density..),col="blue",fill="cyan")+geom_density(col="orange",linetype="dashed")+ ggtitle("Histograma de Beta 3")
  se1<<-ggplot(Std_Err,aes(x=beta1))+geom_histogram(aes(y =..density..),col="lawngreen",fill="lemonchiffon")+geom_density(col="navy",linetype="dashed")+ ggtitle("Histograma del error estándar de Beta 1")
  se2<<-ggplot(Std_Err,aes(x=beta2))+geom_histogram(aes(y =..density..),col="lawngreen",fill="lemonchiffon")+geom_density(col="navy",linetype="dashed")+ ggtitle("Histograma del error estandar de Beta 2")
  se3<<-ggplot(Std_Err,aes(x=beta3))+geom_histogram(aes(y =..density..),col="lawngreen",fill="lemonchiffon")+geom_density(col="navy",linetype="dashed")+ ggtitle("Histograma del error estándar de Beta 3")
  #Las pruebas de normalidad las guardamos en objetos permanentes
  pb1<<-shapiro.test(betas$beta1)
  pb2<<-shapiro.test(betas$beta2)
  pb3<<-shapiro.test(betas$beta3)
}

Se realiza el ejercicio para n = 500 y tamaño = 10,30,50,100,1000,10000. En este solucionario solamente mostraremos el caso con t=1000. el codigo se oculta porque como la funcion MiRegresion imprime todos los resultados, entonces la gran simulacion imprimiría 100 mil veces los resultados y el html se haría demasiado largo. Pero el código es simplemente GranSimul(beta = beta, s_e = s_e, s_xx = s_xx, m_x = m_x, t = 1000, n = 500)

Imprimimos los gráficos

print(h1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(h2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(h3)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(se1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(se2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

print(se3)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Y finalmente imprimimos los resultados de la prueba de normalidad de Shapiro

print(pb1)
## 
##  Shapiro-Wilk normality test
## 
## data:  betas$beta1
## W = 0.99675, p-value = 0.4142
print(pb2)
## 
##  Shapiro-Wilk normality test
## 
## data:  betas$beta2
## W = 0.99812, p-value = 0.8649
print(pb3)
## 
##  Shapiro-Wilk normality test
## 
## data:  betas$beta3
## W = 0.99826, p-value = 0.9008

Los que tengan un P valor por debajo del nivel de significancia no son normales.