Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
110,596 characters · 32 sections · 88 citation commands
Black-box Bayesian inference for economic agent-based models
Simulation models are increasingly employed across the sciences. Their primary advantage is the flexibility they afford modellers: models may be specified mechanistically and at the microscopic level without concern for analytic tractability of the full system behaviour. This flexibility is embraced in the agent-based modelling paradigm wherein individual agents can be described by realistic and heterogeneous decision and interaction rules. Examining the macroscopic properties of the system is then done simply through simulation, freeing the modeller from the burden of reasoning about the global behaviour of the system themselves and facilitating the study of radically more complex systems and emergent phenomena than is otherwise possible.
In general, \glspl{abm} take as input some parameter vector $\boldsymbol{\theta}$, and return a (possibly multivariate) time-series $\mathbf{x}$. Typically, \glspl{abm} are stochastic simulators, with the property that running the simulator repeatedly with a fixed parameter $\boldsymbol{\theta}$ will produce different outputs with each run. The probability (density) with which the simulator produces a particular output $\mathbf{x}$, conditional on an input parameter vector $\boldsymbol{\theta}$, is given by its likelihood function, written $p(\mathbf{x} \mid \boldsymbol{\theta})$, which in some sense encapsulates the statistical properties of the \gls{abm}. We denote simulation of a dataset $\mathbf{x}$ from the \gls{abm} as a draw from its probability density function, $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$.
An essential prerequisite for the successful application of \glspl{abm} in real-world decision making scenarios is the ability to perform parameter estimation, that is, to tune $\boldsymbol{\theta}$ such that the simulated data are good matches to observed data, which we denote with $\mathbf{y}$. Traditional statistical workhorses like maximum likelihood estimation or Bayesian inference rely on pointwise evaluations of the likelihood function. In practice, $p(\mathbf{x} \mid \boldsymbol{\theta})$ cannot be obtained or evaluated in a reasonable time for arbitrary simulation models, and is thus only implicitly defined by the software specifying the behaviour of the simulator. This presents a barrier to the use of \glspl{abm} in real-world use.
To overcome this, a number of approaches to performing statistical inference have been developed in which exact density evaluations are replaced by evaluations of approximate densities or cost functions that are constructed using simulations from the model. Some of these \gls{sbi} methods have been explored within the \gls{abm} community. Perhaps the most prevalent of them is \gls{smd}, in which an estimate $\hat{\boldsymbol{\theta}}$ is obtained by minimising some loss function $f(\mathbf{y}, \boldsymbol{\theta})$ between the observed data $\mathbf{y}$ and simulated data $\mathbf{x} \sim p\left(\cdot \mid \boldsymbol{\theta}\right)$ over some search space $\boldsymbol{\Theta}$:
This general class of estimators includes the maximum likelihood estimator -- corresponding to choosing, for example, $f\left(\mathbf{y}, \boldsymbol{\theta}\right) = -\log{p\left(\mathbf{y} \mid \boldsymbol{\theta}\right)}$ or, where the likelihood is intractable, some estimate thereof diggle1984monte, kukacka2017estimation. It also includes the \gls{msm} franke2009applying -- in which $f$ takes the form
where $g(\mathbf{y})$ denotes a set of moments derived from $\mathbf{y}$, $\hat{g}_{\boldsymbol{\theta}}$ denotes the same set of moments derived from $R \geq 1$ simulations at $\boldsymbol{\theta}$, $W$ is a suitably chosen weight matrix, and $'$ denotes the transpose. As a final example, it includes \gls{ii} gourieroux1993indirect, which follows a similar approach to \gls{msm} but replaces the moments with estimated parameters of a tractable auxiliary model. Some further loss functions have been proposed more recently in the context of \gls{abm} estimation PLATT2020103859.
A major drawback of \gls{smd} and optimisation-based approaches more generally is that only parameter point estimates are produced, while in many contexts it is desirable to obtain meaningful uncertainty quantification regarding appropriate values for $\boldsymbol{\theta}$. Bayesian inference is an alternative inferential paradigm which naturally provides this meaningful notion of uncertainty, and has seen some recent but otherwise limited attention within the \gls{abm} community grazzini2017bayesian, lux2018estimation, PLATT2020103859, platt2021bayesian, lux2021bayesian. Here, the primary object of interest is the parameter posterior distribution $p( \boldsymbol{\theta} \mid \mathbf{y} )$, obtained via Bayes theorem:
In Equation (ref), $p\left(\boldsymbol{\theta}\right)$ is referred to as the prior distribution over parameters $\boldsymbol{\theta} \in \boldsymbol{\Theta}$, which encodes the experimenter's initial beliefs regarding appropriate values for $\boldsymbol{\theta}$ before data is observed, while the simulation model itself appears via its associated data likelihood function $p\left(\mathbf{y} \mid \boldsymbol{\theta}\right)$. Importantly, Bayesian inference procedures seek to accurately approximate the entire distribution $p\left(\boldsymbol{\theta} \mid \mathbf{y}\right)$, rather than any single point estimate derived from this distribution. In this way, accurate and meaningful uncertainty estimates may also be obtained which correctly account for the experimenter's initial beliefs, as well as the evidence provided by the data $\mathbf{y}$.
A well-known investigation into Bayesian estimation methods for \glspl{abm} is grazzini2017bayesian, which explores various means of constructing a surrogate likelihood $\tilde{p}\left(\mathbf{y} \mid \boldsymbol{\theta}\right)$ to account for the intractability of the true likelihood $p\left(\mathbf{y} \mid \boldsymbol{\theta}\right)$. Two immediate drawbacks of the methods grazzini2017bayesian discuss are that they (a) are difficult to reconcile with the fact that \glspl{abm} are frequently dynamical models generating time-series output, and (b) entail a huge computational burden, since each likelihood evaluation in the employed posterior sampling algorithm requires at least one simulation from the model, and it is typically necessary to make many hundreds of thousands of such evaluations. Recent work by platt2021bayesian suggests estimating the model's transition density via a mixture density network. While the approach offers some improvements through the use of a more flexible density estimator and the incorporation of some temporal dependencies, it suffers from similar drawbacks, assuming time-homogeneity and requiring a computationally expensive estimation step for every single sample from the approximate posterior distribution. In summary, there is a need for parameter inference methods that fulfil two main criteria:
To this end, we seek with this paper to analyse the utility of two classes of parameter inference methods that have seen significant activity within the computational statistics and probabilistic machine learning literature in recent years: neural posterior estimation papamakarios2016fast, lueckmann2017flexible, Greenberg2019, and neural density ratio estimation thomas2021lfire, pmlr-v119-hermans20a, Durkan2020. Such methods have been employed successfully in a variety of applied domains, including high-energy physics brhemer2018constraining, cosmology alsing2019fast, and neuroscience gonccalves2020training. In addition, as we will demonstrate below, they work flexibly with potentially multivariate time-series data, requiring minimal model assumptions, and typically generate more accurate parameter inferences with a significantly reduced simulation budget in comparison to common alternatives, such as those reported in grazzini2017bayesian.
We further seek to establish a higher standard of assessment criteria to be used in future benchmarking experiments for approximate Bayesian parameter inference for economic simulation models through the use of common integral probability metrics and inference validation procedures used widely in the statistics and machine learning communities. In Bayesian inference, the object of interest is the posterior density, whose entire shape is important since this captures degrees of beliefs regarding the unknown $\boldsymbol{\theta}$. However, the literature on Bayesian estimation of \glspl{abm} has largely forgone attempts to examine the degree of overlap between the estimated and ground-truth posterior densities or the extent to which the inference procedure has learnt the dynamics of the data, often in favour of a simple comparison of certain point estimates of $\boldsymbol{\theta}$ grazzini2017bayesian, platt2021bayesian, shiono2021estimation. Since the objective of Bayesian inference is to derive a posterior probability density function which naturally quantifies uncertainty, it is important to adopt assessment criteria that quantify the degree of closeness between the estimated and target posterior densities in benchmarking experiments.
In summary, our contributions with this article are to provide:
Bayesian inference has seen relatively little attention from the economic \gls{abm} community as a means to parameter estimation. In this section, we summarise the key literature in this space. Importantly, we emphasise throughout that a common feature of these approaches is that the task of inference -- that is, of constructing the posterior distribution -- is inherently linked to the act of simulating from the \gls{abm} and building an additional approximate generative model (sometimes also called an emulator), in the sense that an approximation to the likelihood function is constructed. We will see in later chapters that this contrasts with a new generation of estimation methods in which inference is decoupled from the act of simulating and can be performed in a discriminative manner, typically resulting in more efficient inferences.
Bayesian estimation of economic \glspl{abm} has gained popularity with the work of grazzini2017bayesian, in which the authors discuss three approaches to approximate Bayesian parameter inference. In this section, we briefly outline the three approaches, and show that each method can seen as explicitly or implicitly approximating the unknown likelihood function associated with the \gls{abm} by
where $K_{\epsilon}$ is a (possibly unnormalised) kernel function, with parameter(s) $\epsilon$, providing a measure of \enquote{distance} between observed data $\mathbf{y}$ and simulated data $\mathbf{x}$. This permits an approximation $\hat{p}\left(\boldsymbol{\theta} \mid \mathbf{y}\right)$ of the posterior density as
Such approaches -- targeting (ref) -- are comprehensively subsumed under \enquote{\gls{abc}} which, as noted by grazzini2017bayesian, is a technique that has received significant attention within the computational statistics community over the last two decades as a means to constructing approximate posterior densities for models with intractable likelihood functions tavare1997inferring, pritchard1999population, beaumont2002approximate, sunnaaker2013approximate, gbiabc. The approach has been successfully applied within a variety of fields, most prominently in ecology, epidemiology, and systems biology toni2009approximate, liepe2014framework, beaumont2010approximate, ju2021sequential. While the approaches discussed by grazzini2017bayesian can be cast in this general form, a judicious choice of a particular kernel $K_\epsilon$ based on reasonable assumptions is crucial. In particular, different kernel choices are called for in different scenarios. Some possibilities are outlined below.
The simplest approach discussed by grazzini2017bayesian involves assuming that the \gls{abm} has entered a statistical equilibrium, such that the observations $y_t \in \mathbb{R}^d$ in the time-series $\mathbf{y} := (y_1, y_2, \dots, y_T) \in \mathbb{R}^{ d\times T }$ are fluctuations around some stationary value $m^{*}$:
where $(\varepsilon_t)$ are iid noise terms with density $g_{\epsilon}$ and parameters $\epsilon$. Although more complex distributions are available, $g_{\epsilon}$ may for example be a zero-mean Gaussian distribution, and $\epsilon$ the elements of the covariance matrix.
Under these assumptions, recalling that the elements of the time-series are assumed to be independent, the parameters $m^{*}$ and $\epsilon$ can be estimated by the means and covariances of the elements of the time-series to give $\hat{m}^*$ and $\hat{\epsilon}$. More precisely, first fixing $\boldsymbol{\theta}$ and drawing a sample $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$, one can calculate the estimates $\hat{m}^*(\mathbf{x})$ and $\hat{\epsilon}(\mathbf{x})$ using, for example, maximum likelihood estimation. Adopting as the measure of distance the probability of the true data being observed under this approximated process yields the choice of kernel
Finally by taking $R \geq 1$ iid simulated datasets\footnote{While $R$ can be chosen to be any natural number, it is usually efficient to take $R=1$ when averages of estimators are concerned and run the Markov chain for an accordingly higher number of iterations bornn2017use, sherlock2017pseudo.} $\mathbf{x}^{(r)} \sim p\left(\mathbf{x} \mid \boldsymbol{\theta}\right)$, $r = 1, \dots, R$, and calculating their associated estimators $\hat{m}^{*}\left(\mathbf{x}^{(r)}\right)$ and $\hat{\epsilon}\left(\mathbf{x}^{(r)}\right)$, the likelihood at $\boldsymbol{\theta}$, i.e. Equation (ref), is approximated with the Monte Carlo average
where we write $\hat{p}_{\boldsymbol{\theta}}$ instead of $\hat{p}(\cdot \mid \boldsymbol{\theta})$. Alternatively, the $R$ simulations may be pooled to generate single estimates $\hat{m}^{*}\left(\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(R)}\right)$ and $\hat{\epsilon}\left(\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(R)}\right)$ at $\boldsymbol{\theta}$. This is closely related to synthetic likelihood approaches wood2010statistical, price2018bayesian which use similar Gaussian distributions, but usually rely on summary statistics.
In either case, the resulting approximate likelihood $\hat{p}_{\boldsymbol{\theta}}$ from Equation (ref) can then be used downstream for parameter estimation, either directly via e.g. maximum likelihood or through the corresponding approximate posterior, simulated e.g. through \gls{mcmc}. This approach suffers from clear limitations, however. Firstly, it treats the data points in the observed time-series as independent -- that is, as lacking a natural ordering -- which destroys important information when the observed $\mathbf{y}$ and simulated $\mathbf{x}$ are time-series. Secondly, choosing an appropriate and sufficiently flexible family of parametric densities $g_{\epsilon}$ to construct the likelihood approximation is non-trivial, with poor choices leading to erroneous Bayesian inference.
An alternative approach, which partially addresses the second limitation described in Section (ref), is to forgo the assumption of a parametric family of densities and instead use a non-parametric method for density estimation. grazzini2017bayesian describe the use of \gls{kde} for this purpose. Here, the data points in the time-series are once again assumed to be independent and fluctuating about some stationary value $m^*$ as in the parametric approach described above. Then, an estimate of the likelihood function is obtained by applying \gls{kde} to $R \geq 1$ iid simulations of length $S$, $\mathbf{x}^{(r)} := (x^{(r)}_1, \dots, x^{(r)}_S) \sim p\left(\mathbf{x} \mid \boldsymbol{\theta}\right)$, $r = 1, \dots, R$, providing an unbiased estimate of the approximated likelihood function in Equation (ref) as
where $\hat{p}_{\epsilon}\left(y_t \mid \boldsymbol{\theta}, \mathbf{x}^{(r)}\right)$ is the estimate of the conditional density $p\left(y_t \mid \boldsymbol{\theta}, \mathbf{x}^{(r)}\right)$ obtained via \gls{kde}:
for some probability kernel $\kappa_{\epsilon}$ with bandwidth parameter(s) $\epsilon$. In particular, the authors employ a Gaussian kernel with bandwidth chosen using Silverman's method silverman1986density, such that
Alternatively, as in Section (ref), the $R$ simulations may be pooled and a single \gls{kde} model fit to the combined dataset to obtain $\hat{p}(\mathbf{y} \mid \boldsymbol{\theta})$. While these non-parametric approaches to density estimation are arguably more flexible than the assumption of a parametric family, they are well known to suffer from the curse of dimensionality, limiting their applicability to low dimensional data only, even under the unrealistic assumption that the $y_t$ are assumed independent.
The authors discuss a third common approach to simulation-based Bayesian inference for intractable simulation models, namely classical \gls{abc}. A typical \gls{abc} algorithm proceeds by proposing a parameter value $\boldsymbol{\theta}$ from some proposal distribution, for example the prior density, and determining whether to accept or reject $\boldsymbol{\theta}$ on the basis of some meaningful notion of distance $d(\mathbf{x}, \mathbf{y})$ between a simulation $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$ and the observed data $\mathbf{y}$. A common choice for the distance $d$ is the Euclidean distance between suitable summary statistics of the data; that is, $d(\mathbf{x}, \mathbf{y}) := \lVert \mathbf{s}(\mathbf{x}) - \mathbf{s}(\mathbf{y}) \rVert$ where $\mathbf{s}$ is a vector of, for example, the mean, standard deviation, and lag-1 autocorrelation of the data. While many versions of \gls{abc} exist, the simplest -- Rejection \gls{abc} -- involves repeatedly proposing candidate parameter values $\boldsymbol{\theta} \sim p(\boldsymbol{\theta})$ and rejecting in the event that $d\left(\mathbf{x}, \mathbf{y}\right) > \epsilon$, where $\mathbf{x} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$. Rejection \gls{abc} can therefore be cast in the form of Equation (ref) by taking
providing an estimate of the (unnormalised) approximate likelihood function in Equation (ref) as
As discussed by grazzini2017bayesian, \gls{abc} enjoys a number of desirable properties, in particular (under mild regularity conditions) consistency with the true posterior distribution as $\epsilon \to 0$, provided no information loss is incurred by the choice of distance $d$. However, a major challenge to its use in practice is designing such a distance function: the success or failure of the approach often hinges on the engineering of summary statistics that are close to sufficient\footnote{A Bayesian sufficient statistic $\mathbf{s}(\mathbf{x})$ for a model $p(\mathbf{x}\mid\boldsymbol{\theta})$ is one for which $p(\boldsymbol{\theta} \mid \mathbf{s}(\mathbf{x})) = p(\boldsymbol{\theta} \mid \mathbf{x})$. Low-dimensional sufficient statistics for arbitrary probabilistic models are generally unobtainable pitman1936sufficient, koopman1936sufficient, darmois1935sufficient.}. A number of works have investigated strategies for summary statistic selection in \gls{abc} Fearnhead2012, blum2013comparative, wiqvist2019partially, Chen2020 and there has been some success in alternative, summary-free distances Park2016, Bernton2019, dyer2021approximate. It is however in general not possible to derive a low-dimensional sufficient statistic for arbitrary simulation models, and a poor choice of summary statistics can lead to an unacceptable loss of information from the original data, resulting in a posterior approximation of little value.
In Algorithm (ref), we show how the approximate posterior (ref), with any of the three possibilities for $K_{\epsilon}$ discussed by grazzini2017bayesian, may be targeted using a variant of the popular \gls{mh} algorithm Hastings1970MonteCS. If the likelihood estimate $\hat{p}_{\boldsymbol{\theta}}(\mathbf{y})$ is an unbiased, non-negative estimate of some desired target $p_{\boldsymbol{\theta}}(\mathbf{y})$, then Algorithm (ref) is an example of a pseudo-marginal \gls{mh} algorithm andrieu2009pseudo. Note that while the sampling procedure in Algorithm (ref) exactly targets the posterior distributions described in grazzini2017bayesian, many alternative posterior sampling algorithms exist beaumont2009adaptive.
We once again emphasise that a major disadvantage to applying the methods described by grazzini2017bayesian to \glspl{abm} is that estimating the posterior involves simulating $R \geq 1$ times from the \gls{abm} at each proposed parameter value $(\boldsymbol{\theta}_i)_{i=1,\ldots,n}$. Since \glspl{abm} are typically expensive to simulate, and $n$ is required to be many hundreds of thousands in order to accurately estimate the targeted posterior, such approaches can rapidly lead to a prohibitively large number of required simulations.
A class of models closely related to approximate Bayesian computation is the class of so-called latent variable models. Here, the real data, $\mathbf{y}$, is assumed to be a noisy observation of an unobserved (latent) process $\mathbf{x}$. Thus, in contrast to previous approaches, some discrepancy between $\mathbf{y}$ and $\mathbf{x}$ is to be expected.
Under the most basic scenario, the data obtained from the simulator are identically $x_t \sim p(x \mid \boldsymbol{\theta})$ and $p(\mathbf{x} \mid \boldsymbol{\theta}) = \prod_{t=1}^T p(x_t \mid \boldsymbol{\theta})$. The observed data is then assumed to have an error distribution, $g_\varepsilon$, say, such that the contribution made by one observation $y_t$ to the likelihood function is thus
and $p(\mathbf{y} \mid \boldsymbol{\theta}) = \prod_{t=1}^T p(y_t \mid \boldsymbol{\theta})$. We draw the attention of the read to the conceptual similarity to (ref); however, (ref) is no longer an approximation but is the exact likelihood under the assumption of measurement error. (Similarly, \gls{abc} can be seen as exact if $K_\varepsilon$ describes the observational error, see e.g. wilkinson2013approximate.) A simple unbiased estimator of (ref) is
for some $R\geq 1$\footnote{Note that opposed to earlier cases the choice $R=1$ is generally not optimal for products of unbiased estimators, see doucet2015efficient, sherlock2015efficiency, schmon2021large.}, which can then also be used in Algorithm (ref), for example, to target the associated approximate posterior. The independence assumption underlying the simulator data is, however, not realistic for \glspl{abm} and other time-series models, and incorporating the temporal aspect of data requires further structural assumptions, such as those of hidden Markov models.
Particle filters, an instance of a broader class of \gls{smc} methods ju2021sequential relax the independence assumption imposed on the output of the simulation model and have previously been employed for Bayesian parameter inference for economic \glspl{abm} lux2018estimation, lux2021bayesian. However, such methods also involve making assumptions about the structure of the model -- specifically, the assumption of a state space structure with a particular error distribution. Furthermore, previous works have noted the significant computational burden required for \gls{smc} methods malleson2020simulating, lux2021bayesian, which can arise due to the fact that a number of particles that grows exponentially with the model dimensions is often required to avoid failure modes such as particle collapse.
As an alternative to kernel density estimation, neural density estimators have previously been employed to perform Bayesian estimation of \glspl{abm}. platt2021bayesian presents a method which diverges from the approaches of grazzini2017bayesian in two main regards. Firstly, in contrast to the previous independence assumption, it assumes an autoregressive structure for the time-series model to better capture temporal correlations in the data. In particular, the approach assumes a time-series model that is Markov of order $L$, that is, $p(x_t \mid x_{1:t-1}, \boldsymbol{\theta}) = p(x_t \mid x_{t-L:t-1}, \boldsymbol{\theta})$, where $x_{1:t} = (x_1, \dots, x_t)$. In words, the distribution of the state at time $t$ depends only on the previous $L$ states (and $\boldsymbol{\theta}$). In light of this assumption, the full likelihood can be factorised as
Secondly, as a consequence of the autoregressive model structure, the method employs a conditional density estimator $q_\phi$ to approximate the transition density function $p(x_t \mid x_{t-L:t-1}, \boldsymbol{\theta})$. The resulting likelihood approximation can be written
where the contribution of the first $L$ terms is ignored. A mixture density network bishop1994mixture assumes the role of the conditional density estimator $q_{\phi(\boldsymbol{\theta})}$ with parameters $\phi(\boldsymbol{\theta})$ trained for each $\boldsymbol{\theta}$. Denoting the $L$ states preceding $y_t$ as $\mathbf{y}_{<t}$ for brevity, the particular form of the mixture density network is
where $\alpha_k$ are the mixing coefficients, and $\mu_k, \Sigma_k$ are the mean and covariance matrix respectively, all of which are parameterised by neural networks and trained via maximum likelihood on $R$ iid model samples $\mathbf{x}^{(r)}\sim p(\mathbf{x} \mid \boldsymbol{\theta})$, $r=1,\ldots, R$. The resulting likelihood estimate can then once again be used also with Algorithm (ref), for example, to target the associated approximate posterior distribution.
While the model structure described above accounts for the sequential nature of simulations generated by \glspl{abm} to some extent, the computational burden associated with the simulation of sufficient training data and the training of a new conditional density estimator at each $\boldsymbol{\theta}$ renders it largely infeasible shiono2021estimation. In other work, shiono2021estimation examines a particular neural approach, with an invertible architecture and the ability to accommodate time-series data, for the specific purpose of estimating the posterior density for a New-Keyensian \gls{abm}. In contrast to both platt2021bayesian and shiono2021estimation, we take in this work a broad perspective on the subject of neural simulation-based inference for \glspl{abm}, by thoroughly contextualising the problem of Bayesian parameter estimation for economic \glspl{abm} and by incorporating both benchmarking tasks and an exploration of suitable assessment criteria.
\glsresetall
In this section, we provide an overview of the two main \gls{sbi} methods we make use of in this paper: \gls{npe} and \gls{nre}, both of which -- as their names suggest -- employ neural networks\footnote{Note that neural networks are not essential here -- the procedures we describe require only flexible function approximators. We frame the discussion around neural networks primarily because they are a convenient choice of flexible function approximators in many settings.} to obtain an estimate of the posterior density $p(\boldsymbol{\theta} \mid \mathbf{y})$. We motivate the use of these methods for \glspl{abm} by framing them as discriminative approaches to simulation-based parameter inference, and contrasting them with the generative approaches described in Section (ref) which seek to approximate the model's distribution (likelihood) using some form of probabilistic model. We then present the core elements of these methods, and further provide a discussion on how they may be leveraged flexibly and automatically to accommodate inference tasks involving high-dimensional and potentially multivariate time-series data, as is often the case in economic applications.
The methods we endorse in this article -- \gls{npe} and \gls{nre} -- can be motivated by the fact that they are:
\paragraph{Simulation efficiency} The algorithms described in Section (ref) share a common pattern. For a fixed parameter $\boldsymbol{\theta}$, iid simulations $\mathbf{x}^{(r)} \sim p(\mathbf{x} \mid \boldsymbol{\theta})$, $r = 1,\ldots,R$, are sampled to produce a proxy likelihood $\hat{p}(\mathbf{x} \mid \boldsymbol{\theta})$. Then, to generate $n$ approximate posterior samples via \gls{mcmc}, this procedure is performed at $n$ different values for $\boldsymbol{\theta}$, where $n$ must typically be a large number -- often a few hundred thousand, if not orders or magnitude larger -- to ensure a low Monte Carlo error. Consequently, at least $nR$ simulations from the \gls{abm} are required in total for inference under these algorithms. Since \glspl{abm} can be very expensive to simulate, and the act of simulating from the \gls{abm} remains the primary bottleneck in Bayesian estimation procedures for \glspl{abm}, this simulation demand can quickly become infeasible.
In contrast, the methods for which we advocate here differ by eliminating the need to simulate when sampling from the posterior. Instead, they decouple the act of simulating from the task of constructing the posterior by employing powerful function approximators to learn -- in essence -- global posterior density estimators on the basis of a limited number of simulations from across the parameter space of the \gls{abm}; that is, they learn functions $h : \mathcal{Y} \times \boldsymbol{\Theta} \to \mathbb{R}$ which approximate the posterior density $p(\boldsymbol{\theta} \mid \mathbf{y})$ across the space of all possible values for $\boldsymbol{\theta}$ and $\mathbf{y}$. This has the potential to significantly reduce the simulation burden associated with approximate Bayesian inference procedures, because the pointwise estimates of $p(\boldsymbol{\theta} \mid \mathbf{y})$ can borrow strength from, and share information between, one another; in contrast, the algorithms described in Section (ref) consider each pointwise evaluation of $p(\boldsymbol{\theta} \mid \mathbf{y})$ as standalone density estimation tasks, such that no information can be shared between them. In this way, a large number of approximate posterior samples can be generated from an algorithm trained on what can in practice be a far smaller number of model samples than is required for the algorithms described in Section (ref).
Finally, the approaches we endorse here are equipped with a further efficiency benefit: amortisation. This phrase captures the fact that the global density estimators can be used to generate samples from an approximate posterior $\hat{p}(\boldsymbol{\theta} \mid \mathbf{y})$ for any data $\mathbf{y}$ without the need to further simulate from the \gls{abm}, which results from the fact that they are trained on the full space of possible values for $\boldsymbol{\theta}$ and $\mathbf{y}$. In contrast, applying the methods described in Section (ref) to a new dataset $\mathbf{y}'$ would entail another $nR$ simulations from the \gls{abm}, multiplying the already high simulation costs.
\paragraph{Black-box inference methods} We saw in Sections (ref) and (ref) that many inference approaches come with restrictive assumptions, for instance, stationarity of the simulated and observed time-series. In general, it is useful to dispense with these assumptions, since they can be difficult to verify and limit the applicability of these methods for arbitrary simulators. Instead, it can be preferable to employ black-box methods that make minimal assumptions about the model and are therefore generically applicable to arbitrary simulators. Doing so enables the modeller to concentrate resources on model design and implementation, rather than on developing bespoke inference algorithms for each new simulator. Furthermore, assumptions such as stationarity are known to be particularly poorly suited to certain economic simulation models such as \glspl{abm}, since these models are known and are even designed to produce non-equilibrium dynamics. Employing inference procedures that are able to handle such dynamics is therefore essential to the task of estimating generic \glspl{abm}.
\paragraph{Discriminative approaches to parameter inference} Discriminative tasks in machine learning are typically simpler than generative tasks, since generative tasks address larger problems than pure discrimination. This is intuitive: for example, it is typically easier for humans and computers alike to distinguish between images of cats and dogs than to generate them. Analogously, it is generally a simpler task to discriminate between complex time-series data than it is to generate such time-series. The approaches described in Section (ref) adopt a generative approach: they each -- either explicitly or implicitly -- seek to derive an approximation to the simulator's likelihood function using a probabilistic model, and thus seek to model the (probability density function of the) simulation output itself. Formally, this may be understood as learning the (stochastic) map $\boldsymbol{\theta} \mapsto \mathbf{x}$. \Gls{npe} and \gls{nre}, in contrast, do not seek to model the simulation output: instead, they map from instances of the simulation output to certain target values which we will describe in more detail below. This can be formally thought of as learning the map $\mathbf{x} \mapsto \boldsymbol{\theta}$. Such an approach to parameter estimation therefore embodies a fundamental departure from the approaches described in Section (ref). It has the potential to be particularly beneficial for \glspl{abm}, which are known to be able to produce complex, non-equilibrium dynamics that are especially difficult to model and generate.
In this section, we provide an introduction to neural conditional density estimators and their use in posterior estimation tasks. The core idea underlying this class of simulation-based Bayesian inference techniques is to use such conditional density estimators to model the parameter posterior density directly. While various conditional density estimators can be used for this purpose, e.g. mixture density networks bishop1994mixture, papamakarios2016fast, we focus here on the case of normalising flows tabak2010density, tabak2013family, rezende due to their widespread use in various density estimation tasks, including in the areas of image generation kingma2018glow and physics noe2019boltz.
Normalising flows involve transforming a simple base distribution into a more complicated one, often with the use of neural networks. Consider a random variable $\mathbf{U} \sim p_{\mathbf{U}}$, where $p_{\mathbf{U}}$ is a probability distribution chosen to be a “simple” base distribution, in the sense that it is easy to both generate samples from $p_{\mathbf{U}}$ and to evaluate $p_{\mathbf{U}}(\mathbf{u})$ for any $\mathbf{u}$. Now consider the transformed random variable $\mathbf{X} = g(\mathbf{U})$, where $g$ is a differentiable, invertible function with differentiable inverse $f := g^{-1}$. Then, by the change of variables formula, we have
where $J_f$ is the Jacobian of $f$. Sampling from $p_{\mathbf{X}}$ then simply involves sampling a value $\mathbf{u}$ from the base distribution $p_{\mathbf{U}}$ and immediately obtaining $\mathbf{x} = g(\mathbf{u})$. Evaluating $p_{\mathbf{X}}(\mathbf{x})$ also then simply involves evaluation of the right-hand side of Equation (ref). The above may be extended easily to a composition, or flow, $g = g_{n} \circ g_{n-1} \circ \dots \circ g_{1}$ of differentiable and invertible functions $g_{i}$ with inverses $f_i$, for which we now have that
Given samples from $p_{\mathbf{X}}$, the question of how to choose a base distribution $p_{\mathbf{U}}$ and a series of functions $g_1,\ldots,g_n$ such that the distribution of the samples is modelled accurately arises. The idea behind normalising flows is to learn this sequence of transformations, that is, to have each $f_{i}$ be a flexible transformation $f_{\phi_i}$ with trainable parameters $\phi_i$. The resulting density estimator $q_{\phi}$ can be written
The simple base distribution is then often taken to be a standard normal distribution\footnote{For this reason, the composition $f = f_{1} \circ \dots \circ f_{n-1} \circ f_{n}$ is often referred to as the normalising flow, since it is a flow of transformations resulting (typically) in a Gaussian base distribution.} of the same dimension as the target distribution, and the flexible transformations $f_{\phi_i}$ are typically neural networks with a structure designed to guarantee the required invertibility and differentiability. The trainable parameters $\phi = \left\{ \phi_i \colon 1 \leq i \leq n \right\}$ are optimised by maximising the log-likelihood of the data $\mathbf{x}^{(r)} \overset{iid}{\sim} p_{\mathbf{X}}, r = 1, \ldots, R$:
Normalizing flows are typically implemented in machine learning libraries that support automatic differentiation such as pytorch paszke2019pytorch or TensorFlow abadi2016tensorflow allowing the effective use of gradient-based optimisation techniques, such as Adam kingma2014adam. We provide a schematic of the sampling (flow, right-pointing arrows) and density evaluation (normalising flow, left-pointing arrows) processes in Figure (ref), in which the simple, left-most distribution is morphed over successive steps into the more complex, right-most distribution.
Normalising flows may also be extended to the case of conditional density estimation papamakarios2016fast, lueckmann2017flexible, papamakarios2019sequential, Greenberg2019. In this way, we can use a normalising flow to estimate the posterior density $p(\boldsymbol{\theta} \mid \mathbf{x})$ associated with a simulation model by following the same procedure as above and finding the neural network parameters as
where $(\mathbf{x}^{(r)}, \boldsymbol{\theta}^{(r)}) \overset{iid}{\sim} p(\mathbf{x}, \boldsymbol{\theta}), r=1, \ldots, R.$ Importantly, this framework allows us to learn a single conditional density estimator for the posterior density function across data $\mathbf{x}$, rather than having to undergo the expensive procedure of training a separate density estimator for each point evaluation of the posterior density as in e.g. platt2021bayesian.
\glsreset{nre}
In this section, we describe an alternative common approach to \gls{sbi} termed \gls{nre}. The core idea underlying this method contrasts with \gls{npe} in that while \gls{npe} seeks to model the posterior density directly, \gls{nre} instead seeks to estimate the following ratio:
This ratio is frequently referred to as the likelihood-to-evidence ratio. In the event that Bayesian inference is the goal, this then permits evaluation of the posterior density as
In the following, we discuss how this ratio may be efficiently estimated with well-calibrated probabilistic classifiers trained on simulated data-parameter pairs.
We outline the approach to approximate Bayesian inference via density ratio estimation as described by pmlr-v119-hermans20a. Here, a probabilistic binary classifier is trained to distinguish between two sets of training examples:
The difference between the two cases is that the $\mathbf{x}$ are generated by the $\boldsymbol{\theta}$ to which they are paired in the former set of examples, while in the latter set of examples the $\mathbf{x}$ bear no relation to the $\boldsymbol{\theta}$ they are paired with.
The function of a probabilistic binary classifier trained on such data is to model the class probability $d(\mathbf{x}, \boldsymbol{\theta}) := p(c = 1 \mid \mathbf{x}, \boldsymbol{\theta}) \in \left[0, 1\right]$; hard classification labels (i.e. decisions regarding the predicted value of $c$) are obtained when the continuous-valued $d(\mathbf{x}, \boldsymbol{\theta})$ is combined with a decision rule, for example the rule that $c = 1$ should be predicted whenever $d(\mathbf{x}, \boldsymbol{\theta}) > 0.5$. One can show pmlr-v119-hermans20a that the optimal estimate of $d(\mathbf{x}, \boldsymbol{\theta})$ is the value
Thus, a good probabilistic classifier trained to distinguish between the two possible types of pairs $(\mathbf{x}, \boldsymbol{\theta})$ will learn a good estimate $\hat{d}(\mathbf{x}, \boldsymbol{\theta})$ of this ratio. Such a probabilistic classifier will allow us to evaluate the posterior density in this way: by noticing that one can rearrange Equation (ref) as
where $r^{*}(\mathbf{x}, \boldsymbol{\theta})$ is the corresponding estimate of the likelihood-to-evidence ratio $p(\mathbf{x} \mid \boldsymbol{\theta})/p(\mathbf{x})$. In practice, of course, only an approximation $\hat{d}(\mathbf{x}, \boldsymbol{\theta})$ will be obtained to the ratio in Equation (ref), yielding a correspondingly imperfect estimate $\hat{r}(\mathbf{x}, \boldsymbol{\theta})$ of the likelihood-to-evidence ratio. Nonetheless, it is known that training expressive, high-capacity neural network classifiers via the cross-entropy loss
can yield classifiers with good probability estimates, and thus good estimates of $r^{*}(\mathbf{x}, \boldsymbol{\theta})$.
In more recent work, Durkan2020 demonstrate that the above approach to density ratio estimation can be generalised to the problem of training a probabilistic classifier to identify the correct $(\mathbf{x}, \boldsymbol{\theta}^{(i)})$ pair from a batch $\lbrace{(\mathbf{x}, \boldsymbol{\theta}^{(b)})}\rbrace_{b=1}^{B}$ of $B$ pairs that otherwise contain only “incorrect” pairs, where “correct” and “incorrect”, respectively, correspond to having been drawn from the joint distribution $p(\mathbf{x}, \boldsymbol{\theta})$ and from the product of the marginal distributions $p(\mathbf{x})p(\boldsymbol{\theta})$. In this case, the optimal estimate of the correct class probability is
Thus, by comparison with Equation (ref), training a neural network $f_{\phi}$ on the loss function
for each $(\mathbf{x}, \boldsymbol{\theta}) \sim p(\mathbf{x} \mid \boldsymbol{\theta})p(\boldsymbol{\theta})$ will induce $f_{\phi}(\mathbf{x}, \boldsymbol{\theta})$ to learn the value $f_{\phi}(\mathbf{x}, \boldsymbol{\theta}) = \log(p(\boldsymbol{\theta} \mid \mathbf{x})/p(\boldsymbol{\theta}))$, thus recovering an estimate of the desired density ratio. This can be extended to a batch of $R$ “correct” data-parameter pairs as
where the terms in the denominator are labelled with $b_r$ to account for the possibility that the contrasting (“incorrect”) set of parameters may be different for different “correct” pairs $(\mathbf{x}^{(r)}, \boldsymbol{\theta}^{(r)})$. The authors further demonstrate that this can yield more accurate density ratio estimators. For this reason, we adopt this approach throughout this article.
After successful training, the posterior density estimator $q_{\hat\phi}(\boldsymbol{\theta} \mid \mathbf{x})$ obtained from \gls{npe} is a parametric approximation of the true posterior distribution, which can then be used to generate iid samples $\boldsymbol{\theta} \sim q_{\hat\phi}(\boldsymbol{\theta} \mid \mathbf{x})$ or to evaluate the posterior density on a pointwise basis in the ways described in Section (ref). While ratio estimation similarly permits a pointwise evaluation of the posterior distribution as $p(\boldsymbol{\theta})\exp{(f_{\hat\phi}(\mathbf{x}, \boldsymbol{\theta}))}$, it requires a further inference step using, for example, Metropolis--Hastings to generate samples from the same posterior. However, the upfront training of the ratio estimator eliminates the need for expensive sampling from the \gls{abm}, reducing the run-time significantly in comparison to alternative approaches such as those described in Section (ref).
The procedures described in Sections (ref) and (ref) are framed as single density (ratio) estimation tasks. This means that a single dataset of $R$ data-parameter pairs $ (\mathbf{x}^{(r)}, \boldsymbol{\theta}^{(r)}), r = 1, \ldots, R$ is created in the following way:
Subsequently, \gls{nre} and \gls{npe} algorithms are trained on the whole dataset to find either the likelihood-to-evidence ratio or the posterior directly following, respectively, (ref) or (ref). Density estimators trained in this way are referred to as amortised density estimators, which reflects the fact that they may be used to construct posterior distributions for any data $\mathbf{x}$ without further simulation from the \gls{abm} or the need to train multiple density estimators.
The disadvantage of this one-stage approach is that all training simulations from the \gls{abm} are generated by parameters drawn from the prior distribution. In many cases, however, interest lies only in the posterior for the particular observation $\mathbf{y}$, which can be concentrated on specific subregions of the entire space covered by the prior density. It can then be preferable to narrow down the search space of good candidate values for the parameter $\boldsymbol{\theta}$ to subregions of the parameter space that could most plausibly have generated $\mathbf{y}$. By doing so, the density (ratio) estimator will be presented with a less varied range of dynamics between which it must learn to distinguish, allowing it to develop a more refined approximation of the density (ratio) in regions of high posterior density. This can facilitate more rapid learning and potentially reduce the number of training examples that must be simulated by the \gls{abm}.
Significant effort has thus been extended towards the constructions of \enquote{round}-based approaches to training density (ratio) estimators for \gls{sbi} papamakarios2019sequential, Greenberg2019, pmlr-v119-hermans20a, Durkan2020. The idea is to split the total budget of $R$ simulations into subsets of size $N_1, N_2, \ldots$ such that $\sum_i N_i = R$. In the first step, $N_1$ data points are created as before and a posterior approximation $q_{\hat\phi}\big(\boldsymbol{\theta} \mid \mathbf{x}\big)$ is constructed. Subsequently, in round $i \geq 2$, we create new data $(\mathbf{x}^{(r)}, \boldsymbol{\theta}^{(r)}), r = 1, \ldots, N_i$ using
The posterior $q_{\phi}\big(\boldsymbol{\theta} \mid \mathbf{x}\big)$ may now be retrained using (a combination of the old and) the new samples, and the process can be repeated for as many rounds as necessary. An identical process can be followed for density ratio estimation, with the exception that parameters in round $i$ are drawn using the likelihood-to-evidence ratio obtained from round $i-1$ and, for example, Metropolis--Hastings.
While the round-based approach is appealing, it does not come without additional challenges. In particular, the joint distribution of the example pairs $(\mathbf{x}^{(r)}, \boldsymbol{\theta}^{(r)})$ is no longer $p(\mathbf{x}, \boldsymbol{\theta})$ for data sampled after the first round. In practice, this needs to be accounted for during training by the use of additional weighing factors, such as those appearing in Equation (ref) (cf. Equation (ref)). We refer the interested reader to the following papers for further details on this matter: papamakarios2016fast, lueckmann2017flexible, Greenberg2019.
Training schemes for \gls{npe} and \gls{nre} which make use of this round-based training approach are referred to as \gls{snpe} and \gls{snre}, respectively, in which the prefix “sequential” reflects the round-based training design of the now non-amortised density (ratio) estimator.
As described above, \gls{npe} and \gls{nre} take as input either $\mathbf{x}$ alone or $(\mathbf{x}, \boldsymbol{\theta})$ pairs. Given that $\mathbf{x}$ is a (possibly multivariate) time-series for many \glspl{abm} and is thus a high-dimensional object, it can be beneficial to incorporate useful inductive biases\footnote{Inductive biases are usually realised as constraints on the form of the network architecture in the case of neural networks. For example, when certain symmetries are known to be present in the data, this inductive bias may be incorporated with the use of (partially) exchangeable neural networks NIPS2017_f22e4747, wiqvist2019partially. Imposing such constraints helps to restrict the space over which appropriate functions are searched for and facilitates effective learning and generalisation from training data.} into the network architecture that account for the sequential nature of $\mathbf{x}$ and reduce its dimensionality. To this end, the density (ratio) estimator can be prefixed with a so-called embedding network with trainable parameters $\varphi$, whose function is to consume the original high-dimensional dataset $\mathbf{x}$ and express this as a lower-dimensional summary statistic vector $\mathbf{s}_{\varphi}(\mathbf{x})$. The parameters of the embedding network and of the density (ratio) estimator may then be learned concurrently using the same loss function, offering a means to automatically learn descriptive low-dimensional features of the raw input data during the same posterior or density ratio estimation task in an end-to-end fashion.
This quality is a further attractive feature of \gls{npe} and \gls{nre}. Finding low-dimensional summary statistics of the data is necessary in many approaches to \gls{sbi}\footnote{This is the prototypical approach to \gls{abc}, for example, although methods that obviate the need to summarise high-dimensional data have been proposed and explored in recent years Park2016, Bernton2019, dyer2021approximate.}, and the question of how to summarise data naturally arises in response to this. Popular solutions that have been explored to-date include: the experimenter themselves devising a collection of hand-crafted summary statistics that they believe adequately describes the data, which is an arduous task that can lead to a difficult-to-quantify loss of information; or learning summary statistics by constructing a separate learning task to the existing problem of density (ratio) estimation Fearnhead2012, Chen2020, imposing an additional burden on the experimenter. \Gls{npe} and \gls{nre}, in contrast, offer a convenient approach to incorporating inductive biases and to learning summary statistics in an end-to-end fashion without the additional complications associated with alternative \gls{sbi} techniques.
To summarise and conclude this section, we follow Greenberg2019 and Durkan2020 and provide in Algorithms (ref) and (ref) procedures for training neural posterior and density ratio estimators, respectively, for \gls{sbi} over multiple rounds. In both cases, we assume that the density (ratio) estimator includes any embedding network used to learn summary statistics concurrently with the density (ratio) estimate, and thus consider $\phi$ to contain the parameters of both the estimator and the embedding network.
In this section, we present experiments in which we compare the ability of \gls{npe} and \gls{nre} to estimate parameter posterior distributions for economic simulation models with tractable likelihood functions against the non-parametric density estimation method described in grazzini2017bayesian, which we term kde\ and outline in Section (ref). We provide details of the neural network architecture and training hyperparameters in Appendix (ref).
Each example presented in this section will be equipped with a tractable transition density, and therefore a tractable likelihood function. Such models are thus amenable to Bayesian inference via standard means, such as \gls{mcmc}, to obtain an approximate ground-truth posterior density\footnote{While the true posterior density is targeted with such a sampling scheme, the ground-truth obtained in this fashion remains an approximation due to the Monte Carlo error associated with the finite sample size.}.
To assess the accuracy of the estimated posteriors in this case, we compare the full ground-truth posterior density with the full posterior estimated with each of the implemented simulation-based approaches. This contrasts with previous studies of Bayesian parameter estimation for \glspl{abm}, in which point estimates alone are often used to assess the quality of the tested inference procedures. For example, platt2021bayesian uses a prior-weighted Euclidean distance between the estimated posterior mean and generating parameters (that is, the $\boldsymbol{\theta}$ that generated the pseudo-observation $\mathbf{y}$) as a performance metric. However, such metrics provide a very limited and at times misleading view on the outcome of the inference process, because (a) \enquote{closesness} is not measured in the geometry of the target distribution, as visualised in (ref); (b) it is possible that an approximate posterior can produce a posterior mean close to the generating parameter but simultaneously under- or over-estimate the width of the distribution or otherwise yield poor uncertainty quantification, and (c) finite datasets are not guaranteed to yield posterior densities that concentrate on the generating parameter.
To compare the ground-truth and simulation-based posteriors, we compute two integral probability metrics which each correspond to a notion of dissimilarity between the ground-truth and estimated posteriors:
\paragraph{Wasserstein distance}
This is a distance measure derived from optimal transport theory wass and used widely in various machine learning contexts, e.g. generative adversarial networks wgan. For some distance $\rho_0$ on $\boldsymbol{\Theta}$ and two probability measures $\mu$ and $\mu'$, the $p$-Wasserstein metric between two sets of samples $\lbrace{\boldsymbol{\theta}_i}\rbrace_{i=1}^{n} \overset{iid}{\sim} \mu$ and $\lbrace{\boldsymbol{\theta}'_i}\rbrace_{i=1}^{m} \overset{iid}{\sim} \mu'$ is computed as
where $\Gamma_{n,m}$ is the set of $n\times m$ matrices with non-negative entries, columns summing to $m^{-1}$, and rows summing to $n^{-1}$. Throughout, we use the Euclidean distance for $\rho_0$.
\paragraph{Maximum mean discrepancy (MMD)} This metric is once again a metric on probability distributions that draws from the theory of reproducing kernel Hilbert spaces gretton2006kernel, mmd and is used widely within the machine learning and simulation-based inference community to assess the dissimilarity between probability distributions papamakarios2019sequential, lueckmann2021benchmarking. Here, under a suitable choice of kernel\footnote{That is: a symmetric, positive semi-definite function on $\boldsymbol{\Theta}$.} $k$ chosen by the experimenter, the discrepancy between two probability distributions $P$ and $Q$ is taken to be
where $\mathcal{H}$ is the \gls{rkhs} associated with $k$ and $\mathbb{E}_{\boldsymbol{\theta}\sim P}[k(\boldsymbol{\theta}, \cdot)]$ is the so-called mean embedding of $P$ via kernel $k$. When only samples $\lbrace{\boldsymbol{\theta}_i}\rbrace_{i=1}^{n} \overset{iid}{\sim} P$ and $\lbrace{\boldsymbol{\theta}'_i}\rbrace_{i=1}^{m} \overset{iid}{\sim} Q$ are available, an unbiased estimator of this metric can be computed as
Throughout, we use a Gaussian kernel as $k$,
where the scale parameter $\sigma^2 = \text{median}\lbrace{\| \tilde{\boldsymbol{\theta}}_i - \tilde{\boldsymbol{\theta}}_j \|^2_2\rbrace}$, following Briol2019, and $\lbrace{\tilde{\boldsymbol{\theta}}_i}\rbrace_{i=1}^{n}$ are the samples drawn from the ground-truth posteriors.
We consider a variant of the model proposed by BROCK19981235, which has previously been used in \gls{abm} calibration experiments PLATT2020103859. The model dynamics can be expressed as the following system of coupled equations:
where $R, \beta, \sigma$ are parameters. We follow PLATT2020103859 and assume that $H=4, R = 1.0, \sigma = 0.04, g_1 = b_1 = b_4 = 0$ and $g_4 = 1.01$. $\beta$ will be $120$ or $10$, depending on the experiment, and we note below which value is used. We consider the task of estimating the posterior $p\left(\boldsymbol{\theta} \mid \mathbf{y}\right)$, where $\boldsymbol{\theta} = \left(g_2, b_2, g_3, b_3\right)$, $\mathbf{y} := (y_1, y_2, \dots, y_T) \sim p(\mathbf{x} \mid \boldsymbol{\theta}^{*})$ is the pseudo-observation, $T=100$, and $\boldsymbol{\theta}^{*}$ is the parameter setting used to generate $\mathbf{y}$. We consider two experimental setups which we describe further in the following two subsections. For each experiment, we note that by rewriting the above system of equations, we are able to find the transition density for observation $y_{t+1}$ as
where
In this way, we are able to obtain approximate ground truth posteriors with standard \gls{mcmc} techniques such as \gls{mh}.
For both \gls{npe} and \gls{nre}, we use a round-based training approach i.e. \gls{snpe} and \gls{snre}: we train over 10 rounds and generate 1000 simulations in each round. For kde, we take $R = 1$ and sample from the posterior with \gls{mh} (see Appendix (ref) for further details).
\paragraph{Summarising data} As discussed in Section (ref), there exists a number of approaches to representing high-dimensional time-series data $\mathbf{x}$ as a low-dimensional vector $\mathbf{s}(\mathbf{x})$ of summary statistics to facilitate \gls{sbi}. In the experiments for the Brock & Hommes model, we demonstrate \gls{npe} and \gls{nre} using summary statistics obtained in two different ways:
We take $\beta = 120$ and consider the task of estimating the posterior density $p(\boldsymbol{\theta} \mid \mathbf{y})$, where $\boldsymbol{\theta}^{*} := (g_2^{*}, b_2^{*}, g_3^{*}, b_3^{*}) = (0.9, 0.2, 0.9, -0.2)$. We use the following uniform priors: $g_2, b_2, g_3 \sim \mathcal{U}\left(0, 1\right)$, while $b_3 \sim \mathcal{U}(-1,0)$.
In Figure (ref), we show the posteriors obtained using different posterior estimation methods: Figure (ref) shows the approximate ground-truth posterior density obtained with the true likelihood function and \gls{mh}; Figure (ref) shows the posterior obtained with kde and \gls{mh}; Figures (ref) and (ref) show the posterior estimated via \gls{snpe} with naive and learned summary statistics, respectively; and Figures (ref) and (ref) show the posterior obtained via \gls{snre} and naive and learned summary statistics, respectively, which were also sampled using \gls{mh}. Here, we use an embedding network consisting of two stacked Elman recurrent units with hidden state of size 32, followed by a single linear layer of size 16. In each figure, the marginals and joint bivariate densities are located on the diagonal and upper diagonal, respectively, while the red lines/dots locate the mean of the true posterior.
We see from the approximate ground-truth in Figure (ref) that the marginal posteriors are sharply peaked on the generating parameters for $b_2, g_3$, and $b_3$, while the ground truth marginal for $g_2$ is shifted towards higher values. While these features are recovered reasonably well for Figures (ref)--(ref) (the most notable exception being the posteriors for $g_2$), we see that they are recovered poorly with kde. This performance gap is also manifested in the Wasserstein\ and mmd\ metrics reported in Table (ref), where lower values indicate better estimates of the posterior. In summary, kde\ both requires a far larger simulation budget to estimate a single posterior, while simultaneously generating worse estimates of that posterior. In contrast, \gls{snpe} and \gls{snre} is able to achieve superior estimates of the ground-truth posterior distribution, despite the 10-fold reduction in the simulation budget they have been afforded.
We now take $\beta = 10$ and consider the task of estimating the posterior density $p(\boldsymbol{\theta} \mid \mathbf{y})$, where $\mathbf{y} := (y_1, y_2, \dots, y_T) \sim p(\mathbf{x} \mid \boldsymbol{\theta}^{*})$, $T=100$, and $\boldsymbol{\theta}^{*} := (g_2^{*}, b_2^{*}, g_3^{*}, b_3^{*}) = (-0.7, -0.4, 0.5, 0.3)$. In this case, we use the following priors: $g_2, b_2 \sim \mathcal{U}(-1,0)$ and $g_3, b_3 \sim \mathcal{U}(0,1)$.
In Figure (ref), we show the posteriors obtained using different posterior estimation methods: Figure (ref) shows the approximate ground-truth posterior density obtained with the true likelihood function and \gls{mh}; Figure (ref) shows the posterior obtained with kde and \gls{mh}; Figures (ref) and (ref) show the posterior estimated via \gls{snpe} with naive and learned summary statistics, respectively; and Figures (ref) and (ref) show the posterior obtained via \gls{snre} and naive and learned summary statistics, respectively, which were also sampled with \gls{mh}. In this experiment, we now use an embedding network consisting of two stacked \glspl{gru} with hidden state of size 32, followed by a single linear layer of size 16. This results in an embedding network with approximately 10,000 trainable parameters.
We see from the approximate ground-truth in Figure (ref) that the marginal posteriors for $b_2$ and $b_3$ remain sharply peaked on the generating parameters, but that the ground truth marginals for $g_2$ and $g_3$ are diffuse and sloped. The diffuseness of the posterior with respect to these two dimensions is also visually apparent in the joint density plots in the upper diagonal. Upon inspection of the estimated posteriors, we see that the overall shape is recovered well by \gls{snpe} and \gls{snre} with learned summary statistics, but less accurately by kde\ and \gls{snpe} and \gls{snre} when using the naive hand-crafted summary statistics described in Section (ref). This is once again reflected by the Wasserstein\ and mmd\ scores, reported in Table (ref), in which we see significantly decreased values for \gls{snpe} and \gls{snre} with learned summary statistics with respect to the alternatives, despite a 10-fold decrease in the simulation budget they have been afforded. It is nonetheless noteworthy that even with the naive hand-crafted summary statistics, \gls{snpe} and \gls{snre} achieve comparable and slightly favourable performance than kde, again with a 10-fold decrease in the allotted simulation budget.
In Section (ref), we demonstrated that \gls{npe} and \gls{nre} can be used to generate Bayesian parameter posteriors that better match the ground truth posterior than kde\ for the univariate Brock & Hommes model, despite the fact that the latter method entailed 10 times as many simulations as the former two methods based on neural networks. In this section, we seek to demonstrate that this remains the case even for models that generate multivariate time-series as output.
To this end, we consider a \gls{mvgbm}, which is used in a variety of applications in financial time-series modelling. The model is a stochastic differential equation in which each of the $d$ components, labelled $i \in \lbrace{1, 2, \dots, d\rbrace}$, of the path $\mathbf{X}_t = \left(X^1_t, X^2_t, \dots, X^d_t\right) \in \mathbb{R}^d_{+}$ evolve according to
where the $b_i$ are drift coefficients, the $\sigma_{ij}$ are volatility coefficients, and $W^i_t$ is a Brownian motion.
We consider the case of $d=3$ and the task of estimating the posterior for the parameters $\boldsymbol{\theta} = (b_1, b_2, b_3)$ given an observation $\mathbf{y} \sim p\left(\mathbf{x} \mid \boldsymbol{\theta}^{*}\right)$ of $T=100$ points spaced equally with spacing $\Delta t = 1/(T-1)$, where $\boldsymbol{\theta}^{*} = \left(0.2, -0.5, 0.0\right)$. We take priors $b_i \sim \mathcal{U}(-1, 1)$ for each $i = 1, 2, 3$. We note that this model once again permits both exact simulations and samples from the exact posterior since the transition density is once again tractable\footnote{Assuming that $\sigma$ is of full rank.} and can be written as
where
In our experiments, we take
In Figure (ref) we show the approximate ground truth posterior obtained with \gls{mh} and the transition density described above, while in Figures (ref), (ref), and (ref) we show the posteriors obtained with kde and \gls{mh}, \gls{npe} with learned summary statistics, and \gls{nre} with \gls{mh} also with learned summary statistics, respectively. The corresponding simulation budgets are $10^6$, $10^3$, and $10^3$, respectively. To learn the summary statistics, we once again use the embedding network described in Section (ref) and learn summary statistics and the density (ratio) estimator concurrently.
We see that the approximate ground-truth posteriors are relatively diffuse, with peaks approximately coinciding with the true posterior mean, shown with red lines/dots. The shape and degree of diffuseness is captured accurately by \gls{npe} and \gls{nre}. In contrast, the posterior obtained with kde\ is insufficiently diffuse and biased, and thus a significantly worse estimate of the ground-truth posterior. These observations are corroborated by the corresponding Wasserstein\ and mmd\ metrics, reported in Table (ref) along with the corresponding simulation budgets. In summary, \gls{npe} and \gls{nre} achieve significantly more accurate posterior estimates here with a 1000-fold decrease in the simulation budget, and are able to flexibly accommodate multivariate time-series as input for the inference problem.
As previously mentioned, a major benefit of the Bayesian inferential paradigm is the fact that uncertainty quantification is a built-in feature captured via the posterior distribution. Its utility relies, however, on the ability to capture the correct degree of diffuseness in the posterior, such that the uncertainty quantification in the recovered posterior is meaningful. This presents a challenge in approximate Bayesian inference, since by definition the experimenter does not know the correct shape of the posterior when exact inference is intractable. In this section, we address this issue by outlining one widely used approach to verifying the accuracy of approximate Bayesian inference pipelines in the absence of any ground-truth posterior densities.
Simulation-based calibration is a general purpose method for validating approximate Bayesian inference pipelines. The core idea is to use the fact that the data-averaged posterior should be identical to the prior distribution $p(\boldsymbol{\theta})$. That is, a perfect Bayesian inference pipeline should uphold the following equality:
A deviation from this equality signifies some error in the posterior sampling procedure; thus, it has been proposed by talts2020validating that testing how close our Bayesian workflow is to satisfying this equality is a test of the accuracy of the Bayesian pipeline. This may be performed by repeating the following steps $P$ times:
It can be shown talts2020validating that the rank statistics obtained as above should follow a discrete Uniform distribution on $\left\{0, \dots, L\right\}$ for any joint distribution $p(\mathbf{y} \mid \boldsymbol{\theta})p(\boldsymbol{\theta})$ for one-dimensional $\boldsymbol{\theta}$, while for $d$-dimensional $\boldsymbol{\theta}$ each of the $d$ components' rank histograms should be uniform. Inspecting the discrepancy between the true distribution of rank statistics and the desired Uniform distribution thus gives an indication of the accuracy of the Bayesian pipeline; furthermore and conversely, particular patterns of deviation from the desired uniformity are interpretable and provide insight into the specific way in which the obtained posteriors are inaccurate talts2020validating. Performing this procedure, and inspecting the resultant rank histograms, is referred to as performing \gls{sbc}.
While this procedure is in principle possible for kde, the computational burden is immense even for cheap simulation models, and becomes infeasible for large-scale \gls{abm}s. This is because the number of simulations required to generate $n$ samples from the posterior would, in general, be at least\footnote{In practice, this number may be larger due to the necessity of performing trial runs to estimate the parameters of the proposal distribution.} $n R$, such that the total number of simulations required to perform \gls{sbc} for kde\ increases to $P n R$. Even for a moderately sized \gls{sbc} task with $P \simeq 10^3$ and the most optimistic $R=1$, the total number of simulations required can easily reach the order of $10^8$, since it is typically necessary to take $n$ to be many hundreds of thousands to obtain a reasonably accurate estimate of the posterior. In fact, larger values of $n$ are likely to be necessary, since \gls{abm}s are usually highly parameterised and thus can have large parameter spaces.
In contrast, such an approach to verifying the accuracy of posteriors derived from \gls{npe} and \gls{nre} is feasible due to the fact that such approaches involve the prior training of a global density (ratio) estimator, obviating any simulations during the inference (posterior sampling) phase. Thus the strength of \gls{npe} and \gls{nre} derives from not only their enhanced performance and improved simulation-efficiency, but also from the fact that such methods permit accuracy checks such as \gls{sbc} that are otherwise infeasible for alternative parameter estimation techniques when the simulator is expensive, as is often the case for large-scale \glspl{abm}.
To further demonstrate the ability of \gls{npe} and \gls{nre} to recover accurate posterior estimates and informative summary statistics in a highly automated manner, we perform \gls{sbc} using a posterior and density ratio estimator trained on the Franke & Westerhoff model franke2012structural, which has previously been used in benchmarking experiments for \gls{abm} estimation methods platt2021bayesian. In particular, we consider the model version referred to as the “Wealth & Predisposition” model which, similarly to the Brock & Hommes model (see Section (ref)), may be written as a system of coupled equations which models the asset price dynamics that result from a heterogeneous system of traders:
where the $\epsilon^j_t \sim \mathcal{N}(0,1)$ are standard normal random variables and $j = \lbrace{f, c\rbrace}$ labels fundamentalist and chartist traders, respectively. We take as output from the model the time series of log returns $r_t = p_t - p_{t-1}$, and assume the task of training a global, amortised density and density ratio estimator via \gls{npe} and \gls{nre}, respectively, using $10^4$ simulations of length $T=100$ from the simulator defined by Equations (ref)--(ref) with the following parameter settings: $\mu = 0.01$, $\beta = 1$, $\phi = 1$, $\chi = 0.9$, $\alpha_0 = 2.1$, $\sigma_f = 0.752$. The parameter vector for which we seek the posterior and density ratio estimators is taken to be $\boldsymbol{\theta} = \left(\alpha_w, \eta, \sigma_c\right)$, and we obtain iid samples from the posterior associated with the density ratio estimator via a sampling-importance-resampling scheme described in Appendix (ref). (Samples may be drawn iid from the \gls{npe} posterior by construction.)
We show in Figure (ref) the rank histograms obtained via the simulation-based calibration procedure in Section (ref) with \gls{npe} (Figure (ref)) and \gls{nre} (Figure (ref)). To construct these histograms, we take $P = 5,000$ and the following uniform priors: $\alpha_w \sim \mathcal{U}\left(0, 15000\right)$, $\eta \sim \mathcal{U}(0, 1)$, and $\sigma_c \sim \mathcal{U}\left(0, 5\right)$. The black line and gray band denotes the expected height and the expected variation in the heights, respectively, of the bars at this $P$ and with this number of bins. We see that the bars tend to lie within the expected range, although some notable deviations exist: for example, the rank histograms for $\alpha_w$ for both \gls{npe} and \gls{nre} exhibit a minor bias towards larger rank values, which suggests that the marginal posteriors for $\alpha_w$ in both cases is slightly biased towards lower values on average than the true posteriors talts2020validating. However, we emphasise that a major benefit of \gls{npe} and \gls{nre} is that they allow the modeller to make statements of this sort in the first place, whereas the number of simulations required to make such statements in the case of more traditional techniques such as kde\ would quickly become prohibitively large.
\glsresetall
In this paper, we investigated the use of two recently developed approaches to Bayesian parameter estimation for intractable simulation models: \gls{npe} for approximating the posterior density directly, and \gls{nre} for approximating the likelihood-to-evidence ratio. We motivated the use of these parameter inference methods as simulation-efficient black-box discriminative approaches to Bayesian estimation for complex time-series simulators such as agent-based models in economics, which contrast with existing methods that (a) entail an often prohibitively large simulation burden, (b) can make restrictive assumptions about the form of the simulator, for example by ignoring important structural features such as temporal dependencies, and (c) are inherently generative, and thus assume a task that is typically more difficult than discrimination alone. We argue that this latter point is likely to be particularly pertinent for the case of agent-based models in the social sciences, since they are known to be able to generate complex and stochastic non-equilibrium dynamics that can be difficult to model and generate.
We further reviewed existing alternatives, and benchmarked \gls{npe} and \gls{nre} against the most popular of these. In these examples, we see that \gls{npe} and \gls{nre} generally achieve superior recovery of the approximate ground-truth posterior distributions, despite requiring simulation budgets that are orders of magnitude lower than traditional alternatives. In addition, we demonstrated that a verification of the accuracy of the density (ratio) estimators is possible with the introduced methods via \gls{sbc}, due to the amortisation of the estimators. In contrast, such checks would necessitate an unreasonably large number of simulations for existing methods in the literature for Bayesian estimation of economic simulation models, for example the methods discussed by grazzini2017bayesian. For these reasons, we argue that such simulation-efficient black-box Bayesian inference as \gls{npe} and \gls{nre} may enable economists and social scientists alike to more readily exploit the potential benefits of the agent-based modelling paradigm.