## S3 method for class 'fixest'
vcov(
  object,
  vcov = NULL,
  se = NULL,
  cluster,
  ssc = NULL,
  attr = FALSE,
  forceCovariance = FALSE,
  keepBounded = FALSE,
  nthreads = getFixest_nthreads(),
  ...
)

# Load panel data
data(base_did)

# Simple estimation on a panel
est = feols(y ~ x1, base_did)

# ======== #
# IID VCOV #
# ======== #

# By default the VCOV assumes iid errors:
se(vcov(est))

# You can make the call for an iid VCOV explicitly:
se(vcov(est, "iid"))

#
# Heteroskedasticity-robust VCOV
#

# By default the VCOV assumes iid errors:
se(vcov(est, "hetero"))

# => note that it also accepts vcov = "White" and vcov = "HC1" as aliases.

# =============== #
# Clustered VCOVs #
# =============== #

# To cluster the VCOV, you can use a formula of the form cluster ~ var1 + var2 etc
# Let's cluster by the panel ID:
se(vcov(est, cluster ~ id))

# Alternative ways:

# -> cluster is implicitly assumed when a one-sided formula is provided
se(vcov(est, ~ id))

# -> using the argument cluster instead of vcov
se(vcov(est, cluster = ~ id))

# For two-/three- way clustering, just add more variables:
se(vcov(est, ~ id + period))

# -------------------|
# Implicit deduction |
# -------------------|
# When the estimation contains FEs, the dimension on which to cluster
# is directly inferred from the FEs used in the estimation, so you don't need
# to explicitly add them.

est_fe = feols(y ~ x1 | id + period, base_did)

# Clustered along "id"
se(vcov(est_fe, "cluster"))

# Clustered along "id" and "period"
se(vcov(est_fe, "twoway"))


# =========== #
# Panel VCOVs #
# =========== #

# ---------------------|
# Newey West (NW) VCOV |
# ---------------------|
# To obtain NW VCOVs, use a formula of the form NW ~ id + period
se(vcov(est, NW ~ id + period))

# If you want to change the lag:
se(vcov(est, NW(3) ~ id + period))

# Alternative way:

# -> using the vcov_NW function
se(vcov(est, vcov_NW(unit = "id", time = "period", lag = 3)))

# -------------------------|
# Driscoll-Kraay (DK) VCOV |
# -------------------------|
# To obtain DK VCOVs, use a formula of the form DK ~ period

se(vcov(est, DK ~ period))

# If you want to change the lag:
se(vcov(est, DK(3) ~ period))

# Alternative way:

# -> using the vcov_DK function
se(vcov(est, vcov_DK(time = "period", lag = 3)))

# -------------------|
# Implicit deduction |
# -------------------|
# When the estimation contains a panel identifier, you don't need
# to re-write them later on

est_panel = feols(y ~ x1, base_did, panel.id = ~id + period)

# Both methods, NM and DK, now work automatically
se(vcov(est_panel, "NW"))
se(vcov(est_panel, "DK"))


# =================================== #
# VCOVs robust to spatial correlation #
# =================================== #

data(quakes)
est_geo = feols(depth ~ mag, quakes)

# ------------|
# Conley VCOV |
# ------------|
# To obtain a Conley VCOV, use a formula of the form conley(cutoff) ~ lat + lon
# with lat/lon the latitude/longitude variable names in the data set
se(vcov(est_geo, conley(100) ~ lat + long))

# Alternative way:

# -> using the vcov_DK function
se(vcov(est_geo, vcov_conley(lat = "lat", lon = "long", cutoff = 100)))

# -------------------|
# Implicit deduction |
# -------------------|
# By default the latitude and longitude are directly fetched in the data based
# on pattern matching. So you don't have to specify them.
# Furhter, an automatic cutoff is deduced by default.

# The following works:
se(vcov(est_geo, "conley"))


# ======================== #
# Small Sample Corrections #
# ======================== #

# You can change the way the small sample corrections are done with the argument ssc.
# The argument ssc must be created by the ssc function
se(vcov(est, ssc = ssc(adj = FALSE)))

# You can add directly the call to ssc in the vcov formula.
# You need to add it like a variable:
se(vcov(est, iid ~ ssc(adj = FALSE)))
se(vcov(est, DK ~ period + ssc(adj = FALSE)))