Luis Ernesto Juárez Linares
19 de Mayo de 2019
library(stargazer)
library(haven)
hprice1 <- read_dta("C:/Users/luisn/Downloads/hprice1.dta")
head(hprice1,n=6)## # A tibble: 6 x 10
## price assess bdrms lotsize sqrft colonial lprice lassess llotsize lsqrft
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 300 349. 4 6126 2438 1 5.70 5.86 8.72 7.80
## 2 370 352. 3 9903 2076 1 5.91 5.86 9.20 7.64
## 3 191 218. 3 5200 1374 0 5.25 5.38 8.56 7.23
## 4 195 232. 3 4600 1448 1 5.27 5.45 8.43 7.28
## 5 373 319. 4 6095 2514 1 5.92 5.77 8.72 7.83
## 6 466. 414. 5 8566 2754 1 6.14 6.03 9.06 7.92modelo_estimado<- lm(price~assess+bdrms+lotsize+colonial+llotsize,data = hprice1)
stargazer(modelo_estimado,type = "text",tixtle = "modelo estimado")##
## ===============================================
## Dependent variable:
## ---------------------------
## price
## -----------------------------------------------
## assess 0.940***
## (0.072)
##
## bdrms 8.620
## (6.791)
##
## lotsize 0.001
## (0.001)
##
## colonial 10.031
## (10.580)
##
## llotsize -13.357
## (17.813)
##
## Constant 68.090
## (146.133)
##
## -----------------------------------------------
## Observations 88
## R2 0.832
## Adjusted R2 0.822
## Residual Std. Error 43.364 (df = 82)
## F Statistic 81.224*** (df = 5; 82)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
##
## ===============
## modelo estimado
## ---------------library(stargazer)
x_mat<- model.matrix(modelo_estimado)
stargazer(head(x_mat,n=6), type = "text")##
## =====================================================
## (Intercept) assess bdrms lotsize colonial llotsize
## -----------------------------------------------------
## 1 1 349.100 4 6,126 1 8.720
## 2 1 351.500 3 9,903 1 9.201
## 3 1 217.700 3 5,200 0 8.556
## 4 1 231.800 3 4,600 1 8.434
## 5 1 319.100 4 6,095 1 8.715
## 6 1 414.500 5 8,566 1 9.056
## -----------------------------------------------------xx_mat<- t(x_mat)%*%x_mat
stargazer(xx_mat, type = "text")##
## ============================================================================================
## (Intercept) assess bdrms lotsize colonial llotsize
## --------------------------------------------------------------------------------------------
## (Intercept) 88 27,784.800 314 793,748 61 783.649
## assess 27,784.800 9,563,053.000 102,507.500 278,300,049.000 19,578.900 250,005.600
## bdrms 314 102,507.500 1,182 2,933,767 228 2,802.953
## lotsize 793,748 278,300,049.000 2,933,767 16,165,159,010 555,967 7,457,452.000
## colonial 61 19,578.900 228 555,967 61 544.060
## llotsize 783.649 250,005.600 2,802.953 7,457,452.000 544.060 7,004.230
## --------------------------------------------------------------------------------------------sn<- solve(diag(sqrt(diag(xx_mat))))
stargazer(sn, type = "text")##
## ======================================
## 0.107 0 0 0 0 0
## 0 0.0003 0 0 0 0
## 0 0 0.029 0 0 0
## 0 0 0 0.00001 0 0
## 0 0 0 0 0.128 0
## 0 0 0 0 0 0.012
## --------------------------------------\(|X^T X|\) Normalizada
options(scipen = 999999)
xx_norm<- (sn%*%xx_mat)%*%sn
stargazer(xx_norm, type = "text")##
## ===================================
## 1 0.958 0.974 0.666 0.833 0.998
## 0.958 1 0.964 0.708 0.811 0.966
## 0.974 0.964 1 0.671 0.849 0.974
## 0.666 0.708 0.671 1 0.560 0.701
## 0.833 0.811 0.849 0.560 1 0.832
## 0.998 0.966 0.974 0.701 0.832 1
## -----------------------------------Autovalores de \(|X^T X|\)
lambdas<- eigen(xx_norm, symmetric = TRUE)
stargazer(lambdas$values, type = "text")##
## ====================================
## 5.193 0.492 0.237 0.049 0.028 0.0005
## ------------------------------------calculo de \(k(x) = \sqrt(Lmax/Lmin)\)
k<- sqrt(max(lambdas$values)/min(lambdas$values))
print(k)## [1] 106.4903Se observa un problema de colinealidad, ya que el valor es k>20
library(mctest)
eigprop(x= x_mat[,-1])##
## Call:
## eigprop(x = x_mat[, -1])
##
## Eigenvalues CI Intercept assess bdrms lotsize colonial llotsize
## 1 5.1934 1.0000 0.0000 0.0014 0.0012 0.0037 0.0079 0.0000
## 2 0.4923 3.2479 0.0000 0.0001 0.0015 0.2963 0.0615 0.0000
## 3 0.2365 4.6860 0.0003 0.0113 0.0041 0.0322 0.8545 0.0002
## 4 0.0491 10.2888 0.0048 0.4576 0.0148 0.0101 0.0021 0.0027
## 5 0.0283 13.5505 0.0012 0.2138 0.9074 0.0093 0.0696 0.0017
## 6 0.0005 106.4903 0.9936 0.3158 0.0711 0.6484 0.0044 0.9954
##
## ===============================
## Row 5==> bdrms, proportion 0.907364 >= 0.50
## Row 6==> lotsize, proportion 0.648440 >= 0.50
## Row 3==> colonial, proportion 0.854491 >= 0.50
## Row 6==> llotsize, proportion 0.995427 >= 0.50sn<- solve(diag(sqrt(diag(x_mat))))
stargazer(sn, type = "text")##
## ===========================
## 1 0 0 0 0 0
## 0 0.053 0 0 0 0
## 0 0 0.577 0 0 0
## 0 0 0 0.015 0 0
## 0 0 0 0 1 0
## 0 0 0 0 0 0.332
## ---------------------------xx_norm<- (sn%*%xx_mat)%*%sn
stargazer(xx_norm, type = "text")##
## ====================================================================
## 88 1,481.988 181.288 11,703.180 61 260.414
## 1,481.988 27,206.410 3,156.693 218,862.700 1,044.301 4,431.284
## 181.288 3,156.693 394 24,973.880 131.636 537.771
## 11,703.180 218,862.700 24,973.880 3,514,165.000 8,197.286 36,538.780
## 61 1,044.301 131.636 8,197.286 61 180.796
## 260.414 4,431.284 537.771 36,538.780 180.796 773.473
## --------------------------------------------------------------------zn <- scale(x_mat[,-1])
stargazer(zn, type = "text")##
## ==========================================
## assess bdrms lotsize colonial llotsize
## ------------------------------------------
## 1 0.350 0.513 -0.284 0.662 -0.340
## 2 0.375 -0.675 0.087 0.662 0.543
## 3 -1.029 -0.675 -0.375 -1.495 -0.641
## 4 -0.881 -0.675 -0.434 0.662 -0.866
## 5 0.035 0.513 -0.287 0.662 -0.349
## 6 1.036 1.702 -0.045 0.662 0.277
## 7 0.546 -0.675 -0.002 0.662 0.367
## 8 -0.163 -0.675 -0.276 0.662 -0.315
## 9 -0.836 -0.675 -0.297 -1.495 -0.378
## 10 -0.624 -0.675 -0.602 -1.495 -1.719
## 11 -0.018 0.513 -0.297 0.662 -0.378
## 12 1.057 1.702 -0.194 0.662 -0.082
## 13 1.241 -0.675 0.316 0.662 0.932
## 14 -0.382 -0.675 -0.252 -1.495 -0.242
## 15 -0.296 -0.675 -0.246 0.662 -0.225
## 16 -0.869 0.513 -0.533 0.662 -1.318
## 17 -0.125 0.513 -0.304 -1.495 -0.402
## 18 -0.106 -0.675 -0.186 0.662 -0.063
## 19 -0.514 -0.675 -0.332 0.662 -0.491
## 20 0.108 0.513 -0.041 0.662 0.284
## 21 -0.225 -0.675 -0.347 0.662 -0.540
## 22 0.032 -0.675 -0.120 0.662 0.104
## 23 -0.226 -0.675 -0.297 -1.495 -0.377
## 24 -1.130 0.513 -0.374 -1.495 -0.635
## 25 -0.798 -0.675 0.040 0.662 0.452
## 26 -0.227 -0.675 -0.286 -1.495 -0.343
## 27 -0.507 -0.675 -0.227 -1.495 -0.172
## 28 0.463 -0.675 -0.044 0.662 0.279
## 29 1.703 4.079 -0.061 0.662 0.241
## 30 0.415 0.513 0.074 0.662 0.519
## 31 -1.028 0.513 -0.414 0.662 -0.786
## 32 0.727 0.513 0.596 -1.495 1.317
## 33 -0.959 -0.675 -0.320 0.662 -0.452
## 34 -0.670 0.513 -0.259 0.662 -0.264
## 35 0.411 0.513 -0.002 0.662 0.367
## 36 -1.083 0.513 -0.543 -1.495 -1.369
## 37 1.434 0.513 0.184 0.662 0.718
## 38 2.123 1.702 0.650 0.662 1.382
## 39 -0.276 0.513 -0.258 0.662 -0.259
## 40 -0.500 -1.864 -0.014 -1.495 0.343
## 41 -0.391 -0.675 -0.266 0.662 -0.284
## 42 3.564 1.702 1.888 0.662 2.469
## 43 -0.445 0.513 -0.194 0.662 -0.081
## 44 -1.087 -0.675 -0.365 -1.495 -0.604
## 45 0.401 1.702 -0.234 0.662 -0.192
## 46 -0.668 -0.675 -0.117 0.662 0.112
## 47 0.087 -0.675 -0.788 -1.495 -3.671
## 48 1.676 0.513 -0.089 -1.495 0.176
## 49 -0.618 -0.675 -0.312 0.662 -0.424
## 50 -0.383 0.513 -0.232 0.662 -0.186
## 51 -0.019 -0.675 -0.234 0.662 -0.192
## 52 -0.377 -1.864 0.614 -1.495 1.339
## 53 -1.228 -0.675 -0.381 0.662 -0.660
## 54 -0.989 -0.675 -0.295 0.662 -0.373
## 55 -0.497 -0.675 -0.060 0.662 0.243
## 56 -0.351 0.513 -0.334 0.662 -0.497
## 57 -0.892 0.513 -0.336 0.662 -0.505
## 58 -0.301 0.513 -0.245 0.662 -0.224
## 59 -0.179 -0.675 -0.291 0.662 -0.360
## 60 -0.012 0.513 -0.342 -1.495 -0.525
## 61 -0.260 -0.675 -0.143 -1.495 0.048
## 62 -0.308 0.513 -0.348 0.662 -0.543
## 63 -0.652 2.890 -0.361 0.662 -0.589
## 64 1.744 1.702 0.670 0.662 1.406
## 65 0.719 0.513 -0.098 0.662 0.156
## 66 2.391 0.513 0.290 0.662 0.891
## 67 0.218 0.513 -0.055 0.662 0.254
## 68 2.092 0.513 0.598 0.662 1.319
## 69 1.272 0.513 0.181 -1.495 0.713
## 70 -0.549 -0.675 -0.267 0.662 -0.288
## 71 -0.161 -0.675 0.249 -1.495 0.827
## 72 -0.682 -0.675 -0.297 0.662 -0.378
## 73 4.122 1.702 2.160 -1.495 2.641
## 74 -0.414 -0.675 -0.488 0.662 -1.098
## 75 0.764 -1.864 1.148 -1.495 1.898
## 76 -0.663 -0.675 -0.344 -1.495 -0.529
## 77 -0.215 0.513 8.223 0.662 4.654
## 78 0.459 -0.675 -0.083 0.662 0.191
## 79 -0.415 0.513 -0.302 0.662 -0.395
## 80 -0.692 -0.675 0.965 -1.495 1.725
## 81 -1.189 0.513 -0.462 0.662 -0.984
## 82 -0.648 -0.675 -0.379 0.662 -0.653
## 83 -0.094 0.513 -0.111 0.662 0.126
## 84 0.027 -0.675 -0.291 0.662 -0.361
## 85 -0.591 -0.675 -0.314 -1.495 -0.431
## 86 -0.605 -0.675 -0.263 -1.495 -0.276
## 87 -0.879 -1.864 -0.261 -1.495 -0.270
## 88 -0.669 0.513 -0.400 0.662 -0.731
## ------------------------------------------Normalizar la matriz X (se muestran las primeras 6 filas)
library(stargazer)
zn<-scale(x_mat[,-1])
stargazer(head(zn,n = 6), type = "text")##
## =========================================
## assess bdrms lotsize colonial llotsize
## -----------------------------------------
## 1 0.350 0.513 -0.284 0.662 -0.340
## 2 0.375 -0.675 0.087 0.662 0.543
## 3 -1.029 -0.675 -0.375 -1.495 -0.641
## 4 -0.881 -0.675 -0.434 0.662 -0.866
## 5 0.035 0.513 -0.287 0.662 -0.349
## 6 1.036 1.702 -0.045 0.662 0.277
## -----------------------------------------Calcular la matriz R
library(stargazer)
n<-nrow(zn)
R<-(t(zn)%*%zn)*(1/(n-1))
stargazer(R,type = "text",digits = 4)##
## ================================================
## assess bdrms lotsize colonial llotsize
## ------------------------------------------------
## assess 1 0.4825 0.3281 0.0829 0.5717
## bdrms 0.4825 1 0.1363 0.3046 0.1695
## lotsize 0.3281 0.1363 1 0.0140 0.8079
## colonial 0.0829 0.3046 0.0140 1 0.0386
## llotsize 0.5717 0.1695 0.8079 0.0386 1
## ------------------------------------------------Calcular \(|R|\)
determinante_R<-det(R)
print(determinante_R)## [1] 0.1419755Estadistico \(\chi_{FG}^2\)
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] 164.9525Valor Critico
gl<-m*(m-1)/2
VC<-qchisq(p = 0.95,df = gl)
print(VC)## [1] 18.30704Como \(\chi_{FG}^2 \geq V.C.\) se rechaza H0, por lo tanto hay evidencia de colinealidad en los regresores
library(psych)
FG_test<-cortest.bartlett(x_mat[,-1])
print(FG_test)## $chisq
## [1] 164.9525
##
## $p.value
## [1] 0.000000000000000000000000000003072151
##
## $df
## [1] 10Matriz de Correlación de los regresores del modelo (Como se obtuvo con anterioridad)
print(R)## assess bdrms lotsize colonial llotsize
## assess 1.00000000 0.4824739 0.32814633 0.08293582 0.5716654
## bdrms 0.48247394 1.0000000 0.13632563 0.30457549 0.1694902
## lotsize 0.32814633 0.1363256 1.00000000 0.01401865 0.8078552
## colonial 0.08293582 0.3045755 0.01401865 1.00000000 0.0386421
## llotsize 0.57166539 0.1694902 0.80785523 0.03864210 1.0000000Inversa de la matriz de correlación \(R^-1\)
inversa_R<-solve(R)
print(inversa_R)## assess bdrms lotsize colonial llotsize
## assess 2.1535576 -0.9010888 0.8347216 0.1520968 -1.7586001
## bdrms -0.9010888 1.5104833 -0.3640744 -0.4021983 0.5687703
## lotsize 0.8347216 -0.3640744 3.2049651 0.1130042 -3.0089889
## colonial 0.1520968 -0.4021983 0.1130042 1.1142184 -0.1531266
## llotsize -1.7586001 0.5687703 -3.0089889 -0.1531266 4.3456744VIF’s para el modelo estimado
VIFs<-diag(inversa_R)
print(VIFs)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674Obtención de los VIF’s, a través de la librería “car”
library(car)
VIFs_car<-vif(modelo_estimado)
print(VIFs_car)## assess bdrms lotsize colonial llotsize
## 2.153558 1.510483 3.204965 1.114218 4.345674Obtención de los VIF’s, a través de la librería “mctest”
library(mctest)
mc.plot(x = x_mat[,-1],y = hprice1$price,vif = 2,)