Covariance between a normal variable (x1) and a sum including it (x1 + x2)

Question

I'm interested in the joint distribution of two variables, $x_1$ and $x_2$.

$$ x_1 \sim Normal(0, \sigma^2_1);\\ \epsilon \sim Normal(0, \sigma^2_{\epsilon});\\ x_2 = x_1 + \epsilon; $$

as a bivariate normal distribution

$$ \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \sim Normal( 0, \Sigma); \\ $$

$$ \Sigma = \begin{bmatrix} \sigma^2_1 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma^2_2 \end{bmatrix} $$

$$ \sigma^2_2 = \sigma^2_1 + \sigma^2_{\epsilon} $$

I can simulate from this distribution quite easily, but would like the closed-form solution for the correlation coefficient $\rho$. I can see that the $\rho$ should be proportional to $\sigma^2_1$: with higher $\sigma^2_1$, the signal present in $x_1$ is stronger and less influenced by the noise $\epsilon$.

Simulations show that the covariance is actually very close to the variance $\sigma^2_1$...

library(tidyverse)
vals = 2^seq(1, 8)
vals = seq(1, 10, 1)
d = expand.grid(var1=vals, var2=vals)

d$cov = pmap_dbl(d, function(var1, var2){
  # Generate random variables and check covariance between x1 and x3 (= x1 + x2)
  n = 1000
  x1 = rnorm(n, 0, sqrt(var1))
  x2 = rnorm(n, 0, sqrt(var2))
  x3 = x1 + x2
  cov(x1, x3)
})
ggplot(d, aes(var1, cov, color=factor(var2))) +
  geom_path() +
  geom_abline(intercept=0, slope=1, size=1, linetype='dashed') +
  labs(x='Var(x1)', y='Cov(x1, x2)', color='Var(x2)') +
  coord_equal()

jet.colors <- c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")
ggplot(d, aes(var1, var2, fill=cov)) +
  geom_tile() +
  scale_fill_gradientn(colors=jet.colors) +
  labs(x='Var(x1)', fill='Cov(x1, x2)', y='Var(x2)') +
  coord_equal()

...but not quite right either, presumably as I need to incorporate $\sigma^2_{\epsilon}$ into the equation somewhere.

Is there a known analytic solution to this?

Answer 1

Suppose $x_1$ and $\epsilon$ are independent, then

$$\mathrm{Cov}\left(\matrix{x_1\\ \epsilon}\right) = \left (\matrix{\sigma_1^2 & 0 \\0&\sigma_\epsilon^2}\right)$$ $$\left(\matrix{x_1\\ x_2}\right) = \left (\matrix{1 & 0 \\1 &1}\right)\left(\matrix{x_1\\ \epsilon}\right) $$

So $$\mathrm{Cov}\left(\matrix{x_1\\ x_2}\right) = \left (\matrix{1 & 0 \\1 &1}\right)\left (\matrix{\sigma_1^2 & 0 \\0&\sigma_\epsilon^2}\right)\left (\matrix{1 & 1 \\0 &1}\right) = \left (\matrix{\sigma_1^2 & \sigma_1^2 \\\sigma_1^2 &\sigma_1^2+\sigma_\epsilon^2}\right)$$

So the correlation coefficient $ρ = \frac {\sigma_1^2}{\sqrt{\sigma_1^2 (\sigma_1^2+\sigma_\epsilon^2)}} = \frac {\sigma_1}{\sqrt{\sigma_1^2+\sigma_\epsilon^2}}$

Answer 2

Thank you to @a_statistician for the correct analytical solution. I've since realised why my simulation results haven't agreed with this analytic result (cov(x1, x2) == var(x1)): the effect of sampling error is surprisingly strong here, even with large N, and so the correlation between $x_1$ and $\epsilon$ weren't necessarily 0, and the sample variances of $x_1$ and $\epsilon$ weren't necessarily their population values of $\sigma^2_1$ and $\sigma^2_{\epsilon}$.

Here's the corrected simulation code.

library(tidyverse)
orthogonalise = function(x1, x2){
  X = data.frame(x1, x2) %>% as.matrix
  pca = prcomp(X)
  pX = X %*% pca$rotation
  return(list(x1=pX[,1], x2=pX[,2]))
}
exact.orthogonal.samples = function(var1, var2){
  n = 1000
  x1 = rnorm(n, 0, 1)
  x2 = rnorm(n, 0, 1)
  # Remove correlation
  ox = orthogonalise(x1, x2)
  x1 = ox$x1
  x2 = ox$x2
  # Set variance
  x1 = x1 * (sqrt(var1)/sd(x1))
  x2 = x2 * (sqrt(var2)/sd(x2))
  return(list(x1=x1, x2=x2))
}

vals = 2^seq(1, 8)
vals = seq(1, 10, 1)
d = expand.grid(var1=vals, var2=vals)
d$cov = res = pmap_dbl(d, function(var1, var2){
  # Generate random variables and check covariance between x1 and x2 (= x1 + e)
  X = exact.orthogonal.samples(var1, var2)
  x1 = X$x1
  e = X$x2
  x2 = x1 + e
  cov(x1, x2)
})


jet.colors <- c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000")
ggplot(d, aes(var1, var2, fill=cov)) +
  geom_tile() +
  scale_fill_gradientn(colors=jet.colors) +
  labs(x='Var(x1)', fill='Cov(x1, x2)', y='Var(x2)') +
  coord_equal()