Geostatistical simulation techniques are used to quantify uncertainty of spatial attributes of interest describing mineral deposits, petroleum reservoirs, hydrogeological horizons, environmental contaminants and so on. The majority of existing methods consider second-order spatial statistics and Gaussian processes, while the more advanced multiple point based simulation approaches are algorithmic and do not consistently account for the high-order spatial relations in data. Recently, simulation techniques for complex and non-Gaussian, spatially distributed variables have been developed, based on high-order spatial cumulants, and make no assumptions on data distribution or require data transformations. In this paper, the previous developments are extended and a new approach for the joint simulation of multiple correlated variables using high-order spatial statistics is proposed. The technique is based on a new algorithm described here for the decorrelation of correlated variables into factors, using the so-termed diagonal domination condition of high-order cumulants. The decorrelated factors are then simulated using high-order simulation and back-transformed into the initial correlated variables. The decorrelation using diagonal domination of high-order statistics is tested with a dataset from a multi-element iron ore deposit and then a fully known two-dimensional dataset with two correlated variables is used to demonstrate the practical intricacies of the proposed method.