Multivariate Pareto type I distribution was introduced by Mardia (1962) as a joint probability distribution of several nonnegative random variables \(X_1, \cdots, X_k\). Its probability density function is given by
$$f(x_1, \cdots, x_k) = \frac{[ \prod_{i=1}^{k} \theta_i] a(a+1) \cdots (a+k-1)}{(\sum_{i=1}^{k} \theta_i x_i - k + 1)^{a+k}},$$
where \(x_i > 1 / \theta_i, a >0, \theta_i>0, i=1,\cdots, 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, \cdots, x_k) = \left( \sum_{i=1}^{k} \theta_i x_i - k + 1 \right)^{-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 \(X_1, \cdots, X_k\) from Mardia's multivariate Pareto type I distribution can be generated through linear transformation of multivariate Lomax random variables \(Y_1, \cdots, Y_k\) by letting \(X_i = Y_i + 1/\theta_i, i = 1, \cdots, k\); see Nayak (1987).