Streamlining Simulation Experiments with Agent-Based Models in Demography

,


Introduction
. Traditionally, demography has been a highly empirical discipline of the social sciences, focusing predominantly on long-term changes of population-level processes (Burch ). Still, in the recent years, the field has witnessed an increasing dissatisfaction with the inadequate level of theoretical explanation o ered for demographic phenomena (Hobcra ) -a discussion echoing similar concerns raised in other areas of the social sciences, such as in sociology (Hedström & Ylikoski ). Some answers to these challenges can be o ered by the tools and techniques of social simulation, in particular, by agent-based modelling (Silverman ). This is especially important, as by using simulation-based methods, we can gain insights into the mechanisms of the underlying processes (the why of demographic phenomena), which cannot be examined by other tools, be it macro-level models of population dynamics, or survey-based studies. .
The suggestions to use agent-based models for studying demographic patterns, rather than just applying the ubiquitous hypothesis-driven approaches, have already gained some traction in the population science community over the last decade and a half. The use of agent-based computational models in demography dates back to the seminal book edited by Francesco Billari and Alexia Prskawetz (Billari & Prskawetz ). This research area has gained momentum especially over the past decade, with some of the most recent advances .
The relative youth of agent-based modeling within demography and the fact that it is not yet an accepted mainstream tool means that no one single approach to these tasks exists, and, given that the suitability of methods depends on the model and the research questions to be asked, is likely to remain this way. Early approaches were understandably rather ad hoc -for instance, in Billari et al. ( ), a model of marriage and social pressure, results are reported for a set of default parameters, and are also provided for a small number of additional scenarios.
. Similarly, the simulation in Hills & Todd ( ) models partnership formation as a search and match process where agents aim to find agents similar to themselves, but relaxed their criteria as to what constitutes a good match over time. Their experiments are aimed at examining the theoretical consequences of greater population heterogeneity and cultural diversity on the matching process by varying the relevant parameters one at a time. They are able to qualitatively match US marriage curves (although no detail as to how this was achieved is given), and they also recreate the divorce rate using the same parameter settings, as an attempt at model validation. .
Later examples are more systematic in their approaches to calibrate and validate the agent-based models by simulation experiments. Most o en, the parameter set which minimizes some distance to observed quantities is sought. For instance, the model in Aparicio Diaz et al. ( ) examines the e ect of social network pressure on transition to first birth in Austria, and replicates changes in the timing of first birth seen in Austria in the s. A metric, measuring the distance between simulated and real age-specific fertility rates, is used to assess model performance. A larger set of experiments are attempted in this work, with combinations of parameters evaluated over a grid, and a 'null' model in which the network e ects are turned o is also investigated. Presumably, the stated default parameters are those on the grid for which this distance is minimized. However, because of the practice of varying two parameters while keeping the others fixed, large areas of the parameter space are le unexplored. .
Other work (Fent et al. ) advances this agenda by examining the e ect of social networks on the success of family policies. The authors modeled preferred family size as dependent on opinions of social network members. This time, social network growth is endogenous and dependent on agent similarity, degree of relatedness, and number of shared network contacts. A 'full factorial' grid search over parameters leads to a total of ,   (Coale ), showing the dynamic relationship between son preference, social pressure relating to family size, and di usion of abortion technologies. The model is calibrated by fitting a regression meta-model to a Latin Hypercube sample of points, and then using these predictions to minimize the predicted root-mean-squared error of their model relative to observations.

