Options:

1. The skellam distribution is the distribution of the difference in Poisson variables. You should be able construct a distribution and then get the middle 95% by evaluating inverse cdf (0.025) and icdf(0.975)

https://en.m.wikipedia.org/wiki/Skellam_distribution

2) assume normality, which is pretty reasonable given the sample size. Mean is u1 - u2, and variance is u1 + u2. This is the simplest, but doesn't account for the block correlation.

3) run a Clustered bootstrap: similar to a simulation except you redraw from the sample and draw the clusters together. Then take quantiles of the difference between the two outcomes. This is the most accurate.

See cluster data block bootstrap on wikipedia. It sounds tricky, but it's literally just a couple lines of Python.

Cluster data describes data where many observations per unit are observed. This could be observing many firms in many states, or observing students in many classes https://en.m.wikipedia.org/wiki/Bootstrapping_(statistics)