# Monitoria de Econometria 2 - 01/04/2022
# Monitor: Alan Leal/ Professor: Raphael Corbi

##### Pacotes:
if(!require(tidyverse)){install.packages("tidyverse");require(tidyverse)}
if(!require(wooldridge)){install.packages("wooldridge");require(wooldridge)}
if(!require(mvtnorm)){install.packages("mvtnorm");require(mvtnorm)}
if(!require(multiwayvcov)){install.packages("multiwayvcov");require(multiwayvcov)}

##### Matrizes, dataframes e outros objetos no R
a <- 5 # inteiro
b <- "String"
# Matriz:
a <- 1:10 %>% matrix(.,ncol=1); seq(1,10,2); rep(10,5)
b <- 31:40 %>% matrix(.,ncol=1)
# Cbind e rbind
# cbind une vetores na direção da coluna
cbind(a,b)
#rbind une vetores na direção da linha
rbind(a,b)
# c() coloca num vetor os elementos
c(1,2,3,4,5)
# matrix cria matrizes no R
matrix(cbind(a,b,rep(c(1,2,3,4,5),2)),ncol=3)
# diag cria uma matriz diagonal cuja diagonal principal é preenchida pelo vetor argumento
diag(a)

##### Operações de matrizes 
# Soma:
a+b
# Produto de matrizes
a%*%t(b); det(a%*%t(b));solve(a%*%t(b))
a*b # multiplicação elemento por elemento


##### Funções no R
# Uma função precisa ter nome e argumentos, mas não necesariamente precisa retornar algum objeto:
eq2 <- function(a,b,c){
  # a, b e c precisam ser números reais
  # b^2-4*a*c maior ou igual a 0
  r1 <- (-b+sqrt(b^2-4*a*c))/(2*a)
  r2 <- (-b-sqrt(b^2-4*a*c))/(2*a)
  return(cbind(r1,r2))
}

# Chamar a função:
eq2(1,0,-16)

##### Loops 
for (i in 1:10){
  print(20*i)
}

##### Estimativa de um modelo de regressão linear via MQO
# Supondo erros iid e clustered
# Lendo a base de dados Census2000 do pacote Wooldridge:
data('census2000')

# Vamos estimar uma regressão linear explicando o log do salário semanal em função de:
# Educação
# Experiência
# Experiência^2

# Vamos construir as matrizes necessárias:
Y <- matrix(census2000$lweekinc,ncol=1)
n <- dim(Y)[1]
X <- matrix(cbind(rep(1,n),census2000$educ,census2000$exper,census2000$expersq),ncol=4)
k <- dim(X)[2]
# O vetor beta será dado por:
beta <- solve(t(X)%*%X)%*%t(X)%*%Y
# Agora, vamos calcular o erro-padrão assumindos erros iid:
u <- Y-X%*%beta
sigma2 <- as.numeric((t(u)%*%u)/(n-k))
Var_beta <- sigma2*solve(t(X)%*%X)
sd_beta <- sqrt(diag(Var_beta))
Var_beta %>% diag() %>% sqrt() # Outra forma de obter a mesma coisa
# %>% é do dplyr e facilita a passagem de argumentos para funções

# Calcular a estatística t e o p-valor:
tcalc <- beta/sd_beta
# E o p-valor:
p_valor <- 2 * (1 - pt(q = abs(tcalc), df = n-k))
# Comparando os resultados da regressão original:
summary(lm(lweekinc~educ+exper+expersq,data=census2000))

# Eu quero criar uma função que leia um vetor coluna Y e uma matriz de variáveis explicativas X
# E que me retorne o vetor de coeficientes, os erros-padrãos desses coeficientes, as estatísticas t e os p-valores
mqo_my <- function(y,X){
  beta <- solve(t(X)%*%X)%*%t(X)%*%Y
  # Agora, vamos calcular nosso erro-padrão do coeficiente beta:
  n <- dim(X)[1]
  k <- dim(X)[2]
  sigma2 <- as.numeric((t(Y-X%*%beta))%*%(Y-X%*%beta))/(n-k)
  Var_beta <- sigma2*solve(t(X)%*%X)
  sd_beta <- sqrt(diag(Var_beta))
  # Teste t e p-valor:
  # O teste t tem por estatísitica: (beta[i]-0)/sd[i]:
  tcalc <- beta/sd_beta
  # Por fim, o p-valor tem por cálculo:
  p_valor <- 2 * (1 - pt(q = abs(tcalc), df = n-k))
  return(cbind(beta,sd_beta,tcalc,p_valor))
}

