Some challenges of calibrating differentiable agent-based models

Arnau Quera-Bofarull * 1 Joel Dyer * 1 2 Anisoara Calinescu 1 Michael Wooldridge 1

Abstractlevel, resulting in exact derivatives.
Despite recent progress, the challenges involved in building
and benefitting from differentiable ABMs remain underexplored, and there exists little guidance to practitioners
interested in implementing and exploiting differentiable
ABMs. The aim of this paper is therefore to discuss some
central challenges in applying AD to ABMs.

Agent-based models (ABMs) are a promising approach to modelling and reasoning about complex
systems, yet their application in practice is impeded by their complexity, discrete nature, and the
difficulty of performing parameter inference andoptimisation tasks. This in turn has sparked interest in the construction of differentiable ABMs
as a strategy for combatting these difficulties, yet
a number of challenges remain. In this paper,
we discuss and present experiments that highlightsome of these challenges, along with potential
solutions.

2. Challenges
2.1. Discrete randomness
The issue of differentiating through discrete structures is inherent in ABMs, which simulate discrete events, transitions,and interactions that are incompatible with conventional
AD. Initial efforts to implement AD within ABMs have primarily centred on transforming the ABMs’ discrete control
flow structure with continuous approximations (Andelfinger, 2021). Furthermore, the use of the Gumbel-Softmax(GS) reparametrisation trick (Jang et al., 2017) allows forthe differentiation of discrete randomness, and has been deployed effectively in epidemiological ABMs (Chopra et al.,2023). However, this approach does not provide an ideal
solution. While it allows for gradient calculations, GS does
not guarantee unbiased or low-variance gradients (Huijben
et al., 2022). Developing unbiased and lower variance methods such as StochasticAD (Arya et al., 2022) in the Julia
programming language (Bezanson et al., 2017) is currently
an active field of research, but we limit the scope of our
discussion here to GS-based methods for discrete ABMs.

1. Introduction
Agent-based models (ABMs, see Appendix A for a brief
overview) have gained considerable popularity across a
range of disciplines, due to their ability to accurately simulate complex systems at a granular level. While thesemodels offer unique advantages, their complexity presents
significant challenges, for example in terms of parameter
calibration (see e.g. Dyer et al., 2022a;b). For such tasks,
multiple factors contribute to their difficulty, including theintractability of the ABM’s likelihood function, and theoften black-box and non-differentiable nature of the ABM.
These drawbacks of ABMs have motivated research into theconstruction of differentiable ABMs (Chopra et al., 2023;
Quera-Bofarull, 2023), for example through the use of differentiable programming and by exploiting automatic differentiation (AD) frameworks. AD – a methodological cornerstone in machine learning, largely underpinning the successof deep learning paradigms due to its ability to accurately
compute derivatives within models – circumvents issues
present in alternative approaches to model differentiation by
applying the chain rule of differentiation at a computational

Despite the potential lack of robustness of GS, GS-based differentiable ABM implementations have shown great successin improving the calibration (Chopra et al., 2023) and sensitivity analyses (Quera-Bofarull et al., 2023) of ABMs. In theExperiments section below, we further show that gradientsobtained using the GS trick are robust enough to enable fast
and accurate Bayesian inference (see subsection 3.2).

*
Equal contribution 1 Department of Computer Science, University of Oxford 2 Institute for New Economic Thinking, Oxford. Correspondence to: Arnau Quera-Bofarull <arnau.querabofarull@cs.ox.ac.uk>, Joel Dyer <joel.dyer@cs.ox.ac.uk>.

2.2. Reverse- vs. Forward-mode AD
In Reverse-mode AD (RMAD), a computation graph must
be stored that records all operations performed within themodel, such that the gradients of the model outputs withrespect to the input parameters can be obtained. This contrasts with Forward-mode AD (FMAD), where the gradients

Published at the Differentiable Almost Everything Workshop of the
40 th International Conference on Machine Learning, Honolulu,
Hawaii, USA. July 2023. Copyright 2023 by the author(s).

1

10−1
Loss

Path-wise
horizon=0
Path-wise
horizon=1
Path-wise
horizon=2
Path-wise
horizon=100
Score

10−2

10−3

where p(· | θ) is the ABM’s likelihood function. The gradients of a differentiable ABM can then be exploited bygradient-assisted methods for minimising the objective in(1), by finding a Monte Carlo estimate of the expression

0

200

400
600
Epoch

800

∇ω Ez∼pω [L(z)] .

For differentiable ABMs, a Monte Carlo estimate of thisgradient can be obtained using the path-wise derivative via
reparametrisation tricks (Mohamed et al., 2020). In suchcases, derivatives of the form ∂xt /∂ωi will contribute to
the estimate. To properly benefit from access to the differentiable ABM’s gradients in these settings, it is criticalthat low-variance, low-bias Monte Carlo estimates of (3)are available. However, as we will demonstrate in subsection 3.1, naively estimating these gradients by accounting
for both (a) the explicit dependency of each xt on ωi , and(b) the implicit dependency on ωi , mediated by the x1:t−1that result from the recursive structure of ABMs, can result in unusable gradient estimates with prohibitively large
variances. Consequently, modifications to vanilla AD canbecome necessary, as illustrated in subsection 3.1.

