library(wooldridge)
data(hprice1)
head(force(hprice1),n=5)
##   price assess bdrms lotsize sqrft colonial   lprice  lassess llotsize   lsqrft
## 1   300  349.1     4    6126  2438        1 5.703783 5.855359 8.720297 7.798934
## 2   370  351.5     3    9903  2076        1 5.913503 5.862210 9.200593 7.638198
## 3   191  217.7     3    5200  1374        0 5.252274 5.383118 8.556414 7.225482
## 4   195  231.8     3    4600  1448        1 5.273000 5.445875 8.433811 7.277938
## 5   373  319.1     4    6095  2514        1 5.921578 5.765504 8.715224 7.829630

1. estime el siguiente modelo:

price = ˆα + ˆα1(lotsize) + ˆα2(sqrft) + ˆα3(bdrms) + ∈

library(stargazer)

options(scipen = 9999)
Modeloautocorrelacion <- lm(formula = price~lotsize+sqrft+ bdrms, data = hprice1)
stargazer(Modeloautocorrelacion, title = "Modelo Estimado", type = "text", digits = 6)
## 
## Modelo Estimado
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                price           
## -----------------------------------------------
## lotsize                     0.002068***        
##                             (0.000642)         
##                                                
## sqrft                       0.122778***        
##                             (0.013237)         
##                                                
## bdrms                        13.852520         
##                             (9.010145)         
##                                                
## Constant                    -21.770310         
##                             (29.475040)        
##                                                
## -----------------------------------------------
## Observations                    88             
## R2                           0.672362          
## Adjusted R2                  0.660661          
## Residual Std. Error     59.833480 (df = 84)    
## F Statistic          57.460230*** (df = 3; 84) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

2. Verifique si los residuos del modelo son independientes entre sí (no autocorrelación), a través de:

a) Prueba de Durbin Watson.

Usando la libreria de “lmtest”

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

No se rechaza la H0, porque el valor p = 0.6218 es mayor que 0.05

Usando la libreria de “car”

library(car)

durbinWatsonTest(Modeloautocorrelacion, simulate= TRUE, reps=1000)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.05900522      2.109796   0.612
##  Alternative hypothesis: rho != 0

En este modelo no hay evidecia de la presencia de autocorrelación

B) Prueba del Multiplicador de Lagrange (verifique autocorrelación de primer y segundo orden).

library(dplyr)
library(tidyr)
library(kableExtra)


Ui <- Modeloautocorrelacion$residuals
Mat_x <- model.matrix(Modeloautocorrelacion)
MatInformación <- Mat_x[,-1]  
Endogena <- hprice1$price
cbind(Ui, Endogena, MatInformación) %>%
  as.data.frame()%>%
  mutate(Lag1 = dplyr::lag(Ui, 1),
         Lag2 = dplyr::lag(Ui, 2))%>%
  replace_na(list(Lag1=0,Lag2=0)) -> Data.PruebaBG
#kable(head(x = Data.PruebaBG, n = 7))

  ## Cálculo de la regresión auxiliar:
Regresion.AuxiliarBG <- lm(Ui~lotsize+sqrft+bdrms + Lag1 + Lag2, data = Data.PruebaBG)
ResumenBG <- summary(Regresion.AuxiliarBG)

  ## Cálculo del estadístico LMbg:
R2 <- ResumenBG$r.squared
n <- nrow(Data.PruebaBG)
LMbg <- n*R2

    ## Sacamos los grados de libertad:
gl <- 2 #xq se esta verificando autocorrelación de orden dos (Lag1+Lag2+...+Lagn= gl)

    ## Sacamos el P Value:
P_Value <- 1-pchisq(q = LMbg, df = gl)

    ## Sacamos Valor Critico:
VC <- qchisq(p = 0.05,df = gl,lower.tail = FALSE)
#VC <- qchisq(p = 0.95,df = gl,lower.tail = TRUE)

    ## Presentamos los resultados:
library(stargazer)
SalidaBG <- c(LMbg, VC,P_Value)
names(SalidaBG) <- c("LMbg", "Valor Crítico","P Value")
stargazer(SalidaBG, title = "Resultado de la prueba de Breusch Godfrey", type = "text", digits = 6 )
## 
## Resultado de la prueba de Breusch Godfrey
## ===============================
## LMbg     Valor Crítico P Value 
## -------------------------------
## 3.033403   5.991465    0.219435
## -------------------------------

Usando libreria lmtest

Autocorrelación de primer orden:

library(lmtest)
PruebaMultiplicadorLagrange1 <- bgtest(Modeloautocorrelacion, order = 1)
PruebaMultiplicadorLagrange1
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  Modeloautocorrelacion
## LM test = 0.39362, df = 1, p-value = 0.5304

Como el p-value (0.05304) > 0.05 no se rechaza la H0, por lo tanto se puede concluir que los residuos que tiene el modelo no siguen autocorrelación de orden “uno”

Autocorrelación de segundo orden:

library(lmtest)
PruebaMultiplicadorLagrange2 <- bgtest(Modeloautocorrelacion, order = 2)
PruebaMultiplicadorLagrange2
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  Modeloautocorrelacion
## LM test = 3.0334, df = 2, p-value = 0.2194

Como el p-value (0.2194) > 0.05 no se rechaza H0, por lo tanto se puede concluir que los residuos que tiene el modelo no siguen autocorrelación de orden “dos”