Consider the following expression:
$$ D = \left[\sum_{m=1}^M \left(\frac{1}{\sqrt{R}} \sum_{r=1}^R x^{p_1}_{R(m-1)+r}\right)^2\right] \left[\sum_{m=1}^M \left(\frac{1}{\sqrt{R}} \sum_{r=1}^R x^{p_2}_{R(m-1)+r}\right)^2\right] - \left[\sum_{m=1}^M \left(\frac{1}{\sqrt{R}} \sum_{r=1}^R x^{p_1}_{R(m-1)+r}\right)\left(\frac{1}{\sqrt{R}} \sum_{r=1}^R x^{p_2}_{R(m-1)+r}\right) \right]^2 $$
where all $x$ are i.i.d random variables (zero-mean Gaussian for simplicity). I am seeking $\mathbb{E}[ D^q ]$ (e.g. $q=1$ and $q=2$).
The expectation operator $\mathbb{E}$ is linear and can be applied by expanding the entire sum, grouping the same $x$ together and replacing them with their moments.
To understand better what I mean consider the much simplified expression with $R=1$:
$$ D = \left(\sum_{m=1}^M x^{2 p_1}_{m} \right) \left(\sum_{m=1}^M x^{2p_2}_{m}\right) - \left(\sum_{m=1}^M x^{p_1+p_2}_{m} \right)^2 = \sum_{m_1=1}^M \sum_{m_2=1}^M \left( x^{2 p_1}_{m_1}x^{2 p_2}_{m_2} - x^{p_1+p_2}_{m_1}x^{p_1+p_2}_{m_2} \right) $$
The two sums can be easily analyzed for the 2 cases $m_1=m_2=m$ and $m_1\neq m_2$ to yield:
$$ \mathbb{E}(D)=M(M-1) [ \mu(2p_1)\mu(2p_2) - \mu^2(p_1+p_2) ] $$
This approach generalizes to higher order sums using symmetric functions. However, for $R>1$ the expression is not symmetric any more.
Another approach I follwed was to find a pattern in the expanded sum (e.g., using Mathematica). However, this becomes infeasible. For example, $\mathbb{E}[D]$ expands to a 6-fold sum and $\mathbb{E}[D^2]$ to a 12-fold sum. If $M=R=4$ the expanded sum has already $4^{12}=16777216$ terms - even through they collapse to just a handfull terms when replacing them with their moments. However, for finding a patterns I would need to sweep $M$ and $R$ at least between one and 10 or so.
Any thoughts on how to tackle this sum?