# Roteiro da aula de painel
# Prof. Raphael B. Conbi/ Monitor: Alan Leal

# Lendo os pacotes necessários:
if(!require(tidyverse)){install.packages("tidyverse");require(tidyverse)}
if(!require(estimatr)){install.packages("estimatr");require(estimatr)}
if(!require(plm)){install.packages("plm");require(plm)}
if(!require(stargazer)){install.packages("stargazer");require(stargazer)}

# Agora, lendo a base de dados usada aqui:
dairy <- read_csv("https://pages.stern.nyu.edu/~wgreene/Econometrics/dairy.csv") %>% arrange(FARM,YEAR)

# Parte I) MQO empilhado, efeitos fixos e demean.
# X1- log de qtd de vacas
# X2 - log de área de terra
# X3 - log do trabalho
# X4 - log da alimentação
mqo_a <- lm(MILK~X1+X2+X3+X4,data=dairy)
ef_a <- lm_robust(MILK~X1+X2+X3+X4,# usando estimatr
                  data = dairy, 
                  fixed_effect = ~FARM, 
                  se_type = "stata")
ef_b <- plm(MILK~X1+X2+X3+X4,# usando plm
            data=dairy,
            effect = "individual",
            model="within",
            index = c("FARM","YEAR"))
# Por fim, fazendo o MQO retirado de suas médias:
dairy_demean <- dairy %>% 
  mutate(MILK=MILK-ave(MILK,FARM),
         X1=X1-ave(X1,FARM),
         X2=X2-ave(X2,FARM),
         X3=X3-ave(X3,FARM),
         X4=X4-ave(X4,FARM))
mqo_demean <- lm(MILK~-1+X1+X2+X3+X4,data=dairy_demean)
# O Cunningham não corrige os erros padrão dessa regressão retirada da média. Isso é errado!
# Por isso, é melhor trabalhar com a regressão de efeitos fixos programada no programa.
# Agora, vamos comparar todas essas estimações:
stargazer(mqo_a,ef_b,mqo_demean,type="text")

# Parte II) MQO de variáveis dummies equivale ao  modelo de efeitos fixos:
mqo_dummies <- lm(MILK~X1+X2+X3+X4+as.factor(FARM),data=dairy)
stargazer(ef_b,mqo_dummies,type="text")

# Parte III) Primeiras diferenças para dois períodos equivale aos efeitos fixos para dois períodos (para 1997 e 1998)
dairy_fd <- dairy %>% 
  filter(YEAR>=97) %>% group_by(FARM) %>% 
  mutate(MILK=c(NA,diff(MILK)),
         X1=c(NA,diff(X1)),
         X2=c(NA,diff(X2)),
         X3=c(NA,diff(X3)),
         X4=c(NA,diff(X4))) %>% na.omit()

frst_dif <- lm(MILK~-1+X1+X2+X3+X4,data=dairy_fd)
# Agora, calculando o modelo de efeitos fixos apenas para 1997 e 1998:
dairy_9798 <- dairy %>% filter(YEAR>=97)
ef_b_1 <- plm(MILK~-1+X1+X2+X3+X4,# usando plm
              data=dairy_9798,
              effect = "individual",
              model="within",
              index = c("FARM","YEAR"))
stargazer(ef_b_1,frst_dif,type="text")