**3. Generalized Cauchy Distribution**

**Proposition 2.** *If X* ∼ *ESN* (*<sup>ε</sup>*1) *and Y* ∼ *ESN* (*<sup>ε</sup>*2) *are independent random variables, then Z d* = *X Y has what we call a generalized epsilon–skew–Cauchy (GESC) distribution with parameters ε*1 *and ε*2 *and with density given by*

$$f\_Z\left(z; \varepsilon\_1, \varepsilon\_2\right) = \frac{1}{2\pi \left[\left(\frac{1}{1+\varepsilon\_2}\right)^2 + \left(\frac{z}{1+\varepsilon yn(z)\varepsilon\_1}\right)^2\right]} + \frac{1}{2\pi \left[\left(\frac{1}{1-\varepsilon\_2}\right)^2 + \left(\frac{z}{1-\varepsilon yn(z)\varepsilon\_1}\right)^2\right]},$$

*where z* ∈ R*,* |*<sup>ε</sup>*1| < 1*,* |*<sup>ε</sup>*2| < 1 *and we write Z* ∼ *GESC*(*<sup>ε</sup>*1,*ε*2)*.*

**Proof.** Arguing as in Proposition 1, we have that

$$f\_{\mathbf{Z}}\left(z; \varepsilon\_{1}, \varepsilon\_{2}\right) = \int\_{-\infty}^{\infty} \frac{|w|}{2\pi} \exp\left\{-\frac{1}{2} \left[\frac{z^{2}}{\left(1 - \operatorname{sgn}\left(zw\right)\varepsilon\_{1}\right)^{2}} + \frac{1}{\left(1 - \operatorname{sgn}\left(w\right)\varepsilon\_{2}\right)^{2}}\right] w^{2}\right\} \,\mathrm{d}w.$$

By considering separately the cases in which *z* ≥ 0 and *z* < 0, result follows directly.

From Proposition 2 two special cases are obtained directly,

1. If *ε*2 = 0 a generalized Cauchy (GC) distribution is obtained. In this case we write *Z* ∼ *GC*(*ε*), and its pdf is given by

$$f\_{Z}\left(z;\varepsilon\right) = \frac{1}{2\pi\left[1 + \left(\frac{z}{1 + \varsigma\eta n(z)\varepsilon}\right)^{2}\right]} + \frac{1}{2\pi\left[1 + \left(\frac{z}{1 - \varsigma\eta n(z)\varepsilon}\right)^{2}\right]},$$

where *z* ∈ R, 0 < *ε* < 1 and *Z d* = *X*/*Y* where *X* ∼ *ESN* (*ε*) and *Y* ∼ *N*(0, 1) are independent random variables.

2. If *ε*1 = *ε*2 = *ε* an epsilon–skew–Cauchy distribution is obtained. In this case we write *Z* ∼ *ESC*2(*ε*), and its pdf is given by

$$f\_{Z}\left(z;\varepsilon\right) = \frac{1}{2\pi\left[\left(\frac{1}{1+\varepsilon}\right)^{2} + \left(\frac{z}{1+\operatorname{sgn}(z)\varepsilon}\right)^{2}\right]} + \frac{1}{2\pi\left[\left(\frac{1}{1-\varepsilon}\right)^{2} + \left(\frac{z}{1-\operatorname{sgn}(z)\varepsilon}\right)^{2}\right]},$$

where *z* ∈ R, |*ε*| < 1 and *Z d* = *X*/*Y* where *X* ∼ *ESN* (*ε*) and *Y* ∼ *ESN* (*ε*) are independent random variables.

### **4. General Bivariate Mudholkar-Hutson Distribution**

Define *Z*1, *Z*2, *Z*3 to be i.i.d. standard normal random variables. For *i* = 1, 2, 3 define

$$\mathcal{U}\_i = \begin{cases} -\alpha\_i & \text{with probability } \gamma\_i\\ \beta\_i & \text{with probability } 1 - \gamma\_i \end{cases}$$

where *α*1, *α*2, *α*3, *β*1, *β*2, *β*3 are positive numbers and 0 < *γ*1, *γ*2, *γ*3 < 1. So, the parameters *αi* and *βi* indicate the propensity in which the discrete random variable takes negative and positive values, respectively. The parameters *γi* control how often negative and positive values are taken by *Ui*.

