homocedasticidad y autocorrelacion
correr el modelo
library(readr)
library(stargazer)
ejemplo_regresion <- read_csv("C:/Users/Kathya Hernandez/Downloads/ejemplo_regresion.csv")
modelo_lineal<-lm(Y~X1+X2,data = ejemplo_regresion)
stargazer(modelo_lineal,title = "modelo estimado",type = "html")
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
|
ejemplo de la prueba white manual
library(stargazer)
u_i<-modelo_lineal$residuals
data_prueba_white<-as.data.frame(cbind(u_i,ejemplo_regresion))
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)
n<-nrow(data_prueba_white)
R_2<-sumario$r.squared
LM_w<-n*R_2
gl=2+2+1
p_value<-1-pchisq(q = LM_w,df = gl)
VC<-qchisq(p = 0.95,df = gl)
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 = "html",digits = 6)
Resultados de la prueba de White
|
|
|
LMw
|
Valor Crítico
|
p value
|
|
|
|
3.690182
|
11.070500
|
0.594826
|
|
|
Como 0.5948>0.05, no se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica
Uso de la libreria lmtest
library(lmtest)
prueba_white<-bptest(modelo_lineal,~I(X1^2)+I(X2^2)+X1*X2,data = ejemplo_regresion)
print(prueba_white)
##
## studentized Breusch-Pagan test
##
## data: modelo_lineal
## BP = 3.6902, df = 5, p-value = 0.5948
Como 0.5948>0.05 No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica.
autocorrelacion
correlacion de primer orden
prueba Durbin Watson
# usando lmtest
library(lmtest)
dwtest(modelo_lineal,alternative = "two.sided",iterations = 1000)
##
## Durbin-Watson test
##
## data: modelo_lineal
## DW = 1.9483, p-value = 0.5649
## alternative hypothesis: true autocorrelation is not 0
# usando libreria car
library(car)
durbinWatsonTest(modelo_lineal,simulate = TRUE,reps = 1000)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.04366918 1.948305 0.544
## 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
autocorrelacion de orden superior
**prueba del multiplicador de lagrange [Breusch-Godfrey]
manual
library(dplyr)
library(tidyr)
library(kableExtra)
cbind(u_i,ejemplo_regresion) %>%
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
|
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 = "html",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”.
usando la libreria lmtest
library(lmtest)
bgtest(modelo_lineal,order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_lineal
## 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”
El test BG puede usarse también para verificar la autocorrelación de 1° orden:
library(lmtest)
bgtest(modelo_lineal,order = 1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_lineal
## LM test = 0.051063, df = 1, p-value = 0.8212
Como pvalue>0.05 No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de 1° orden
Estimadores HAC
corriendo el modelo
options(scipen = 999999)
library(foreign)
datos_regresion <- read.dta("C:/Users/Kathya Hernandez/Downloads/crime.dta")
modelo_crime<-lm(crime~poverty+single,data=datos_regresion)
print(modelo_crime)
##
## Call:
## lm(formula = crime ~ poverty + single, data = datos_regresion)
##
## Coefficients:
## (Intercept) poverty single
## -1368.189 6.787 166.373
realizar las pruebas de heterocedasticidad y autocorrelacion
heterocedasticidad
Prueba de White (prueba de Breusch Pagan)
library(lmtest)
white_test<-bptest(modelo_crime,~I(poverty^2)+I(single^2)+poverty*single,data = datos_regresion)
print(white_test)
##
## studentized Breusch-Pagan test
##
## data: modelo_crime
## BP = 10.73, df = 5, p-value = 0.057
hay evidencia de heterocedasticidad ya que Pvalue<0.05 (por que el pvalue es muy pequeño)
Prueba del Multiplicador de Lagrange (Breusch Godfrey)
Autocorrelación de 2º Orden:
library(lmtest)
prueba_LM<-bgtest(modelo_crime,order = 2)
print(prueba_LM)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_crime
## LM test = 0.27165, df = 2, p-value = 0.873
no hay evidencia de autocorrelacion de orden 2 ya que pvalue>0.05
Autocorrelación de 1º orden (prueba de Durbin Watson)
library(car)
durbinWatsonTest(model = modelo_crime)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.07014421 2.040007 0.96
## Alternative hypothesis: rho != 0
no hay evidencia de autocorrelacion de orden 1 ya que pvalue>0.05
Estimación Robusta (uso del estimador HAC)
corrigiendo la homocedasticidad
sin corregir
options(scipen = 99999)
library(lmtest)
#Sin corregir:
coeftest(modelo_crime)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 187.2052 -7.3085 0.00000000247861 ***
## poverty 6.7874 8.9885 0.7551 0.4539
## single 166.3727 19.4229 8.5658 0.00000000003117 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
corrigiendo (usando un estimador HAC)
options(scipen = 99999)
library(lmtest)
library(sandwich)
#Corregido
#HC0 Corrige Sólo Heterocedasticidad, use HC1 para corregir también Autocorrelación de Primer Orden
estimacion_omega<-vcovHC(modelo_crime,type = "HC0")
coeftest(modelo_crime,vcov. = estimacion_omega)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 276.4111 -4.9498 0.00000956181 ***
## poverty 6.7874 10.6010 0.6403 0.5251
## single 166.3727 25.4510 6.5370 0.00000003774 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estimador HAC de NeweyWest
Nota:(si la correlación es de orden 1, es más simple usar el estimador HC1)
Si se hubiera detectado la correlación de segundo orden se tendría que corregir así:
library(lmtest)
library(sandwich)
#Corregido:
estimacion_omega<-NeweyWest(modelo_crime,lag = 2)
coeftest(modelo_crime,vcov. = estimacion_omega)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 303.8466 -4.5029 0.00004279768 ***
## poverty 6.7874 10.5943 0.6407 0.5248
## single 166.3727 25.9154 6.4198 0.00000005708 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estimación Robusta
options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_crime_robust<-lmrob(crime~poverty+single,data=datos_regresion)
# print(summary(modelo_crime_robust))
stargazer(modelo_crime,modelo_crime_robust,type = "html",title = "comparativa")
comparativa
|
|
|
|
Dependent variable:
|
|
|
|
|
|
crime
|
|
|
OLS
|
MM-type
|
|
|
|
linear
|
|
|
(1)
|
(2)
|
|
|
|
poverty
|
6.787
|
11.466
|
|
|
(8.989)
|
(9.263)
|
|
|
|
|
|
single
|
166.373***
|
176.569***
|
|
|
(19.423)
|
(23.223)
|
|
|
|
|
|
Constant
|
-1,368.189***
|
-1,539.640***
|
|
|
(187.205)
|
(235.765)
|
|
|
|
|
|
|
|
Observations
|
51
|
51
|
|
R2
|
0.707
|
0.795
|
|
Adjusted R2
|
0.695
|
0.787
|
|
Residual Std. Error (df = 48)
|
243.610
|
191.864
|
|
F Statistic
|
57.964*** (df = 2; 48)
|
|
|
|
|
Note:
|
p<0.1; p<0.05; p<0.01
|
ejemplo simulacion
#Bias Proportion
Um<-function(pronosticado,observado){
library(DescTools)
((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado)
}
#Variance Proportion
Us<-function(pronosticado,observado){
library(DescTools)
((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc<-function(pronosticado,observado){
library(DescTools)
(2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U<-function(pronosticado,observado){
library(DescTools)
RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}
script de simulacion
options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(50)
numero_de_muestras<-1000
muestras<- BostonHousing$medv %>%
createDataPartition(p = 0.8,
times = numero_de_muestras,
list = TRUE)
Modelos_Entrenamiento<-vector(mode = "list",
length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance<-vector(mode = "list",
length = numero_de_muestras)
#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- BostonHousing[muestras[[j]], ]
Datos_Prueba<- BostonHousing[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = medv~.,data=Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]]<-data.frame(
R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv),
RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv),
MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv),
MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv)*100,
THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv,type = 1),
Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv),
Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv),
Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
Datos_Entrenamiento$medv)
)
Resultados_Performance[[j]]<-data.frame(
R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$medv),
RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$medv)*100,
THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$medv,
type = 1),
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
)
}
resultados de la simulacion
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "html",
digits = 3)
Medidas de Performance Datos del Modelo
|
|
|
Statistic
|
N
|
Mean
|
St. Dev.
|
Min
|
Pctl(25)
|
Pctl(75)
|
Max
|
|
|
|
R2
|
1,000
|
0.743
|
0.013
|
0.713
|
0.734
|
0.751
|
0.794
|
|
RMSE
|
1,000
|
4.653
|
0.141
|
4.177
|
4.565
|
4.759
|
4.948
|
|
MAE
|
1,000
|
3.265
|
0.095
|
2.905
|
3.204
|
3.332
|
3.512
|
|
MAPE
|
1,000
|
16.387
|
0.464
|
14.813
|
16.085
|
16.718
|
17.691
|
|
THEIL
|
1,000
|
0.096
|
0.003
|
0.087
|
0.095
|
0.099
|
0.102
|
|
Um
|
1,000
|
0.000
|
0.000
|
0
|
0
|
0
|
0
|
|
Us
|
1,000
|
0.074
|
0.004
|
0.058
|
0.072
|
0.077
|
0.085
|
|
Uc
|
1,000
|
0.928
|
0.004
|
0.918
|
0.925
|
0.931
|
0.945
|
|
|
bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "html",
digits = 3)
Medidas de Performance Simulación
|
|
|
Statistic
|
N
|
Mean
|
St. Dev.
|
Min
|
Pctl(25)
|
Pctl(75)
|
Max
|
|
|
|
R2
|
1,000
|
0.723
|
0.056
|
0.452
|
0.690
|
0.764
|
0.840
|
|
RMSE
|
1,000
|
4.862
|
0.575
|
3.465
|
4.450
|
5.226
|
6.961
|
|
MAE
|
1,000
|
3.411
|
0.281
|
2.633
|
3.216
|
3.598
|
4.492
|
|
MAPE
|
1,000
|
17.197
|
1.618
|
12.875
|
16.066
|
18.262
|
23.137
|
|
THEIL
|
1,000
|
0.101
|
0.012
|
0.073
|
0.092
|
0.108
|
0.148
|
|
Um
|
1,000
|
0.011
|
0.016
|
0.000
|
0.001
|
0.015
|
0.205
|
|
Us
|
1,000
|
0.081
|
0.066
|
0.00000
|
0.027
|
0.122
|
0.333
|
|
Uc
|
1,000
|
0.918
|
0.066
|
0.667
|
0.875
|
0.971
|
1.010
|
|
|