Multivariate Lomax (Pareto type II) distribution was introduced by Nayak (1987) as a joint probability distribution of several skewed positive random variables \(X_1, X_2, \cdots, X_k\). Its probability density function is given by
$$f(x_1, x_2, \dots, x_k) = \frac{[ \prod_{i=1}^{k} \theta_i] a(a+1) \cdots (a+k-1)}{(1+\sum_{i=1}^{k} \theta_i x_i)^{a+k}},$$
where \(x_i > 0, a>0, \theta_i>0, i=1,\dots,k\).
Cumulative distribution function \(F(x_1, \dots, x_k)\) is obtained by the following formula related to survival function \(\bar{F}(x_1, \dots, x_k)\) (Joe, 1997)
$$F(x_1, \dots, x_k) = 1 + \sum_{S \in \mathcal{S}} (-1)^{|S|} \bar{F}_S(x_j, j \in S),$$
where the survival function is given by
$$\bar{F}(x_1, \dots, x_k) = ( 1+\sum_{i=1}^{k} \theta_i x_i )^{-a}.$$
Equicoordinate quantile is obtained by solving the following equation for \(q\) through the built-in one dimension root finding function uniroot:
$$\int_{0}^{q} \cdots \int_{0}^{q} f(x_1, \cdots, x_k) dx_k \cdots dx_1 = p,$$
where \(p\) is the given cumulative probability.
Random numbers from multivariate Lomax distribution can be generated by simulating \(k\) independent exponential random variables having a common environment parameter following gamma distribution with shape parameter \(a\) and scale parameter \(1\); see Nayak (1987).