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ó.