mqo_my(Y,X)


##### MQo com erros por cluster em alguma variável
u <- cbind(Y-X%*%beta,census2000$puma)
ids <- unique(census2000$puma)
omega <- matrix(0,ncol = dim(X)[2],nrow=dim(X)[2])
X2 <- matrix(c(rep(1,n),census2000$educ,census2000$exper,census2000$expersq,census2000$puma),ncol=5)
for (i in 1:length(unique(census2000$puma))){
  x_restrito <- matrix(X2[X2[,5]==ids[i],],ncol=5);x_restrito <- matrix(x_restrito[,-5],ncol=4)
  u_restrito <- matrix((u[u[,2]==ids[i],]),ncol=1);u_restrito <- matrix(u_restrito[1:(length(u_restrito)/2),],ncol=1)
  omega1 <- t(x_restrito)%*%u_restrito%*%t(u_restrito)%*%x_restrito
  omega <- omega+omega1
}
# Calculando agora a variância com cluster em puma:
var_beta_cluster <- solve(t(X)%*%X)%*%omega%*%solve(t(X)%*%X)
sd_beta <- sqrt(diag(var_beta_cluster))
# Teste t e p-valor:
# O teste t tem por estatísitica: (beta[i]-0)/sd[i]:
tcalc <- beta/sd_beta
# Por fim, o p-valor tem por cálculo:
p_valor <- 2 * (1 - pt(q = abs(tcalc), df = n-k))
# Comparando com o resultado da regressão original:
ols_cluster <- lm(lweekinc~educ+exper+expersq,data=census2000)
cluster_R <- cluster.vcov(ols_cluster,cluster = census2000$puma) %>% diag() %>% sqrt()
# Comparando nossos resultados e os obtidos automaticamente pelo R, temos:
cbind(sd_beta,cluster_R)




set.seed(12345)# SEMPRE USEM ESSE COMANDO QUANDO RETIRANDO NÚMEROS DE DISTRIBUIÇÕES ALEATÓRIAS

##### Simulações sobre o (não-)uso de erros-padrão por cluster
#Nessa última parte da aula, basicamente a mesma implementação existente em Cunningham, vamos verificar as rejeições falsas de nulas para o intercepto, no caso de dados nos quais as observações se agrupem no formato de cluster. 

gen_cluster <- function(param = c(.1, .5), n = 1000, n_cluster = 50, rho = .5) {
  # Function to generate clustered data
  # Required package: mvtnorm
  
  # individual level
  Sigma_i <- matrix(c(1, 0, 0, 1 - rho), ncol = 2)
  values_i <- rmvnorm(n = n, sigma = Sigma_i)
  
  # cluster level
  cluster_name <- rep(1:n_cluster, each = n / n_cluster)
  Sigma_cl <- matrix(c(1, 0, 0, rho), ncol = 2)
  values_cl <- rmvnorm(n = n_cluster, sigma = Sigma_cl)
  
  # predictor var consists of individual- and cluster-level components
  x <- values_i[ , 1] + rep(values_cl[ , 1], each = n / n_cluster)
  
  # error consists of individual- and cluster-level components
  error <- values_i[ , 2] + rep(values_cl[ , 2], each = n / n_cluster)
  
  # data generating process
  y <- param[1] + param[2]*x + error
  
  df <- data.frame(x, y, cluster = cluster_name)
  return(df)
}

