more when the spreads are narrow and the volume is high. These agents have not
learned intuitive behaviour.

Speculating as to why this is the case: when $s=1$, being 60% of the spread
states, we have that these states get returned to more often and the agent can
learn optimal actions in these states better than in the other states.
Moreover, the agents are learning to trade less in all the states over the
training period so it could be the case that the agents are visiting these
states more often and learning to reduce their volume. While the wider spread
states are visited less often and are going to be learning to reduce the
volume at a slower rate, which is what we see in the plots. This effect does
occur in the smaller state space agents but it is somewhat mitigated. Once
again the agents with the larger state space should be trained for longer so
that the agent can visit each state-action pair more often. The plot for the
third agent with the larger state space can be seen in Figure 41(b), for
similar results.

Figure 13: To compare the learning dynamics of the agents with the larger
state space to the agents with the smaller state space we plot the average
action for the high (resp. low) narrow (resp. wide) spread states. The agents
with the larger state space have not been able to learn the intuitive actions
as they trade significantly more when the volumes are low and the spreads are
wide. As noted previously, due to the large state space the agents are not
able to return as often to the existing states and update the actions
therefore the agents find it more difficult to learn in the episodes. With a
much longer training period, we speculate that these agents would be able to
learn intuitive behaviour. The plot for the third agent with the larger state
space can be seen in Figure 41(b), for similar results.

To understand how the agents have learned to trade for different time and
inventory states we chose to see how the agents’ actions differ when they
trade at the beginning, with a large amount of inventory, as compared to when
they are trading at the end with less inventory. States with a lot of, and a
little time, remaining are classified respectively as states that have $t\geq
4$ and $t\leq 2$. States with a lot of, and a little inventory remaining have
$i\geq 4$ and $i\leq 2$ respectively. Figure 14 plots the results for the two
agents with the smaller state space, the third agent’s results can be found in
Figure 42(a) in the appendix. It is clear that the agent has learned to trade
more inventory when it is early in the simulation and the agent has a large
amount of inventory when compared with trading later when it has less
inventory. This effect is present no matter how much initial inventory the
agent was given. This finding corroborates the results seen in Figure 10.

This effect may be due to the agents trading higher voloumes initially in
order to guarantee the execution of all the inventory, but then reducing the
volume later to reduce the price impact. It must also be noted that the
average action is rapidly decreasing when there is less time and inventory and
if the episodes were increased it is likely the agent would continue to
decrease the amount of volume traded in these states.

Figure 14: To determine the average behaviour in different time and inventory
state combinations we plot the average action taken by the greedy policy when
the agent has lots of inventory at the beginning of the simulation and when
the agent has less inventory at the end of the simulation. The agent has
learned to trade more when it is early in the simulation and the agent has a
large amount of inventory. This effect is present when there is a low initial
inventory (left) and when there is a high initial inventory (right). Figure
42(a), plots the results for the other smaller state space agent for similar
results.

The time and inventory interaction results for the two agents with the larger
state space are shown in Figure 15, with the third’s results being shown in
Figure 42(b) in the appendix. States with a lot, and a little time, remaining
are classified as states that respectively have $t\geq 7$ and $t\leq 4$; and
states with a lot, and a little inventory remaining have $i\geq 7$ and $i\leq
4$, respectively. Once again we see that both agents have learned to decrease
the amount of inventory traded as the agents are trained. However, the same
effect, of trading more initially when the agent has a large inventory is not
seen in the agents with the larger state space. The agents are trading roughly
the same amounts in each state combination.

This reason for this discrepancy between the state spaces can be seen by
plotting the marginal effects of time and inventory on the average actions.
Figures 43 and 44 in the appendix plot these results. Here we see that all the
agents have learned to trade, on average, more when there is less time
remaining and trade more when there is less inventory remaining. However, in
the smaller state space, we have that the time effect is much less than the
inventory effect, which translates to the agents trading more when more time
and inventory remains. However, for the agents with the larger state space,
the time and inventory effects are roughly similar, therefore the effects are
cancelled out in the interaction plot.

Figure 15: Plots of the average action were taken by two of the larger state
space agents for a different time and inventory state combinations. The effect
seen in the smaller state space (Figure 14) where the agents traded more
initially when the agent has lots of time and inventory remaining is not seen
in this figure. The agents are trading roughly the same amount in each state
combination. The reason for this discrepancy between the state spaces can be
seen by plotting the marginal effects of time and inventory on the average
actions. Figures 43 and 44 show that all the agents have learned to trade, on
average, more when there is less time remaining and trade more when there is
more inventory remaining. However, in the large state space, the effects
cancel out whereas in the smaller state space the inventory effect dominates.
The third larger state space agent’s similar results are shown in Figure
42(b).

###  6.4 Stylised fact dynamics

Every 100 episodes, including the first and the last, all the raw messages of
the simulation were recorded. The stylised fact moments were computed at each
of these points so that we could see how adding in the RL agents affected the
stylised facts of the ABM, and if these effects changed as the agents learned
to trade.

Four key moments were computed on the both the JSE, and the ABM, and are
plotted in Figure 16 along with their associated confidence intervals. This
includes the moments for each of the RL agents over the episodes. The GPH
statistic was roughly the same for the ABM and the JSE, and adding an RL agent
to the simulation has not made a significant difference. This is also true for
the GARCH moment, except for one erroneous point at the start of training for
one of the agents. For the Hill estimator, the confidence intervals were close
but did not overlap and the moments for the ABM’s with the RL agents have Hill
estimators that fall into the ABM confidence interval.

The Hurst exponent moment was the only moment that showed slight differences
between the ABM with and without the RL agent. For the RL simulations, the
Hurst exponent is lower than that of the ABM and consistently falls below the
lower confidence interval of the ABM. The only ones to fall into the ABM’s
confidence interval are the RL agents that had the smallest parent orders.
Indicating that the larger the amount of inventory traded by the RL agent, the
more this increases the amount of mean reversion in the simulation. As the
episodes increase there is a slight downward trend in the Hurst exponent for
the agents with larger parent orders. Figure 16 only shows four out of the
eight moments. The results for the other moments can be seen in Figure 45, and
they are similar to the results of the GPH, Hill, and GARCH estimators.

Figure 16: Compares the ABM, and JSE moments to the moments computed on data
generated by the simulations of the ABM with the RL agents. During the
training of the RL agents, the moments were computed every 100 episodes so
that as the agents trained we could also see if this generated an effect on
the moments. The GPH statistic was roughly the same for the ABM and the JSE,
and adding an RL agent to the simulation has not made a significant
difference. This is also true for the GARCH moment, except for one erroneous
point at the start of training for one of the agents. For the Hill estimator,
the confidence intervals were close but did not overlap and the moments for
the ABMs with the RL agents have Hill estimators that fall into the ABM
confidence interval. The Hurst exponent moment was the only moment that showed
slight differences between the ABM with and without the RL agent. For the RL
simulations, the Hurst exponent is lower than that of the ABM and consistently
falls below the lower confidence interval of the ABM. The only ones to fall
into the original ABM’s confidence interval are the RL agents that had the
smallest parent order, indicating that the larger the amount of inventory
traded by the RL agents the more this increases the amount of mean reversion
in the simulation. This plot only shows four out of the eight moments. The
results for the other moments can be seen in Figure 45, and they are similar
to the results of the GPH, Hill, and GARCH estimators.

##  7 Stability, realism and learnability

###  7.1 e-ABM

Even though the objective surface was difficult to optimise the calibration
was able to find reasonable parameters. As argued by Fabretti [18], a
successful calibration would mean that we need to show that the simulated data
comes from the same, or at least similar distribution as the empirical data.
We did this by comparing the moments of the simulated price path to the
empirical price path and found that, except for the KS statistic, the moments
of the simulated price path are similar to those found on the JSE (See Table
4). We saw that four out of the eight moments had overlapping confidence
intervals. While there are clear differences between the distributions, the
simulated data reasonably resembles the empirical data.

While acknowledging the challenges with a stylised fact based comparison, we
found that the model was able to correctly reproduce many of the stylised
facts observed in the market. Namely, Fat-tails, volatility clustering,
significant negative autocorrelations for only the first few lags in the log-
returns, decay in the trade-sign autocorrelation, and similarity in the
extreme-value distributions of the ABM and the JSE. Even though the model was
able to replicate the decay in the autocorrelation of the trade-signs, we did
find that the decay did not exhibit the long-memory property that was found on
the JSE. This can be seen in Figure 23. We know that this decay is due to two
effects: herding and order-splitting.

The e-ABM did not have any order-splitting agents therefore it is not
surprising to see that this stylised fact has not been fully reproduced.
Additionally, there is a large herding effect caused mainly by the chartists
all submitting approximately the same orders at the same time, and it was also
observed that they change their trades from buys to sells very slowly. The RL
agent is an order-splitting agent by construction, so to investigate whether
this agent affected the decay of the trade-sign autocorrelation, we plot the
trade-sign autocorrelation decay for all six of the RL agents in Figure 17. We
see that adding a single learning agent that performs order-splitting, does
not guarantee the correct decay in the trade-sign autocorrelation. However,
agents plotted in green, blue and pink, all illustrate that power-law decay is
possible using this model. After a few initial lags, the decay of
autocorrelation is approximately linear on the log-linear scale indicating
that the decay is moving closer to a power-law relationship.

Figure 17: Plots the trade-sign autocorrelation for the ABM with all the
different RL agents added. The rows represent the state space sizes and the
columns represent the initial inventories. We see that adding a single
learning agent that performs order-splitting, does not guarantee the correct
decay in the trade-sign autocorrelation. However, agents plotted in green,
blue and pink, all illustrate that adding a strategic order-splitting agent
can reduce the trade-sign autocorrelations and move the model closer to power-
law decay. After a few initial lags, the decay of autocorrelation is
approximately linear on the log-linear scale indicating that the decay is
moving closer to a power-law relationship. This figure demonstrates that
adding more RL agents and explicitly calibrating for the trade-sign
autocorrelation effect may be able to produce the correct decay.

The absolute log-returns are compared in Figure 18 which shows those for the
e-ABM (Figure 18(a)), those found in the real world measurements from the JSE
(Figure 18(b)) as compared to the simulations of the e-ABM with a single RL
agent over different initial inventory and state space sizes (Figure 18(c)).
We can see from this figure that the autocorrelations on the JSE decay far
slower than the simulated moments, suggesting that there is higher volatility
clustering found in the real market when compared with the simulated markets.
We do not see significant differences between the autocorrelations of the
e-ABM and the e-ABMs with the single RL agent added.

(a)

(b)

(c)

Figure 18: Plots the autocorrelation in the absolute log-returns for the
calibrated e-ABM in Figure 18(a) and the JSE in Figure 18(b), as well as for
all the models that include the RL agents in Figure 18(c). Comparing the ACF
for the JSE to rest, we see that the ACF found on the JSE decays much slower
indicating that the JSE exhibits a stronger volatility clustering signature.
We note that the addition of agents engaging in strategic order-splitting can
boost the autocorrelation in the absolute returns as seen in Figure 18(c) on
the right.

Changing the stochastic time model of Jericevich et al. [34] to an event-based
model does change the stylised facts. The first difference was that the e-ABM
model was able to increase the variance in the returns and reduce the amount
of density at zero returns. This made the log-return distribution more similar
to the empirical distribution. The e-ABM model reduced the difference between
the tail-index parameters of the simulated and empirical returns observed in
the e-ABM when compared to Jericevich et al. [34]. There are clear differences
in the autocorrelation of the returns. Jericevich et al. [34] found
significant, positive and negative lags, for a larger number of lags, when
compared with the e-ABM model. We only observed significant negative
correlations for the first few lags. The log-return autocorrelations in the
e-ABM model are thus more similar to the empirical data. The reason for this
is likely because we have multiple fundamentalist agents who all trade mean
reversion strategies, rather than just the one that was implemented in
Jericevich et al. [34]. Both models were able to show persistence in the
absolute log-return autocorrelations indicating volatility clustering, fat-
tails in the log-returns, as well as correct depth profiles. Additionally,
Jericevich et al. [34] was able to produce more accurate decay in the trade-
sign autocorrelation.

Price impact curves for the e-ABMs that include the RL agents are shown in
Figure 19 with the separation between buyer initiated (Figure 19(a)) and
seller initiated (Figure 19(b)) trades. The linear upwards trend in the price
impact on the log-log scale is indicative of a power law price impact
function, where the price impact increases with the volume of the trade.
Comparing these figures to Figure 3, we see that adding in the RL agents does
not change the power law price impact seen in the underlying e-ABM. Even
though the RL agents are all selling agents the power law price impact is
demonstrated for both buyer and seller initiated trades. Additionally, we
observe the price impact increase as the inventory traded increases. These
price impact curves are consistent with previous research [28, 32, 41, 34].

(a)

(b)

Figure 19: Plots the price impact curves for the e-ABMs that include the RL
agents, for buyer initiated (Figure 19(a)) and seller initiated (Figure 19(b))
trades. The linear upwards trend in the price impact on the log-log scale is
indicative of a power law price impact function, where the price impact
increases with the volume of the trade. Comparing these figures to Figure 3,
we see that adding in the RL agents does not change the slope of the power law
price impact function relative to that seen in the underlying e-ABM. The
observed price impact increases as the inventory traded increases and these
price impact curves remain consistent with previous research [28, 32, 41, 34].

###  7.2 Reinforcement learning

Each RL agent learned to increase its profit and decrease the variance in the
profit from an initial random agent to a reasonably optimised agent, and this
is shown in Figure 4. The liquidity effect caused by increasing the size of
the parent order is evident, as the agents with larger initial inventories
increase the frequency and severity of the liquidity shocks. This is
consistent with the power-law price impact effect seen in the underlying ABM.

