Second-order Control of Complex Systems with Correlated Synthetic Data

The generation of synthetic data is an essential tool to study complex systems, allowing for example to test models of these in precisely controlled settings, or to parametrize simulation models when data is missing. This paper focuses on the generation of synthetic data with an emphasis on correlation structure. We introduce a new methodology to generate such correlated synthetic data. It is implemented in the field of socio-spatial systems, more precisely by coupling an urban growth model with a transportation network generation model. We also show the genericity of the method with an application on financial time-series. The simulation results show that the generation of correlated synthetic data for such systems is indeed feasible within a broad range of correlations, and suggest applications of such synthetic datasets.


Introduction
Developing methods to study complex systems, such as simulation models or data-mining techniques, often requires testbeds and benchmarks to ensure expected properties. The use of synthetic data, in the sense of statistical populations generated randomly under constraints of proximity of patterns to a studied system, is a widely used methodology tackling this issue. This approach is used in several disciplines related to complex systems such as therapeutic evaluation [1], territorial science [2,3], machine learning [4] or bio-informatics [5].
Generation of synthetic datasets can consist in data disaggregation by producing a microscopic population with fixed macroscopic properties [6]. The creation of synthetic populations for microsimulation models is a typical example where empirical statistical distributions are reproduced [7]. In data extensive contexts, several methods have been developed and improved for a better reproduction of margin distributions [8].
Synthetic datasets can also be generated at the same scale than the targeted real dataset, with a broad range of realism levels and corresponding constraints on the generated data [9]. For example, [10] show that some datamining techniques such as decision trees can be inverted to produce datasets capturing complex non-linear patterns.
The constraints of proximity to reality of synthetic dataset will depend on expected applications. They range for example from a strong statistical fit on given indicators, to weaker assumptions of similarity on aggregated patterns. In the case of systems where emergence plays a central role, a microscopic property does not directly imply given macroscopic patterns, and synthetic datasets may have to capture some of these. This approach therein becomes part of the complex systems simulation toolbox. Indeed, with the rise of new computational paradigms [11], data (simulated, measured or hybrid) shape our understanding of complex systems. Methodological tools for data-mining, modeling and simulation, including the generation of synthetic data, are therefore crucial to be developed.

Proposed approach
This literature review on different aspects of synthetic data generation unveils at least two gaps: (i) a lack of attention on controlling covariance structures when generating synthetic data; and (ii) an absence of such methods applied to the study of socio-spatial systems at aggregated scales. As spatio-temporal dependencies structures are essential in driving the dynamics of such systems [29,30], the combination of these two aspects appears as an unexplored research problem.
We propose in this paper to study the generation of correlated synthetic data, and more particularly in the case of socio-spatial systems. We introduce here a generic methodology taking into account the dependance structure for the generation of synthetic datasets, more precisely by controlling the average of correlation matrices. It is suited to be applied on cases where microscopic data is not available and system similarity is expected on aggregated indicators.
We investigate thus the question of how to generate correlated synthetic data at aggregated levels, where constraints on macroscopic indicators are fulfilled and correlation structure is controlled. We focus on this problem in the particular case of socio-spatial systems, but keep in mind the genericity of the approach.
Our contribution is twofold: (i) we implement a generation of spatial synthetic data for socio-spatial systems, which to the best of our knowledge has never been done in that context; (ii) the method introduced is generic, and we illustrate it with an application to financial time-series.
The rest of the paper is organized as follows. The generic method to generate correlated synthetic data is first formally described. We then apply it to a generative model of territorial configurations, composed by the sequential coupling of a reaction-diffusion model for population density with a road network generation model, and study the 2 produced correlation patterns. We also illustrate in a following section illustrate the genericity of our method by applying it to financial time-series, which are an other example of highly complex signals for which correlations are crucial.

