Solución del Parcial 2

Solución del Parcial 2 Econometría

PARTE I TEÓRICA

1) para un modelo con 2 regresores, la matriz de correlacion es de [1,0.96;0.96,1] el VIF es de: 12.7551

options(scipen = 999999)
regresor<-matrix(data = c(1,0.96,0.96,1),nrow = 2,ncol = 2,byrow = FALSE)
VIF<-diag(solve(regresor))
print(VIF)

## [1] 12.7551 12.7551

2) prueba robusta para evaluar la normalidad de los residuos, independientemente del tamaño de la muestra: Prueba Jarquer Bera puede utilizarse con cualquier tamaño muestral

3) En una prueba FG, para un modelo con 5 regresores, y 60 obserbaciones con un alfa de 4.3%, el estadistico de prueba dio 40,el valor del determinante de la matriz de correlacion del modelo es de: El determinante es de 0.488499

det<-40/-(59-1-(2*4+5)/6)
determinante<-exp(det)
print(determinante)

## [1] 0.488499

4) Si el VIF para una variable es de 2.5, el coeficiente de correlacion entre las variables y el resto de regresores es de: 0.6

#VIF=1/1-R^2
#2.5=1/1-R^2
1-1/2.5

## [1] 0.6

5) Sea 10, 15, -10,-15,4,-4, los residuos de un modelo, el estadistico de prueba para el contraste de normalidad de KS es: 0.1374

options(scipen = 999999)
#matriz de residuos
Residuos<-matrix(data = c(10,15,-10,-15,4,-4),nrow = 6,ncol = 1,byrow = FALSE) 
library(nortest)
lillie.test(Residuos)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  Residuos
## D = 0.1374, p-value = 0.9763

6) En una prueba de FG, para un modelo con 5 regresores y un alfa de 4.3% el valor critico es de: 18.79093

m_<-5
gl_<-m_*(m_-1)/2
# Valor crítico
VC<-qchisq(0.043,gl_,lower.tail = FALSE)
print(VC)

## [1] 18.79093

7) Sea 10, 15, -10,-15,4,-4, los residuos de un modelo, el estadistico de prueba para el contraste de normalidad de JB es: 0.51072

library(normtest)
jb.norm.test(Residuos)

## 
##  Jarque-Bera test for normality
## 
## data:  Residuos
## JB = 0.51072, p-value = 0.596

8) para una tolerancia de 0.05 el VIF es de: 20

VIF<- 1/0.05
print(VIF)

## [1] 20

9) Sea 10, 15, -10,-15,4,-4, los residuos de un modelo, el estadistico de prueba para el contraste de normalidad de Shapiro Willk es: 0.96164

Residuos_SW<-matrix(data = c(10, 15, -10, -15, 4, -4), nrow = 6, ncol= 1)
print(Residuos_SW)

##      [,1]
## [1,]   10
## [2,]   15
## [3,]  -10
## [4,]  -15
## [5,]    4
## [6,]   -4

# Prueba Shapiro Wilk
shapiro.test(Residuos_SW)

## 
##  Shapiro-Wilk normality test
## 
## data:  Residuos_SW
## W = 0.96164, p-value = 0.8323

10) Para un VIF=2.5, la tolerancia es de: 0.4

VIF3<-2.5
Tolerancia<-1/VIF3
print(Tolerancia)

## [1] 0.4

PARTE II PRÁCTICA

Ejercicio_1: Sea el conjunto de datos, indicados en el enlace, tomados en 24 meses correspondientes a los gastos de comercialización (C) de una empresa, el nivel de ventas (V), su coste de personal (P) y los costes de materias primas (M); se trata de estimar el nivel de ventas a partir de las restantes variables.

1. Estime el modelo de regresión lineal, correspondiente y verifique el supuesto de normalidad, usando todas las pruebas vistas en clase. Comente sus resultados.

PRUEBAS DE NORMALIDAD

Cargar los datos

library(readxl)
ventas_empresa <- read_excel("C:/Users/Grupo Betel/Downloads/ventas_empresa.xlsx")
modelo_ventas<-lm(formula = V~C+P+M,data = ventas_empresa)

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus)
library(stargazer)
fit_normal<-fitdist(data = modelo_ventas$residuals,distr = "norm")
plot(fit_normal)

Prueba de Normalidad de Jarque - Bera

library(normtest) 
jb.norm.test(modelo_ventas$residuals) 

## 
##  Jarque-Bera test for normality
## 
## data:  modelo_ventas$residuals
## JB = 1.4004, p-value = 0.2775

qqnorm(modelo_ventas$residuals)
qqline(modelo_ventas$residuals)

hist(modelo_ventas$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

Comentario para la prueba JB. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Kolmogorov - Smirnov

library(nortest)  
lillie.test(modelo_ventas$residuals) 

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_ventas$residuals
## D = 0.13659, p-value = 0.2935

qqnorm(modelo_ventas$residuals)
qqline(modelo_ventas$residuals)

hist(modelo_ventas$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

Comentario para la prueba KS. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Shapiro - Wilk

shapiro.test(modelo_ventas$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_ventas$residuals
## W = 0.95315, p-value = 0.3166

Comentario para la prueba WS. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Calcular Wn

U<- 0.0038915*(log(24))^3 - 0.083751*(log(24))^2 - 0.31082*(log(24)) -1.5861
print(U)

## [1] -3.294879

S<- exp(0.0030302*(log(24))^2-0.082676*(log(24))-0.4803)
print(S)

## [1] 0.4904442

WN<- (log(1-0.95315)-(U))/S
print(WN)

## [1] 0.4772707

Prueba de MULTICOLINEALIDAD

2. Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. Comente sus resultados.

Cargando los datos

library(haven)
hprice1 <- read_excel("C:/Users/Grupo Betel/Downloads/ventas_empresa.xlsx")
head(hprice1,n=6)

## # A tibble: 6 x 4
##       V     C     P     M
##   <dbl> <dbl> <dbl> <dbl>
## 1   607   197   173   110
## 2   590   208   152   107
## 3   543   181   150    99
## 4   558   194   150   102
## 5   571   192   163   109
## 6   615   196   179   114

Estimando el modelo

library(stargazer)
modelo_estimado<-lm(formula = V~C+P+M,data = ventas_empresa)
stargazer(modelo_estimado,type = "text",title = "modelo estimado")

## 
## modelo estimado
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  V             
## -----------------------------------------------
## C                            0.923***          
##                               (0.223)          
##                                                
## P                            0.950***          
##                               (0.156)          
##                                                
## M                            1.298***          
##                               (0.431)          
##                                                
## Constant                    107.444***         
##                              (18.057)          
##                                                
## -----------------------------------------------
## Observations                    24             
## R2                             0.980           
## Adjusted R2                    0.977           
## Residual Std. Error       9.506 (df = 20)      
## F Statistic           323.641*** (df = 3; 20)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Calcular la matriz |Xt X| de forma MANUAL

library(stargazer)
X_mat<-model.matrix(modelo_estimado)
stargazer(head(X_mat,n=6),type="text")

## 
## =========================
##   (Intercept)  C   P   M 
## -------------------------
## 1      1      197 173 110
## 2      1      208 152 107
## 3      1      181 150 99 
## 4      1      194 150 102
## 5      1      192 163 109
## 6      1      196 179 114
## -------------------------

XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "text")

## 
## ===================================================
##             (Intercept)     C         P        M   
## ---------------------------------------------------
## (Intercept)     24        5,308     4,503    2,971 
## C              5,308    1,187,852 1,007,473 664,534
## P              4,503    1,007,473  859,157  564,389
## M              2,971     664,534   564,389  372,387
## ---------------------------------------------------

Cálculo de la matriz de normalización

library(stargazer)
options(scipen = 999999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "text")

## 
## =======================
## 0.204   0     0     0  
## 0     0.001   0     0  
## 0       0   0.001   0  
## 0       0     0   0.002
## -----------------------

# Tambien se puede hacer a traves de 
Sn_1<- diag(1/sqrt(diag(XX_matrix)))

# Matriz Normalizada
library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "text",digits = 4)

## 
## ===========================
## 1      0.9941 0.9917 0.9938
## 0.9941   1    0.9973 0.9992
## 0.9917 0.9973   1    0.9978
## 0.9938 0.9992 0.9978   1   
## ---------------------------

Calculando el indice de condicion

#Autovalores de |XtX| Normalizada:
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")

## 
## =======================
## 3.987 0.010 0.003 0.001
## -----------------------

#Cálculo de κ(x)
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)

## [1] 71.16349

Comentario: la multicolinealidad es severa ya que es mayor a 30.

Uso de la libreria “mctest”, para obtener el indice de condición

library(mctest)
library(stargazer)
library(haven)
hprice1 <- read_excel("C:/Users/Grupo Betel/Downloads/ventas_empresa.xlsx")
modelo_estimado<-lm(formula = V~C+P+M,data = ventas_empresa)
source(file = "C:/Users/Grupo Betel/Downloads/correccion_eigprop.R")
my_eigprop(mod= modelo_estimado)

## 
## Call:
## my_eigprop(mod = modelo_estimado)
## 
##   Eigenvalues      CI (Intercept)      C      P      M
## 1      3.9869  1.0000      0.0007 0.0001 0.0003 0.0001
## 2      0.0095 20.4852      0.8776 0.0049 0.0877 0.0075
## 3      0.0028 37.8141      0.1183 0.1594 0.8478 0.0636
## 4      0.0008 71.1635      0.0034 0.8356 0.0642 0.9288
## 
## ===============================
## Row 4==> C, proportion 0.835554 >= 0.50 
## Row 3==> P, proportion 0.847805 >= 0.50 
## Row 4==> M, proportion 0.928751 >= 0.50

El índice de condición es igual a 71.1635 y cae dentro del rango mayor a 30, por lo que se tiene evidencia de multicolinealidad severa en los regresores del modelo de ventas.

CALCULANDO LA PRUEBA DE FG

# Cargando datos
library(haven)
ventas_empresa<- read_excel("C:/Users/Grupo Betel/Downloads/ventas_empresa.xlsx")

# Estimando el modelo
modelo_estimado<-lm(formula = V~C+P+M,data = ventas_empresa)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_estimado)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "text")

## 
## ======================
##     C      P      M   
## ----------------------
## 1 -0.983 -0.587 -0.975
## 2 -0.536 -1.430 -1.187
## 3 -1.634 -1.510 -1.753
## 4 -1.105 -1.510 -1.541
## 5 -1.186 -0.988 -1.046
## 6 -1.024 -0.346 -0.692
## ----------------------

Calculando la matriz R

library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
#También se puede calcular R a través de cor(X_mat[,-1])
R_1<- cor(X_mat[,-1])
stargazer(R,type = "text",digits = 4)

## 
## ======================
##     C      P      M   
## ----------------------
## C   1    0.8205 0.9312
## P 0.8205   1    0.8579
## M 0.9312 0.8579   1   
## ----------------------

Calculando el determinante de R

determinante_R<-det(R)
print(determinante_R)

## [1] 0.03459107

Aplicando la prueba de Farrer Glaubar (Bartlett)

# Estimando el Estadistico
m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)

## [1] 71.20805

Calculando el Valor Crítico

gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)

## [1] 7.814728

# Graficar
library(fastGraph)
shadeDist(xshade = chi_FG, ddist = "dchisq", parm1 = gl, lower.tail = FALSE, sub= paste("VC:", VC))

FG USANDO MCTEST

library(mctest)
mctest(modelo_estimado)

## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0346         0
## Farrar Chi-Square:        71.2080         1
## Red Indicator:             0.8711         1
## Sum of Lambda Inverse:    20.9196         1
## Theil's Method:            0.5430         1
## Condition Number:        105.2299         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

FG USANDO de la libreria PSYCH

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test)

## $chisq
## [1] 71.20805
## 
## $p.value
## [1] 0.000000000000002352605
## 
## $df
## [1] 3

# Graficar
library(fastGraph)
shadeDist(xshade = FG_test$chisq, ddist = "dchisq", parm1 = FG_test$df, lower.tail = FALSE, sub= paste("VC:", "FG", FG_test$chist))

Con un estadístico FG de 71.20805 mayor al valor critico que es 7.81472, se rechaza la hipótesis nula ya que la prueba evidencia que existe alta multicolinealidad en los regresores del modelo.

CALCULANDO LOS VIF DE FORMA MANUAL

# Matriz de Correlación de los regresores del modelo
# Cargando datos
library(haven)
hprice1 <- read_excel("C:/Users/Grupo Betel/Downloads/ventas_empresa.xlsx")

# Estimando el modelo
modelo_estimado<-lm(formula = V~C+P+M,data = ventas_empresa)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_estimado)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "text",digits = 6)

## 
## ============================
##      C        P        M    
## ----------------------------
## C    1     0.820452 0.931240
## P 0.820452    1     0.857916
## M 0.931240 0.857916    1    
## ----------------------------

Inversa de la matriz de correlación R

inversa_R<-solve(R)
print(inversa_R)

##            C          P         M
## C  7.6314513 -0.6223111 -6.572822
## P -0.6223111  3.8389114 -2.713943
## M -6.5728221 -2.7139427  9.449210

VIF’s para el modelo estimado

VIFs<-diag(inversa_R)
print(VIFs)

##        C        P        M 
## 7.631451 3.838911 9.449210

VIF CON LIBRERIA CAR

library(car)
VIFs_car<-vif(modelo_estimado)
print(VIFs_car)

##        C        P        M 
## 7.631451 3.838911 9.449210

VIF’s, a través de la librería “mctest”

library(mctest)
mc.plot(modelo_estimado, vif = 2)

El VIF muestra que las variables independientes superan el umbral que se ha establecido, siendo los gastos de comercialización (C) y los costes de materia prima (M) los regresores que más inflan la varianza y aportan un nivel alto de multicolinealidad en el modelo de ventas.

Ejercicio_2: Se tienen los datos para trabajadores hombres,en el archivo adjunto, con ellos estime un modelo donde educ es años de escolaridad, como variable dependiente, y como regresores sibs (número de hermanos), meduc (años de escolaridad de la madre) y feduc (años de escolaridad del padre)