What could be seen as a surprising result is that none of the agents were able
to significantly improve on the amount of profit received over trading an
initial single order with the value: $X_{0}\times m_{0}$. In previous
implementations of RL to solve the optimal execution problem the agents are
trained offline with the assumption that the order book is resilient. Meaning
that when the agent trades and pushes the price down, this does not affect the
next price observed. This assumption has to be made because the agents are
trained offline using static data sets, and the data was historic and used to
generate a historic simulation of changes in the prices. This was because in a
static data set the next observed price, that in the real world would be
caused by that trade, is not easily simulated in a way that captures all the
other agents’ reactions to that specific RL agent’s specific trade. In this
ABM, we train the RL agent online using a dynamic data set, and we have that
the trade from an RL agent will cause other agents to react to this trading
decision. We speculate that this may be one of the reasons why we see that the
RL agents are then unable to increase the profit more than $X_{0}\times m_{0}$
as there is a permanent contribution to price impact. This can be tested by
taking a static simulated price path, such as the one used in prior work,
training the same RL agent with the same implementation as seen here, and
determining if it can improve upon $X_{0}\times m_{0}$. Consider the
$\epsilon$ decay function (Figure 1) and the convergence of the Q-matrix (left
of Figure 6) and notice that Q-matrix differences stop improving suddenly at a
similar time to when $\epsilon$ is at its lowest ($\epsilon=0.01$). Agents
with the highest $X_{0}$ show this effect most prominently. The rate of
decrease in the one-step difference of the Q-matrix is much faster before
$\epsilon$ is at its lowest indicating that the agents’ Q-matrix was
converging rather quickly.

It is reasonable to speculate that having more training episodes and a slower
decay in $\epsilon$, will promote greater convergence in the Q-matrix which
will indicate that the agent is learning a better policy, which could increase
the profit beyond $X_{0}\times m_{0}$. This should be the case for the larger
state space agent but Figure 6 suggests this could also be the case for the
smaller state space agents. Furthermore, we saw that the policies are changing
very slowly and to see significant behaviour changes we will have to
significantly increase the training time for the RL agents.

In all the learning dynamics we saw that the agents initially increased the
the inventory it traded, as shown in Figures (12, 13, 41(a), 41(b), 14, 15,
42(a), 42(b), 43, and 44). This is likely due to the initial reward signal,
and as the agent learned it, they started to significantly decrease the
inventory traded. This is evidence that the agent is learning about price
impact and that it must reduce the inventory traded if it is to maximise its
profit. When considering the smaller state space agent, the decreases for the
most part have not converged, demonstrating that the agents have not converged
to stable policies and are still trying to improve the return. This is
additional evidence that the agents could be trained for a longer period, and
that this could change their behaviour and improve their return.

All RL agents learned to trade more when they have more inventory and trade
more when they are at the end of the simulation, as evidenced by Figures 43
and 44. These are intuitive actions because if you have a large amount of
inventory you will want to increase your rate of trading to ensure that you
fully execute the parent order. This is also true for when there is little
time remaining then on average you will want to increase the inventory traded
to increase the chances of fully executing the parent order.

The agents’ dynamical behaviour change from the start to the end of the
simulation and this is shown in Figures 14 and 15. These confirm that the
agent will start with high time and inventory states, and move to the lower
states as time passes and they trade inventory. For the smaller state space
agent, we see large differences in the behaviour of the agent at the start and
end of the simulation, as the agent has learned to trade significantly more
initially and reduces the amount traded at the end of the simulation. However,
the larger state space agent does not exhibit a change in behaviour as the
simulation progresses. The smaller state space agents have learned the correct
trading behaviour when the initial inventory is high, but it is still learning
this effect for the lower inventory levels. However, since convergence in the
policy has not occurred it is speculated that the agent will learn this if the
episodes are increased. The larger state space agents have not learned the
correct behaviour.

This shows that smaller and larger state space agents have drastically
different behaviours. The remaining question is why is this not shown as
different profits in the return function? A potential reason for this is that
each RL agent’s profits are being washed out by the other agents in the ABM
because they are reacting to their trades, and even though the agents’
behaviours are different they are not different enough to overcome this
effect.

Figures 16 and 45, demonstrate that adding an RL agent to the ABM does not
cause the moments to differ too far away from the underlying ABM’s moments.
Illustrating that the underlying ABM’s moments are robust to the addition of a
single RL agent.

###  7.3 Limitations and future work

The model of Jericevich et al. [34] faced several challenges, here we found
similar issues. From the depth profiles, shown in Figure 24(b), we see that
most of the volume is placed at the best bid/ask, which does conform well to
what is observed empirically [34, 55]. This is the result of using a Gamma
distribution to generate the price placements [34]. When the volumes are
generated from a power-law distribution we find that either the LTs are unable
to generate price moves due to the large volumes at the best bid/ask, or they
quickly clear out one side of the LOB, generating liquidity shocks and extreme
price movements. This can make the ABM fragile and difficult to optimise as
different runs with the same parameters then produce very different price
paths. This effect is exacerbated by the increased number of orders submitted
by the agents. It was observed that a single order could have a volume of up
to $1\times 10^{7}$, for either LPs or LTs.

Looking at the empirical evidence [9] we could generate price placement as a
power-law, and then generate the volume conditional on the price. The
conditional distribution of the volume based on the price is said to be
independent of the price for the first 20-50 ticks and then decay as a power-
law Bouchaud et al. [9]. More in-depth analysis of the price placement, limit
order volume distributions, and trade distributions are necessary, especially
for the JSE if we wish to be drawn from realistic distributions that match the
market we are going to calibrate to. This did not detract from the overall
intention of this work, which was aimed to demonstrate the impact of strategic
order-splitting.

As discussed previously, we have observed that the decay in the trade-signs
does not completely match what is observed in the market. Indicating that the
agents are submitting similar trades, and for longer periods of time, than
what they would do in the real market. This rule should be changed to give
greater heterogeneity to the agents so that they don’t submit similar orders
at similar times, and they should be programmed to change their orders at a
quicker rate. However, here we show that by adding order-splitting agents,
such as optimally executing RL agents, we can produce more realistic decay
rates.

A further limitation to this model that can be improved, is to increase the
number of episodes for training the RL agents. We have demonstrated that for
the smaller state space agents, the agents are still decreasing the average
inventory traded and there is evidence that increasing the number of episodes
would make consecutive Q-matrices more similar, thus indicating the
convergence in Q-matrix. Moreover, we have observed the policies changing at a
slow rate and significantly increasing the episodes will increase the chances
of the agent improving the return beyond that of a single large order of
value: $m_{0}\times X_{0}$.

Given that the market is made up of agents with bounded rationality it would
be interesting to determine how the model changes when we increase the number
of learning agents. This would be the first investigation into adding learning
agents into an ABM and could provide additional features to the model as that
could increase its realism. It may be able to determine what effects we see in
the market are due to adaptation and learning vs those that emerge just
through the interaction of simple rule based agents.

Adding additional adaptive agents into the market may decrease the stationary
properties of the simulation making it harder for each individual RL agent to
learn. This Multi-agent RL (MARL) environment could be a good test bed to see
how adding additional agents affects the learning dynamics of the RL agents,
when compared with just adding a single agent.

The RL agents weren’t able to improve upon $m_{0}\times X_{0}$, however, more
complex RL agents would be able to learn more complex behaviour that could
allow them to increase the return. Another application of this model could be
to compare the performance of more complex RL agents in an online manner.
However, we have demonstrated that the more complex the agent the longer you
will have to train it for, therefore, this will have to be taken into account
when adding complexity to the agents.

##  8 Conclusion

This work combines a realistic but minimally intelligent agent-based model
with a learning agent with bounded rationality in event time. This required
moving the market clearing mechanism into one that is event-driven so that the
minimally intelligent (zero-intelligence) agents can interact with the
learning agent in a truly asynchronous continuous-double auction market [35]
without the need for batched market clearing events to be externally imposed
in calendar-time for implementation synchronisation. The key challenge for
this type of implementation is the requirement of moving the agent-based model
(ABM) into an event-driven agent-based model (e-ABM) framework that is both
reactive and reasonably realistic (Section 2).

The e-ABM was built on two class of agents: liquidity takers who only submit
market orders, and liquidity suppliers who only submit limit-orders. The
liquidity takers were split in turn into two sub-classes: fundamentalists and
chartists. The fundamentalist started each trading session, such as a trading
day, with different notions of fair, or fundamental value, and traded
accordingly; while the chartists considered a moving average of past price
changes in order to trigger trading decisions. In order to adapt the hybrid
ABM the e-ABM model includes a spread condition as a cost to trading, and the
software implementation was moved into reactive framework (Sections 3.1,
3.2.1, and 3.2.2). Our implementation adapts the model of Jericevich et al.
[34] and again explicitly excludes direct calibration to neither the order-
flow, or trade-sign autocorrelations, nor the shape of the price impact.

We were able to replicate key stylised facts of the market as demonstrated in
Section 5 and D. The real market observations were based on trade and order-
book data extracted from the Johannesburg Stock Exchange (JSE) as used by
Jericevich et al. [34]171717Naspers trade-and-quote data for the 08-07-2019.
The dataset used to extract the day of trading data can be found here.. The
recovery of the stylised facts are evidenced by Figures 20 through 22 where,
respectively, the log-returns exhibit fat-tails, similar extreme value
distributions to those observed in the market, volatility clustering, as well
as replicating the ACF of the log-returns. The model was also able to capture
the decay in the auto-correlation of the trade-signs, but with has a slower
decay rate when compared with the JSE as evidenced by Figure 23. The JSE
exhibited the long-memory property described by Lillo and Farmer [40],
however, this was not replicated by the e-ABM. The decay of the auto-
correlation in the absolute values of the log-returns was slower in the real
market than the model. However, the model was able to replicate reasonable
price impact functions and depth profiles shown in Figures 3 and 24(b). It was
speculated that the introduction of a learning agent would increase the
autocorrelations in the absolute values of the price fluctuations (sometime
called the leverage effect), and reduce the order-flow autocorrelations, both
due to the introduction of strategic order-splitting.

The learning agent with bounded rationality was represented by simple tabular
Q-learning algorithm that was learning to solve an optimal execution selling
problem using only market orders. The algorithm considered two market states
based on features representing the limit-order book, the spread and the best
bid volume, and then two states representing the current execution agents
states, the remaining time available in the trading session, and the agents
remaining inventory. The RL agent attempted to learn how to modify an existing
TWAP trading schedule to optimally execute its sell order.

After the calibration of the e-ABM the single RL optimal execution agent with
its sell order was added into the market to interact with the e-ABM. The RL
agent then attempts to modify its event-based TWAP strategy. We tested this
agent with different parent order sizes and different state space sizes to
observe whether any agents’ profits could be found to increase and converge,
to improve a random initial agent configuration. It was found that all agents
were able to increase the profits over their initial random agent,
demonstrating that each agent has learned. However, it was noted that there
was an upper bound on the profits obtained. No agent was able to learn how to
improve upon $m_{0}\times X_{0}$, meaning that the agents were able to reduce
the implementation shortfall, but they were not able to increase profits
enough to observe a consistent and significant negative implementation
shortfall. Adding RL agents with increasing order sizes caused the ABM to have
significant liquidity effects as expected.

There were clear differences in the dynamics for agents with different state
space sizes. The smaller state space agents had the number of states they
visited converge much faster than the larger state space agents. They quickly
finished the state discovery process while the larger state space agent was
still rapidly discovering new states, which limited their ability to learn.
The larger state space agents showed less switching of actions in each state
when compared with the small state space agents, which was due to the agents
not returning to individual states as often to update their state-action
values.

This showed that the larger state space agents had difficulty learning in such
a small number of episodes. This was also evident in the behaviour learned. If
a large state space agent had a large amount of inventory at the beginning of
the simulation versus a smaller amount at the end, the agent did not learn to
change behaviour while the smaller state space agent learned to trade more
initially and then reduce the rate of trading at the end. There was also
evidence that the smaller state space agent was able to start learning to
trade intuitively using the spread and volume states. Given this, all agents
learned to reduce the number of trades and guarantee the execution of the
entire inventory, and they learned some structure during the training process
as evidenced by the change in the action volume trajectories over the training
period.

The statistical analysis work is based on the comparison of the calibrated
e-ABM via its simulated moments against the same e-ABM in the presence of the
single RL optimal execution selling agent with varying initial order sizes.
The e-ABM was calibrated using a simulated method of moments where the
optimisations are implemented using Nelder-Mead with Threshold Acceptance [25]
approach, using adaptive parameters [20]. This method of calibration was
chosen to balance computational speed with stability F.4 and was informed by
prior work [25, 20, 34]. The calibrated parameters are provided with
indicative confidence intervals in Table 3 following [34]. We see that the
most significant parameters are, in decreasing order:
$\sigma_{f},\kappa,N^{c}_{{}_{\mathrm{LT}}},\nu,N^{f}_{{}_{\mathrm{LT}}},$ and
$\delta$.

Testing the effect the RL agents had on the moments of the ABM, we plotted the
moments of the e-ABM that included the RL agents over the training period and
compared them to the ABM’s moments. We saw that, with the exception of the
Hurst exponent, the moments were reasonably robust to the addition of a single
RL agent even with increasing initial order sizes.

Future work will include both optimal buying and selling execution agents.
This will again be carried out by testing the effect on the moments when
adding an increasing number of RL agents to the model. The intention is to
gain insight into the ability of the underlying e-ABM model to support an
increasing number of learning agents, both optimally and optimally buying,
before becoming saturated, and hence to potentially expose the domain in which
the interplay between the RL agents can come to dominate market dynamics.