Method Formalization
The domain-specific methods described above are too broad to be summarized within a same formalism. We therefore introduce here a generic and model-agnostic framework, focused on the control of correlations structures in synthetic data.
Let X I a multidimensional stochastic process (which can be indexed e.g. with time in the case of time-series, but also with space, or any other indexation). We assume to have a real dataset X = (X i,j ), which is interpreted as a set of realizations of the stochastic process. We propose to generate a statistical populationX =X i,j such that 1. A given criteria of proximity to data is verified, i.e. given a precision ε and some aggregated indicator f , we have 2. The level of correlation is controlled, i.e. given a matrix R representing the correlation structure (any symmetric matrix with coefficients in [−1, 1] and a unity diagonal), we have the estimated covariance matrix given bŷ where the standard deviation diagonal matrix Σ is estimated on the synthetic population.
The second requirement will generally be conditional to parameter values determining generation procedure, either generation models being simple or complex (R itself is a parameter). Formally, we can also understand synthetic processes as parametric familiesX i [ α].
We propose to apply the methodology on very different examples, both typical of complex systems: territorial systems and financial high-frequency time-series. We illustrate the flexibility of the method, and claim to help building interdisciplinary bridges by methodology transposition and reasoning analogy. In the first case, morphological calibration of a population density distribution model allows to respect real data proximity. Correlations of urban form with transportation network measures are empirically obtained by exploration of coupling with a network morphogenesis model. The control is in this case indirect and the feasible space of correlations is empirically determined. In the second case, proximity to data is the equality of signals at a fundamental frequency, to which higher frequency synthetic components with controlled correlations are superposed.

Correlated population density and road network
We now apply the method to territorial systems of human settlements, in the particular case here of population distribution in correlation with road network. In this application, simulation appears as a crucial step to implement the method.

Territorial configuration model
We propose in our case to generate territorial systems summarized in a simplified way as a spatial population density d( x) and a transportation network n( x). Correlations we aim to control are correlations between urban morphological measures and network measures. The question of interactions between territories and networks is already well-studied [31] but remains highly complex and difficult to quantify [32]. A dynamical modeling of implied processes should shed light on these interactions [33], and [34] has investigated the concept of co-evolution within such models. We develop here in that context a simple coupling (i.e. without any feedback loop) between a population density distribution model and a network morphogenesis model.

Density model
We use a model D similar to aggregation-diffusion models [35] to generate a discrete spatial distribution of population density. A generalization of the basic model is proposed in [36], providing a calibration on morphological objectives (entropy, hierarchy, spatial auto-correlation, mean distance) against real values computed on the set of 50 km sized grid extracted from European density grid [37]. We recall here rapidly the processes included in the model. An square grid of width W , initially empty, is represented by population (P i (t)) 1≤i≤W 2 . At each time step, until the total population reaches a fixed parameter P m , • total population is increased of a fixed number N G (growth rate), following a preferential attachment such that P[P i (t + 1) = P i (t) + 1|P (t + 1) = P (t) + 1] = (P i (t)/P (t)) α (P i (t)/P (t)) α (3) • a fraction β of population is diffused to four closest neighbors is operated n d times The two opposite processes of urban concentration and urban sprawl are captured by the model, what allows reproducing with a good precision a large number of existing morphologies regarding macroscopic urban form indicators.

