Como leer los datos

#Si son .csv USAMOS LA LIBRERIA(readr) EJ: read.csv("")
#Si son .dta USAMOS LA LIBRERIA(haven) EJ: read_dta("")
#Si son .RData USAMOS load("")

Todas las librerías a utilizar

library(readr)
library(dplyr)
library(stargazer)
library(normtest) #Prueba de Jarque-Bera
library(nortest) #Prueba KS y SW
library(haven) #para leer .dta
library(mctest) #Obtener indice de condicion
library(psych) #Para prueba de Farrar Glaudar
library(car) #Prueba de VIF y durbin watson
library(lmtest) #Prueba de White 
library(tidyr)
library(kableExtra)
library(foreign)
library(sandwich)

Pruebas de Normalidad

1.0 Leyendo datos

#si se desea hacer paso por paso ver "https://rpubs.com/Carlos_Martinez/501521" las 3 pruebas.
library(readr)
dataframe <- read.csv("C:/Users/ejhar/Desktop/econometria/ej2.csv")
colnames(dataframe) <- c("x1", "x2", "y")
head(dataframe, n=5)
##     x1   x2    y
## 1 3.92 7298 0.75
## 2 3.61 6855 0.71
## 3 3.32 6636 0.66
## 4 3.07 6506 0.61
## 5 3.06 6450 0.70

1.1 ¿Como estimar el modelo?

#Usamos lm
library(stargazer)
options(scipen=9999)
mod_lineal<-lm(formula = y~x1+x2, data= dataframe) #Así.
stargazer(mod_lineal, title = "Ejemplo de Regresion Multiple",
          type = "text", digits = 8)
## 
## Ejemplo de Regresion Multiple
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  y             
## -----------------------------------------------
## x1                         0.23719750***       
##                            (0.05555937)        
##                                                
## x2                        -0.00024908***       
##                            (0.00003205)        
##                                                
## Constant                   1.56449700***       
##                            (0.07939598)        
##                                                
## -----------------------------------------------
## Observations                    25             
## R2                          0.86529610         
## Adjusted R2                 0.85305030         
## Residual Std. Error    0.05330222 (df = 22)    
## F Statistic         70.66057000*** (df = 2; 22)
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

1.1.1 Prueba de Jarque Bera

library(normtest)
jb.norm.test(mod_lineal$residuals)
## 
##  Jarque-Bera test for normality
## 
## data:  mod_lineal$residuals
## JB = 0.93032, p-value = 0.477
#En caso que se pida graficar, mirar "https://rpubs.com/Carlos_Martinez/491076"

  • En caso de la prueba de Jarque-Bera la condicion de no rechazar de la Ho se puede evaluar por medio del p-value en la cual no se rechasa si p > \(\alpha\) y cuando el estadístico JB<V.C. JB(0.93032) < V.C(5.9915) p (0.4935) > \(\alpha\)(0.05)

  • NO SE RECHAZA LA HIPOTESIS NULA POR LO QUE HAY EVIDENCIA QUE LOS RESIDUOS SIGUEN UNA DISTRIBUCION NORMAL

1.1.2 Prueba de Kolmogorov Smirnov

