3.1 Estimation

Take a simulation model depending on \(K\) parameters \(\theta = \left[ \theta_1,\theta_2,\dots,\theta_K \right]\) to output \(M\) simulated summary statistics \(S(\theta) = \left[ S_1(\theta), S_2(\theta),\dots,S_M(\theta) \right]\). Given real-world observed summary statistics \(S^*\) we want to find the parameter vector \(\theta^*\) that generated them.
Our method proceeds as follows:

1. Repeatedly run the model each time supplying it a random vector \(\hat \theta\)
2. Collect for each simulation its output as summary statistics \(S(\hat \theta)\)
3. Run \(K\) separate regularized regressions, \(r_i\), one for each parameter. The parameter is the dependent variable and all the candidate summary statistics are independent variables: \[ \left\{\begin{matrix} \theta_1 = r_1(S_1,S_2,\dots,S_M)\\ \theta_2 = r_2(S_1,S_2,\dots,S_M)\\ \vdots \\ \theta_n = r_n(S_1,S_2,\dots,S_M) \end{matrix}\right. \tag{3.1} \]
4. Plug in the “real” summary statistics \(S^*\) in each regression to predict the “real” parameters \(\theta^*\)

In step 1 and 2 we use the model to generate a data-set: we repeatedly input random parameters \(\hat \theta\) and we collect the output \(\hat S_1,\dots,\hat S_M\) (see Table 3.1 for an example). In step 3 we take the data we generated and use it to train regressions (one for each model input) where we want to predict each model input \(\theta\) looking only at the model output \(S_1,\dots,S_M\). The regressions map observed outputs back to the inputs that generated them.
Only in step 4 we use the real summary statistics (which we have ignored so far) to plug in the regressions we built in step 3. The prediction the regressions make when we feed in the real summary statistics \(S^*\) are the estimated parameters \(\theta^*\).

The key implementation detail is to use regularized regressions. We used elastic nets (Friedman, Hastie, and Tibshirani 2010) because they are linear and they combine the penalties of LASSO and ridge regressions. LASSO penalties automatically drop summary statistics that are uninformative. Ridge penalties lower the weight of summary statistics that are strongly correlated to each other.
We can then start with a large number of summary statistics (anything we think may matter and that the model can reproduce) and let the regression select only the useful ones. The regressions could also be made non-linear and non-parametric but this proved not necessary for our examples.

Compare this to simulated minimum distance. There, the mapping between model outputs \(S(\cdot)\) and inputs \(\theta\) is defined by the minimization of the distance function: \[ \theta^* = \arg_{\theta}\min (S(\theta)-S^*) W (S(\theta)-S^*) \tag{3.2} \] When using a simulated minimum distance approach, both the summary statistics \(S(\cdot)\) and their weights \(W\) need to be chosen in advance. Then one can estimates \(\theta^*\) by explicitly minimizing the distance. In effect the \(\arg \min\) operator maps model output \(S(\cdot)\) to model input \(\theta\).
Instead, in our approach, the regularized regression not only selects and weights summary statistics but, more importantly, also maps them back to model input. Hence, the regression bypasses (takes the place of) the usual minimization.

3.2 Model Selection

We observe summary statistics \(S^*\) from the real world and we would like to know which candidate model \(\{m_1,\dots,m_n\}\) generated them.
We can solve this problem by training a classifier against simulated data.

1. Repeatedly run each model
2. Collect for each simulation its generated summary statistics \(S(m_i)\)
3. Build a classifier predicting the model by looking at summary statistics and train it on the data-set just produced: \[ i \sim g(S_1,S_2,\dots,S_M) \tag{3.3} \]
4. Plug in “real” summary statistics \(S^*\) in the classifier to predict which model generated it

The classifier family \(g(\cdot)\) can also be non-linear and non-parametric. We used a multinomial regression with elastic net regularization(Friedman, Hastie, and Tibshirani 2010); its advantages is that the regularization selects summary statistics and the output is a classic logit formula that is easily interpretable.

Model selection cannot be done on the same data-set used for estimation. If the parameters of a model need to be estimated, the estimation must be performed on on either a separate training data-set or the inner loop of a nested cross-validation(Hastie, Tibshirani, and Friedman 2009).

4.1 Fit Lines

4.1.1 Regression-based methods fit straight lines as well as ABC or Simulated Minimum Distance

In this section we present a simplified one-dimensional parametrization problem and show that regression methods perform as well as ABC and simulated minimum distance.

We observe 10 summary statistics (in this case these are just direct observations) \(S^*=(S^*_0,\dots,S^*_9)\) and we assume they were generated by the model \(S_i=\beta i + \epsilon\) where \(\epsilon \sim \mathcal{N}(0,1)\). Assuming the model is correct, we want to estimate the \(\beta^*\) parameter that generated the observed \(S^*\). Because the model is a straight line, we could solve the estimation by least squares. Pretend however this model is a computational black box where the connection between it and drawing straight lines is not apparent.

We run 1,000 simulations with \(\beta \sim U[0,2]\) (the range is arbitrary and could represent either a prior or the feasible interval), collecting the 10 simulated summary statistics \(\hat S(\beta_i)\). We then train a linear elastic net regression of the form: \[ \beta = a + \sum_{i=0}^9 b_i \hat S_i \tag{4.1} \] Table 4.1 shows the \(b\) coefficients of the regression. \(S_0\) contains no information about \(\beta\) (since it is generated as \(S_0 = \epsilon\)) and the regression correctly ignores it by dropping \(b_0\) (setting the coefficient to 0). The elastic net weights larger summary statistics more because the effect of \(\beta\) is larger there compared to the \(\epsilon\) noise.

Table 4.1: Coefficients generated by fitting an elastic net against a training sample of size 1000. Notice that \(b_0\) has been automatically dropped since in this case \(S_0\) contains no information about \(\beta\)

Term	Estimate
(Intercept)	0.0297673
\(b_1\)	0.0010005
\(b_2\)	0.0065579
\(b_3\)	0.0056334
\(b_4\)	0.0109605
\(b_5\)	0.0134782
\(b_6\)	0.0190075
\(b_7\)	0.0203925
\(b_8\)	0.0238694
\(b_9\)	0.0281099

We test this regression by making it predict, out-of-sample, the \(\beta\) of 1000 more simulations looking only at their summary statistics. Figure 4.1 and table 4.2 compare the estimation quality with those achieved by the standard rejection ABC and MCMC ABC (following Marjoram et al. 2003; Wegmann, Leuenberger, and Excoffier 2009; both methods implemented in Jabot et al. 2015). All achieve similar error rates.

Figure 4.1: A side by side comparison between the ‘real’ and discovered \(\beta\) as generated by either the regression method or two standard ABC algorithms: rejection sampling or MCMC. The x-axis represents the real \(\beta\), the y-axis its estimate. Each point represents the estimation against a different set of summary statistics. The red line is the 45 degree line: the closer the points are to it the better the estimation. All methods perform equally well.

Table 4.2: All three methods achieve on average similar estimation errors, although the regression-based method is much quicker

Estimation	Average RMSE	Runtime	Runtime(scaled)
Regression Method	0.0532724	0.001515357 mins	1.000
Rejection ABC	0.0823950	0.6868809 mins	453.280
MCMC ABC	0.0835164	4.166347 mins	2749.417

We can also compare the regression against the standard simulated minimum distance approach. Simulated minimum distance search for the parameters that minimize: \[ \beta^* = \arg_{\beta}\min (S(\beta)-S^*) W (S(\beta)-S^*) \tag{3.2} \] The steeper its curvature, the easier it is to minimize numerically.
In figure 4.2 we compare the curvature defined by \(W=I\) and \(W= \left(\text{diag}(\Sigma)\right)^{-1}\) (where \(I\) is the identity matrix and \(\Sigma\) is the covariance matrix of summary statistics) against the curvature implied by our regression \(r\): \(|r(S^*)-r(S)|\) (that is, the distance between the estimation made given the real summary statistics, \(r(S^*)\), and the estimation made at other simulated summary statistics \(r(S)\)). This problem is simple enough that the usual weight matrices produce distance functions that are easy to numerically minimize.

Figure 4.2: Curvature of the fitness function as implied by the regression versus weighted Euclidean distance between summary statistics. We take the 1000 testing simulation runs, we pick one in the middle (\(\beta \approx 1\)) as the ‘real’ one and we compute the summary statistics distances between all the other simulations to it. The red line represents the real \(\beta\), each dot is the distance to the real \(\beta\) for that simulation run and its associated \(\beta\).

We showed here that regression-based methods perform as well as ABC and simulated minimum distance in a simple one-dimensional problem.

4.1.2 Regression-based methods fit broken lines better than ABC or simulated minimum distance

In this section we present another simplified one-dimensional parametrization but show that regression methods outperform both ABC and simulated minimum distance because of better selection of summary statistics.

We observe 10 summary statistics \(S=(S_0,\dots,S_9)\) but we assume they were generated by the “broken line” model: \[ S_i=\left\{\begin{matrix} \epsilon & i < 5\\ \beta i + \epsilon & i\geq5 \end{matrix}\right. \tag{4.2} \]
where \(\epsilon \sim \mathcal{N}(0,1)\).
We want to find the \(\beta\) that generated the observed summary statistics. Notice that \(S_0,\dots,S_4\) provide no information on \(\beta\).

As before we run the model 1000 times to train an elastic net of the form: \[ \beta = a + \sum_{i=0}^9 b_i S_i \tag{4.1} \] Table 4.3 shows the coefficients found. The regression correctly drops \(b_0,\dots,b_4\) and weighs the remaining summary statistics more the higher they are.

We run the model another 1000 times to test the regression ability to estimate their \(\beta\). As shown in Figure 4.3 and Table 4.4 the regression outperforms both ABC alternatives. This is due to the elastic net ability to prune unimportant summary statistics.

Table 4.3: Coefficient estimates after 1000 training observations; notice that all the uninformative summary statistics have been automatically dropped.

Term	Estimate
(Intercept)	0.0360036
\(b_5\)	0.0189285
\(b_6\)	0.0202486
\(b_7\)	0.0222020
\(b_8\)	0.0288608
\(b_9\)	0.0273448

Figure 4.3: A side by side comparison between the ‘real’ and discovered \(\beta\) as generated by either the regression method or two standard ABC algorithms: rejection sampling or MCMC. The x-axis represents the real \(\beta\), the y-axis its estimate. Each point represents the estimation against a different set of summary statistics. The red line is the 45 degree line: the closer the points are to it the better the estimation. Notice that rejection sampling (ABC) is inaccurate when \(\beta\) is very small (near 0) or very large (near 2) while MCMC has a larger average error than regression-based methods

Table 4.4: Regression errors achieve lower estimation errors in a fraction of the time. This is because they automatically ignore all summary statistics that do not inform \(\beta\)

Estimation	Average RMSE	Runtime	Runtime(scaled)
Regression Method	0.0579736	0.0007997513 mins	1.0000
Rejection ABC	0.1278042	0.692666 mins	866.1018
MCMC ABC	0.1343737	4.470746 mins	5590.1706

In simulated minimum distance methods, choosing the wrong \(W\) can compound the problem of poor summary statistics selection. As shown in Figure 4.4 the default choice \(W= \left(\text{diag}(\Sigma)\right)^{-1}\) detrimentally affects the curvature of the distance function, by adding noise. This is because summary statistics \(S_0,\dots,S_4\), in spite of containing no useful information, are the values with the lowest variance and are therefore given larger weights.

Figure 4.4: Curvature of the fitness function as implied by the regression versus weighted Euclidean distance between summary statistics. We take the 1000 training simulation runs, we pick one in the middle (\(\beta \approx 1\)) as the ‘real’ one and we compute the summary statistics distances between all the other simulations to it. The red line represents the real \(\beta\), each dot is the distance to the real \(\beta\) for that simulation run and its associated \(\beta\). For this model the default choice of weighing matrix (the inverse of the variance diagonal) makes numerical minimization harder.

We showed here how, even in a simple one dimensional problem, the ability to prune uninformative summary statistics will result in better estimation (as shown by RMSE in table 4.4) .

4.1.3 Model selection through regression-based methods is informative and accurate

In this section we train a classifier to choose which of the two models we described above (‘straight-’ and ‘broken-line’) generated the observed data. We show that the classifier does a better job than choosing the model that minimizes prediction errors.

We observe 10 summary statistics \(S^*=(S^*_0,\dots,S^*_9)\) and we would like to know if they were generated by the “straight-line” model: \[ S_i = \beta i + \epsilon \tag{4.3} \] or the “broken-line” model: \[ S_i=\left\{\begin{matrix} \epsilon & i < 5\\ \beta i + \epsilon & i\geq5 \end{matrix}\right. \tag{4.2} \]
where \(\epsilon \sim N(0,1)\) and \(i \in [0,9]\). Assume we have already estimated both models and \(\beta=1\).

We run each model 1000 times and then train a logit classifier: \[ \text{Pr}(\text{broken line model}) = \frac{1}{1 + \exp(-[a +\sum b_i S_i])} \tag{4.4} \] Notice that \(S_0\) and \(S_5,\dots,S_9\) are generated by the same rule in both models: both draw a straight line through those points. Their values will then not be useful for model selection. We should focus instead exclusively on \(S_1,\dots,S_4\) since that is the only area where one model draws a straight line and the other does not.
The regularized regression discovers this, by dropping \(b_0\) and \(b_5,\dots,b_9\) as Table 4.5 shows.

Table 4.5: Coefficient estimates after 2000 training observations; notice how only the low summary statistics have a coefficient associated to them. This is because they are the only part of the line that is ‘broken’: \(S_5,\dots,S_9\) are generated the same way for both models and are therefore of no use for differentiating the two models.

Term	Estimate
(Intercept)	-11.1578231
\(b_1\)	0.5733553
\(b_2\)	1.1016989
\(b_3\)	2.0715734
\(b_4\)	1.8777375

Figure 4.5 tabulates the success rate from 2000 out-of-sample simulations by the classifier and compares it with picking models by choosing the one that minimizes the distance between summary statistics. The problem is simple enough that choosing \(W=I\) achieves almost perfect success rate; using the more common variance-based weights does not.

Figure 4.5: A comparison between training a classifier to do model selection versus choosing the model whose distance to summary statistics is smaller. All weighing schemes except \(W=I\) underperform compared to the classifier. Model selection performed on 2000 testing runs (not used for training the classifier)

4.2 Fit RBC

4.2.1 Fit against cross-correlation matrix

In this section we parametrize a simple macro-economic model by looking at the time series it simulates. We show that we can use a large number of summary statistics efficiently and that we can diagnose estimation strengths and weaknesses the same way we diagnose regression outputs.

We observe 150 steps² of 5 quarterly time series: \(Y,r,I,C,L\). Assuming we know they come from a standard real business cycle (RBC) model (the default RBC implementation of the gEcon package in R, see Klima, Podemski, and Retkiewicz-Wijtiwiak 2018, see also Appendix A), we estimate the 6 parameters (\(\beta,\delta,\eta,\mu,\phi,\sigma\)) that generated them.

Summarize the time series observed with (a) the \(t=-5,\dots,+5\) cross-correlation vectors of \(Y\) against \(r,I,C,L\), (b) the lower-triangular covariance matrix of \(Y,r,I,C,L\), (c) the squares of all cross-correlations and covariances, (d) the pair-wise products between all cross-correlations and covariances. This totals to 2486 summary statistics: too many for most ABC approaches.
This problem is however still solvable by regression. We train 6 separate regressions, one for each parameter, as follows: \[ \left\{\begin{matrix} \beta &= \sum a_{\beta,i} S_i + \sum b_{\beta,i} S_{i}^2 + \sum c_{\beta,i,j} S_i S_j \\ \gamma &= \sum a_{\gamma,i} S_i + \sum b_{\gamma,i} S_{i}^2 + \sum c_{\gamma,i,j} S_i S_j \\ &\vdots \\ \sigma &= \sum a_{\sigma,i} S_i + \sum b_{\sigma,i} S_{i}^2 + \sum c_{\sigma,i,j} S_i S_j \\ \end{matrix}\right. \tag{4.5} \] Where \(a,b,c\) are elastic net coefficients.

We train the regressions over 2000 RBC runs then produce 2000 more to test its estimation. Figure 4.6 and Table 4.6 show that all the parameters are estimated correctly but \(\eta\) and \(\delta\) less precisely than the others. We judge the quality of the estimation through the predictivity coefficient (see Salle and Yıldızoğlu 2014; also known as modelling efficiency as in Stow et al. 2009) which measures the improvement of using an estimate \(\hat \beta\) to discover the real \(\beta^*\) as opposed to just using the hindsight average \(\bar \beta\): \[ \text{Predictivity} = \frac{ \sum \left( \beta^*_i - \bar \beta \right)^2 - \sum \left( \beta^*_i - \hat \beta_i \right)^2 }{\sum \left( \beta^*_i - \bar \beta \right)^2} \tag{4.6} \] Intuitively, this is just an extension of the definition to \(R^2\) to out-of-sample predictions. A predictivity of 1 implies perfect estimation, 0 implies that averaging would have performed as well as the estimate.

Figure 4.6: Estimated versus real parameters for 2000 RBC runs in the testing data-set (not used for training the regressions). Each dot is a pair real-estimate parameter from a single RBC run. Red line is the 45 degree line: the closer the dots are to it the better estimated that parameter is. The parameters \(\delta\) and \(\eta\) are not as well estimated as the other 4

Table 4.6: Predictivity for each variable when estimating parameters of the RBC model looking at cross-correlations and covariance; \(\delta\) and \(\eta\) are less precisely estimated than the other parameters. All results obtained on a testing data-set (not used for training) of 2000 RBC runs. Bias is defined as average distance between real and estimated parameter: \(\theta^*-\theta\), RMSE is the root mean square error: \(\sqrt{\sum \frac{\left(\theta^*-\theta\right)^2}{2000}}\). Variable bounds are the interval 20% below and above the default RBC parameters in the original gEcon implementation.

Variable	Variable Bounds	Average Bias	Average RMSE	Predictivity
\(\beta\)	[0.891,0.999]	0.0003105	0.0000351	0.9635677
\(\delta\)	[0.225,0.275]	-0.0000426	0.0000010	0.4997108
\(\eta\)	[1.8,2.2]	-0.0009051	0.0062672	0.5168588
\(\mu\)	[0.27,0.33]	0.0001315	0.0000774	0.7419450
\(\sigma\)	[0.01,0.03]	-0.0000391	0.0000071	0.7941670
\(\phi\)	[0.855,0.999]	-0.0001568	0.0001203	0.9265618

We showed here that we can estimate the parameters of a simple RBC model by looking at their cross-correlations. However this method is better served by looking at multiple auxiliary models as we show in the next section.

4.2.2 Fit against auxiliary fits

In this section we parametrize the same simple macro-economic model by looking at the same time series it simulates but using a more diverse array of summary statistics. We show how that improves the estimation of the model.

Again, we observe 150 steps of 5 quarterly time series: \(Y,r,I,C,L\). Assuming we know they come from a standard RBC model, we want to estimate the 6 parameters (\(\beta,\delta,\eta,\mu,\phi,\sigma\)) that generated them.

In this section we use a larger variety of summary statistics. Even though they derive from auxiliary models that convey much of the same information, we are able to use them efficiently for estimation. We summarize output by (i) coefficients of regressing \(Y\) on \(Y_{t-1},I_{t},I_{t-1}\), (ii) coefficients of regressing \(Y\) on \(Y_{t-1},C_{t},C_{t-1}\), (iii) coefficients of regressing \(Y\) on \(Y_{t-1},r_{t},r_{t-1}\), (iv) coefficients of regressing \(Y\) on \(Y_{t-1},L_{t},L_{t-1}\), (v) coefficients of regressing \(Y\) on \(C,r\) (vi) coefficients of fitting AR(5) on \(Y\), (vii) the (lower triangular) covariance matrix of \(Y,I,C,r,L\).
In total there are 48 summary statistics to which we add all their squares and cross-products.

Again we train 6 regressions against 2000 RBC observations and test their estimation against 2000 more.
Figure 4.7 and Table 4.7 shows that all the parameters are estimated correctly and with higher predictivity coefficient. In particular \(\delta\) and \(\eta\) are better predicted than by looking at cross-correlations.

Figure 4.7: Estimated versus real parameters for 2000 RBC runs in the testing data-set (not used for training the regressions). Each dot is a pair of real-estimate parameters from a single RBC run. The red line is the 45 degree line: the closer the dots are to it the better estimated that parameter is. All parameters are well estimated

Table 4.7: Predictivity for each variable when estimating parameters of the RBC model looking at cross-correlations and covariance; All parameters are well estimated. All results obtained on a testing data-set (not used for training) of 2000 RBC runs. Bias is defined as average distance between real and estimated parameter: \(\theta^*-\theta\), RMSE is the root mean square error: \(\sqrt{\sum \frac{\left(\theta^*-\theta\right)^2}{2000}}\). Variable bounds are the interval 20% below and above the default RBC parameters in the original gEcon implementation.

Variable	Variable Bounds	Average Bias	Average RMSE	Predictivity
\(\beta\)	[0.891,0.999]	0.0000154	0.0000023	0.9975601
\(\delta\)	[0.225,0.275]	-0.0000308	0.0000003	0.8285072
\(\eta\)	[1.8,2.2]	0.0001992	0.0012552	0.9040992
\(\mu\)	[0.27,0.33]	0.0000126	0.0000357	0.8793255
\(\sigma\)	[0.01,0.03]	-0.0000638	0.0000061	0.8100463
\(\phi\)	[0.855,0.999]	0.0002740	0.0000244	0.9862208

4.3 Heuristic selection in the face of input uncertainty

In this section we select between 10 variations of an agent-based model where some inputs and parameters are unknowable. By training a classifier we are able to discover which summary statistics are robust to the unknown parameters and leverage this to outperform standard model selection by prediction error.

Assume all fishermen in the world can be described by one of 10 decision-making algorithms (heuristics). We observe 100 fishermen for five years and we want to identify which heuristic they are using. The heuristic selection is however complicated by not knowing basic biological facts about the environment the fishermen exploit.
This is a common scenario in fisheries: we can monitor boats precisely through electronic logbooks and satellites but more than 80% of global catch comes from unmeasured fish stocks(Costello et al. 2012).
Because heuristics are adaptive, the same heuristic will generate very different summary statistics when facing different biological constraints. We show here that we can still identify heuristics even without knowing how many fish there are in the sea.

We can treat each heuristic as a separate model and train a classifier that looks at summary statistics to predict which heuristic generated them. The critical step is to train the classifier on example runs where the uncertain inputs (the unknown biological constraints) are randomized. This way the classifier automatically learns which summary statistics are robust to biological changes and look only at those when selecting the right heuristic.

We use POSEIDON(Bailey et al. 2018), a fishery agent-based model to simulate the fishery. We observe the fishery for five years and we collect 17 summary statistics: (i) aggregate five year averages on landings, effort (hours spent fishing), distance travelled from port, number of trips per year and average hours at sea per trip, (ii) a discretized(3 by 3 matrix) heatmap of geographical effort over the simulated ocean, (iii) the coefficients and \(R^2\) of a discrete choice model fitting area fished as a function of distance from port and habit (an integer describing how many times that particular fisher has visited that area before in the past 365 days ).
These summary statistics are a typical information set we might expect to obtain in a developed fishery.

Within each simulation, all fishers use the same heuristic. There are 10 possible candidates,all described more in depth in a separate paper³:

• Random: each fisher each trip picks a random cell to trawl
• Gravitational Search: agents copy one another as if planets attracted by gravity where mass is equal to profits made in the previous trip (see Rashedi, Nezamabadi-pour, and Saryazdi 2009)
• \(\epsilon\)-greedy bandit: agents have a fixed 20% chance of exploring a new random area otherwise exploiting(fishing) the currently most profitable discovered location. Two parametrizations of this heuristic are available, depending on whether fishers discretize the map in a 3 by 3 matrix or a 9 by 9 one.
• “Perfect” agent: agents knows perfectly the profitability of each area \(\Pi(\cdot)\) and chooses to fish in area \(i\) by SOFTMAX: \[ \text{Probability}(i) = \frac{e^{\Pi(i)}}{\sum_j e^{\Pi(j)}} \tag{4.7} \]
• Social Annealing: agents always return to the same location unless they are making less than fishery’s average profits(across the whole fishery), in which case they try a new area at random
• Explore-Exploit-Imitate: like \(\epsilon\)-greedy agents, they have a fixed probability of exploring. When not exploring each fisher checks the profits of two random “friends” and copies their location if either are doing better. Three parametrizations of this heuristics are possible: exploration rate at 80% or 20% with exploration range of 5 cells, or 20% exploration rate with exploration range of 20 cells.
• Nearest Neighbour: agents keep track of the profitability of all the areas they have visited through a nearest neighbour regression and always fish the peak area predicted by the regression.

All runs are simulations of the single-species “baseline” fishery described in the ESM of Bailey et al. (2018) (see also Appendix B).
We run the model 10,000 times in total, 1000 times for each heuristic. We also randomize in each simulation three biological parameters: (i) carrying capacity of the each cell map, (ii) diffusion speed of fish, (iii) the Malthusian growth parameter of the stock.

We train a multi-logit classifier on the simulation outputs as shown in section 3.2: \[ \text{Pr}(\text{Heuristic }i) = \frac{e^{b_i' S}}{1+\sum_j e^{b_j' S}} \tag{4.8} \] where \(b_1,\dots,b_{10}\) are vectors of elastic net coefficients.

We build the testing data-set by running the model 10,000 more times. We compare the quality of model selection against selection by weighted prediction error.
First, in Figure 4.8 the prediction error is computed by re-running simulations with the correct random biological parameters. In spite of having a knowledge advantage over the classifier, selection by prediction error performs worse. In Figure 4.9 we do not assume to know the random biological parameters. Prediction error (which now tries to minimize the distance between simulations which differ in both heuristics and biological parameters) performs much worse.

Figure 4.8: Comparison of success rate in correctly estimating which heuristics generated a given summary statistics \(S\). The classifier success rate has solid borders. In this set of simulations model selection by minimum prediction error is advantaged by knowing the correct biological parameters while the classifier is a single global model trained over all combinations of biological parameters.

Figure 4.9: Comparison of success rate in correctly estimating which heuristics generated a given summary statistics \(S\). The classifier success rate has solid borders. In this set of simulations model selection by minimum prediction error has no additional knowledge about the biological parameters.

Model selection by prediction error is hard if we are not certain about all model inputs. However it is trivial to train a classifier by feeding it observations from different parameter configurations such that only the summary statistics “robust” to such parameter noise are used for model selection.

5.1 Advantages of a regression-based approach

Regression are commonly-used and well-understood. We leverage that knowledge for parameter estimation.
By looking at out-of-sample prediction we can immediately diagnose parameter identification problems: the lower the accuracy, the weaker the identification (compare this to diagnosis of numerical minimization outputs as in Canova and Sala 2009) . Moreover, if the regression is linear, its coefficients show which summary statistic \(S_i\) inform which model parameter \(\theta_i\).

A secondary advantage of building a statistical model is the ease with which additional input uncertainty can be incorporated as shown in section 4.3. If we are uncertain about an input and cannot estimate it directly (because of lack of data or under-identification), it is straightforward to simulate for multiple values of the input, train a statistical model and then check whether it is possible to estimate reliably the other parameters. In a simulated minimum distance setting we would have to commit to a specific guess of the uncertain input before running a minimization or minimize multiple times for each guess (without any assurance that the weighing matrix ought to be kept constant between guesses).

5.2 Weaknesses

This method has two limitations. First, users cannot weight summary statistics by importance. As a consequence, the method cannot be directed to replicate a specific summary statistic. Imagine a financial model of the housing market which we may judge primarily by its ability of matching real house prices. When estimating its parameters the regression may however ignore house prices entirely and focus on another seemingly minor summary statistic (say, total number of bathrooms). While that may be an interesting and useful way of parametrizing the model, it does not guarantee a good match with observed house prices.
Nevertheless, in this case one option would be to transition to a weighted regression and weigh more the observations whose generated summary statistic \(S_i\) were close to the value of \(S_i^*\) we are primarily interested in replicating. As such, the model would be forced to include component deemed important, even if not the most efficient in terms of explanatory power. Comparing this fit with the unweighted one may also be instructive.

The second limitation is that the each parameter is estimated separately. Because each parameter is given a separate regression, this method may underestimate underlying connections (dependencies) between two or more parameters when such connections are not reflected in changes to summary statistics. This issue is mentioned in the ABC literature by Beaumont (2010) although it is also claimed that in practice it does not seem to be an issue. It may however be more appropriate to use a simultaneous system estimation routine(Henningsen and Hamann 2007) or a multi-task learning algorithm(Evgeniou and Pontil 2004).

5.3 Possible Extensions

Our method contains a fundamental inefficiency: we build a global statistical model to use only once when plugging in the real summary statistics \(S^*\). Given that the real summary statistics \(S^*\) are known in advance, a direct search ought to be more efficient.
Our method makes up for this by solving multiple problems at once as it implicitly selects summary statistics and weights them while estimating the model parameters.

In machine learning, transduction substitutes regression when we are interested in a prediction for a known \(S^*\) only, without the need to build a global statistical model(see Chapter 10 of Cherkassky and Mulier 2007; Beaumont, Zhang, and Balding 2002 weights simulations depending on the distance to known \(S^*\) in essentially the same approach). Switching to transduction would however remove the familiarity advantage of using regressions as well as the ease to test its predictive power with randomly generated summary statistics.

Another avenue for improvement is to use non-linear, non-parametric regressions. Blum and Francois (2010) used neural networks, for example. Non-parametrics are more flexible but they may be more computationally demanding and their final output harder to interpret.

It may also be possible to use regressions as an intermediate step in a more traditional simulated minimum distance problem. We would still build a regression for each parameter but then using them as the distance function to minimize by standard numerical methods.

5.4 Conclusion

We presented here a general method to fit agent-based and simulation models. We perform indirect inference using prediction rather than minimization, and, by using regularized regressions, we avoid three major problems of estimation: selecting the most valuable of the summary statistics, defining the distance function, and achieving reliable numerical minimization. By substituting regression with classification we can further extend this approach to model selection. This approach is relatively easy to understand and apply, and at its core uses statistical tools already familiar to social scientists. The method scales well with the complexity of the underlying problem, and as such may find broad application.