Is the sample correlation always positively correlated with the sample variance?

The sample correlation $$r$$ and the sample standard deviation of $$X$$ (call it $$s_X$$) seem to be positively correlated if I simulate bivariate normal $$X$$, $$Y$$ with a positive true correlation (and seem to be negatively correlated if the true correlation between $$X$$ and $$Y$$ is negative). I found this somewhat counterintuitive. Very heuristically, I suppose it reflects the fact that $$r$$ represents the expected increase in Y (in units of SD(Y)) for a one-SD increase in X, and if we estimate a larger $$s_X$$, then $$r$$ reflects the change in Y associated with a larger change in X.

However, I would like to know if $$Cov(r, s_x) >0$$ for $$r>0$$ holds in general (at least for the case in which X and Y are bivariate normal and with large n). Letting $$\sigma$$ denote a true SD, we have:

$$Cov(r, s_X) = E [ r s_X] - \rho \sigma_x$$

$$\approx E \Bigg[ \frac{\widehat{Cov}(X,Y)}{s_Y} \Bigg] - \frac{Cov(X,Y)}{\sigma_Y}$$

I tried using a Taylor expansion on the first term, but it depends on $$Cov(\widehat{Cov}(X,Y), s_Y)$$, so that’s a dead end. Any ideas?

EDIT

Maybe a better direction would be to try to show that $$Cov(\widehat{\beta}, s_X)=0$$, where $$\widehat{\beta}$$ is the OLS coefficient of Y on X. Then we could argue that since $$\widehat{\beta} = r \frac{s_Y}{s_X}$$, this implies the desired result. Since $$\widehat{\beta}$$ is almost like a difference of sample means, maybe we could get the former result using something like the known independence of the sample mean and variance for a normal RV?

• It would be unchanged. Hmm. I'm afraid I don't yet see the relevance, though. – half-pass Apr 1 at 19:45
• I should probably also note that while I wish this were a homework question, it's not... :) – half-pass Apr 1 at 19:46
• Ah, I didn't read the question carefully enough. My apologies. – jbowman Apr 1 at 19:50
• The first equality in your calculation is not correct. $s_x = \sqrt{s^2_x}$ is consistent for the standard deviation, but is not unbiased: en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation – Andrew M Apr 1 at 19:51
• It's extremely close to unbiased for large n, though -- the rule-of-thumb correction factor for a normal RV is (n - 1.5) vs. (n-1). – half-pass Apr 1 at 19:54

TL;dr

The off-diagonal entries of the sample covariance will generally be correlated with the diagonal entries because $$E(XY^3) - E(XY)E(Y^2) = 0$$ only when special conditions of the mixed 4th order moments hold. When $$(X,Y)$$ is bivariate Gaussian, these conditions hold only when $$X$$ is independent of $$Y$$.

Details

