In this work we propose a new approach to build multivariate pure jump processes. We introduce linear and nonlinear dependence, without restrictions on marginal properties, by imposing a multi-factorial structure separately on both positive and negative jumps. Such a new approach provides higher flexibility in calibrating nonlinear dependence than in other comparable Lévy models in the literature. Using the notion of multivariate subordinator, this modeling approach can be applied to the class of univariate Lévy processes which can be written as the difference of two subordinators. A common example in the financial literature is the variance gamma process, which we extend to the multivariate (multi-factorial) case. The model is tractable and a straightforward multivariate simulation procedure is available. An empirical analysis documents an accurate multivariate fit of stock index returns in terms of both linear and nonlinear dependence. An example of multi-asset option pricing emphasizes the importance of the proposed multivariate approach.