1. Estime el modelo de regresión lineal, correspondiente y verifique el supuesto de normalidad, usando todas las pruebas vistas en clase. Comente sus resultados.

PRUEBAS DE NORMALIDAD

Cargar los datos

library(readxl)
load("C:/Users/Grupo Betel/Downloads/wage2.RData")
modelo_educ<- lm(formula = educ~sibs+meduc+feduc,data = wage2) 

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus)
library(stargazer)
fit_normal<-fitdist(data = modelo_educ$residuals,distr = "norm")
plot(fit_normal)

Prueba de Normalidad de Jarque - Bera

library(normtest) 
jb.norm.test(modelo_educ$residuals) 

## 
##  Jarque-Bera test for normality
## 
## data:  modelo_educ$residuals
## JB = 35.655, p-value < 0.00000000000000022

qqnorm(modelo_educ$residuals)
qqline(modelo_educ$residuals)

hist(modelo_educ$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

Comentario para la prueba JB. se rechaza H0, por lo tanto no hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Kolmogorov - Smirnov

library(nortest)  
lillie.test(modelo_educ$residuals) 

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_educ$residuals
## D = 0.089992, p-value = 0.000000000000003394

qqnorm(modelo_educ$residuals)
qqline(modelo_educ$residuals)

hist(modelo_educ$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

Comentario para la prueba KS. se rechaza H0, por lo tanto no hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Shapiro - Wilk

shapiro.test(modelo_educ$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_educ$residuals
## W = 0.96692, p-value = 0.00000000001058

Comentario para la prueba WS. se rechaza H0, por lo tanto no hay evidencia que los residuales poseen una distribución normal

Calcular Wn

U<- 0.0038915*(log(935))^3 - 0.083751*(log(935))^2 - 0.31082*(log(935)) -1.5861
print(U)

## [1] -6.385615

S<- exp(0.0030302*(log(935))^2-0.082676*(log(935))-0.4803)
print(S)

## [1] 0.4049237

WN<- (log(1-0.95315)-(U))/S
print(WN)

## [1] 8.210955

Prueba de MULTICOLINEALIDAD

2. Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. Comente sus resultados.

Cargando los datos

library(haven)
load("C:/Users/Grupo Betel/Downloads/wage2.RData")

Estimando el modelo

library(stargazer)
modelo_educ<- lm(formula = educ~sibs+meduc+feduc,data = wage2)
stargazer(modelo_educ,type = "text",title = "Modelo de educación")

## 
## Modelo de educación
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                educ            
## -----------------------------------------------
## sibs                         -0.094***         
##                               (0.034)          
##                                                
## meduc                        0.131***          
##                               (0.033)          
##                                                
## feduc                        0.210***          
##                               (0.027)          
##                                                
## Constant                     10.364***         
##                               (0.359)          
##                                                
## -----------------------------------------------
## Observations                    722            
## R2                             0.214           
## Adjusted R2                    0.211           
## Residual Std. Error      1.987 (df = 718)      
## F Statistic           65.198*** (df = 3; 718)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Calcular la matriz |Xt X| de forma MANUAL

library(stargazer)
X_mat<-model.matrix(modelo_educ)
stargazer(head(X_mat,n=6),type="text")

## 
## ==============================
##   (Intercept) sibs meduc feduc
## ------------------------------
## 1      1       1     8     8  
## 2      1       1    14    14  
## 3      1       1    14    14  
## 4      1       4    12    12  
## 5      1       10    6    11  
## 7      1       1     8     8  
## ------------------------------

XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "text")

## 
## ============================================
##             (Intercept)  sibs  meduc  feduc 
## --------------------------------------------
## (Intercept)     722     2,064  7,802  7,404 
## sibs           2,064    9,552  20,967 19,949
## meduc          7,802    20,967 90,078 83,895
## feduc          7,404    19,949 83,895 83,806
## --------------------------------------------

Cálculo de la matriz de normalización

library(stargazer)
options(scipen = 999999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "text")

## 
## =======================
## 0.037   0     0     0  
## 0     0.010   0     0  
## 0       0   0.003   0  
## 0       0     0   0.003
## -----------------------

# Tambien se puede hacer a traves de 
Sn_1<- diag(1/sqrt(diag(XX_matrix)))

# Matriz Normalizada
library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "text",digits = 4)

## 
## ===========================
## 1      0.7859 0.9674 0.9518
## 0.7859   1    0.7148 0.7051
## 0.9674 0.7148   1    0.9656
## 0.9518 0.7051 0.9656   1   
## ---------------------------

Calculando el indice de condicion

#Autovalores de |XtX| Normalizada:
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")

## 
## =======================
## 3.558 0.376 0.042 0.025
## -----------------------

#Cálculo de κ(x)
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)

