ejercicio de libro de Damodar Gujarati,pag.335

“estimacion de la Función de consumo para Estados Unidos,1947-2000”

Estimando el modelo

options(scipen = 9999)
library(stargazer)
library(readxl)
funcion_de_Consumo2_EEUU <- read_excel("Practicas_ECMA/datos_Excel/Datos_Gujarati/funcion_de_Consumo2_EEUU.xlsx")
modelo_EEUU<-lm(formula =lC~lYd+lW+I,data = funcion_de_Consumo2_EEUU)
stargazer(modelo_EEUU,title = "Función de consumo para Estados Unidos,1947-2000",type = "html",digits = 8)
Función de consumo para Estados Unidos,1947-2000
Dependent variable:
lC
lYd 0.80487290***
(0.01749786)
lW 0.20127000***
(0.01759260)
I -0.00268906***
(0.00076193)
Constant -0.46771110***
(0.04277803)
Observations 54
R2 0.99955970
Adjusted R2 0.99953320
Residual Std. Error 0.01193376 (df = 50)
F Statistic 37,832.61000000*** (df = 3; 50)
Note: p<0.1; p<0.05; p<0.01 
load("C:/Users/manue/Desktop/Practicas_ECMA/datos_Rdata/Datos_Gujarati/Función_de_consumo_EEUU.RData")
print(Media_de_la_variable_Y)
     [,1]

[1,] 7.826093

print(Desviacion_estandar_de_Y)
      [,1]

[1,] 0.5523683

Los resultados demuestran que todos los coeficientes estimados son muy significativos desde el punto de vista estadístico, pues sus valores p son muy pequeños. Los coeficientes estimados se interpretan como sigue: la elasticidad del ingreso es≈0.80, lo que indica que, cuando las demás variables se mantienen constantes, si el ingreso aumenta 1%, la media del gasto de consumo aumenta alrededor de 0.8%. El coeficiente de riqueza es≈0.20, lo que significa que si la riqueza aumenta 1%, la media del consumo se incrementa sólo 0.2%, de nuevo cuando las demás variables se mantienen constantes. El coeficiente de la variable tasa de interés indica que,a medida que la tasa de interés aumenta un punto porcentual, el gasto de consumo disminuye0.26%, ceteris paribus.

Pruebas de Normalidad MRLM

1.Prueba de Normalidad de Jarque - Bera

library(normtest)
jb.norm.test(modelo_EEUU$residuals)
## 
##  Jarque-Bera test for normality
## 
## data:  modelo_EEUU$residuals
## JB = 1.8055, p-value = 0.2655

1.8055<5.9915 Como(JB<VC) NO rechaza Ho es decir los residuos tienen una distribucion normal.

0.2625>0.05 Como(P>alfa) NO rechaza Ho es decir los residuos tienen una distribucion normal.

2.Prueba de Normalidad de Kolmogorov - Smirnov

library(nortest)
lillie.test((modelo_EEUU$residuals))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  (modelo_EEUU$residuals)
## D = 0.085441, p-value = 0.4197

0.085441<0.1191(KS<VC) NO rechaza Ho es decir los residuos tienen una distribucion normal.

0.4197>0.05(P>alfa) NO rechaza Ho es decir los residuos tienen una distribucion normal.

3.Prueba de Normalidad de Shapiro - Wilk