library(nortest)
lillie.test(mod_lineal$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  mod_lineal$residuals
## D = 0.082345, p-value = 0.9328
#Para calcular n ver "http://www.real-statistics.com/statistics-tables/kolmogorov-smirnov-table/"
  • En caso de la prueba de Kolmogorov-Smirnov para un nivel de significancia del 5% y una muesta n=25 el V.C. = 0.26404 la condicion de no rechazar de la \(H_o\) es que el estadistico D < V.C., ademas tambien se puede evaluar por medio del p-value en la cual la condicion de no rechazo es p-value > \(\alpha\)

D (0.082345) < V.C.(0.159) p-value (0.9328) > \(\alpha\)(0.05)

  • NO SE RECHAZA LA HIPOTESIS NULA POR LO QUE HAY EVIDENCIA QUE LOS RESIDUOS SIGUEN UNA DISTRIBUCION NORMAL

1.1.3 Prueba de Shapiro Wilk

#se usa la misma libreria que la anterior prueba de KS
shapiro.test(mod_lineal$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod_lineal$residuals
## W = 0.97001, p-value = 0.6453

En caso de la prueba de Shapiro-Wilk para un nivel de significanciadel 5% el V.C. = 1.644854 la condicion de no rechazar de la \(H_o\) es que el estadistico SW < V.C., ademas tambien se puede evaluar por medio del p-value en la cual la condicion de no rechazo es p-value > \(\alpha\)

SW (0.97001) < V.C.(1.644854) p-value (0.6453) > \(\alpha\)(0.05)

PRUEBAS DE COLINEALIDAD

2.0 Calcular indice de condición.

Mide la sensibilidad de las estimaciones mínimo cuadráticos ante pequeños cambios en los datos.

Es necesario para esta prueba hacer lo siguiente:

  • Estimar el modelo.
  • Calcular |X’X|
  • Normalizar |X’X|
  • Calcular indice de condicion
#modelo estimado
model_es<-lm(y~x1+x2,data = dataframe)
stargazer(model_es,type = "text",title = "modelo estimado")
## 
## modelo estimado
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  y             
## -----------------------------------------------
## x1                           0.237***          
##                               (0.056)          
##                                                
## x2                          -0.0002***         
##                              (0.00003)         
##                                                
## Constant                     1.564***          
##                               (0.079)          
##                                                
## -----------------------------------------------
## Observations                    25             
## R2                             0.865           
## Adjusted R2                    0.853           
## Residual Std. Error       0.053 (df = 22)      
## F Statistic           70.661*** (df = 2; 22)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Calculo de |x’x|

#De un datafraeme estimado, así calcular X
library(stargazer)
X_mat<-model.matrix(model_es)
stargazer(head(X_mat,n=6),type="text")
## 
## =========================
##   (Intercept)  x1    x2  
## -------------------------
## 1      1      3.920 7,298
## 2      1      3.610 6,855
## 3      1      3.320 6,636
## 4      1      3.070 6,506
## 5      1      3.060 6,450
## 6      1      3.110 6,402
## -------------------------
#más fácil, imposible. :v
XX_matrix<-t(X_mat)%*%X_mat

Ahora, una vez tenemos la sigma matriz debemos normalizarla, tan simple como…

library(stargazer)
options(scipen = 999)
Sn<-solve(diag(sqrt(diag(XX_matrix)))) #sacamos raiz a la diagonal principal de la sigma matriz.
stargazer(Sn,type = "text")
## 
## ===================
## 0.200   0      0   
## 0     0.051    0   
## 0       0   0.00003
## -------------------
XX_norm<-(Sn%*%XX_matrix)%*%Sn #Sólo seguimos la fórmula estandar de normalización.
stargazer(XX_norm,type = "text",digits = 4)
## 
## ====================
## 1      0.9893 0.9909
## 0.9893   1    0.9988
## 0.9909 0.9988   1   
## --------------------

Ahora tan sencillo como calcular los autovalores, bueno con el comando “eigen” lo es.

#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)
stargazer(lambdas$values,type = "text")
## 
## =================
## 2.986 0.013 0.001
## -----------------
K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)
## [1] 51.03081

si k < 20 -> Multicolinealidad leve, no se considera un problema. si 20 < k < 30 -> Multicolinealidad moderada. si k > 30 -> Multicolinealidad severa.

3. Prueba de Farrar Glaubar

3.1 Aprendamos a calcular R y su determinante

En este caso se normalizará X y no |X’X| y luego segimos la formula, al final lo que nos interesa es el determinante de R

library(stargazer)
Zn<-scale(X_mat[,-1]) #Normalizando X
stargazer(head(Zn,n=6),type = "text")
## 
## ===============
##     x1     x2  
## ---------------
## 1 0.115  0.055 
## 2 -0.421 -0.387
## 3 -0.922 -0.605
## 4 -1.354 -0.735
## 5 -1.371 -0.791
## 6 -1.285 -0.839
## ---------------
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1)) #Formula de calculo.
#También se puede calcular R a través de cor(X_mat[,-1])
stargazer(R,type = "text",digits = 4)
## 
## ================
##      x1     x2  
## ----------------
## x1   1    0.9410
## x2 0.9410   1   
## ----------------
#VERIFICAR QUE R ESTA BIEN SI LA DIAGONAL ES IDENTIDAD.
#Determinante, calculo.
determinante_R<-det(R)
print(determinante_R)
## [1] 0.1145205

3.2 Estadistico \(X_{FG}^2\)

