I am working on implementing a Fisher Exact Test for some unemployment and wage data. The idea is to describe two populations (one receiving an assistance program (the "treatment") and one not receiving it) via a summary statistic.

The first statistic that I am using is the difference in average wage-level post-treatment between the two populations.

To calculate a p-value for such a set-up, the idea is to make a large number $N$ of random permutations of the assignment vector (the list of 0's and 1's that indicates whether a given observation belongs to the control population or not). I calculate the summary statistic for each of these random assignments, to get a distribution of summary statistics in counterfactual experiments where the control group had been different.

Since I also have the summary statistic for the observed assignment vector (the actual vector of 0's and 1's for the experiment that was truly observed), I can then count the number of simulated summary statistics which are more extreme than the observed summary statistic.

This proportion out of all of my simulated trials serves as the estimated p-value for my summary statistic.

My question is as follows: is there a standard way to get a standard error for such an estimated p-value? Obviously, I can calculate the Monte Carlo standard error for the summary statistics in my simulations, but since all of the simulations are used for just one single p-value calculation, it''s not clear how to get the standard error for that.

I had the following thought about what the standard error might be in this case.

In usual Monte Carlo, we have some function $f$ that we compute at each of the simulated draws $x_{i}$ (where here, the $x_{i}$ are understood to be the assignment vectors). If I define:

$$ f(x_{i}) = \mathbb{1}_{ |stat(x_{i})| > |stat(x_{obs})| } $$

then it seems that the p-value I calculate is given by

$$ \hat{p} = \frac{1}{N}\sum_{i=1}^{N}f(x_{i}) = \frac{\textrm{# more extreme}}{\textrm{total samples}}$$

And following the usual Monte Carlo formulas, would it then make sense to write the variance of the estimate as:

$$ Var(\hat{p}) = \frac{1}{N}\sum_{i=1}^{N}[ f(x_{i}) - \hat{p} ]^{2} $$

and then take the square root to obtain the standard error?

The reason this confuses me is that for each $i$, $f(x_{i})$ will be binary, either the computed statistic was more extreme in that iteration or it wasn't. It seems like it would be error prone to sum up a bunch of binary things like that to estimate the variance in a p-value, but that could just be my unfamiliarity with this method.

Can anyone confirm that this is right? Also, if I exposed any other ignorance about what I am doing here, corrective comments are appreciated.