#===============================================================================================#
# Metody Ekonometryczne
# Zadanie 5. Monte Carlo ilustrujace efektywnosc estymatora WLS
# Michal Gradzewicz
#===============================================================================================#

rm(list=ls())
library(car)
library(lmtest)
library(AER)
library(stargazer)

N <- 250
beta_0_true <- 1.5
beta_1_true <- 1
kappa <- 1
sigma <- 1

steps <- 1000

#DGP

x <- matrix(0, nrow = N, ncol = steps)
for (i in seq(steps)) {
  x[ , i] <- rnorm(N, 2, 2)
}

y <- matrix(0, nrow = N, ncol = steps)
for (i in seq(steps)) {
  for (j in seq(N)) {
    y[j , i] <- beta_0_true + beta_1_true*x[j, i] + rnorm(1, 0, exp(x[j,i] * kappa) * sigma^2)
  }
}

#(i)
beta_0 <- vector(length = steps)
beta_1 <- vector(length = steps)
for (i in seq(steps)) {
  reg <- lm(y[, i] ~ x[, i])
  beta_0[i] <- coef(reg)[1]
  beta_1[i] <- coef(reg)[2]
}
plot(density(beta_0))
mean(beta_0)
sd(beta_0)

plot(density(beta_1))
mean(beta_1)
sd(beta_1)

#(ii)
beta_0 <- vector(length = steps)
beta_1 <- vector(length = steps)
beta_0_wgls <- vector(length = steps)
beta_1_wgls <- vector(length = steps)
beta_0_fwgls <- vector(length = steps)
beta_1_fwgls <- vector(length = steps)
for (i in seq(steps)) {
  reg <- lm(y[, i] ~ x[, i])
  beta_0[i] <- coef(reg)[1]
  beta_1[i] <- coef(reg)[2]
  # (i)
  reg_wgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(x[, i]*kappa)))
  beta_0_wgls[i] <- coef(reg_wgls)[1]
  beta_1_wgls[i] <- coef(reg_wgls)[2]
  # (ii)
  aux.lm <- lm(log(reg$residuals^2) ~ x[ , i])
  reg_fwgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(aux.lm$fitted.values)))
  beta_0_fwgls[i] <- coef(reg_fwgls)[1]
  beta_1_fwgls[i] <- coef(reg_fwgls)[2]
}

plot(density(beta_0_fwgls), col = "red")
lines(density(beta_0_wgls), col = "blue")
lines(density(beta_0), col = "green")

print("beta_0:")
c(ols = mean(beta_0), fwgls= mean(beta_0_fwgls), wgls = mean(beta_0_wgls))
c(ols = sd(beta_0), fgls = sd(beta_0_fwgls), wgls = sd(beta_0_wgls))

plot(density(beta_1_fwgls), col = "red")
lines(density(beta_1_wgls), col = "blue")
lines(density(beta_1), col = "green")

print("beta_1:")
c(ols = mean(beta_1), fwgls= mean(beta_1_fwgls), wgls = mean(beta_1_wgls))
c(ols = sd(beta_1), fgls = sd(beta_1_fwgls), wgls = sd(beta_1_wgls))

#(iii)

for (kappa in c(0.5, 1.5)) {
  #DGP
  x <- matrix(0, nrow = N, ncol = steps)
  for (i in seq(steps)) {
    x[ , i] <- rnorm(N, 2, 2)
  }
  
  y <- matrix(0, nrow = N, ncol = steps)
  for (i in seq(steps)) {
    for (j in seq(N)) {
      y[j , i] <- beta_0_true + beta_1_true*x[j, i] + rnorm(1, 0, exp(x[j,i] * kappa) * sigma^2)
    }
  }
  
  beta_0 <- vector(length = steps)
  beta_1 <- vector(length = steps)
  beta_0_wgls <- vector(length = steps)
  beta_1_wgls <- vector(length = steps)
  beta_0_fwgls <- vector(length = steps)
  beta_1_fwgls <- vector(length = steps)
  for (i in seq(steps)) {
    reg <- lm(y[, i] ~ x[, i])
    beta_0[i] <- coef(reg)[1]
    beta_1[i] <- coef(reg)[2]
    # (i)
    reg_wgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(x[, i]*kappa)))
    beta_0_wgls[i] <- coef(reg_wgls)[1]
    beta_1_wgls[i] <- coef(reg_wgls)[2]
    # (ii)
    aux.lm <- lm(log(reg$residuals^2) ~ x[ , i])
    reg_fwgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(aux.lm$fitted.values)))
    beta_0_fwgls[i] <- coef(reg_fwgls)[1]
    beta_1_fwgls[i] <- coef(reg_fwgls)[2]
  }

  cat("kappa = ", kappa, "\n")
  print("beta_0:")
  print(c(ols = mean(beta_0), fwgls= mean(beta_0_fwgls), wgls = mean(beta_0_wgls)))
  print(c(ols = sd(beta_0), fgls = sd(beta_0_fwgls), wgls = sd(beta_0_wgls)))
  print("beta_1:")
  print(c(ols = mean(beta_1), fwgls= mean(beta_1_fwgls), wgls = mean(beta_1_wgls)))
  print(c(ols = sd(beta_1), fgls = sd(beta_1_fwgls), wgls = sd(beta_1_wgls)))
}


#(iv)

kappa <- 1

for (N in c(30, 100, 10000)) {
  #DGP
  x <- matrix(0, nrow = N, ncol = steps)
  for (i in seq(steps)) {
    x[ , i] <- rnorm(N, 2, 2)
  }
  
  y <- matrix(0, nrow = N, ncol = steps)
  for (i in seq(steps)) {
    for (j in seq(N)) {
      y[j , i] <- beta_0_true + beta_1_true*x[j, i] + rnorm(1, 0, exp(x[j,i] * kappa) * sigma^2)
    }
  }
  
  beta_0 <- vector(length = steps)
  beta_1 <- vector(length = steps)
  beta_0_wgls <- vector(length = steps)
  beta_1_wgls <- vector(length = steps)
  beta_0_fwgls <- vector(length = steps)
  beta_1_fwgls <- vector(length = steps)
  for (i in seq(steps)) {
    reg <- lm(y[, i] ~ x[, i])
    beta_0[i] <- coef(reg)[1]
    beta_1[i] <- coef(reg)[2]
    # (i)
    reg_wgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(x[, i]*kappa)))
    beta_0_wgls[i] <- coef(reg_wgls)[1]
    beta_1_wgls[i] <- coef(reg_wgls)[2]
    # (ii)
    aux.lm <- lm(log(reg$residuals^2) ~ x[ , i])
    reg_fwgls <- lm(y[, i] ~ x[, i], weights = 1/sqrt(exp(aux.lm$fitted.values)))
    beta_0_fwgls[i] <- coef(reg_fwgls)[1]
    beta_1_fwgls[i] <- coef(reg_fwgls)[2]
  }
  
  cat("N = ", N, "\n")
  print("beta_0:")
  print(c(ols = mean(beta_0), fwgls= mean(beta_0_fwgls), wgls = mean(beta_0_wgls)))
  print(c(ols = sd(beta_0), fgls = sd(beta_0_fwgls), wgls = sd(beta_0_wgls)))
  print("beta_1:")
  print(c(ols = mean(beta_1), fwgls= mean(beta_1_fwgls), wgls = mean(beta_1_wgls)))
  print(c(ols = sd(beta_1), fgls = sd(beta_1_fwgls), wgls = sd(beta_1_wgls)))
}