1000

Figure 1. Training loss as a function of epochs for each gradientestimation method and gradient horizon value. Solid line shows
10-step moving average; dashed lines are without averaging.

are computed during the forward simulation. There are two
important computational considerations when comparing
FMAD vs. RMAD. The first is that the computational timeassociated with FMAD scales with the number of model
inputs, while that of RMAD scales with the number ofmodel outputs. In machine learning, the latter option ismore prevalent, since machine learning models often have
many more inputs than outputs. However, the computationgraph that RMAD must store in (often GPU) memory canbe extremely large, hindering the possibility of differentiating through large models. This is particularly pertinent
for ABMs: the size of the computation graph grows with
the number of agents and time-steps, which can pose a challenge to the use of RMAD for ABMs with a large numberof agents and time-steps.

3. Experiments
In this section, we present experiments on the use ofgradient-assisted calibration methods for two ABMs, where
each experiment serves to highlight different combinationsof the challenges described in Section 2.
While there exist many different gradient-assisted calibration methods, we focus on a variational approach toBayesian parameter inference termed Generalised Variational Inference (GVI, Knoblauch et al., 2022) – a
likelihood-free Bayesian inference approach that has previously been used to calibrate the parameters θ ∈ Rd of
a differentiable ABM (Quera-Bofarull et al., 2023). Here,a variational procedure targets a “generalised” posterior(Bissiri et al., 2016)

To address this, in subsection 3.2 we discuss a differentiating
strategy that alternates between FMAD and RMAD whencalibrating ABMs, and we apply it to an epidemiologicalsimulation involving 8 million agents.
2.3. Monte Carlo gradient estimation

πw,y (θ) ∝ e−w·ℓ(y,θ) π(θ),

Since ABMs are typically stochastic models, it can often be
the case that practitioners are interested in performing an
optimisation problem of the form
min Ez∼pω [L(z)] ,ω∈Ωmin Ex∼p(·|θ) [D(x, y)] ,(4)

where π(θ) is a prior distribution, ℓ(y, θ) is a loss function
capturing the compatibility between the observed data yand the behaviour of the ABM at parameter vector θ, andw > 0 is a hyperparameter. To target this posterior, wetrain a normalising flow qϕ with trainable parameters ϕto minimise the Kullback-Liebler divergence KL(qϕ ∥πw,y )from qϕ to πw,y , yielding the minimising parameters
ϕw,y,π = arg min w Eqφ [ℓ(y, θ)] + KL (qϕ ∥π) . (5)(1)

where pω ∈ {pω′ : ω ′ ∈ Ω} is a probability distribution
on some domain Z indexed by a parameter ω belonging
to some set Ω, and L : Z → R is a loss function. For
example, certain parameter calibration procedures that seek
to identify the parameters θ in some set Θ that minimisesome discrepancy D(·, y) between the model output x andreal-world data y can be cast in the form
θ∈Θ

(3)

ϕ

Further details are provided in Appendix D. Code to reproduce the results and perform GVI on differentiable ABMs
can be found at https://github.com/arnauqb/
blackbirds.

(2)
2

H=0

H=1

H=2

H = 100

Score Function

140
120

Counts

100
80
60
40
20
0

10−3

10−2

10−1

10−3

10−2

10−1

10−3

10−2

10−1

10−3

10−2

10−1

10−3

10−2

10−1

Standard Deviation of Gradient Estimate

Figure 2. Histograms of the standard deviations in estimates of ∂Eqφ ℓ(y, θ)/∂ϕj for different gradient horizons and the score-based
method across all j for the Brock & Hommes model.

3.1. The Brock & Hommes model

vanilla pathwise derivative (red curve). Indeed, we seein this case that the access to the simulator’s gradients appears to offer no improvement over the score-based gradient,shown with the purple curve and obtained as

The Brock & Hommes model (Brock & Hommes, 1998) is aheterogeneous agent model for the price xt ∈ R of an asset
over time 1 ≤ t ≤ T . At each time step, the agents in themodel subscribe to one of a set of J > 1 trading strategies,
each of which is characterised by a trend-following parameter gj and bias parameter bj , j ∈ {1, . . . , J}. Following
Dyer et al. (2022a), we note that the price xt may be written
as deterministic transformations ft of the input parametersθ = (g1 , . . . , gJ , b1 , . . . , bJ ), auxiliary parameters α, andstandard Normal random variables:
xt = ft (ϵ1 , . . . , ϵt , θ, α),ϵt ∼ N (0, 1).

∇ϕ Eqφ [ℓ(y, θ)] = Eqφ [ℓ(y, θ)∇ϕ log qϕ (θ)] .