## [1] 11.90937

Comentario: la multicolinealidad es severa ya que es mayor a 30.

Uso de la libreria “mctest”, para obtener el indice de condición

library(mctest)
library(stargazer)
library(haven)
load("C:/Users/Grupo Betel/Downloads/wage2.RData")
modelo_educ<- lm(formula = educ~sibs+meduc+feduc,data = wage2)
source(file = "C:/Users/Grupo Betel/Downloads/correccion_eigprop.R")
my_eigprop(mod= modelo_educ)

## 
## Call:
## my_eigprop(mod = modelo_educ)
## 
##   Eigenvalues      CI (Intercept)   sibs  meduc  feduc
## 1      3.5576  1.0000      0.0033 0.0194 0.0031 0.0046
## 2      0.3756  3.0778      0.0015 0.7200 0.0107 0.0184
## 3      0.0417  9.2337      0.3235 0.1056 0.0813 0.8786
## 4      0.0251 11.9094      0.6717 0.1549 0.9049 0.0984
## 
## ===============================
## Row 2==> sibs, proportion 0.720032 >= 0.50 
## Row 4==> meduc, proportion 0.904919 >= 0.50 
## Row 3==> feduc, proportion 0.878599 >= 0.50

El índice de condición es igual a 11.90937 y cae dentro del rango menor a 20, por lo que se tiene evidencia de que la multicolinealidad en este caso es leve y no se considera un problema.

CALCULANDO LA PRUEBA DE FG

# Cargando datos
library(haven)
load("C:/Users/Grupo Betel/Downloads/wage2.RData")

# Estimando el modelo
modelo_educ<- lm(formula = educ~sibs+meduc+feduc,data = wage2)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_educ)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "text")

## 
## ======================
##    sibs  meduc  feduc 
## ----------------------
## 1 -0.826 -0.992 -0.682
## 2 -0.826 1.129  1.133 
## 3 -0.826 1.129  1.133 
## 4 0.507  0.422  0.528 
## 5 3.173  -1.699 0.225 
## 7 -0.826 -0.992 -0.682
## ----------------------

Calculando la matriz R

library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
#También se puede calcular R a través de cor(X_mat[,-1])
R_1<- cor(X_mat[,-1])
stargazer(R,type = "text",digits = 4)

## 
## =============================
##        sibs    meduc   feduc 
## -----------------------------
## sibs     1    -0.2913 -0.2269
## meduc -0.2913    1    0.5765 
## feduc -0.2269 0.5765     1   
## -----------------------------

Calculando el determinante de R

determinante_R<-det(R)
print(determinante_R)

## [1] 0.6075382

Aplicando la prueba de Farrer Glaubar (Bartlett)

# Estimando el Estadistico
m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)

## [1] 358.3897

Calculando el Valor Critico

gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)

## [1] 7.814728

# Graficar
library(fastGraph)
shadeDist(xshade = chi_FG, ddist = "dchisq", parm1 = gl, lower.tail = FALSE, sub= paste("VC:", VC))

FG USANDO MCTEST

library(mctest)
mctest(modelo_educ)

## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.6075         0
## Farrar Chi-Square:       358.3897         1
## Red Indicator:             0.3952         0
## Sum of Lambda Inverse:     4.1666         0
## Theil's Method:            0.3575         0
## Condition Number:         11.2768         0
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

FG USANDO de la libreria PSYCH

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test)

## $chisq
## [1] 358.3897
## 
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000227501
## 
## $df
## [1] 3

# Graficar
library(fastGraph)
shadeDist(xshade = FG_test$chisq, ddist = "dchisq", parm1 = FG_test$df, lower.tail = FALSE, sub= paste("VC:", "FG", FG_test$chist))

Con un estadistico FG de 358.3897 mayor al valor crítico que es 7.814728 se rechaza la hipótesis nula ya que la prueba evidencia que existe multicolinealidad severa en los regresores del modelo de ventas

CALCULANDO LOS VIF DE FORMA MANUAL

# Matriz de Correlación de los regresores del modelo
# Cargando datos
library(haven)
load("C:/Users/Grupo Betel/Downloads/wage2.RData")

# Estimando el modelo
modelo_educ<- lm(formula = educ~sibs+meduc+feduc,data = wage2)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_educ)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "text",digits = 6)

## 
## ===================================
##         sibs      meduc     feduc  
## -----------------------------------
## sibs      1     -0.291255 -0.226889
## meduc -0.291255     1     0.576495 
## feduc -0.226889 0.576495      1    
## -----------------------------------

Inversa de la matriz de correlación R

inversa_R<-solve(R)
print(inversa_R)

##             sibs      meduc       feduc
## sibs  1.09894956  0.2641069  0.09708311
## meduc 0.26410686  1.5612542 -0.84013197
## feduc 0.09708311 -0.8401320  1.50635875

