## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, 
                      eval = TRUE,
                      comment = "#>")
Sys.setenv(lang = "en")

library(fixest)
setFixest_nthreads(1)

## -----------------------------------------------------------------------------
library(fixest)
data(trade)
# OLS estimation
gravity = feols(log(Euros) ~ log(dist_km) | Destination + Origin + Product + Year, trade)
# Two-way clustered SEs
summary(gravity, vcov = "twoway")
# Two-way clustered SEs, without small sample correction
summary(gravity, vcov = "twoway", ssc = ssc(adj = FALSE, cluster.adj = FALSE))

## ----eval = TRUE, include = FALSE---------------------------------------------
is_plm = requireNamespace("plm", quietly = TRUE)

if(!is_plm){
    knitr::opts_chunk$set(eval = FALSE)
    cat("Evaluation of the next chunks requires 'plm' which is not installed.")
} else {
    knitr::opts_chunk$set(eval = TRUE)
    
    library(plm)
    library(sandwich)
}


## ----eval = !is_plm, include = !is_plm----------------------------------------
#  # NOTE:
#  # Evaluation of the next chunks requires the package 'plm' which is not installed.
#  # The code output is not reported.

## -----------------------------------------------------------------------------
library(sandwich)
library(plm)

data(Grunfeld)

# Estimations
res_lm    = lm(inv ~ capital, Grunfeld)
res_feols = feols(inv ~ capital, Grunfeld)

# Same standard-errors
rbind(se(res_lm), se(res_feols))

# Heteroskedasticity-robust covariance
se_lm_hc    = sqrt(diag(vcovHC(res_lm, type = "HC1")))
se_feols_hc = se(res_feols, vcov = "hetero")
rbind(se_lm_hc, se_feols_hc)

## -----------------------------------------------------------------------------

# Estimations
est_lm    = lm(inv ~ capital + as.factor(firm) + as.factor(year), Grunfeld)
est_plm   = plm(inv ~ capital + as.factor(year), Grunfeld, index = c("firm", "year"), model = "within")
# we use panel.id so that panel VCOVs can be applied directly
est_feols = feols(inv ~ capital | firm + year, Grunfeld, panel.id = ~firm + year)

#
# "iid" standard-errors
#

# By default fixest clusters the SEs when FEs are present,
#  so we need to ask for iid SEs explicitly.
rbind(se(est_lm)["capital"], se(est_plm)["capital"], se(est_feols, vcov = "iid"))

# p-values:
rbind(pvalue(est_lm)["capital"], pvalue(est_plm)["capital"], pvalue(est_feols, vcov = "iid"))


## -----------------------------------------------------------------------------
# Clustered by firm
se_lm_firm    = se(vcovCL(est_lm, cluster = ~firm, type = "HC1"))["capital"]
se_plm_firm   = se(vcovHC(est_plm, cluster = "group"))["capital"]
se_stata_firm = 0.06328129    # vce(cluster firm)
se_feols_firm = se(est_feols) # By default: clustered according to firm

rbind(se_lm_firm, se_plm_firm, se_stata_firm, se_feols_firm)

## -----------------------------------------------------------------------------
# How to get the lm version
se_feols_firm_lm = se(est_feols, ssc = ssc(fixef.K = "full"))
rbind(se_lm_firm, se_feols_firm_lm)

# How to get the plm version
se_feols_firm_plm = se(est_feols, ssc = ssc(fixef.K = "none", cluster.adj = FALSE))
rbind(se_plm_firm, se_feols_firm_plm)

## -----------------------------------------------------------------------------
#
# Newey-west
#

se_plm_NW   = se(vcovNW(est_plm))["capital"]
se_feols_NW = se(est_feols, vcov = "NW")

rbind(se_plm_NW, se_feols_NW)

# we can replicate plm's by changing the type of SSC:
rbind(se_plm_NW, 
      se(est_feols, vcov = NW ~ ssc(adj = FALSE, cluster.adj = FALSE)))

#
# Driscoll-Kraay
#

se_plm_DK   = se(vcovSCC(est_plm))["capital"]
se_feols_DK = se(est_feols, vcov = "DK")

rbind(se_plm_DK, se_feols_DK)

# Replicating plm's
rbind(se_plm_DK, 
      se(est_feols, vcov = DK ~ ssc(adj = FALSE, cluster.adj = FALSE)))


## ----eval = TRUE, include = FALSE---------------------------------------------
is_lfe = requireNamespace("lfe", quietly = TRUE)
is_lfe_plm = is_lfe && is_plm
if(is_lfe){
    # avoids ugly startup messages popping + does not require the use of the not very elegant suppressPackageStartupMessages
    library(lfe)
}

## ----eval = !is_lfe_plm, include = !is_lfe_plm, echo = FALSE------------------
#  if(!is_lfe){
#    cat("The evaluation of the next chunks of code requires the package 'lfe' which is not installed")
#  } else {
#    cat("The evaluation of the next chunks of code requires the package 'plm' (for the data set) which is not installed.",
#      "\nThe code output is not reported.")
#  }

## ----eval = is_lfe_plm, warning = FALSE---------------------------------------
library(lfe)

# lfe: clustered by firm
est_lfe = felm(inv ~ capital | firm + year | 0 | firm, Grunfeld)
se_lfe_firm = se(est_lfe)

# The two are different, and it cannot be directly replicated by feols
rbind(se_lfe_firm, se_feols_firm)

# You have to provide a custom VCOV to replicate lfe's VCOV
my_vcov = vcov(est_feols, ssc = ssc(adj = FALSE))
se(est_feols, vcov = my_vcov * 199/198) # Note that there are 200 observations

# Differently from feols, the SEs in lfe are different if year is not a FE:
# => now SEs are identical.
rbind(se(felm(inv ~ capital + factor(year) | firm | 0 | firm, Grunfeld))["capital"],
      se(feols(inv ~ capital + factor(year) | firm, Grunfeld))["capital"])

# Now with two-way clustered standard-errors
est_lfe_2way  = felm(inv ~ capital | firm + year | 0 | firm + year, Grunfeld)
se_lfe_2way   = se(est_lfe_2way)
se_feols_2way = se(est_feols, vcov = "twoway")
rbind(se_lfe_2way, se_feols_2way)

# To obtain the same SEs, use cluster.df = "conventional"
sum_feols_2way_conv = summary(est_feols, vcov = twoway ~ ssc(cluster.df = "conv"))
rbind(se_lfe_2way, se(sum_feols_2way_conv))

# We also obtain the same p-values
rbind(pvalue(est_lfe_2way), pvalue(sum_feols_2way_conv))

## -----------------------------------------------------------------------------
setFixest_ssc(ssc(adj = FALSE))

## -----------------------------------------------------------------------------
setFixest_vcov(no_FE = "iid", one_FE = "iid", 
               two_FE = "hetero", panel = "driscoll_kraay")