library(readr)
ejemplo_regresion <- read_csv("C:/Users/Eduardo/Downloads/ejemplo_regresion.csv")
head(ejemplo_regresion, n=6)
## # A tibble: 6 x 3
## X1 X2 Y
## <dbl> <dbl> <dbl>
## 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.7
## 6 3.11 6402 0.72
library(stargazer)
modelo_ejemplo<- lm(formula = Y~X1+X2, data = ejemplo_regresion)
stargazer(modelo_ejemplo, title = "Ejemplo Regresion Heterocedasticidad", type= "text")
##
## Ejemplo Regresion Heterocedasticidad
## ===============================================
## 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
Verificando el supuesto de Heterocedasticidad (es donde la varianza del error no es la misma para todos los niveles de las variables explicativas, consistia en que nuestros interbalos de confianza no eran robustos)
Prueba de White (creando una regresion auxiliar)
library(stargazer)
resid<-modelo_ejemplo$residuals
data_regre_auxi<-as.data.frame(cbind(resid,ejemplo_regresion))
regresion_auxiliar<-lm(I(resid^2)~X1+X2+I(X1^2)+I(X2^2)+X1*X2,data = data_regre_auxi)
resumen<-summary(regresion_auxiliar)
R_2<-resumen$r.squared
n<-nrow(data_regre_auxi)
LM_w<-n*R_2
gl<-2+2+1
VC<-qchisq(p=0.95,df=gl)
p_value<-1-pchisq(q=LM_w,df=gl)
salida_prueba<-c(LM_w,VC,p_value)
names(salida_prueba)<-c("LMw","Valor Critico","P Value")
stargazer(salida_prueba,title = "Prueba White",type = "html",digits = 6)
Prueba White
|
|
|
LMw
|
Valor Critico
|
P Value
|
|
|
|
3.690182
|
11.070500
|
0.594826
|
|
|
Prueba White con libreria lmtest
library(lmtest)
prueba_white<-bptest(modelo_ejemplo,~I(X1^2)+I(X2^2)+X1*X2,data=ejemplo_regresion)
print(prueba_white)
##
## studentized Breusch-Pagan test
##
## data: modelo_ejemplo
## BP = 3.6902, df = 5, p-value = 0.5948
# Concluimos que no hay Heterocedasticidad
Autocorrelacion (lo que se trata de evaluar es que si los residuos del modelo tienen asociacion con sus mismos valores resagados o retrdados, como que si el error presente es multiple del error pasado )
Prueba de Durbin – Watson con libreria lmtest
library(lmtest)
dwtest(modelo_ejemplo,alternative = "two.sided",iterations = 1000)
##
## Durbin-Watson test
##
## data: modelo_ejemplo
## DW = 1.9483, p-value = 0.5649
## alternative hypothesis: true autocorrelation is not 0
# Al tener un p-value de este tamaño podemos ver que no hay autocorrelacion y el estadistico de prueba es DW = 1.9484 y podemos decir que no existe autocorrelacion
Prueba de Durbin – Watson con libreria car
library(car)
durbinWatsonTest(modelo_ejemplo,simulate = TRUE,reps = 1000)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.04366918 1.948305 0.576
## Alternative hypothesis: rho != 0
# Cloncluimos que no hay autocorrelacion de primer orden
Autocorrelacion de Orden Superior (orden “m”)
library(dplyr)
library(tidyr)
library(stargazer)
library(kableExtra)
cbind(resid,ejemplo_regresion)%>% as.data.frame()%>%
mutate(Lag_1=dplyr::lag(resid,1),
Lag_2=dplyr::lag(resid,2)) %>%
replace_na(list(Lag_1=0,Lag_2=0))-> data_prueba_BG
kable(head(data_prueba_BG,n=6))
|
resid
|
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
|
regresion_auxi_BG<-lm(resid~X1+X2+Lag_1+Lag_2,data = data_prueba_BG)
resumen_BG<-summary(regresion_auxi_BG)
n<-nrow(ejemplo_regresion)
gl<-2
R_2_BG<-resumen_BG$r.squared
LM_BG<-n*R_2_BG
VC<-qchisq(p=0.05,lower.tail = FALSE,df=gl)
p_value_BG<-pchisq(q=LM_BG,lower.tail = FALSE,df=gl)
salida_BG<-c(LM_BG,VC,p_value_BG)
names(salida_BG)<-c("LM BG","Valor Critico","P value BG")
stargazer(salida_BG,title = "Prueba de BG",type = "html",digits = 6)
Prueba de BG
|
|
|
LM BG
|
Valor Critico
|
P value BG
|
|
|
|
3.305189
|
5.991465
|
0.191552
|
|
|
# No se rechaza la H0
Usando lmtest
library(lmtest)
bgtest(modelo_ejemplo,order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_ejemplo
## LM test = 3.3052, df = 2, p-value = 0.1916