remove(list=ls())
#setwd("/Users/michael/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt08/R")
#setwd("/Volumes/MichaelsHD/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt08/R")
setwd("/Volumes/Daten/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt08/R")
set.seed(23514)

#Daten einlesn
dat1 <- matrix(scan("./data8a.dat"),ncol=2,byrow=T)
dat2 <- matrix(scan("./data8b.dat"),ncol=2,byrow=T)
dat3 <- matrix(scan("./data8c.dat"),ncol=2,byrow=T)

#erstes Modell
x1 <- dat1[,1]
t1 <- dat1[,2]
plot(t1,x1,main="data8a.dat")

mod1 <- lm(x1 ~ t1)
#residuen plotten
plot(mod1$residuals,type="l")
#vermutlich kein quadratischer term noetig, ueberpruefen:
mod11 <- lm(x1 ~ t1 + I(t1*t1))
summary(mod1)
summary(mod11)
#-> kein quadratischer Term noetig

library(lmtest)
#autokorreliert?
dwtest(mod1,alternative="two.sided")
dwtest(mod11,alternative="two.sided")
#Teststatistik nahe 4 -> negative Autokorrelation
qqnorm(mod1$residuals)
qqline(mod1$residuals)
tmpeps <- rnorm(101)
qqnorm(tmpeps)
qqline(tmpeps)
#modell x1 ~ t1 mit negativ autokorrelierten Eps

#zusammen plotten
op <- par(mfrow=c(2,2))
plot(t1,x1,main="data8a.dat")
plot(mod1$residuals,type="l",main="Residuen bei lm(x1 ~ t1)")
qqnorm(mod1$residuals,main="QQ-Plot der Residuen")
qqline(mod1$residuals)
qqnorm(tmpeps, main="QQ-Plot bei Normalverteilung")
qqline(tmpeps)
par(op)

##zweites Modell
x2 <- dat2[,1]
t2 <- dat2[,2]
plot(t2,x2,main="data8b.dat")

#erster Ansatz
mod2 <- lm(x2 ~ t2)
#Residuen plotten
plot(mod2$residuals,type="l")
#evtl. quadratisches Modell
mod21 <- lm(x2 ~ t2 + I(t2*t2))
plot(mod21$residuals,type="l")
#sieht besser aus, evtl. beta2=0?
mod22 <- lm(x2 ~ I(t2*t2) )
summary(mod2)
summary(mod21)
summary(mod22)
#also beta2=0

qqnorm(mod2$residuals)
qqline(mod2$residuals)
qqnorm(mod21$residuals)
qqline(mod21$residuals)
qqnorm(mod22$residuals)
qqline(mod22$residuals)
#test auf autokorrelation schlaegt nicht an
dwtest(mod22,alternative="two.sided")

#Modell x2 ~ t2*t2 mit unkorrelierten eps

#zusammen plotten
op <- par(mfrow=c(2,2))
plot(t2,x2,main="data8b.dat")
plot(mod2$residuals,type="l",main="Residuen bei lm(x2 ~ t2)")
plot(mod21$residuals,type="l",main="Residuen bei lm(x2 ~ t2 + I(t2*t2))")
qqnorm(mod22$residuals,main="QQ-Plot der Residuen bei mod22")
qqline(mod22$residuals)
par(op)

##drittes Modell
x3 <- dat3[,1]
t3 <- dat3[,2]
plot(t3,x3,main="plot8c.dat",ylim=c(0,20))
#punktwolke sieht "seltsam" aus.
mod3 <- lm(x3~t3)
summary(mod3)
mod31 <- lm(x3 ~ t3 + I(t3^2))
summary(mod31)
plot(mod3$residuals,type="l")
plot(mod31$residuals,type="l")
qqnorm(mod3$residuals)
qqline(mod3$residuals)
dwtest(mod3,alternative="two.sided")
dwtest(mod31,alternative="two.sided")
#also positiv korreliert, versuchen log-Modell
mod33 <- lm(log(x3) ~ t3)
summary(mod33)
mod34 <- lm(log(x3) ~ t3 + I(t3^2))
summary(mod34)
#quadratischer term wohl nicht noetig
plot(mod33$residuals,type="l")
qqnorm(mod33$residuals)
qqline(mod33$residuals)
dwtest(mod33,alternative="two.sided")
#Modell log-normal, also
#x3 ~ exp(t3 + eps) mit unkorrelierten eps (heteroskedastisch)

#zusammen plotten
op <- par(mfrow=c(2,2))
plot(t3,x3,main="plot8c.dat",ylim=c(0,20))
plot(mod3$residuals,type="l",main="Residuen bei lm(x3 ~ t3)")
plot(mod33$residuals,type="l",main="Residuen bei lm(log(x3) ~ t3)")
qqnorm(mod33$residuals,main="QQ-Plot der Residuen bei log-Modell")
qqline(mod33$residuals)
par(op)