Network model
On top of the population density model, we generate a planar transportation network with a model N at a similar scale. Several processes can be taken into account to simulate network growth [38]. Other model types could be used as well, such as biological self-generated networks [39], local network growth based on geometrical constraints optimization [40], or a more complex model based on multi-dimensional network percolation [41] which would allow the creation of loops for example. [38] generates networks in the frame of a modular architecture, in which the choice of the network generation heuristic can be adapted to a specific need (as e.g. proximity to real data, constraints on output indicators, variety of generated forms, etc.).
We choose here an heuristic based on spatial interaction potential breakdown, which corresponds in practice to a network answering to the strongest demand patterns. The algorithm assumes realistic thematic assumptions: a connected initial network and the creation of links based on spatial interactions.
The heuristic network generation procedure is the following: 1. A fixed number N c of centers that will be first nodes of the network are distributed given density distribution. Their spatial distribution follows a similar law than the aggregation process, i.e. the probability to be distributed in a given patch is (Pi/P ) α (Pi/P ) α . Population is then attributed according to Voronoi areas of centers, such that a center cumulates population of patches within its triangulation extent.
2. Centers are connected deterministically through a percolation between closest clusters: as long as the network is not fully connected, the two closest connected components, in the sense of the minimal euclidian distance between their vertices, are connected with the link realizing this distance. It yields a tree-shaped network at this stage.
3. The network is modified by adding links following a spatial interaction potential breaking. More precisely, a generalized gravity potential between two centers i and j is defined by is a weight to determine the role of populations, γ gives the shape of the hierarchy across population values, r g is a characteristic interaction distance and d 0 is a distance shape parameter. 4. A fixed number K · N L of potential new links is taken among the couples having greatest euclidian distance potential (K = 5 is fixed).
5. Among potential links, N L are effectively realized, that are the one with smallest rateṼ At this stage only the gap between euclidian and network distance is taken into account:Ṽ ij does indeed not depend on populations and is increasing with d N at constant d ij .
6. Planarity of the network is imposed by creating nodes at possible intersections formed by new links.
The nature and range of correlations produced by this model coupling, as a function of model parameters, are to be determined by simulation experiments.

Parameter space
The parameter space for the coupled model is constituted first by density generation parameters α D = (P m /N G , α, β, n d ).
We study for the sake of simplicity the rate between population and growth rate P m /N G instead of both varying, i.e. the number of steps needed to generate the distribution. These are completed by network generation parameters α N = (N C , k h , γ, r g , d 0 ). We write α = ( α D , α N ).

Indicators
Urban form and network structure are quantified by numerical indicators in order to modulate correlations between these. Morphology is defined as a vector M = (r, d, ε, a) giving spatial auto-correlation (Moran index), mean distance, entropy and hierarchy (see [42] for a precise definition of these indicators). Network measures G = (c, l, s, δ) are with network denoted (V, E) • Average centrality c defined as average betweeness-centrality (normalized in [0, 1]) on all links.
• Average path length l given by 1 • Average network speed [43] which corresponds to network performance compared to direct travel, defined as We study the cross-correlation matrix Cov M , G between morphology and network. We estimate it on a set of with the standard unbiased estimator, given by Eq. 5 below.ρ The covariance is estimated by Eq. 6, where variables are indexed by t over T realizations.
The variance is estimated by Eq. 7.V Null model In order to provide a minimal benchmark of our correlated data generation method, we also introduce a null model to control if the produced correlation are not intrinsic to the specification of indicators for example. The procedure to generate null configuration is the following: (i) generate a random population density, by randomly selecting a proportion r L between random pairs of nodes; (iv) planarize the resulting network by adding nodes at link intersections. In this model, population density and network are either totally independent, or linked through network node density only. We thus expect the corresponding correlation to behave as a baseline of how correlations between indicators behave when no particular care is given to including interaction processes.

Generating correlated synthetic data
The coupling of generative models is done both in a formal and operational way. We indeed loosely couple independent implementations. The OpenMOLE software [44] for model exploration offers a proper framework for this. Its modular workflow language allows to compose model tasks and integrate these into diverse numerical experiments. For the population density generation, we use the scala implementation provided by [36]. The network generation model is implemented in NetLogo [45], which offers a good compromise between performance and interactive model validation and exploration. The two models are coupled with a specific OpenMOLE script. Source code is available at https://github.com/JusteRaimbault/CityNetwork/tree/master/Models/Synthetic.

