# Does the (normalized) product of two independent binomial variables converges in distribution to a normal variable?

(I asked this question on MSE 10 days ago, but I got no answer.)

Let $$X$$ and $$Y$$ be two independent identically distributed binomial random variables with parameters $$n \in \mathbb{N}$$ and $$p \in (0,1)$$. Let $$Z := XY$$ be their product.

Is it true or false that $$\tilde{Z} := (Z - \mathbf{E}[Z]) / \sqrt{\mathbf{Var}[Z]}$$ converges in distribution to a standard normal random variable (as $$n \to \infty$$) ?

At a first glance, I would be tempted to write $$X = \sum_{i=1}^n A_i$$ and $$Y = \sum_{i=1}^n B_i$$, where $$A_i$$ and $$B_i$$ are Bernoulli random variables, and then to apply the central limit theorem to $$Z = \sum_{i=1}^n \sum_{j = 1} A_i B_j$$... but $$A_i B_j$$ are not independent...

Thanks for any help

P.S.1 For $$p=1/2$$, it is easy to check that $$\mathbf{E}[Z] = n^2 / 4$$ and $$\mathbf{Var}[Z] = n^3 / 8 + n^2 / 16$$. Moreover, expanding $$(XY - n^2/4)^k$$ with the binomial theorem and using the formula for the moments of the binomial distribution, I got that

$$\mathbf{E}\left[\tilde{Z}^k\right] = \frac1{(n^3 / 8 + n^2 / 16)^{k/2}} \sum_{j=0}^k \binom{k}{j} \left(\sum_{i=0}^j {j \brace i} (n)_{i} (1/2)^i\right)^2 (-n^2/4)^{k-j}$$

where $${j \brace i}$$ are Stirling numbers of second kind and $$(n)_{i}$$ is a falling factorial. With this formula, I verified that the first 20 moments of $$\tilde{Z}$$ tend to the moments of a standard normal variable. However, I still do not know how to prove this for all moments.

P.S.2 I think that one cannot prove the claim only using the fact that the (normalized) $$X$$ and $$Y$$ converge in distribution to normal variables. In fact, it is known that the product of two independent normal variables is not a normal variable.

• Your example is for $p=1/2$, right? Sep 10 at 15:36
• @BrendanMcKay Yes, it for $p=1/2$. Sorry I forgot that Sep 10 at 18:02

$$\newcommand{\R}{\mathbb R}\newcommand\ep\epsilon\newcommand\tsi{\tilde\sigma}$$Yes, of course. This follows by the multivariate (here, bivariate) so-called delta method.
Indeed, we may assume that $$\begin{equation*} X=\sum_{i=1}^n X_i,\quad Y=\sum_{i=1}^n Y_i, \end{equation*}$$ where $$X_1,Y_1,\dots,X_n,Y_n$$ are independent identically distributed (iid) Bernoulli random variables (r.v.'s) with parameter $$p\in(0,1)$$.
For each $$i\in[n]:=\{1,\dots,n\}$$, let $$\begin{equation*} V_i:=(X_i-p,Y_i-p), \end{equation*}$$ so that $$V_1,\dots,V_n$$ are iid zero-mean random vectors in $$\R^2$$. Then $$\begin{equation*} XY=n^2 f(\bar V)+n^2p^2, \tag{1}\label{1} \end{equation*}$$ where $$\bar V:=\frac1n\,\sum_{i=1}^n V_i$$ and for $$(x,y)\in\R^2$$ $$\begin{equation*} f((x,y)):=f(x,y):=(x+p)(y+p)-p^2, \end{equation*}$$ so that $$f(0,0)=0$$ and for $$L:=f'(0,0)$$ we have $$\begin{equation*} L(x,y)=px+py \end{equation*}$$ and $$\begin{equation*} |f(x,y)-L(x,y)|=|xy|\le\tfrac12\,\|(x,y)\|^2, \end{equation*}$$ where $$\|(x,y)\|:=\sqrt{x^2+y^2}$$. So, condition (2.1) of this paper holds (for any real $$\ep>0$$ and $$M_\ep=1$$).
Moreover, $$v_3:=E\|V_1\|^3<\infty$$ and $$\begin{equation*} \tsi:=\sqrt{EL(V_1)^2}=\sqrt{p^2 E(X_1-p+Y_1-p)^2}=\sqrt{2p^3q}>0, \end{equation*}$$ where $$q:=1-p$$. So, by \eqref{1} and Theorem 2.9 of the same paper, $$\begin{equation*} \frac{XY-n^2p^2}{n^2\sqrt{2p^3q/n}} =\frac{f(\bar V)}{\tsi/\sqrt n}\to Z\sim N(0,1) \tag{2}\label{2} \end{equation*}$$ (in distribution as $$n\to\infty$$).
Note also that $$EXY=n^2p^2$$ and $$Var\,XY=2n^3p^3q+n^2p^2q^2=2n^3p^3q(1+O(1/n))\\ \sim2n^3p^3q=(n^2\sqrt{2p^3q/n})^2.$$ Thus, by \eqref{2}, $$\begin{equation*} \frac{XY-EXY}{\sqrt{Var\,XY}}\to Z, \end{equation*}$$ as desired.
Moreover, it follows from cited Theorem 2.9 that $$\begin{equation*} \sup_{z\in\R}\Big|P\Big(\frac{XY-EXY}{\sqrt{Var\,XY}}\le z\Big)-P(Z\le z)\Big|\le \frac{C_p}{\sqrt n} \end{equation*}$$ for some real $$C_p>0$$ depending only on $$p\in(0,1)$$ and all natural $$n$$.