There is an asymptotic result that can be shown here by examining the limiting distribution of $$\sqrt n$$-times the sample covariance (by the CLT, it's going to be multivariate normal) and then applying the delta method. This unfortunately means we'll have to take detour through a derivation of the distribution of the sample covariance$$^1$$ since I can't find any good references to it online. Alternately, if you are willing to assume normality, then knowledge of the covariance of the Wishart distribution would let you skip directly to section 2.

1 The asymptotic distribution of the sample covariance

Let $$V_1, \dotsc, V_n$$ be an iid sample from a bivariate distribution $$V_i = \begin{pmatrix} X_i \\ Y_i \end{pmatrix}$$ with finite fourth moments, and let $$\text{Cov}(V_i) = \begin{pmatrix} \sigma^2 & \rho \sigma \tau \\ \rho \sigma \tau & \tau^2 \end{pmatrix} = \Sigma.$$ Without loss of generality and to avoid some annoying additional bookkeeping we'll assume $$E(V_i) = \mathbf{0}$$.

Then by the linearity of expectation and the weak law of large numbers, the sample covariance $$S_n = \frac{1}{n-1} \sum_{i=1}^n (V_i - \bar V_n) (V_i - \bar V_n)^T = \frac{1}{n-1}\sum_{i=1} V_i V_i^T - \frac{n}{n-1} \bar V_n \bar V_n^T$$ is unbiased and consistent for $$\Sigma$$, and in fact $$\sqrt{n} (S_n - \Sigma) \rightarrow_d N(0, \Lambda).$$

The exercise thus passes to determining $$\Lambda$$. For a symmetric matrix $$\mathbf{A} = \begin{pmatrix} a & b \\ b & c \end{pmatrix}$$, let $$\tilde{\mathbf{A}} = (a, b, c)^T$$ be the "vectorization" of its upper triangle. Now consider a single element of the average that enters into the leading term (the scatter matrix) of $$S_n$$: $$\tilde Z_i = \widetilde{V_i V_i^T} = \begin{pmatrix} X_i^2 \\ X_i Y_i \\ Y_i^2 \end{pmatrix}.$$ Clearly by the zero-mean assumption, already $$E(Z_i) = \tilde \Sigma$$ and by considering the powers of $$X$$ and $$Y$$ that appear in $$\tilde{Z}_i \tilde{Z}_i^T$$ we can just write $$\text{Cov}(\tilde Z_i) = E(\tilde Z_i \tilde Z_i^T) - E(\tilde Z_i) E(\tilde Z_i)^T = \begin{pmatrix} \kappa_{40} \sigma^4 & \kappa_{31} \sigma^2 \tau & \kappa_{22} \sigma^2 \tau^2 \\ \kappa_{31} \sigma^2 \tau & \kappa_{22} \sigma^2 \tau^2 & \kappa_{13} \sigma \tau^3 \\ \kappa_{22} \sigma^2 \tau^2 & \kappa_{13} \sigma \tau^3 & \kappa_{04} \tau^4 \end{pmatrix} - \tilde \Sigma \tilde \Sigma^T.$$

Here $$\kappa_{ij} = E \left[ \left( \frac{X_i}{\sigma} \right)^i \left( \frac{Y_i}{\tau} \right)^j \right]$$ indicates the $$ij$$th mixed standardized moment (about the mean, but we assumed mean zero at the onset).

Alternately, we have the factorization $$\text{Cov}(\tilde Z_i) = D(\sigma, \tau) [ K - R(\rho) R(\rho)^T ] D(\sigma, \tau), \quad (1)$$ where $$D(\sigma, \tau) = \text{diag}(\sigma^2, \sigma \tau, \tau^2)$$, $$R(\rho) = (1, \rho, 1)^T$$ and $$K = \begin{pmatrix} \kappa_{04} & \kappa_{31} & \kappa_{22} \\ \kappa_{31} & \kappa_{22} & \kappa_{13} \\ \kappa_{22} & \kappa_{13} & \kappa_{04} \end{pmatrix}.$$

Thus we have that $$Z_{11}$$ and $$Z_{12}$$, representing the sample variance of $$X$$ and the covariance of $$X,Y$$ are correlated unless $$\rho = \kappa_{31}$$. When $$V_i$$ is multivariate normal, this occurs only when $$\rho = 0$$.

2 The correlation coefficient

Now consider the transformation $$g(x, y, z) = (x, \frac{y}{\sqrt{z}\sqrt{x}})$$ on $$\tilde{S_n}$$. This provides the bivariate distribution of the sample correlation coefficient and the sample variance of x. By the delta method and asymptotic normality of $$S_n$$, $$\sqrt{n}( g(\tilde{S_n}) - (\rho, \sigma^2)^T ) \rightarrow N(0, \mathbf{J}(\tilde \Sigma)^T \tilde \Lambda \mathbf{J}(\tilde \Sigma)),$$ where $$\mathbf{J}(\tilde \Sigma) = [\nabla g_1^T, \nabla g_2^T]^T$$ is the jacobian of $$g$$.

I find (though you probably want to check my algebra..) that the gradient of the second component of $$g$$ is $$\nabla g_2 (\sigma^2, \rho \sigma \tau, \tau^2) = \left( -\frac{\rho}{2\sigma^2}, \frac{1}{\sigma \tau}, -\frac{\rho}{2 \tau^2} \right)^T,$$
So $$J(\sigma, \rho, \tau) = \begin{pmatrix} 1 & -\frac{\rho}{2\sigma^2} \\ 0 & \frac{1}{\sigma \tau} \\ 0 & -\frac{\rho}{2 \tau^2} \end{pmatrix}.$$

Putting it all together with the factorization in equation (1) yields

$$J(\sigma, \rho, \tau)^T D(\sigma, \tau) [ K - R(\rho) R(\rho)^T ] D(\sigma, \tau) J(\sigma, \rho, \tau).$$

Plugging in some easy to use numbers, say $$\sigma = \tau = 1$$ and $$\rho = .5$$, we'd have for $$J(\sigma, \rho, \tau)^T D(\sigma, \tau) [ K - R(\rho) R(\rho)^T ] D(\sigma, \tau)J(\sigma, \rho, \tau) = \begin{pmatrix} -1/4 & 1 & -1/4 \\ 1 & 0 & 0 \end{pmatrix} \mathbf I \Omega \mathbf I \begin{pmatrix} -1/4 & 1 \\ 1 & 0 \\ -1/4 & 0 \end{pmatrix} = \mathbf{Q},$$ where $$\Omega = K - R(\rho) R(\rho)^T$$ is generally some dense matrix. Courtesy of Mathematica, I expanded this product in terms of entries in $$K$$ and recount below $$Q_{12}$$ $$n \times Q_{12} = n \times \text{Cov}(r, s^2_x) = \kappa_{31} -\frac{\kappa_{04} + \kappa_{22}}{4} \quad (2)$$ which is an opaque expression in terms of the mixed moments, but certainly doesn't seem like it's going to be zero, generally.

3 Specializing to the normal case

Isserlis theorem provides a way to derive the mixed moments of a Gaussian. Again assuming $$\sigma = \tau = 1$$ and $$\rho = .5$$, we'd have $$\kappa_{31} = 3/2, \kappa_{04} = 3, \kappa_{22} = 3/2$$, thus $$Q_{12} = 3/2 - (3 + 3/2)/4 = 3/8 > 0$$, as you observe.

4 Simulation and Example

Below find a simulation verifying equation (1). For $$n=100$$ and $$n=1000$$ (in red and blue, respectively) iid observations from a multivariate normal, I derive the covariance of $$\sqrt{n} \tilde S_n$$ by bootstrap. The covariance between $$S_{xy}$$ and $$S_{xx}$$ is plotted on y axis as $$\rho$$ varies from $$-.9$$ to $$.9$$. The theoretical value from equation (1) and using facts about the 4th order moments of the bivariate Gaussian is plotted in a dashed black line. A fun exercise would be to try to find a family of copula that for any value of $$\rho$$ would render $$\text{Cov}(S_{xy}, S_{xx}) = 0$$...

library(mvtnorm)
library(tidyverse)
library(boot)
params = expand.grid(sx = 1, sy = 1, n = c(100, 1000), rho = seq(-.9, .9, by = .1), replicate = 1:10) %>% mutate(k04 = 3*sx^4, k31 = 3*sx*rho*sx*sy, q12 = k31 - rho*sx*sy)

Sn_tilde = function(dat, idx){
Sn = cov(dat[idx,,drop =FALSE])*sqrt(length(idx))
Sn[upper.tri(Sn, diag = TRUE)]
}

out = params %>% group_by_all() %>% do({
x = with(., rmvnorm(n = .$$n, sigma = matrix(c(sx^2, rho*sx*sy, rho*sx*sy, sy^2), nrow = 2))) colnames(x) = c('X', 'Y') b = boot(x, Sn_tilde, R = 500) cov_Sn = cov(b$$t)
rownames(cov_Sn) = colnames(cov_Sn) = c('Sxx', 'Sxy', 'Syy')
as_tibble(cov_Sn, rownames = 'j')
})

ggplot(filter(out,  j == 'Sxx'), aes(x = rho, y = Sxy, color = factor(n))) + geom_point(size = .5, alpha = .5) + geom_smooth(method = 'lm') + geom_line(data = filter(params, replicate == 1, n == 100), aes(y = q12), lty = 2, color = 'black') + theme_minimal() + ylab('Cov(Sxy, Sxx)')


$$^1$$ This heavily uses Michael Perlman's lecture notes on probability and mathematical statistics, which I really wish were available as a electronically so I could replace mine when they wear out...

• Thank you! However, I think there may be a misstep somewhere: in fact it seems that $Cov(r, s_x) \to 0$, not 3/8, empirically (although $Corr(r, s_x)$ does not because $Var(r, s_x) \to 0$ and $Var(r, s_x) \to 0$). – half-pass Apr 4 at 11:17
• I am going to ask a new question about this since I also don't know how to show $Cov(r, s_x) \to 0$. – half-pass Apr 4 at 11:45
• (+1) Very interesting post. It appears that for bivariate $N(0,1)$, expression $(2)$ evaluates to $3\rho - 1 - 0.5 \rho^2$. This leads to the result that if $\rho<0.35 \implies \text{Cov}(r, s^2_x) <0$ while if $\rho>0.35 \implies \text{Cov}(r, s^2_x) >0$. – Alecos Papadopoulos Apr 4 at 22:53
• @half-pass: The pair $(r, s_x)$ needs to be scaled up by $\sqrt n$ to have a (non-degenerate) limiting distribution. If you want to examine the correlation per se, you could use the result in section 1 and just modify the $g$-function in section 2 accordingly. – Andrew M Apr 5 at 15:44
• @AlecosPapadopoulos: expression 2 is already specialized to the case that $\tau^2 = \sigma^2 = 1$ and $\rho = .5$. If all you care about is the sign of the association between $s_x$ and $r$, could just examine the [1,2] entry in the difference $K - R(\rho)R(\rho)^T$ in equation 1 using facts about the mixed moments of a bivariate normal to plug in for $K$ as a function of $\rho$. – Andrew M Apr 5 at 15:47

Edit: This answer is incorrect. I'm not sure whether it's better to leave it here for the record, or to just delete it.

Yes, it does hold asymptotically regardless of the distribution of X and Y. I was on the right track with the Taylor expansion: • I look at eq.$(2.1)$, second line. We have a probability limit minus a fixed quantity. If we have consistency, the probability limit is equal to the fixed quantity. Therefore the 2nd line of eq. $(2.1)$ appears to be equal to zero... which should be expected, since the probability limit of $\text{Cov}(r,s_x)$ is equal to $\text{Cov}(\rho, \sigma_x)$. But both $\rho$ and $\sigma_x$ are constants, so their covariance is zero. It appears that the result obtained depends critically 1) on ignoring the Taylor remainder (whose sign we don't know) and (CONTD) – Alecos Papadopoulos Apr 2 at 23:01
• (CONTD) ...and 2) on applying the consistency property selectively on certain terms but not on others. Are you sure that these are valid manipulations? – Alecos Papadopoulos Apr 2 at 23:01
• Thanks for the pushback. I made that manipulation in the second line because, for general RVs U and V, plim E[U] E[V] = plim E[U] plim E[V], giving me the second term. But for the first term, plim E[ UV ] != plim E[U] plim E[V]. Therefore, I think I've applied consistency throughout, just on different steps. – half-pass Apr 3 at 17:23
• I am not sure, though... – half-pass Apr 3 at 17:25
• Okay -- this is indeed incorrect! However, I don't think it's because of the first manipulation but rather the fact that I ignored that Var(s_Y) -> 0 as well, leading to a tautology. Unfortunately I cannot downvote my own answer. – half-pass Apr 3 at 23:45

It will depend on the joint distribution. For the example you mention, the bivariate (zero-mean) Normal distribution is characterized by the $$\rho, \sigma_x, \sigma_y$$. It follows that one can have all possible combinations of values of these three parameters, implying that no relation between $$\rho$$ and the standard deviations can be established.

For other bivariate distributions, the correlation coefficient may be fundamentally a function of the standard deviations (essentially both will be functions of more primitive parameters), in which case one can examine whether a monotonic relation exists.

• I understand that the three parameters can have arbitrary relationships for the BVN distribution, but I don't think it follows that the sample estimates of these are asymptotically independent. – half-pass Apr 1 at 21:31