Results
The study of the density model alone is developed in [36]. It is in particular calibrated on European density grid data, on 50km width square areas with 500m resolution for which real indicator values have been computed on whole Europe. Furthermore, a grid exploration of model behavior yields feasible output space in reasonable parameters bounds (roughly α ∈ [0. 5,2], N G ∈ [500, 3000], P m ∈ [10 4 , 10 5 ], β ∈ [0, 0.2], n d ∈ {1, . . . , 4}). The reduction of indicators space to a two dimensional plan through a Principal Component Analysis (variance explained with two components 80%) allows to isolate a set of output points that covers reasonably precisely real point cloud. It confirms the ability of the model to reproduce morphologically the set of real configurations.
With a fixed population density, the conditional exploration of network generation model parameter space suggest a good flexibility on global indicators G, together with good convergence properties. In order to apply the synthetic data generation method in relation with the thematical question of interactions between networks and territories, the exploration has been oriented towards the study of cross-correlations.
Given the large relative dimension of the parameter space, an exhaustive grid exploration is not possible. We use a Latin Hypercube sampling procedure with bounds given above for α D and for α N , we take N C ∈ [50, 120], 4,20]. For number of model replications for each parameter point, less than 50 are enough to obtain confidence intervals at 95% on indicators of width less than standard deviations. For correlations a hundred give confidence intervals (obtained with Fisher method) of size around 0.4, we take thus n = 80 for experiments. The null model is simulated also with n = 80, for random and density-based node distributions, and r We show in Fig. 1 examples of generated territorial configurations. This visualization and some values of associated correlations already suggest that the method application yields a broad spectrum of generated correlation patterns. We obtain for example low density configurations, in aggregated or dispersed settings (top left, resp. bottom left panel), inducing very different types of networks. Similarly, areas with population centers which are closer to urban areas (Top right and bottom right panel), can also produce different network shapes. In the latest case, increasing the role of hierarchy through γ and k h leads from a negative correlation between average distance d and centrality c to a positive correlation. This corresponds to a transition from processes where population dispersal decrease centrality (redundant networks) to inverse processes (centralized networks).
Regarding the generation of correlated synthetic data in itself, several results presented in Fig. 2 are worth noting. First of all, the statistical distributions of correlation coefficients (histograms, top left panel of Fig. 2) between morphology and network indicators are not systematically simple and some are bimodal. For example, the correlation ρ[a, l] between hierarchy of the population distribution a and mean path length in the network l has a mode around 0 and an other around 0.6. This means that in a certain regime, these tend to decorrelate in average, while in an other regime they are strongly correlated. The latest correspond to configurations with a high Moran index and a high hierarchy, which means that more centralized urban configurations constrain the network path length through this correlation.
Second, still based on distributions in Fig. 2, but also on heatmaps for amplitude and maximal correlation (top right panel), we observe that it is possible to modulate up to a relatively high level of correlation for all indicators, since the maximal absolute correlation varies between 0.6 and 0.9. The amplitude of correlations ranges between 0.9 and 1.6, allowing thus a broad spectrum of values.
As the cross-correlation matrix is of dimension 16, we proceed to a principal component analysis on all generated correlation matrices (one matrix per row) to visualize the covered space in two dimensions. This PCA yields 38% of variance for the first component and 68% of cumulated variance for the second. We visualize the corresponding point cloud in the principal plan, with transformed confidence intervals (bottom left panel of Fig. 2) and with particular points (bottom right panel). The point cloud in the principal plan has a large extent but is not uniform: it is not possible to modulate in any direction any coefficient as they stay themselves correlated because of underlying  generation processes. A more refined study at higher orders (correlation of correlations) would be necessary to precisely understand degrees of freedom in the generation of correlations. However, the covered area remains broad and confirms a rather flexible output space for generated correlations. When comparing to the null model runs (black dots and error bars), we find as expected that null model correlations are around 0 (all confidence intervals covering the origin), and that a part of the generated point cloud falls in the same area. An other important part of points fall outside the range of the null model in a statistically significant way. These are the interesting points for a possible application of the synthetic dataset, and we show thus that the method produces non-trivial and significant correlation patterns. When evaluating the proximity of indicator values to real points (Equation 1 in the abstract description of the method), which is given by the color level in the bottom right panel of Fig. 2, we note that the points with the highest level of correlation are also the ones which are closest to real data (red points). The points which are the farthest from real configurations are the uncorrelated ones, which also coincide with the null model. This suggests that in the frame of model hypotheses, real configurations exhibit high correlations between network properties and urban form. [46] confirms this fact by studying real effective correlations.
Finally, some examples of configurations taken on particular points in the principal plan, highlighted in blue in Fig. 2 and described above (Fig. 1), show that similar population density profiles can yield very different correlation profiles. This confirms the flexibility of the method and the possibility to control on correlation structure.