shapiro.test(modelo_EEUU$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_EEUU$residuals
## W = 0.96726, p-value = 0.1457

Calcular Wn

W<-0.96726
mew<-0.0038915*I(log(54))^3-0.083751*I(log(54))^2-0.31082*I(log(54))-1.5861
des<-exp(0.0030302*I(log(54))^2-0.083751*(log(54))-0.4803)
Wn<-(log(1-W)-mew)/des
print(Wn)
## [1] 1.059484

1.059484<1.644854 Como(Wn<VC) NO rechaza Ho es decir los residuos tienen una distribucion normal.

0.1457>0.05 Como(P>alfa) NO rechaza Ho es decir los residuos tienen una distribucion normal.

Pruebas de Multicolinealidad.

Indice de condicion( calculo manual)

Calcular la matriz |X X|

# matriz Xmat
Xmat<-model.matrix(modelo_EEUU)
print(Xmat)
##    (Intercept)      lYd        lW          I
## 1            1 6.942350  8.550012 -10.350940
## 2            1 6.993933  8.571825  -4.719804
## 3            1 6.999057  8.631834   1.044063
## 4            1 7.083975  8.658609   0.407346
## 5            1 7.112327  8.713756  -5.283152
## 6            1 7.144249  8.739355  -0.277011
## 7            1 7.191053  8.757094   0.561137
## 8            1 7.203406  8.824241  -0.138476
## 9            1 7.268084  8.877974   0.261997
## 10           1 7.314753  8.905876  -0.736124
## 11           1 7.339213  8.897721  -0.260683
## 12           1 7.348394  8.970810  -0.574630
## 13           1 7.392524  9.010432   2.295943
## 14           1 7.417460  9.030227   1.511181
## 15           1 7.450080  9.101850   1.296432
## 16           1 7.497485  9.115100   1.395922
## 17           1 7.534496  9.152298   2.057616
## 18           1 7.604347  9.210680   2.026599
## 19           1 7.664347  9.265095   2.111669
## 20           1 7.716283  9.261227   2.020251
## 21           1 7.758120  9.333626   1.212616
## 22           1 7.803108  9.404707   1.054986
## 23           1 7.833719  9.364970   1.732154
## 24           1 7.874739  9.363065   1.166228
## 25           1 7.917646  9.418404  -0.712241
## 26           1 7.963564  9.510439  -0.155737
## 27           1 8.030182  9.478913   1.413839
## 28           1 8.023520  9.381668  -1.042571
## 29           1 8.041896  9.444175  -3.533585
## 30           1 8.084408  9.507238  -0.656766
## 31           1 8.119905  9.531431  -1.190427
## 32           1 8.168345  9.578484   0.113048
## 33           1 8.196602  9.638219   1.704210
## 34           1 8.204672  9.678151   2.298496
## 35           1 8.227135  9.678153   4.703847
## 36           1 8.240570  9.699688   4.449027
## 37           1 8.270499  9.737719   4.690972
## 38           1 8.344648  9.771484   5.848332
## 39           1 8.377425  9.855785   4.330504
## 40           1 8.408850  9.929644   3.768031
## 41           1 8.430000  9.963439   2.819469
## 42           1 8.473053 10.013775   3.287061
## 43           1 8.498316 10.071533   4.317956
## 44           1 8.520029 10.047810   3.595025
## 45           1 8.523772 10.087899   1.802757
## 46           1 8.554354 10.103084   1.007439
## 47           1 8.568133 10.130318   0.624790
## 48           1 8.593636 10.135337   2.206002
## 49           1 8.619587 10.219747   3.333143
## 50           1 8.644302 10.290388   3.083201
## 51           1 8.674966 10.394031   3.120000
## 52           1 8.727227 10.479736   3.583909
## 53           1 8.751474 10.586364   3.245271
## 54           1 8.785570 10.549745   3.575970
## attr(,"assign")
## [1] 0 1 2 3
# matriz XXmat
XXmat<-t(Xmat)%*%Xmat
print(XXmat)
##             (Intercept)       lYd        lW         I
## (Intercept)    54.00000  428.4718  512.6252  65.44629
## lYd           428.47179 3416.0687 4083.5998 568.75115
## lW            512.62518 4083.5998 4882.7449 671.47742
## I              65.44629  568.7511  671.4774 478.60612
## sn calculo de la matriz de normalizacion
sn<-diag(1/sqrt(diag(XXmat)))
print(sn)
##           [,1]       [,2]       [,3]       [,4]
## [1,] 0.1360828 0.00000000 0.00000000 0.00000000
## [2,] 0.0000000 0.01710948 0.00000000 0.00000000
## [3,] 0.0000000 0.00000000 0.01431093 0.00000000
## [4,] 0.0000000 0.00000000 0.00000000 0.04570996
## XXmat Normalizada 
XX_norm<-(sn%*%XXmat)%*%sn
print(XX_norm)
##           [,1]      [,2]      [,3]      [,4]
## [1,] 1.0000000 0.9976129 0.9983229 0.4070981
## [2,] 0.9976129 1.0000000 0.9998799 0.4448052
## [3,] 0.9983229 0.9998799 1.0000000 0.4392485
## [4,] 0.4070981 0.4448052 0.4392485 1.0000000
## autovalores
library(stargazer)
#autovalores
lambdas<-eigen(XX_norm,symmetric = TRUE)$values
print(lambdas)
## [1] 3.24479087364 0.75342466252 0.00172912759 0.00005533624
# Indice de condicion
K<-sqrt(max(lambdas)/min(lambdas))
print(K)
## [1] 242.1523

En el caso de que κ(x)≥30 la multicolinealidad es severa. Y dado que el valor obtenido en en el indice de condicion es de 242.1523, por lo tanto existe una alta multicolinealidad entre los regresores del modelo.

Uso de la libreria “mctest”, para obtener el indice de condición

library(mctest)
source("C:/Users/manue/Desktop/Practicas_ECMA/datos_Rdata/correccion_eigprop.R")
my_eigprop(mod=modelo_EEUU)
## 
## Call:
## my_eigprop(mod = modelo_EEUU)
## 
##   Eigenvalues       CI (Intercept)    lYd     lW      I
## 1      3.2448   1.0000      0.0001 0.0000 0.0000 0.0157
## 2      0.7534   2.0753      0.0001 0.0000 0.0000 0.6119
## 3      0.0017  43.3191      0.5335 0.0216 0.0046 0.3658
## 4      0.0001 242.1523      0.4663 0.9784 0.9954 0.0067
## 
## ===============================
## Row 4==> lYd, proportion 0.978426 >= 0.50 
## Row 4==> lW, proportion 0.995357 >= 0.50 
## Row 2==> I, proportion 0.611885 >= 0.50

Prueba de Farrar-Glaubar ( manual)

library(fastGraph)
m<-ncol(Xmat[,-1]) # cantidad de variables explicativas k-1
n<-nrow(Xmat)
determinante_R<- det(cor(Xmat[,-1])) # determinanre de la matriz de correlacion
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)
## [1] 206.8664
## Valor Critico
gl<-m*(m-1)/2
VC<-qchisq(0.05,gl,lower.tail = FALSE)
print(VC)
## [1] 7.814728
shadeDist(xshade = chi_FG,ddist = "dchisq",parm1 = gl,lower.tail = FALSE,sub=paste("VC:",VC,"FG:",chi_FG))

