library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(lmtest)
library(orcutt)
library(readxl)
library(prais)## Loading required package: sandwich## Loading required package: pcse##
## Attaching package: 'pcse'## The following object is masked from 'package:sandwich':
##
## vcovPC# Lectura de la Base de Datos
df_econometria <- read_excel("~/Desktop/DataScience/bases de datos/ bd_ahorrovoluntario.xlsx")# Convertir en Series de Tiempo
datos <- df_econometria[1:154, ]
datos.ts = ts(datos, start=c(2010, 1), frequency=12) # Adjusted start to include month
# Creating individual time series
fecha.ts = ts(datos[,1], start=c(2010, 1), frequency=12)
av.ts = ts(datos[,2], start=c(2010, 1), frequency=12)
ren.ts = ts(datos[,3], start=c(2010, 1), frequency=12)
sal.ts = ts(datos[,4], start=c(2010, 1), frequency=12)
trab.ts = ts(datos[,5], start=c(2010, 1), frequency=12)
mtrab.ts = ts(datos[,6], start=c(2010, 1), frequency=12)
cetes.ts = ts(datos[,7], start=c(2010, 1), frequency=12)
igae.ts = ts(datos[,8], start=c(2010, 1), frequency=12)
inf.ts = ts(datos[,9], start=c(2010, 1), frequency=12)
desem.ts = ts(datos[,10], start=c(2010, 1), frequency=12)model1 = lm(av.ts ~ ren.ts + sal.ts + mtrab.ts + trab.ts + cetes.ts + igae.ts + inf.ts + desem.ts)
summary(model1)##
## Call:
## lm(formula = av.ts ~ ren.ts + sal.ts + mtrab.ts + trab.ts + cetes.ts +
## igae.ts + inf.ts + desem.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4703.1 -313.8 24.7 316.1 4964.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.532e+03 3.530e+03 1.567 0.11923
## ren.ts -3.308e+02 1.106e+02 -2.990 0.00328 **
## sal.ts -5.813e+00 6.673e+00 -0.871 0.38511
## mtrab.ts 1.266e+04 1.708e+04 0.741 0.45972
## trab.ts -1.251e-02 1.708e-02 -0.733 0.46489
## cetes.ts -8.534e+01 8.045e+01 -1.061 0.29058
## igae.ts -4.987e+01 2.588e+01 -1.927 0.05591 .
## inf.ts -2.009e+02 2.279e+02 -0.881 0.37951
## desem.ts 1.877e+02 2.425e+02 0.774 0.44009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1058 on 145 degrees of freedom
## Multiple R-squared: 0.2265, Adjusted R-squared: 0.1838
## F-statistic: 5.308 on 8 and 145 DF, p-value: 7.548e-06Hacemos dos pruebas para ver si nuestra regresión presenta autocorrelación.
dwtest(model1)##
## Durbin-Watson test
##
## data: model1
## DW = 2.8534, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0\[H_0: p < 0.05; H_1: p > 0.05\]
Un valor de 2.8534 sugiere una posible autocorrelación negativa. Por otra parte, un valor p de 1 indica que rechazamos la hipótesis nula. Es decir, que existe autocorrelación.
bgtest(model1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model1
## LM test = 28.844, df = 1, p-value = 7.844e-08\[H_0: p > 0.05; H_1: p < 0.05\] El valor p nos dice de que efectivamente, hay autocorrelación. Por otra parte, el estadístico LM mide la presencia de autocorrelación en los residuos de un modelo de regresión. Ambas pruebas nos dan suficiente evidencia de la autocorrelación en la regresión.
model2= cochrane.orcutt(model1)
summary(model2)## Call:
## lm(formula = av.ts ~ ren.ts + sal.ts + mtrab.ts + trab.ts + cetes.ts +
## igae.ts + inf.ts + desem.ts)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5994e+03 2.5312e+03 2.212 0.02853 *
## ren.ts -3.4923e+02 6.9374e+01 -5.034 1.413e-06 ***
## sal.ts -7.3990e+00 4.2383e+00 -1.746 0.08299 .
## mtrab.ts -1.2126e+04 1.5395e+04 -0.788 0.43219
## trab.ts 1.2277e-02 1.5394e-02 0.798 0.42646
## cetes.ts -7.6041e+01 5.1878e+01 -1.466 0.14489
## igae.ts -4.5611e+01 1.9083e+01 -2.390 0.01814 *
## inf.ts -1.3981e+02 1.6166e+02 -0.865 0.38854
## desem.ts 2.0163e+02 1.7727e+02 1.137 0.25724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 948.5934 on 144 degrees of freedom
## Multiple R-squared: 0.4309 , Adjusted R-squared: 0.3993
## F-statistic: 13.6 on 8 and 144 DF, p-value: < 1.348e-14
##
## Durbin-Watson statistic
## (original): 2.85343 , p-value: 1e+00
## (transformed): 1.90112 , p-value: 1.304e-01Al parecer, nuestro modelo sigue teniendo correlación, que se observa por la alta R-Cuadrada pero la nula significancia de distintos modelos, por lo que necesitaremos realizar otra corrección y observar si nuestro modelo mejora.
model3 <- prais_winsten(model1, data = datos.ts, index = fecha.ts)## Iteration 0: rho = 0
## Iteration 1: rho = -0.4317
## Iteration 2: rho = -0.4619
## Iteration 3: rho = -0.4633
## Iteration 4: rho = -0.4634
## Iteration 5: rho = -0.4634
## Iteration 6: rho = -0.4634model3##
## Call:
## prais_winsten(formula = model1, data = datos.ts, index = fecha.ts)
##
## Coefficients:
## (Intercept) ren.ts sal.ts mtrab.ts trab.ts cetes.ts
## 5.593e+03 -3.492e+02 -7.396e+00 -1.211e+04 1.226e-02 -7.639e+01
## igae.ts inf.ts desem.ts
## -4.546e+01 -1.421e+02 2.008e+02
##
## AR(1) coefficient rho: -0.4634dwtest(model3)##
## Durbin-Watson test
##
## data: model3
## DW = 2.8534, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0bgtest(model3)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model3
## LM test = 28.844, df = 1, p-value = 7.844e-08#usando los grados de libertad del modelo para comparar LM con chi2
df <- 1
alpha <- 0.05
chi_critical <- qchisq(1 - alpha, df)
chi_critical_r <- round(chi_critical, 2)
print(paste("chi2:", chi_critical_r))## [1] "chi2: 3.84"