Let $(X_i)^n_{i=1}$ be independent random variables from the same distribution (with CDF and PDF denoted as F and f). The kth order statistic $X_{(k)}$ is $Beta(k,n-k-1)$ - the PDF is:
$$f_{X_{(k)}}(t) = \frac{n!}{(n-k)!(k-1)!}[F(t)]^{k-1}[1-F(t)]^{n-k}f(t)$$
which is easy to prove.
My task demands to use, for $1 \leq k < m \leq n$, the joint distribution $(X_{(k)}, X_{(m)})$. Why the formula for it is:
$$f_{(X_{(k)}, X_{(m)})}(x,y) = C\cdot F^{k-1}(x)f(x)[F(y)-F(x)]^{m-k-1}f(y)[1-F(y)]^{n-m}\mathbb{1}_{\{x \leq y\}}(x,y),$$
where $C = \frac{n!}{(k-1)!(m-k-1)!(n-m)!}$?
I have found the formula for this (many sources, ie. en.wikipedia, expectation-of-product-of-two-order-statistics) but not how to get it. When I try to compute it, I receive a completely different formula...
My approach was as for the $f_{X_{(k)}}(t)$: to compute the CDF, then differentiate it, and reduce it. But the series is quite complex: $$F_{X_{(k)},X_{(m)}}(x,y) = \sum^{m-1}_{a=k}\sum^n_{b=m} {{m-1}\choose{k}}{n\choose m}[F(x)]^a[1-F(x)]^{m-1-a}[F(y)]^b[1-F(y)]^{n-b}$$
Is there a simpler way? Maybe directly from the PDF function $f_{X_{(k)}}$?
Here is the full proof.
To define the joint distribution of $(X_{(k)}, X_{(m)})$, for $1\leq k< m\leq n$, consider first the following picture of the event $(x_k < X_{(k)} \leq x_k + \delta x_k,x_m < X_{(m)} \leq x_m + \delta x_m)$:
$X_i\leq x_k$ for $k-1$ of the $X_i$'s, $x_k< X_i\leq x_k + \delta x_k$ for exactly one of the $X_i$'s, $x_k+\delta x_k < X_i \leq x_m$ for $m-k-1$ of the $X_i$'s, $x_m < X_i \leq x_m+\delta x_m$ for exactly one of the $X_i$'s, and $X_i> x_m + \delta x_m$ for the remaining $n-m$ of the $X_i$'s. Considering $\delta x_k$ and $\delta x_m$ to be both small, we may write
\begin{align*} P(x_k< &X_{(k)}\leq x_k+\delta x_k, x_m< X_{(m)}\leq x_m+\delta x_m)\\ & = \frac{n!}{(k-1)!(m-k-1)!(n-m)!}F(x_k)^{k-1}\times\\ &\quad\left\{F(x_m) - F(x_k+\delta x_k)\right\}^{m-k-1}\left\{1 - F(x_m+\delta x_m)\right\}^{n-k}\times\\ &\quad\left\{F(x_k+\delta x_k) - F(x_k)\right\}\left\{F(x_m+\delta x_m) - F(x_m)\right\}+\\ &\quad O((\delta x_k)^2\delta x_m) + O(\delta x_k(\delta x_m)^2), \end{align*}
where $ O((\delta x_k)^2\delta x_m)$ , $O(\delta x_k(\delta x_m)^2)$ are higher order terms which correspond to the probabilities of the event of having more than one $X_i$ in the interval $(x_k,x_k+\delta x_k]$ and at least one $X_i$ in the interval $(x_m,x_m+\delta x_m]$, and of the event of having one $X_i$ in $(x_k,x_k+\delta x_k]$ and more than one $X_i$ in $(x_m,x_m+\delta x_m]$, respectively.
From the above equation, we can derive the joint density function of $X_{(k)}, X_{(m)}$ to be
\begin{align*} f_{X_{(k)},X_{(m)}}&(x_k,x_m) = \lim_{\delta x_k\to0,\delta x_m\to 0}\left\{\frac{P(x_k< X_{(k)}\leq x_k+\delta x_k, x_m< X_{(m)}\leq x_m+\delta x_m)}{\delta x_k\delta x_m}\right\}\\ &=\frac{n!}{(k-1)!(m-k-1)!(n-m)!}F(x_k)^{k-1} F(x_k)^{k-1}\{F(x_m)-F(x_k)\}^{m-k-1}\times\\ &=\{1-F(x_m)\}^{n-m}f(x_k)f(x_m),\quad -\infty< x_k < x_m<\infty, \end{align*} Q.E.D.