VIF’s para el modelo estimado

VIFs<-diag(inversa_R)
print(VIFs)

##     sibs    meduc    feduc 
## 1.098950 1.561254 1.506359

VIF CON LIBRERIA CAR

library(car)
VIFs_car<-vif(modelo_educ)
print(VIFs_car)

##     sibs    meduc    feduc 
## 1.098950 1.561254 1.506359

VIF’s, a través de la librería “mctest”

library(mctest)
mc.plot(modelo_educ, vif = 2)

El VIF plot muestra que ninguno de los regresores supera el umbral de 2 y por lo tanto no existe evidencia de multicolinealidad severa, sino que en este caso, se trata de una colinealidad leve

Ejercicio_3: El sueldo inicial medio (salary) para los recién graduados de la Facultad de Economía se determina mediante una función lineal: log(salary)=f(SAT,GPA ,log(libvol),log(cost),rank) Donde LSAT es la media del puntaje LSAT del grupo de graduados, GPA es la media del GPA (promedio general) del grupo, libvol es el número de volúmenes en la biblioteca de la Facultad de Economía, cost es el costo anual por asistir a dicha facultad y rank es una clasificación de las escuelas de Economía (siendo rank 1 la mejor)

1. Estime el modelo de regresión lineal, correspondiente y verifique el supuesto de normalidad, usando todas las pruebas vistas en clase. Comente sus resultados.

PRUEBAS DE NORMALIDAD

Cargar los datos

library(readxl)
load("C:/Users/Grupo Betel/Downloads/LAWSCH85.RData")
modelo_sueldo<- lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data = LAWSCH85) 

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus)
library(stargazer)
fit_normal<-fitdist(data = modelo_sueldo$residuals,distr = "norm")
plot(fit_normal)

Prueba de Normalidad de Jarque - Bera

library(normtest) 
jb.norm.test(modelo_sueldo$residuals) 

## 
##  Jarque-Bera test for normality
## 
## data:  modelo_sueldo$residuals
## JB = 0.36511, p-value = 0.827

qqnorm(modelo_sueldo$residuals)
qqline(modelo_sueldo$residuals)

hist(modelo_sueldo$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

### Comentario para la prueba JB. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Kolmogorov - Smirnov

library(nortest)  
lillie.test(modelo_sueldo$residuals) 

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_sueldo$residuals
## D = 0.054571, p-value = 0.4123

qqnorm(modelo_sueldo$residuals)
qqline(modelo_sueldo$residuals)

hist(modelo_sueldo$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 

Comentario para la prueba KS. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Prueba de Normalidad de Shapiro - Wilk

shapiro.test(modelo_sueldo$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_sueldo$residuals
## W = 0.99282, p-value = 0.7235

Comentario para la prueba WS. no se rechaza H0, por lo tanto hay evidencia que los residuales poseen una distribución normal

Calcular Wn

U<- 0.0038915*(log(156))^3 - 0.083751*(log(156))^2 - 0.31082*(log(156)) -1.5861
print(U)

## [1] -4.7903

S<- exp(0.0030302*(log(156))^2-0.082676*(log(156))-0.4803)
print(S)

## [1] 0.4401989

WN<- (log(1-0.95315)-(U))/S
print(WN)

## [1] 3.928896

Prueba de MULTICOLINEALIDAD

2. Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. Comente sus resultados.

Cargando los datos

library(haven)
load("C:/Users/Grupo Betel/Downloads/LAWSCH85.RData")

Estimando el modelo

library(stargazer)
modelo_sueldo<- lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data = LAWSCH85)
stargazer(modelo_sueldo,type = "text",title = "Modelo de sueldo")

## 
## Modelo de sueldo
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               lsalary          
## -----------------------------------------------
## LSAT                           0.005           
##                               (0.004)          
##                                                
## GPA                          0.248***          
##                               (0.090)          
##                                                
## llibvol                      0.095***          
##                               (0.033)          
##                                                
## lcost                          0.038           
##                               (0.032)          
##                                                
## rank                         -0.003***         
##                              (0.0003)          
##                                                
## Constant                     8.343***          
##                               (0.533)          
##                                                
## -----------------------------------------------
## Observations                    136            
## R2                             0.842           
## Adjusted R2                    0.836           
## Residual Std. Error      0.112 (df = 130)      
## F Statistic          138.230*** (df = 5; 130)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Calcular la matriz |Xt X| de forma MANUAL

library(stargazer)
X_mat<-model.matrix(modelo_sueldo)
stargazer(head(X_mat,n=6),type="text")

## 
## ===========================================
##   (Intercept) LSAT  GPA  llibvol lcost rank
## -------------------------------------------
## 1      1      155  3.150  5.375  9.029 128 
## 2      1      160  3.500  5.545  8.851 104 
## 3      1      155  3.250  6.050  9.703  34 
## 4      1      157  3.200  5.796  9.774  49 
## 5      1      162  3.380  5.805  9.030  95 
## 6      1      161  3.400  5.740  9.030  98 
## -------------------------------------------

XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "text")

## 
## ==================================================================================
##             (Intercept)    LSAT        GPA       llibvol      lcost       rank    
## ----------------------------------------------------------------------------------
## (Intercept)     136       21,557     450.110     783.372    1,277.008    10,847   
## LSAT          21,557     3,419,799  71,440.370 124,336.700 202,521.800  1,697,437 
## GPA           450.110   71,440.370  1,494.950   2,599.264   4,228.169  34,980.460 
## llibvol       783.372   124,336.700 2,599.264   4,536.406   7,362.770  60,522.090 
## lcost        1,277.008  202,521.800 4,228.169   7,362.770  12,010.100  100,735.700
## rank          10,847     1,697,437  34,980.460 60,522.090  100,735.700  1,190,245 
## ----------------------------------------------------------------------------------

Cálculo de la matriz de normalización

library(stargazer)
options(scipen = 999999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "text")

## 
## ===================================
## 0.086   0     0     0     0     0  
## 0     0.001   0     0     0     0  
## 0       0   0.026   0     0     0  
## 0       0     0   0.015   0     0  
## 0       0     0     0   0.009   0  
## 0       0     0     0     0   0.001
## -----------------------------------

# Tambien se puede hacer a traves de 
Sn_1<- diag(1/sqrt(diag(XX_matrix)))

# Matriz Normalizada
library(stargazer)
XX_norm<-(Sn%*%XX_matrix)%*%Sn
stargazer(XX_norm,type = "text",digits = 4)

## 
## =========================================
## 1      0.9996 0.9982 0.9973 0.9992 0.8526
## 0.9996   1    0.9991 0.9983 0.9993 0.8413
## 0.9982 0.9991   1    0.9981 0.9979 0.8293
## 0.9973 0.9983 0.9981   1    0.9975 0.8236
## 0.9992 0.9993 0.9979 0.9975   1    0.8425
## 0.8526 0.8413 0.8293 0.8236 0.8425   1   
## -----------------------------------------

Calculando el indice de condicion

#Autovalores de |XtX| Normalizada:
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")

## 
## =====================================
## 5.735 0.260 0.002 0.002 0.0004 0.0002
## -------------------------------------

#Cálculo de κ(x)
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)

