Ejercicio 1. Sea el conjunto de datos, indicados en el enlace de abajo, tomados en 24 meses correspondientes a los gastos de comercialización (C) de una empresa, el nivel de ventas (V), su coste de personal (P) y los costes de materias primas (M); se trata de estimar el nivel de ventas a partir de las restantes variables.
Estimación del modelo
library(readxl)
library(stargazer)
ventas_empresa_1_ <- read_excel("C:/Users/familia/Downloads/ventas_empresa (1).xlsx")
modelo_ventas<-lm(formula= V~C+P+M, data = ventas_empresa_1_)
stargazer(modelo_ventas, title = "Modelo Laboratorio 1", digits = 6, type = "text")
##
## Modelo Laboratorio 1
## ===============================================
## Dependent variable:
## ---------------------------
## V
## -----------------------------------------------
## C 0.922567***
## (0.222733)
##
## P 0.950177***
## (0.155845)
##
## M 1.297786***
## (0.430729)
##
## Constant 107.443500***
## (18.057490)
##
## -----------------------------------------------
## Observations 24
## R2 0.979817
## Adjusted R2 0.976789
## Residual Std. Error 9.505570 (df = 20)
## F Statistic 323.641500*** (df = 3; 20)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
1. Verifique los supuestos de Heterocedastidad y Autocorrelación para el modelo propuesto.
Prueba de White
library(stargazer)
resid<-modelo_ventas$residuals
data_reg_aux<-as.data.frame(cbind(resid,ventas_empresa_1_))
reg_aux<-lm(formula = I(resid^2)~C+P+M+I(C^2)+I(P^2)+I(M^2)+C*P+C*M+M*P,
data = data_reg_aux)
resumen<-summary(reg_aux)
R_2<-resumen$r.squared
n<-nrow(data_reg_aux)
LM_w<-n*R_2
gl<-3+3+3
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 Crítico", "P value")
stargazer(salida_prueba, title="Prueba de White", type="text", digits=6)
##
## Prueba de White
## ===============================
## LMw Valor Crítico P value
## -------------------------------
## 7.122650 16.918980 0.624351
## -------------------------------
Como pvalue=0.6243 >0.05 No se rechaza la H0, por lo tanto hay evidencia que la varianza de los residuos es homocedástica.
Prueba de White con libreria lmtest
library(lmtest)
prueba_white<-bptest(modelo_ventas,~I(C^2)+I(P^2)+I(M^2)+C*P+C*M+M*P, data = ventas_empresa_1_)
print(prueba_white)
##
## studentized Breusch-Pagan test
##
## data: modelo_ventas
## BP = 7.1227, df = 9, p-value = 0.6244
Como pvalue=0.6243 >0.05 No se rechaza la H0, por lo tanto hay evidencia que la varianza de los residuos es homocedástica.
Autocorrelación de 1er Orden
Prueba de Durbin-Watson
Usando libreria “lmtest”
library(lmtest)
dwtest(modelo_ventas, alternative = "two.sided", iterations = 1000)
##
## Durbin-Watson test
##
## data: modelo_ventas
## DW = 1.2996, p-value = 0.05074
## alternative hypothesis: true autocorrelation is not 0
Usando libreria “car”
library(car)
durbinWatsonTest(modelo_ventas, simulate = TRUE, reps = 1000)
## lag Autocorrelation D-W Statistic p-value
## 1 0.3013888 1.299572 0.046
## Alternative hypothesis: rho != 0
En ambos casos, al utilizar la tabla del estadístico Durbin-Watson 1.101<DW<1.656 cae en la zona no concluyente y el pvalue es cercano al nivel de significancia, por lo tanto se aplicará el test de BG para verificar la autocorrelación de primer orden.
Test de BG para verificar autocorrelación de primer orden
library(lmtest)
bgtest(modelo_ventas, order = 1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_ventas
## LM test = 2.5963, df = 1, p-value = 0.1071
Como pvalue>0.05 no se rechaza H0, por lo tanto se concluye que los residuos del modelo, no siguen autocorrelación de 1er orden.
Autocorrelación de orden superior (orden “m”)
library(dplyr)
library(tidyr)
library(stargazer)
library(kableExtra)
cbind(resid,ventas_empresa_1_) %>%
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
|
V
|
C
|
P
|
M
|
Lag_1
|
Lag_2
|
|
10.673678
|
607
|
197
|
173
|
110
|
0.000000
|
0.000000
|
|
7.372511
|
590
|
208
|
152
|
107
|
10.673678
|
0.000000
|
|
-2.435532
|
543
|
181
|
150
|
99
|
7.372511
|
10.673678
|
|
-3.322264
|
558
|
194
|
150
|
102
|
-2.435532
|
7.372511
|
|
-9.913932
|
571
|
192
|
163
|
109
|
-3.322264
|
-2.435532
|
|
8.704039
|
615
|
196
|
179
|
114
|
-9.913932
|
-3.322264
|
reg_aux_BG<-lm(resid~C+P+M+Lag_1+Lag_2, data = data_prueba_BG)
resumen_BG<-summary(reg_aux_BG)
n1<-nrow(data_prueba_BG)
gl1<-2
R_2_BG<-resumen_BG$r.squared
LM_BG<-n1*R_2_BG
VC<-qchisq(p=0.95, df=gl1)
p_value_BG<-pchisq(q=LM_BG, lower.tail = FALSE, df=gl1)
salida_BG<-c(LM_BG,VC,p_value_BG)
names(salida_BG)<-c("LM_BG","Valor Crítico", "Valor P")
stargazer(salida_BG, title = "Prueba de BG", type = "text", digits = 6)
##
## Prueba de BG
## ===============================
## LM_BG Valor Crítico Valor P
## -------------------------------
## 3.840869 5.991465 0.146543
## -------------------------------
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 libreria “lmtest”
library(lmtest)
bgtest(modelo_ventas, order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_ventas
## LM test = 3.8409, df = 2, p-value = 0.1465
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”.
2. En caso de encontrar evidencia de violación de los supuestos, planteados en el literal anterior, corrija a través de un estimador HAC apropiado, el modelo propuesto. Comentario: El modelo propuesto no viola ninguno de los supuestos, no hay evidencia de Heterocedasticidad y Autocorrelación por lo tanto no se procede a corregir a través de los estimadores HAC.
Ejercicio 2. Se tienen los datos para trabajadores hombres,en el archivo adjunto, con ellos estime un modelo donde educ es años de escolaridad, como variable dependiente, y como regresores sibs (número de hermanos), meduc (años de escolaridad de la madre) y feduc (años de escolaridad del padre)
Estimación del modelo
library(stargazer)
load("D:/Arturo - Rstudio/wage2.RData")
modelo_escolar<-lm(formula = educ~sibs+meduc+feduc, data= wage2)
stargazer(modelo_escolar, title = "Modelo Laboratorio 2", type="text")
##
## Modelo Laboratorio 2
## ===============================================
## Dependent variable:
## ---------------------------
## educ
## -----------------------------------------------
## sibs -0.094***
## (0.034)
##
## meduc 0.131***
## (0.033)
##
## feduc 0.210***
## (0.027)
##
## Constant 10.364***
## (0.359)
##
## -----------------------------------------------
## Observations 722
## R2 0.214
## Adjusted R2 0.211
## Residual Std. Error 1.987 (df = 718)
## F Statistic 65.198*** (df = 3; 718)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
1. Verifique los supuestos de Heterocedastidad y Autocorrelación para el modelo propuesto.
Prueba white con libreria lmtest
library(lmtest)
prueba_white1<-bptest(modelo_escolar,~I(sibs^2)+I(meduc^2)+I(feduc^2)+sibs*meduc+sibs*feduc+meduc*feduc, data = wage2)
print(prueba_white1)
##
## studentized Breusch-Pagan test
##
## data: modelo_escolar
## BP = 15.537, df = 9, p-value = 0.0772
Como pvalue=0.0772 >0.05 No se rechaza la H0, por lo tanto hay evidencia que la varianza de los residuos es homocedástica.
Autocorrelación de 1er Orden
Prueba de Durbin-Watson
Usando libreria “lmtest”
library(lmtest)
dwtest(modelo_escolar, alternative = "two.sided", iterations = 1000)
##
## Durbin-Watson test
##
## data: modelo_escolar
## DW = 1.8989, p-value = 0.1705
## alternative hypothesis: true autocorrelation is not 0
Usando libreria “car”
library(car)
durbinWatsonTest(modelo_escolar, simulate = TRUE, reps = 1000)
## lag Autocorrelation D-W Statistic p-value
## 1 0.05018452 1.898938 0.17
## Alternative hypothesis: rho != 0
En ambos casos, no hay evidencia de Autocorrelación de 1er orden en los residuos del modelo, ya que el Pvalue>0.05, por lo tanto no se rechaza la H0.
Autocorrelación de orden superior (orden “m”)
Usando libreria “lmtest”
library(lmtest)
bgtest(modelo_escolar, order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_escolar
## LM test = 4.5747, df = 2, p-value = 0.1015
Como el pvalue>0.05 No se rechaza la H0 por lo tantom los residuos del modelo no siguen autocorrelación de orden “2”
2. En caso de encontrar evidencia de violación de los supuestos, planteados en el literal anterior, corrija a través de un estimador HAC apropiado, el modelo propuesto. Comentario: El modelo propuesto no viola ninguno de los supuestos, por lo tanto no hay evidencia de Heterocedasticidad y Autocorrelación por lo tanto no se procede a corregir a través de los estimadores HAC.
Ejercicio 3. El sueldo inicial medio (salary) para los recién graduados de la Facultad de Economía se determina mediante una función lineal: log(salary)=f(SAT,GPA,log(libvol),log(cost),rank). Donde LSAT es la media del puntaje LSAT del grupo de graduados, GPA es la media del GPA (promedio general) del grupo, libvol es el número de volúmenes en la biblioteca de la Facultad de Economía, cost es el costo anual por asistir a dicha facultad y rank es una clasificación de las escuelas de Economía (siendo rank 1 la mejor)
Estimación del Modelo
library(stargazer)
load("D:/Arturo - Rstudio/LAWSCH85.RData")
modelo_salario<-lm(formula = lsalary~LSAT+GPA+llibvol+lcost+rank, data = LAWSCH85)
stargazer(modelo_salario, title = "Modelo de Salario", type = "text")
##
## Modelo de Salario
## ===============================================
## Dependent variable:
## ---------------------------
## lsalary
## -----------------------------------------------
## LSAT 0.005
## (0.004)
##
## GPA 0.248***
## (0.090)
##
## llibvol 0.095***
## (0.033)
##
## lcost 0.038
## (0.032)
##
## rank -0.003***
## (0.0003)
##
## Constant 8.343***
## (0.533)
##
## -----------------------------------------------
## Observations 136
## R2 0.842
## Adjusted R2 0.836
## Residual Std. Error 0.112 (df = 130)
## F Statistic 138.230*** (df = 5; 130)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
1. Verifique los supuestos de Heterocedastidad y Autocorrelación para el modelo propuesto.
Heterocedasticidad
Prueba de White con libreria “lmtest”
library(lmtest)
prueba_white2<-bptest(modelo_salario,~I(LSAT^2)+I(GPA^2)+I(llibvol^2)+
I(lcost^2)+I(rank^2)+LSAT*GPA+LSAT*llibvol+LSAT*lcost+
LSAT*rank+GPA*llibvol+GPA*lcost+GPA*lcost+GPA*rank+
llibvol*lcost+llibvol*rank+lcost*rank, data = LAWSCH85)
print(prueba_white2)
##
## studentized Breusch-Pagan test
##
## data: modelo_salario
## BP = 34.295, df = 20, p-value = 0.0242
Como pvalue=0.0242<0.05 se rechaza la H0, por lo tanto hay evidencia que la varianza de los residuos es heterocedástica.
Autocorrelación de 1er orden
Prueba de Durbin-Watson
Usando librería “lmtest”
library(lmtest)
dwtest(modelo_salario, alternative = "two.sided", iterations = 1000)
##
## Durbin-Watson test
##
## data: modelo_salario
## DW = 1.7058, p-value = 0.07519
## alternative hypothesis: true autocorrelation is not 0
Usando libreria “car”
library(car)
durbinWatsonTest(modelo_salario, simulate = TRUE, reps = 1000)
## lag Autocorrelation D-W Statistic p-value
## 1 0.144458 1.705846 0.092
## Alternative hypothesis: rho != 0
En ambos casos, no hay evidencia de Autocorrelación de 1er orden en los residuos del modelo, ya que el Pvalue>0.05, por lo tanto no se rechaza la H0.
Autocorrelación de orden superior (orden “m”)
Usando libreria “lmtest”
library(lmtest)
bgtest(modelo_salario, order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_salario
## LM test = 3.2116, df = 2, p-value = 0.2007
Como el pvalue>0.05 No se rechaza la H0 por lo tanto los residuos del modelo no siguen autocorrelación de orden “2”
2. En caso de encontrar evidencia de violación de los supuestos, planteados en el literal anterior, corrija a través de un estimador HAC apropiado, el modelo propuesto. Existen evidencia de violación en el supuesto de heterocedasticidad, por lo tanto se procede a corregirlo a través de un estimador HAC
Modelo Sin corregir
options(scipen = 999999)
library(lmtest)
coeftest(modelo_salario)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.34322596 0.53251920 15.6675 < 0.00000000000000022 ***
## LSAT 0.00469647 0.00401049 1.1710 0.243722
## GPA 0.24752388 0.09003704 2.7491 0.006826 **
## llibvol 0.09499321 0.03325435 2.8566 0.004988 **
## lcost 0.03755380 0.03210608 1.1697 0.244270
## rank -0.00332459 0.00034846 -9.5408 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estimador Robusta (uso del estimador HAC)
options(scipen = 99999)
library(lmtest)
library(sandwich)
estimacion_omega1<-vcovHC(modelo_salario,type = "HC0")
coeftest(modelo_salario,vcov. = estimacion_omega1)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.34322596 0.50982819 16.3648 < 0.00000000000000022 ***
## LSAT 0.00469647 0.00447644 1.0492 0.2960540
## GPA 0.24752388 0.08861505 2.7932 0.0060073 **
## llibvol 0.09499321 0.02703852 3.5133 0.0006095 ***
## lcost 0.03755380 0.03258921 1.1523 0.2512966
## rank -0.00332459 0.00030126 -11.0356 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Comparación
options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_salario_robust<-lmrob(lsalary~LSAT+GPA+llibvol+lcost+rank, data= LAWSCH85)
stargazer(modelo_salario, modelo_salario_robust, type = "text", title = "Comparativa")
##
## Comparativa
## =================================================================
## Dependent variable:
## ----------------------------------
## lsalary
## OLS MM-type
## linear
## (1) (2)
## -----------------------------------------------------------------
## LSAT 0.005 0.003
## (0.004) (0.005)
##
## GPA 0.248*** 0.283***
## (0.090) (0.105)
##
## llibvol 0.095*** 0.089***
## (0.033) (0.029)
##
## lcost 0.038 0.045
## (0.032) (0.042)
##
## rank -0.003*** -0.003***
## (0.0003) (0.0003)
##
## Constant 8.343*** 8.464***
## (0.533) (0.529)
##
## -----------------------------------------------------------------
## Observations 136 136
## R2 0.842 0.845
## Adjusted R2 0.836 0.839
## Residual Std. Error (df = 130) 0.112 0.090
## F Statistic 138.230*** (df = 5; 130)
## =================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
Comentario: Al elaborar la corrección mediante el estimador HAC, ahora podemos tener mejores resultados al elaborar pronósticos.