Should bootstrapping and collecting sample means from a series of binomial distributions result in standard normal? I’m trying to better understand several statistical concepts (bootstrapping, central limit theorem, and confidence intervals) by applying them to a binomial distribution (you can think of it as a coin flip, for example). 
I’ll explain my expectations/sanity checks, and then hopefully someone can tell me why my expectations don’t match the results I’m getting in my simulation in code.
Procedure:


*

*Given a probability of success (e.g. $p=0.5$), I obtain the results of 1000 bernoulli trials.

*I do this 1000 times, which means I get 1000 resamples.  From each resample, I compute the sample mean (p=% of positives from the 1000 bernoulli trials), and the sample variance ($n*p*(1-p)$).

*I can then construct a “bootstrap distribution” which is the distribution of the sample means.

*construct test statistic for $S_n$, standardized and normalized: $$\frac{S_n - p}{\sigma/\sqrt{1000}}$$, where $p$ is the "true" probability of success.

*Determine the desired $\alpha$ (type I error rate), and get corresponding $Z$ values:  $Z_{\alpha/2}$ and $Z_{1-\alpha/2}$ (for $\alpha=0.05$ I determined these to be $+/- 1.96$.  

*Determine if the test statistic from step 4) is within the $Z$-values: $Z_{\alpha/2} \le \frac{S_n - p}{\sigma/\sqrt{1000}} \le Z_{1-\alpha/2}$

*If so, then we correctly accept the null hypothesis that the probability of success is $p$.  If not, then we reject the null hypothesis, and call it a type I error.

*Fix $p$ (true probability of success), and $\alpha$ (type I error rate), and repeat steps 1) through 7) 1000 times, and we would expect 50 type I errors.  


Expectation (with justification):


*

*I expect the bootstrap distribution to have the following:
a. mean ($S_n$) = mean of sample means of the resamples (by the law of large numbers).
b. variance (($\sigma^2)/1000$) = mean of sample variances (by the law of large numbers) divided by 1000 (number of resamples).
c. Normally distributed (by central limit theorem, see below)

*I expect the central limit theorem to apply here because (see “classic CLT” described here):
a. The sample means are i.i.d

*If the two expectations above are true, then I would expect 95% of the sample means to fall within the confidence interval for the normal distribution parameterized by $S_n$ and $\sigma^2/1000$

*I expect the percentage of type I errors to be 5% (50 out of 1000).  


    import numpy as np
    from scipy import stats

    def run_test(error_count):
        n=1000 # number of bernoulli trials. one set of bernoulli trials is a resample
        size=1000  # number of resamples
        theta=0.5
        sample_means = []
        sample_variances = []

        # 1) 
        bootstrap_resamples = np.random.binomial(n, theta, size)

        # 2)
        for resample in bootstrap_resamples:
            sample_mean = resample/float(n)
            sample_variance = float(n)*sample_mean*(1-sample_mean)
            sample_means.append(sample_mean)
            sample_variances.append(sample_variance)

        # 3)
        # sample_means (see above) 

        # 4) 
        S_n = np.mean(sample_means)
        sigma_squared = np.mean(sample_variances)
        sigma = np.sqrt(sigma_squared)
        test_statistic  = (S_n - theta)/(sigma/np.sqrt(1000))

        # 5) 
        alpha=0.05
        if not stats.norm.ppf(alpha/2.0) <  test_statistic < stats.norm.ppf(1-(alpha/2.0)):
            error_count +=1
        return error_count


    error_count = 0
    for i in range(1000):
        error_count = run_test(error_count)

    print error_count

The above code consistently returns 0 errors.  In other words, print error_count returns 0, when I would expect it to return approximately 50.  
In most cases it might be good to have 0 errors, but I'm trying to validate that the $\alpha$ level influences the number of errors, and I'm not seeing that here.  
I have a suspicion that the test statistic is computed incorrectly, but I can't figure out what it is...Do you perhaps see an issue in my computation?  
Thank you.
 A: Let's first describe with mathematical notation what happens in the simulation (when dividing by n) and then lay out the statistics that would lead one to expect the result. 
bootstrap_resamples = np.random.binomial(n, theta, size)

Generates $S=1000$ binomial random variables, $Y_s$ each consisting of the sum of $N=1000$ Bernoulli trials with success probability $\theta$, $Y_s = \sum_j^{N} X_{s, j}$
sample_mean = resample/float(n)
sample_variance = sample_mean*(1-sample_mean) / n

Calculates $\hat{\theta}_s = \frac{1}{N} \sum_j^{N} X_{s, j}$ and $\hat{\text{v}}_s = \frac{\hat{\theta_s}(1-\hat{\theta_s})}{N}$
  S_n = np.mean(sample_means)
  sigma_squared = np.mean(sample_variances)
  sigma = np.sqrt(sigma_squared)
  test_statistic  = (S_n - theta)/(sigma/np.sqrt(1000))

Calculates $\hat{\theta} = \frac{1}{S} \sum_s^S \hat{\theta}_s$ and $\hat{\text{v}} = \frac{1}{S} \sum_s^S \hat{\text{v}_s}$ then $\text{stat} = \frac{ \hat{\theta} - \theta}{ \sqrt{ S^{-1} \hat{\text{v}} } }$
       if not stats.norm.ppf(alpha/2.0) <  test_statistic < stats.norm.ppf(1-(alpha/2.0)):
        error_count +=1

Finally, compute the indicator variable $Q = \text{stat} \not \in (-1.96, 1.96)$
Now, what is the sampling distribution of $Q$? If $\text{stat}$ is approximately standard normal, then it approximately equals $5\%$. This is in fact the case because a central limit theorem applies.
Looking at $\hat{\theta} = \frac{1}{S \cdot N} \sum_s \sum_j X_{s,j}$ that is just an average of $S\cdot N$ Bernoulli variables. If we subtract its mean divide by its standard deviation, it will have mean zero and variance 0. It consists of a sum of iid random variables, so will be approximately normal in large samples and be within (-1.96, 1.96) $95\%$ of the time.  We have $E[ \hat{\theta} ] = \theta$ and $Var( \hat{\theta} ) = \text{V} =  \frac{\theta(1-\theta)}{S\cdot N} $. 
So $$\frac{ \hat{\theta} - \theta }{ \sqrt{\text{V}}}$$ is approximately standard normal by a CLT. In the simulation, the approximation will be really good because the sum is over a million random variables. If you want to focus only on understanding the CLT, divide by $v$ in the simulation. And you could have just looked at the mean of 1000 bernoulli trials instead of the mean over 1000 means. 
However, the denominator of $\text{stat}$ isn't $\sqrt{V}$ and yet the simulation still yields $Q$ at around $5\%$ anyways. Why? Because the denominator of $\text{stat}$ converges in probability to the root of $V$ by a law of large numbers. And there's another theorem (https://en.wikipedia.org/wiki/Slutsky%27s_theorem) that says asymptotically, that's as good as if it it actually was $V$. To see the denominator works out:
The denominator of $\text{stat}$ is the root of $\frac{1}{S \cdot S} \sum_s^S \hat{\text{v}_s} = \frac{1}{S \cdot S} \sum_s^S \frac{\hat{\theta_s}(1-\hat{\theta_s})}{N} = \frac{1}{S \cdot N} \left( \frac{1}{S} \sum_s \hat{\theta_s}  - \frac{1}{S}  \sum_s \hat{\theta_s}^2 \right)$. The quantities inside the brackets will converge to $\theta$ and $\theta^2$ by a law of large numbers, so the denominator will converge to the root of $V$ and Slutsky's theorem kicks in. 