## Reproducibility

The code used to implement, calibrate, and perform the sensitivity analysis
can be found on the following GitHub branch:
https://github.com/matthewdicks98/MDTG-MALABM/tree/Calibrated-ABM. The code
used to implement each RL agent inside the e-ABM can be found on the following
GitHub branch: https://github.com/matthewdicks98/MDTG-MALABM/tree/RL-ABM.

## Acknowledgements

We thank Andrew Paskaramoorthy for helpful comments and critique. We thank
Ivan Jericevich for support and advice with respect to the matching engine and
hybrid agent-based model implementation.

## References

  * Almgren and Chriss [2001] Almgren, R., Chriss, N., 2001\.  Optimal execution of portfolio transactions.  Journal of Risk 3, 5–40.  doi:10.21314/JOR.2001.041.
  * Aloud et al. [2017] Aloud, M., Fasli, M., Tsang, E., Dupuis, A., Olsen, R., 2017.  Modeling the high-frequency fx market: An agent-based approach.  Computational Intelligence 33, 771–825.  doi:10.1111/coin.12114.
  * Barto and Mahadevan [2003] Barto, A.G., Mahadevan, S., 2003\.  Recent advances in hierarchical reinforcement learning.  Discrete event dynamic systems 13, 41–77.  doi:10.1023/A:1022140919877.
  * Bellman [1954] Bellman, R., 1954.  The theory of dynamic programming.  Bulletin of the American Mathematical Society 60, 503–515.  doi:10.1090/S0002-9904-1954-09848-8.
  * Bertsekas and Tsitsiklis [1996] Bertsekas, D.P., Tsitsiklis, J.N., 1996\.  Neuro-dynamic programming.  Athena Scientific.  doi:10.1007/978-0-387-74759-0_440.
  * Bertsimas and Lo [1998] Bertsimas, D., Lo, A.W., 1998\.  Optimal control of execution costs.  Journal of Financial Markets 1, 1–50.  doi:10.1016/S1386-4181(97)00012-8.
  * Bezanson et al. [2021] Bezanson, J., Edelman, A., Karpinski, S., Shah, V.B., 2021\.  Asynchronous Programming.  URL: https://docs.julialang.org/en/v1/manual/asynchronous-programming/#. accessed: 2022-03-15.
  * Bouchaud et al. [2003] Bouchaud, J.P., Gefen, Y., Potters, M., Wyart, M., 2003\.  Fluctuations and response in financial markets: the subtle nature ofrandom’price changes.  Quantitative finance 4, 176\.  doi:10.1080/14697680400000022.
  * Bouchaud et al. [2002] Bouchaud, J.P., Mézard, M., Potters, M., 2002.  Statistical properties of stock order books: empirical results and models.  Quantitative finance 2, 251–256.  doi:10.1088/1469-7688/2/4/301.
  * Cartea et al. [2015] Cartea, Á., Jaimungal, S., Penalva, J., 2015.  Algorithmic and high-frequency trading.  Cambridge University Press.
  * Chang et al. [2021] Chang, P., Pienaar, E., Gebbie, T., 2021.  The Epps effect under alternative sampling schemes.  Physica A: Statistical Mechanics and its Applications , 126329URL: https://www.sciencedirect.com/science/article/pii/S0378437121006026, doi:https://doi.org/10.1016/j.physa.2021.126329.
  * Clauset et al. [2009] Clauset, A., Shalizi, C.R., Newman, M.E., 2009.  Power-law distributions in empirical data.  SIAM review 51, 661–703.  URL: http://www.jstor.org/stable/25662336, doi:10.1137/070710111.
  * Cont [2001] Cont, R., 2001.  Empirical properties of asset returns: stylized facts and statistical issues.  Quantitative Finance 1, 223–236.  doi:10.1080/713665670.
  * Crafa [2021] Crafa, S., 2021.  From agent-based modeling to actor-based reactive systems in the analysis of financial networks.  Journal of Economic Interaction and Coordination 16, 649–673.  URL: https://doi.org/10.1007/s11403-021-00323-8, doi:10.1007/s11403-021-00323-8.
  * Dieci and He [2018] Dieci, R., He, X.Z., 2018.  Heterogeneous agent models in finance.  Handbook of computational economics 4, 257–328.  URL: https://www.sciencedirect.com/science/article/pii/S157400211830008X, doi:10.1016/bs.hescom.2018.03.002.
  * Dijkstra [1968] Dijkstra, E.W., 1968.  Cooperating sequential processes, in: The origin of concurrent programming. Springer, pp. 65–138.  doi:10.1007/978-1-4757-3472-0_2.
  * Doya [1996] Doya, K., 1996.  Efficient nonlinear control with actor-tutor architecture, in: Mozer, M., Jordan, M., Petsche, T. (Eds.), Advances in Neural Information Processing Systems, MIT Press.  URL: https://proceedings.neurips.cc/paper/1996/file/eb86d510361fc23b59f18c1bc9802cc6-Paper.pdf.
  * Fabretti [2013] Fabretti, A., 2013.  On the problem of calibrating an agent based model for financial markets.  Journal of Economic Interaction and Coordination 8, 277–293.  doi:10.1007/s11403-012-0096-3.
  * Farmer et al. [2005] Farmer, J.D., Patelli, P., Zovko, I.I., 2005.  The predictive power of zero intelligence in financial markets.  Proceedings of the National Academy of Sciences 102, 2254–2259.  URL: https://www.pnas.org/content/102/6/2254, doi:10.1073/pnas.0409157102.
  * Gao and Han [2012] Gao, F., Han, L., 2012.  Implementing the nelder-mead simplex algorithm with adaptive parameters.  Computational Optimization and Applications 51, 259–277.  doi:10.1007/s10589-010-9329-3.
  * Garcia and Ndiaye [1998] Garcia, F., Ndiaye, S.M., 1998\.  A learning rate analysis of reinforcement learning algorithms in finite-horizon, in: ICML ’98: Proceedings of the Fifteenth International Conference on Machine Learning, Morgan Kaufmann Publishers Inc., San Francisco, CA, USA. pp. 215–223.
  * Gaunersdorfer and Hommes [2007] Gaunersdorfer, A., Hommes, C., 2007\.  A nonlinear structural model for volatility clustering, in: Long memory in economics. Springer, pp. 265–288.  doi:10.1007/978-3-540-34625-8_9.
  * Geweke and Porter-Hudak [1983] Geweke, J., Porter-Hudak, S., 1983\.  The estimation and application of long memory time series models.  Journal of time series analysis 4, 221–238.  doi:10.1111/j.1467-9892.1983.tb00371.x.
  * Gilles [2006] Gilles, D., 2006.  Asynchronous simulations of a limit order book.  Master’s thesis. University of Manchester.
  * Gilli and Winker [2003] Gilli, M., Winker, P., 2003\.  A global optimization heuristic for estimating agent based models.  Computational Statistics & Data Analysis 42, 299–312.  doi:10.1016/S0167-9473(02)00214-1.
  * Goosen and Gebbie [2020] Goosen, K., Gebbie, T., 2020\.  Calibrating high-frequency trading data to agent-based models using approximate bayesian computation.  Master’s thesis. University of cape town.  URL: http://hdl.handle.net/11427/33699, doi:10.25375/uct.12894005.v1.
  * Gould et al. [2013] Gould, M.D., Porter, M.A., Williams, S., McDonald, M., Fenn, D.J., Howison, S.D., 2013\.  Limit order books.  Quantitative Finance 13, 1709–1742.  doi:10.1080/14697688.2013.803148.
  * Harvey et al. [2016] Harvey, M., Hendricks, D., Gebbie, T., Wilcox, D., 2016\.  Deviations in expected price impact for small transaction volumes under fee restructuring.  Physica A: Statistical Mechanics and its Applications 471, 416–426.  doi:10.1016/j.physa.2016.11.042.
  * Hendricks [2017] Hendricks, D., 2017.  Using real-time cluster configurations of streaming asynchronous features as online state descriptors in financial markets.  Pattern Recognition Letters 97, 21–28.  doi:10.1016/j.patrec.2017.06.026.
  * Hendricks and Wilcox [2014] Hendricks, D., Wilcox, D., 2014\.  A reinforcement learning extension to the almgren-chriss framework for optimal trade execution, in: 2014 IEEE Conference on Computational Intelligence for Financial Engineering & Economics (CIFEr), IEEE. pp. 457–464.  doi:10.1109/CIFEr.2014.6924109.
  * Huberman and Stanzl [2004] Huberman, G., Stanzl, W., 2004\.  Price manipulation and quasi-arbitrage.  Econometrica 72, 1247–1275.  doi:https://doi.org/10.1111/j.1468-0262.2004.00531.x.
  * Jericevich et al. [2020] Jericevich, I., Chang, P., Gebbie, T., 2020.  Comparing the market microstructure between two South African exchanges.  arXiv:2011.04367.
  * Jericevich et al. [2021a] Jericevich, I., Chang, P., Gebbie, T., 2021a.  Simulation and estimation of a point-process market-model with a matching engine.  arXiv:2105.02211.
  * Jericevich et al. [2021b] Jericevich, I., Chang, P., Gebbie, T., 2021b.  Simulation and estimation of an agent-based market-model with a matching engine.  arXiv:2108.07806.
  * Jericevich et al. [2022] Jericevich, I., Sing, D., Gebbie, T., 2022.  Cointossx: An open-source low-latency high-throughput matching engine.  SoftwareX 19, 101136\.  URL: https://www.softxjournal.com/article/S2352-7110(22)00087-5/pdf, doi:10.1016/j.softx.2022.101136.
  * Johannesburg Stock Exchange [2013] Johannesburg Stock Exchange, 2013.  New Equity Market Trading and Information Solution: Volume 00 – Trading and Information Overview.  Technical Report. JSE.  V2.03.
  * Kyle [1985] Kyle, A.S., 1985.  Continuous auctions and insider trading.  Econometrica 6, 1315–1335.  doi:doi.org/10.2307/1913210.
  * LeBaron [2000] LeBaron, B., 2000.  Agent-based computational finance: Suggested readings and early research.  Journal of Economic Dynamics and Control 24, 679–702.  URL: https://www.sciencedirect.com/science/article/pii/S0165188999000226, doi:10.1016/S0165-1889(99)00022-6.
  * Lee and Ready [1991] Lee, C.M.C., Ready, M.J., 1991\.  Inferring trade direction from intraday data.  The Journal of Finance 46, 733–746.  doi:10.1111/j.1540-6261.1991.tb02683.x.
  * Lillo and Farmer [2007] Lillo, F., Farmer, J.D., 2007\.  The long memory of the efficient market.  Studies in Nonlinear Dynamics & Econometrics 8, 1–1.  doi:10.2202/1558-3708.1226.
  * Lillo et al. [2003] Lillo, F., Farmer, J.D., Mantegna, R., 2003.  Master curve for price-impact function.  Nature 421, 129–30.  doi:10.1038/421129a.
  * Mandelbrot [1963] Mandelbrot, B.B., 1963.  The variation of certain speculative prices.  The Journal of Business 36, 394–419.  URL: http://www.jstor.org/stable/2350970.
  * Mandeş [2015] Mandeş, A., 2015.  Microstructure-based order placement in a continuous double auction agent based model.  Algorithmic Finance 4, 105–125.  doi:10.3233/AF-150049.
  * Mandeş and Winker [2017] Mandeş, A., Winker, P., 2017\.  Complexity and model comparison in agent based modeling of financial markets.  Journal of Economic Interaction and Coordination 12, 469–506.  doi:10.1007/s11403-016-0173-0.
  * Microsoft Azure [2022] Microsoft Azure, 2022.  Linux virtual machines.  Overview of Linux VMs in Azure (PDF Documentation 3/10/2021).  URL: https://docs.microsoft.com/en-us/azure/virtual-machines/linux/overview. accessed: 2022-05-24.
  * Nelder and Mead [1965] Nelder, J.A., Mead, R., 1965\.  A simplex method for function minimization.  The Computer Journal 7, 308–313.  doi:https://doi.org/10.1093/comjnl/7.4.308.
  * Nevmyvaka et al. [2006] Nevmyvaka, Y., Feng, Y., Kearns, M., 2006.  Reinforcement learning for optimized trade execution, in: Proceedings of the 23rd international conference on Machine learning, pp. 673–680.  doi:10.1145/1143844.1143929.
  * Ning et al. [2018] Ning, B., Lin, F.H.T., Jaimungal, S., 2018.  Double deep q-learning for optimal execution.  arXiv:1812.06600.
  * Nuyts [2010] Nuyts, J., 2010.  Inference about the tail of a distribution: Improvement on the hill estimator.  International Journal of mathematics and mathematical sciences 2010.  doi:10.1155/2010/924013.
  * Obizhaeva and Wang [2013] Obizhaeva, A.A., Wang, J., 2013\.  Optimal trading strategy and supply/demand dynamics.  Journal of Financial Markets 16, 1–32.  doi:10.1016/j.finmar.2012.09.001.
  * Perold [1988] Perold, A.F., 1988.  The implementation shortfall: Paper versus reality.  Journal of Portfolio Management 14, 4.  doi:10.3905/jpm.1988.409150.
  * Platt [2020] Platt, D., 2020.  A comparison of economic agent-based model calibration methods.  Journal of Economic Dynamics and Control 113, 103859.  URL: https://www.sciencedirect.com/science/article/pii/S0165188920300294, doi:10.1016/j.jedc.2020.103859.
  * Platt and Gebbie [2016] Platt, D., Gebbie, T., 2016\.  The problem of calibrating an agent-based model of high-frequency trading.  arXiv:1606.01495.
  * Platt and Gebbie [2018] Platt, D., Gebbie, T., 2018\.  Can agent-based models probe market microstructure?  Physica A: Statistical Mechanics and its Applications 503, 1092–1106.  URL: https://www.sciencedirect.com/science/article/pii/S0378437118309956, doi:10.1016/j.physa.2018.08.055.
  * Potters and Bouchaud [2003] Potters, M., Bouchaud, J.P., 2003\.  More statistical properties of order books and price impact.  Physica A: Statistical Mechanics and its Applications 324, 133–140.  doi:10.1016/S0378-4371(02)01896-4.
  * Sing [2017] Sing, D., 2017.  Cointossx software.  URL: https://github.com/dharmeshsing/CoinTossX.
  * Sing and Gebbie [2017] Sing, D., Gebbie, T., 2017.  JSE Matching Engine Simulator.  Master’s thesis. University of the Witwatersrand.  URL: http://wiredspace.wits.ac.za/handle/10539/25136.
  * Sutton and Barto [2018] Sutton, R.S., Barto, A.G., 2018\.  Reinforcement learning: An introduction.  MIT press.
  * Toth et al. [2015] Toth, B., Palit, I., Lillo, F., Farmer, J.D., 2015\.  Why is equity order flow so persistent?  Journal of Economic Dynamics and Control 51, 218–239.  doi:10.1016/j.jedc.2014.10.007.
  * Von Neumann and Morgenstern [1944,2007] Von Neumann, J., Morgenstern, O., 1944,2007.  Theory of games and economic behavior.  Princeton Classic Editions, Princeton university press.  doi:10.2307/2280142.
  * Watkins and Dayan [1992] Watkins, C.J., Dayan, P., 1992\.  Q-learning.  Machine learning 8, 279–292.  doi:10.1007/BF00992698.
  * Watkins [1989] Watkins, C.J.C.H., 1989.  Learning from delayed rewards.  Ph.D. thesis. King’s College, Cambridge United Kingdom.
  * [63] Wilcox, D., Gebbie, T., .  Hierarchical causality in financial economics.  Available at SSRN: URL: https://ssrn.com/abstract=2544327, doi:http://dx.doi.org/10.2139/ssrn.2544327.
  * Winker et al. [2007] Winker, P., Gilli, M., Jeleskovic, V., 2007.  An objective function for simulation based inference on exchange rate data.  Journal of Economic Interaction and Coordination 2, 125–145.  doi:10.1007/s11403-007-0020-4.