## [1] 186.7153

Comentario: la multicolinealidad es severa ya que es mayor a 30.

Uso de la libreria “mctest”, para obtener el indice de condición

library(mctest)
library(stargazer)
library(haven)
load("C:/Users/Grupo Betel/Downloads/LAWSCH85.RData")
modelo_sueldo<- lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data = LAWSCH85)
source(file = "C:/Users/Grupo Betel/Downloads/correccion_eigprop.R")
my_eigprop(mod= modelo_sueldo)

## 
## Call:
## my_eigprop(mod = modelo_sueldo)
## 
##   Eigenvalues       CI (Intercept)   LSAT    GPA llibvol  lcost   rank
## 1      5.7351   1.0000      0.0000 0.0000 0.0000  0.0001 0.0000 0.0021
## 2      0.2604   4.6930      0.0000 0.0000 0.0002  0.0004 0.0001 0.2884
## 3      0.0021  52.4800      0.0058 0.0030 0.0007  0.8411 0.1155 0.1357
## 4      0.0018  55.7648      0.0002 0.0010 0.3355  0.1095 0.1756 0.0161
## 5      0.0004 123.2068      0.4254 0.0588 0.4407  0.0423 0.6610 0.4700
## 6      0.0002 186.7153      0.5686 0.9371 0.2229  0.0066 0.0478 0.0877
## 
## ===============================
## Row 6==> LSAT, proportion 0.937119 >= 0.50 
## Row 3==> llibvol, proportion 0.841136 >= 0.50 
## Row 5==> lcost, proportion 0.661004 >= 0.50

El indice de condicion es de 186.7153 un valor muy por encima de 30}, lo cual representa evidencia de multicolinealidad severa en el modelo.

CALCULANDO LA PRUEBA DE FG

# Cargando datos
library(haven)
load("C:/Users/Grupo Betel/Downloads/LAWSCH85.RData")

# Estimando el modelo
modelo_sueldo<- lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data = LAWSCH85)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_sueldo)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
stargazer(head(Zn,n=6),type = "text")

## 
## =====================================
##    LSAT   GPA   llibvol lcost   rank 
## -------------------------------------
## 1 -0.763 -0.809 -0.910  -0.955 0.983 
## 2 0.325  0.965  -0.508  -1.426 0.494 
## 3 -0.763 -0.302  0.685  0.829  -0.932
## 4 -0.328 -0.556  0.085  1.016  -0.627
## 5 0.759  0.357   0.107  -0.952 0.311 
## 6 0.542  0.458  -0.048  -0.952 0.372 
## -------------------------------------