In addition assume that all six random variables *Z*1, *Z*2, *Z*3, *U*1, *U*2 and *U*3 are independent, and define 

$$f(X,Y) = \left(\frac{\left\langle L\_1 \left| Z\_1 \right|, \left| L\_2 \left| Z\_2 \right| \right| \right\rangle}{\left\langle L\_3 \left| Z\_3 \right|, \left| L\_3 \left| Z\_3 \right| \right| \right.}\right). \tag{2}$$

The model (2) is highly flexible since it allows for different behavior in each of the four quadrants of the plane. From (2) it may be recognized that only fractional moments of *X* and *Y* exist.

Note that if we define (*<sup>W</sup>*1, *<sup>W</sup>*2) = |*Z*1| |*<sup>Z</sup>*3|, |*<sup>Z</sup>*2| |*<sup>Z</sup>*3|, it is readily verified that

$$f\_{\mathcal{W}\_1, \mathcal{W}\_2}(w\_1, w\_2) = \frac{2}{\pi} \left(1 + w\_1^2 + w\_2^2\right)^{-3/2} I\left(w\_1 > 0, w\_2 > 0\right),\tag{3}$$

in which case we say that (*<sup>W</sup>*1, *<sup>W</sup>*2) has a standard bivariate half-Cauchy distribution.

Using (3) and conditioning on *U*1, *U*2 and *U*3 we obtain the density of (*<sup>X</sup>*,*<sup>Y</sup>*) in the albeit complicated form:

*f*(*<sup>x</sup>*, *y*) = 2*π γ*1*γ*2*γ*3 *α*23 *α*1*α*2 1 + *α*3 *α*1 *x*2 + *α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( + 2*π <sup>γ</sup>*1*γ*2(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *α*1*α*2 1 + *β*3 *α*1 *x*2 + *β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 2*π <sup>γ</sup>*1(<sup>1</sup> − *<sup>γ</sup>*2)*<sup>γ</sup>*3 *α*23 *<sup>α</sup>*1*β*2 1 + *α*3 *α*1 *x*2 + *α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 2*π <sup>γ</sup>*1(<sup>1</sup> − *<sup>γ</sup>*2)(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *<sup>α</sup>*1*β*2 1 + *β*3 *α*1 *x*2 + *β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 2*π* (1 − *<sup>γ</sup>*1)*<sup>γ</sup>*2*γ*3 *α*23 *β*1*α*2 1 + *α*3 *β*1 *x*2 + *α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 2*π* (1 − *<sup>γ</sup>*1)*<sup>γ</sup>*2(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *β*1*α*2 1 + *β*3 *β*1 *x*2 + *β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 2*π* (1 − *<sup>γ</sup>*1)(<sup>1</sup> − *<sup>γ</sup>*2)*<sup>γ</sup>*3 *α*23 *β*1*β*2 1 + *α*3 *β*1 *x*2 + *α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 2*π* (1 − *<sup>γ</sup>*1)(<sup>1</sup> − *<sup>γ</sup>*2)(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *β*1*β*2 1 + *β*3 *β*1 *x*2 + *β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( . (4)

Some special cases which might be considered include the following:

1. Mudholkar and Hutson type. For this we set: *αi* = 1 + *εi*, *βi* = 1 − *εi* and *γi* = (1 + *εi*) /2 for *i* = 1, 2, 3.

In this case the density (4) simplifies somewhat to become:

*f*(*<sup>x</sup>*, *y*) = 14*π* (1 + *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>+*ε*3 1+*ε*1 *x*2 + <sup>1</sup>+*ε*3 1+*ε*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( + 14*π* (1 − *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>−*ε*<sup>3</sup> 1+*ε*1 *x*2 + <sup>1</sup>−*ε*<sup>3</sup> 1+*ε*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 14*π* (1 + *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>+*ε*3 1+*ε*1 *x*2 + <sup>1</sup>+*ε*3 1−*ε*<sup>2</sup> *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 14*π* (1 − *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>−*ε*<sup>3</sup> 1+*ε*1 *x*2 + <sup>1</sup>−*ε*<sup>3</sup> 1−*ε*<sup>2</sup> *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 14*π* (1 + *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>+*ε*3 1−*ε*<sup>1</sup> *x*2 + <sup>1</sup>+*ε*3 1+*ε*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 14*π* (1 − *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>−*ε*<sup>3</sup> 1−*ε*<sup>1</sup> *x*2 + <sup>1</sup>−*ε*<sup>3</sup> 1+*ε*2 *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 14*π* (1 + *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>+*ε*3 1−*ε*<sup>1</sup> *x*2 + <sup>1</sup>+*ε*3 1−*ε*<sup>2</sup> *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 14*π* (1 − *<sup>ε</sup>*3)<sup>3</sup> 1 + <sup>1</sup>−*ε*<sup>3</sup> 1−*ε*<sup>1</sup> *x*2 + <sup>1</sup>−*ε*<sup>3</sup> 1−*ε*<sup>2</sup> *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( . (5)

A further specialization of the density (5) can be considered, as follows. 2. Homogenous Mudholkar and Hutson type. For this we set: *αi* = 1 + *ε*, *βi* = 1 − *ε* and *γi* = (1 + *ε*) /2 for *i* = 1, 2, 3.

This homogeneity results in a little simplification of (4), thus:

*f* (*<sup>x</sup>*, *y*;*ε*) = 14*π* )(1 + *ε*)<sup>3</sup> + (1 − *ε*)<sup>3</sup> 1 + *x*2 + *y*<sup>2</sup>−3/2 *<sup>I</sup>*(*x* > 0, *y* > 0)\* + 14*π* (1 − *ε*)<sup>3</sup> 1 + <sup>1</sup>−*<sup>ε</sup>* 1+*ε x*2 + <sup>1</sup>−*<sup>ε</sup>* 1+*ε y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 14*π* (1 + *ε*)<sup>3</sup> 1 + *x*2 + <sup>1</sup>+*ε* 1−*ε y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 14*π* (1 − *ε*)<sup>3</sup> 1 + <sup>1</sup>−*<sup>ε</sup>* 1+*ε x*2 + *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 14*π* (1 + *ε*)<sup>3</sup> 1 + <sup>1</sup>+*ε* 1−*ε x*2 + *y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 14*π* (1 − *ε*)<sup>3</sup> 1 + *x*2 + <sup>1</sup>−*<sup>ε</sup>* 1+*ε y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 14*π* (1 + *ε*)<sup>3</sup> 1 + <sup>1</sup>+*ε* 1−*ε x*2 + <sup>1</sup>+*ε* 1−*ε y*<sup>2</sup> <sup>−</sup>3/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( . (6)

It is easy to see that the parameter *ε* is not identifiable in (6) because *f* (*<sup>x</sup>*, *y*;*ε*) = *f* (*<sup>x</sup>*, *y*; <sup>−</sup>*<sup>ε</sup>*). An adjustment to ensure identifiability involves introducing the constraint *ε* ≥ 0.

3. Equal weights. In this case we assume that *α*1, *α*2, *α*3, *β*1, *β*2, *β*3 are positive numbers and *γi* = 1/2 for *i* = 1, 2, 3.

*f* (*<sup>x</sup>*, *y*) = 14*π α*23 *α*1*α*2 1 + *α*3 *α*1 *x*2 + *α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 + *β*23 *β*1*β*2 1 + *β*3 *β*1 *x*2 + *β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2( ×*I* (*x* > 0, *y* > 0) + 14*π β*23 *α*1*α*2 1 + *β*3 *α*1 *x*2 + *β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 + *α*23 *β*1*β*2 1 + *α*3 *β*1 *x*2 + *α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2( ×*I* (*x* < 0, *y* < 0) + 14*π α*23 *<sup>α</sup>*1*β*2 1 + *α*3 *α*1 *x*2 + *α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2 + *β*23 *β*1*α*2 1 + *β*3 *β*1 *x*2 + *β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2( ×*I* (*x* > 0, *y* < 0) + 14*π α*23 *β*1*α*2 1 + *α*3 *β*1 *x*2 + *α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>3/2 + *β*23 *<sup>α</sup>*1*β*2 1 + *β*3 *α*1 *x*2 + *β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>3/2( ×*I* (*x* < 0, *y* > <sup>0</sup>). (7)

The pdf (7) is not identifiable because the values of *αi* can be interchanged with those of *βi* and *f*(*<sup>x</sup>*, *y*) does not change. Moreover, multiplying all of the *α*'s and *β*'s by a constant does not change *f*(*<sup>x</sup>*, *y*). So, one way to ge<sup>t</sup> identifiability in the model (7) is to set *αi* = *βi* (*i* = 1, 2, 3) and *α*3 equal to 1. In that case, Equation (7) takes the form

$$f\left(x,y\right) = \frac{1}{2\pi} \frac{1}{\alpha\_1 \alpha\_2} \left[1 + \left(\frac{x}{\alpha\_1}\right)^2 + \left(\frac{y}{\alpha\_2}\right)^2\right]^{-3/2}.$$

.

However, this is then recognizable as being simply a scaled version of the standard bivariate Cauchy density (compare with Equation (3)).

We now consider the marginal densities for the random variable (*<sup>X</sup>*,*<sup>Y</sup>*) defined by (2). From (2) we have *X* = *<sup>U</sup>*1|*<sup>Z</sup>*1| *<sup>U</sup>*3|*<sup>Z</sup>*3| and *Y* = *<sup>U</sup>*2|*<sup>Z</sup>*2| *<sup>U</sup>*3|*<sup>Z</sup>*3| and the density of *W*1 = |*<sup>Z</sup>*1| |*<sup>Z</sup>*3| = *Z*1 *Z*3 is a standard half-Cauchy density, i.e., *fW*1 (*<sup>w</sup>*1) = 2*π* (1 + *w*21)−1.

Consequently, the density of *X* is of the form

$$\begin{array}{rcl}f\_{\mathbf{X}}\left(\mathbf{x}\right) &=& \frac{2}{\pi}\gamma\_{1}\gamma\_{3}\frac{a\_{3}}{a\_{1}}\left[1+\left(\frac{a\_{3}}{a\_{1}}\mathbf{x}\right)^{2}\right]^{-1}I\left(\mathbf{x}>0\right)+\frac{2}{\pi}\gamma\_{1}(1-\gamma\_{3})\frac{\beta\_{3}}{a\_{1}}\left[1+\left(\frac{\beta\_{3}}{a\_{1}}\mathbf{x}\right)^{2}\right]^{-1}I\left(\mathbf{x}<0\right) \\ &+\frac{2}{\pi}(1-\gamma\_{1})\left\{\gamma\_{3}\frac{a\_{3}}{\overline{\rho}\_{1}}\left[1+\left(\frac{a\_{3}}{\overline{\rho}\_{1}}\mathbf{x}\right)^{2}\right]^{-1}I\left(\mathbf{x}<0\right)+(1-\gamma\_{3})\frac{\beta\_{3}}{\overline{\rho}\_{1}}\left[1+\left(\frac{\beta\_{3}}{\overline{\rho}\_{1}}\mathbf{x}\right)^{2}\right]^{-1}I\left(\mathbf{x}>0\right)\right\},\end{array}$$

which is a mixture of half-Cauchy densities. The density of *Y* is then of the same form with *α*1, *β*1 replaced by *α*2, *β*2.

To ge<sup>t</sup> a more direct bivariate version of the original Mudholkar-Hutson density in which the *α*'s, *β*'s and *γ*'s were related to each other, we go back to the original representation, i.e., Equation (2), but now we set *α*1 = (1 + *<sup>ε</sup>*1), *β*1 = (1 − *<sup>ε</sup>*1), *γ*1 = (1 + *<sup>ε</sup>*1) /2, *α*2 = (1 + *<sup>ε</sup>*2), *β*2 = (1 − *<sup>ε</sup>*2), *γ*2 = (1 + *<sup>ε</sup>*2) /2 and *α*3 = *β*3 = *γ*3 = 1. So the density is of the form

$$\begin{split} f\left(\mathbf{x},y\right) &= \ \frac{1}{2\pi} \left[1+\left(\frac{\mathbf{x}}{1+\varepsilon\_{1}}\right)^{2}+\left(\frac{\mathbf{y}}{1+\varepsilon\_{2}}\right)^{2}\right]^{-3/2}I\left(\mathbf{x}>0,y>0\right) \\ &\quad + \frac{1}{2\pi} \left[1+\left(\frac{\mathbf{x}}{1+\varepsilon\_{1}}\right)^{2}+\left(\frac{\mathbf{y}}{1-\varepsilon\_{2}}\right)^{2}\right]^{-3/2}I\left(\mathbf{x}>0,y<0\right) \\ &\quad + \frac{1}{2\pi} \left[1+\left(\frac{\mathbf{x}}{1-\varepsilon\_{1}}\right)^{2}+\left(\frac{\mathbf{y}}{1+\varepsilon\_{2}}\right)^{2}\right]^{-3/2}I\left(\mathbf{x}<0,y>0\right) \\ &\quad + \frac{1}{2\pi} \left[1+\left(\frac{\mathbf{x}}{1-\varepsilon\_{1}}\right)^{2}+\left(\frac{\mathbf{y}}{1-\varepsilon\_{2}}\right)^{2}\right]^{-3/2}I\left(\mathbf{x}<0,y<0\right) \end{split}$$

with marginal density

$$f\_{\mathcal{X}}(\mathbf{x}) = \frac{1}{\pi} \left[ 1 + \left( \frac{\mathbf{x}}{1 + \varepsilon\_1} \right)^2 \right]^{-1} I\left(\mathbf{x} > 0\right) + \frac{1}{\pi} \left[ 1 + \left( \frac{\mathbf{x}}{1 - \varepsilon\_1} \right)^2 \right]^{-1} I\left(\mathbf{x} < 0\right) \dots$$

A more general version of the density with *γ*3 = 1, is of the form: (without loss of generality set *α*3 = 1)

$$\begin{split} f(\mathbf{x},y) &= \quad \frac{2}{\pi} \left\{ \gamma\_1 \gamma\_2 \frac{1}{a\_1 a\_2} \left[ 1 + \left( \frac{\mathbf{x}}{a\_1} \right)^2 + \left( \frac{\mathbf{y}}{a\_2} \right)^2 \right]^{-3/2} I(\mathbf{x} > 0, y > 0) \right\} \\ &+ \frac{2}{\pi} \left\{ \gamma\_1 (1 - \gamma\_2) \frac{1}{a\_1 \beta\_2} \left[ 1 + \left( \frac{\mathbf{x}}{a\_1} \right)^2 + \left( \frac{\mathbf{y}}{\beta\_2} \right)^2 \right]^{-3/2} I(\mathbf{x} > 0, y < 0) \right\} \\ &+ \frac{2}{\pi} \left\{ (1 - \gamma\_1) \gamma\_2 \frac{1}{\beta\_1 a\_2} \left[ 1 + \left( \frac{\mathbf{x}}{\beta\_1} \right)^2 + \left( \frac{\mathbf{y}}{a\_2} \right)^2 \right]^{-3/2} I(\mathbf{x} < 0, y > 0) \right\} \\ &+ \frac{2}{\pi} \left\{ (1 - \gamma\_1)(1 - \gamma\_2) \frac{1}{\beta\_1 \beta\_2} \left[ 1 + \left( \frac{\mathbf{x}}{\beta\_1} \right)^2 + \left( \frac{\mathbf{y}}{\beta\_2} \right)^2 \right]^{3/2} I(\mathbf{x} < 0, y < 0) \right\} \end{split}$$

with corresponding marginal

$$f\_{\mathcal{X}}\left(\mathbf{x}\right) = \frac{2}{\pi}\gamma\_1\frac{1}{a\_1}\left[1 + \left(\frac{\mathbf{x}}{a\_1}\right)^2\right]^{-1}I\left(\mathbf{x} > 0\right) + \frac{2}{\pi}(1-\gamma\_1)\frac{1}{\beta\_1}\left[1 + \left(\frac{\mathbf{x}}{\beta\_1}\right)^2\right]^{-1}I\left(\mathbf{x} < 0\right)\dots$$

Needless to say we can consider analogous models in which, instead of assuming that (*<sup>W</sup>*1, *<sup>W</sup>*2) = |*Z*1| |*<sup>Z</sup>*3| , |*<sup>Z</sup>*2| |*Z*3| has a density of the form (3) we assume that it has another bivariate density with support R<sup>+</sup> × R+, e.g., bivariate normal restricted to R<sup>+</sup> × R+, or bivariate Pareto, etc.

Another general bivariate MH model which includes the bivariate skew-Cauchy distribution given by (2) is the bivariate skew *t* model. This model can be obtained replacing |*<sup>Z</sup>*3| by *V*3 ∼ *χ*2*ν* in Equation (2). Thus, if we assume that all six random variables *Z*1, *Z*2 iid∼ *N* (0, <sup>1</sup>), *V*3 ∼ *<sup>χ</sup>*2*ν*; *U*1, *U*2 and *U*3 are independent, and define

$$\chi(X,Y) = \left(\frac{\mathcal{U}\_1 \, |Z\_1|}{\mathcal{U}\_3 \sqrt{V\_3/\nu}}, \frac{\mathcal{U}\_2 \, |Z\_2|}{\mathcal{U}\_3 \sqrt{V\_3/\nu}}\right),\tag{8}$$

then, because (*<sup>W</sup>*1, *<sup>W</sup>*2) = |*<sup>Z</sup>*1| √*V*3/*ν* , |*<sup>Z</sup>*2| <sup>√</sup>*V*3/*ν* has density

$$f\_{\mathcal{W}\_1, \mathcal{W}\_2}(w\_1, w\_2) = \frac{2}{\pi} \left( 1 + \frac{w\_1^2 + w\_2^2}{\nu} \right)^{-(\nu + 2)/2} I\left(w\_1 > 0, w\_2 > 0\right),$$

the density of (*<sup>X</sup>*,*<sup>Y</sup>*) is

*f* (*<sup>x</sup>*, *y*) = *f* (*<sup>x</sup>*, *y* |*<sup>u</sup>*1, *u*2, *u*3 ) *f* (*<sup>u</sup>*1, *u*2, *<sup>u</sup>*3) d*<sup>u</sup>*1d*<sup>u</sup>*2d*<sup>u</sup>*3 = *u*23 *u*1*u*2 *fW*1,*W*2 *u*3 *u*1 *x*, *u*3 *u*2 *y u*1, *u*2, *u*3 *f* (*<sup>u</sup>*1, *u*2, *<sup>u</sup>*3) d*<sup>u</sup>*1d*<sup>u</sup>*2d*<sup>u</sup>*3 = 2*π u*23 *u*1*u*2 1 + 1*ν u*3 *u*1 *x*2 + 1*ν u*3 *u*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *f* (*<sup>u</sup>*1, *u*2, *<sup>u</sup>*3) d*<sup>u</sup>*1d*<sup>u</sup>*2d*<sup>u</sup>*3 = 2 *π γ*1*γ*2*γ*3 *α*23 *α*1*α*2 1 + 1*ν α*3 *α*1 *x*2 + 1*ν α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( + 2*π <sup>γ</sup>*1*γ*2(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *α*1*α*2 1 + 1*ν β*3 *α*1 *x*2 + 1*ν β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 2*π <sup>γ</sup>*1(<sup>1</sup> − *<sup>γ</sup>*2)*<sup>γ</sup>*3 *α*23 *<sup>α</sup>*1*β*2 1 + 1*ν α*3 *α*1 *x*2 + 1*ν α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 2*π <sup>γ</sup>*1(<sup>1</sup> − *<sup>γ</sup>*2)(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *<sup>α</sup>*1*β*2 1 + 1*ν β*3 *α*1 *x*2 + 1*ν β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 2*π* (1 − *<sup>γ</sup>*1)*<sup>γ</sup>*2*γ*3 *α*23 *β*1*α*2 1 + 1*ν α*3 *β*1 *x*2 + 1*ν α*3 *α*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* < 0, *y* > 0)( + 2*π* (1 − *<sup>γ</sup>*1)*<sup>γ</sup>*2(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *β*1*α*2 1 + 1*ν β*3 *β*1 *x*2 + 1*ν β*3 *α*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* > 0, *y* < 0)( + 2*π* (1 − *<sup>γ</sup>*1)(<sup>1</sup> − *<sup>γ</sup>*2)*<sup>γ</sup>*3 *α*23 *β*1*β*2 1 + 1*ν α*3 *β*1 *x*2 + 1*ν α*3 *β*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* < 0, *y* < 0)( + 2*π* (1 − *<sup>γ</sup>*1)(<sup>1</sup> − *<sup>γ</sup>*2)(<sup>1</sup> − *<sup>γ</sup>*3) *β*23 *β*1*β*2 1 + 1*ν β*3 *β*1 *x*2 + 1*ν β*3 *β*2 *y*<sup>2</sup> <sup>−</sup>(*ν*+<sup>2</sup>)/2 *<sup>I</sup>*(*x* > 0, *y* > 0)( . (9)

The pair of variables in (8) allows one to model a wider variety of paired data sets than that given in (2) because it also can model light and heavy "tails" in a different way for each quadrant of the coordinate axis. Additionally, it is possible to compute the *<sup>r</sup>*−order moments.

**Proposition 3.** *The expected value of* (*<sup>X</sup>*,*<sup>Y</sup>*) *in (8) is given by*

$$\mathbb{E}\left(X,Y\right) = \frac{\sqrt{\nu}}{\sqrt{\pi}} \left(\beta\_1 \left(1-\gamma\_1\right) - a\_1\gamma\_1, \beta\_2 \left(1-\gamma\_2\right) - a\_2\gamma\_2\right) \left(\frac{1-\gamma\_3}{\beta\_3} - \frac{\gamma\_3}{a\_3}\right) \frac{\Gamma\left(\nu/2 - 1/2\right)}{\Gamma\left(\nu/2\right)}\lambda$$

*provided that ν* > 1*.* **Proof.** For *i* = 1, 2 it is immediate that E (*Ui*) E (|*Zi*|) = (*βi* (1 − *<sup>γ</sup>i*) − *<sup>α</sup>iγi*) √2/ √ *π* and E *U*−<sup>1</sup> 3 = (1 − *<sup>γ</sup>*3) /*β*3 − *γ*3/*<sup>α</sup>*3.

On the other hand, since *V*3 ∼ *χ*2 *ν* then

$$\mathbb{E}\left(\frac{\sqrt{\nu}}{\sqrt{V\_3}}\right) = \int\_0^\infty \frac{\sqrt{\nu}}{\sqrt{\infty}2^{\nu/2}\Gamma\left(\nu/2\right)} e^{-\mathbf{x}/2} d\mathbf{x} \stackrel{\nu \ge 1}{=} \frac{\sqrt{\nu}\Gamma\left(\nu/2 - 1/2\right)}{\sqrt{2}\Gamma\left(\nu/2\right)}.$$

Therefore, from (8) we have

$$\mathbb{E}\left(X,Y\right) = \left(\mathbb{E}\left(\mathcal{U}\_1\right)\mathbb{E}\left(|Z\_1|\right), \mathbb{E}\left(\mathcal{U}\_2\right)\mathbb{E}\left(|Z\_2|\right)\right) \mathbb{E}\left(\frac{1}{\mathcal{U}\_3}\right) \mathbb{E}\left(\frac{\sqrt{\nu}}{\sqrt{\mathcal{V}\_3}}\right)$$

and the result is obtained straightforward.

Following the proof of the previous proposition it is possible to obtain the *<sup>r</sup>*−order moments provided that *ν* > *r*.

In applications, it will usually be appropriate to augmen<sup>t</sup> these models by the introduction of location, scale and rotation parameters, i.e., to consider

$$\left(\widetilde{X}, \widetilde{Y}\right) = \mu + \Sigma^{1/2} \left(X, Y\right)\_{\mathcal{A}}$$

where *μ* ∈ (− <sup>∞</sup>, ∞) 2 and **Σ** is positive definite.