##  Appendix A e-ABM parameters

Table 5: Describes the parameters of the ABM, and detail whether the parameters were kept fixed or were allowed to vary in the calibration. If the variable was kept fixed its value is also shown. For a detailed description of why certain parameters were kept fixed and others were allowed to vary see F.3. |  |  |
---|---|---|---
Parameters |  Description | Type | Value
$T$ |  Total simulation time, not including the initialisation of the LOB. | Fixed | 25s
$\phi$ |  Time to cancellation of an order. | Fixed | 1s
$N_{{}_{\mathrm{LP}}}$ |  The number of liquidity providers (HFTs). | Fixed | 30
$N^{c}_{{}_{\mathrm{LT}}}$ |  The number of chartist liquidity takers.. | Free | -
$N^{f}_{{}_{\mathrm{LT}}}$ |  The number of fundamentalist liquidity takers. | Free | -
$\delta$ |  The mid-price cut-off $\%$ to set minimum volumes. | Free | -
$\kappa$ |  Limit order placement intensity for HFTs. | Free | -
$\nu$ |  Shape parameter setting the order-imbalance (equation 1). | Free | -
$\sigma_{f}$ |  Fundamental value variance | Free | -
$\lambda_{min}$ |  The minimum chartist forgetting factor | Fixed | 0.0005
$\lambda_{max}$ |  Maximum chartist forgetting factor | Fixed | 0.05

##  Appendix B Reinforcement learning

A Reinforcement Learning (RL) agent interacts with the environment at discrete
time steps, $t=1,2,3,...$ 181818We restrict ourselves to considering only the
discrete-time case as this is what our agent will observe. However, we note
that continuous-time extensions to this formulation are available (see [5,
17]).. At each time step the agent observes the state of the environment
$S_{t}\in\mathcal{S}$ where $\mathcal{S}$ is the set of possible states, and
takes an action $A_{t}\in\mathcal{A}(S_{t})$ where $\mathcal{A}(S_{t})$ is the
set of actions the agent can take in state $S_{t}$. At the following time step
the agent observes the reward $R_{t+1}\in\mathcal{R}\subset\mathbb{R}$
191919We use $R_{t+1}$ instead of $R_{t}$ to denote the reward received from
taking action $A_{t}$ in state $S_{t}$ to indicate that $R_{t+1}$ and
$S_{t+1}$ are observed at the same time. This convention was taken from [58].,
for taking action $A_{t}$ in state $S_{t}$, and it observes the new state
$S_{t+1}$.

To select the actions at each time step the agent uses a policy, which maps
the state to the probabilities of selecting each possible action [58]. The
policy is denoted by $\pi_{t}$ where $\pi_{t}(a|s)$ is the probability that
$A_{t}=a$ when $S_{t}=s$ [58].

The goal of the agent is to find the policy such that it maximises the total
reward received. Or more formally, at each time step $t$ the agent wishes to
maximise the expected return, $G_{t}$, which is a function of the reward
sequence. In the simplest case this is just the sum of the rewards:
$G_{t}=R_{t+1}+R_{t+2}+R_{t+3}+\cdots+R_{T}$, where T is the final time step
[58]. This is the expected return function that our agent will be using.

This formulation of the return maps nicely to problems that have a final time
step, where the agent-environment interaction can be broken up into sub-
sequences, such as repeated plays of a game [58]. For example, blackjack, or
in our case repeated executions of a large parent order. These repeated tasks
are called episodes and tasks that can be broken up into episodes are called
episodic tasks. At the final time step, of an episode, the agent moves into a
terminal state, after which the state of the environment gets reset to some
starting state.

Some tasks can’t be broken up into episodes, where $T=\infty$. A trading agent
that just aims to maximise its profit over its lifetime would be an example of
such an agent. These tasks are called continuing tasks [58]. Since $T=\infty$,
we have that the expected return could be infinite. This is remedied by adding
discounting, which gives us the expected discounted return:
$G_{t}=\sum_{k=0}^{\infty}\gamma^{k}R_{t+k+1}$. Here $0\leq\gamma\leq 1$ is
the discounting rate, and it controls how much the agent values immediate vs.
farsighted returns. As $\gamma$ approaches zero the agent becomes more myopic
and values immediate rewards more heavily, while as $\gamma$ approaches one
the agent becomes more farsighted and starts to value future returns more
heavily. Now, an issue occurs because we have two tasks and two reward
functions. To recover a single notation that describes both tasks we use the
convention given to us by [58], where we consider episode termination to be
the case of entering an absorbing state that transitions only to itself and
has a reward of zero. This allows us to describe both tasks return functions.

Given action $a$ and state $s$, the probability of observing a reward $r$ and
moving to state $s^{\prime}$ is given by:

$\displaystyle
p(s^{\prime},r|s,a)=\mathbb{P}(S_{t+1}=s^{\prime},R_{t+1}=r|S_{t}=s,A_{t}=a).$
(12)

This defines the one-step dynamics of the model for the MDP. From this we can
find the state-transition probabilities:

$\displaystyle p(s^{\prime}|s,a)=\sum_{r\in\mathcal{R}}p(s^{\prime},r|s,a),$
(13)

as well as the expected reward for taking action $a$ in state $s$:

$\displaystyle
r(s,a)=\sum_{r\in\mathcal{R}}r\sum_{s^{\prime}\in\mathcal{S}}p(s^{\prime},r|s,a).$
(14)

Using these definitions we can derive equations for the state-value function
$v_{\pi}(s)$ and the action-value function $q_{\pi}(s,a)$ for policy $\pi$,
where $v_{\pi}(s)$ is the value of state $s$ under policy $\pi$, which is
defined as:

$\displaystyle
v_{\pi}(s)=\mathbb{E}_{\pi|S_{t}=s}[G_{t}]=\mathbb{E}_{\pi|S_{t}=s}\left[\sum_{k=0}^{\infty}\gamma^{k}R_{t+k+1}\right]=\sum_{a}\pi(a|s)\sum_{s^{\prime},r}p(s^{\prime},r|s,a)\left[r+\gamma
v_{\pi}(s^{\prime})\right].$ (15)

Here the action-value function, which is the value of taking action $a$ in
state $s$, and following $\pi$ thereafter, is:

$\displaystyle q_{\pi}(s,a)$
$\displaystyle=\mathbb{E}_{\pi|\begin{subarray}{c}S_{t}=s\\\
A_{t}=a\end{subarray}}[G_{t}]=\mathbb{E}_{\pi|\begin{subarray}{c}S_{t}=s\\\
A_{t}=a\end{subarray}}\left[\sum_{k=0}^{\infty}\gamma^{k}R_{t+k+1}\right]=\sum_{r,s^{\prime}}p(s^{\prime},r|s,a)\left[r+\gamma\sum_{a^{\prime}}\pi(a^{\prime}|s^{\prime})q_{\pi}(s^{\prime},a^{\prime})\right].$
(16)

Value functions give us a partial ordering of policies. We say that policy
$\pi$ is better than or equal to a policy $\pi^{\prime}$ if its expected
return is greater than or equal to the expected return of $\pi^{\prime}$ for
all states [58]. The optimal policy is the policy such that no other policy
has a greater expected return for any of the states. There could be more than
one of these policies, however, we denote the optimal policy as $\pi_{*}$.
Using this ordering of policies we can define the optimal state-value function
as: $v_{*}(s)=\max_{\pi}v_{\pi}(s)\leavevmode\nobreak\ \forall
s\in\mathcal{S}$, and the optimal action-value function as:
$q_{*}(s,a)=\max_{\pi}q_{\pi}(s,a)$ $\forall s\in\mathcal{S}$ and $\forall
a\in\mathcal{A}$. The optimal policy can then be found by maximising over
$q_{*}(s,a)$ using that if $a=\mathrm{argmax}_{a\in\mathcal{A}}q_{*}(s,a)$
then $\pi^{*}(s,a)=1$, otherwise $\pi^{*}(s,a)=0$. Therefore if we know
$q_{*}(s,a)$ we know the optimal policy. This is key since most RL algorithms,
and the Q-learning algorithm that we will implement will aim to optimise the
action-value function directly. Therefore, we will need to recover the optimal
policy from the optimal action-value function.

We do not have the full model of the environment therefore, we can’t use
equations 15 and 16. However, we can learn the optimal state-value and action-
value functions from samples of sequences of states, actions, and rewards,
that can be generated through actual or simulated interaction with the
environment. The two main methods for learning these functions through
experience are Monte Carlo and Temporal Difference (TD) learning. When using
Monte Carlo methods to estimate the value functions, the value functions are
only updated at the end of every episode. This can make learning inefficient.
Moreover, Monte Carlo methods can’t be used to learn in a step-by-step online
way. Given the learning inefficiencies and the long runtimes of a simulation,
it is likely that we won’t be able to have enough episodes for the RL agent to
accurately estimate the value functions, using Monte Carlo methods. Therefore,
we will only consider TD learning.

TD methods allow the agent to learn in each episode and does not require
waiting for the episode to complete like Monte Carlo. The TD algorithm that we
will be using is the Q-learning algorithm first proposed by [62]. In this
particular case, the update to the value function can be done after a single
step in an episode.

Q-learning is an off-policy TD control algorithm which is in contrast to the
SARSA on-policy control algorithm. This means that the policy that generates
the actions is not the same as the target policy. Using this method simplifies
the analysis of the algorithm and reduces the conditions needed for guaranteed
convergence to $q_{*}$.

At a particular time t in the an arbitrary episode we give the simplest form
of Q-learning, the one-step Q-learning update equation is as follows:

$\displaystyle Q(S_{t},A_{t})=$ $\displaystyle Q(S_{t},A_{t})+\alpha
R_{t+1}+\alpha\left[{\gamma\max_{a}Q(S_{t+1},a)-Q(S_{t},A_{t})}\right].$ (17)

Here $\alpha$ is the learning rate. Provided each state is visited infinitely
often and each action is taken in each state infinitely often we have that $Q$
converges to $q_{*}$ with probability 1, for appropriately decaying step
sizes202020The conditions on the step sizes are $\sum_{t}\alpha_{t}=\infty$
and $\sum_{t}\alpha_{t}^{2}<\infty$. [61].

##  Appendix C System implementation

Table 6: Describes the hardware specifications on which the simulations were run and shows indicative runtimes for a single simulation of the ABM, the calibrations and sensitivity analysis, as well as the training time for a single RL agent and the total amount of time need to train all the RL agents. The sensitivity analysis and calibration were run on Virtual Machine 1 and the RL experiments were run on the upgraded Virtual Machine 2. The approximate runtimes for a single run of the simulation as well as the sensitivity analysis, calibration, training time for a single RL agent and the time it took to train all RL agents. Hardware | Specifications | Primary Task | Runtime
---|---|---|---
Virtual Machine 1 | Standard A4m V2 | Sensitivity Analysis | 30 hours
Operating system | 64-bit Ubuntu 20.04.4 LTS | Calibration | 8.3 hours
CPU | Intel(R) Xeon(R) Platinum 8171M |  |
| @ 2.60GHz with 4 vCPUs |  |
RAM | 32GiB |  |
Virtual Machine 2 | Standard F16s V2 | 1 simulation run | 26 seconds212121The actual simulation time is 25 seconds but due to the Just-In-Time (JIT) compilation of Julia and the initialisation of the LOB a single simulation run takes approximately 26 seconds.
Operating system | 64-bit Ubuntu 20.04.4 LTS | Training time for 1 RL agent | 8.3 hours
CPU | Intel(R) Xeon(R) Platinum 8272CL | Total RL training time | 50 hours
| @ 2.60GHz with 16 vCPUs |  |
RAM | 32GiB |  |