.
Other papers have also looked at the behavior of simulations over a grid of points, including Klabunde ( ), which examined circular migration between Mexico and the USA, modeling agents as deciding whether to migrate using a discrete choice framework, i.e, the theory of planned behavior (Klabunde & Willekens ). Following Werker & Brenner ( ), some parameters of the model were fixed empirically using Mexican Migration Project data, whilst others were based on behavior rules, the parameters for which were found through a grid search of possible combinations, with a goodness of fit metric forming the criterion of choice. In contrast, and focusing on the phenomenon of assortative mating, Grow & Van Bavel ( ) develop an agent-based model examining how increased educational attainment amongst women has led to changes in the marriage market. This simulation is calibrated against empirical data using regression meta-models fit on central composite designs, and validated against a hold-back set of data from other countries. .
The above examples show that experimentation with agent-based models is important in addressing the sorts of questions typically of interest to demographers. However, conducting these experiments requires significant e orts. Furthermore, agent-based models of demographic processes are not always produced with reproducibility in mind, and simulation code and descriptions of the experiments conducted are o en not provided or incomplete, with notable exceptions including Grow & Van Bavel ( ) and Klabunde ( ). In part, this may be because of the di iculty in easily specifying and sharing sets of experiments. Provision of a simple and flexible way of specifying, sharing and running experiments with agent-based computational models of demographic processes would therefore further agent-based modeling within demography.
. The remainder of the paper is structured as follows. In Section , the needs for flexible support and for documenting simulation experiments with agent-based models are outlined, followed by an introduction of a Scala layer for specifying and executing simulation experiments (SESSL -Simulation Experiment Specification via a Scala Layer). In Section , the Modeling Language for Linked Lives (ML ) is introduced. The test model, a simple prototype model targeting migration networks, is described in Section . Section describes the use of ML and SESSL to define and execute experiments on this simple migration model, used here for the transparency of presentation. In Section , further experiments are described and executed for a more complex and realistic model of social care. The final section provides a concluding discussion. All models and experiments in the paper are available online (https://git.informatik.uni-rostock.de/mosi/jasss ).

Simulation Experiments with SESSL
. The Simulation Experiment Specification on a Scala Layer (SESSL) (Ewald & Uhrmacher ) is a simulationsystem-agnostic layer between user and simulation system -in other words, a layer that can work with di erent simulation environments. It has been developed to support the compact specification and execution of a wide variety of simulation experiments allowing users to add code that is executed during the experiment. Thus, proficient users can directly add missing features "on-the-fly". As a domain-specific language embedded in the programming language Scala (https://www.scala-lang.org), SESSL o ers many extension points where usersupplied functionality can be injected. .

The current features of SESSL include
• simulation-based multi-objective optimization, via a binding for Opt J (Lukasiewycz et al. ), • statistical/simulation-based model-checking (Sen et al. ), • di erent experimental design methods, including full factorial design, Latin Hypercube sampling (Mckay et al. ), or central composite design, • replication criteria, such as maximal relative width of a confidence interval of simulation output, • structured output of simulation results to CSV files, and • bindings for di erent simulation packages. .
Although SESSL specifications are valid, executable Scala code, the resulting experiment specifications do not resemble typical program code. The salient aspects of simulation experiments, such as the model input parameter configuration or the observation of model outputs, are specified in a declarative style. But vanilla Scala functions, for example to post-process the observed output, can still seamlessly be integrated if needed. Larger fragments of code can be packaged in a reusable binding. SESSL bindings are particularly useful to integrate third-party tools and make their features available in experiment specifications. Bindings to external tools are mostly translating elements of a SESSL specification to invocations of external tools. Thus, they are slim and easily implemented. When setting up an experiment, users can choose from the available bindings to enrich their experiment with features and connect it to one or more external tools. SESSL relies on Apache Maven (https://maven.apache.org) for the automatic resolution of dependencies and download of so ware artifacts for bindings and third-party so ware.
. val prey = observe(count(" * /Prey")); val predator = observe(count(" * /Predator")) observeAt(range(0.0, 0.1, 11)) set("a" <~0.014, "b" <~0.6, "d" <~0.7, "k" ). A sequential probability ratio test (Wald ) with specific statistical parameters is chosen in line as the statistical test procedure to check the property. Lines -finally specify that the result of the check shall be output a er the experiment finishes. example, a binding to a simulation package may provide traits for parallel execution, observation, or report generation; a binding for statistical model-checking may provide a trait that checks the model output against a behavior specification given in temporal logics. Statistical model checking is an established technique that makes behavioral properties -in other approaches, also called 'hypotheses' (Lorig et al. ) -that shall be checked formally, explicit, and provides means for automatically testing these properties (Agha & Palmskog ).
. Figure illustrates the working of SESSL with an example of statistical model-checking executed on a Lotka-Volterra 'predator and prey' model of interdependent dynamics of two populations, which has been defined in ML-Rules, a domain-specific modeling language for modeling and simulation in a di erent domain, namely systems biology (Maus et al. ). To perform a statistical model-checking experiment, the property to check must be explicitly defined. The experiment in Figure checks for "coexisting state" (a commonly checked property of Lotka-Volterra models (Droz & Pȩkalski )) that states that both prey and predators will not become extinct (see Peng et al. for details).
. Traits and bindings make the specification of simulation experiments in SESSL flexible and agile. Due to the automatic management of so ware artifacts, executing and repeating SESSL experiments is straightforward across machines and platforms. Based on SESSL, simulation experiments such as the statistical model checking shown above can even be automatically generated and executed to test whether specific behavior properties of models do still hold when models are extended (Peng et al. ) or composed (Peng et al. ). At the same time, SESSL specifications are succinct and readable, allowing easy sharing and communication of experiment specifications. This makes SESSL an excellent tool for experiment specification in domains that require diverse experiment design methods and non-standard approaches, such as agent-based computational demography. ) is a domain specific modeling language specifically designed for agent-based demographic models with dynamic social networks. While most agentbased models are executed in a step-wise manner, and commonly used modeling environments such as Net-Logo (Wilensky ) pursue such a discrete-time approach, time in ML is continuous. This allows to capture the temporal behavior more realistically, and avoids many of the inherent problems of time-stepped simulation(Willekens ). The ML simulator so ware is implemented in Java (https://www.java.com) and distributed using Apache Maven. The ML Sandbox, an editor for ML models, is available for download in the ML repository (https://git.informatik.uni-rostock.de/mosi/ml ).

Key features of ML .
The primary entities of all ML models are agents, representing individual persons as well as higher-level demographic actors such as households. To distinguish the di erent kinds of entities that are modeled as agents, each agent has a type that determines its properties and behavior. Thus, each ML model begins with the definition of agent types. Figure (lines to ) shows such an agent type definition. There, an agent type Person is defined. Every agent of this type has two attributes: their sex, which can have one of the two values "female" and "male", and their income, which is a real number. As aging of agents is an important part of many demographic models, e.g., with respect to fertility and mortality, all agents have an attribute age in addition to their explicitly defined attributes. During the simulation, the age of agents is initialized when an agent is "born", and updated automatically as time passes. Finally, all agents are either alive or dead. Agents are alive when they are created, and might die during the simulation. Dead agents do no longer act independently, but other agents might still be influenced by them. .
To facilitate the modeling of dynamically evolving social networks, these networks are an explicit part of an ML model. Agents in ML are interconnected with dynamic links, which model their relationships. For example, parents and their children might be connected with a link. The definition of the kinds of links agents of a certain type might have is the second part of every ML model, as Figure (line ) demonstrates for the example of parents and children. It reads as follows: Every person is linked to zero or more other persons, their children (read from le to right). Every person is also linked to exactly two other persons, their parents (read from right to le ). Links are always directed (meaning the parent-direction is di erent from the children-direction), bidirectional (there has to be a parent-direction for the child-direction), and dynamic (links may change during the simulation, e.g., parents might get more children). .
The behavior of agents is governed by stochastic rules. For each agent type, a set of behavioral rules must be defined, and those rules apply to all agents of that type that are currently alive. Rules are defined as guard-ratee ect triples, stating who is subject to the rule (guard), when the rule is applied (rate), and what happens when the rule is applied (e ect). Figure lines to shows a simple mortality rule as an example. The rule applies to all male persons (line ). Its timing is governed by the rate expression ( ). Here, the rate is given by the map maleMortalityRate, a special data structure ML uses to hold time-series data. In this case the map might hold the age-dependent force of mortality for di erent cohorts as it might have been estimated from data. The  Figure ). The experiment object, defined in the SESSL binding for ML , is invoked by a vanilla Scala main function. Subsequently, an experiment object in the ML simulation package is created, which starts simulation runs and observes them. A er n simulation runs are finished, the n trajectories that were observed during the runs are checked against the property defined in the experiment. The checking is handled by the StatisticalModelChecking trait that lives in the SESSL verification package. The process is repeated until enough simulation runs have been observed to decide whether the property is satisfied or not. rate expression gives the (possibly time dependent) parameter of an exponential distribution. Finally, the rule ends with the e ect (line ). In the case of this mortality rule, the agent simply dies. During the simulation for each pair of agent and rule a waiting time is drawn from the according waiting time distribution, and the agent-rule pair with minimal waiting time is executed. In addition to stochastically timed rules, which are the most commonly used, ML supports a range of deterministic ways to schedule an event, e.g., to schedule events when an agent reaches a certain age or to schedule periodic events.

The SESSL-ML binding .
To exploit the features of SESSL for experimentation with ML models, we implemented a SESSL binding for ML . The binding covers the basic features of SESSL experiments, such as choosing a model file and a simulation algorithm, setting simulation stop conditions and replication numbers, configuring parallel execution, and the observations made during the simulation. Similar as in other SESSL bindings for simulation packages, all specifications are translated to API (application programming interface; the interface through which other so ware communicates with a piece of so ware) calls of the ML package. A sketch of the call hierarchy of an experiment is shown in Figure . In Figure , an experiment identical to the one shown in Figure is specified using SESSL-ML instead of ML-Rules. .
In addition to traits for the generic experiment aspects, we also implemented experiment aspects in the binding that are specific to ML . To realize the binding for ML we adopted a new concept of model parameterization.
Previous SESSL bindings executed models with scalar model inputs. However, as ML models are aimed at describing demographic phenomena, many model parameters are time series data or distributions, which ML represents as parameter maps. For example, for individuals the risk of dying depends on their age (see Figure  line ). The pattern of age-dependent event rates is typical for applications in demography. In the return migration model the age distribution at migration is realized as a parameter map. Consequently, we implemented a trait ParameterMaps that allows reading in parameter maps from .csv files. val prey = observe(agentCount("Prey")); val predator = observe(agentCount("Predator")) observeAt(range(0.0, 0.1, 11)) set("a" <~0.014, "b" <~0.6, "d" <~0.7, "k" <~0.002) initializeWith ( In lines , , , and , drop-in replacements were available to switch from using ML-Rules to ML . Line had to be added to construct the initial state for the model.

Illustration: A Model of Return Migration
. As a first case study, we present a novel but prototypical and relatively simple agent-based model of return migration. To illustrate the main features of our experimental approach in a succinct and transparent way, we use this simple model rather than a more complex one, which is discussed next. Primarily the model consists of four processes: the inflow of migrants into the system; the formation of a social network between the migrants; the employment and earnings of the migrants, depending on the social network; and the decision of migrants to return to their country of origin. The model is kept simple on purpose, to capture only the main features of the migration processes of interest, enabling to illustrate the process of the experimental design with a flexible simulation environment o ered by ML and SESSL. Following the approach of Cio i-Revilla (Cio i-Revilla ) the model, as presented here, is merely a relatively simple initial model version M 0 . The results of the experiments we conduct with this initial model suggest a way forward for refinement of the model in an iterative process, by which through identifying the model parameters of key importance, and possibly collecting additional information allowing for identifying them better, we would be able to design subsequent model versions (Courgeau et al. ). For more realistic applications, a comprehensive review of recent advancements in agent-based modeling of migration is provided by Klabunde & Willekens ( ). .
The primary entities of our models are migrant agents. Migrants are characterized by a number of attributes that influence their economic success and their behavior regarding their return migration decision (Figure ). Migration streams predominantly consist of younger migrants, and so the simulated ages of new migrants are drawn from a distribution reflecting this. The distribution was constructed by fitting a -parameter model based on that developed by Rogers ( ) to UK migration data from the O ice of National Statistics; it is provided to the simulation as an input file. Migrant age is significant in the model as once a predetermined retirement age is reached ( for the runs described here), agents are removed from the simulation. In a later model version other processes, e.g., employment and income, might also depend on age. .
Simulated migrants are also endowed with a randomly determined 'skill level' which determines their earnings potential. This reflects migrants di ering ability to succeed in the host country, which is o en incorporated into similar economic models (Borjas & Bratsberg ; Dustmann et al. ). Skill levels are drawn from a normal distribution centered on and truncated from below at . The variance of this distribution is a variable param- .
On initialization, there are no migrants present in the system. Migrants enter the system stochastically with exponentially distributed inter-arrival times. The rate depends upon the di erence between the current average earnings and migrant's target earnings, in line with economic theory, which suggests, that migrant flow depends on earnings di erentials (Harris & Todaro ). In the ML implementation this process is governed by a stochastic rule (Figure ). When the rule is executed, a new migrant agent is created (line ) and its attributes are initialized according to the respective distributions (lines -and -), and the new migrant is introduced into the network (lines -). .
Network ties are determined on entry to the system, with an initial 'contact' chosen at random for each new migrant from those agents already in the system, and further ties made from amongst this contact's neighbors. This reflect the tendency for migrants' moves to be facilitated by someone they already know in the host country (Liu ). When an agent leaves the system (either because they have satisfied one of the return migration conditions described below, or they have reached the retirement age), their network ties are removed. .
Migrants also participate in labor markets a er migration. The model of the labor market used is relatively crude. Employment opportunities come along at a fixed rate (although the model allows for unemployment to increase as more migrants enter the system). Wages are assumed to respond to labor supply, so that more migrants in the system results in lower wages. More specifically, labor unit prices respond quadratically to supply. Of course, the precise relationship between migrant wages and number of immigrants depends on the many factors including local labor market conditions and production technology; however, for an initial choice, diminishing returns from labor would appear to be an appropriate assumption, and similar assumptions have been made elsewhere (Epstein ). However, skill-based heterogeneity in earnings amongst migrants is included as a multiplier to this basic relationship as shown in ( ). Future models could include explicit price setting and bargaining based upon the information and preferences of both employer(s) and migrants.
In Equation ( ), y i is the earnings of the i th migrant; N is the total number of migrants present in the system; and s i is the skill of the migrant in question.
. Total income or utility is a ected not only by earned income but also through social capital associated with the extent of a migrant's network. In line with work by Portes ( ), the social capital associated with migrant networks could realize itself in a number of forms, including access to credit and housing, sharing of information regarding the host country, loaning of durable goods, or provision of childcare. Although these are generally provided in non-monetary forms, it is presumed here that they act to increase net income by reducing expenditure on certain goods; the benefit it provides is fungible (ibid). In empirical work, it is o en assumed that the relevant quantity determining network e ects is the total number of migrants present in a given location (Epstein , e.g.). However, in this case, the specific number of network ties each migrant has determines their where g is the size of the agent's network. .
Agents are presumed to follow one of two strategies: either they intend to settle permanently in their new country of residence; or they wish to save until they reach some target value which will allow them to live more comfortable in their country of origin. This distinction between temporary and permanent migrants in incoming migration streams has support in empirical migration literature (Constant & Massey ; Bijwaard , e.g). These strategies define the rules used to describe when agents migrate back to their home country.

.
In the case of permanent 'settler' migrants, migrants return if they earn less than some minimum value C (Figure , line ). Their counterparts, 'saving' migrants, instead return once they reach some target value T (Figure , line ), which is drawn from a normal distribution (although an income target for 'savers' is also defined, but at a lower value, as they are willing to accept lower current consumption in favor of later reward upon return to their country of origin). If the minimum income conditions are met, migrants do not leave immediately; instead, they remain in the country for a short period z in case the situation changes, before finally returning. .
The relative frequency of 'savers' in the immigrant distribution is determined by parameter θ. Di erent proportions of the two migrant types may be expected to produce di erent behavior in the overall system, both because of their individual characteristics and due to potential emergent properties resulting from their interdependence in the network. However, in order to assess the relative influence of the di erent elements of the model on the proportion of 'savers', a more thorough and systematic experimentation is required. Figure : The four rules that describe how agents leave the system: when they retire (line ), when their fail to reach their income target for too long (line ), when they are unemployed for too long (line ), and when they have fulfilled their savings target (line ). As the e ect of all four events is the same, i.e., the removal of the agent, the actual e ect was implemented as a procedure that is called in each of the rules (line -).

Experiments with the Migration Model
. Thanks to the binding between SESSL and ML , all features of SESSL are now ready to use for agent-based models specified in ML , including the migration model presented in the previous section. The variety of experiment designs, of optimization methods and statistical model checking, form an expressive portfolio for analyzing the behavior of agent-based models in demography. However, so far one class of experimental methods has not been supported: meta-modeling, which increasingly receives attention, also in the area of agent-based demography (Bijak et al. ; Silverman et al. a). Meta-models, or response surfaces, or surrogate models, provide means to approximate the input-output function that is realized by the underlying simulation model. They are in fact statistical models of the underlying simulation models. As such, they are used in many application fields (see, for example, the Managing Uncertainty in Complex Models project, http://www.mucm.ac.uk). The meta-models are particularly useful, whenever the calculation of individual simulation runs, and thus parameter scans or optimizations, imply significant computational e ort and expense. .
To demonstrate how SESSL can be extended to this new class of analytical methods for agent-based models, two examples of meta-models are derived based on two experimental designs. Firstly, a regression-based metamodel is fit to a parameter scan relying on central composite design (Kleijnen ), and secondly, a Gaussian Process emulator is used to analyze a Latin hypercube sample obtained in a second set of experiments (Kennedy & O'Hagan ). In both cases, the links between the di erent designs and analysis models are natural, as central composite designs assume relationships between variables which can be captured through the low-order polynomials used in regression meta-models. In contrast, Gaussian processes are semi-parametric, meaning that no particular functional form for the relationship between inputs and outputs is assumed, and thus the more flexible Latin hypercube design is preferable. More complex stochastic models might even ask for additional methods to deal with the high-dimensionality of agent-based models e ectively, e.g., Salemi et al. ( ). .
With the two meta-models we employ two di erent strategies for SESSL to be integrated into an experimentation workflow. For the linear regression meta-model, the whole process of fitting the meta-model is part of the experiment specified in SESSL. This includes the execution of the central composite design as well as the actual fitting of the regression model. However, for the Gaussian process meta-model we demonstrate a different approach. While SESSL is used to conduct the simulation runs according to the Latin Hypercube design, the meta-model is fitted with an external tool. To this end, the SESSL experiment produces CSV files containing the data produced by the simulation. Those can then be processed externally, in our case using the R statistical programming language (R Development Core Team ). This approach has the advantage of added flexibility, as the simulationist is not restricted to the methods provided by SESSL itself. However, the advantages SESSL o ers for documentation and reproducibility are partly lost. While the documentation and reproducibility is still ensured by SESSL for the experimental setup, it is not for the fitting of the regression model.

Example : Central composite design and regression meta-models .
Experimental design allows to gain information about the relationships between experimental inputs and outputs by controlling the location of input points. A central composite design provides for the identification of second-order polynomial relationships in the data. It consists of the combination of three sets of points (factorial, axial and central), as detailed in Figure (a). The first such set is generally a two-level factorial (grid-based) design, where each parameter is assigned a high and low value, and all combinations of these parameter values are identified. Secondly, a set of axial points are identified, whereby all parameters but one are set to their central value, and the remaining parameter is given extreme values (beyond those used for the factorial design). Finally, a central point is defined, where all parameters are set to their central value. Generally, repeated runs are used at the central point to gain information about the inherent stochasticity of the system under study (in the case of stochastic simulation, this is the part of variability in output caused by pseudo-random number generation, and unrelated to parametric changes). Although Figure displays only two dimensions, the design generalizes naturally to higher dimensions.
. Parameters controlling the degree of heterogeneity in 'skill' τ ; the proportion of savers and settlers in the population θ; the number of initial network ties m and the extent to which network ties e ect overall earnings γ were varied for this case study. Initially, the aim was to understand the e ect of these parameters on the steady-state number of total migrants in the system. Although time is continuous in the ML simulations, avoiding problems of simultaneity, in this case counts of migrants were observed every year for ease of interpretation. Figure  provides the SESSL code needed to define the experiment described, with the design itself specified at line , and the regression meta-model at line .

.
The estimates of standardized coe icients are given in Figure . Considering the main (linear) e ects, it can be seen that θ and γ have the largest e ects on the average number of migrants in the system. The quadratic and interaction e ects underline the importance of these variables, while the heterogeneity in skill τ and number of network ties m appear to have little e ect.

Example : Latin hypercube sample and Gaussian process emulator
. A similar experiment was defined varying the same parameters and monitoring the same output variable, but utilizing a di erent design and meta-model; a Latin hypercube design was used to fit a Gaussian Process emulator to the data. A Latin hypercube sample is a space-filling design which divides each dimension into regularly val regr = fitLinearModel(result)("theta", "gamma", "tau", "m", "theta" * "theta", "theta" * "gamma", "theta" * "tau", "gamma" * "gamma", "gamma" * "tau", "tau" * "tau", "m" * "gamma", "m" * "theta", "m" * "tau", "m" * "m")(migrants) println ( Figure . sized bins equal to the number of desired observations, and chooses points pseudo-randomly in such a way that there is exactly one observation in each bin for every dimension. In contrast to a factorial or grid-based design, it does not repeat values of any single parameter (Urban & Fricker ), as can be seen from the tick marks along the axis in Figure . As with the Central Composite design, the Latin Hypercube design generalizes easily to higher dimensions and can be specified in SESSL with a simple command, as shown in Figure . .
Gaussian process emulators work on the assumption that values close to each other in input space will also be close to each other in terms of output. More specifically, outputs are considered to be jointly multivariate normal, with the degree of correlation determined by a covariance function operating on inputs. Full details can be found in Kennedy & O'Hagan ( ). Gaussian process are flexible in that they don't make assumptions about the possible shapes of the output as a function of inputs, with the exception that some degree of smoothness is assumed, to be estimated from the data. .
Once fitted, the emulator can be used to produce predictive probability distributions for output values at any combination of input points, incorporating both the uncertainty about the mean value of the simulation at that point, and uncertainty about the degree of underlying simulation stochasticity (Kennedy & O'Hagan ). In this example, we have used the implementation of the method presented in Hilton & Bijak ( ). Predictions at mean values across the range of each covariate are displayed in Figure . The results are broadly consistent with the results obtained from the meta-model, although local movements in the mean surface are captured, as can be seen to some extent in the predictions over the range of the settler-saver-proportion θ in the Figure. .
A more comprehensive view can be gained by conducting a global sensitivity analysis of our initial model M 0 ; this examines the proportion of the total variance that can be attributed to each variable and combination of variables, averaged over prior distribution on these variables. The mathematical details can be found in Kennedy & O'Hagan ( ), and a summary in Silverman et al. ( a). In our example, as can be seen from the Figure , in model M 0 the majority of the variance is associated with θ and γ, in agreement with results from the meta-model coe icients displayed in Figure , although a significant part is also related to m.

.
These results indicate that the next phase of the model-building process should concentrate on obtaining more information about θ and γ, rather than focus on τ or m. This can be possibly achieved by conducting a surveybased or experimental empirical study pertaining to these parameters, including the propensity for saving (θ)

A Model of Social Care
. To demonstrate that our approach is also applicable to larger scale models, we executed a further simulation experiment on a linked lives model of social care by Noble et al. ( ), that aims at capturing the e ects of the aging population and changing family structures on the cost of state-funded social care in the United Kingdom (UK). The model is concerned with the demographic processes that a ect the supply and demand of social care. In the model, the population of the UK is represented by agents with a scaling factor of 1 : 10, 000, i.e. one agent represents 10, 000 people in the real world. For all agents life-course transitions that are important for social care, e.g., fertility, mortality, partnership formation and internal migration, are simulated. Agents are born as dependent children of their parents. With age , they reach adulthood. At this point they enter the workforce and become a taxpayer. They may marry another agent, move to a di erent part of the UK, and get children of their own. When they reach the retirement age, they retire and leave the workforce. Anytime during their life, their health status might degrade, which increases their need for social care. Finally, they die, which removes them from the population. For the purpose of internal migration, the model divides the UK into a grid of towns, with consist of multiple houses. Agents may migrate to a di erent house in the same town or to a di erent town multiple times in their life. We have reimplemented the model, originally implemented in Python, in ML (Warnke et al. ). .
The central output variable of the model is the total cost of social care per taxpayer. Every agent needs a certain amount of social care, measured in hours per week, depending on their care need category -which may include zero. A part of this care need can be fulfilled informally by relatives, the rest needs to be paid for by the state. Individuals who need no or only little care themselves can provide some hours for informal care each week. However, they will not deliver care to anybody, but only to persons living in the same household and to their parents, as long as the parents live in the same town. The care they can deliver is distributed among the persons they are willing to care for. Therefore, household structures and mobility a ect the amount of informal social care that is actually delivered. .
Given the high complexity of this model, in this section we present only the key high-level features and findings related to the conducted experiments. In particular, we show a simulation-based optimization experiment, a further experimental method supported by SESSL. The term simulation-based optimization refers to the search for those values of model input parameters that minimize or maximize an objective, which is a function of the simulation output (Amaran et al. ). Simulation-based optimization is commonly employed for model calibration, by minimizing the di erence between the simulated model output and the calibration target. While meta-models can also be used for optimization purposes (Barton & Meckesheimer ) they cannot yield satisfying results for all optimization problems, e.g., with discontinuous objective functions, many input variables or discrete inputs (Nguyen et al. ). There is no 'one-size-fits-all' optimization method, and a plethora of algorithms exist for di erent classes of optimization problems (see Amaran et al. ( ) for an overview). In SESSL simulation-based optimization is realized through a binding to the Opt J (Lukasiewycz et al. ) library of optimization algorithms, allowing the user to choose between many di erent methods.

.
In the experiment (Figure ), we want to minimize the cost of social care per taxpayer by varying the retirement age. Therefore, we embedded the experiment specification (line -) in a minimize-block, which is provided by the SESSL-Opt J binding. As the calculation of the cost of social care is very complex (see Section . ), it can not be realized with SESSL-ML . However, we can exploit the extendibility of SESSL, by implementing a custom observation directly using the observation-interface of the ML simulator. We have wrapped that in the trait HealthCareCostObservation (lines and ). Similarly, we have implemented a custom initial state builder (line ), which creates the initial state exactly as in the Python implementation of the model. In both cases, SESSL allows for making direct calls to the more powerful, but di icult to use, API of the ML simulator when necessary, while keeping the easy to use features of SESSL when possible.
. In line the model parameter for the age of retirement is set to the value chosen by the optimization algorithm.
In lines -the objective function of the optimization is calculated. The optimization itself is then parameterized a er the experiment specification, where a range for the retirement age parameter and an optimization algorithm are chosen. The optimization yields a minimum of the cost of social care by setting the retirement age to about , with a weekly cost per taxpayer of £ . . This is in line with the findings of Silverman et al. ( b) that care cost per taxpayer actually increases when the retirement age is increased beyond , as the increased number of taxpayers will be o set by the reduction of available informal care in households with elderly members.

Discussion
. The review about benefits and limitations of the ODD protocol (Overview, Design concepts, and Details), which was published in to standardize the description of agent-based models (Grimm et al. ), recommends to enhance ODD by a separate section on simulation experiments. Information about experiments done with a model supports assessing the range for which a model might be valid, interpreting and questioning published simulation results, which is essential in valuing and reusing models. In addition, compact, declarative specifications of simulation experiments facilitate the generation and adaptation of simulation experiments in general, and with agent-based models in particular. They contribute to a more e icient and systematic engineering of simulation experiments (Teran-Somohano et al. ). .

SESSL (Simulation Experiment Specification via a Scala
Layer) is a domain-specific language for supporting the specification and execution of simulation experiments. It presents a simulation-system-agnostic layer between simulation system and user and supports a wide variety of simulation experiments, including specific experiment designs, parameter scans, simulation-based optimization, and statistical model checking. Due to the new presented binding between SESSL and ML (Modeling Language for Linked Lives), the entire methodological portfolio provided in SESSL can now be used for experimenting with agent-based models specified in ML and for documenting these simulation experiments in a non-ambiguous manner. .
As an internal domain specific language, SESSL facilitates adaptations and extensions for particular application domains. For demographic models, we realized and integrated a new experimental design method, i.e., central composite design, and two meta-modeling methods, i.e., regression meta-models and a Gaussian process emulator. We combined them with further methods SESSL already supported to analyze a simple agent-based return migration model and a more complex model of social care. The resulting simulation experiments can be succinctly described, executed, and adapted in SESSL and show that in this model the parameters for migrant strategies distribution and network ties are crucial for the system's behavior. The binding between SESSL and ML demonstrates the potential of domain specific languages for conducting and documenting simulation experiments with agent-based models in a convenient, replicable, and at the same time, flexible manner. ) and Burch ( ), these developments can substantially enhance the explanatory power of agent-based models, especially in such a strongly empirical discipline as demography. On the one hand, they can help connect the micro and macro levels of analysis Billari ( ), and provide a template for designing theoretical microfoundations explaining populationlevel phenomena in a rigorous and systematic way. At the same time, the proposed approach remains within the spirit of the Manifesto of Computational Social Science (Conte et al. ), arguing for a greater empirical relevance of social simulation models, for example with respect to modelling the human decision processes involved (Klabunde et al. ), and also helping unravel the mechanisms underpinning the social processes under study. The approach proposed in this paper provides a blueprint -as well as e icient so ware toolsfor addressing these research challenges in a systematic fashion through an iterative process of design, execution, analysis and documentation of computer experiments. Further developments in this area could include integrating a full statistical analysis, which would additionally enable the estimation of the free model parameters given the available data and assessing the uncertainty of the variables of interest. The ultimate goal of this process is to enhance the analytical possibilities o ered by social simulations by embedding them within the wider framework of experimental design -and by bringing in this way the empirical and theoretical aspects of demographic enquiries closer together.