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


tabulka <- read.csv2 (file = "voda.csv")
summary (tabulka)


model <- lm (Y ~ x, data = tabulka)
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 (Y ~ x, data = tabulka)
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 


# Odstraneni AR(1)
Y.new <- tabulka$Y[-1] - theta2 * tabulka$Y[-n]
x.new <- tabulka$x[-1] - theta2 * tabulka$x[-n]
model.new <- lm (Y.new ~ x.new)
r.new <- model.new$residuals
n <- length (r.new)
# vykresleni novych residui
r1 <- r.new[-n]
r2 <- r.new[-1]
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.new, type = "b", xlab="index", ylab = "r", main = "residual plot")
abline (h = 0, lty = 2)
# Durbin-Watson
DW <- dwtest (Y.new ~ x.new)
D <- DW$statistic
D
# asymptoticky test
theta1 <- as.numeric( (t(r1) %*% r2) / (t(r1) %*% r1))
theta1
U <- abs (sqrt(n) * theta1)
c (U, qnorm (1 - alpha/2)) 
U > qnorm (1 - alpha/2)  # nezamitame hypotezu H: theta = 0
DW$p.value # Podle p-hodnoty D-W testu take nezamitame 



# 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 %*% tabulka$Y
beta

model$coefficients

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


# Reseni pomoci funkce gls
model2 <- gls (Y ~ x, correlation = corAR1 (form = ~ 1), data = tabulka)
summary (model2, correlation = TRUE)


# grafy predikovanych regresnich primek pro 3 modely
xx <- seq (300, 1650, by = 1)
Y <- predict (model, data.frame (x = xx))
Y1 <- cbind (rep (1, times = length (xx)), xx) %*% beta
Y2 <- predict (model2, data.frame (x = xx))
plot (range (xx), c (50, 400), type = "n", xlab = "x", ylab = "Y")
points (tabulka$x, tabulka$Y)
lines (xx, Y, col = "red")
lines (xx, Y1, col = "green")
lines (xx, Y2, col = "blue")