Inconsistent MNW (2023) wildcluster bootstrapped p-values with default engine ("R")

[amarbler]

Hi,

I am trying to reproduce the MNW (2023) wildcluster bootstrapped p-values. While I can replicate their results using fwildclusterboot::boottest() with engine = "WildBootTests.jl", the default engine (engine = "R") produces slightly different p-values.

library(data.table)
library(fixest)
library(fwildclusterboot)

url = "https://journaldata.zbw.eu/dataset/2a462673-c91b-4ac5-9e46-6137e4fe64f6/resource/9c452211-4d03-498c-8bc6-3fd177613613/download/mnw-files.zip"

# set path and download and unzip MNW data
dir <- " "

file <- paste0(dir, "MNW23_rep.zip")
download.file(url, file)
system(paste0("unzip ", file, " -d ", dir)) 

data = paste0(dir, "min_wage_teen_hours.csv") |> data.table::fread()

reg = fixest::feols(hours2 ~ mw + black + female + i(age) + i(educ) 
                    | state + year, 
                    data = data)

coeftable(reg, 
          cluster = c("state"), 
          ssc = ssc("fixef.K" = "nested"))[1, ]

#>    Estimate  Std. Error     t value    Pr(>|t|) 
#> -0.15389068  0.06231056 -2.46973654  0.01697872 

# WCR-C Julia
set.seed(9326) # according to MNW (2023)
wcr_c_J = boottest(reg, 
                    B = 999999, # according to MNW (2023)
                    param = "mw", 
                    fe = "state",
                    clustid = "state",
                    nthreads = parallel::detectCores()-1,
                    bootstrap_type = "11",
                    engine = "WildBootTests.jl")

wcr_c_J_p = round(wcr_c_J$p_val, 4)

# WCR-C R
dqrng::dqset.seed(9326)
wcr_c_R = boottest(reg, 
                   B = 999999, 
                   param = "mw", 
                   fe = "state",
                   clustid = c("state"),
                   nthreads = parallel::detectCores()-1,
                   bootstrap_type = "11")

wcr_c_R_p = round(wcr_c_R$p_val, 4)
#
# WCR-S Julia
set.seed(9326)
wcr_s_J = boottest(reg, 
                   B = 999999, 
                   param = "mw", 
                   clustid = "state",
                   fe = "state",
                   nthreads = parallel::detectCores()-1,
                   bootstrap_type = "31",
                   engine = "WildBootTests.jl")
wcr_s_J_p = round(wcr_s_J$p_val, 4)
#
# WCR-S R
dqrng::dqset.seed(9326)
wcr_s_R = boottest(reg, 
                   B = 999999, 
                   param = "mw", 
                   fe = "state",
                   clustid = "state",
                   nthreads = parallel::detectCores()-1,
                   bootstrap_type = "31")
wcr_s_R_p = round(wcr_s_R$p_val, 4)

data.table("source" = c("MNW23", "Julia", "R"),
           "WCR_C" = c("0.0362", wcr_c_J_p, wcr_c_R_p),
           "WCR_S" = c("0.0374", wcr_s_J_p, wcr_s_R_p))

#>    source  WCR_C  WCR_S
#>    <char> <char> <char>
#> 1:  MNW23 0.0362 0.0374
#> 2:  Julia 0.0363 0.0373
#> 3:      R 0.0358 0.0368

With the Julia engine's t-statistics almost matching those from fixest::feols() (using ssc(fixef.K = "nested")):

data.table("source" = c("fixest", "Julia", "R"),
           "t-stat" = c(tstat(reg)["mw"], wcr_c_J$t_stat, wcr_c_R$t_stat))

#>    source    t-stat
#>    <char>     <num>
#> 1: fixest -2.469737
#> 2:  Julia -2.469739
#> 3:      R -2.469774

Could these differences in p-values result from the use of different 'original' test statistics or differences in the small sample correction applied by the two engines? If the small sample correction is indeed the source of the discrepancy, would it be feasible to directly obtain K from fixest to ensure consistency?