Parcial-2 Econometria UES GT-02
KEVIN ADONAY LEMUS PEREZ CARNET:LP16027
19/6/2020
PARTE TEÓRICA
#1. para un modelo con 2 regresores, la matriz de correlacion es de [1,0.96;0.96,1]
#Respuesta =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.
#Respuesta = La prueba jarquer Bera(en el video lo expresa)
#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:
#Respuesta = 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:
#Respuesta = 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:
#Respuesta = 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:
#Respuesta = 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:
#Respuesta = 0.51072
library(normtest)
jb.norm.test(Residuos)##
## Jarque-Bera test for normality
##
## data: Residuos
## JB = 0.51072, p-value = 0.593#8. para una tolerancia de 0.05 el VIF es de:
#Respuesta = 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: #Respuesta = 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:
#Respuesta = 0.4
VIF3<-2.5
Tolerancia<-1/VIF3
print(Tolerancia)## [1] 0.4Solución_Ejercicio_1
#Sea el conjunto de datos, indicados en el enlace de abajo, 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.
#Importar datos
library(readxl)
nivel_ventas<- read_excel("C:/Users/kevin/Desktop/segundo parcial de econometria/ventas_empresa_1.xlsx")
print(nivel_ventas)## # A tibble: 24 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
## 7 606 203 169 113
## 8 593 200 166 113
## 9 582 198 159 115
## 10 646 221 206 119
## # ... with 14 more rowscorriendo el modelo de regresion
options(scipen = 9999)
modelo_nivel_ventas<-lm(formula = V~C+P+M, data = nivel_ventas)
summary(modelo_nivel_ventas)##
## Call:
## lm(formula = V ~ C + P + M, data = nivel_ventas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.279 -6.966 1.555 6.721 14.719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 107.4435 18.0575 5.950 0.00000808 ***
## C 0.9226 0.2227 4.142 0.000505 ***
## P 0.9502 0.1558 6.097 0.00000586 ***
## M 1.2978 0.4307 3.013 0.006872 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.506 on 20 degrees of freedom
## Multiple R-squared: 0.9798, Adjusted R-squared: 0.9768
## F-statistic: 323.6 on 3 and 20 DF, p-value: < 0.00000000000000022library(fitdistrplus)
library(stargazer)
library(MASS)
library(survival)
fit_normal<-fitdist(data = modelo_nivel_ventas$residuals,distr = "norm")
plot(fit_normal)
summary(fit_normal)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 0.0000000000000004254771 1.771258
## sd 8.6773587050283431665321 1.252469
## Loglikelihood: -85.91174 AIC: 175.8235 BIC: 178.1796
## Correlation matrix:
## mean sd
## mean 1 0
## sd 0 1Prueba de normalidad jarque-Bera
library(normtest) #Carga los comandos para las pruebas de normalidad
jb.norm.test(modelo_nivel_ventas$residuals) #Ejecuta la prueba de Jarque -Bera##
## Jarque-Bera test for normality
##
## data: modelo_nivel_ventas$residuals
## JB = 1.4004, p-value = 0.272qqnorm(modelo_nivel_ventas$residuals)
qqline(modelo_nivel_ventas$residuals)
#Respuesta = No se rechaza la hipotesis nula, ya que los residuos muestran un comportamiento normal
Prueba de Normalidad de Kolmogorov - Smirnov
library(nortest) #Carga los comandos para las pruebas de normalidad
lillie.test(modelo_nivel_ventas$residuals) #Ejecuta la prueba KS, con la corrección de Lilliefors##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_nivel_ventas$residuals
## D = 0.13659, p-value = 0.2935qqnorm(modelo_nivel_ventas$residuals)
qqline(modelo_nivel_ventas$residuals)
hist(modelo_nivel_ventas$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 
#Respuesta = No se rechaza la Hipotesis nula, ya que los residuos muestran un comportamiento normal # Prueba de Normalidad de Shapiro - Wilk
shapiro.test(modelo_nivel_ventas$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_nivel_ventas$residuals
## W = 0.95315, p-value = 0.3166Normalizando W
U<- 0.0038915*((log(24))^3) - 0.083751*((log(24))^2) - 0.31082*(log(24)) -1.5861
print(U)## [1] -3.294879S<- exp(0.0030302*((log(24))^2)-0.082676*(log(24))-0.4803)
print(S)## [1] 0.4904442WN<- (log(1-0.95315)-(U))/S
print(WN)## [1] 0.4772707#Wn= 0.4772707 Respuesta
#Respuesta: #Se Concluye que no se debe rechazar la Hipotesis nula, ya que los residuos muestran un comportamiento normal
2. Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. .
#cargando los datos
library(haven)
library(readxl)
hprice1 <- read_excel("C:/Users/kevin/Desktop/segundo parcial de econometria/ventas_empresa_1.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#estimacion del modelo
options(scipen = 9999999)
library(stargazer)
Modelo_ventas_empresa<-lm(formula = V~C+P+M,data = nivel_ventas)
stargazer(Modelo_ventas_empresa,type = "html",title = "modelo ventas")| 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 |
Encontrando indice de condicion
utilizando libreria mctest para obtener el indice de condicion
library(mctest)
library(stargazer)
library(haven)
hprice1 <- read_excel("C:/Users/kevin/Desktop/segundo parcial de econometria/ventas_empresa_1.xlsx")
Modelo_ventas_empresa<-lm(formula = V~C+P+M,data = nivel_ventas)
source(file ="C:/Users/kevin/Desktop/segundo parcial de econometria/correccion_eigprop.R" )
my_eigprop(mod= Modelo_ventas_empresa)##
## Call:
## my_eigprop(mod = Modelo_ventas_empresa)
##
## 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#Respuesta = La multicolinealidad es severa ya que es mayor a 30. # Calculando el Indice de condición #Estimando nuevamente el modelo
options(scipen = 9999999)
library(stargazer)
Modelo_ventas_empresa<-lm(formula = V~C+P+M,data = nivel_ventas)
stargazer(Modelo_ventas_empresa,type = "html",title = "ventas_empresa")| 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 |
Encontrando indice de condicion
#calculo manual
#smat
xmat<-model.matrix(Modelo_ventas_empresa)
print(xmat)## (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
## 7 1 203 169 113
## 8 1 200 166 113
## 9 1 198 159 115
## 10 1 221 206 119
## 11 1 218 181 120
## 12 1 213 192 123
## 13 1 207 191 122
## 14 1 228 217 131
## 15 1 249 190 133
## 16 1 225 221 135
## 17 1 237 189 133
## 18 1 236 192 128
## 19 1 231 193 134
## 20 1 260 233 135
## 21 1 254 196 139
## 22 1 239 199 138
## 23 1 248 202 146
## 24 1 273 240 153
## attr(,"assign")
## [1] 0 1 2 3#xxmat
XXmat<-t(xmat)%*%xmat
print(XXmat)## (Intercept) C P M
## (Intercept) 24 5308 4503 2971
## C 5308 1187852 1007473 664534
## P 4503 1007473 859157 564389
## M 2971 664534 564389 372387#Sn matriz de normalizacion
Sn<-diag(1/sqrt(diag(XXmat)))
print(Sn)## [,1] [,2] [,3] [,4]
## [1,] 0.2041241 0.000000000 0.000000000 0.000000000
## [2,] 0.0000000 0.000917527 0.000000000 0.000000000
## [3,] 0.0000000 0.000000000 0.001078857 0.000000000
## [4,] 0.0000000 0.000000000 0.000000000 0.001638712# XXmat_Normalizada
XXmat_norm<-(Sn%*%XXmat)%*%Sn
print(XXmat_norm)## [,1] [,2] [,3] [,4]
## [1,] 1.0000000 0.9941322 0.9916538 0.9938018
## [2,] 0.9941322 1.0000000 0.9972774 0.9991686
## [3,] 0.9916538 0.9972774 1.0000000 0.9978035
## [4,] 0.9938018 0.9991686 0.9978035 1.0000000# calculo de los autovalores
lambas<-eigen(XXmat_norm,symmetric = TRUE)$values
print(lambas)## [1] 3.9869237681 0.0095007154 0.0027882470 0.0007872695# indice de condicion
K<-sqrt(max(lambas)/min(lambas))
print(K)## [1] 71.16349Prueba FG
library(fastGraph)## Warning: package 'fastGraph' was built under R version 3.5.3#m =variables explicativas
m<-ncol(xmat[,-1]) #cantidad de variables explicativas K-1
n<-nrow(xmat)
#determinante de la matriz de correlacion
determinante_R<-det(cor(xmat[,-1]))
Chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(Chi_FG)## [1] 71.20805# valor critico
gl<-m*(m-1)/2
VC<-qchisq(0.05,gl,lower.tail = FALSE)
print(VC)## [1] 7.814728#encontrando p-value
p_value<-pchisq(Chi_FG,gl,lower.tail = FALSE)
print(p_value)## [1] 0.000000000000002352605shadeDist(xshade=Chi_FG,ddist="dchisq",parm1=gl,lower.tail= FALSE,sub=paste("VC:",VC,"FG",Chi_FG))
#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.
FG con mctest
library(mctest)
mctest(Modelo_ventas_empresa)##
## 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 testVIF
# manualmente
VIF<-diag(solve(cor(xmat[,-1])))
print(VIF)## C P M
## 7.631451 3.838911 9.449210VIF con libreria car
#{r,message=FALSE,warning=FALSE} #library(car) #VIF_car<-vif(Modelo_ventas_empresa) #print(VIF_car) # # VIF con libreria MCTEST
library(mctest)
mc.plot(Modelo_ventas_empresa,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.
Solución_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.
options(scipen = 999999)
#Carga de objetos en formato .Rdata
load("C:/Users/kevin/Desktop/segundo parcial de econometria/wage2_1.rdata")
modelo_edu<-lm(formula = educ~sibs+meduc+feduc,data = wage2)options(scipen = 9999999)
library(stargazer)
myformula<-as.formula(educ~sibs+meduc+feduc)
modelo_edu<-lm(formula = myformula,data =wage2)
stargazer(modelo_edu,type = "html",title = "modelo_edu")| 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 |
library(fitdistrplus)
library(stargazer)
library(MASS)
library(survival)
fit_normal<-fitdist(data = modelo_edu$residuals,distr = "norm")
plot(fit_normal)
summary(fit_normal)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 0.00000000000000002159635 0.07374528
## sd 1.98153989082280057587582 0.05214573
## Loglikelihood: -1518.231 AIC: 3040.462 BIC: 3049.626
## Correlation matrix:
## mean sd
## mean 1 0
## sd 0 1Prueba de normalidad jarque-Bera
library(normtest) #Carga los comandos para las pruebas de normalidad
jb.norm.test(modelo_edu$residuals) #Ejecuta la prueba de Jarque -Bera##
## Jarque-Bera test for normality
##
## data: modelo_edu$residuals
## JB = 35.655, p-value < 0.00000000000000022#Respuesta = No se acepta la hipotesis nula, por lo tanto no hay evidencia que los residuales poseen una distribución normal
Prueba de Normalidad de Kolmogorov - Smirnov
library(nortest) #Carga los comandos para las pruebas de normalidad
lillie.test(modelo_edu$residuals) #Ejecuta la prueba KS, con la corrección de Lilliefors##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_edu$residuals
## D = 0.089992, p-value = 0.000000000000003394qqnorm(modelo_edu$residuals)
qqline(modelo_edu$residuals)
hist(modelo_edu$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 
#Respuesta = No se acepta la hipotesis nula, por lo tanto no hay evidencia que los residuales poseen una distribución normal
Prueba de Normalidad de Shapiro - Wilk
shapiro.test(modelo_edu$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_edu$residuals
## W = 0.96692, p-value = 0.00000000001058Normalizando W
U<- 0.0038915*((log(935))^3) - 0.083751*((log(935))^2) - 0.31082*(log(935)) -1.5861
print(U)## [1] -6.385615S<- exp(0.0030302*((log(935))^2)-0.082676*(log(935))-0.4803)
print(S)## [1] 0.4049237WN<- (log(1-0.95315)-(U))/S
print(WN)## [1] 8.210955#Wn=8.210955
#Respuesta = No se acepta la Hipotesis nula, por lo tanto no hay evidencia que los residuales poseen una distribución normal
2. Utilizando todas las herramientas vistas en clase, evalue la situación de colinealidad del modelo.
Encontrando indice de condicion
#calculo manual
#smat
xmat<-model.matrix(modelo_edu)
stargazer(head(xmat,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
## ------------------------------#xxmat
XXmat<-t(xmat)%*%xmat
print(XXmat)## (Intercept) sibs meduc feduc
## (Intercept) 722 2064 7802 7404
## sibs 2064 9552 20967 19949
## meduc 7802 20967 90078 83895
## feduc 7404 19949 83895 83806#Sn matriz de normalizacion
Sn<-diag(1/sqrt(diag(XXmat)))
print(Sn)## [,1] [,2] [,3] [,4]
## [1,] 0.03721615 0.00000000 0.00000000 0.000000000
## [2,] 0.00000000 0.01023182 0.00000000 0.000000000
## [3,] 0.00000000 0.00000000 0.00333189 0.000000000
## [4,] 0.00000000 0.00000000 0.00000000 0.003454319# XXmat_Normalizada
XXmat_norm<-(Sn%*%XXmat)%*%Sn
print(XXmat_norm)## [,1] [,2] [,3] [,4]
## [1,] 1.0000000 0.7859482 0.9674488 0.9518319
## [2,] 0.7859482 1.0000000 0.7147921 0.7050768
## [3,] 0.9674488 0.7147921 1.0000000 0.9655820
## [4,] 0.9518319 0.7050768 0.9655820 1.0000000# calculo de los autovalores
lambas<-eigen(XXmat_norm,symmetric = TRUE)$values
print(lambas)## [1] 3.55762739 0.37556335 0.04172605 0.02508320# indice de condicion
K<-sqrt(max(lambas)/min(lambas))
print(K)## [1] 11.90937#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.
utilizando libreria mctest para obtener el indice de condicion
library(mctest)
library(stargazer)
library(haven)
load("C:/Users/kevin/Desktop/segundo parcial de econometria/wage2_1.rdata")
modelo_estimado_2<- lm(formula = educ~sibs+meduc+feduc,data = wage2)
source(file = "C:/Users/kevin/Desktop/segundo parcial de econometria/correccion_eigprop.R")
my_eigprop(mod= modelo_edu)##
## Call:
## my_eigprop(mod = modelo_edu)
##
## 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.50Prueba de FG
library(fastGraph)
#m =variables explicativas
m<-ncol(xmat[,-1]) #cantidad de variables explicativas K-1
n<-nrow(xmat)
#determinante de la matriz de correlacion
determinante_R<-det(cor(xmat[,-1]))
Chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(Chi_FG)## [1] 358.3897# valor critico
gl<-m*(m-1)/2
VC<-qchisq(0.05,gl,lower.tail = FALSE)
print(VC)## [1] 7.814728#encontrando p-value
p_value<-pchisq(Chi_FG,gl,lower.tail = FALSE)
print(p_value)## [1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000227501shadeDist(xshade=Chi_FG,ddist="dchisq",parm1=gl,lower.tail= FALSE,sub=paste("VC:",VC,"FG",Chi_FG))
#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.
FG con mctest
library(mctest)
mctest(modelo_edu)##
## 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 testVIF
# manualmente
VIF<-diag(solve(cor(xmat[,-1])))
print(VIF)## sibs meduc feduc
## 1.098950 1.561254 1.506359VIF con libreria car
#{r,message=FALSE,warning=FALSE} #library(car) #VIF_car<-vif(modelo_edu) #print(VIF_car) #
VIF con libreria MCTEST
library(mctest)
mc.plot(modelo_edu,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
Solución_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.
#Cargando el modelo
options(scipen = 999999)
#Carga de objetos en formato .Rdata
load("C:/Users/kevin/Desktop/segundo parcial de econometria/LAWSCH85.rdata")
modelo_salario<-lm(formula = log(salary)~LSAT+GPA+log(libvol)+log(cost)+rank,data = LAWSCH85)
summary(modelo_salario)##
## Call:
## lm(formula = log(salary) ~ LSAT + GPA + log(libvol) + log(cost) +
## rank, data = LAWSCH85)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.301356 -0.084982 -0.004359 0.077935 0.288614
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.3432262 0.5325192 15.667 < 0.0000000000000002 ***
## LSAT 0.0046965 0.0040105 1.171 0.24372
## GPA 0.2475238 0.0900371 2.749 0.00683 **
## log(libvol) 0.0949932 0.0332544 2.857 0.00499 **
## log(cost) 0.0375539 0.0321061 1.170 0.24427
## rank -0.0033246 0.0003485 -9.541 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1124 on 130 degrees of freedom
## (20 observations deleted due to missingness)
## Multiple R-squared: 0.8417, Adjusted R-squared: 0.8356
## F-statistic: 138.2 on 5 and 130 DF, p-value: < 0.00000000000000022options(scipen = 9999999)
library(stargazer)
myformula<-as.formula(lsalary~LSAT+GPA+llibvol+lcost+rank)
modelo_salario1<-lm(formula = myformula,data =LAWSCH85)
stargazer(modelo_salario1,type = "html",title = "modelo salario")| 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 |
library(fitdistrplus)
library(stargazer)
library(MASS)
library(survival)
fit_normal<-fitdist(data = modelo_salario1$residuals,distr = "norm")
plot(fit_normal)
summary(fit_normal)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean -0.00000000000000000007262062 0.009424193
## sd 0.10990402695661732890286544 0.006661428
## Loglikelihood: 107.3325 AIC: -210.6649 BIC: -204.8396
## Correlation matrix:
## mean sd
## mean 1.0000000000000000000 -0.0000000000002230343
## sd -0.0000000000002230343 1.0000000000000000000Prueba de normalidad jarque-Bera
library(normtest) #Carga los comandos para las pruebas de normalidad
jb.norm.test(modelo_salario1$residuals) #Ejecuta la prueba de Jarque -Bera##
## Jarque-Bera test for normality
##
## data: modelo_salario1$residuals
## JB = 0.36511, p-value = 0.8275#Respuesta = No se rechaza la hipotesis nula, por lo tanto hay evidencia que los residuales poseen una distribución normal
Prueba de Normalidad de Kolmogorov - Smirnov
library(nortest) #Carga los comandos para las pruebas de normalidad
lillie.test(modelo_salario1$residuals) #Ejecuta la prueba KS, con la corrección de Lilliefors##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_salario1$residuals
## D = 0.054571, p-value = 0.4123qqnorm(modelo_salario1$residuals)
qqline(modelo_salario1$residuals)
hist(modelo_salario1$residuals,main = "Histograma de los residuos",xlab = "Residuos",ylab = "frecuencia") 
#Respuesta = No se rechaza la Hipotesis nula, ya que los residuos muestran un comportamiento normal
Prueba de Normalidad de Shapiro - Wilk
shapiro.test(modelo_salario1$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_salario1$residuals
## W = 0.99282, p-value = 0.7235Normalizando W
U<- 0.0038915*((log(156))^3) - 0.083751*((log(156))^2) - 0.31082*(log(156)) -1.5861
print(U)## [1] -4.7903S<- exp(0.0030302*(log(156))^2-0.082676*(log(156))-0.4803)
print(S)## [1] 0.4401989WN<- (log(1-0.95315)-(U))/S
print(WN)## [1] 3.928896#Wn= 3.928896
#Respuesta: #Se concluye que no se rechaza la Hipotesis nula, ya que los residuos muestran un comportamiento normal
#2. Utilizando todas las herramientas vistas en clase, evalue la situación de colinealidad del modelo. Comente sus resultados
Encontrando indice de condicion
#calculo manual
#smat
xmat<-model.matrix(modelo_salario)
stargazer(head(xmat,n=6),type = "text")##
## ===================================================
## (Intercept) LSAT GPA log(libvol) log(cost) 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
## ---------------------------------------------------#xxmat
XXmat<-t(xmat)%*%xmat
print(XXmat)## (Intercept) LSAT GPA log(libvol) log(cost) rank
## (Intercept) 136.0000 21557.00 450.110 783.3715 1277.008 10847.00
## LSAT 21557.0000 3419799.00 71440.370 124336.7464 202521.842 1697437.00
## GPA 450.1100 71440.37 1494.950 2599.2640 4228.169 34980.46
## log(libvol) 783.3715 124336.75 2599.264 4536.4063 7362.770 60522.09
## log(cost) 1277.0080 202521.84 4228.169 7362.7704 12010.096 100735.74
## rank 10847.0000 1697437.00 34980.460 60522.0931 100735.741 1190245.00#Sn matriz de normalizacion
Sn<-diag(1/sqrt(diag(XXmat)))
print(Sn)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.08574929 0.000000000 0.00000000 0.00000000 0.000000000 0.0000000000
## [2,] 0.00000000 0.000540754 0.00000000 0.00000000 0.000000000 0.0000000000
## [3,] 0.00000000 0.000000000 0.02586346 0.00000000 0.000000000 0.0000000000
## [4,] 0.00000000 0.000000000 0.00000000 0.01484718 0.000000000 0.0000000000
## [5,] 0.00000000 0.000000000 0.00000000 0.00000000 0.009124872 0.0000000000
## [6,] 0.00000000 0.000000000 0.00000000 0.00000000 0.000000000 0.0009166041# XXmat_Normalizada
XXmat_norm<-(Sn%*%XXmat)%*%Sn
print(XXmat_norm)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 0.9995823 0.9982420 0.9973380 0.9991966 0.8525542
## [2,] 0.9995823 1.0000000 0.9991485 0.9982590 0.9993057 0.8413471
## [3,] 0.9982420 0.9991485 1.0000000 0.9981161 0.9978511 0.8292662
## [4,] 0.9973380 0.9982590 0.9981161 1.0000000 0.9974980 0.8236445
## [5,] 0.9991966 0.9993057 0.9978511 0.9974980 1.0000000 0.8425432
## [6,] 0.8525542 0.8413471 0.8292662 0.8236445 0.8425432 1.0000000# calculo de los autovalores
lambas<-eigen(XXmat_norm,symmetric = TRUE)$values
print(lambas)## [1] 5.7351306271 0.2604004362 0.0020823558 0.0018442637 0.0003778105
## [6] 0.0001645067# indice de condicion
K<-sqrt(max(lambas)/min(lambas))
print(K)## [1] 186.7153#Respuesta = La multicolinealidad es severa ya que es mayor a 30. # utilizando libreria mctest para obtener el indice de condicion
library(mctest)
library(stargazer)
library(haven)
load("C:/Users/kevin/Desktop/segundo parcial de econometria/LAWSCH85.rdata")
modelo_salario1<- lm(formula = educ~sibs+meduc+feduc,data = wage2)
source(file = "C:/Users/kevin/Desktop/segundo parcial de econometria/correccion_eigprop.R")
my_eigprop(mod= modelo_salario1)##
## Call:
## my_eigprop(mod = modelo_salario1)
##
## 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.50prueba de FG
library(fastGraph)
#m =variables explicativas
m<-ncol(xmat[,-1]) #cantidad de variables explicativas K-1
n<-nrow(xmat)
#determinante de la matriz de correlacion
determinante_R<-det(cor(xmat[,-1]))
Chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(Chi_FG)## [1] 391.509# valor critico
gl<-m*(m-1)/2
VC<-qchisq(0.05,gl,lower.tail = FALSE)
print(VC)## [1] 18.30704#encontrando p-value
p_value<-pchisq(Chi_FG,gl,lower.tail = FALSE)
print(p_value)## [1] 0.000000000000000000000000000000000000000000000000000000000000000000000000000006031952shadeDist(xshade=Chi_FG,ddist="dchisq",parm1=gl,lower.tail= FALSE,sub=paste("VC:",VC,"FG",Chi_FG))
#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.
FG con mctest
library(mctest)
mctest(modelo_salario1)##
## 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 testVIF
# manualmente
VIF<-diag(solve(cor(xmat[,-1])))
print(VIF)## LSAT GPA log(libvol) log(cost) rank
## 3.635214 3.369004 2.110801 1.573583 3.124106VIF con libreria car
#{r,message=FALSE,warning=FALSE} #library(car) #VIF_car<-vif(modelo_salario1) #print(VIF_car) #
VIF con libreria MCTEST
library(mctest)
mc.plot(modelo_salario1,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. # el vif con libreria car se ha dejado como comentario en los tres parciales ya que no es posible instalar la libreria car en la version R #3.5.1