# Simulate a dataset with clusters and fit OLS
# Calculate cluster-robust SE when cluster_robust = TRUE
cluster_sim <- function(param = c(.1, .5), n = 1000, n_cluster = 50,
                        rho = .5, cluster_robust = FALSE) {
  # Required packages: mvtnorm, multiwayvcov
  df <- gen_cluster(param = param, n = n , n_cluster = n_cluster, rho = rho)
  fit <- lm(y ~ x, data = df)
  b1 <- coef(fit)[2]
  if (!cluster_robust) {
    Sigma <- vcov(fit)
    se <- sqrt(diag(Sigma)[2])
    b1_ci95 <- confint(fit)[2, ]
  } else { # cluster-robust SE
    Sigma <- cluster.vcov(fit, ~ cluster)
    se <- sqrt(diag(Sigma)[2])
    t_critical <- qt(.025, df = n - 2, lower.tail = FALSE)
    lower <- b1 - t_critical*se
    upper <- b1 + t_critical*se
    b1_ci95 <- c(lower, upper)
  }
  return(c(b1, se, b1_ci95))
}

# Function to iterate the simulation. A data frame is returned.
run_cluster_sim <- function(n_sims = 1000, param = c(.1, .5), n = 1000,
                            n_cluster = 50, rho = .5, cluster_robust = FALSE) {
  # Required packages: mvtnorm, multiwayvcov, dplyr
  df <- replicate(n_sims, cluster_sim(param = param, n = n, rho = rho,
                                      n_cluster = n_cluster,
                                      cluster_robust = cluster_robust))
  df <- as.data.frame(t(df))
  names(df) <- c('b1', 'se_b1', 'ci95_lower', 'ci95_upper')
  df <- df %>% 
    mutate(id = 1:n(),
           param_caught = ci95_lower <= param[2] & ci95_upper >= param[2])
  return(df)
}

# Distribution of the estimator and confidence intervals
sim_params <- c(.4, 0)   # beta1 = 0: no effect of x on y
sim_nocluster <- run_cluster_sim(n_sims = 10000, param = sim_params, rho = 0)
hist_nocluster <- ggplot(sim_nocluster, aes(b1)) +
  geom_histogram(color = 'black') +
  geom_vline(xintercept = sim_params[2], color = 'red')
print(hist_nocluster)

ci95_nocluster <- ggplot(sample_n(sim_nocluster, 100),
                         aes(x = reorder(id, b1), y = b1, 
                             ymin = ci95_lower, ymax = ci95_upper,
                             color = param_caught)) +
  geom_hline(yintercept = sim_params[2], linetype = 'dashed') +
  geom_pointrange() +
  labs(x = 'sim ID', y = 'b1', title = 'Randomly Chosen 100 95% CIs') +
  scale_color_discrete(name = 'True param value', labels = c('missed', 'hit')) +
  coord_flip()
print(ci95_nocluster)

sim_nocluster %>% summarize(type1_error = 1 - sum(param_caught)/n())

#Data with clusters
sim_params <- c(.4, 0)   # beta1 = 0: no effect of x on y
sim_cluster_ols <- run_cluster_sim(n_sims = 10000, param = sim_params)
hist_cluster_ols <- hist_nocluster %+% sim_cluster_ols
print(hist_cluster_ols)

#Confidence interval
ci95_cluster_ols <- ci95_nocluster %+% sample_n(sim_cluster_ols, 100)
print(ci95_cluster_ols)
sim_cluster_ols %>% summarize(type1_error = 1 - sum(param_caught)/n())

#clustered robust
sim_params <- c(.4, 0)   # beta1 = 0: no effect of x on y
sim_cluster_robust <- run_cluster_sim(n_sims = 10000, param = sim_params,
                                      cluster_robust = TRUE)

hist_cluster_robust <- hist_nocluster %+% sim_cluster_ols
print(hist_cluster_robust)

#Confidence Intervals
ci95_cluster_robust <- ci95_nocluster %+% sample_n(sim_cluster_robust, 100)
print(ci95_cluster_robust)

sim_cluster_robust %>% summarize(type1_error = 1 - sum(param_caught)/n())