First and second moments of deep nesting of the Binomial-Binomial hierarchical model? I am interested in the Binomial-Binomial hierarchical model, where the number of trials itself follows a binomial distribution. I would like to know the expected value (first central moment, $\mu_1$) and variance (second central moment, $\mu_2$) of this model (if not possible to derive the Probability Mass Function), and, whether it is possible to derive those quantities for any arbitrary $i^{th}$ nesting of these models.

Question: Is there a closed-form expression for $\mu_{1}(X_{i})$ and $\mu_2(X_{i})$?
 A: It's actually pretty straightforward.  First, the expectation of an arbitrarily deep nesting of binomial distributions with probability parameters $p_i, i \in \{1,\cdots,K\}$ and "top" number of trials parameter $n$ is just $n \Pi_{i=1}^Kp_i$.  This should be intuitively clear, as expectation is a linear operator.  For the variance, life is slightly more complicated, so I'll walk you through a simple characteristic function analysis that works for distributions of sums of i.i.d. random variables, where the sum is of $n$ terms and $n$ follows some other distribution.  
A brief review of characteristic functions:  The ch.f. of a distribution $F$ is defined as $\phi(it) = \int e^{itx}dF(x)$, where $i = \sqrt{-1}$.  By expanding $e^{itx}$ inside the integral, which I won't do, we can see that taking the derivatives of $\phi(it)$ will give us the noncentral moments of the distribution: $\mathbb{E}[X^k]= (-i)^k \phi^{(k)}(0)$, where $\phi^{(k)}(0)$ is the $k^{th}$ derivative of $\phi$ with respect to $t$ evaluated at $t=0$.  We can recover the central moments by using the known relationships between them and the noncentral moments.
Let $x$ follow a distribution with ch.f. $\phi(it)$.  In your case $x$ is Bernoulli with parameter $p$. The sum of $n$ i.i.d. $x$s, label it $y$, has ch.f. $\phi^n(it)$, and is (obviously) conditional on $n$.   If $n$ follows a distribution $f$ with ch.f. $\theta(it)$, we can find the ch.f. of $y$ unconditional on $n$, label that $\eta(it)$, from:
$\eta(it) = \sum_{n=0}^\infty \phi^n(it) f(n) = \sum_{n=0}^\infty \exp\{n \log(\phi(it))\}f(n) = \theta(\log(\phi(it)))$
We just substitute $\log(\phi(it))$ into $\theta(it)$ wherever we see $it$ in the latter.  Now we can take the derivatives of $\eta(it)$ with respect to $t$ and find the moments (after the second, it becomes a real pain) in terms of the moments of $n$ and $x$.
I won't do that explicitly here, it's good practice if you want to work on your calculus, but, jumping to the answer:
$\sigma^2_y = \text{E}_n \sigma^2_x + \text{E}^2_x \sigma^2_n$
This formula pops up in inventory control when calculating safety stocks, as it is the variance of leadtime demand where per-period demand is i.i.d. and leadtime is random.
Obviously you can calculate this iteratively for any number of nested binomial distributions, starting with the penultimate one and working your way down the chain.
