I recently tackled a problem and I arrived at something of the following form,
$$ \frac{1}{n^2} \sum_{i=0}^{n-1}\sum_{j=0}^{n-1} f(T^ix, T^jy), $$
where $T$ is a measure preserving transformation and I am interested in the limit as $n$ tends to infinity. In my case it is the shift operator in a probability space.
In the uni-variate case Birkhoff's ergodic theorem states that for a measure preserving transformation $T$,in a measurable space $(X, \mathscr{B})$, with a $\sigma$-finite measure $\mu$ $$\frac{1}{n} \sum_{i=0}^{n-1} f(T^ix)$$ converges a.e. to a function $f*\in L^1$, with $f^* = \frac{1}{\mu(X)}\int fd\mu$. Is there an equivalent result for the multivariate case? Will the first equation converge to it's average?
Rephrased in an other way, the question is about the convergence of $$ \frac 1{n^2}\sum_{i=1}^n\sum_{j=1}^nf\left(X_i,X_j\right), $$ where $\left(X_i\right)_{i\geqslant 1}$ is a strictly stationary sequence. In order to hope something, we have to assume that $f\left(X_i,X_j\right)$ is integrable for all $i,j$. The diagonal part $\sum_{i=1}^nf\left(X_i,X_i\right)$ has, by the classical ergodic theorem, a negligible contribution hence we are reduced to study the asymptotic behavior of $$ \frac 1{n^2}\sum_{1\leqslant i<j\leqslant n}f\left(X_i,X_j\right). $$ (the part with $j>i$ follows by using $g\colon (u,v)\mapsto f(v,u)$). The keyword is then $U$-statistics. A version of the ergodic theorem could be $$ \frac 1{n^2}\sum_{1\leqslant i<j\leqslant n}\left(f\left(X_i,X_j\right)-\mathbb E\left[f\left(X_i,X_j\right)\right]\right)\to 0 \mbox{ a.s.} $$ but the question seems to be open.
There are some results when the product of marginal distribution functions is continuous (see
Borovkova, S.; Burton, R.; Dehling, H. Consistency of the Takens estimator for the correlation dimension. Ann. Appl. Probab. 9 (1999), no. 2, 376--390. doi:10.1214/aoap/1029962747. https://projecteuclid.org/euclid.aoap/1029962747).