library (car)
library (nlme)		# gls
library (lmtest)	#	D-W test
library (nortest)

Y <- embed (LakeHuron)
x <- 1875:1972
x <- (x-mean(x))/sd(x)
plot(x,Y,pch=20)
summary(mod <- lm(Y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6)+I(x^7)))
res <- mod$residuals
n <- length(res)
# graficky
plot(res[-n],res[-1],xlab="e_1,...,e_n-1",ylab="e_2,...,e_n",main="residual plot")
model.r <- lm (res[-1] ~ res[-n])
abline (model.r$coefficients)

# odhad theta
(theta1 <- as.numeric((t(res[-n])%*%res[-1])/(t(res[-n])%*%res[-n])))
# Durbin-Watson
#install.packages("lmtest")
library(lmtest)
(ts <- dwtest(mod))
D <- ts$statistic
(theta2 <- 1-D/2)

# asymptoticky test
alpha <- 0.05
abs(sqrt(n)*theta1)
qnorm(1-alpha/2)  # Zamitame hypotezu H: theta = 0

# Odstraneni
Ynew <- Y[-1]-theta1*Y[-n]
xnew <- x[-1]-theta1*x[-n]
summary(modnew <- lm(Ynew~xnew+I(xnew^2)+I(xnew^3)+I(xnew^4)+I(xnew^5)+I(xnew^6)+I(xnew^7)))
summary(modnew <- lm(Ynew~xnew+I(xnew^2)+I(xnew^3)+I(xnew^5)+I(xnew^7)))
resnew <- modnew$residuals
# vykresleni novych residui
plot(resnew[-n+1],resnew[-1],xlab="e_1,...,e_n-1",ylab="e_2,...,e_n",main="residual plot")
# testovani pomoci DW
dwtest(modnew)
# vykresleni
plot(xnew,Ynew)
curve(predict(modnew,data.frame(xnew=x)),add=T,col=2)

# podle vety
ind <- sapply (1:n, function (i) {sapply (1:n, function (j) {abs(i-j)})})
W <- theta1^ind
X <- as.matrix(model.matrix(mod))
Winv <- solve(W)
beta <- solve(t(X)%*%Winv%*%X)%*%t(X)%*%Winv%*%Y
res <- Y-X%*%beta
(theta <- as.numeric((t(res[-n])%*%res[-1])/(t(res[-n])%*%res[-n])))

# asymptoticky test
alpha <- 0.05
abs(sqrt(n)*theta)
qnorm(1-alpha/2)  # Zamitame hypotezu H: theta = 0
# vykresleni
plot(x,Y)
lines(x,X%*%beta,col=2)
curve(predict(mod,data.frame(x=x)),add=T)
legend("bottomleft",legend=c("Puvodni","Novy"),col=c("black","red"),lty=c(1,1))

# normalita residui
shapiro.test(residuals(modnew))

#install.packages("nortest")
shapiro.test(x)
ad.test(x)
cvm.test(x)
pearson.test(x)
sf.test(x)
lillie.test(x)