Como χ2FG≥V.C. se rechaza H0, por lo tanto hay evidencia de colinealidad en los regresores.

Prueba de Farrar-Glaubar con mctest

library(mctest)
mctest(modelo_EEUU)
## 
## 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.0175         0
## Farrar Chi-Square:       206.8664         1
## Red Indicator:             0.7601         1
## Sum of Lambda Inverse:    72.2108         1
## Theil's Method:            0.3298         0
## Condition Number:        918.0045         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

Prueba de Farrar-Glaubar Uso de la libreria psych

library(psych)
library(fastGraph)
FG_test<-cortest.bartlett(Xmat[,-1])
VC_1<-qchisq(0.05,FG_test$df,lower.tail = FALSE)
print(FG_test)
## $chisq
## [1] 206.8664
## 
## $p.value
## [1] 0.00000000000000000000000000000000000000000001384868
## 
## $df
## [1] 3
shadeDist(xshade = FG_test$chisq,ddist = "dchisq",parm1 = FG_test$df,lower.tail = FALSE,sub=paste("VC:",VC_1,"FG:",FG_test$chisq))

Cálculando los VIF para el modelo estimado(forma manual)

Matriz de Correlación de los regresores del modelo (Como se obtuvo con anterioridad):

VIF<-diag(solve(cor(Xmat[,-1])))
print(VIF)
##      lYd       lW        I 
## 35.02073 35.56237  1.62765

Obtención de los VIF’s, a través de la librería “car”

library(car)
VIFs_car<-vif(modelo_EEUU)
print(VIFs_car)
##      lYd       lW        I 
## 35.02073 35.56237  1.62765

Obtención de los VIF’s, a través de la librería “mctest”

library(mctest)
mc.plot(modelo_EEUU,vif = 2,)

existe una alta colinealidad entre las variables de ingreso personal disponible y el nivel de riqueza de los individuos.

Pruebas de Heterocedásticidad y Autocorrelación

prueba de White

library(stargazer)
u_i<-modelo_EEUU$residuals
data_prueba_white<-as.data.frame(cbind(u_i,funcion_de_Consumo2_EEUU))
regresion_auxiliar<-lm(I(u_i^2)~lYd+lW+I+I(lYd^2)+I(lW^2)+I(I^2)+lYd*lW+lYd*I+lW*I,data = data_prueba_white)
sumario<-summary(regresion_auxiliar)
n<-nrow(data_prueba_white)
R_2<-sumario$r.squared
LM_w<-n*R_2
gl=3+3+3
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
13.391770 16.918980 0.145665

Como 0.145665>0.05 No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica.

De igual forma como LMw<VC No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica.

Uso de la librería “lmtest”