##  Appendix D Stylised facts

Table 7: Defines and describes the stylised facts that we aim to reproduce with our ABM. Stylised Fact |  Description
---|---
Fat/Heavy tails |  The unconditional distribution of the log-returns has been shown to have power-law behaviour in the tails [42], with a finite tail-index between 2 and 5 [13].
Return autocorrelation |  The log-returns exhibit an absence of significant autocorrelations [13]. In high-frequency time scales the return series observes negative significant autocorrelation for the first few lags which are usually given as evidence for a fast mean-reverting component[13].
Trade-sign autocorrelation |  The trade sign time-series has been shown to have persistence in their autocorrelations [40, 8, 59]; implying that buy (resp. sell) orders are followed by more buy (resp. sell) orders. Herding and order-splitting have been used to explain this observed phenomena; with order-splitting being the main reason for the persistence in order flow [59].
Volatility Clustering |  Here large (resp. small) price variations are often followed by large (resp. small) price variations. The persistence of the autocorrelation of the absolute log-returns are signatures of volatility clustering and long-range dependence in the volatility [13].
Depth profile |  The Paris Bourse and the NASDAQ have an observed hump in their shape with a slight increase in volume as the price moves away from the best bid/offer, with it reaching a maximum a few ticks away from the best bid/offer before a steep decline [9, 27]222222The depth profile measures the average volume observed at different prices on the ask and bid sides.. The Standard and Poor’s Depository Receipts (SPY) showed a maximum at the best bid/offer [55].
Price impact |  Price impact can be measured as the change in bid/ask price caused by a trade of a given volume. As the volume increases the observed price changes increase linearly on a log-log scale indicating a power-law relationship [28, 32, 41].

Figure 20: Plots the log-return distributions for the ABM (green) and the JSE
(purple). Using QQ-plots the insets compare the distributions to the normal
distribution fitted to the data using MLE and reveal that the log-return
distributions have fatter tails than the normal distribution, demonstrating
the ABM’s ability to reproduce the required stylised fact.

Figure 21: The extreme distributions taken from the upper 95th percentile as
well as the lower 5th percentile for the ABM (green) and JSE (purple) are
plotted. This figure demonstrates that both return distributions have near
power-law behaviour in their tails. If their distributions were exactly a
power-law we would expect that the data would be linear on a log-log scale.
However, we observe a slight downward curve which indicates only near power-
law behaviour. Comparing the empirical-versus-simulated returns we see that
both distributions observe this behaviour, giving credibility to the returns
generated by the ABM. We also note that the ABM’s distribution is fatter
tailed when compared with the JSE’s distribution.

Figure 22: Plots the autocorrelations in the log-returns and the inset plots
the autocorrelations in the absolute log-returns, for the ABM (green) and the
JSE (purple). The significant negative autocorrelation in the log-returns for
the first few lags is present in both the plots. The absolute log-returns for
both the simulated and empirical both have persistence in their
autocorrelations, meaning that both time series exhibit evidence for
volatility clustering. However, the decay in the autocorrelations is far
quicker for the simulated returns.

Figure 23: The autocorrelations in the trade-signs are plotted for the ABM
(green) and the JSE (purple). The model’s autocorrelation function (ACF) does
not exhibit the same behaviour as the JSE. The ABM does capture the stylised
fact of persistent decay of the trade-signs. However, the decay is far slower
than that found in the JSE data. The insets plot the trade-sign decay on a
log-linear scale where the difference in the rate of decay is also clearly
demonstrated.

(a) The power-law decay in the order-flow autocorrelation show the long-memory
property Lillo and Farmer [40]. Here we plot the trade-sign autocorrelations
for the JSE on a log-log scale and observe an approximately linear decrease to
recover similar results. Using the MLE of the power-law described in equation
9, we show that $\alpha=1.286$, which means that the distribution is not
integrable and exhibits the long-memory property.

(b) The depth profile is computed as the average volume observed at seven
price levels throughout the simulation [34]. The maximum volume occurring at
the best bid/ask and the rapid decrease in volume as the price increases is
consistent with empirical data [55] as well as previous ABM results [34].

##  Appendix E Sensitivity Analysis

Table 8: Defines and describes the moments that are used to characterise the price paths. These moments are used in the sensitivity analysis to see how changing the parameters changes the price paths and they are used to calibrate the ABM to market data found on the JSE. Moments |  Description
---|---
$\mu$, $\sigma$ |  The mean and standard deviation measure the shape of the log-return distribution. The mean is expected to be close to zero.
KS |  The Kolmogorov-Smirnov statistic is used to compare the maximum absolute difference between the CDFs of the empirical and simulated log-returns. This is used to test if the samples are drawn from the same distribution.
GPH |  Provides an estimate for the fractional integration parameter $d\in[-0.5,0.5]$ in absolute log-returns time series. [23] estimates this by estimating the memory parameter in the log-periodogram. This provides a measure for the long-range dependence in the log-returns.
ADF |  The augmented Dicky-Fuller statistic tests for the presence of a unit root in the time series of the log-returns. The more negative the test statistic the less evidence there is for the presence of a unit root.
Garch(1,1): $\alpha+\beta$ |  The parameters of the GARCH(1,1) model provide a useful summary for the short-range dependence in the squared returns. [64] found that the sum of the parameters in the GARCH(1,1) model is a more robust statistic when compared with the individual parameters.
”Improved” Hill Estimator (HE) |  The Hill estimator is used to estimate the tail-index of a power-law distribution which can be used to infer the power-law behaviour in the tails of experimental distributions [49]. [49] proposed an improved Hill estimator to remedy the issues with the standard estimator.
Hurst Exponent (H) |  Provides a measure of the long-memory of the log-returns. If $H\in[0,1]$ is significantly far from 0.5 we have evidence that the log-returns do not follow a random walk. Furthermore, if $H\in[0,0.5)$ then the process is said to be mean-reverting and if $H\in(0.5,1]$ the process is said to be trending.

Due to computational constraints, six parameters were chosen to be varied in
the sensitivity analysis. These are the free parameters defined in Table 5 and
were informed by prior work [34]. The parameters that were held constant
during the sensitivity analysis as well as their values are also shown. The
$\lambda_{min}$ and $\lambda_{max}$ parameters were chosen in such a way that
a chartist agent has an equal probability of having a long or short trading
horizon, and the time to cancellation, of the LPs limit-orders, was
arbitrarily tuned to ensure that liquidity is not being removed too quickly
from the order book, which would cause frequent liquidity crashes, and having
too slow a cancellation time, which would cause almost no cancellations. This
would cause behaviour that does not represent that of a high-frequency agent,
which the LPs are trying to model. The total simulation time, after
initialisation, $T$ was set at 25 seconds because it was the time that
produced order counts most similar to the day of trading on the JSE that we
will use to calibrate the model.

Comparing the free parameters in my model with that of [34], it can be seen
that I am varying the number of chartists and the number of fundamentalists
while keeping the number of high-frequency agents constant. However, [34] kept
the number of chartists and fundamentalists constant, at one agent each, and
varied the number of high-frequency agents to determine how the ratio of
liquidity providers to liquidity takers,
$N=\frac{N_{{}_{\mathrm{LP}}}}{N_{{}_{\mathrm{LT}}}}$ , changed the simulated
price paths. During initial testing of the model, it was observed that having
different values for $N^{f}_{{}_{\mathrm{LT}}}$ and $N^{c}_{{}_{\mathrm{LT}}}$
changed the behaviour of the simulated price paths. For example, when there
are more chartists ($N^{c}_{{}_{\mathrm{LT}}}=12$) than fundamentalists
($N^{f}_{{}_{\mathrm{LT}}}=5$) the log-return distribution has a smaller
variance and but has a lower tail-index when compared with more
fundamentalists and less chartists. A possible reason for this is there are
fewer large orders placed by fundamentalists and smaller orders placed by the
chartists on account of equations 4 and 5. The initial lags in the trade-sign
autocorrelation are larger when there are more chartists and fundamentalists,
which is possibly attributed to greater herding behaviour. Furthermore, as the
number of fundamentalists increases more large orders are placed, and this
caused wider spreads, which would give rise to more extreme price moves.

Due to these observed changes, it was decided that a better parameter set
would be to keep the number of high-frequency agents constant and vary both
the number of chartists and fundamentalists. Since we are still changing the
ratio of liquidity takers to liquidity providers we will also be changing N
and therefore, we will still be modelling how a changing N affects the moments
as was done in [34].

When implementing the sensitivity analysis we considered 4 different values
for each parameter, and for each value, we varied all other parameters and
estimated the moments of all the simulated price paths, which leads to a total
of 4096 parameter combinations. Therefore, for each parameter value and each
moment, there were 1024 values to be aggregated over. Another benefit of this
method is that it allowed us to create parameter surfaces for the 2 way
interactions for each of the parameters. These surfaces allow us to visualise
how the interactions of two parameters affect the statistical properties of
the simulated price paths and the objective function. These effects can be
seen in Figures 26 and 28.

It must be noted that when implementing the sensitivity analysis, for each of
the parameter combinations the moments were estimated for a single run of the
simulation with a single seed set which ensured the same fundamental prices
and the same chartist forgetting factors. This was done due to time and
computational cost constraints. Therefore, these results should be considered
indicative rather than presenting the exact surfaces and moment values as
there will be uncertainty around these values that have not been reduced
through Monte Carlo replications over multiple seeds.

To ensure that the simulations for each of the parameter combinations remain
independent of each other we had to start and stop trading sessions with
CoinTossX before and after each simulation was run. This means, the LOB was
cleared before the start of the next session and there were no orders from the
previous session that could be matched with orders from the current session.

Due to a large number of parameter combinations after each simulation, the
Julia garbage collector was run. Helping to free up memory from the Julia
program and free up memory in the JVM created by the Julia program. This was
recommended in the JavaCall documentation. They note that the pressure on the
Julia and java heaps are different and the Julia garbage collector may not be
run even if the java heap is full and memory could be freed with garbage
collection. Therefore, they recommend periodically manually running the Julia
garbage collector to ensure this issue does not cause an error.

###  E.1 Marginals

Figure 25, shows the distributions of each of the moments for different values
of each of the parameters as well as the objective function. These
distributions consist of holding one parameter value fixed and varying all
other parameters giving us 1024 samples for each of the distributions. For
almost all of the parameters, the distributions of the moments across their
values do not change considerably, and there is a large overlap of the
distributions. Indicating that a change to a single parameter value is not
going to cause significant changes to the model. One parameter-moment
combination where we do see a change in distribution when changing the
parameter values is between the number of fundamentalist traders and the Hurst
exponent. As we increase the number of fundamentalists the distribution of the
Hurst exponent shifts to smaller values where most of the distribution lies
further below 0.5. As described in Table 8, this represents more mean-
reversion in the time series. This is consistent with the fact that the
fundamentalists follow a mean-reverting trading strategy.

Due to large outliers and large variations in the objective function, for each
of the parameters, it is difficult to see a significant difference in the
distributions of the objective function across the range of parameter values.
The only observation we can make is that certain parameter values may be
subject to producing large outliers. For example, smaller values of the
$\kappa$ parameter produce larger outliers.

Figure 25: Plots the distributions of each of the moments and the objective
function for different values of each of the parameters. The columns and rows
represent the values of the moments and parameters respectively. These
distributions consist of holding one parameter value fixed and varying all
other parameters giving us 1024 samples for each of the distributions.

###  E.2 Objective surfaces

A better way to visualise the changes in the objective function as the
parameters change is to look at the behaviour of the objective function for an
interaction of two variables. Figure 26 plots the micro-price objective
function when holding two variables fixed and averaging over all other
combinations of the other parameters. Recall that the sensitivity analysis was
run once for a single seed so there were no Monte Carlo replications. The
surfaces are also course grained and made up of a 4x4 grid. Therefore, these
surfaces are just an indicative representation of the objective function and
are not to be interpreted as the true objective function. Despite this these
surfaces let us see some indicative patterns.

In almost all the plots we have that the objective function decreases from
areas of extremely high values too much lower values, giving evidence that
optimisation is possible. However, optimisation of the objective function will
not be trivial since there is evidence for local minima where the optimisation
can get stuck. For example, consider the surface for parameters
$N^{f}_{{}_{\mathrm{LT}}}$ and $\nu$. The optimisation algorithm could
possible get stuck in the local minima at $N^{f}_{{}_{\mathrm{LT}}}=9$ and
$\nu=7$ rather than reaching the actual minimum at
$N^{f}_{{}_{\mathrm{LT}}}=3$ and $\nu=2$. Another reason that the objective
function will be difficult to minimise is that the minima that the optimiser
converges to may depend heavily on the initial state of the optimiser. For an
extreme example, consider the interaction surface between
$N^{f}_{{}_{\mathrm{LT}}}$ and $\delta$. Depending on which side of the saddle
you are on will cause you to find different minima. These two issues make it
extremely hard to calibrate the ABM. Calibration difficulties are not new and
for detailed discussions on the difficulty of calibrating ABMs, we refer the
reader to [53, 52].

