Estimating variance of dgp that measured as Bernoulli Parameter I have some parameters $p_1,\ldots,P_n$ distributed iid. with mean and variance $\mu, \sigma^2$.
Then, for each $t\in\{1,\ldots,n\}$ I have $\{X_{1, t}, \dots X_{N, t}\}$ realizations of the Bernoulli experiment with probability $p_t$ (iid conditional on $p_t$), with whom I estimate 
$$ \tilde p_t = \frac{1}{N} \sum_{i=1}^N X_{i, t} $$
Then, I estimate $\sigma^2$ using $n$ many estimates
$$ \tilde\sigma^2 = \frac{1}{n-1} \sum_{t=1}^n (\tilde p_t - \bar{\tilde{p_t}})^2$$
Is this estimator unbiased? In particular, I'm worried that for small $N$ and high $\mu$, my estimate is biased downwards. My intuition is that "In that case, under small $N$, realizations $X_1$ unlikely to reflect the 0 event. Hence $\tilde p_t$ will be biased towards $1$, and having too little variance. 
 A: The estimator proposed in OP (hereafter "original estimator") is biased - it overestimates (on average) the variance. Intuitively, the variation in $\tilde{p}_t$ consist  of both the variation $p_t$ and the extra variation in the estimates. In the remainder of the answer, this is shown by using the law of total variance. Then, by further observing the law-of-total-variance decomposition, a bias correction is derived. Finally, bias of the original estimator and unbiasedness of the corrected estimator are verified in a computer simulation (Python 3 code). 
Note that the (mistaken) intuition for underestimation in OP seems to be based on $\tilde{p}_t$ being somehow "biased towards $1$". However, $\tilde{p}_t$ is an unbiased estimator of $p_t$. 
 Proof of the bias
First, note that the proposed estimator is the sample variance of $\tilde{p}_t$ which (due to the $n-1$ correction) is known to be an unbiased estimator of the variance of $\tilde{p}_t$. That is, 
\begin{equation}
\mathrm{E}\left[\textrm{Original estimator}\right] = \mathrm{Var}\left[\tilde{p}_t\right].
\end{equation}
Now, apply the law of total variance:
\begin{equation}
= \mathrm{Var}\left[\mathrm{E}\left[\tilde{p}_t \mid p_t\right]\right] + \mathrm{E}\left[\mathrm{Var}\left[\tilde{p}_t \mid p_t\right]\right]
\end{equation}
Use the fact that $\tilde{p}_t$ is an unbiased estimator of $p_t$, i.e., $\mathrm{E}(\tilde{p}_t \mid p_t) = p_t$:
\begin{align}
&= \mathrm{Var}[p_t]+  \mathrm{E}[\mathrm{Var}[\tilde{p}_t \mid p_t]] \\
&= \sigma^2 + \mathrm{E}[Var[\tilde{p}_t  \mid p_t]].
\end{align}
The second term is positive since the estimator $\tilde{p}_t$ is not a constant given the parameter. Thus, it has been shown that
\begin{equation}
\mathrm{E}[\textrm{original estimator}] > \sigma^2
\end{equation}
Deriving an unbiased estimator
So, the bias of the original estimator is $\mathrm{E}[\mathrm{Var}[\tilde{p}_t  \mid p_t]]$. Let us investigate this term close:
\begin{equation}
\mathrm{E}\left[\mathrm{Var}\left[\tilde{p}_t  \mid p_t\right]\right] = \mathrm{E}\left[\mathrm{Var}\left[X_{i,t} \mid p_t\right]/N\right]  
\end{equation}
For the right hand side, we can use the sample variance of $X_{1:N, t}$ that is an unbiased estimator of the variance - so:
\begin{equation}
\mathrm{E}\left[\frac{\sum_{i=1}^N (X_{i,t} - \tilde{p}_t)^2}{N-1}  \mid p_t \right] =  \mathrm{Var}\left[X_{i,t} \mid p_t\right]
\end{equation}
taking expectations and multiplying by $1/N$
\begin{equation}
\mathrm{E}\left[\frac{1}{N}\,\frac{\sum_{i=1}^N (X_{i,t} - \tilde{p}_t)^2}{N-1}\right] = \mathrm{E}\left[\frac{1}{N}\,\mathrm{Var}(X_{i,t} \mid p_t)\right].
\end{equation}
The expectation does not change  by taking average over all $t$, so we have also
\begin{equation}
\mathrm{E}\left[\frac{\sum_{t=1}^n\sum_{i=t}^N (X_{i,t} - \tilde{p}_t)^2}{n\,N\,(N-1)}\right] 
= \mathrm{E}\left[\mathrm{Var}\left[\tilde{p}_t  \mid p_t\right]\right].
\end{equation}
(for unbiasedness could as well have used just one arbitrary $t$ but I assume this decreases the variance of the estimator, no proof though)
We have now obtained
\begin{equation}
E\left[\textrm{original estimator} - \frac{\sum_{t=1}^n\sum_{i=t}^N (X_{i,t} - \tilde{p}_t)^2}{n\,N\,(N-1)}\right] = \sigma^2.
\end{equation}
(Looks like this new estimator could be simplified furher but I did not immediately see anything)
Simulation verification
Let us simulate $1000$ replications where $p_t \sim \mathrm{Beta}(1, 2)$, $n=5$, $N=3$ (these choices have no significance, just picked some small numbers and some distribution whose variance I know and that is easy simulate from, but which does not feel likely to be some special case. The last criterion excludes for example, uniform distribution and $N=n=2$).
Python 3 code which I hope is self-explaining - some places should be vectorized but I expect for loops to be easier to understand for people not using Python. 
import numpy as np
def simulate_data(N,n, alpha, beta):

    p = np.random.beta(alpha, beta, size=n)

    X = np.zeros((N,n))
    for t in range(n):
        for i in range(N):
            X[i,t] = np.random.binomial(n=1, p=p[t])
    return X

def estimator(X):
    N = X.shape[0]
    n = X.shape[1]
    tilde_p_t = np.zeros(n)

    for t in range(n):
        tilde_p_t[t] = np.mean(X[:, t])


    bar_tilde_p_t = np.mean(tilde_p_t)

    return np.sum((tilde_p_t - bar_tilde_p_t)**2) / (n-1)

def correction_term(X):
    N = X.shape[0]
    n = X.shape[1]

    var_X_t = np.zeros(n)
    for t in range(n):
        var_X_t[t] = np.sum((X[:,t] - np.mean(X[:,t]))**2) / (N-1)

    return np.mean(var_X_t) / N

np.random.seed(1)
replications = 10000

estimates = []
corrected_estimates = []
alpha = 1
beta = 2
true_variance = (alpha*beta)/((alpha+beta)**2 * (alpha+beta+1))

for i in range(replications):
    X = simulate_data(N=3, n=5, alpha=alpha, beta=beta)
    estimates.append(estimator(X))
    corrected_estimates.append(estimator(X) - correction_term(X))

print("Theoretical variance")
print(true_variance)
print("Mean of estimates with original estimator")
print(np.mean(estimates))
print("Mean of corrected estimates")
print(np.mean(corrected_estimates))

For the results I get:
Theoretical variance
0.05555555555555555
Mean of estimates with original estimator
0.111303333333
Mean of corrected estimates
0.05575

