# Load packagesa
library(wooldridge) # The data are stored here
library(sandwich)
library(nlme)       # The gls function is here

data("wage1") # Load data
str(wage1)    # Structure

# Simple Exploratory data analysis  on the variable of interest
summary(wage1[,1:4])
par(mfrow=c(1,4))
for(i in 1:4) hist(wage1[,i], main = names(wage1)[i], breaks = 30)

library(PerformanceAnalytics)
chart.Correlation(wage1[,c(2,3,4,1)])

# Fit the linear model and explore the residuals
ols <- lm(wage ~ educ + exper + tenure, data = wage1)
summary(ols)

plot(ols, 1)

# Breusch - Pagan Test
library(lmtest)
bptest(ols)
help(bptest)


residuals(ols) # y_i - hat y_i
# Bartlett Test

e_std <- rstandard(ols)
G <- 5
min(fitted(ols))
max(fitted(ols))
cuts <- quantile(fitted(ols), probs = seq(0, 1, length.out = G + 1))
grp <- cut(fitted(ols), breaks = cuts, include.lowest = TRUE)
table(grp)

par(mfrow=c(1,1))
plot(fitted(ols), e_std, xlab = "Fitted values",
     ylab = "Standardized residuals")
abline(v = cuts[2:G], lty = 2, col = "red")


bartlett.test(e_std ~ grp) # bartlett test is in the stats library

# By hand 
ng  <- tapply(e_std, grp, length)
S2g <- tapply(e_std, grp, var)
n   <- sum(ng)
S2p <- sum((ng - 1) * S2g) / (n - G)

TB <- ((n - G) * log(S2p) - sum((ng - 1) * log(S2g))) /
  (1 + (sum(1/(ng - 1)) - 1/(n - G)) / (3 * (G - 1)))
TB
pchisq(TB, df = G - 1, lower.tail = FALSE)


#  Fit the GLS assuming that the variance 
# of the errors is proportional to $|\hat y_i|^{2\delta}$
pred <- predict(ols)
glsFit <- gls(wage ~ educ + exper + tenure, data = wage1,
              weights = varPower(form = ~ pred), method = "ML")
summary(glsFit)


# Compare the coefficient (and std err.) OLS vs GLS
cbind(coef(ols), coef(glsFit))
cbind(summary(ols)$coefficients[,2], sqrt(diag(glsFit$varBeta)))


# Plot of residuals 
par(mfrow=c(1,2))
plot(fitted(ols), residuals(ols), xlab = "Fitted values (OLS)",
     ylab = "Residuals", main = "Heteroscedasticity in OLS")
abline(h = 0, col = "red", lty = 2)

plot(fitted(ols), rstandard(ols), xlab = "Fitted values (OLS)",
     ylab = "Standardized Residuals", main = "Heteroscedasticity in OLS")
abline(h = 0, col = "red", lty = 2)

par(mfrow=c(1,1))
plot(fitted(glsFit), resid(glsFit, type="normalized"),     xlab = "Fitted values (GLS)",
     ylab = "Standardized residuals", main = "Variance stabilized in GLS")
abline(h = 0, col = "red", lty = 2)


# Logarithmic transformation
ols_log <- lm(log(wage) ~ educ + exper + tenure, data = wage1)

bptest(ols_log)

e_std <- rstandard(ols_log)
G <- 5
cuts <- quantile(fitted(ols_log), probs = seq(0, 1, length.out = G + 1))
grp <- cut(fitted(ols_log), breaks = cuts, include.lowest = TRUE)
bartlett.test(e_std ~ grp)

# Residuals 
plot(ols_log,1)

pred <- predict(ols_log)
gls_log <- gls(log(wage) ~ educ + exper + tenure, data = wage1,
               weights = varPower(form = ~ pred), method = "ML")


par(mfrow=c(1,2))
plot(fitted(ols_log), rstandard(ols_log), xlab = "Fitted values (OLS)",
     ylab = "Standardized residuals", main = "OLS")
abline(h = 0, col = "red", lty = 2)


plot(fitted(gls_log), resid(gls_log, type="normalized"),     xlab = "Fitted values (GLS)",
     ylab = "Standardized residuals", main = "GLS")
abline(h = 0, col = "red", lty = 2)

# Compare the coefficient (and std err.) OLS vs GLS
cbind(coef(ols_log), coef(gls_log))
cbind(summary(ols_log)$coefficients[,2], sqrt(diag(gls_log$varBeta)))



y ~ a0 + a1*x1 + a3*x3 + a4*x4

y ~ a0 + a1*x1 + a2*x2



# Comparison in terms of RSE
cbind(summary(ols)$sigma, summary(glsFit)$sigma)
cbind(summary(ols_log)$sigma, summary(gls_log)$sigma)


# Information criteria
cbind(AIC(ols), AIC(glsFit))
cbind(AIC(ols_log), AIC(gls_log))

AIC(ols_log) + sum(2*log(wage1$wage))
AIC(gls_log) + sum(2*log(wage1$wage))


# Multicollinearity (an example)

library(ISLR2)
data(Credit)
par(mfrow=c(1,2))
with(Credit, plot(Limit, Age))
with(Credit, plot(Limit, Rating))

# Fit simple linear models 
lmAge <- lm(Balance ~ Age, data = Credit)
summary(lmAge)$coefficients
lmLimit <- lm(Balance ~ Limit, data = Credit)
summary(lmLimit)$coefficients
lmRating <- lm(Balance ~ Rating, data = Credit)
summary(lmRating)$coefficients

# Age and Limit
lmAge <- lm(Balance ~ Age, data = Credit)
summary(lmAge)$coefficients
lmLimit <- lm(Balance ~ Limit, data = Credit)
summary(lmLimit)$coefficients
lmAgeLimit <- lm(Balance ~ Age + Limit, data = Credit)
summary(lmAgeLimit)$coefficients

# Rating and Limit 
lmRating <- lm(Balance ~ Rating, data = Credit)
summary(lmRating)$coefficients
lmLimit <- lm(Balance ~ Limit, data = Credit)
summary(lmLimit)$coefficients
lmLimitRating <- lm(Balance ~ Limit + Rating, data = Credit)
summary(lmLimitRating)$coefficients

# The reason lies in the high correlation 
cor(Credit[, c("Age", "Limit", "Rating")])

# VIF
lmFull <- lm(Balance ~ Age + Limit + Rating, data = Credit)
library(car)
vif(lmFull)

vif(lmAgeLimit)