Figure 26: Plots the micro-price objective surfaces for all the two-way
interactions of the parameters. For each parameter combination, two of the
parameters are held fixed and we average over the micro-price objective
function for all other parameters. The figure demonstrates that there are
regions in the objective function that have higher and lower values and
therefore parameters can be differentiated by their objective function value
and the function can be optimised. It also shows that there are local minima
in the objective surface that will have to be taken into account when
performing optimisation.

###  E.3 Correlation matrix

The correlation matrix shown in Figure 27, depicts the correlations between
the parameters, the moments, and the objective function. There are only two
large correlations namely standard deviation with the objective function
(0.918) and the Hurst exponent with the number of fundamentalists (-0.699).
The first gives us that the as standard deviation in the log-returns increases
we have that there is an increase in the objective function, indicating that
the log-returns of the JSE data have a lower variance when compared with the
ABM. The second shows that as we increase the number of fundamentalists we
have that the Hurst exponent decreases. This effect can also be seen by
looking at the marginal distributions. Figure 27 also shows a mild negative
correlation (-0.505) between the number of value traders and the KS statistic,
indicating that as the number of value traders increases the distributions of
the log-returns become more similar.

Considering the bottom row of the correlation matrix we see that none of the
parameters are correlated with the objective function however, we have just
seen with the objective surfaces that the value of the objective function does
depend on the values of the parameters. Therefore, it is good to note that
while the correlations are good at finding the linear dependencies it does not
tell us about the non-linear dependencies in the model that are present. With
this in mind interaction plots of the parameters showing how they affect each
of the moments have been generated.

Figure 27: The correlation matrix shown depicts the correlations between the
parameters, the moments, and the objective function. There are only two large
correlations namely standard deviation with the objective function (0.918) and
the Hurst exponent with the number of fundamentalists (-0.699). A mild
negative correlation (-0.505) between the number of value traders and the KS
statistic is also present.

###  E.4 Moment surfaces

To show how the interaction of parameters affects the moments of the price
paths, moment surfaces were produced in the same way as was done for the
objective function. If the mean and standard deviation are not considered we
have six moments and six parameters, which means we will have 90 surfaces.
This is far too many surfaces to present here so we have considered just the
surfaces for each of the Hurst exponent, the Hill estimator, and the KS
statistic, for parameter combinations of $N^{c}_{{}_{\mathrm{LT}}}$ and
$N^{f}_{{}_{\mathrm{LT}}}$ as well as $\kappa$ and $\nu$. These can be found
in Figure 28. Considering the Hurst exponent surface for the
$N^{c}_{{}_{\mathrm{LT}}}$ and $N^{f}_{{}_{\mathrm{LT}}}$ parameters we can
see that the Hurst exponent decrease drastically as the number of
fundamentalists increase. There is also a small decrease in the Hurst exponent
as the number of chartists increases but this is a much smaller effect and
should not be interpreted as a trend. From this plot, there is evidence that
the mean reversion of the fundamentalists is dominating the momentum from
chartists. To illustrate some of the non-linear effects not captured by the
correlations, consider the Hill estimator surface for the
$N^{c}_{{}_{\mathrm{LT}}}$ and $N^{f}_{{}_{\mathrm{LT}}}$ parameters as well
as the $\kappa$ and $\nu$ parameters.

Figure 28: Micro-price moment surfaces show how the estimated moments change
on average when varying two of the parameters. Due to a large number of moment
surfaces we have just considered the surfaces for each of the Hurst exponent,
the Hill estimator, and the KS statistic, for parameter combinations
$N^{c}_{{}_{\mathrm{LT}}}$ and $N^{f}_{{}_{\mathrm{LT}}}$ as well as $\kappa$
and $\nu$. Looking at the Hurst exponent surface for the
$N^{c}_{{}_{\mathrm{LT}}}$ and $N^{f}_{{}_{\mathrm{LT}}}$ parameters we can
see that the Hurst exponent decreases drastically as the number of
fundamentalists increases, which is consistent with the fact that the
fundamentalists employ a mean reverting strategy.

##  Appendix F Calibration

Initial ABM research manually selected parameters, for the ABM, such that
these parameters were able to produce the required log-return stylised facts
[18]. There were no checks to see if the distribution of the simulated data
came from the same distribution as that of the empirical data [18]. This
presents a calibration problem in which the model must be calibrated to the
empirical data [53].

The three most common methods used to calibrate financial ABMs are (1) maximum
likelihood estimation, (2) Bayesian inference and, (3) method-of-moment using
simulated minimum distance (MM-SMD). Maximum likelihood is reserved for those
models involving closed-form solutions usually for daily sampled data [53],
approximating closing auctions and not intraday trading. Therefore, this
method was deemed unsuitable for calibrating our model. The second method,
Bayesian inference, constructs a posterior distribution of the parameters
using either the kernel density estimation (KDE) or parametric approaches.
Parametric methods have been said to often not be flexible enough to
accurately characterise the distribution [53]. However, Bayesian estimation
has been shown to outperform MM-SMD in a wide range of circumstances using the
KDE method [52], but this method is extremely computationally intensive for
large ABMs such as the one described in Section 2, and therefore was ruled out
as a calibration method. The method chosen is MM-SMD, because it is more
computationally tractable than Bayesian inference, it does not require a
closed-form solution, and it means that we are consistent with the methods
used in [34], whose model we try to replicate in event-time.

MM-SMD does have some pitfalls that need to be noted. The MM-SMD method has
been criticised for producing parameter degeneracies and causing haphazard
effects on the objective function [53]. These haphazard effects limit the
ability of the optimiser to find identifiable parameters due to the optimiser
converging to local optima. From the indicative objective surfaces in the
sensitivity analysis, we see that there are likely local optima in our
objective function, which could present the same calibration problems found in
[53]. The method has also been criticised for the arbitrary selection of
moments that will only represent a limited number of aggregate properties of
the data, and therefore may not capture important dynamics [53]. Despite these
shortcomings, MM-SMD presents a fast, simple and reliable (when there are more
moments than parameters) method for calibrating an ABM. However, due to these
shortcomings, these parameters should be considered indicative rather than
robust [34].

###  F.1 Calibration data

Since we are simulating a continuous-double-auction (CDA) all data relating to
the opening and closing auctions needed to be discarded. Therefore, the only
activity that occurred between 9:00 and 16:50 was considered. Moreover, it was
noticed that the first minute of trading produced erroneous data and it was
decided to remove the first minute of trading since it would make the model
more difficult to calibrate. It is also vital to remove any events related to
intraday volatility auctions and the impact of various futures close-outs.
Given that the future close-outs occur only in March, June, September and
December we do not need to consider this case. To remove most of the unwanted
trades such as after-hour trades (LT), correction of previous day’s trades
(LC) and auction uncrossing price trades (IP) I only considered automated
trades (AT).

Often there are larger trades that get executed against a set of smaller
trades. This means that single trade, at a given timestamp, gets split up into
smaller trades that change the best bid and best ask as it moves up and down
the order book. This event is represented in the data as multiple trades when
it is a single trade. Also, some timestamps have multiple quotes. To deal with
these events we performed trade and quote compacting. This is the process of
changing the data to better represent the fact that a single trade has
occurred, or to ensure that there is only a single quote per timestamp. For
quotes with the same timestamp, the most recent quote in the sequence of
quotes is kept and the rest are removed. For trades with the same timestamp,
and the same order type, the volume of the trade is calculated to be the
aggregated volume for that timestamp and the price of the trade is the volume
weighted average price.

###  F.2 Selection of moments

The first step in the construction of MM-SMD is the selection of the moments.
Winker et al. [64] recommend selecting a large number of moments taking into
account a large number of characteristics of the empirical data, which will
increase the probability of identification of the model parameters. The
moments that have been chosen are described in Table 8. These moments have
been chosen as they measure several stylised facts that we aim to reproduce. A
description of these stylised facts can be found in Table 7. We differ from
[34] as we do not include kurtosis as a moment. The fat tails of the log-
return distribution call into question the existence of such higher order
moments and kurtosis has been shown to not be a robust moment [64]. Jericevich
et al. [34] stated that it was the least reliable moment as large kurtosis
values were present when there were extreme price changes initially but the
price remained unchanged for the rest of the simulation. To measure the fat
tails of the return distribution the ”improved” Hill estimator was used to
measuring the tail-index of the upper 5% of the right tail of the log-return
distribution [64]. Removing the kurtosis moment is also valid as it still
ensures that we have more moments than parameters. The standard deviation of
the log-returns was kept as it is a robust statistic [64] and is far more
likely to exist when compared with the kurtosis. During the initial testing of
the model, the maximum likelihood estimate of the tail-index of the Pareto
distribution using the upper and lower 5% of the tails of the log-return
distribution was above 2 for all simulations, however, the estimate never went
above 4 again giving evidence for the existence of the second moment and the
nonexistence of the fourth.

The moments were estimated on the micro-price log-returns as Jericevich et al.
[34] found that using the mid-prices was unreliable when trying to identify
the distributional properties.

###  F.3 Simulated minimum distance

The calibration problem can be defined as the optimisation of a stochastic
approximation of an objective function
$\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}f(\boldsymbol{\theta})$ where
$\boldsymbol{\theta}$ is a parameter vector, $\boldsymbol{\Theta}$ is the
space of all possible parameters and f is the objective function.

The parameters that we will be optimising are described in Table 5. We have
that $\boldsymbol{\theta}=\\{N^{c},N^{f},\delta,\kappa,\nu,\sigma_{f}\\}$ are
the free parameters of the model. There are several parameters in the model
that were fixed to ensure that the fewest and most influential parameters were
the only ones that varied. Keeping the number of free parameters as low as
possible reduces the dimensionality of the surface and makes it easier to
optimise. The maximum and minimum forgetting factors of the uniform
distribution, $\lambda_{min}$ and $\lambda_{max}$, that were sampled were kept
fixed. These parameters aimed to establish a range of forgetting factors where
each chartist, when sampling from the uniform distribution, had an equal
chance of having a long or short trading horizon. The smallest forgetting
factor was set low enough such that the moving average tracked the mid-price
almost perfectly and the largest forgetting factor kept the moving average
reasonably constant across the simulation. Since these values have already
fulfilled their role there is no need to calibrate them. The time to trade,
after initialisation, was set at 25 seconds because it was deemed to give the
most similar order counts when compared with the day of data we calibrated to
(see Section 4.1). The time to cancellation parameter for the high-frequency
agent limit orders was set at 1 second as this was the fastest time to
cancellation that did not remove liquidity too quickly to cause consistent
liquidity crashes. The initial mid-price was arbitrarily set at 10000 and the
minimum volume for the limit orders, $x_{m}$, was reduced by half when
compared to [34] to account for the increase in the number of limit orders due
to the high-frequency agents trading in every event loop.

The MM-SMD method considers: the true moments, estimated moments and then
moments estimated from the simulated paths. The true moments,
$\boldsymbol{m}=(m_{1},...,m_{d})$, are the true moments of the $d$ selected
moments of the log-returns or the features. These values can’t be known and
they must be estimated from a particular realisation of a price path. The
estimated moments $\boldsymbol{m^{e}}=(m_{1}^{e},...,m_{d}^{e})$ are estimated
from the empirical intraday price path collected from the data source.
Finally, moments
$\left[\boldsymbol{m^{s}|\theta}\right]=(m_{1}^{s},...,m_{d}^{s})$ can be
estimated from the simulated price path of the model given a particular set of
parameters $\boldsymbol{\theta}$. Using the standard method of moments
condition
$\mathbb{E}\left[\left(\boldsymbol{m^{s}|\theta}\right)\right]=\boldsymbol{m}$.
If we know $\boldsymbol{m}$ then
$\mathbb{E}\left[\left(\boldsymbol{m^{s}|\theta}\right)-\boldsymbol{m}\right]=0$
so that from estimating the expectation over $I$ simulation runs we have that
$\frac{1}{I}\sum^{I}_{i=1}\left[\left(\boldsymbol{m_{i}^{s}|\theta}\right)-\boldsymbol{m}\right]=0$.
However, we don’t have the true moments $\boldsymbol{m}$, so instead, we use
the estimated moments $\boldsymbol{m^{e}}$. We define an error measure in
terms of the distance between the simulated moments and the estimated moments:

$\displaystyle\hat{G}(\boldsymbol{\theta})=\frac{1}{I}\sum^{I}_{i=1}\left[\left(\boldsymbol{m_{i}^{s}|\theta}\right)-\boldsymbol{m^{e}}\right]$
(18)

Here $I=5$ for the entire calibration process following [53]. The Monte Carlo
replications are key to reducing the variance of the stochastic approximation
of the objective function. Moreover, for each simulation run, we have that the
seed was also set to $i$. The seed controls the fundamental prices and the
forgetting factors for each of the fundamentalists and chartists respectively.
Averaging over multiple seeds is key as we want to try to average out the
effects of the random fundamental prices and forgetting factors. Five Monte
Carlo replications are low and this value was used due to computational
constraints. Please see Table 6 that shows the runtimes to get an idea of the
computational burden. In the ideal setting, each seed should be averaged over
far more than 5 times to reduce the variance caused by the randomness in each
simulation as well as the random fundamental prices and forgetting factors.

We can construct the full objective function as a weighted sum-of-squares of
the errors defined in $G$:
$f(\boldsymbol{\theta})=\hat{G}(\boldsymbol{\theta})^{T}\boldsymbol{W}\hat{G}(\boldsymbol{\theta})$.
Here $\boldsymbol{W}$ is a matrix of weights. Using the method of Winker et
al. [64] we choose
$\boldsymbol{W}=\mathbb{C}\mathrm{ov}\left[\boldsymbol{m^{e}}\right]^{-1}$.
The intuition behind the choice is that it takes into account the uncertainty
of the empirically estimated moments, $\boldsymbol{m^{e}}$, by assigning
larger weights to moments with greater uncertainty [53].