m<-ncol(X_mat[,-1])
n_2<-nrow(X_mat[,-1])
chi_FG<--(n_2-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)
## [1] 48.75753

3.3 Calculamos Valor critico

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

La condicion es la siguiente:

  • Si X > VC

Se rechaza la Ho

  • En este caso X es mayor, por lo que hay evidencia de colinealidad en los regresores

3.3.1 COMO PLUS, CALCULAR FARRAR CON libreria(psych)

# colocamos la matriz x dentro y ya.
library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
print(FG_test)
## $chisq
## [1] 48.75753
## 
## $p.value
## [1] 0.000000000002896446
## 
## $df
## [1] 1

la conclusion sera igual…. como p-value es muy pequeño ademas, es obvio que rechazamos Ho.

PRUEBA DE LOS FIV

Cuantifica la intensidad de multicolinealidad en un analisis de regresion normal de minimos cuadrados.

4.1 Calcular VIF para modelor estimado

extraemos de R y el solve(R)

inversa_R<-solve(R)
print(inversa_R)
##           x1        x2
## x1  8.732060 -8.216861
## x2 -8.216861  8.732060
#Si quieres saber como calcularla revisa de nuevo el punto 3.1 sn*|x'x|*sn

4.1.1 Para el modelo estimado

Los VIFs se calculan de la diagonal principal de la inversa de R

VIFs<-diag(inversa_R)
print(VIFs)
##      x1      x2 
## 8.73206 8.73206

Si VIF>2 presenta colinealidad Si VIF>5 son altamente colineales

4.2 Uso de la librería “car” y “mctest”

#A traves de car
library(car)
VIF_s<-vif(model_es)
print(VIF_s)
##      x1      x2 
## 8.73206 8.73206
#A traves de mctest
library(mctest)
mc.plot(x=X_mat[,-1], y= dataframe$y, vif = 2) #en Y seleccionamos la variable dependiente.

Si VIF>2 presenta colinealidad Si VIF>5 son altamente colineales

A.5 Heterocedasticidad

Se dan casos de heterocedasticidad cuando la varianza del termino de error no es constante

** Criterios de decisión **

A.5.1 calculo de la prueba de white paso a paso

library(stargazer)
u_i<-model_es$residuals #necesitamos los residuos
data_prueba_white<-as.data.frame(cbind(u_i,dataframe)) #solo juntamos en en una sola la data con los residuos 
regresion_auxiliar<-lm(I(u_i^2)~x1+x2+I(x1^2)+I(x2^2)+x1*x2,data = data_prueba_white)
sumario<-summary(regresion_auxiliar)
n2<-nrow(data_prueba_white)
R_2<-sumario$r.squared
LM_w<-n2*R_2
gl_w=2+2+1
p_value<-1-pchisq(q = LM_w,df = gl_w)
VC<-qchisq(p = 0.95,df = gl_w)
salida_white<-c(LM_w,VC,p_value)
names(salida_white)<-c("LMw","Valor Crítico","p value")
stargazer(salida_white,title = "Resultados de la prueba de White",type = "text",digits = 6)
## 
## Resultados de la prueba de White
## ===============================
## LMw      Valor Crítico p value 
## -------------------------------
## 3.690182   11.070500   0.594826
## -------------------------------

Como p-value (0.5948) > 0.05 No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedastica.

A.5.1.1 Calculo de la prueba de white libreria(lmtest)

