#### ejercicio de series ######################### set.seed(7) sigmau<-0.7 ### es el desvío. ues<-rnorm(30,0,sigmau) ro<-0.8 sigma2error<-sigmau^2/(1-ro^2) # varianza del error sigma2error epsilon<-rnorm(1,0,sqrt(sigma2error)) for(i in 2:30) epsilon[i]<-ro*epsilon[i-1]+ues[i] epsilon x1<-c(1.034, 1.069, 0.950, 0.820, 0.758, 1.013, 0.877,1.112, 1.015, 1.105,0.992, 0.976, 1.145, 0.881,1.278, 0.987, 1.049, 0.817, 1.083, 0.939, 1.005, 1.061, 0.916, 1.072, 1.052, 1.087, 0.945, 0.905, 0.932, 0.980) x2<-c(3.02, 3.22, 2.90, 2.93, 2.93, 3.17, 2.97, 2.93, 3.04, 2.99, 3.13, 3.02, 2.92, 2.84, 2.91, 3.09,2.96, 3.08, 2.94, 2.85, 2.95, 2.86, 3.11, 3.16, 3.04, 3.03, 3.06, 3.10, 2.83, 3.08) y<-1+2*x1+x2+epsilon ajuste<-lm(y~x1+x2,x=T) # muy mala estimación!! (sabemos la verdad) ¿qué está pasando? #quizás no debería ir la intersept por cómo fueron construidas x1 y x2, la columna de 1s puede ser dependiente summary(ajuste) ajuste$x ajuste$x[,2]+ajuste$x[,3] plot(x1,y) # es decir con más de dos variables no tiene por qué haber plot(x2,y) # estructura lineal entre la respuesta y las covariables por separado predichos<-ajuste$fitted.values residuos<-y-predichos # El modelo fue construido homocedástico, normal y respuesta fue contruida para que # dependa linealmente de x1 y x2 salvo un error plot(predichos,residuos) #no hay esctructura lineal abline(0,0) plot(x1,residuos) #idem abline(0,0) plot(x2,residuos) #idem abline(0,0) resst<-rstandard(ajuste) plot(predichos,resst) #el modelo fue construido homocedástico, no vemos crecimiento en las varianzas abline(0,0) plot(1:30,residuos) # aparece gráficamente la correlación, parecería haber algún patrón entre errores abline(0,0) plot(1:30,residuos,type="b") abline(0,0) #estadístico de DW a mano durbin<-function(residuos) { m<-length(residuos) v<-rep(0,m-1) for(i in 2:m) { v[i]<-(residuos[i]-residuos[i-1])^2 } numeradord<-sum(v) denominadord<-sum(residuos^2) numeradord/denominadord } durbin(residuos) #ver la tabla correspondiente library(car) durbinWatsonTest(ajuste, alternative=c("two.sided")) library(lmtest) dwtest(ajuste, alternative=c("two.sided")) runsTest(residuos) #rechazamos en ambos tests #Normalidad qqnorm(residuos) shapiro.test(residuos) qqnorm(resst) shapiro.test(resst) ################################## transformo los datos rhogorro=function(res){ n=length(res) e1=res[1:(n-1)] e2=res[2:n] rho=sum(e1*e2)/sum(res^2) list(rho.est=rho) } rhoestimado<-rhogorro(residuos)$rho.est rhoestimado #obviamente da lo mismo que la estimación que hace durbinWatsonTest: durbinWatsonTest(ajuste, alternative=c("two.sided"))$r y2=y[2:30] y1=y[1:(30-1)] ytransf=y2-rhoestimado*y1 x11=x1[1:(30-1)] x12=x1[2:30] x21=x2[1:(30-1)] x22=x2[2:30] x1transf=x12-rhoestimado*x11 x2transf=x22-rhoestimado*x21 ajuste2<-lm(ytransf~x1transf+x2transf, x=T) summary(ajuste2) intestimado <- summary(ajuste2)$coeff[1]/(1-rhoestimado) intestimado residuos2<-ajuste2$residuals predichos2<-ajuste2$fitted.values plot(predichos2,residuos2) #no hay esctructura lineal abline(0,0) resst2<-rstandard(ajuste2) plot(predichos2,resst2) abline(0,0) plot(1:29,residuos2) abline(0,0) plot(1:29,residuos2,type="b") abline(0,0) ## mejoró mucho visualmente !!! #estadístico de DW a mano durbin(residuos2) #ver la tabla correspondiente durbinWatsonTest(ajuste2, alternative=c("two.sided")) dwtest(ajuste2, alternative=c("two.sided")) runsTest(residuos2) #Normalidad qqnorm(residuos2) shapiro.test(residuos2) qqnorm(resst2) shapiro.test(resst2) rhoestimado<-durbinWatsonTest(ajuste2, alternative=c("two.sided"))$r rhoestimado ################################## Veamos qué pasa si repetimos este proceso una segunda vez y<-ytransf x1<-x1transf x2<-x2transf y2=y[2:29] y1=y[1:(29-1)] ytransf=y2-rhoestimado*y1 x11=x1[1:(29-1)] x12=x1[2:29] x21=x2[1:(29-1)] x22=x2[2:29] x1transf=x12-rhoestimado*x11 x2transf=x22-rhoestimado*x21 ajuste2<-lm(ytransf~x1transf+x2transf, x=T) summary(ajuste2) intestimado <- summary(ajuste2)$coeff[1]/(1-rhoestimado) intestimado residuos2<-ajuste2$residuals predichos2<-ajuste2$fitted.values plot(predichos2,residuos2) #no hay esctructura lineal abline(0,0) resst2<-rstandard(ajuste2) plot(predichos2,resst2) abline(0,0) plot(1:28,residuos2) abline(0,0) plot(1:28,residuos2,type="b") abline(0,0) ## mejoró mucho visualmente !!! #estadístico de DW a mano durbin(residuos2) #ver la tabla correspondiente durbinWatsonTest(ajuste2, alternative=c("two.sided")) dwtest(ajuste2, alternative=c("two.sided")) runsTest(residuos2) #Normalidad qqnorm(residuos2) shapiro.test(residuos2) qqnorm(resst2) shapiro.test(resst2) rhoestimado<-durbinWatsonTest(ajuste2, alternative=c("two.sided"))$r rhoestimado # paramos acá