The matrix $\boldsymbol{W}$ is estimated using a moving block bootstrap with a
window of size 2000 [34]. The moving block bootstrap is used to preserve the
autocorrelations in the time series. The total number of moving block
bootstrap samples that were obtained was 1000. Therefore, we have 1000 samples
in which the empirical moments could be estimated, and 1000 samples of each
moment to estimate $\boldsymbol{W}$. Since we are inverting a covariance
matrix the condition number needs to be checked to determine the invertibility
of the matrix. The condition number of
$\mathbb{C}\mathrm{ov}\left[\boldsymbol{m^{e}}\right]$ is
$1.9536096786628117\times 10^{13}$. This is concerning however it is far less
than [34] who’s condition number was $8.375811152908104\times 10^{19}$.
Therefore the matrix may only just be invertible but it is at least comparable
and better than that of [34].

###  F.4 Optimisation

The objective function described in Section F.3 is computationally expensive
to evaluate. A single run of the simulation lasts 25 seconds and because the
number of Monte Carlo replications was set at 5, this means that it takes
approximately 125 seconds to evaluate the objective function. Therefore, the
Nelder-Mead (NM) optimisation algorithm [46] was chosen as it reduces the
number of objective function evaluations by choosing efficient search steps
[25]. The roughness of the stochastic approximation of the objective function
will cause lots of local minima. To avoid the rapid descent into local minima
the NM optimisation algorithm is combined with the threshold accepting
algorithm (TA) to give us Nelder-Mead with threshold accepting (NMTA) used in
[25, 34]. The NM algorithm has been shown to become inefficient in higher
dimensions and therefore the NM algorithm implementation with adaptive
parameters was chosen to increase the efficiency of the optimisation [20].

Given $n$ free parameters that need to be optimised the NM algorithm creates a
simplex of $n+1$ vertices, where each vertex represents a potential solution
to the optimisation. At each iteration, the simplex is transformed by a
reflection and either an expansion, contraction or shrinkage operator. This is
repeated until either the convergence criteria are met or until the
optimisation is halted by a time limit or a fixed number of iterations. As
noted previously this algorithm will converge quickly to local optima if the
objective surface is rough.

The TA algorithm selects a single point on the objective function and perturbs
it to create a new point. This point is then evaluated and if it has a better
objective function value then it is accepted as the new solution. Furthermore,
if it has a worse objective function value by no more than some threshold
$\tau$ of the current solution then the new solution replaces the current
solution [64]. The accepting of worse solutions allows this algorithm to
escape from local optima. The threshold, $\tau$, is decreased throughout the
optimisation so that convergence can occur. Each value for $\tau$ is used for
a set number of steps to allow the algorithm to search the parameter space
using each threshold several times [34].

Combining these two algorithms gives the NMTA algorithm. At each iteration
either a single step in the TA algorithm is executed or a single step in the
NM algorithm is executed. A step of either the TA or NM algorithm is performed
with probability $\xi$ and $1-\xi$ respectively with NM being chosen with a
higher probability. The TA algorithm is applied to the vertex with the lowest
objective function value and the size of the perturbation is dependent on the
mean value, computed over the simplex, of the parameter being chosen [34].
Moreover, for each of the operations in the simplex search the objective
function values are also compared using the same TA thresholding to ensure
that worse solutions are selected in the NM algorithm.

Initial calibration attempts caused parameters of the model to fall outside
reasonable bounds. For, example some vertices of the simplex would contain
negative values for the number of liquidity takers. This is unrealistic and
bounds on the parameters would need to be placed so that the simplex does not
search areas of the objective function that have no grounds in reality. Table
9 describes the parameter bounds and the justification for the bounds.

Table 9: Describes the bounds on the parameters of the model and gives a justification as to why the particular bounds were chosen. Parameter Bounds |  Justification
---|---
$N^{c}_{{}_{\mathrm{LT}}}$, $N^{f}_{{}_{\mathrm{LT}}}\in\left[1,20\right]$ |  The lower bound is used to ensure that there are always liquidity takers. The upper bound was set such that there is a case for more liquidity takers than providers but it was not set too high that the program would drop messages from the matching engine.
$\delta\in\left[0,10\right]$ |  Assumes the fundamental prices and the moving averages would never be more than 10 times that of the mid-price. This is a reasonable assumption due to the bounds on the variance on the fundamental prices and the method used to sample the forgetting factors. During calibration, this parameter never went outside this range.
$\kappa\in\left[0,\infty\right)$ |  Was allowed to vary on the positive real line.
$\nu\in\left[1.5,\infty\right)$ |  The lower bound was adjusted from 1 to 1.5 as values too close to 1 would cause volumes so high that there were errors when converting the sampled volume from a float to an integer. These volumes were deemed unrealistic.
$\sigma_{f}\in\left[0,0.05\right]$ |  Ensures that the fundamental prices are within approximately 10% of the mid-price. It was found that increasing this further did not make a difference as this value was rarely truncated and the prices did not move more than 10% over the course of the calibrated simulations.

Due to the bounds on the parameters the initial vertex was set at the mean
values for each of the parameters and an affine transformation was applied to
this vertex to give the initial simplex. The simulation was run for 100
iterations as it was determined that the simplex vertices converged to a state
where the parameters in the vertices had stopped changing in a meaningful way
at this point. The objective function was also shown to have converged to an
optima in this time.

###  F.5 Calibration convergence results

The plot on the left of Figure 29 depicts the convergence of the objective
function for the best vertex over the number of iterations (blue) as well as
the estimate for the scaled standard deviation of the objective function
values over the entire simplex (purple). The purple line can be thought of as
the convergence criteria of the objective function because as this decreases
this means that the simplex values are becoming more similar in their
objective function values and are therefore likely closer together on the
objective function. The plot on the right shows the convergence of the
objective function for each of the vertices on the simplex. The objective
function value for the best vertex does seem to have converged to an optima
because there has only been a very small improvement in the last 32
iterations. The last 35 iterations had no thresholding applied. The
convergence criteria have also decreased as the simplex vertices seem to have
more similar objective function values. This is also shown in the right-hand
side plot of the figure. Since the objective surface is stochastic and most
likely rough, there will always be noise in the objective function values even
if they are very close in the parameter space. Therefore, just looking at the
objective function convergence will not provide the full picture of
convergence. The key question is whether the vertices in the simplex are
similar enough where any of the simplex transformations or threshold accepting
perturbations will cause the simplex to continue to move around the parameter
space.

Figure 29: The left plot shows the convergence for the best vertex over the
number of iterations (blue) as well as the estimate for the scaled standard
deviation of the objective function values of the entire simplex (purple). The
purple line can be thought of as the convergence criteria of the objective
function because as this decreases this means that the simplex vertices’
objective function values are becoming more similar and are therefore likely
closer on the objective surface. The plot on the right shows the convergence
of the objective function for each of the vertices on the simplex. The
objective function value for the best vertex does seem to have converged to an
optima because there has only been a very small improvement in the last 32
iterations. The last 35 iterations had no thresholding applied. The
convergence criteria have also decreased as the simplex values seem to have
more similar objective function values. This is also shown in the right-hand
side plot of the figure.

Figure 30 shows the parameter value for each of the vertices in the simplex
(dashed lines) and the parameter value for the centroid of the simplex (solid
line). These plots show that throughout the optimisation we have that each
parameter value on each vertex has converged to the same place on the
objective surface because by the end of the calibration all the simplex
vertices have essentially the same parameter values. This means that if the
optimisation process was run for further iterations there would be no
significant changes in the vertices of the simplex. The simplex search
transformations are based on the difference in the parameter values of the
simplex and the centroid and because they are so similar there would not be
significant changes to the vertices and the algorithm would become a threshold
accepting algorithm. The threshold accepting steps are run with probability
$\xi$ which is much smaller than $(1-\xi)$ and therefore will only be run very
few times and the optimisation will become very inefficient to continue the
search as the simplex will change very slowly if at all for the rest of the
algorithm. For these reasons, the optimisation was halted at 100 iterations.
It can also be noted that there is variance in the objective function values
of each vertex in the simplex, in the last few iterations, even though the
parameter values are not significantly different. Highlighting the stochastic
and rough nature of the objective surface and its difficulty to optimise.

Figure 30: To determine if the parameters converged we plot trace plots for
the parameters. The parameter value for each of the vertices in the simplex
(dashed lines) and the parameter value for the centroid of the simplex (solid
lines) are plotted in these trace plots. Throughout the optimisation, we have
that the simplex has converged on the objective surface because by the end of
the calibration all the simplex vertices have essentially a single parameter
value. This means that if the optimisation process was run for further
iterations there would be no significant changes in the vertices of the
simplex.

The right-hand plot in Figure 29 has large spikes in the value of the
objective function. It was found that these spikes are due to the shrink
operator of the simplex algorithm. For a deterministic smooth function, the
shrinkage operator will shrink the simplex closer to the minimum as it assumes
it is closer to an optima. However, for a stochastic approximation of a
function shrinking the simplex could produce vertices that have much worse
objective function values due to the variance of the objective function and
its rough surface, especially if the simplex vertices still have significant
differences.

##  Appendix G Learning agent results

###  G.1 Best bid volume state spaces

Table 10: Describes how the states, for the best bid volume ($v_{k}^{b}$), are
determined as well as the approximate probability of being in each of those
states, based on the observed best bid volume at some event $k$ and the
simulated historical distributions. Table 10(a) and 10(b) shows how the volume
states are assigned for the agents with the smaller and larger state spaces
respectively.

(a) Smaller volume state space assignment: ($n_{T},n_{I},n_{S},n_{V}=5$). State | Volume Range | Prob.
---|---|---
$v$ | $\cdot<v_{k}^{b}\leq\cdot$ |
1 | 0 | 31 | 0.2
2 | 31 | 266 | 0.2
3 | 266 | 1453 | 0.2
4 | 1453 | 5209 | 0.2
5 | 5209 | $\cdot$ | 0.2

(b) Larger volume state space assignment:($n_{T},n_{I},n_{S},n_{V}=10$). State | Volume Range | Prob.
---|---|---
$v$ | $\cdot<v_{k}^{b}\leq\cdot$ |
1 | 0 | 11 | 0.1
2 | 11 | 31 | 0.1
3 | 31 | 93 | 0.1
4 | 93 | 266 | 0.1
5 | 266 | 636 | 0.1
6 | 636 | 1453 | 0.1
7 | 1453 | 2930 | 0.1
8 | 2930 | 5209 | 0.1
9 | 5209 | 2322 | 0.1
10 | 12322 | $\cdot$ | 0.1

###  G.2 Spread state spaces

Table 11: Describes how the states, for the spread ($s_{k}$), are determined
as well as the approximate probability of being in each of those states, based
on the observed spread at some event $k$ and the simulated historical
distributions. Table 11(a) and 11(b) shows how the spread states are assigned
for the agents with the smaller and larger state spaces respectively.

(a) Smaller spread state space assignment: ($n_{T},n_{I},n_{S},n_{V}=5$). State | Spread Range | Prob.
---|---|---
$s$ | $\cdot<s_{k}^{b}\leq\cdot$ |
1 | 0 | 1 | 0.6
2 | 1 | 2 | 0.123
3 | 2 | 3 | 0.062
4 | 3 | 7 | 0.091
5 |  7 | $\cdot$ | 0.124

(b) Larger spread state space assignment:($n_{T},n_{I},n_{S},n_{V}=10$). State | Spread Range | Prob.
---|---|---
$s$ | $\cdot<s_{k}^{b}\leq\cdot$ |
1 | 0 | 1 | 0.6
2 | 1 | 2 | 0.123
3 | 2 | 3 | 0.062
4 | 3 | 4 | 0.031
5 | 4 | 5 | 0.031
6 | 5 | 6 | 0.031
7 | 6 | 8 | 0.031
8 | 8 | 12 | 0.031
9 | 12 | 19 | 0.031
10 |  19 | $\cdot$ | 0.031

###  G.3 Convergence results

Figure 31: Trace plot for the actions taken in each state during training for
agents with $X_{0}=21500$ and $n_{T},n_{I},n_{S},n_{V}=5$ (left) and
$n_{T},n_{I},n_{S},n_{V}=10$ (right).

Figure 32: Trace plot for the actions taken in each state during training for
agents with $X_{0}=86000$ and $n_{T},n_{I},n_{S},n_{V}=5$ (left) and
$n_{T},n_{I},n_{S},n_{V}=10$ (right).

###  G.4 Behaviour analysis

(a) Plots the policy for the agent with $X_{0}=21500$ and
$n_{T},n_{I},n_{S},n_{V}=5$

(b) Plots the action volume trajectory for the agent with $X_{0}=21500$ and
$n_{T},n_{I},n_{S},n_{V}=5$.

Figure 33: Plots the policy (left) and the action volume trajectory (right)
for the RL agent with $X_{0}=21500$ and $n_{T},n_{I},n_{S},n_{V}=5$. From the
change in the action volume trajectory, it is clear that the agent has learned
to change its behaviour from an initial random agent.

(a) Plots the policy for the agent with $X_{0}=86000$ and
$n_{T},n_{I},n_{S},n_{V}=5$.

(b) Plots the action volume trajectory for the agent with $X_{0}=86000$ and
$n_{T},n_{I},n_{S},n_{V}=5$.

Figure 34: Plots the policy (left) and the action volume trajectory (right)
for the RL agent with $X_{0}=86000$ and $n_{T},n_{I},n_{S},n_{V}=5$. From the
change in the action volume trajectory, it is clear that the agent has learned
to change its behaviour from an initial random agent.

