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.
I don't know if there is a smart way to solve the problem, but the standard way is the following:
For $1\le j< k\le n$, the joint density function should be $$ f_{X_{j},X_{_k}}=\frac{n!}{(j-1)!\,(k-j-1)!\,(n-k)!}\,x^{j-1}\,(y-x)^{k-1-j}\,(1-y)^{n-k} $$ where $0<x<y<1.$
Now to compute $E[X_{j}X_{_k}]=\int_0^1\mathbb{d}y\int_0^y xy\cdot f(x,y)\,\mathbb{d}x$, you need to be familiar with Beta function. Let $$C=\frac{n!}{(j-1)!\,(k-j-1)!\,(n-k)!},$$ so that $$\begin{eqnarray} E[X_{j}X_{_k}] &=&C\int_0^1\mathbb{d}y\int_0^y xy\cdot x^{j-1}(y-x)^{k-1-j}(1-y)^{n-k}\,\mathbb{d}x \\ &=&C\int_0^1\mathbb{d}y\int_0^y \left[\frac{x}{y}\right]^{j}\left[1-\frac{x}{y}\right]^{k-1-j}y^k(1-y^{n-k})\,\mathbb{d}x\\ &=&C\cdot B(j+1,k-j)\int_0^1y^{1+k}(1-y)^{n-k}\,\mathbb{d}y \\ &=&C \cdot B(j+1,k-j)\cdot B(k+2,n-k+1)\\ &=& \frac{j\cdot (k+1)}{(n+1)(n+2)}. \end{eqnarray}$$
Finally, $$ \begin{eqnarray}\operatorname{Cov}[X_{j}X_{_k}]&=&E[X_{j}X_{_k}]-E[X_{j}]\cdot E[X_{_k}]\\ &=& \frac{j\cdot (k+1)}{(n+1)(n+2)}-\frac{j\cdot k}{(n+1)^2} \\ &=& \frac{j\cdot (n+1-k)}{(n+1)^2(n+2)}. \end{eqnarray}$$