library(AER)
## Warning: package 'AER' was built under R version 3.6.3
## Loading required package: car
## Warning: package 'car' was built under R version 3.6.3
## Loading required package: carData
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.6.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.6.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 3.6.3
## Loading required package: survival
library(prais)
## Warning: package 'prais' was built under R version 3.6.3
library(orcutt)
## Warning: package 'orcutt' was built under R version 3.6.3
library(sandwich)
data("USConsump1993",package ="AER")
d=as.data.frame(USConsump1993)
plot(d$income, type="l")

head(d)
## income expenditure
## 1 6284 5820
## 2 6390 5843
## 3 6476 5917
## 4 6640 6054
## 5 6628 6099
## 6 6879 6365
#3. En Macroeconomia existe la hipotesis del ingreso permanente que relaciona al consumo con el ingreso. Utiliza la variable de gasto como variable proxy del consumo y la variable ingreso para correr el siguiente modelo:
reg=lm(expenditure~income,data=USConsump1993)
summary(reg)
##
## Call:
## lm(formula = expenditure ~ income, data = USConsump1993)
##
## Residuals:
## Min 1Q Median 3Q Max
## -294.52 -67.02 4.64 90.02 325.84
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -65.795821 90.990824 -0.723 0.474
## income 0.915623 0.008648 105.874 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 153.6 on 42 degrees of freedom
## Multiple R-squared: 0.9963, Adjusted R-squared: 0.9962
## F-statistic: 1.121e+04 on 1 and 42 DF, p-value: < 2.2e-16
#Por la teoría del ingreso permanente por los resultados que arroja, morimos con dinero por eso es negativo. Con gasto negativo
#4. Grafica los residuales y explica que observas. ¿Son estos errores independientes?
plot(reg$residuals,type="l")

#5. Estima la prueba de Durbin-Watson a mano e indica si la regresión de consumo e ingreso presenta autocorrelación de primer orden. Verifica tu respuesta con el comando dwtest de R.
dwtest(reg)[1]
## $statistic
## DW
## 0.4607776
dwtest(reg)
##
## Durbin-Watson test
##
## data: reg
## DW = 0.46078, p-value = 3.274e-11
## alternative hypothesis: true autocorrelation is greater than 0
#de 0 a 2 autocorrelación positiva, de 2 a 4 negativa y 2 no hay.
#6. Estima la prueba del Multiplicador de Lagrage para autocorrelacion de orden 1 a mano y verifica tu respuesta con el comando de R para esta prueba.
bgtest(reg)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: reg
## LM test = 24.901, df = 1, p-value = 6.034e-07
#7. Utilizando la prueba de DW, encuentra el posible valor del parémetro rho del error. Estima cuasiprimeras diferencias y grafica los residuales. Explica que encuentras. Has de nuevo la prueba DW.
rho=1-(as.numeric(dwtest(reg)[1])/2)
rho
## [1] 0.7696112
n=44
cuasiincome=d$income[2:n]-(rho*d$income[1:n-1])
plot(cuasiincome)

plot(d$income)

cuasiexpenditure=d$expenditure[2:n]-(rho*d$expenditure[1:n-1])
plot(cuasiexpenditure)

plot(d$expenditure)

reg2=lm(cuasiexpenditure~cuasiincome)
summary(reg2)
##
## Call:
## lm(formula = cuasiexpenditure ~ cuasiincome)
##
## Residuals:
## Min 1Q Median 3Q Max
## -239.766 -56.881 5.583 61.416 235.762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -39.55245 62.76270 -0.63 0.532
## cuasiincome 0.92645 0.02425 38.20 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.8 on 41 degrees of freedom
## Multiple R-squared: 0.9727, Adjusted R-squared: 0.972
## F-statistic: 1459 on 1 and 41 DF, p-value: < 2.2e-16
plot(reg2$residuals,type="l")

dwtest(reg2)
##
## Durbin-Watson test
##
## data: reg2
## DW = 1.9314, p-value = 0.3513
## alternative hypothesis: true autocorrelation is greater than 0
#Mejoramos mucho la autocorrelación pero no la eliminamos por completo.
#8. Utiliza la funcion del estimador de Cochrane y Orcutt y encuentra el valor de rho. ¿Cuál es el tiempo de convergencia?. Estima cuasiprimeras diferencias y grafica los residuales. Has de nuevo la prueba DW. Explica que encuentras.
tab=cochrane.orcutt(reg)
rho1=tab$rho
rho1
## [1] 0.7828935
cuasiincome1=d$income[2:n]-(rho1*d$income[1:n-1])
cuasiexpenditure1=d$expenditure[2:n]-(rho1*d$expenditure[1:n-1])
reg3=lm(cuasiexpenditure1~cuasiincome1)
summary(reg3)
##
## Call:
## lm(formula = cuasiexpenditure1 ~ cuasiincome1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -239.196 -57.904 4.692 59.147 234.966
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.9771 62.7304 -0.589 0.559
## cuasiincome1 0.9264 0.0256 36.181 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.8 on 41 degrees of freedom
## Multiple R-squared: 0.9696, Adjusted R-squared: 0.9689
## F-statistic: 1309 on 1 and 41 DF, p-value: < 2.2e-16
plot(reg3$residuals,type="l")