Correlated financial time-series
We also apply the method to a totally different type of system, namely financial complex systems. Financial time-series are heterogeneous, multi-scalar and non-stationary [47]. Correlations are broadly explored in that field. For example, Random Matrix Theory allows to distinguish signal from noise for a correlation matrix computed for a large number of asset with low-frequency signals, typically with a time scale of a day [48]. Similarly, Complex Network Analysis on networks constructed from correlations introduced methods such as Minimal Spanning Tree [49] or more refined topologically constrained network generation methods [50]. These provide reconstructions of economic sectors structure. At high frequencies, the precise estimation of interdependence parameters assuming models for asset dynamics has been extensively studied from a theoretical point of view aimed at refinement of models and estimators [51]. Theoretical results must be tested on synthetic datasets as they ensure a control of most parameters in order to check that a predicted effect is indeed observable all things being otherwise equal. Empirical confirmation of the improvement of estimators is obtained on a synthetic dataset at a fixed correlation level.
We consider a network of assets (X i (t)) 1≤i≤N sampled at high-frequency (typically 1s). We use a multi-scalar framework (used e.g. in wavelet analysis approaches [52] or in multi-fractal signal processing [53]) to interpret observed signals as the superposition of components at different time scales. We thus write X i = ω X ω i . We denote by T ω i = ω ≤ω X ω i the filtered signal at a given frequency ω. A typical problem in the study of complex systems is the prediction of a trend at a given scale. It can be viewed as the identification of regularities and their distinction from components considered as random. For the sake of simplicity, we represent such a process as a trend prediction model at a given temporal scale ω 1 , formally an estimator M ω1 : (T ω1 i (t )) t <t →T i ω1 (t) which aims to minimize error on the real trend T ω1 i −T ω1 i . In the case of autoregressive multivariate estimators, the performance will depend among other parameters on respective correlations between assets. It is thus interesting to apply the method to the evaluation of performance as a function of correlation at different scales. We assume a Black-Scholes dynamic for assets [54], i.e. dX = σ · dW , with W Wiener process. Such a dynamic model allows an easy modulation of correlation levels.

