Comparing clustering of standard errors between felm and feols functions I'm using the lfe and fixest packages to run regressions with high-dimensional fixed effects. For these regressions, I would like to cluster the standard errors by several dimensions (eg. product, destination and time). However, I'm confused about the syntax and how it differs between the felm and feols commands. Would the clustering in the following two models be equivalent?
EDIT: I ran the two models and found that m2 has larger standard errors than m1.
m1 <- felm(y ~ x1+ x2 | fe1 + fe2|0|product + destination + time, data=df) #with lfe package
summary(m1)
m2 <- feols(y ~ x1+ x2 | fe1 + fe2, data=df) #with fixest package 
summary(m2, cluster=~product + destination + time)


 A: Actually there is no single way to compute the standard-errors. The way they are computed in fixest and how they compare to lfe are explained in this vignette.
There were also a couple of (minor) bugs in the SEs in fixest version < 0.6.0 making the SEs look slightly different.
Here's a comparison related to your example with toy data:
library(fixest) ; library(lfe)
data(trade)

est_felm  =  felm(log(Euros) ~ log(dist_km) | Origin + Destination | 0 | Origin + Destination + Year, trade)
est_feols = feols(log(Euros) ~ log(dist_km) | Origin + Destination, trade)

# Same SEs but different p-values
coeftable(est_felm)
#>               Estimate Cluster s.e.   t value     Pr(>|t|)
#> log(dist_km) -2.072132    0.1516212 -13.66651 2.525297e-07
coeftable(est_feols, cluster = ~ Origin + Destination + Year)
#>               Estimate Std. Error   t value     Pr(>|t|)
#> log(dist_km) -2.072132  0.1516212 -13.66651 2.024366e-42

# Same SEs and p-values (t.df is explained in the vignette)
coeftable(est_feols, cluster = ~ Origin + Destination + Year,
          dof = dof(t.df = "min"))
#>               Estimate Std. Error   t value     Pr(>|t|)
#> log(dist_km) -2.072132  0.1516212 -13.66651 2.525297e-07

