Carga de Datos
Enuanciado: 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.
Estime el modelo de regresión lineal, correspondiente y verifique el supuesto de normalidad, usando todas las pruebas vistas en clase. EL CASO DEL JB Y DEL SW REPRESENTE SUS RESULTADOS DE FORMA GRÁFICA Y COMENTE AL RESPECTO.
library(stargazer)
modelo_estimado<- lm(V~C+P+M, data = ventas_empresa)
stargazer(modelo_estimado, title = "Modelo para estimar el nivel de Ventas", type = "text")
##
## Modelo para estimar el nivel de Ventas
## ===============================================
## Dependent variable:
## ---------------------------
## V
## -----------------------------------------------
## C 0.923***
## (0.223)
##
## P 0.950***
## (0.156)
##
## M 1.298***
## (0.431)
##
## Constant 107.444***
## (18.057)
##
## -----------------------------------------------
## Observations 24
## R2 0.980
## Adjusted R2 0.977
## Residual Std. Error 9.506 (df = 20)
## F Statistic 323.641*** (df = 3; 20)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
Prueba de Normailidad usando la prueba Jaque Bera
library(tseries)
salida_JB<-jarque.bera.test(modelo_estimado$residuals)
salida_JB
##
## Jarque Bera Test
##
## data: modelo_estimado$residuals
## X-squared = 1.4004, df = 2, p-value = 0.4965
Rechazar Hipotesis Nula H0 si: Criterio de decisión:
P-value=0.4965 α=0.05 0.4965>0.05
Dado que el p-valor es de 0.4965 mayor que 0.05, existe evidencia estadística suficiente para no rechazar la hipótesis nula y se concluye que los residuos se distribuyen normalmente. Por lo tanto, el modelo cumple con el supuesto de normalidad bajo la prueba de Jarque-Bera.
Representacion Grafica
library(fastGraph)
alpha_sig<-0.05
JB<-salida_JB$statistic
gl<-salida_JB$parameter
VC<-qchisq(1-alpha_sig,gl,lower.tail = TRUE)
shadeDist(JB,ddist = "dchisq",
parm1 = gl,
lower.tail = FALSE,xmin=0,
sub=paste("VC:",round(VC,2)," ","JB:",round(JB,2)))
Rechazar Hipotesis Nula H0 si: Criterio de decisión:
Rechazar H0 sí JB≥VC
Usando el P-value
Interpretacion P-value=0.4965 α=0.05 0.4965>0.05
Dado que el p-valor es de 0.4965 mayor que 0.05, existe evidencia estadística suficiente para no rechazar la hipótesis nula y se concluye que los residuos se distribuyen normalmente. Por lo tanto, el modelo cumple con el supuesto de normalidad bajo la prueba de Jarque-Bera.
Verificacion de Normalidad Usando la Prueba Kolmogorow Smirnov
Usando Nortest
library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals)
prueba_KS
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo_estimado$residuals
## D = 0.13659, p-value = 0.2935
Rechazar Hipotesis Nula H0 si: Criterio de decisión: Rechazar H0 sí KS≥VC Usando el P-value - Rechazar H0 sí Pvalue≤α Los residuos no se distribuyen normalmente Interpretacion P-value=0.2935 α=0.05 0.2935>0.05 Dado que el p-valor es de 0.2935 mayor que 0.05, existe evidencia estadística suficiente para no rechazar la hipótesis nula y se concluye que los residuos se distribuyen normalmente. Por lo tanto, el modelo cumple con el supuesto de normalidad bajo la prueba de Kolmogorv Sminorv.
Verificacion del supuesto de Normalidad Usando la Prueba Shapiro Wilk Usando la liberia stats (Precargada al inicar R)
salida_SW<-shapiro.test(modelo_estimado$residuals)
print(salida_SW)
##
## Shapiro-Wilk normality test
##
## data: modelo_estimado$residuals
## W = 0.95315, p-value = 0.3166
Rechazar Hipotesis Nula H0 si: Criterio de decisión: Rechazar H0 sí SW≥VC Usando el P-value - Rechazar H0 sí Pvalue≤α Los residuos no se distribuyen normalmente Interpretacion P-value=0.3166 α=0.05 0.2935>0.05 Dado que el p-valor es de 0.3116 mayor que 0.05, existe evidencia estadística suficiente para no rechazar la hipótesis nula y se concluye que los residuos se distribuyen normalmente. Por lo tanto, el modelo cumple con el supuesto de normalidad bajo la prueba de Shapiro Wilk.
Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. EN EL CASO DE LA LA PRUEBA DE FG, REPRESENTE SUS RESULTADOS GRÁFICAMENTE Y COMENTE AL RESPECTO
#2. Diagnostico de Multicolinealidad
Calculo Manual
library(stargazer)
x_mat<-model.matrix(modelo_estimado)
stargazer(head(x_mat, n=6), type = "text")
##
## =========================
## (Intercept) C P M
## -------------------------
## 1 1 197 173 110
## 2 1 208 152 107
## 3 1 181 150 99
## 4 1 194 150 102
## 5 1 192 163 109
## 6 1 196 179 114
## -------------------------
xx_matrix<-t(x_mat) %*% x_mat
stargazer(xx_matrix, type = "text")
##
## ===================================================
## (Intercept) C P M
## ---------------------------------------------------
## (Intercept) 24 5,308 4,503 2,971
## C 5,308 1,187,852 1,007,473 664,534
## P 4,503 1,007,473 859,157 564,389
## M 2,971 664,534 564,389 372,387
## ---------------------------------------------------
library(stargazer)
options(scipen = 9999)
Sn<-solve(diag(sqrt(diag(xx_matrix))))
stargazer(Sn, type = "text")
##
## =======================
## 0.204 0 0 0
## 0 0.001 0 0
## 0 0 0.001 0
## 0 0 0 0.002
## -----------------------
X^tX Normalizada: Calcular la Matriz Normalizada
library(stargazer)
xx_norm<-(Sn%*%xx_matrix)%*%Sn
stargazer(xx_norm, type = "text", digits = 4)
##
## ===========================
## 1 0.9941 0.9917 0.9938
## 0.9941 1 0.9973 0.9992
## 0.9917 0.9973 1 0.9978
## 0.9938 0.9992 0.9978 1
## ---------------------------
Autovalores de X^tX
library(stargazer)
#autovalores
lamdas<-eigen(xx_norm, symmetric = TRUE)
stargazer(lamdas$values, type = "text")
##
## =======================
## 3.987 0.010 0.003 0.001
## -----------------------
CALCULO DE K(X)=√(λmax/λmin)
k<-sqrt(max(lamdas$values)/min(lamdas$values))
print(k)
## [1] 71.16349
Sí κ(x) es inferior o igual a 20, la multicolinealidad es leve, no se considera un problema. Para 20<κ(x)<30, la multicolinealidad se considera moderada. En el caso de que κ(x)≥30 la multicolinealidad es severa.
INTERPRETACION Como k(x)>20 con un valor de 71.1635 la multicolinealidad no se considera un problema
Calculo del indice de condicion utilizando la libreria mctest
library(mctest)
x_mat<-model.matrix(modelo_estimado)
mctest(mod = modelo_estimado)
##
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf,
## theil = theil, cn = cn)
##
##
## Overall Multicollinearity Diagnostics
##
## MC Results detection
## Determinant |X'X|: 0.0346 0
## Farrar Chi-Square: 71.2080 1
## Red Indicator: 0.8711 1
## Sum of Lambda Inverse: 20.9196 1
## Theil's Method: 0.5430 1
## Condition Number: 71.1635 1
##
## 1 --> COLLINEARITY is detected by the test
## 0 --> COLLINEARITY is not detected by the test
Como se puede observar donde dice “condition Number” es el indice y podemos concluir que: * Como k(x)>30 con un valor de 71.1635 la multicolinealidad se considera un problema.
Calculo del Indice de Condicion usando la libreria “Otsrr”.
library(olsrr)
ols_eigen_cindex(model=modelo_estimado)
## Eigenvalue Condition Index intercept C P
## 1 3.9869237681 1.00000 0.0007213226 0.00009632167 0.0002714936
## 2 0.0095007154 20.48523 0.8775733849 0.00494748750 0.0877003286
## 3 0.0027882470 37.81406 0.1182991896 0.15940219484 0.8478049007
## 4 0.0007872695 71.16349 0.0034061029 0.83555399599 0.0642232771
## M
## 1 0.00008216652
## 2 0.00754342318
## 3 0.06362314121
## 4 0.92875126909
CALCULO MANUAL PARA REALIZAR LA PRUEBA DE FARRAR-GLAUBAR Caculo lRl
library(stargazer)
Zn<-scale(x_mat[, -1])
stargazer(head(Zn, n=6), type = "text")
##
## ======================
## C P M
## ----------------------
## 1 -0.983 -0.587 -0.975
## 2 -0.536 -1.430 -1.187
## 3 -1.634 -1.510 -1.753
## 4 -1.105 -1.510 -1.541
## 5 -1.186 -0.988 -1.046
## 6 -1.024 -0.346 -0.692
## ----------------------
Calcular la Matriz R
library(stargazer)
n<-nrow(Zn)
R<-(t(Zn)%*%Zn)*(1/(n-1))
stargazer(R, type = "text", digits = 4)
##
## ======================
## C P M
## ----------------------
## C 1 0.8205 0.9312
## P 0.8205 1 0.8579
## M 0.9312 0.8579 1
## ----------------------
Calcular lRl
determinante_R<-det(R)
print(determinante_R)
## [1] 0.03459107
Aplicando la prueba de Farrer Glaubar (Barlett)
m<-ncol(x_mat[, -1])
n<-nrow(x_mat[, -1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)
## [1] 71.20805
Encontrar el valor critico V.C.
gl<-m*(m-1)/2
VC<-qchisq(p=0.95, df=gl)
print(VC)
## [1] 7.814728
Regla de desicion: Como X^2FG es de 71.20805 > que el valor critico de 7.814728 se rechaza H0, por lo tanto hay evidencia de colinealidad en los regresores.
** Calculando del estadistico utilizando la libreria “mctest and pysich”
Calculo de FG usando mctest
library(mctest)
mctest::omcdiag(mod = modelo_estimado)
##
## Call:
## mctest::omcdiag(mod = modelo_estimado)
##
##
## Overall Multicollinearity Diagnostics
##
## MC Results detection
## Determinant |X'X|: 0.0346 0
## Farrar Chi-Square: 71.2080 1
## Red Indicator: 0.8711 1
## Sum of Lambda Inverse: 20.9196 1
## Theil's Method: 0.5430 1
## Condition Number: 71.1635 1
##
## 1 --> COLLINEARITY is detected by the test
## 0 --> COLLINEARITY is not detected by the test
y nuevamente tenemos que el Chi-Square es de 71.2080 superior al valor critico de 7.8 encontrado por lo tanto se detecta como bien lo dice el diagnostico colinealidad en el modelo.
Calculo de FG usando la pysych
library(psych)
FG_test<-cortest.bartlett(x_mat[, -1])
print(FG_test)
## $chisq
## [1] 71.20805
##
## $p.value
## [1] 0.000000000000002352605
##
## $df
## [1] 3
Resultado de forma grafica
library(fastGraph)
# Usamos tus variables calculadas: chi_FG (71.2080) y gl (3)
shadeDist(xshade = chi_FG,
ddist = "dchisq",
parm1 = gl,
lower.tail = FALSE,
main = "Prueba de Farrar-Glauber (Bartlett)",
sub = paste("Chi-cuadrado =", round(chi_FG, 4)))
Nuevamente se confirma que es de 31.38122 el estadistico FG mayor que el valor critico y que hay suficiente videncia para rechazar la hipotesis nula concluyendo que hay colinealidad en los regresores.
Usando la libreria “lmtest” para verificar si la varianza residual es homocedastica a traves de la prueba de White (incluya los terminos cruzados).
library(stargazer)
resid<-modelo_estimado$residuals
data_prueba_white<-as.data.frame(cbind(resid, modelo_estimado))
regresion_axuliar<-lm(formula=I(resid^2)~C+P+M+I(C^2)+I(P^2)+I(M^2)+
C*P + C*M + P*M, data = ventas_empresa)
sumario<-summary(regresion_axuliar)
n<-nrow(data_prueba_white)
R_2<-sumario$r.squared
LM_w<-n*R_2
gl=3*2+choose(3,2)
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 Critico","p value")
stargazer(Salida_white, title = "Resultados de la prueba White", type = "text",digits = 2)
##
## Resultados de la prueba White
## ==========================
## LMw Valor Critico p value
## --------------------------
## 7.12 16.92 0.62
## --------------------------
Rechazar Hipotesis Nula H0 si: Criterio de decisión:
Rechazar H0 sí LMW≥VC
Rechazar H0 sí Pvalue≤α 7.12<16.92 En este caso utilizare el estadistico LMw. El resultado da que el valor critico de 16.92 es mayor que LMw de 7.12, rechazando la hipotesis Nula (H0). Por lo tanto se concluye que hay evidencia de que la varianza de los residuos es homocedastica.
Calculamos el estadisctico para realizar la prueba White usando la libreria lmtest
library(lmtest)
Wtest<-bptest(modelo_estimado, ~I(C^2)+I(P^2)+I(M^2)+
C*P + C*M + P*M, data = ventas_empresa)
print(Wtest)
##
## studentized Breusch-Pagan test
##
## data: modelo_estimado
## BP = 7.1227, df = 9, p-value = 0.6244
Rechazar Hipotesis Nula H0 si: Criterio de decisión:
Rechazar H0 sí BP≥VC
Rechazar H0 sí Pvalue≤α
Utlizando el P_value para la prueba de hipotesis se tiene que 0.6244>0.05 Por lo tanto se concluye que: hay evidencia de que la varianza de los residuos es homocedastica.
Resultados de forma gráfica a través de la librería fastGraph
gl=3*2+choose(3,2)
VC<-qchisq(p = 0.95,df = gl)
library(fastGraph)
shadeDist(Wtest$statistic,ddist = 'dchisq',col=c('blue','green'),parm1 =VC,lower.tail = FALSE)
Verifique el supuesto de No Autocorrelación, (verifique primero
y segundo orden). USE LA PRUEBA DEL MULTIPLICADOR DE LAGRANGE,
REPRESENTE GRÁFICAMENTE SUS RESULTADOS Y COMENTE AL
RESPECTO.
Utilizamos la libreria lmtest
library(lmtest)
bgtest(modelo_estimado, order = 2)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_estimado
## LM test = 3.8409, df = 2, p-value = 0.1465
Rechazar Hipotesis Nula H0 si: Criterio de decisión:
Utlizando el P_value para la prueba de hipotesis se tiene que 0.1465>0.05
Ahora usamos BG para berificar la autocorrelacion de 1° orden
library(lmtest)
bgtest(modelo_estimado, order=1)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_estimado
## LM test = 2.5963, df = 1, p-value = 0.1071
Prueba de hipotesis: Regla de desicion
Rechazar H0 si pvalue≤α
No rechazar H0 si pvalue>=α
Como pvalue>0.05 No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden “1”