Calculando la matriz R

library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
#También se puede calcular R a través de cor(X_mat[,-1])
R_1<- cor(X_mat[,-1])
stargazer(R,type = "text",digits = 4)

## 
## ===============================================
##          LSAT     GPA   llibvol  lcost   rank  
## -----------------------------------------------
## LSAT       1    0.7727  0.6348  0.4545  -0.7184
## GPA     0.7727     1    0.5857  0.1730  -0.7034
## llibvol 0.6348  0.5857     1    0.3290  -0.6990
## lcost   0.4545  0.1730  0.3290     1    -0.4452
## rank    -0.7184 -0.7034 -0.6990 -0.4452    1   
## -----------------------------------------------

Calculando el determinante de R

determinante_R<-det(R)
print(determinante_R)

## [1] 0.05208986

Aplicando la prueba de Farrer Glaubar (Bartlett)

# Estimando el Estadístico
m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)

## [1] 391.509

Calculando el Valor Crítico

gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)

## [1] 18.30704

# Graficar
library(fastGraph)
shadeDist(xshade = chi_FG, ddist = "dchisq", parm1 = gl, lower.tail = FALSE, sub= paste("VC:", VC))

FG USANDO MCTEST

library(mctest)
mctest(modelo_sueldo)

## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0521         0
## Farrar Chi-Square:       391.5090         1
## Red Indicator:             0.5819         1
## Sum of Lambda Inverse:    13.8127         0
## Theil's Method:           -0.3680         0
## Condition Number:        181.9505         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

FG USANDO de la libreria PSYCH

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test)

## $chisq
## [1] 391.509
## 
## $p.value
## [1] 0.000000000000000000000000000000000000000000000000000000000000000000000000000006031929
## 
## $df
## [1] 10

# Graficar
library(fastGraph)
shadeDist(xshade = FG_test$chisq, ddist = "dchisq", parm1 = FG_test$df, lower.tail = FALSE, sub= paste("VC:", "FG", FG_test$chist))

Con un estadístico FG de 391.509 mayor al valor critico que es 18.3079 se rechaza la hipótesis nula y se evidencia que existe alta multicolinealidad en los regresores del modelo.

CALCULANDO LOS VIF DE FORMA MANUAL

# Matriz de Correlación de los regresores del modelo
# Cargando datos
library(haven)
load("C:/Users/Grupo Betel/Downloads/LAWSCH85.RData")

# Estimando el modelo
modelo_sueldo<- lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank,data = LAWSCH85)

# Calcular la matriz x del modelo
X_mat<-model.matrix(modelo_sueldo)

# Normalizar la matriz X
library(stargazer)
Zn<-scale(X_mat[,-1])
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R,type = "text",digits = 6)

## 
## =========================================================
##           LSAT       GPA     llibvol    lcost     rank   
## ---------------------------------------------------------
## LSAT        1     0.772659  0.634751  0.454493  -0.718443
## GPA     0.772659      1     0.585743  0.173031  -0.703416
## llibvol 0.634751  0.585743      1     0.328961  -0.699046
## lcost   0.454493  0.173031  0.328961      1     -0.445248
## rank    -0.718443 -0.703416 -0.699046 -0.445248     1    
## ---------------------------------------------------------

Inversa de la matriz de correlación R

inversa_R<-solve(R)
print(inversa_R)

##               LSAT        GPA     llibvol       lcost      rank
## LSAT     3.6352138 -2.0934988 -0.52380567 -0.96473975 0.3433791
## GPA     -2.0934988  3.3690041 -0.10955061  0.94287625 1.2089864
## llibvol -0.5238057 -0.1095506  2.11080158  0.02215876 1.0320295
## lcost   -0.9647398  0.9428763  0.02215876  1.57358251 0.6862486
## rank     0.3433791  1.2089864  1.03202954  0.68624864 3.1241059

VIF’s para el modelo estimado

VIFs<-diag(inversa_R)
print(VIFs)

##     LSAT      GPA  llibvol    lcost     rank 
## 3.635214 3.369004 2.110802 1.573583 3.124106

VIF CON LIBRERIA CAR

library(car)
VIFs_car<-vif(modelo_sueldo)
print(VIFs_car)

##     LSAT      GPA  llibvol    lcost     rank 
## 3.635214 3.369004 2.110802 1.573583 3.124106

VIF’s, a través de la librería “mctest”

library(mctest)
mc.plot(modelo_sueldo, vif = 2)

El VIF plot muestra que las variables independientes que superan el umbral de VIF que se ha establecido, produce que la multicolinealidad del modelo de salario sea severa.