Data generation
We can straightforward generateX i such that Var X ω1 i = t Σ · R · Σ (with Σ are estimated standard deviations and R is a fixed correlation matrix) and verifying X ω≤ω0 i =X ω≤ω0 i . This means in practice that the data proximity indicator is the identity of components at a lower frequency than a fundamental frequency ω 0 < ω 1 . We use therefore the simulation of Wiener processes with fixed correlation. Indeed, if dW 1 |= dW |= 1 (and σ 1 < σ 2 indicatively, assets being interchangeable), then is such that ρ(dW 1 , dW 2 ) = ρ 12 . Signals for other components can be constructed the same way by Gram orthonormalization. We isolate the component at the desired frequency ω 1 by filtering the signal, i.e. using signals constructed with Eq. 8 such thatX ω1 where F ω0 is a low-pass filter with cut-off frequency ω 0 . We reconstruct then the hybrid synthetic signals by taking The method is tested on an example with two assets from foreign exchange market (EUR/USD and EUR/GBP), on a six month period from June 2015 to November 2015. Data was obtained from http://www.histdata.com/. The data cleaning procedure, starting from original series sampled at a frequency around 1s, consists in a first step to the determination of the minimal common temporal range (missing sequences being ignored, by vertical translation of series, i.e. S(t) := S(t) · S(tn) S(tn−1) when t n−1 , t n are extremities of the "hole" and S(t) value of the asset, what is equivalent to keep the constraint to have returns at similar temporal steps between assets). We study then log-prices and log-returns [47], defined by X(t) := log S(t) S0 and ∆X(t) = X(t) − X(t − 1). Raw data are filtered at a maximal frequency ω m = 10min (which will be the maximal frequency for following treatments) for concerns of computational efficiency. As time-series are then sampled at 3 · ω m to avoid aliasing, a day of size 86400 for 1s sampling is reduced to a much smaller size of 432. We use a non-causal gaussian filter of total width ω. We fix the fundamental frequency ω 0 = 24h and we propose to construct synthetic data at frequencies ω 1 = 30min, 1h, 2h. We show in Fig. 3 an example of signal structure at the scales ω m and ω 1 = 30min, compared with the non-filtered raw signal.
It is crucial to consider the interference between ω 0 and ω 1 frequencies in the reconstructed signal: the correlation which is indeed estimated is Assuming to be in the reasonable limit σ 1 σ 0 (fundamental frequency small enough), that Cov ∆X ω1 i , ∆X ω j = 0 for all i, j, ω 1 > ω and that returns are centered at any scale, we can develop the previous expression to compute the correction on effective correlation due to interferences. We obtain at the first order the expression of effective correlation given by what corresponds to the correlation that we can effectively simulate in synthetic data.
Correlations are in practice estimated with a Pearson estimator, the covariances being corrected for bias, i.e. following Eq. 5, Eq. 6 and Eq. 7.
The generated synthetic data are then used to test a toy model. We propose in particular to investigate the predictive power of a very simple linear model. The tested predictive model M ω1 is a simple ARMA for which parameters p = 2, q = 0 are fixed (as we do not create lagged correlation, we do not expect large orders of auto-regression as these kind of processes have short memory for real data; furthermore smoothing is not necessary as data are already filtered). It is however applied in an adaptive way, in the sense that given a time window T W , we estimate for any t the model on [t − T W + 1, t] in order to predict signals at t + 1.
Experiments are implemented in the R language, using in particular the MTS [55] library for time-series models. Cleaned data and source code are available on an open git repository at https://github.com/JusteRaimbault/ SyntheticAsset. Fig. 4 shows the effective correlations computed on synthetic data. For standard parameter values (for example ω 0 = 24h, ω 1 = 2h and ρ = −0.5), we find ρ 0 0.71 et ε i 0.3 what yields |ρ e − ρ| 0.05. We observe a good agreement between observed ρ e and values predicted by Equation 11 in the interval ρ ∈ [−0.5, 0.5]. On the contrary, for larger absolute values, a deviation increasing with |ρ| and as ω 1 decreases : it confirms the intuition that when frequency decreases and becomes closer to ω 0 , interferences between the two components are not negligible anymore and invalidate independence assumptions for example.
Application to the study of a predictive model performance The predictive model described above is then applied to synthetic data, in order to study its average performance as a function of correlation between signals. Results for ω 1 = 1h, 1h30, 2h are shown in Fig. 5. The a priori counter-intuitive result of a maximal performance at vanishing correlation for one of the assets confirms the role of T , local polynomial smoothing to ease reading). It is interesting to note the U-shape for EUR/USD, due to interference between components at different scales. Correlation between simulated noises deteriorates predictive power. The study of lagged correlations (here ρ[∆X EURUSD (t), ∆X EURGBP (t − τ )]) on real data clarifies this phenomenon: the fourth graph show an asymmetry in curves at any scale compared to zero lag (τ = 0) what leads fundamental components to increase predictive power for the dollar, amelioration then perturbed by correlations between simulated components. Dashed lines show time steps (in equivalent τ units) used by the ARMA at each scale, what allows to read the corresponding lagged correlation on fundamental component.
13 synthetic data to better understand system mechanisms : the study of lagged correlations shows an asymmetry in the real data that we can understand at a daily scale as an increased influence of EUR/GBP on EUR/USD with a rough two hours lag. The existence of this lag allows a "good" prediction of EUR/USD thanks to fundamental component. This predictive power is perturbed by added noises in a way that increases with their correlation. The more noises correlated are, the more he model will take them into account and will make false predictions because of the Markovian character of simulated brownian (note that the model used has theoretically no predictive power at all on pure brownians).
This case study on a toy-model illustrates the relevance of using simulated synthetic data. Further developments can be directed towards the simulation of more realistic data (presence of consistent lagged correlation patterns, more realistic models than Black-Scholes) and apply it on more operational predictive models.