Drawing inspiration from the literature on backpropagationthrough-time (e.g. truncated back-propagation in the context of RNN training, see Sutskever (2013)), we consider
pruning a subset of the paths in the computation graph that
contribute to each of the ∂xt /∂θ i as a possible solution
to this problem. We achieve this by invoking an appropriate stop gradient operation (e.g. .detach() inpytorch, Paszke et al. (2019)) on terms xt′ that (a) contribute directly/explicitly to the evaluation of xt and (b) for
which t′ < t − H for some “gradient horizon” H ≥ 0.
As evident from Figure 1, we observe that a finite gradient horizon can dramatically improve the gradient-assisted
training of qϕ . In this experiment, the best performance wasobserved while using a gradient horizon of H = 0.

(6)Further details are provided in Appendix B. Thus, providedℓ is chosen to be a differentiable function of θ, this enables
us to employ gradient-based approaches to minimising theobjective (5) that exploit the reparameterisation trick.
Fixing g1 = b1 = b4 = 0 and g4 = 1.01, we consider
the task of calibrating parameters g2 , g3 , b2 , b3 given synthetic data y = (y1 , . . . , yT ) generated from the model at(g2 , g3 , b2 , b3 ) = (0.9, 0.9, 0.2, −0.2) with T = 100. We
follow Cherief-Abdellatif & Alquier (2020) and target thegeneralised posterior (4) given by the choiceℓ(y, θ) = MMD2 (PT , Pθ ),(8)

We posit that this is a manifestation of a bias-variance tradeoff in the Monte Carlo gradient estimation step: pruning
a subset of paths in the computation graph with the use ofa finite gradient horizon may introduce some bias in, but
can significantly reduce the variance of, Monte Carlo estimates of the gradient ∇ϕ Eqφ [ℓ(y, θ)] when employing
the pathwise derivative. This hypothesis is supported by
Figure 2, which shows histograms of the standard deviationof the estimates of ∂Eqφ [ℓ(y, θ)] /∂ϕj across all j for gradient horizons H ∈ {0, 1, 2, 100}, and for the score-based
estimator. There, we see that the histogram shifts towards
larger values as H increases. Further results supporting thishypothesis are given in Appendix B.2.

(7)where MMD2 (PT , Pθ ) is the maximum mean discrepancy between the empirical measure of returns PT =(y1 , . . . , yT ) and the distribution Pθ of returns implied by
the simulator at parameters θ. Using a Gaussian kernel
within the MMD computation, the operations comprisingevaluation of ℓ(y, θ) are also all differentiable and deterministic, enabling evaluation of the term ∇ϕ Eθ∼qφ ℓ(y, θ)in (5) using the reparameterisation trick (see Appendix B.1).

3.2. The JUNE model

Despite our ability to compute the partial derivatives∂xt /∂ϕi exactly, the posterior estimator qϕ struggles totrain with a gradient-assisted approach to minimising (5).
This can be seen in Figure 1, in which the objective function
decreases slowly with the number of epochs when trainedwith AdamW (Loshchilov & Hutter, 2017) and using the

The J UNE model (Aylett-Bullock et al., 2020) is a largescale epidemiological ABM of England based on a realistic
synthetic population constructed from the English census.
Calibrating the original implementation required the construction of a surrogate model due to its high computational
3

101

Cumulative number of infections

Memory allocation[MB / (1,000 agents)]

102

Reverse diff.
Forward diff.

100

0

100
200
Number of time-steps

300

104
103

Trained Flow
Prior

10

20
30
Time-step

Untrained Flow
Ground Truth

40

50

since these evaluations are embarrassingly parallelisable,
the impact on performance can be minimal.
With the above in mind, we set up an experiment with theLondon’s population (8.1 million people) in G RADABMJ UNE. We generate a synthetic time-series of daily infections for 50 days using some assumed parameters that we
aim to recover through our calibration process. The parameters that we vary are the contact intensities at 10 different
locations, as well as the number of initial cases. Further
details of the experimental setup are shown in Appendix C.

3.2.1. R EDUCING MEMORY CONSUMPTION THROUGHFORWARD - MODE AD
G RADABM-J UNE can simulate the entire English population at a scale of 1:1 — 53 million people at the time of the2011 census. As discussed in subsection 2.2, differentiating
through this model using RMAD is challenging due to thehigh memory demand to store the computation graph. To
perform GVI for this model, we implement a hybrid AD
technique: we use FMAD to obtain the Jacobian Jθ of the
ABM outputs with respect to the ABM parameters, andcombine it with RMAD through the flow qϕ , yielding∂Eqφ ℓ(y, θ)∈ R1×d .
∂θ

105

Figure 4. Cumulative infections over time. Black: synthetic
ground truth; green, orange, blue: runs obtained from samples
from the prior, untrained flow, and trained flow, respectively.

cost (Vernon et al., 2022). The G RADABM-J UNE model(Quera-Bofarull, 2023) is a differentiable implementationof J UNE which employs the GS reparameterisation trick todifferentiate through discrete randomness. Compared to its
non-differentiable counterpart, G RADABM-J UNE has beenused to more efficiently generate parameter point estimates,
as well as sensitivity analyses (Quera-Bofarull et al., 2023).

Jθ (Eqφ ℓ(y, θ)) =106

0Figure 3. Memory consumption of running AD in the JUNE modelper 1,000 agents, for the Forward and Reverse mode of AD.

∇ϕ Eqφ [ℓ(y, θ)] = Jθ (Eqφ [ℓ(y, θ)]) · ∇ϕ θ,

107

We apply the GVI procedure to calibrate the J UNE modelwith 11 free parameters. The flow converges after approximately 3,000 model evaluations, highlighting the potential
for simulation-efficient calibration with gradient-assisted
methods. In Figure 4, we show a comparison of runs obtained by sampling the ABM parameters from the trained
flow, the untrained flow, and the prior. This demonstrates
that the trained flow generates parameters that result inclose agreement between the simulator and ground truth
data, while providing useful uncertainty quantification.

with (9)(10)

Here, (10) is the Jacobian obtained through FMAD and∇ϕ θ ∈ Rd×F is the gradient, obtained with RMAD, of thed ABM parameters generated by the normalising flow withparameters ϕ ∈ RF .

4. Conclusions and discussion
This study examines some challenges that arise from theapplication of vanilla AD to ABMs, such as overcoming theinherent discreteness of ABMs and the variance and highcomputational requirements of passing gradients throughlarge simulators. We have shown that these challenges canbe overcome to some extent with different modifications
to vanilla AD. As supporting evidence, we successfully
calibrate differentiable implementations of the Brock &Hommes and J UNE models with these modifications, thelatter involving over eight million agents and discrete randomness. In this way, this study helps to pave the way
towards robust calibration of large-scale agent-based models.

In Figure 3, we plot the memory costs of employing FMAD
and RMAD to compute the ABM’s Jacobian. We see thatthe cost of FMAD is independent of the number of timesteps, since no computation graph is stored. In contrast, thecost of RMAD scales linearly with the number of time-steps
and agents. Simulating the entire English population for 300
time-steps would require 5TB of memory, while doing sowith FMAD would require merely 18GB, regardless of thenumber of time-steps. The increase in computational timeof FMAD comes at an increase of computational cost since
it requires d evaluations of the model for θ ∈ Rd ; however,4References

Chopra, A., Rodrı́guez, A., Subramanian, J., Quera-Bofarull,A., Krishnamurthy, B., Prakash, B. A., and Raskar, R.
Differentiable agent-based epidemiology. In Proceedingsof the 2023 International Conference on Autonomous
Agents and Multiagent Systems, AAMAS ’23, pp. 1848–
1857, Richland, SC, 2023. International Foundation forAutonomous Agents and Multiagent Systems. ISBN 9781-4503-9432-1.

Andelfinger, P. Differentiable Agent-Based Simulation
for Gradient-Guided Simulation-Based Optimization.
arXiv:2103.12476 [cs, eess], March 2021.
Arya, G., Schauer, M., Schäfer, F., and Rackauckas, C. Automatic Differentiation of Programs with Discrete Randomness, October 2022.

Dyer, J., Cannon, P., Farmer, J. D., and Schmon, S. Blackbox Bayesian inference for economic agent-based models.
arXiv preprint arXiv:2202.00625, 2022a.

Aylett-Bullock, J., Cuesta-Lazaro, C., Quera-Bofarull, A.,
Icaza-Lizaola, M., Sedgewick, A., Truong, H., Curran, A.,
Elliott, E., Caulfield, T., Fong, K., Vernon, I., Williams, J.,
Bower, R., and Krauss, F. June: Open-source individualbased epidemiology simulation. Royal Society Open Science, 8(7):210506, 2020. doi: 10.1098/rsos.210506.

Dyer, J., Cannon, P., Farmer, J. D., and Schmon, S. M.
Calibrating Agent-based Models to Microdata with Graph
Neural Networks. In ICML 2022 Workshop AI for AgentBased Modelling, 2022b.

Aylett-Bullock, J., Cuesta-Lazaro, C., Quera-Bofarull, A.,
Katta, A., Pham, K. H., Hoover, B., Strobelt, H., Jimenez,
R. M., Sedgewick, A., Evers, E. S., Kennedy, D., Harlass,S., Maina, A. G. K., Hussien, A., and Luengo-Oroz, M.
Operational response simulation tool for epidemics within
refugee and IDP settlements: A scenario-based case studyof the Cox’s Bazar settlement. PLOS Computational
Biology, 17(10):e1009360, October 2021. ISSN 15537358. doi: 10.1371/journal.pcbi.1009360.

Huijben, I. A. M., Kool, W., Paulus, M. B., and van Sloun,
R. J. G. A Review of the Gumbel-max Trick and its
Extensions for Discrete Stochasticity in Machine Learning. IEEE Transactions on Pattern Analysis and MachineIntelligence, 45(2):1353–1371, 2022.
Jang, E., Gu, S., and Poole, B. Categorical Reparameterization with Gumbel-Softmax. arXiv:1611.01144 [cs, stat],August 2017.

Baptista, R., Farmer, J. D., Hinterschweiger, M., Low, K.,
Tang, D., and Uluc, A. Staff working paper no. 619
macroprudential policy in an agent-based model of theuk housing market. 2016.

Knoblauch, J., Jewson, J., and Damoulas, T.
An
optimization-centric view on bayes’ rule: Reviewing andgeneralizing variational inference. Journal of Machine
Learning Research, 23(132):1–109, 2022.

Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.
Julia: A fresh approach to numerical computing. SIAM
review, 59(1):65–98, 2017.

Loshchilov, I. and Hutter, F. Decoupled weight decay regularization. In International Conference on Learning
Representations, 2017.

Bissiri, P. G., Holmes, C., and Walker, S. G. A generalframework for updating belief distributions. Journal of
the Royal Statistical Society: Series B (Statistical Methodology), 78(5):1103, 2016.

Mohamed, S., Rosca, M., Figurnov, M., and Mnih, A.
Monte Carlo Gradient Estimation in Machine Learning.
The Journal of Machine Learning Research, 21(1):5183–
5244, 2020.

Brock, W. A. and Hommes, C. H. Heterogeneous beliefs androutes to chaos in a simple asset pricing model. Journal of
Economic Dynamics and Control, 22(8):1235–1274, July
1998. ISSN 0165-1889. doi: 10.1016/S0165-1889(98)00011-6.

Papamakarios, G., Pavlakou, T., and Murray, I. Masked
autoregressive flow for density estimation. In Guyon, I.,
Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), Advances in Neural
Information Processing Systems, volume 30. Curran Associates, Inc., 2017.

Cherief-Abdellatif, B.-E. and Alquier, P. MMD-Bayes:
Robust Bayesian Estimation via Maximum Mean Discrepancy.
In Zhang, C., Ruiz, F., Bui, T., Dieng, A. B., and Liang, D. (eds.), Proceedings of
The 2nd Symposium on Advances in Approximate
Bayesian Inference, volume 118 of Proceedings of Machine Learning Research, pp. 1–21. PMLR, 08 Dec
2020. URL https://proceedings.mlr.press/
v118/cherief-abdellatif20a.html.

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J.,
Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga,
L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang,
L., Bai, J., and Chintala, S. PyTorch: An Imperative
Style, High-Performance Deep Learning Library. Curran
Associates Inc., Red Hook, NY, USA, 2019.
5

Quera-Bofarull, A. Arnauqb/GradABM-JUNE: AAMAS
zenodo release. Zenodo, February 2023.
Quera-Bofarull, A., Chopra, A., Aylett-Bullock, J., CuestaLazaro, C., Calinescu, A., Raskar, R., and Wooldridge,
M. Don’t simulate twice: One-shot sensitivity analyses via automatic differentiation. In Proceedings of the2023 International Conference on Autonomous Agentsand Multiagent Systems, AAMAS ’23, pp. 1867–1876,
Richland, SC, 2023. International Foundation for Autonomous Agents and Multiagent Systems. ISBN 978-14503-9432-1.
Quera-Bofarull, A., Chopra, A., Calinescu, A., Wooldridge,
M., and Dyer, J. Bayesian calibration of differentiableagent-based models. ICLR Workshop on AI for Agentbased Modelling, 2023.
Stimper, V., Liu, D., Campbell, A., Berenz, V., Ryll, L.,
Schölkopf, B., and Hernández-Lobato, J. M. Normflows: A PyTorch Package for Normalizing Flows. arXiv
preprint arXiv:2302.12014, 2023.
Sutskever, I. Training Recurrent Neural Networks. University of Toronto Toronto, ON, Canada, 2013.
Vernon, I., Owen, J., Aylett-Bullock, J., Cuesta-Lazaro, C.,
Frawley, J., Quera-Bofarull, A., Sedgewick, A., Shi, D.,
Truong, H., Turner, M., Walker, J., Caulfield, T., Fong, K.,
and Krauss, F. Bayesian emulation and history matchingof JUNE. Philosophical Transactions of the Royal Society
A: Mathematical, Physical and Engineering Sciences,380(2233):20220039, October 2022. doi: 10.1098/rsta.
2022.0039.

6

A. Agent-based models
Agent-based modelling is the name given to a broad approach to modelling complex systems that consist of multiple discrete,autonomous, and heterogeneous interacting components – the “agents” of the system. Examples of such complex systems
include the housing market (Baptista et al., 2016), in which a large collection of renters, homeowners, financial institutions
etc. interact and take actions which affect, for example, the availability of housing and mortgage rates. An agent-based
approach to modelling such a system would model the system at the level of these individual agents in the system, oftenwith the intention of observing how aggregate, macroscopic properties of the system emerge from the microscopic details of
the system.
While this is often a natural approach to modelling systems of this kind, the inherently discrete nature of the model’s
components and dynamics give rise to difficulties in applying gradient-based optimisation and calibration techniques. We
expand on these difficulties in section 2.

B. The Brock & Hommes model
The dynamics of the Brock and Hommes model are often expressed as the following system of coupled equations:

xt =





J
X

1 
(gj xt−1 + bj ) nj,t + σϵt  , ϵt ∼ N (0, 1),
R j=1

nj,t = PJ

exp (βUj,t−1 )j ′ =1 exp (βUj ,t−1 )

,

(11)
(12)

′

Uj,t−1 = (xt−1 − Rxt−2 ) (gj xt−3 + bj − Rxt−2 ) ,(13)

where R, β, σ are auxiliary parameters. We fix J = 4, R = 1.01, σ = 0.04, g1 = b1 = b4 = 0, g4 = 1.01, and β = 120 for
the experiment presented in the main body of the paper. By rewriting the above system of equations, we are able to find thetransition density for observation xt+1 asp(xt+1 | x1:t , θ, α) = N f (xt−2:t , θ, α), σ 2 /R2



(14)

where
J

f (xt−2:t , θ, α) =1 X
exp [β (xt − Rxt−1 ) (gj xt−2 + bj − Rxt−1 )](gj xt + bj )PJ
R j=1 j ′ =1 exp [β (xt − Rxt−1 ) (gj ′ xt−2 + bj ′ − Rxt−1 )](15)and α = (R, β, σ). The model is taken to be initialised with x−2 = x−1 = x0 = 0.
B.1. The asset prices as differentiable and deterministic transformations of input noise
We claim in the main body that we may rewrite the xt as deterministic transformations of standard Normal random variables.
By exploiting the autoregressive structure of the model, we explicitly provide these forms for x1 and x2 below to demonstrate
this claim. Throughout, ϵt ∼ N (0, 1) are iid standard Normal random variables. They are as follows:J
1 X bj+ σϵ1  := f1 (ϵ1 , θ, α)x1 =R j=1 J
h
P
i
bj ′′
βb
J
J
exp Rjj ′′ =1 J + σϵ1
1 X
σ
h
P
i + ϵ2 := f2 (ϵ1:2 , θ, α).
x2 =
P
βb
b
′
′′
J
J
j
j
R j=1
R
j ′ =1 exp
j ′′ =1 J + σϵ1
R

(16)

(17)

Repeating this process, we find that the xt may all be expressed in the form xt = ft (ϵ1:t , θ, α) for a deterministic mappingft : Rt × R2J × R3 → R.
7

Takingℓ(y, θ) = MMD2 (PT , Pθ )

(18)
′

′

= Ex,x′ ∼Pθ [k(x, x )] + Ey,y′ ∼PT [k(y, y )] − 2Ex∼Pθ ,y∼PT [k(x, y)]
≈

X
X
1
1
2
k(xt , xt′ ) +k(yt , yt′ ) − 2
T (T − 1)
T
(T
−
1)
T
′
′
t̸=t

t̸=t

T
X

k(xt , yt′ )

(19)
(20)

t,t′ =1

with a Gaussian kernel k, the loss ℓ(y, θ) is a deterministic and differentiable transformation of the noise drawn from thebase distribution ρ of the normalising flow and of the separate noise source with distribution ν given as input to the simulator.
This permits us to estimate the gradient of the first term in (5) as∇ϕ Eqφ [ℓ(y, θ)] = ∇ϕ Eu∼ρ [ℓ(y, θ ϕ (u))](21)= Eu∼ρ [∇ϕ ℓ(y, θ ϕ (u))]
X1= Eu∼ρ ∇ϕ k (ft (ϵ1:t , θ ϕ (u), α), ft′ (ϵ1:t′ , θ ϕ (u), α)) −T (T − 1)
′

(22)

t̸=t


T
2 X
∇ϕ k(ft (ϵ1:t , θ ϕ (u), α), yt′ )
T2 ′

(23)

t,t =1

where in the first line we use the Law of the Unconscious Statistician and assume throughout that the order of derivatives
and integrals can be exchanged freely.
B.2. Further experimental results for the Brock & Hommes model
B.2.1. C ALIBRATION RESULTS WITH GRADIENT HORIZON H = 0
In Figure 5, we show the generalised posterior approximated by the converged normalising flow and with gradient horizon
H = 0, which achieved a loss close to 0. Since the objective function is lower-bounded by 0, this posterior can be taken to
be a good approximation to the generalised posterior it targets.
B.2.2. F URTHER EVIDENCE IN SUPPORT OF THE BIAS - VARIANCE TRADE - OFF IN THE REPARAMETERISED M ONTEC ARLO GRADIENT ESTIMATOR AT DIFFERENT GRADIENT HORIZONS
To further test the hypothesis that the Monte Carlo gradient estimators at different gradient horizons result in a bias-variance
trade-off that can result in favourable performance at finite gradient horizons, we inspect the empirical distribution of theestimates ηN of the gradient (23) based on N Monte Carlo samples,


T
N
X
2 X
1 X
1
(24)
∇ϕ k (xt , yt′ )
∇ϕ k (xt , xt′ ) − 2
ηN :=N n=1 T (T − 1)
T′′t,t =1t̸=twhere here we take a diagonal Gaussian distribution over R4 as the posterior estimator qϕ . In this experiment, therefore,
ϕ = (µ1 , . . . , µ4 , σ1 , . . . , σ4 ), where µi and σi are the mean and standard deviation in each dimension of this choice of qϕ .
When implemented in the form given by Equations (14) and (15) – as is necessary to avoid the cumbersome task of finding
the explicit form of ft (ϵ1:t , θ, α) for each t – the xt depend on θ both explicitly and implicitly via xt−3:t−1 . Thus in general
we have
∞

X
∂xt
=
∂θ i
l=1

X

Length l paths
v=(v0 ,v1 ,...,vl−1 ,vl )with v0 =xt ,vl =θ i

l−1
Y

∂vi
,
∂v
i+1
i=0

(25)

where we abuse notation by taking the derivative on the left-hand side to mean “holding only (θ 1 , . . . , θ i−1 , θ i+1 , . . . , θ d )constant” while the partial derivatives on the right-hand side are partial derivatives in the true sense of the term (or
8

0
0

2

0

2
−

2

0

2
−

−

1.
5
0.
0
1.
5
3.
0
−
2
0
2
4

−

2

b3

2

−

2

b2

2

−

g3
2 0 2 4

Flow
Untrained Flow

g2

g3

b2

b3

Figure 5. The inferred generalised posterior for the Brock & Hommes model using a gradient horizon H = 0.

equivalently the xt on the left-hand side is viewed only as a function of θ, while on the right-hand side they are viewed
as functions of both θ and x1:t ). Choosing a gradient horizon of H > 0 then amounts to retaining paths of length greaterthan 1 if their first edge corresponds to an edge connecting xt to any node in XH := {xt−h : h ∈ HH }, where H0 = ∅ and
Hj = {1, . . . , j} when j > 0. In this way, the derivative (25) is then taken instead as
∞

∂xt X
∂xt
=
+
∂θ i
∂θ i
l=2

X

Length l paths
v=(v0 ,v1 ,...,vl−1 ,vl )with v0 =xt ,v1 ∈XH ,vl =θ i

l−1
Y

∂vi
,
∂vi+1
i=0

(26)

where the same abuse of notation is once again used. This elimination of terms from the summation can be expected toreduce the variance since for two random variables X0 , X1 , it is the case that
Var(X0 + X1 ) = Var(X0 ) + Var(X1 ) + 2Cov(X0 , X1 ),(27)

which can be greater than Var(X0 ) if Var(X1 ) + 2Cov(X0 , X1 ) > 0. It may also be preferable to stricter truncation of thecomputation graph – for example, by pruning all paths beyond a certain length – as it retains information on long-range
dependencies while still potentially reducing variance.
We plot in Figure 6 boxplots for the distribution of ηN obtained with N = 5 and with different values of H, obtained at a
fixed value for ϕ (the results were qualitatively similar for different ϕ we tried, and so we show only the results from one
settings). We also show the same boxplots for the gradient estimate obtained with the score-based estimator,N∇ϕ Eqφ [ℓ(y, θ)] = Eqφ [ℓ(y, θ)∇ϕ log qϕ (θ)] ≈1 Xℓ(y, θ (n) )∇ϕ log qϕ (θ (n) ).
N n=1

(28)

Orange (green) dashed lines show the mean (median) of the distributions. The blue crosses show the mean of the distributionof the gradient estimate (28) obtained using N = 1000, which provide a good estimate of the target value (since thescore-based estimator is unbiased). We see from this that, generally speaking, the variance of these estimates increaseas H increases, while the bias in the estimates do not degrade substantially. This highlights the possibility that using a9

Gradient Estimate ηN

finite gradient horizon can be beneficial when performing Monte Carlo gradient estimation for differentiable time series
simulation models, such as ABMs, when reparameterisation is possible. Further work will be required to establish thegeneral applicability and suitability of this technique.
H=0

H=1

H=2

H = 100

Score functionµ 1 µ 2 µ 3 µ 4 σ1 σ2 σ3 σ4µ 1 µ 2 µ 3 µ 4 σ1 σ2 σ3 σ4µ 1 µ 2 µ 3 µ 4 σ1 σ2 σ3 σ4µ 1 µ 2 µ 3 µ 4 σ1 σ2 σ3 σ4µ 1 µ 2 µ 3 µ 4 σ1 σ2 σ3 σ4

0.4
0.2
0.0
−0.2
−0.4

Parameters

Figure 6. Boxplots for the empirical distribution of Monte Carlo gradient estimates from gradient horizons H = 0, 1, 2, 100, and from thescore-based gradient estimate, (28). Orange (green) dashed lines show the mean (median) of the distributions. Blue crosses show themean of the distribution of the gradient estimate (28) obtained using N = 1000.

C. The JUNE model
The JUNE model (Aylett-Bullock et al., 2020) is an agent-based epidemiological model that generates a synthetic population
at a highly detailed level using the English census data. This model has been applied in various scenarios, including
analyzing the impact of the first and second waves of SARS-CoV2 in England (Vernon et al., 2022) and devising strategies
to control disease transmission in refugee settlements (Aylett-Bullock et al., 2021).
To enhance its performance and enable gradient-based calibration, the JUNE model has been incorporated into the
G RADABM framework (Chopra et al., 2023; Quera-Bofarull et al., 2023). This integration allows for faster execution andmore efficient parameter calibration. The JUNE model offers a wide range of configurable parameters related to diseasetransmission and progression, vaccination, and non-pharmaceutical interventions.
Given a susceptible agent exposed to an infection at location L, the probability of that agent getting infected is given bX
p = 1 − exp −ψs βL ∆t
Ii (t),(29)
i∈gwhere the summation is conducted over all contacts an agent has with infected individuals at the given location L. The term
Ii (t) represents the time-dependent infectious profile of each infected agent, while ∆t is the duration of the interaction.
Additionally, βL corresponds to a location-specific parameter that captures the variation in the nature of interactions acrossdifferent locations.
Since the βL parameters are not directly measurable physical quantities, they are typically calibrated using available dataon the number of cases or fatalities over a specific time period. For the current work we consider the calibration of 11 βLparameters corresponding to the contact intensity at households, companies, schools, universities, pubs and restaurants,
gyms, cinemas, shops, care homes, and residence visits. Additionally, we also calibrate the initial number of infections, I0 ,which are distributed randomly across the population. The synthetic ground truth data is generated by using I0 = 10−3.5 Na
where Na is the number of agents, and βhousehold = βcarehome = 0.6, βschool = βcompany = βuniversity = 0.4, βpub =βshop = βgym = βcinema = βvisit = 0.1.
The normalising flow is trained setting ℓ(y, θ) to be the squared distance between the log10 of the infection time-series.
This choice of loss function keeps the training robust against outliers, since the number of infections can oscillate between
several orders of magnitude. The specific training parameters are described in Appendix D.
C.1. Further experimental results for the JUNE model
We show in Figure 7 the loss as a function of epoch when performing SVI with G RADABM-J UNE. We observe a rapid
convergence after 600 epochs. Since we are sampling 5 Monte Carlo samples to estimate Equation 3, this results in 3,000
10model evaluations. We also make a corner plot of the train and untrained normalising flow which we show in Figure 8,
where the solid black line denotes the prior density. We observe that the flow is very confident about the value of moresensitive parameters such as the initial number of infections and the contact intensity at companies, while it is less certain forvenues which have a low impact in the overall number of infections, such as cinemas. It is worth noting that this calibration
challenge is very underdetermined, that is, it is quite difficult by just observing the overall number of infections over time toinfer the contact intensities at each location. Nonetheless, the flow fits well the synthetic ground truth data.

D. Further experimental details
We use the NORMALIZING - FLOWS library (Stimper et al., 2023) to implement the normalising flows in PyTorch. All models
are trained using the AdamW optimizer (Loshchilov & Hutter, 2017) with a learning rate of 10−3 .
To calibrate the Brock & Hommes and J UNE models, we employ a masked affine autoregressive flow (Papamakarios et al.,2017) with 16 transformations, each parametrized by 2 blocks with 20 hidden units. We also set the regularisation weight to
w = 10−3 for both models and estimate Equation 3 using 5 Monte Carlo samples.

Total Loss

100

10−1

10−2

0

200

400
Epoch

600

800

Figure 7. Loss function per epoch for the J UNE calibration.

11

I0
β household
β company
β school
β university
β pub

12
β grocery
β gym
β cinema
β visitβ care home

Figure 8. Corner plot of samples from the untrained normalising flow, the trained normalising flow, and the prior for the J UNE model
calibration.
0

6

0.

6
2
4
6
8

0.
0.
0.

3

0.

0.
0.

0

0.

3

3

0.

0.

0

0.

3

8

1.
0.

6

0.

−

−

−

−

8

5
0
5
0

0.
1.

6

0.

0.

3

0.

0.

0

0.

3

4

0.

0.

4
0

0.

4

0.

0.

0

0.

4

4
8

0.

0.

0

0.

4

0.

0.

−

5

0.

0.
−

0

0.

1.
− 0
1.
0
−
0.
5

−

−

4

0.

2

0.

0.

0

1

2

−
−

3

−

−
5

0.
0

0.
5

β care home
0.

1.
0

−
0.
3

0.
0

0.
3

β visit
0.
6

−
0.

0.
0

0.
4

β cinema
4

−
0.
4

0.
0

β gym
0.
4

−
0.
8

−
0.

0.
0

0.
4

β grocery
4

0.
8

−
1.
0

−
0.
5

0.
0

0.

β pub
5

0.
2

0.
4

0.
6

0.
8

β university
1.
0

−
0.
3

0.
0

0.
3

β school
0.
6

−
0.
3

0.
0

0.
3

0.

β company
6

0.
2

0.
4

0.
6

0.
8

β household
1.
0