library(lmtest)
prueba_white<-bptest(modelo_EEUU,~I(lYd^2)+I(lW^2)+I(I^2)+lYd*lW+lYd*I+lW*I,data = funcion_de_Consumo2_EEUU)
print(prueba_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_EEUU
## BP = 13.392, df = 9, p-value = 0.1457

Como 0.145665>0.05 No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica.

De igual forma como LMw<VC No se rechaza la H0, por lo tanto hay evidencia de que la varianza de los residuos es homocedástica.

Pruebas de Autocorrelación

Prueba de Durbin – Watson.

Usando libreria “lmtest”

library(lmtest)
dwtest(modelo_EEUU,alternative = "two.sided",iterations = 10000)
## 
##  Durbin-Watson test
## 
## data:  modelo_EEUU
## DW = 1.2892, p-value = 0.001889
## alternative hypothesis: true autocorrelation is not 0

Usando libreria “car”

library(car)
durbinWatsonTest(modelo_EEUU,simulate = TRUE,reps = 10000)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2977614      1.289219   0.003
##  Alternative hypothesis: rho != 0

En ambos casos, se puede rechazar la hipotesis nula por lo tanto Hay evidencia de autocorrelación de primer orden, en los residuos del modelo ya que el pvalue<0.05

Autocorrelación de orden superior

1.Preparación de datos:

library(dplyr)
library(tidyr)
library(kableExtra)
cbind(u_i,funcion_de_Consumo2_EEUU) %>% 
  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 Año lC lYd lW I Lag_1 Lag_2
0.0151790 1947 6.883872 6.942350 8.550012 -10.350940 0.0000000 0.0000000
0.0063945 1948 6.905854 6.993933 8.571825 -4.719804 0.0151790 0.0000000
0.0325784 1949 6.932741 6.999057 8.631834 1.044063 0.0063945 0.0151790
0.0191472 1950 6.994758 7.083975 8.658608 0.407346 0.0325784 0.0063945
-0.0153335 1951 7.009499 7.112327 8.713755 -5.283152 0.0191472 0.0325784
-0.0013297 1952 7.040887 7.144249 8.739354 -0.277011 -0.0153335 0.0191472

2.Calculando la regresión auxiliar y el estadístico LMBP

regresion_auxiliar_BG<-lm(u_i~lYd+lW+I+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
6.447585 5.991465 0.039804

Como pvalue<0.05 SI se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, siguen autocorrelación de orden “2”.

Usando la librería “lmtest”

library(lmtest)
bgtest(modelo_EEUU,order = 2)
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  modelo_EEUU
## LM test = 6.4476, df = 2, p-value = 0.0398

Como pvalue<0.05 SI se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, 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_EEUU,order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_EEUU
## LM test = 5.3122, df = 1, p-value = 0.02118

Como pvalue<0.05 se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, siguen autocorrelación de 1° orden

Estimadores HAC para la correcion de Heterocedasticidad y Autocorrelación.

Estimador HAC de NeweyWest.

en este modelo estamos ante la presencia de autocorrelacion de 1 y 2 orden por lo cual utilizaremos el estimador de NeweyWest para corregirlo ya que no hay evidencia de heterocedasticidad.

Sin corregir

options(scipen = 99999)
library(lmtest)
#Sin corregir:
coeftest(modelo_EEUU)
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value              Pr(>|t|)    
## (Intercept) -0.46771107  0.04277803 -10.9334  0.000000000000007331 ***
## lYd          0.80487293  0.01749786  45.9984 < 0.00000000000000022 ***
## lW           0.20126996  0.01759260  11.4406  0.000000000000001434 ***
## I           -0.00268906  0.00076193  -3.5293             0.0009047 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Corregido (usando un estimador HAC de NeweyWest)

library(lmtest)
library(sandwich)

#Corregido:

estimacion_omega<-NeweyWest(modelo_EEUU,lag = 2)
coeftest(modelo_EEUU,vcov. = estimacion_omega)
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value              Pr(>|t|)    
## (Intercept) -0.46771107  0.05134092 -9.1099 0.0000000000033614123 ***
## lYd          0.80487293  0.01815318 44.3378 < 0.00000000000000022 ***
## lW           0.20126996  0.01694371 11.8787 0.0000000000000003595 ***
## I           -0.00268906  0.00080794 -3.3283              0.001645 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimación Robusta

options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_EEUU_robust<-lmrob(lC~lYd+lW+I,data=funcion_de_Consumo2_EEUU)
# print(summary(modelo_EEUU_robust))
stargazer(modelo_EEUU,modelo_EEUU_robust,type = "html",title = "comparativa")
comparativa
Dependent variable:
lC
OLS MM-type
linear
(1) (2)
lYd 0.805*** 0.811***
(0.017) (0.017)
lW 0.201*** 0.196***
(0.018) (0.017)
I -0.003*** -0.003***
(0.001) (0.001)
Constant -0.468*** -0.469***
(0.043) (0.048)
Observations 54 54
R2 1.000 1.000
Adjusted R2 1.000 1.000
Residual Std. Error (df = 50) 0.012 0.013
F Statistic 37,832.610*** (df = 3; 50)
Note: p<0.1; p<0.05; p<0.01