I am looking for an analytic solution to the variance of the variance of a Bernoulli Random Variable (RV).
The Bernoulli RV probability mass function is defined as:
$$ p(x, \theta) = \theta^x(1-\theta)^{1-x}$$
It's moment generating function is obtained from:
$$M(t) = E[e^{tx}] = \sum_{x=0}^1e^{tx}\theta^x(1-\theta)^{1-\theta} = e^{t0}\theta^0(1-\theta)^{1-0} + e^{t}\theta^1(1-\theta)^{1-1} = (1-\theta) + e^t\theta$$
Thus the moments are:
$$M'(t) =\theta e^t |_{t=0} = \theta = M''(t) ... $$
Now to calculate the variance of the variance of $\theta$ i did the following:
$$var(var(\theta)) = var(E[\theta - E[\theta]]^2) = var(E[\theta^2] - (E[\theta])^2)$$
Now from the MGF we know that: $E[\theta] = \theta = E[\theta^2] = E[\theta^3] = E[\theta^4] $, thus:
$$var(E[\theta^2] - (E[\theta])^2) = var(\theta - \theta^2) = E[(\theta - \theta^2)^2] - (E[\theta - \theta^2])^2 = $$ $$ = E[\theta^2 - 2\theta^3 + \theta^4] - 0^2 = \theta - 2 \theta + \theta = 0$$.
Is there an error in my argument, or is this the correct solution.
In your example $\theta$ is not a Bernoulli random variable, but the parameter of the Bernoulli distribution, i.e., a number between 0 and 1. Hence, everything from "Now to calculate the variance..." onwards makes no sense. Also note $\mathbb{E}(\theta^n)=\theta^n$ for $n=1,2,\ldots$, $\operatorname{Var}(\theta)=0$, and $\operatorname{Var}(\operatorname{Var}(\theta))=\operatorname{Var}(0)=0$.