Figure 35: Plots the policy for the agent with $X_{0}=21500$ and
$n_{T},n_{I},n_{S},n_{V}=10$ Figure 36: Plots the action volume trajectory for
the agent with $X_{0}=21500$ and $n_{T},n_{I},n_{S},n_{V}=10$.

Figure 37: Plots the policy for the agent with $X_{0}=43000$ and
$n_{T},n_{I},n_{S},n_{V}=10$. Figure 38: Plots the action volume trajectory
for the agent with $X_{0}=43000$ and $n_{T},n_{I},n_{S},n_{V}=10$.

Figure 39: Plots the policy for the agent with $X_{0}=86000$ and
$n_{T},n_{I},n_{S},n_{V}=10$. Figure 40: Plots the action volume trajectory
for the agent with $X_{0}=86000$ and $n_{T},n_{I},n_{S},n_{V}=10$.

###  G.5 Learning dynamics

(a) Plots the average action for the high (resp. low) volume and narrow (resp.
wide) spread states over the training episodes. This figure plots the results
for the agent with $X_{0}=43000$ and $n_{T},n_{I},n_{S},n_{V}=5$.

(b) Plots the average action for the high (resp. low) volume and narrow (resp.
wide) spread states over the training episodes. This figure plots the results
for the agent with $X_{0}=43000$ and $n_{T},n_{I},n_{S},n_{V}=10$.

Figure 41: Plots the average action for the high (resp. low) volume and narrow
(resp. wide) spread states over the training episodes. This figure plots the
results for the agents with $X_{0}=43000$. The smaller state space agent is
shown on the left (Figure 41(a)) and the larger state space agent is shown on
the right (Figure 41(b))

(a) Plot the average action taken when the agent has lots of inventory at the
beginning of the simulation and when the agent has less inventory at the end
of the simulation. This figure plots the results for the agent with
$X_{0}=43000$ and $n_{T},n_{I},n_{S},n_{V}=5$.

(b) Plot the average action taken when the agent has lots of inventory at the
beginning of the simulation and when the agent has less inventory at the end
of the simulation. This figure plots the results for the agent with
$X_{0}=43000$ and $n_{T},n_{I},n_{S},n_{V}=10$.

Figure 42: Plots the average action taken when the agent has lots of inventory
at the beginning of the simulation and when the agent has less inventory at
the end of the simulation. This figure plots the results for the agents with
$X_{0}=43000$. The smaller state space agent is shown on the left (Figure
42(a)) and the larger state space agent is shown on the right (Figure 42(b))

Figure 43: Plots the average action taken by the RL agents with
$n_{T},n_{I},n_{S},n_{V}=5$, when there is lots of time remaining and little
time remaining as well as when the agent has a large amount of inventory
remaining and a small amount of inventory remaining. Once again we see that
the agents are learning to reduce the amount they trade as they get trained.
Moreover, we see that all the agents have learned to trade, on average, more
when there is less time remaining and trade more when there is less inventory
remaining. It is clear that the amount of inventory remaining has a more
significant effect when compared with the time remaining, as there is a larger
difference between the average actions taken in the more and less, time and
inventory remaining states.

Figure 44: Plots the average action taken by the RL agents with
$n_{T},n_{I},n_{S},n_{V}=10$, when there is lots of time remaining and little
time remaining as well as when the agent has a large amount of inventory
remaining and a small amount of inventory remaining. Here, we see that the
larger state space agents have also learned to trade less as the number of
episodes increase. In contrast to the smaller states space agents (Figure 43),
the time and inventory effects are roughly similar, therefore the effects are
cancelled out in the interaction plot.

###  G.6 Stylised fact dynamics

Figure 45: Compares the ABM, and JSE moments to the moments computed on data
generated by the simulations of the ABM with the RL agents. The moments shown
here are the ADF statistic, Kolmogorov-Smirnov (KS) statistic, mean, and
standard deviation. Here, we see that the mean and the KS statistic generated
by the RL agents across their learning episodes do not significantly differ
from these same moments generated by the e-ABM. The standard deviation (Std)
moment was zero for e-ABM and the JSE, up to the third significant digit,
however, we do see that some simulations saw larger variations in the log-
returns. The ADF statistic sees more significant differences between the e-ABM
and the e-ABMs with the RL agents. However, most of the moments fall into
either the JSE or e-ABM confidence intervals, and all are significantly
negative indicating that there insufficient evidence for the presence of a
unit root.

##  Appendix H Pseudocode

1

Input : LOB state $\mathcal{L}_{k}$ and simulation parameters
$\boldsymbol{\Theta}$

Output : Limit order $o_{ik}=(p_{ik},\omega_{ik},D_{ik})$

2 Function  _HighFrequencyAgentAction(__)_ :

3 if  _simulation is over_ then return

4 if  _simulation is not initialising_ then cancel orders that have not been
matched in $\phi$ time

5 Sample $u\sim U(0,1)$

6 if  _$u <(\rho_{t}+1$)/2_ then $D_{ik}=$ sell

7 else $D_{ik}=$ buy

8 if  _$D_{ik}=$ sell_ then

9 $\alpha=1-\rho_{t}/\nu$

10 Sample $\eta\sim\Gamma(s_{k},e^{\rho_{t}/\kappa})$

11 $p_{ik}=b_{k}+1+\lfloor\eta\rfloor$

12 $\omega_{ik}\sim f(x;5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

13 else

14 $\alpha=1+\rho_{t}/\nu$

15 Sample $\eta\sim\Gamma(-s_{k},e^{\rho_{t}/\kappa})$

16 $p_{ik}=a_{k}-1-\lfloor\eta\rfloor$

17 $\omega_{ik}\sim f(x;5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

18 end if

19 if  _initialising_ then

20 if  _$p_{k}$ is not inside initial spread_ then Submit order $o_{ik}$

21

22 else

23 Submit order $o_{ik}$

24 end if

25

26 end

Algorithm 1 Liquidity Provider Action Algorithm

1

Input : LOB state $\mathcal{L}_{k}$ and simulation parameters
$\boldsymbol{\Theta}$

Output : Market order $o_{ik}=(0,\omega_{ik},D_{ik})$

2 Function  _ChartistAgentAction(__)_ :

3 if  _simulation is over or initialising_ then return

4 $\bar{m}_{i,k}=\bar{m}_{i,k}+\lambda_{i}(\bar{m}_{i,k}-m_{k})$

5 if  _$\bar{{m}}_{i,k} >m_{k}+1/2s_{k}$_ then

6 $D_{ik}$ = sell

7 else if  _$\bar{{m}}_{i,k} <m_{k}-1/2s_{k}$_ then

8 $D_{ik}$ = buy

9 else

10 return

11 end if

12 $x_{m}=20$

13 if  _$|m_{k}-\bar{m}_{i,k}| >\delta m_{k}$_ then $x_{m}=50$

14 if  _$D_{ik}=$ sell_ then

15 $\alpha=1-\rho_{t}/\nu$

16 $\omega_{ik}\sim f(x_{m};5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

17 else

18 $\alpha=1+\rho_{t}/\nu$

19 $\omega_{ik}\sim f(x_{m};5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

20 end if

21 if  _there are LOs on the contra side and the order won’t cause a
volatility auction_ then

22 Submit order $o_{ik}$

23 else

24 return

25 end if

26

27 end

Algorithm 2 Chartist Agent Action Algorithm

1

Input : LOB state $\mathcal{L}_{k}$ and simulation parameters
$\boldsymbol{\Theta}$

Output : Market order $o_{ik}=(0,\omega_{ik},D_{ik})$

2 Function  _FundamentalistAgentAction(__)_ :

3 if  _simulation is over or initialising_ then return

4 if  _$f_{i} <m_{k}-1/2s_{k}$_ then

5 $D_{ik}$ = sell

6 else if  _$f_{i} >m_{k}+1/2s_{k}$_ then

7 $D_{ik}$ = buy

8 else

9 return

10 end if

11 $x_{m}=20$

12 if  _$|m_{k}-\bar{m}_{i,k}| >\delta m_{k}$_ then $x_{m}=50$

13 if  _$D_{ik}=$ sell_ then

14 $\alpha=1-\rho_{t}/\nu$

15 $\omega_{ik}\sim f(x_{m};5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

16 else

17 $\alpha=1+\rho_{t}/\nu$

18 $\omega_{ik}\sim f(x_{m};5,\alpha)$ Power law distribution with tail-index
parameter $\alpha$ (equation 1)

19 end if

20 if  _there are LOs on the contra side and the order won’t cause a
volatility auction_ then

21 Submit order $o_{ik}$

22 else

23 return

24 end if

25

26 end

Algorithm 3 Fundamentalist Agent Action Algorithm

1

Input : LOB state $\mathcal{L}_{k}$, simulation parameters
$\boldsymbol{\Theta}$, RL agent parameters
$\boldsymbol{\Theta}_{{}_{\mathrm{RL}}}$, and RL event messages from previous
orders, previous state $S_{k-1}$, and previous action $A_{k-1}$

Output : Market order $o_{ik}=(0,\omega_{ik},D_{ik})$

2 Function  _RLAgentAction(__)_ :

3 if  _simulation is over or initialising_ then return

$v_{k}$, $R_{k}$ = ProcessMessages(messages) // previous trades volume traded
and profit received

4 remainingInventory = remainingInventory - $v_{k}$

5 remainingTime = $T_{0}-\mathrm{currentTime}$

6 $s_{k}$, done = GetState($\mathcal{L}_{k}$, remainingTime,
remainingInventory)

7 if  _rl trading time is up and final trade submitted_ then state = terminal
state

8 if  _a previous trade was submitted_ then

9
$Q(S_{k-1},A_{k-1})=Q(S_{k-1},A_{k-1})+\alpha(R_{k}+\gamma\max\limits_{a}Q(S_{k},a)-Q(S_{k-1},A_{k-1}))$

10

11 end if

$A_{k}$ = GetAction($Q,S_{k},\epsilon$) // Sample action from
$\epsilon$-greedy policy

12 if  _rl trading time is up_ then

13 if  _still has inventory remaining_ then

$A_{k}$ = 2 // Submit largest trade possible

14

15 else

$A_{k}$ = 0 // Ensure no trade

16

17 end if

18

19 end if

20 $D_{ik}$ = sell

21 $\omega_{ik}$ = min($A_{k}(X_{0}/N_{dp})$, remainingInventory)

22 if  _there are LOs on the contra side and the order won’t cause a
volatility auction_ then

23 Submit order $o_{ik}$

24 else

25 return

26 end if

27 $S_{k-1}$ = $S_{k}$

28 $A_{k-1}$ = $A_{k}$

29 end

Algorithm 4 Learning Agent Action Algorithm

1

Input : Simulation parameters $\boldsymbol{\Theta}$, RL agent parameters
$\boldsymbol{\Theta}_{RL}$, and CoinTossX trading gateway

2 Function  _Simulate(__)_ :

3 Start the CoinTossX trading session

4 Initialise the LOB $\mathcal{L}_{0}$ with the initial spread $s_{0}$, best
bid and offer $b_{0}$ and $a_{0}$, as well as the mid-price $m_{0}$

5 Connect to the UDP socket with the IP and port of the CoinTossX market data
feed

6 Start the asynchronous task that will listen for and store CoinTossX
messages

7 Initialise the LPs, LTs and the RL agents (if it is an RL simulation)

8 Initialise the LOB using the initial LOB parameters and the LPs

9 Loop  __ // Event loop

10 if  _there are messages from CoinTossX_ then

11 for  _message in messages_ do

12 UpdateLOBState($\mathcal{L}_{k}$,message)

13 if  _there are RL traders_ then store the RL trade messages

14

15 end for

16 Pass updated LOB to all the agents for them to trade off

17

18 else

19 break

20 end if

21

22 EndLoop

23 Close asynchronous task resources, the UDP socket and close the CoinTossX
trading session

24 Write required data to files

25

26 end

Algorithm 5 Simulation Loop Algorithm

1

Input : Calibrated simulation parameters $\boldsymbol{\Theta}$ and RL agent
parameters $\boldsymbol{\Theta}_{rl}$

Output : Training results for an RL agent

2 Function  _TrainRLAgent(__)_ :

3 Initialise $Q_{0}$

4 Login to CoinTossX

5 for  _i in $\\{1,\ldots,\mathrm{number\;\;of\;\;episodes}\\}$_ do

6 Simulate the e-ABM with the RL agent

7 Update $\epsilon$ using the $\epsilon$ decay function in Figure 1

$Q_{i-1}=Q_{i}$ // Store previous Q-matrix for next iteration

8 Manual garbage collection

9

10 end for

11 Logout of CoinTossX

12 Close CoinTossX

13 Write the RL results to files

14

15 end

Algorithm 6 Learning Agent Training Algorithm

1

Input : The number of spread states $n_{{}_{S}}$ and the historical spread
distribution $\boldsymbol{spreads}$

Output : $n_{{}_{S}}$ quantile values used to assign discrete spread states

2 Function  _SpreadStateSpaceCreation(__)_ :

3 splitIntervalFactor = $n_{{}_{S}}$

4 loop  __

5 for  _p in range(0.6,1,splitIntervalFactor)_ do

6 store quantile($\boldsymbol{spreads}$, p)

7 end for

8 if  _number of unique quantiles is $n_{{}_{S}}$_ then

9 break

10 else

11 splitIntervalFactor = splitIntervalFactor + 1

12 end if

13

14 EndLoop

15

16 end

Algorithm 7 Spread State Table Creation Algorithm (Table 11)