Unbiased estimators for the moments of 2 non-independent random variables Let $X$ and $Y$ be two non independent random variables. Suppose to generate $n$ realizations of both variables, and indicate with $(X_i, Y_i)$ their values. Also, let's pose that $X$ and $Y$ are non independent only within a particular realization; in other words, any realization does not affect the result of the others.
What can I do to find a statistical unbiased estimator of central moments
$$m_{a,b} = \mathbb{E}[(X-\mu_x)^a(Y-\mu_y)^b]$$
which uses the $n$ realizations $(X_i, Y_i)$ except from doing by hand all calculations for each $a$ and $b$? Is there some readings I'm missing?
Some additions
First of all, I forgot to write these:
$$\mu_X = \mathbb{E}[X], \mu_Y = \mathbb{E}[Y]$$
 A: You seek the bivariate h-statistic $h_{r, t}$ defined by:
$$E[h_{r, t}] = \mu _{r,t} $$
That is, $h_{r, t}$ is the statistic whose expectation is the bivariate central moment:
$$\mu _{r,t} = E\left[(X-E[X])^r (Y-E[Y])^t\right]$$
Unfortunately, deriving h-statistics by hand is extremely complicated and prone to error, and they are best obtained by computer algorithm, such as the HStatistic function in the mathStatica package (I should add I am an author of the latter). In fact, much of the published literature in this field contain errors, in some cases remaining uncorrected for 50 or more years.
To illustrate, an unbiased estimator of the bivariate central moment $\mu _{3,1}$ is the bivariate h-statistic $h_{3,1}$:

where the solution is expressed in terms of bivariate power sums:
$$s_{a, b}=\sum _{i=1}^n X_i^a Y_i^b $$
which you can then compute given your $(X_i, Y_i)$ data. All done.

To verify that pp above is, in fact, an unbiased estimator, compute $E[$[pp], which is just the 1st RawMoment of pp (and express the solution in terms of CentralMoments):


For more detail
For more detail, see, for instance, Chapter 7 of our book, "Mathematical Statistics with Mathematica". A free download of the chapter is now available here:
http://www.mathStatica.com/book/Rose_and_Smith_2002edition_Chapter7.pdf
Subsequent editions extend from univariate h-statistics to the multivariate case.

OP wrote:  

By the way, I need to have unbiased estimator for $h_{1,2}$, $h_{1,1}$, $h_{2,1}$ ...

No problem ...

And @Ivan  seeks $h_{2,2}$:

