Expected correlation of sample correlations It's stated in a number of references online (eg 1) that the sampling distribution of Pearson's r is approximately $\mathcal N(\rho, \frac{1}{\sqrt{n-3}})$, with Fisher's transformation indicated for large absolute correlations (say greater than 0.4).
I'm taking that as given, but what I'm trying to determine is as follows: taking $r_{x,z}$ and $r_{y,z}$ as respectively the sample estimates of $\rho_{x,z} = \text{Cor}(X,Z)$ and $\rho_{y,z} = \text{Cor}(Y,Z)$, and repeating that for many samples, what is the expected value of $\text{Cor}(r_{x,z},r_{y,z})$? It seems natural that it be equal to $\rho_{x,y} = \text{Cor}(X,Y)$, particularly when considering the extrema where $\rho_{x,y}$ is equal to 1, 0, or -1. This seems to be borne out in my empirical tests, at least where $(X, Y, Z)$ have a multivariate normal distribution (see code below).
I've spent some time trying to prove this but have been unable to do so, so would appreciate any pointers people might have. I'm not entirely convinced it is indeed correct, and even if so, what assumptions (if any) are implicitly relied upon.
Below is the R code I used for doing empirical testing (which I've tested on many parameter variations):
library(MASS)

pxz =  0.10   # population correlation of X and Z
pyz =  0.05   # population correlation of Y and Z
pxy = -0.50   # population correlation of X and Y
N = 100000    # population size
n = 100       # sample size
m = 1000      # number of samples to draw

# generate population data
data = mvrnorm(N, mu=c(0, 0, 0), Sigma=matrix(c(1, pxy, pxz, pxy, 1, pyz, pxz, pyz, 1), 3));
X = data[, 1]
Y = data[, 2]
Z = data[, 3]

# verify the sampling distribution of rxz and ryz
r = t(sapply(1:m, function(., i=sample(N, n)) c(cor(X[i], Z[i]), cor(Y[i], Z[i]))))
tanh(apply(atanh(r), 2, mean)) # ~= c(pxz, pyz)
apply(atanh(r), 2, sd) # ~= 1/sqrt(n-3) * sqrt(N-n)/sqrt(N-1)

# calculate the correlation of m sample correlations 100 times
vals = replicate(100,cor(t(sapply(1:m, function(., i=sample(N, n)) c(cor(X[i], Z[i]), cor(Y[i], Z[i])))))[2])
mean(vals) # ~= pxy

 A: Let's start out with a small correction. The large-sample distribution of the raw correlation coefficient is $$r \sim N\left(\rho, \frac{(1 - \rho^2)^2}{n - 1}\right)$$ (or often just $n$ in the denominator). If we let $z_r$ denote the Fisher's r-to-z transformed correlation coefficient, then the large-sample distribution is $$z_r \sim N\left(z_\rho, \frac{1}{n-3}\right),$$ where $z_\rho$ is the r-to-z transformed true correlation. Note that this is also a large-sample approximation, but this one works very well even for relatively small sample sizes. Also, I am giving the variance in the equations, not the standard deviations.
As for the correlation between correlations, I think the earliest derivation can be found in Pearson and Filon (1989). A more recent summary can be found in Steiger (1980) and several other articles.
In particular, the covariance between two correlation coefficients sharing one variable is $$Cov[r_{jk}, r_{jh}] = \frac{\rho_{kh} (1 - \rho_{jk}^2 - \rho_{jh}^2) - \frac{1}{2}(\rho_{jk}\rho_{jh}) (1 - \rho_{jk}^2 - \rho_{jh}^2 - \rho_{kh}^2)}{n - 1}$$ (again using $n-1$ in the denominator for consistency). So, if you want the correlation, just divide the covariance by the standard deviations, which yields $$Cor[r_{jk}, r_{jh}] = \frac{\rho_{kh} (1 - \rho_{jk}^2 - \rho_{jh}^2) - \frac{1}{2}(\rho_{jk}\rho_{jh}) (1 - \rho_{jk}^2 - \rho_{jh}^2 - \rho_{kh}^2)}{(1-\rho_{jk}^2)(1-\rho_{jh}^2)}.$$
References
Pearson, K., & Filon, L. N. G. (1989). Mathematical contributions to the theory of evolution: IV. On the probable errors of frequency constants and on the influence of random selection on variation and correlation. Philosophical Transactions of the Royal Society of London, Series A, 191, 229-311.
Steiger, J. H. (1980). Tests for comparing elements of a correlation matrix. Psychological Bulletin, 87, 245-251.
