# Quantifying dependence of Cauchy random variables

Given two Cauchy random variables $$\theta_1 \sim \mathrm{Cauchy}(x_0^{(1)}, \gamma^{(1)})$$ and $$\theta_2 \sim \mathrm{Cauchy}(x_0^{(2)}, \gamma^{(2)})$$. That are not independent. The dependence structure of random variables can often be quantified with their covariance or correlation coefficient. However, these Cauchy random variables have no moments. Thus, covariance and correlation do not exist.

Are there other ways of representing the dependence of the random variables? Is it possible to estimate those with Monte Carlo?

Just because they don’t have a covariance doesn’t mean that the basic $$x^t\Sigma^{-1} x$$ structure usually associated with covariances can’t be used. In fact, the multivariate ($$k$$-dimensional) Cauchy can be written as:

$$f({\mathbf x}; {\mathbf\mu},{\mathbf\Sigma}, k)= \frac{\Gamma\left(\frac{1+k}{2}\right)}{\Gamma(\frac{1}{2})\pi^{\frac{k}{2}}\left|{\mathbf\Sigma}\right|^{\frac{1}{2}}\left[1+({\mathbf x}-{\mathbf\mu})^T{\mathbf\Sigma}^{-1}({\mathbf x}-{\mathbf\mu})\right]^{\frac{1+k}{2}}}$$

which I have lifted from the Wikipedia page. This is just a multivariate Student-$$t$$ distribution with one degree of freedom.

For the purposes of developing intuition, I would just use the normalized off-diagonal elements of $$\Sigma$$ as if they were correlations, even though they are not. They reflect the strength of the linear relationship between the variables in a way very similar to that of a correlation; $$\Sigma$$ has to be positive definite symmetric; if $$\Sigma$$ is diagonal, the variates are independent, etc.

Maximum likelihood estimation of the parameters can be done using the E-M algorithm, which in this case is easily implemented. The log of the likelihood function is:

$$\mathcal{L}(\mu, \Sigma) = -{n\over 2}|\Sigma| – {k+1 \over 2}\sum_{i=1}^n\log(1+s_i)$$

where $$s_i = (x_i-\mu)^T\Sigma^{-1}(x_i-\mu)$$. Differentiating leads to the following simple expressions:

$$\mu = \sum w_ix_i/\sum w_i$$

$$\Sigma = {1 \over n}\sum w_i(x_i-\mu)(x_i-\mu)^T$$

$$w_i = (1+k)/(1+s_i)$$

The E-M algorithm just iterates over these three expressions, substituting the most recent estimates of all the parameters at each step.

For more on this, see Estimation Methods for the Multivariate t Distribution, Nadarajah and Kotz, 2008.