Estimate the full model, including the parameter of interest. Keep the t-statistic, using analytically clustered standard errors.
Re-estimate the model, imposing the null hypothesis of no effect. The simplest way to do this is to just re-estimate the model, but omit the parameter of interest. Collect the fitted values and residuals for each observation.
Construct a bootstrap replicate for each cluster. That is, with equal probability, for all observations within a cluster, multiply the residual by +1 or -1. Given the new bootstrap residual, add this to the fitted values to construct a bootstrap replicate value for the dependent variable.
Re-compute the full model, again using analytically clustered standard errors, and keep the resulting t-statistic.
Repeat steps 3 and 4 at least 500 times.
Use the bootstrap distribution of the t-statistic to test whether the original t-statistic is significant.
Anyway, it seemed to me that it would be much faster to re-compute the full model in step 4 if I could tell reghdfe to use the estimated parameters, and values of the fixed effects, from step 1 as seed values. Otherwise it starts from scratch, taking 11 iterations or whatever. I guess the fixed effects are not actually estimated, so I'm not sure how one would do this, but this would be the general idea.
Tentative steps: