Mfrhtyj

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^tSigma^{-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}; {mathbfmu},{mathbfSigma}, k)= frac{Gammaleft(frac{1+k}{2}right)}{Gamma(frac{1}{2})pi^{frac{k}{2}}left|{mathbfSigma}right|^{frac{1}{2}}left[1+({mathbf x}-{mathbfmu})^T{mathbfSigma}^{-1}({mathbf x}-{mathbfmu})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) = -{nover 2}|Sigma| - {k+1 over 2}sum_{i=1}^nlog(1+s_i)$$

where $s_i = (x_i-mu)^TSigma^{-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.

While $text{cov}(X,Y)$ does not exist, for a pair of variates with Cauchy marginals, $text{cov}(Phi(X),Phi(Y))$ does exist for, e.g., bounded functions $Phi(cdot)$. Actually, the notion of covariance matrix is not well-suited to describe joint distributions in every setting, as it is not invariant under transformations.

Borrowing from the concept of copulas (which may also help in defining a joint distribution¹ for $(X,Y)$), one can turn $X$ and $Y$ into Uniform $(0,1)$ variates, by using their marginal cdfs, $Phi_X(X)simmathcal{U}(0,1)$ and $Phi_Y(Y)simmathcal{U}(0,1)$, and look at the covariance or correlation of the resulting variates.

¹For instance, when $X$ and $Y$ are both standard Cauchys,$$Z_X=Phi^{-1}({argtan(X)/pi+1}/2)$$is distributed as a standard Normal, and the joint distribution of $(Z_X,Z_Y)$ can be chosen to be a joint Normal
$$(Z_X,Z_Y) sim mathcal{N}_2(0_2,Sigma)$$This is a Gaussian copula.