library(lmtest)
#bptest de breush pagan
prueba_white<-bptest(model_es, ~I(x1^2)+I(x2^2)+x1*x2, data = dataframe)
print(prueba_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_es
## BP = 3.6902, df = 5, p-value = 0.5948

Como p-value (0.5948) > 0.05 No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedastica.

B.5 Autocorrelacion

Los valores del pasado tienen alta asociacion con los del presente

B.5.1 Autocorrelacion de 1er grado - Durbin Watson

Usando libreria (lmtest)

library(lmtest)
dwtest(model_es,alternative = "two.sided",iterations = 1000)
## 
##  Durbin-Watson test
## 
## data:  model_es
## DW = 1.9483, p-value = 0.5649
## alternative hypothesis: true autocorrelation is not 0

usando libreria (car)

library(car)
durbinWatsonTest(model_es, simulate = TRUE, reps = 1000)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.04366918      1.948305   0.594
##  Alternative hypothesis: rho != 0

En ambos casos, se puede rechazar la presencia de autocorrelación (No se rechaza la H0), ya que el pvalue > 0.05

B.5.2 Autocorrelacion de orden superior - Breush Godfrey

  1. Preparación de datos
library(dplyr)
library(tidyr)
library(kableExtra)
cbind(u_i,dataframe) %>% 
  as.data.frame() %>% 
  mutate(Lag_1=dplyr::lag(u_i,1),
         Lag_2=dplyr::lag(u_i,2)) %>% 
  replace_na(list(Lag_1=0,Lag_2=0))->data_prueba_BG
kable(head(data_prueba_BG,6))
u_i x1 x2 y Lag_1 Lag_2
0.0734697 3.92 7298 0.75 0.0000000 0.0000000
-0.0033412 3.61 6855 0.71 0.0734697 0.0000000
-0.0391023 3.32 6636 0.66 -0.0033412 0.0734697
-0.0621832 3.07 6506 0.61 -0.0391023 -0.0033412
0.0162403 3.06 6450 0.70 -0.0621832 -0.0391023
0.0124247 3.11 6402 0.72 0.0162403 -0.0621832
  1. Calculando la regresión auxiliar y el estadístico LMBP
regresion_auxiliar_BG<-lm(u_i~x1+x2+Lag_1+Lag_2,data = data_prueba_BG)
sumario_BG<-summary(regresion_auxiliar_BG)
R_2_BG<-sumario_BG$r.squared
n<-nrow(data_prueba_BG)
LM_BG<-n*R_2_BG
gl=2
p_value<-1-pchisq(q = LM_BG,df = gl)
VC<-qchisq(p = 0.95,df = gl)
salida_bg<-c(LM_BG,VC,p_value)
names(salida_bg)<-c("LMbg","Valor Crítico","p value")
stargazer(salida_bg,title = "Resultados de la prueba de Breusch Godfrey",type = "text",digits = 6)
## 
## Resultados de la prueba de Breusch Godfrey
## ===============================
## LMbg     Valor Crítico p value 
## -------------------------------
## 3.305189   5.991465    0.191552
## -------------------------------

Como pvalue>0.05 No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden “2”.

B.5.2 Autocorrelación por librerias

Usando lmtest

library(lmtest)
bgtest(model_es,order = 2) #si quisiera que fuera de orden 1 sólo cambio 2 por 1 y se lee igual el resultado.
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  model_es
## LM test = 3.3052, df = 2, p-value = 0.1916

Como pvalue>0.05 No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden “2”

6. Pruebas HAC

Primero se estima el modelo y ya… ## 6.1 Estimacion robusta En caso que se detecte heterocedasticidad

library(lmtest)
library(sandwich)
#Sin corregir:

coeftest(model_es)
## 
## t test of coefficients:
## 
##                 Estimate   Std. Error t value             Pr(>|t|)    
## (Intercept)  1.564496771  0.079395981 19.7050 0.000000000000001817 ***
## x1           0.237197475  0.055559366  4.2693            0.0003126 ***
## x2          -0.000249079  0.000032048 -7.7719 0.000000095087907942 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Corregido:

estimacion_omega<-vcovHC(model_es,type = "HC1")
coeftest(model_es,vcov. = estimacion_omega)
## 
## t test of coefficients:
## 
##                 Estimate   Std. Error t value            Pr(>|t|)    
## (Intercept)  1.564496771  0.092172938 16.9735 0.00000000000003986 ***
## x1           0.237197475  0.049431319  4.7985 0.00008591053735281 ***
## x2          -0.000249079  0.000033313 -7.4769 0.00000017801116050 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.2 Estimacion robusta

En caso de autocorrelacion

library(lmtest)
library(sandwich)

#Corregido:

estimacion_omega<-NeweyWest(model_es,lag = 2)
coeftest(model_es,vcov. = estimacion_omega)
## 
## t test of coefficients:
## 
##                 Estimate   Std. Error  t value              Pr(>|t|)    
## (Intercept)  1.564496771  0.068445738  22.8575 < 0.00000000000000022 ***
## x1           0.237197475  0.040049244   5.9226        0.000005838330 ***
## x2          -0.000249079  0.000024714 -10.0783        0.000000001047 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

NOTA: Si se detecta tan autocorrelacion como heterocedasticidad la ultima prueba vale para corregir ambos casos.