Asymptotic distribution of statistics with proportions I am taking a graduate course in mathematical statistics. 
Here is what I feel must be a straightforward problem in application of standard tools in "asymptotics" (such as application of Slutsky's theorem, Cramer's theorem, and quadratic forms). 
I am struggling to factorize the problem properly and end up running into cumbersome forbidding expression I am not able to develop. There must be a neat way to solve this problem.  
I have little practice solving problems in statistics, any help would be very much welcome!       


Suppose we have $n$ iid samples of a bivariate Bernoulli random variable $(x, y)$ where $x$ and $y$ take values 0 or 1 and are independent. 
Define $\hat{p}_{ij}$ to be the proportion of the samples such that $(x= i, y = j)$. Also define $\hat{p}_{i*} := \hat{p}_{i0} + \hat{p}_{i1}$ (which is identified as the proportion of samples where $x = i$) and similarly $\hat{p}_{*j} := \hat{p}_{0j} + \hat{p}_{1j}$. 
What is the asymptotic distribution of the following statistics: 
$$
S(n) = \sum_{i=0}^{1} \sum_{j=0}^{1}\frac{n(\hat{p}_{ij} - \hat{p}_{i*}\hat{p}_{*j})^2}{\hat{p}_{i*}\hat{p}_{j*}}\,.
$$

Added later: 
Following some comments, here is an elaboration on the definition of the quantities involved:
Let $n$ iid samples of the bivariate Bernoulli $\{(x, y)_k\}_{k=1}^n$, define $$N_{i,j} := \{\text{number of indices } k ~| ~(x= i, y = j)_k\}$$ with this definition at hand the proportion $\hat{p}_{ij}$ is defined to be 
$$
\hat{p}_{ij} := \frac{N_{ij}}{n}\,.
$$
and the asymptotic of $S(n)$ is with respect to the limit $n \to \infty$. 
 A: This will not be an answer, only the explanation of how I would approach this problem. Asymptotics is usually always about the sums, so it is better to express all the quantities in term of sums. The main sum has 4 summands, where each summand depends on $n$, so it is clear that we need to deal with each summand independently. 
Assume that we have a sample $Z_i$, $i=1,...,n$ where $Z_i=(X_i,Y_i)$ and $X_i$ and $Y_i$ are independent Bernoulli random variables. Then we have
$$\hat{p}_{1*}=\frac{\sum_{i=1}^nX_i}{n},\;\; \hat{p}_{0*}=\frac{\sum_{i=1}^n(1-X_i)}{n}.$$
The definitions for $\hat{p}_{*1}$ and $\hat{p}_{*0}$ are respectively the same. 
Now 
$$\hat{p}_{11}=\frac{\sum_{i=1}^nX_iY_i}{n},\;\; \hat{p}_{00}=\frac{\sum_{i=1}^n (1-X_i)(1-Y_i)}{n},$$
and 
$$\hat{p}_{01}=\frac{\sum_{i=1}^n(1-X_i)Y_i}{n},\;\; \hat{p}_{10}=\frac{\sum_{i=1}^n X_i(1-Y_i)}{n}.$$
By substituting these expressions into the initial sum and doing some algebra, the resulting expressions either will look as the selfnormalizing sums of independent random variables, which will give you asymptotic results immediately, or you will get some form of $\chi^2$ statistic for which again you can find some sort of asymptotic results. These are of course some sort of educated guesses.
