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

data (longley)

model <- lm (Employed ~ GNP + Population, data = longley)
summary (model, correlation = TRUE)
r <- model$residuals
n <- length (r)
# graficky 
r1 <- r[-n]
r2 <- r[-1]
cor (r1, r2)
plot (r1, r2, xlab="r_1, ..., r_n-1", ylab = "r_2, ..., r_n", main = "residual plot") 
model.r <- lm (r2 ~ r1) 
abline (model.r$coefficients)
# indexovy graf rezidui
plot (1:n, r, type = "b", xlab="index", ylab = "r", main = "residual plot")
abline (h = 0, lty = 2)

# odhad theta
theta1 <- as.numeric( (t(r1) %*% r2) / (t(r1) %*% r1))
theta1
# Durbin-Watson
DW <- dwtest (Employed ~ GNP + Population, data = longley)
D <- DW$statistic
D
theta2 <- 1 - D/2
theta2
# asymptoticky test
alpha <- 0.05
U <- abs (sqrt(n) * theta1)
c (U, qnorm (1 - alpha/2)) 
U > qnorm (1 - alpha/2)  # Zamitame hypotezu H: theta = 0
DW$p.value # Podle p-hodnoty D-W testu take zamitame 


# Model s AR(1) a reseni zobecnenou metodou nejmensich ctvercu
X <- model.matrix (model)
n <- nrow (X)
W <- diag (n)
W <- theta2^(abs (row (W) - col (W)))
W.inv <- solve (W)
beta <- solve (t(X) %*% W.inv %*% X) %*% t(X) %*% W.inv %*% longley$Employed
beta

model$coefficients

Yhat <- X %*% beta
r <- longley$Employed - Yhat 
n <- length (r)
r1 <- r[-n]
r2 <- r[-1]
cor (r1, r2)


# Reseni pomoci funkce gls
model2 <- gls (Employed ~ GNP + Population, correlation = corAR1 (form = ~ Year), data = longley)
summary (model2, correlation = TRUE)