Contributions
We investigated in this paper the possibility of generating synthetic data at a macrosopic level with a controlled correlation structure. The generic method we introduce can be applied to any complex system, where the proximity to real data is measured on aggregated indicators. The method was designed more particularly for socio-spatial systems. We show in the case of transportation network and territories, by exploring a weak coupled model for population density and road network generation, that varying model parameters yield a broad output space of effective correlations. Two configurations with the same first order indicator values can capture very different underlying correlations. This means that future applications to the study of upstream models to the sensitivity of spatial initial configuration, such as the one done by [27] but in which correlation structure is controlled, should be made possible by our approach.
We postulate that the method can be also applied in other fields where similar constraints can be of interest. Indeed, in the context of financial data, considering data proximity indicators based on low-frequency components of signals, we showed how correlation can be controlled and even analytically predicted to a certain extent. Our work recalls thus the interest in generating hybrid data, and is differentiated from most work where only the microscopic level is taken into account.
As already mentioned, most of simulation models need an initial state generated artificially as soon as model parametrization is not done completely on real data. An advanced model sensitivity analysis implies a control on parameters for synthetic dataset generation, seen as model meta-parameters [27]. In the case of a statistical analysis of model outputs it provides a way to operate a second order statistical control.

Future work
Regarding the application to geographical data, the calibration of the network generation component at given density, on real data for transportation network, is a potential development. It would be relevant typically on road networks given the shape of generated networks, what should be possible using OpenStreetMap open data which have a reasonable quality for Europe [56]. This should theoretically allow to unveil parameter sets reproducing accurately existing configurations both for urban morphology and network shape. By attributing a synthetic dataset similar to a given real configuration, we would be able to compute a sort of intrinsic correlation proper to this configuration.
We studied in the second example stochastic processes in the sense of random time-series, whereas time did not have a role in the first case. We can suggest a strong coupling between the two model components (or the construction of an integrated model) and to observe indicators and correlations at different time steps during the generation. In dynamical spatial models the existence of lagged interdependences in space and time [29] is an important feature of complex dynamics. This would provide a better understanding of the link between static and dynamic correlations.
We were limited to the control of first and second moments of generated data, but we could imagine a theoretical generalization allowing the control of moments at any order. However, as shown by the geographical example, the difficulty of generation in a concrete complex case questions the possibility of higher orders control when keeping a consistent structure model and a reasonable number of parameters. The study of non-linear dependence structures as proposed in [57] is in an other perspective an interesting possible development.
We could also apply specific exploration algorithms to explore more exhaustively the feasible correlation space. Such an algorithm based on Novelty Search has been introduced by [58]. Coupling it with our method would allow establishing the full range of feasible correlations for a given generation model.

Conclusion
We proposed an abstract method to generate synthetic datasets in which correlation structure is controlled, but the empirical data required can be sparse or targeting macroscopic aggregated criteria. Its implementation in two very different fields shows its flexibility and the broad range of possible applications. More generally, it is crucial to favorise such practices of systematic validation of computational models by statistical analysis, in particular for agent-based models for which the question of validation remains an open issue.
Furthermore, our overall approach enters a particular epistemological frame. On the one hand it has a strong multidisciplinary aspect, and on the other hand the importance of empirical component through computational exploration methods make this approach typical of Complex Systems science [59]. The combination of empirical knowledge obtained from data mining, with knowledge obtained by modeling and simulation is generally central to the conception and exploration of multi-scalar heterogeneous models. The method and results presented here is an illustration of such an hybrid paradigm.