This answer takes seriously the assumption that the data-generating distribution is known. Whether a closed-form analytical solution exists depends on the distribution, but it's possible to write down a universal formula for continuous distributions that permits numerical calculation of the variances of the quantile means.
Let's begin by formulating the problem. A random sample of a continuous distribution $F$ (with density $f$) is collected. When sorted, the values may be written as order statistics
$$x_{(1)} \le x_{(2)} \le \cdots \le x_{(n)}.$$
The background for this question concerns partitioning these data into "quantiles." For instance, we might take the first fifth of the order statistics from $x_{(1)}$ through $x_{(n/5)}$ (rounding $n/5$ suitably), the next fifth, and so on, creating five subsets; and characterize each such subset in terms of its mean.
Generally, given two indexes $1\le i \le j \le n,$ define the subsample
$$X_{i;j} = (x_{(i)}, x_{(i+1)}, \ldots, x_{(j)})$$
and let
$$\bar x_{i;j} = \frac{1}{j-i+1} (x_{(i)} + x_{(i+1)} + \cdots + x_{(j)})$$
be its mean. These are statistics--functions of the sample--and as such have a sampling distribution. The question I will address concerns finding the first and second moments of $(\bar x_{i;j}, \bar x_{i^\prime;j^\prime}).$ These could be used to compute the variance-covariance matrix of the quantile means, for instance.
The analysis is straightforward upon observing that the distribution of $x_{(i)}$ has density
$$f_i(x) = \binom{n}{i-1,1,n-i} F(x)^{i-1}(1-F(x))^{n-i} f(x)$$
and the joint distribution of $(x_{(i)}, x_{(j)})$ has density
$$\eqalign{&f_{i,j}(x,y) = \\&\binom{n}{i-1,1,j-i-1,1,n-j} F(x)^{i-1} (F(y)-F(x))^{j-i-1} (1-F(y))^{n-j} f(x) f(y) \mathcal{I}(x\le y).}$$
(See https://stats.stackexchange.com/a/4684/919 for an intuitive explanation of how such a formula can be derived and remembered.) The normalizing constants are multinomial coefficients and can be evaluated as
$$\binom{n}{a,b,\ldots, z} = \exp\left(\log \Gamma(n+1)-\log\Gamma(a+1)-\log\Gamma(b+1)-\cdots - \log\Gamma(z+1)\right).$$
The definitions of moments imply
$$\mu_{i;k} = E[x_{(i)}^k] = \int_\mathbb{R} x^k f_i(x) \mathrm{d}x$$
for the univariate moments of order $k$ and
$$\mu_{i,j;k,l} = E[x_{(i)}^k x_{(j)}^l] = \iint_{\mathbb{R}^2} x^k y^l f_{i,j}(x,y) \mathrm{d}x \mathrm{d}y$$
for the bivariate moments of order $(k,l).$ For $j\gt i$ the covariance is
$$\operatorname{Cov}(x_{(i)}, x_{(j)}) = \mu_{i,j;1,1} - \mu_{i;1}\mu_{j;1}$$
and for $j=i$ it is
$$\operatorname{Cov}(x_{(i)}, x_{(i)})=\operatorname{Var}(x_{(i)}) = \mu_{i,2} - \mu_{i;1}^2.$$
In some special cases--notably, when $F$ is uniform--these integrals can be evaluated, but usually they need numerical integration.
Let, then, $\Sigma_F$ denote the variance-covariance matrix of all the order statistics,
$$\Sigma_F(i,j) = \operatorname{Cov}(x_{(i)}, x_{(j)}).$$
The quantile means $\bar x_{i,j}$ are linear combinations of the order statistics with coefficients $1/(j-i+1)$ for $x_{(k)}$ when $i\le k \le j$ and zero otherwise. Let $\xi_{i,j}$ be this vector of coefficients:
$$\xi_{i,j}(k) = \left\{\matrix{\frac{1}{j-i+1}& i\le k \le j \\ 0& \text{otherwise.}}\right.$$
Because covariance is a quadratic form,
$$\operatorname{Cov}(\bar x_{i,j}, \bar x_{i^\prime, j^\prime}) = \xi_{i,j}^t\,\Sigma_F\,\xi_{i^\prime, j^\prime}.$$
The R code below carries out this analysis for any distribution for which you can compute the values of its CDF $F$ and PDF $f.$ It illustrates it with the standard Normal distribution, then checks its results with a simulation of 10,000 independent samples. The output of that check is a histogram of the relative errors of the $5\times 5 = 25$ entries in the variance-covariance matrix of the five quintiles. That they are all small--not much larger than $1/\sqrt{10000}=0.01$ in size--attests to the correctness of the results.
For large values of $n,$ the computational demands will grow too large. Further analysis of the integrals for any particular $F$ will help develop approximations or asymptotic formulas. Saddle-point approximations may be particularly attractive.
#
# Log multinomial coefficients.
# The arguments are the *bottom* numbers: they must sum to the top value `n`.
#
lmultinom <- function(...) {
i <- unlist(list(...))
lgamma(sum(i)+1) - sum(lgamma(i+1))
}
#
# Univariate moments of order statistics.
# `lpf` returns the log of the CDF. It must recognize a "lower.tail=FALSE"
# argument to return the log of the complementary CDF.
# `ldf` returns the log of the PDF.
# `k` is the order of the moment.
#
e1 <- function(i, n, lpf, ldf, lower, upper, k=1, ...) {
lcnst <- lmultinom(i-1, 1, n-i)
f <- function(x) {
x^k * exp(lcnst + (i-1)*lpf(x) + (n-i)*lpf(x, lower.tail=FALSE) + ldf(x))
}
integrate(f, lower, upper, ...)
}
#
# Bivariate moments of order statistcs.
# `k` and `l` are the orders of the moments.
#
e2 <- function(i, j, n, lpf, ldf, lower, upper, k=1, l=1, ...) {
lcnst <- lmultinom(i-1, 1, j-i-1, 1, n-j)
f <- function(x, y) {
x^k * y^l * exp(lcnst + (i-1)*lpf(x) + (n-j)*lpf(y, lower.tail=FALSE) +
(j-i-1)*log(exp(lpf(y)) - exp(lpf(x))) + ldf(x) + ldf(y))
}
g <- Vectorize(function(y) {
integrate(f, lower, y, y=y, ...)$value
})
integrate(g, lower, upper, ...)
}
#
# EXAMPLE: The standard Normal distribution.
#
lpf <- function(x, ...) pnorm(x, log=TRUE, ...)
ldf <- function(x, ...) dnorm(x, log=TRUE, ...)
lower <- -8; upper <- 8 # Effective range
n <- 20 # The double integrals take a few seconds, but only have to be
# ..performed once and can be stored for efficiency.
#---Compute the means of the order statistics.
y <- sapply(1:n, function(i) e1(i, n, lpf, ldf, lower, upper)$value)
# plot(1:n, y,
# sub=bquote(n==.(n)),
# main="Expected Order Statistics\nStandard Normal Distribution",
# xlab="Order", ylab="Value")
#---Compute the second moments of the order statistics.
z <- matrix(0.0, n, n)
for (i in seq_len(n)) {
for (j in seq_len(n-i)+i) {
z[i,j] <- z[j,i] <- e2(i, j, n, lpf, ldf, lower, upper)$value
}
}
diag(z) <- sapply(1:n, function(i) e1(i, n, lpf, ldf, -8, 8, k=2)$value)
#---Obtain the variance-covariance matrix.
v <- z - outer(y,y)
#
# Compute the variance-covariance matrix for the quintile means.
#
quantiles <- round(n * 0:5/5)
Sigma <- matrix(NA_real_, length(quantiles)-1, length(quantiles)-1)
for (a in seq_len(dim(Sigma)[1])) {
i <- quantiles[a]+1
j <- quantiles[a+1]
x <- rep(1/(j-i+1), j-i+1)
for (b in seq_len(dim(Sigma)[2])) {
i0 <- quantiles[b]+1
j0 <- quantiles[b+1]
x0 <- rep(1/(j0-i0+1), j0-i0+1)
Sigma[a,b] <- Sigma[b,a] <- crossprod(x, v[i:j, i0:j0] %*% x0)
}
}
#------------------------------------------------------------------------------#
#
# Verify with a simulation.
# Each column of `x` is an ordered sample.
#
n.sim <- 1e4
set.seed(17)
x <- apply(matrix(rnorm(n*n.sim), n, n.sim), 2, sort)
#---Compute the quantile means of each sample.
xbar <- sapply(seq_len(dim(Sigma)[1]), function(a) {
i <- quantiles[a]+1
j <- quantiles[a+1]
colMeans(x[i:j, ])
})
#---Estimate the variance-covariance matrix from the simulated quantile means.
Sigma.hat <- cov(xbar)
#---Compare the estimate to the numerical result. We're looking for tiny
# values, on the order of 1/sqrt(n.sim).
#
hist(Sigma.hat / Sigma - 1, main="Relative errors",
xlab=expression(hat(Sigma) / Sigma - 1))