dwtest(reg3)
##
## Durbin-Watson test
##
## data: reg3
## DW = 1.9567, p-value = 0.3832
## alternative hypothesis: true autocorrelation is greater than 0
#el estimador de orcutt nos acerca más al dos.
#9. Investiga cual es la metodología de Prais-Winsten (1954). Utiliza la función del estimador de Prais-Winsten de R y encuentra el valor de rho. ¿Cuál es el tiempo de convergencia? Estima cuasiprimeras diferencias y grafica los residuales. Has de nuevo la prueba DW. Explica que encuentras.
prais_winsten(reg,data=d)
## Iteration 0: rho = 0
## Iteration 1: rho = 0.7921
## Iteration 2: rho = 0.8021
## Iteration 3: rho = 0.8037
## Iteration 4: rho = 0.804
## Iteration 5: rho = 0.804
## Iteration 6: rho = 0.804
## Iteration 7: rho = 0.804
## Iteration 8: rho = 0.804
##
## Call:
## prais_winsten(formula = reg, data = d)
##
## Coefficients:
## (Intercept) income
## -15.1436 0.9142
##
## AR(1) coefficient rho: 0.804
#Nos da otro intercepto totalmente diferente.
rho3=0.804
rho3
## [1] 0.804
cuasiincome2=d$income[2:n]-(rho3*d$income[1:n-1])
cuasiexpenditure2=d$expenditure[2:n]-(rho3*d$expenditure[1:n-1])
reg4=lm(cuasiexpenditure2~cuasiincome2)
summary(reg4)
##
## Call:
## lm(formula = cuasiexpenditure2 ~ cuasiincome2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -238.152 -59.415 3.069 57.763 233.892
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -32.2848 62.6884 -0.515 0.609
## cuasiincome2 0.9260 0.0281 32.951 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.8 on 41 degrees of freedom
## Multiple R-squared: 0.9636, Adjusted R-squared: 0.9627
## F-statistic: 1086 on 1 and 41 DF, p-value: < 2.2e-16
plot(reg4$residuals,type="l")

dwtest(reg4)
##
## Durbin-Watson test
##
## data: reg4
## DW = 1.9933, p-value = 0.4311
## alternative hypothesis: true autocorrelation is greater than 0
#10. Estima los coeficientes de primeras diferencias y grafica los residuales. Has de nuevo la prueba DW. Explica que encuentras.
rho4=1
cuasiincome3=d$income[2:n]-(rho4*d$income[1:n-1])
cuasiexpenditure3=d$expenditure[2:n]-(rho4*d$expenditure[1:n-1])
reg5=lm(cuasiexpenditure3~cuasiincome3)
summary(reg5)
##
## Call:
## lm(formula = cuasiexpenditure3 ~ cuasiincome3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.16 -57.23 -20.96 44.22 209.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.21328 22.78489 1.765 0.085 .
## cuasiincome3 0.72506 0.09026 8.033 5.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 100.1 on 41 degrees of freedom
## Multiple R-squared: 0.6115, Adjusted R-squared: 0.602
## F-statistic: 64.53 on 1 and 41 DF, p-value: 5.971e-10
plot(reg5$residuals,type="l")

dwtest(reg5)
##
## Durbin-Watson test
##
## data: reg5
## DW = 1.8041, p-value = 0.2548
## alternative hypothesis: true autocorrelation is greater than 0
#Primeras diferencias no es lo mejor la prueba de dw nos dice que no desaparece la autocorrelación pero p value disminuye.
#11. Así como existe la matriz de var-cov de estimador de White para heteroscedasticidad, tambien existe la matriz var-cov de Newey-West que es consistente para autocorrelación y heteroscedasticidad. Utiliza el comando NeweyWest de R para estimar los errores estandar de tu regresión consistentes con autocorrelación. Estima esta matriz con 3 lags dentro de las opciones de la función.
coeftest(reg,vcov=NeweyWest(reg,lag=3,prewhite = FALSE,adjust=TRUE))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -65.795821 133.345400 -0.4934 0.6243
## income 0.915623 0.015458 59.2319 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#el intercepto regresó al negativo. Corrige la matriz y luego la usa. Nos ayuda con las pruebas t pero no se corrigió.