Relatively recent techniques for categorical simulations are based on multi-point statistical approaches where a training image is used to derive complex spatial relationships using patterns. However, simulated geological realizations are driven by the training image utilized, while the spatial statistics of the hard data is ignored. This paper presents a data-driven high-order simulation approach based upon the approximation of high-order spatial indicator moments. The high-order spatial statistics are expressed as functions of spatial distances similar to variogram models for two-point methods. It is shown that the higher-order statistics are connected with lower-orders via boundary conditions. Using an advanced recursive B spline approximation algorithm, the high-order statistics are reconstructed from hard data. Finally, conditional distribution is calculated using Bayes rule and random values are simulated sequentially for all unsampled grid nodes. The main advantages of the proposed technique are its ability to simulate without a training image, which reproduces the high-order statistics of hard data, and to adopt the complexity of the model to the information available in the hard data. The approach is tested with a synthetic dataset and compared to a conventional second-order method, sisim, in terms of cross-correlations and high-order spatial~statistics.