Probability – In-depth Interpretation of Total Law of Covariance

covarianceprobabilityvariance-decomposition

let $X,Y,Z$ be random variables defined on the same probability space and let covariance of $X$ and $Y$ be finite, then the law of total covariance / covariance decomposition formula states:
\begin{align}
\text{Cov}(X,Y)=\underbrace{\mathbb{E}\big[\text{Cov}(X,Y\lvert Z)\big]}_{\text{(i)}}+\underbrace{\text{Cov}\big[\mathbb{E}(X\lvert Z),\mathbb{E}(Y\lvert Z)\big]}_{\text{(ii)}}
\end{align}
What is the interpretation of $\text{(i)}$ and $\text{(ii)}$?

My thoughts: in (ii) the two conditional expectations can be seen as random variables themselves, I also know that this is a generalization of the law of total variance / variance decomposition formula which can be shown by setting $X=Y$, where the interpretation is then that of a variation in $Y$ explained by $Z$ and unexplained by $Z$. But what is the correct interpretation in the above covariance formula for (i) and (ii)? Wikipedia offers a brief description which is not very satisfying.

Best Answer

The first term (i): $E[\operatorname{cov}(X,Y|Z)]$

Think of $\operatorname{cov}(X,Y)$ as a function of $Z$. As you examine different values of $Z$, you will correspondingly get a value for $\operatorname{cov}(X,Y)$. The expectation simply averages these different covariances with respect to $Z$.

The second term (ii): $\operatorname{cov}([E[X|Z],E[Y|Z])$

Think of $E[X|Z]$ and $E[Y|Z]$ as functions of $Z$. As you examine different values of $Z$, you will correspondingly get a value of $E[X|Z]$ and a value of $E[Y|Z]$ realized simultaneously. Therefore, for every value of $Z$, you will get an $(X, Y)$ coordinate. This term is simply the covariance of all these coordinate points.

Part 1

Let's solve this from first principles and then reflect on what the solution might teach us.

First Step: Finding the Marginal and Joint Density Functions

The marginal density function of $X$ is

$$f_X(x) = I_{[0,1]}(x),$$

the indicator function for its interval of support $[0,1]$.

The conditional density of $Y$ is

$$f_{Y|x}(y) = \phi(y;x,|x|)= \phi(y;x,x)$$

where $\phi(x;\mu,\sigma)$ is the Normal PDF with mean $\mu$ and standard deviation $\sigma$ (and $|x|=x$ since $0\le x\le 1$). Therefore the joint density function is

$$f_{X,Y}(x,y) = I_{[0,1]}(x) \phi(y;x,x).$$

The marginal density of $Y$ is obtained by integrating out $x$,

$$f_Y(y) = \int f_{X,Y}(x,y) dx = \int_0^1 \phi(y;x,x) dx.$$

On the principle of Just-in-time-computation, I will not evaluate these integrals until that might be needed. (It turns out it is not needed at all!) Note, therefore, that no calculations have yet been performed: everything is just applying definitions.

Warm-up: Reproducing the Moments in the Question

Let's proceed to compute moments. The moments of $X$ are the moments of the uniform distribution, which are easily calculated or looked up. (The mean is $\mu_X^{(1)}=1/2$ and the mean square is $\mu_X^{(2)}=1/3$, for instance.) The $k^\text{th}$ moment of $Y$, $\mu_Y^{(k)}$, can be computed as

$$\mu_Y^{(k)} = \int y^k f_Y(y) dy = \int y^k \int_0^1 \phi(y;x,x) dx dy.$$

Let's interchange the order of integration (Fubini's Theorem) to obtain

$$\mu_Y^{(k)} = \int_0^1 \left(\int y^k \phi(y;x,x) dy\right) dx.$$

The inner integral computes the $k^\text{th}$ moment of a Normal variable $Y$ with mean $x$ and variance $x^2$, so it is immediate that this moment is $x$ when $k=1$ and $x^2 + (x)^2 = 2x^2$ when $k=2$. Whence

$$\mu_Y^{(1)} = \int_0^1 (x) dx = \frac{1}{2}$$

and

$$\mu_Y^{(2)} = \int_0^1 (2x^2) dx = \frac{2}{3}.$$

We conclude the variance of $Y$ is $\mu_Y^{(2)} - \left(\mu_Y^{(1)}\right)^2 = 2/3 - (1/2)^2 = 5/12$.

Solution

Having established that this method of direct integration, using Fubini's Theorem, reproduces the calculations in the question, let's apply it to the covariance:

$$\mathbb{E}(XY) = \int xy f_Y(y) dy = \int_0^1 x \left( \int y \phi(y;x,x) dy \right) dx = \int_0^1 x(x)dx = \frac{1}{3}.$$

Therefore

