Covariance of order statistics (uniform case) Let $X_1, \ldots, X_n$ be uniformly distributed on $[0,1]$ and $X_{(1)}, ..., X_{(n)}$ the corresponding order statistic. I want to calculate $Cov(X_{(j)}, X_{(k)})$ for $j, k \in \{1, \ldots, n\}$.
The problem is of course to calculate $\mathbb{E}[X_{(j)}X_{(k)}]$.
The joint density of $X_{(j)}$ and $X_{(k)}$ is given by $$f_{X_{(j)}, X_{(k)}}=\binom{n}{k}\binom{k}{j-1}x^{j-1}(y-x)^{k-1-j}(1-y)^{n-k}$$ where $0\leq x\leq y\leq 1$. (I used the general formula here.)
Sadly, I see no other way to calculate $\mathbb{E}[X_{(j)}X_{(k)}]$ than by $$\mathbb{E}[X_{(j)}X_{(k)}]=\binom{n}{k}\binom{k}{j-1}\int_0^1\int_0^yxyx^{j-1}(y-x)^{k-1-j}(1-y)^{n-k}\,dx\,dy.$$
But this integral is too much for me. I tried integration by parts, but got lost along the way.
Is there a trick to do it? Did I even get the limits of integration right?
Apart from that, I wonder if there's a smart approach to solve the whole problem more elegantly.

Thank you so much for posting this -- I too looked at this covariance integral (presented in David and Nagaraja's 2003 3rd edition Order Statistics text) and thought that it looked ugly. However, I fear that there may be a few small mistakes in your math on E(X_jX_k), assuming that I'm following you right. The joint density should have (j-1)! in the denominator instead of j! at the outset -- otherwise j! would entirely cancel out the j! in the numerator of Beta(j+1,k-j) instead of ending up with j in the numerator of the solution... Right? 