$$\text{Cov}(X,Y) = \mathbb{E}(XY) - \mathbb{E}(X)\mathbb{E}(Y) = \frac{1}{3} - \left(\frac{1}{2}\right) \left(\frac{1}{2}\right)= \frac{1}{12}.$$

Reflection

It should now be evident that the reversal of the order of integration is really the same thing as the Law of Iterated Expectation (LIE). The strategy was to compute the covariance in terms of the raw moments $\mathbb{E}(XY)$, $\mathbb{E}(X)$, and $\mathbb{E}(Y)$ by exploiting the knowledge of the conditional moments of $Y$, given $X=x$, and integrating those over $x$. Thus we could have invoked the LIE at the outset in the form

$$\mathbb{E}(XY) = \mathbb{E}(\mathbb{E}(XY | X)) = \mathbb{E}(X\,\mathbb{E}(Y | X)) = \mathbb{E}(X(X)) = \mathbb{E}(X^2) = \frac{1}{3}$$

with a similar (but easier) manipulation to find $\mathbb{E}(Y)$. Comparison with the previous integral expression shows that the crucial idea is to factor $X$ out of the conditional expectation:

$$\mathbb{E}(XY | X) = X\,\mathbb{E}(Y | X).$$

That may be the main thing worth remembering from this exercise.

Although these results can be used (in a straightforward manner) to derive an iterated law for covariances, comparable to the law for variances quoted in the question, I find such rules to be a burden on the memory. By identifying and focusing on fundamental principles and ideas, we only need to remember and learn a few things--and we can learn them well. Everything else can then easily be derived as needed.

Part 2

Let $xz = y$, entailing $dy = x dz + z dx$. Since

$$\eqalign{ \phi(y;x,x)dxdy &= \frac{1}{x}\phi\left(\frac{y}{x};1,1\right) dxdy \\ &= \frac{1}{x}\phi\left(z;1,1\right) dx(x dz + z dx) \\ &= x\frac{1}{x}\phi\left(z;1,1\right) dx dz \\ &= \phi\left(z;1,1\right) dx dz, }$$

the joint distribution easily factors as

$$f_{X,Y}(x,y)dxdy = I_{[0,1]}(x) \left[\phi(z;1,1)dxdz\right] = \left[I_{[0,1]}(x)dx\right]\, \left[\phi(z;1,1)dz\right].$$

Since all probability calculations can now be carried out separately for $x$ and $z$, it follows that $X$ and $Z=Y/X$ are independent, QED.

(The manipulations with differentials are interpreted as wedge products. This is fully explained at https://stats.stackexchange.com/a/154298.)

As a check, we may simulate a large number of independent realizations of $(X,Y)$ and plot them. The following carries this out with R:

n <- 1e4
x <- runif(n)
y <- rnorm(n, x, x)
plot(x, y/x, cex=1/3, col="#00000020")

In any vertical slice the distributions of the $y_i/x_i$ look approximately the same. (They are all standard Normal with mean $1$, as we can see from their PDF $\phi(z-1)$.) That is what independence of $Z$ and $X$ means: the distribution of $Z$ does not vary with the value of $X$ (nor does the distribution $X$ vary with the value of $Z$).

Solved – Meaning of Square Root of Covariance / Precision Matrices

I will write matrix square roots of $\Sigma$ as $\Sigma=A A^T$, to be consistent with the Cholesky decomposition which is written as $\Sigma = L L^T$ where $L$ is lowtri (lower triangular). So let $X$ be a random vector with $\DeclareMathOperator{\E}{\mathbb{E}} \DeclareMathOperator{\Cov}{\mathbb{Cov}}\DeclareMathOperator{\Var}{\mathbb{Var}} \E X=\mu$ and $\Var X =\Sigma$. Let now $Z$ be a random vector with expectation zero and unit covariance matrix.

Note that there are (except for the scalar case) infinitely many matrix square roots. If we let $A$ be one of the, then we can find all the others as $A \mathcal{O}$ where $\mathcal{O}$ is any orthogonal matrix, that is, $\mathcal{O} \mathcal{O}^T = \mathcal{O}^T \mathcal{O} =I$. This is known as unitary freedom of square roots.

Let us look at some particular matrix square roots.

First a symmetric square root. Use the spectral decomposition to write $\Sigma = U \Lambda U^T = U\Lambda^{1/2}(U\Lambda^{1/2})^T$. Then $\Sigma^{1/2}=U\Lambda^{1/2}$ and this can be interpreted as the PCA (principal component analysis) of $\Sigma$.
The Cholesky decomposition $\Sigma=L L^T$ and $L$ is lowtri. We can represent $X$ as $X=\mu + L Z$. Multiplying out to get scalar equations, we get a triangular system in $Z$, which in the time series case can be interpreted as a MA (moving average) representation.
The general case $A= L \mathcal{O}$, using the above we can interpret this as a MA representation after rotating $Z$.